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%)?