[quinlana]

Hi,Quinlana;
I have no idea about the last step yet,in your paper workflow,incluing use BWA NOVOALIGN ,HYDRA and BEDTOOLS ,in http://code.google.com/p/hydra-sv/wiki/TypicalWorkflow linkage,I think it is not a good guidelines

Thank you for letting me know. When I find some time I will try to make the approach more intuitive.

[zee]

Perhaps a python/shell script or a Makefile would simplify things. Any plans for this in the future? Perhaps I could help out on this task since we are working on improving use cases for our aligner.

[quinlana]

Hi Zee,
Thanks for the suggestion. However, I am reticent to automate anything about the process as it very much depends on the quality and variability of the user data, the organism, etc. I think it's best for me to merely improve the documentation and give thorough explanations of why things are done the way they are.

Aaron

[quinlana]

Hi all,
I just posted a script on the Hydra site that allows one to convert the Hydra breakpoint calls (in BEDPE format) to BED12 for visualization on IGV, UCSC, etc. I've had several requests for such a tool and finally got around to doing it.

http://code.google.com/p/hydra-sv/downloads/detail?name=bedpeToBed12.py

Best,
Aaron

[wdt]

Aaron,

This is a samtools question but since the o/p is going into Hydra, I thought of posting here.

Following the workflow on http://code.google.com/p/hydra-sv/wiki/TypicalWorkflow, for extracting the discordant reads using

samtools view -uF 2 sample.tier1.bam | \
bamToFastq -bam stdin \
-fq1 sample.tier1.disc.1.fq \
-fq2 sample.tier1.disc.2.fq

Should we see exact same number of reads (and identical pairwise read Ids) in the 1.fq and 2.fq file?

For 1.fq and 2.fq files I have, I don't see corresponding match of read IDs. Do
you require BAM to be sorted based on read id?


Thanks in advance.

[quinlana]

Aaron,

This is a samtools question but since the o/p is going into Hydra, I thought of posting here.

Following the workflow on http://code.google.com/p/hydra-sv/wiki/TypicalWorkflow, for extracting the discordant reads using

samtools view -uF 2 sample.tier1.bam | \
bamToFastq -bam stdin \
-fq1 sample.tier1.disc.1.fq \
-fq2 sample.tier1.disc.2.fq

Should we see exact same number of reads (and identical pairwise read Ids) in the 1.fq and 2.fq file?

For 1.fq and 2.fq files I have, I don't see corresponding match of read IDs. Do
you require BAM to be sorted based on read id?


Thanks in advance.

Assuming you have used BWA for Tier 1, your Tier 1 BAM file should be in "query order". That is, the order of the alignments in the BAM file should be in the order of the input FASTQ files. Using the settings described in the workflow, you should have a single alignment for each end of every pair. Thus, when creating 1.fq and 2.fq, you should have the exact same number of reads in each. If not, I suspect you have either a) used a different aligner for Tier1 or b) not used a "query-ordered" BAM file.

I have updated the workflow to indicate that bamToFastq expects query-ordered BAM files.

Best,
Aaron

[gpcr]

missing reads

Aaron,
I have merged bam file resulting from multiple lanes of paired end alignments. when I extracted fastq from the alignemnt . I have unequal reads in pair1 and 2. when i examined the reads in bam file. I could see for some, there is only one mate (either fwd oir reverse) only aligned and other missing. Do I need to append missing read from the raw read lane? or exclude them from the analysis?.

[epigen]

get discordant pairs

Aaron,
I have merged bam file resulting from multiple lanes of paired end alignments. when I extracted fastq from the alignemnt . I have unequal reads in pair1 and 2. when i examined the reads in bam file. I could see for some, there is only one mate (either fwd oir reverse) only aligned and other missing. Do I need to append missing read from the raw read lane? or exclude them from the analysis?.

If you extract the discordant pairs from your BAM file like that:

samtools view -hb -F 1038 orig.bam > discordant.bam

you get reads that have neither flag 2 (proper pair) nor flag 4 (read itself unmapped) nor flag 8 (mate unmapped) nor flag 1024 (is duplicate).
If the number of reads1 and reads2 is still not equal, maybe your aligner messed up? As Aaron wrote above, you need to have exactly one alignment per read.

By the way, the above also works for coordinate-sorted BAMs. Afterwards you just have to namesort discordant.bam with samtools sort -n option.

[gpcr]

thanks @epigen

[plichel]

breakpoints

Does hydra also utilize partially mapped reads ?
I see a lot of softclipped alignments in my bwa aligned sam file. I am wondering whether this information is used when searching for breakpoints.

Edit1:
It seems that softclipped reads are implicitly used, since their edit distance is often abnormal after clipping. (depending from which end it is clipped, the outer or the inner distance changes...)

To me it seems: simply by looking at softclipped reads it might be possible to detect the exact breakpoint position (at single nucleotide resolution). Why nobody use it ? Do I miss something ?

[epigen]

identify breakpoints with softclipped reads

To me it seems: simply by looking at softclipped reads it might be possible to detect the exact breakpoint position (at single nucleotide resolution). Why nobody use it ? Do I miss something ?

That's a very good question. Maybe softclipped reads don't give enough evidence?

[plichel]

I guess it will greatly depend on coverage and uniqueness of the alignment.
Suppose you have 30x and see indeed at a particular position 30 ore more softclipped reads supporting the breakpoint and the mapper/aligner doesnt report multiple hits. What could lead here to a false positive conclusion ?

[tez]

The CREST algorithm described here:

http://www.nature.com/nmeth/journal/vaop/ncurrent/full/nmeth.1628.html

Uses soft-clipped reads, but requires reads of 75bp or longer (so is out for my SOLiD project). It doesn't appear to use discordance information though, so there may be some benefit in using multiple tools as a part of a workflow.

[idan.gabdank]

I am trying to run with the typical workflow of hydra, but after tier 3 I am getting stuck with the script "pairDiscordants.py".

I am getting:

Traceback (most recent call last):
File "/usr/local/bin/pairDiscordants.py", line 294, in <module>
sys.exit(main())
File "/usr/local/bin/pairDiscordants.py", line 290, in main
pairReads(opts.inFile, opts.numMappings, opts.order, opts.dist, opts.minSpan, opts.minConcRange, opts.maxConcRange, opts.mode, opts.anchorThresh, opts.multiThresh, opts.editSlop)
File "/usr/local/bin/pairDiscordants.py", line 28, in pairReads
printHydraMappings(pairs, editDistance, editSlop)
UnboundLocalError: local variable 'editDistance' referenced before assignment

The input file looks like:

CHROMOSOME_II 9172809 9172833 000_1000_1326_R3/2 0 -
CHROMOSOME_IV 1641291 1641314 000_1000_171_F3/1 0 +

The command was:

> cat result | pairDiscordants.py -i stdin -m hydra -z 800

What am I doing wrong?

[aquinom85]

Does Hydra report zygosity of SVs it calls?