I am trying to understand the effects of sorting on MarkDuplicates and how to interpret the numbers in the featureCounts log.

According to the MarkDuplicates documentation, it handles coordinated-sorted and query-sorted BAM inputs slightly differently:

The program can take either coordinate-sorted or query-sorted inputs, however the behavior is slightly different. When the input is coordinate-sorted, unmapped mates of mapped records and supplementary/secondary alignments are not marked as duplicates. However, when the input is query-sorted (actually query-grouped), then unmapped mates and secondary/supplementary reads are not excluded from the duplication test and can be marked as duplicate reads.

Q1: Should the sorting in the MarkDuplicates step affect featureCounts? It sounds like query-sorted may remove some additional _unmapped_ mates and secondary/supplementary reads. According to the featureCounts documentation, it does not count any multi-mappings by default. And I would imagine quantification only works on mapped reads and not unmapped either, so I thought the featureCounts result should be the same no matter how MarkDuplicates was performed. But when I ran featureCounts on both types of MarkDuplicates sorting outputs, there are some slight difference in the counts files, albeit not a lot. Why is that?

Q2: I am having some trouble seeing where the numbers from the featureCounts logs come from. The BAM file below was generated by MarkDuplicates using coordinate-sorted mode.

samtools flagstat output for BAM input to featureCounts:

35483648 + 0 in total (QC-passed reads + QC-failed reads)
24895287 + 0 primary
10588361 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
20155191 + 0 mapped (56.80% : N/A)
9566830 + 0 primary mapped (38.43% : N/A)
24895287 + 0 paired in sequencing
12334176 + 0 read1
12561111 + 0 read2
9539764 + 0 properly paired (38.32% : N/A)
9539764 + 0 with itself and mate mapped
27066 + 0 singletons (0.11% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

MarkDuplicates log that generated the deduped BAM file summarized above:

LIBRARY:            Unknown Library
UNPAIRED_READS_EXAMINED:    395123 
READ_PAIRS_EXAMINED:        17983810
SECONDARY_OR_SUPPLEMENTARY_RDS: 10588361
UNMAPPED_READS:         15328457
UNPAIRED_READ_DUPLICATES:   368057
READ_PAIR_DUPLICATES:       13213928
READ_PAIR_OPTICAL_DUPLICATES:   0
PERCENT_DUPLICATION:        0.736906
ESTIMATED_LIBRARY_SIZE:     4893981

featureCounts log (relevant parts):

//========================== featureCounts setting ===========================\\
||                                                                            ||
||                 Threads : 1                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : yes                                              ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||   Min overlapping bases : 1                                                ||
||                                                                            ||
||          Chimeric reads : not counted                                      ||
||        Both ends mapped : not required                                     ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file /home/chan/mice_microgravity/genomes/transcripts/ ... ||
||    Features : 663359                                                       ||
||    Meta-features : 43432                                                   ||
||    Chromosomes/contigs : 45                                                ||
||                                                                            ||
|| Process BAM file /home/chan/mice_microgravity/outputs_coordinate/picar ... ||
||    Paired-end reads are included.                                          ||
||    Assign fragments (read pairs) to features...                            ||
||                                                                            ||
||    WARNING: reads from the same pair were found not adjacent to each       ||
||             other in the input (due to read sorting by location or         ||
||             reporting of multi-mapping read pairs).                        ||
||                                                                            ||
||    Read re-ordering is performed.                                          ||
||                                                                            ||
||    Total fragments : 17999315                                              ||
||    Successfully assigned fragments : 2902259 (16.1%)                       ||
||    Running time : 1.30 minutes                                             ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

featureCounts output:

Assigned            2902259
Unassigned_Ambiguity        132370
Unassigned_MultiMapping     5895062
Unassigned_NoFeatures       1602957
Unassigned_Unmapped     7466667
Unassigned_MappingQuality   0
Unassigned_FragmentLength   0
Unassigned_Chimera      0
Unassigned_Secondary        0
Unassigned_Nonjunction      0
Unassigned_Duplicate        0

In the featureCounts log, it says the total fragments is 17999315 (2902259 + 132370 + 5895062 + 1602957 + 7466667, of which 2902259 is assigned), but not sure where this number came from.

I think I see where the unmapped fragments (Unassigned_Unmapped) is from. Since this is paired end data, 7466667 * 2 = 14933334 unmapped alignments = total - mapped - singleton - UNPAIRED_READ_DUPLICATES = 35483648 - 20155191 - 27066 - 368057 (looking at the samtools flagstat and MarkDuplicates outputs)

But not sure about the remaining mapped fragments: 2902259 + 132370 + 5895062 + 1602957 = 10532648 * 2 = 21065296 mapped alignments. But looking at the samtools flagstat output, there's only 20155191 mapped alignments, or counting singletons, 20155191 + 27066 = 20182257 alignments, which is still less than the ones reported by featureCounts. Where is featureCounts getting its total fragment count from?