This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Efficient Collapsing Of Bam Reads By Sequence

Is there an efficient utility out there for collapsing BAM files by sequence? I.e. keep only one of each sequence read (ideally with some constraints on which quality score reads to keep when there are multiple identical sequence reads with distinct quality scores)? thanks.

To clarify, I'd like to be able to only remove duplicates if their sequences are identical - so keep reads with the same alignment position if they have distinct sequences.

rna-seq sequence sam samtools

3 answers

You might also look in to the Broad's ReduceReads program.

http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_compression_reducereads_ReduceReads.html

The downside to using this is that the BAM files generated could only be used with GATK's tools. And (correct me if I'm wrong) but I thought this was part of the GATK v2 code that isn't open.

At least there is a specification.

Well you have the rmdup command of samtools:

Remove potential PCR duplicates: if multiple read pairs have identical external coordinates, 
only retain the pair with highest mapping quality.

If I am not mistaken, the creator of this tool (Li Heng, user lh3) actually recommends using Picard MarkDuplicates instead of samtools rmdup.

Interesting, I did not know that.

indeed, it all depends what the OP really needs - people often use the two concepts interchangeably (identical reads vs reads that map to the same location) though I agree that these are not the same

I've seen several references to rmdup being not recommended by its own creators etc. but if it's not recommended, why is it part of samtools still? If the goal is to remove PCR duplicates, then why isn't Picard's program identical to samtools's rmdup? What does it mean for one to be better than the other, if they are advertised to do the same thing, which is remove PCR duplicate reads (reads on same strand, position, regardless of sequence)?

I had better experience with Picard's MarkDuplicates than with samtools rmdup: http://picard.sourceforge.net/command-line-overview.shtml#MarkDuplicates

However, I did not find exact explanation how it works. I noticed in IGV that it removed exactly same reads (sequences) at exactly same positions and that samtools flagstat got better, but did not deeply examined what else it did though. Hope this helps.

If they remove the same reads, then what made the experience with MarkDuplicates better?

Log in to answer this question.