bwa quality trimming and samtools rmdup

[greigite]

I am trying to use bwa and samtools to map 76 bp reads from multiple bacterial strains back onto a reference sequence, with the ultimate goal of extracting SNP frequencies. To obtain accurate variant frequencies, it is important to me to remove PCR duplicates. It occurred to me that quality trimming reads with BWA using the "-q" flag during alignment could affect how well rmdup works downstream. Say 2 reads are PCR duplicates, but one is rather low-quality and is trimmed to a different length than the other. The alignment start and stop positions would no longer be the same for these duplicates, and they would not be filtered using samtools rmdup, which requires identical external coordinates.

Is this right? Does BWA do hard trimming of reads with the "-q" flag? Or does it ignore low quality bases when calculating the alignment, but still use them to determine alignment coordinates?

While I'm at it, could anyone explain how the BWA quality trimmer works? Please don't say "read the man page," I did, and was very confused by the explanation there:

-q INT 	
Parameter for read trimming. BWA trims a read down to argmax_x{\sum_{i=x+1}^l(INT-q_i)} if q_l<INT
where l is the original read length. [0]

Noob question #3: Does rmdup work for single read data or doesn't it? The samtools man page states:

Samtools’ rmdup does not work for single-end data and does not remove duplicates across chromosomes. Picard is better.

However, there is an option "-s" to use rmdup on single read data, and when I apply it to an alignment, it looks to reduce the size of the resulting bam file. The man page and the command itself don't agree!

[xuer]

did you get the answer? nobody know the answer?

[greigite]

I did not get an answer- now I use Mosaik instead of bwa, it has much better documentation.

[Fedde]

Explanation of BWA read trimming

The BWA trimming feature seems to be explained a little more clearly here: http://seqanswers.com/forums/showthread.php?t=6251 . The real C source code is in the function bwa_trim_read() in the file bwaseqio.c, but I found the comments and variable names of the Perl example referenced in the other thread more clear.

[lbthrice]

Also see the SolexaQA FAQ for an enlightening discussion of the bwa algorithm vs. SolexaQA algorithm.

http://solexaqa.sourceforge.net/questions.htm#bwa

[dan]

I did not get an answer- now I use Mosaik instead of bwa, it has much better documentation.

Where is the real bottleneck ;-)

[mrood]

Duplicate Marking & Trimming

I am trying to use bwa and samtools to map 76 bp reads from multiple bacterial strains back onto a reference sequence, with the ultimate goal of extracting SNP frequencies. To obtain accurate variant frequencies, it is important to me to remove PCR duplicates. It occurred to me that quality trimming reads with BWA using the "-q" flag during alignment could affect how well rmdup works downstream. Say 2 reads are PCR duplicates, but one is rather low-quality and is trimmed to a different length than the other. The alignment start and stop positions would no longer be the same for these duplicates, and they would not be filtered using samtools rmdup, which requires identical external coordinates.

Is this right? Does BWA do hard trimming of reads with the "-q" flag? Or does it ignore low quality bases when calculating the alignment, but still use them to determine alignment coordinates?

While I'm at it, could anyone explain how the BWA quality trimmer works? Please don't say "read the man page," I did, and was very confused by the explanation there:

-q INT 	
Parameter for read trimming. BWA trims a read down to argmax_x{\sum_{i=x+1}^l(INT-q_i)} if q_l<INT
where l is the original read length. [0]

Noob question #3: Does rmdup work for single read data or doesn't it? The samtools man page states:

Samtools’ rmdup does not work for single-end data and does not remove duplicates across chromosomes. Picard is better.

However, there is an option "-s" to use rmdup on single read data, and when I apply it to an alignment, it looks to reduce the size of the resulting bam file. The man page and the command itself don't agree!

I have the same question. I know it has been a long time since you wrote this post, but I am wondering what you ended up doing. I am trying to develop a good pipeline for doing similar analyses.

[dan]

I have the same question. I know it has been a long time since you wrote this post, but I am wondering what you ended up doing. I am trying to develop a good pipeline for doing similar analyses.

Which question? He posted 3. I found a good description of the protocol here:

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_protocol