SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
extract alignment from SAM with a GFF file NicoBxl Bioinformatics 4 08-02-2011 01:45 PM
Extract perfectly mapped reads from SAM/BAM file Graham Etherington Bioinformatics 2 07-21-2011 07:27 AM
Extract Gap and Mismatch From MAF alignment Output peveralldubois Bioinformatics 0 01-14-2011 07:23 PM
BWA: specifying SAM/BAM file header fields before read alignment? nora Bioinformatics 3 12-04-2010 09:11 PM
Filter BAM records by positions using picard guavajuice Bioinformatics 0 04-02-2010 02:45 PM

Reply
 
Thread Tools
Old 12-14-2011, 05:47 AM   #1
CHRYSES
Member
 
Location: Netherlands

Join Date: Dec 2009
Posts: 13
Default Fastest way to extract differing positions from each alignment in a BAM file

Hi,

What would be the fastest way (I have to do this hundreds of millions times) to extract for each aligned read in a BAM file:
1) The positions where the read bases differ from a reference sequence.
2) The PHRED base quality values of these bases. If the difference is an indel, the quality value will, of course, be skipped.

As far as I know, I cannot use mpileup or anything I know of due to memory limitation as this is a very custom amplicon reference analysis, with >500 million coverage per base position on the reference amplicon.

In short, I need to apply an efficient approach to extract all differing positions for each aligned read.

Thanks.

Last edited by CHRYSES; 12-14-2011 at 05:52 AM.
CHRYSES is offline   Reply With Quote
Old 12-14-2011, 06:08 AM   #2
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

samtools and bamtools both provide very fast APIs you can use, however this requires a minimal experience with a programming language or script...
adaptivegenome is offline   Reply With Quote
Old 12-14-2011, 06:13 AM   #3
CHRYSES
Member
 
Location: Netherlands

Join Date: Dec 2009
Posts: 13
Default

Quote:
Originally Posted by genericforms View Post
samtools and bamtools both provide very fast APIs you can use, however this requires a minimal experience with a programming language or script...
Yeah, I tried to get into that, but I am not good at "C" language, I could not follow it. I hope someone else has created/thought of something...

I could go directly into text format (i.e. SAM) and parse it with PERL, but that's really very slooooooow.
CHRYSES is offline   Reply With Quote
Old 12-14-2011, 06:18 AM   #4
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

If you are going to examine every read in order then I suppose you could parse a giant text file. This sort of sequential analysis is hard to speed up unless you parallelize it. If you are parsing a text file you could break it into many parts and then run your PERL script on the many parts in parallel (if you have access to that kind of equipment).

I am not aware of an off the shelf tool. Sorry!

Personally I opt for pthreads and C/C++...
adaptivegenome is offline   Reply With Quote
Old 12-14-2011, 08:11 AM   #5
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Illumina reads are error prone. If you pull every single read with a discrepancy from reference, you are going to pull a lot of noise.

I don't think that a pileup can be generated with only variant positions, but you could grep the pileup to only get lines with alterante letters. The pileup will have the position, all the letters called by all the reads that cross the position, and all the qualities for all the reads that cross the position.
swbarnes2 is offline   Reply With Quote
Old 12-14-2011, 11:28 AM   #6
CHRYSES
Member
 
Location: Netherlands

Join Date: Dec 2009
Posts: 13
Default

Quote:
Originally Posted by swbarnes2 View Post
Illumina reads are error prone. If you pull every single read with a discrepancy from reference, you are going to pull a lot of noise.

I don't think that a pileup can be generated with only variant positions, but you could grep the pileup to only get lines with alterante letters. The pileup will have the position, all the letters called by all the reads that cross the position, and all the qualities for all the reads that cross the position.
Yes, but how can I run a pileup on a single position with 500 million X coverage ? I think i will need to do this read by read...
CHRYSES is offline   Reply With Quote
Reply

Tags
bam parse samtools

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 12:15 PM.


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