Dear Community,

I'd like to clarify unexpected behaviour of featureCounts when counting on ONT cDNA bam files, generated with minimap2 with the splice alignment preset. FeatureCounts version: 2.1.1. When running

featureCounts -L -C -Q 10 --primary -M -O --fraction -t transcript


aired-end : no                                               ||
||        Count read pairs : no                                               ||
||              Annotation : genes.filtered.gtf (GTF)                         ||
||      Dir for temp files : /scratch/local/snakepipes.khA6aeMgrY             ||
||                                                                            ||
||                 Threads : 8                                                ||
||                   Level : meta-feature level                               ||
||      Multimapping reads : counted (fractional)                             ||
||     Multiple alignments : primary alignment only                           ||
|| Multi-overlapping reads : counted                                          ||
||   Min overlapping bases : 1                                                ||
||          Long read mode : yes                                              ||
||                                                                            ||
\\============================================================================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file genes.filtered.gtf ...                                ||
||    Features : 267506                                                       ||
||    Meta-features : 67928                                                   ||
||    Chromosomes/contigs : 439                                               ||
||                                                                            ||
|| Process BAM file sorted.bam...          ||
||    Single-end reads are included.                                          ||
||    Total alignments : 9121703                                              ||
||    Successfully assigned alignments : 3639166 (39.9%)                      ||
||    Running time : 1.37 minutes                                             ||

i.e. ~ 40% of the fragments are counted. In the summary file, I notice 4481575 fragments (~ 49%) are Unassigned:Mapping Quality.

Samtools stats run on the same bam file on the genic regions tells me that <12% of alignments have MAPQ<10.

What could be the reason for the discrepancy between low MAPQ fraction calculated by samtools stats and by featureCounts ?

Best wishes,

Katarzyna