Marking optical or PCR duplicates with picard vs. samtools flagstat

Question

I am trying to identify the best library prep method for noninvasive prenatal test samples, to be sequenced with Illumina Novaseq.

One metric that I am evaluating is the number/ percentage of PCR or optical duplicates. In order to achieve that I have done the following steps:

1. aligned the fastq reads against the reference genome with bwa mem
2. the .sam output from bwa mem was converted to bam, sorted and indexed
3. the duplicates were marked with the command MarkDuplicates from picard

Then if I call samtools flagstat on the sorted bam file which had the duplicates marked with picard I get:

26595942 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
466 + 0 supplementary
1636809 + 0 duplicates
24969064 + 0 mapped (93.88% : N/A)
26595476 + 0 paired in sequencing
13297738 + 0 read1
13297738 + 0 read2
21102678 + 0 properly paired (79.35% : N/A)
24049752 + 0 with itself and mate mapped
918846 + 0 singletons (3.45% : N/A)
462212 + 0 with mate mapped to a different chr
182138 + 0 with mate mapped to a different chr (mapQ>=5)

Which sounds a bit strange to me, as no duplicates were found.

However, the MarkDuplicates command also outputted a metrics file (dups/dupsMetrics.txt) that did identify some duplicates. This was achieved with the following command:

java -jar picard.jar MarkDuplicates I=alignment/sample.sorted.bam O=alignment/sample_markDup.bam M=dups/dupsMetrics.txt

I see that the value under the column PERCENT_DUPLICATION in dups/dupsMetrics.txt is 0.065555 (given the column name, I would expect that the value is indeed 0.06%, and not 6.55%)

What would be the explanation between the discrepancy in the output of samtools flagstat and MarkDuplicates. Is it reasonable to expect that the value displayed by samtools flagstat is rounded to zero, given the low percentage which was observed by MarkDuplicates (0.06%)?

Answer 1

There are duplicates, in this line:

1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356

According to samtools documentation for flagstat:

Provides counts for each of 13 categories based primarily on bit flags in the FLAG field. Each category in the output is broken down into QC pass and QC fail. In the default output format, these are presented as "#PASS + #FAIL" followed by a description of the category.

So you can usually take whatever is in the first column

Answer 2

@StupidWolf's answer is correct -- that first number in the flagstat output is what you want to look at to see the number of reads marked as duplicates. I wanted to add that the number given in the Picard metrics file is in fact 6%, not 0.06%. This is due to an ambiguity that is widespread in Picard metrics files; in many places the program emits a fraction even when the column name is "percent of something". Soo Hee Lee wrote a very short GATK blog post about this way back in the day, see the archived version here.