This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Marking Duplicates With Molecular Tag (Umi)

I have some reads that contain a molecular index so I can know whether they are PCR duplicates. I am going to use the 0x400 flag as specified in the SAM spec to mark them as optical/PCR duplicates.

Should I mark all of the reads in a group (having the same POS and molecular tag) with that flag or should I leave the one with the highest quality (or by whatever metric) unmarked?

I will be sending the result to the GATK SNP-calling pipeline.

markduplicates index
• 6,662 views • link
updated by Ram Ram • written by brentp brentp

I'm not sure I understand:are your SAM records already marked with the flag 0x4 ? or are you looking for a method to set the flag according to the chrom/pos/your-index ?

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

No, they are mapped reads. I want to set the flag 0x400 (1024) to show that they are PCR duplicates.

• 1 views • link
written by brentp brentp

Hi Brent,

Did you ever get this working?

Is the code available for download somewhere? I need it but don't want to re-invent the wheel.

Thank you

• 0 views • link
updated by Ram Ram • written by Michele Busby Michele Busby

Hi Brent,

I know this is old, but in case anyone else needs this, I have code here to add the UMI to the bam file by reading information from an original FASTQ: https://github.com/mbusby/AddUMIsToBam in the RX and QX fields.

The Picard MarkDuplicates, I hear from that team but did not test myself yet, will handle this in its duplicate marking.

• 0 views • link
updated by Ram Ram • written by Michele Busby Michele Busby

1 answer

you have to leave one read with the highest quality and flag the others with the duplicate flag.

here is the code of picard markdup: all the reads are flagged but the best: (https://github.com/nh13/picard/blob/master/src/java/net/sf/picard/sam/MarkDuplicates.java )

for (final ReadEnds end : list) {
    if (end != best) {
        addIndexAsDuplicate(end.read1IndexInFile);
        addIndexAsDuplicate(end.read2IndexInFile);
    }
• 27 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

in the source it is stored in the member "store" 

private short getScore(final SAMRecord rec) {
        short score = 0;
        for (final byte b : rec.getBaseQualities()) {
            if (b >= 15) score += b;
        }
(...)
     pairedEnds.score += getScore(rec);

and then this score is used to get the best pair:

for (final ReadEnds end : list) {
            if (end.score > maxScore || best == null) {
                maxScore = end.score;
                best = end;
            }
        }
• 45 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

cheers, my attention span is too short for java.

• 0 views • link
written by brentp brentp

Log in to answer this question.

More posts like this