Marking optical or PCR duplicates with picard vs. samtools flagstat 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:


aligned the fastq reads against the reference genome with bwa mem
the .sam output from bwa mem was converted to bam, sorted and indexed
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&gt;=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%)?
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