Fastest way to extract differing positions from each alignment in a BAM file

[CHRYSES]

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.

[adaptivegenome]

samtools and bamtools both provide very fast APIs you can use, however this requires a minimal experience with a programming language or script...

[CHRYSES]

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.

[adaptivegenome]

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++...

[swbarnes2]

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.

[CHRYSES]

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...