How does DuplicateSetIterator choose the same representative in the duplicate sets for a read and its mate?

[limwz01]

I saw that DuplicateScoringStrategy allows SUM_OF_BASE_QUALITIES to be set as scoring strategy. However, I am confused because it only takes the base qualities for the record in question and ignores the mate. Wouldn't that mean that for paired-end reads the read chosen in the DuplicateSet may be different from the read chosen in the mate's DuplicateSet? I'm asking because UmiAwareMarkDuplicatesWithMateCigar uses this scoring strategy (inherited from MarkDuplicates). Hope someone could clarify! Thanks

@nh13 @jacarey

[yfarjoun]

If I'm not mistaken, it passes over the file twice. first pass it collects the information it needs to decide on the duplicates and representatives...and then on the second pass it marks and emits the reads.

[limwz01]

hi @yfarjoun , thanks for the response, I looked over the code before and I think it wasn't a two-pass algorithm. could you point me to the piece of code that does the two passes? thanks!

[yfarjoun]

First, on picard/sam/markduplicates/MarkDuplicates.java:270 it calls buildSortedReadEndLists which opens the inputs and iterates over the reads.

Then in the loop controlled by picard/sam/markduplicates/MarkDuplicates.java:338 the iterator (opened in picard.sam.markduplicates.util.AbstractMarkDuplicatesCommandLineProgram.SamHeaderAndIterator) is consumed and the reads marked using the information collected in the first pass.

[limwz01]

right, that's for MarkDuplicates. UmiAwareMarkDuplicatesWithMateCigar inherits from SimpleMarkDuplicatesWithMateCigar which overrides the doWork method with the DuplicateSetIterator approach.

[limwz01]

I just tested it and confirmed that it's a bug. Here's the test SAM file:

@HD	VN:1.5	SO:coordinate
@SQ	SN:chr1	LN:101
@SQ	SN:chr2	LN:101
@SQ	SN:chr3	LN:101
@SQ	SN:chr4	LN:101
@SQ	SN:chr5	LN:101
@SQ	SN:chr6	LN:101
@SQ	SN:chr7	LN:404
@SQ	SN:chr8	LN:202
@RG	ID:1AAXX.1	SM:Hi,Mom!	LB:mylib	PL:ILLUMINA
@PG	ID:MarkDuplicates	PN:MarkDuplicatesGATK	VN:1	CL:MarkDuplicates2 merge1.sam	PP:bwa
@PG	ID:bwa	PN:bwa	VN:1	CL:bwa aln
1AAXX.1.1	83	chr7	1	255	101M	=	302	201	CAACAGAAGCNGGNATCTGTGTTTGTGTTTCGGATTTCCTGCTGAANNGNTTNTCGNNTCNNNNNNNNATCCCGATTTCNTTCCGCAGCTNACCTCCCAAN	/EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	PG:Z:MarkDuplicates	RG:Z:1AAXX.1	MC:Z:101M	RX:Z:AAAAAAAA
1AAXX.1.2	83	chr7	1	255	101M	=	302	201	CAACAGAAGCNGGNATCTGTGTTTGTGTTTCGGATTTCCTGCTGAANNGNTTNTCGNNTCNNNNNNNNATCCCGATTTCNTTCCGCAGCTNACCTCCCAAN	EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	PG:Z:MarkDuplicates	RG:Z:1AAXX.1	MC:Z:101M	RX:Z:AAAAAAAA
1AAXX.1.1	163	chr7	302	255	101M	=	1	-201	NCGCGGCATCNCGATTTCTTTCCGCAGCTAACCTCCCGACAGATCGGCAGCGCGTCGTGTAGGTTATTATGGTACATCTTGTCGTGCGGCNAGAGCATACA	EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	PG:Z:MarkDuplicates	RG:Z:1AAXX.1	MC:Z:101M	RX:Z:AAAAAAAA
1AAXX.1.2	163	chr7	302	255	101M	=	1	-201	NCGCGGCATCNCGATTTCTTTCCGCAGCTAACCTCCCGACAGATCGGCAGCGCGTCGTGTAGGTTATTATGGTACATCTTGTCGTGCGGCNAGAGCATACA	/EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	PG:Z:MarkDuplicates	RG:Z:1AAXX.1	MC:Z:101M	RX:Z:AAAAAAAA

Running it as follows:

picard UmiAwareMarkDuplicatesWithMateCigar --INPUT test.sam --OUTPUT test_md.bam --METRICS_FILE metrics.txt --UMI_METRICS_FILE umi_metrics.txt

Results:

1AAXX.1.1	1107	chr7	1	255	101M	=	302	202	CAACAGAAGCNGGNATCTGTGTTTGTGTTTCGGATTTCCTGCTGAANNGNTTNTCGNNTCNNNNNNNNATCCCGATTTCNTTCCGCAGCTNACCTCCCAAN	/EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	MC:Z:101M	RG:Z:1AAXX.1	MQ:i:255	RX:Z:AAAAAAAA
1AAXX.1.2	83	chr7	1	255	101M	=	302	202	CAACAGAAGCNGGNATCTGTGTTTGTGTTTCGGATTTCCTGCTGAANNGNTTNTCGNNTCNNNNNNNNATCCCGATTTCNTTCCGCAGCTNACCTCCCAAN	EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	MC:Z:101M	RG:Z:1AAXX.1	MQ:i:255	RX:Z:AAAAAAAA
1AAXX.1.1	163	chr7	302	255	101M	=	1	-202	NCGCGGCATCNCGATTTCTTTCCGCAGCTAACCTCCCGACAGATCGGCAGCGCGTCGTGTAGGTTATTATGGTACATCTTGTCGTGCGGCNAGAGCATACA	EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	MC:Z:101M	RG:Z:1AAXX.1	MQ:i:255	RX:Z:AAAAAAAA
1AAXX.1.2	1187	chr7	302	255	101M	=	1	-202	NCGCGGCATCNCGATTTCTTTCCGCAGCTAACCTCCCGACAGATCGGCAGCGCGTCGTGTAGGTTATTATGGTACATCTTGTCGTGCGGCNAGAGCATACA	/EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE	MC:Z:101M	RG:Z:1AAXX.1	MQ:i:255	RX:Z:AAAAAAAA

[limwz01]

oh, the reason the bug triggers is not so trivial. in UmiGraph DuplicateSets using the default comparator is being created. These DuplicateSets will use the TOTAL_MAPPED_REFERENCE_LENGTH scoring strategy, which is fine since it uses both the read and its mate to score, as long as the orders of the records were the same when looking at a read region and the mate region. However, this latter assumption is not satisfied because by the time the position-deduplicated DuplicateSet is passed to the UmiGraph the order has already been altered by the other scoring strategy SUM_OF_BASE_QUALITIES that was used for position-based deduplication.

[limwz01]

I suggest:

1. change back the default scoring strategy to TOTAL_MAPPED_REFERENCE_LENGTH
2. allow users to specify their own symmetric scoring strategies by a custom SAM tag. so that we can calculate our own TOTAL_SUM_OF_BASE_QUALITIES which include the mate, for example
3. correct the bug in UmiGraph where the DuplicateSets are created with the default comparator