SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Obtaining unique sequence tag file from fastQ format ramadatta.88 Introductions 0 09-26-2011 01:25 AM
Converting Solexa FASTQ file to unique sequence tags DrD2009 Bioinformatics 16 08-08-2010 11:30 PM
Solexa - same sequence but unique identifier Layla Bioinformatics 5 11-27-2009 05:08 AM
PubMed: Mining SNPs from DNA Sequence Data; Computational Approaches to SNP Discovery Newsbot! Literature Watch 0 09-22-2009 02:00 AM
In Sequence: Illumina GA’s Deep Coverage Shown to Be Useful for Profiling, Discovery Newsbot! Illumina/Solexa 0 02-26-2008 02:20 PM

Reply
 
Thread Tools
Old 02-09-2011, 06:17 AM   #21
quinlana
Senior Member
 
Location: Charlottesville

Join Date: Sep 2008
Posts: 119
Default

Quote:
Originally Posted by wdt View Post
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
quinlana is offline   Reply With Quote
Old 02-28-2011, 01:44 PM   #22
gpcr
Member
 
Location: usa

Join Date: May 2010
Posts: 18
Default 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?.
gpcr is offline   Reply With Quote
Old 03-01-2011, 10:35 AM   #23
epigen
Senior Member
 
Location: Germany

Join Date: May 2010
Posts: 101
Default get discordant pairs

Quote:
Originally Posted by gpcr View Post
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.
epigen is offline   Reply With Quote
Old 03-01-2011, 04:50 PM   #24
gpcr
Member
 
Location: usa

Join Date: May 2010
Posts: 18
Default

thanks @epigen
gpcr is offline   Reply With Quote
Old 05-18-2011, 11:42 PM   #25
plichel
Junior Member
 
Location: Germany

Join Date: Mar 2010
Posts: 9
Default 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 ?

Last edited by plichel; 05-19-2011 at 08:28 AM.
plichel is offline   Reply With Quote
Old 05-20-2011, 07:15 AM   #26
epigen
Senior Member
 
Location: Germany

Join Date: May 2010
Posts: 101
Default identify breakpoints with softclipped reads

Quote:
Originally Posted by plichel View Post
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?
epigen is offline   Reply With Quote
Old 05-26-2011, 10:43 AM   #27
plichel
Junior Member
 
Location: Germany

Join Date: Mar 2010
Posts: 9
Default

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 ?
plichel is offline   Reply With Quote
Old 07-25-2011, 11:47 PM   #28
tez
Junior Member
 
Location: Australia

Join Date: Jul 2011
Posts: 4
Default

The CREST algorithm described here:
http://www.nature.com/nmeth/journal/...meth.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.
tez is offline   Reply With Quote
Old 09-02-2011, 09:52 AM   #29
idan.gabdank
Junior Member
 
Location: California

Join Date: Jul 2011
Posts: 1
Default

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?
idan.gabdank is offline   Reply With Quote
Old 12-14-2011, 06:46 AM   #30
aquinom85
Research Bioinformaticist
 
Location: Boston

Join Date: Dec 2011
Posts: 19
Default

Does Hydra report zygosity of SVs it calls?
aquinom85 is offline   Reply With Quote
Reply

Tags
hydra, structural-variation

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:13 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO