I have bulk rnaseq data and I am interested in 1 gene, it's transcripts and exons. I ran the rnaseq nf-core pipeline (v3.9) (which uses salmon) to generate gene counts and transcript counts. And then I ran featureCounts on the BAM files to generate exon counts.

I assume that total gene count == sum of all transcript counts == sum of all exon counts

These are the counts for this one gene. Values are mean of raw counts over n samples. Not normalised or scaled.

• Gene count: 11446
• Sum of all transcripts count: 11494
• Sum of all exons count:
  • M+O+f+: 31740
  • M+O-f+: 486
  • M-O+f+: 30775
  • M-O-f+: 484
  • M+O+f-: 31740
  • M-O+f-: 30775
  • M-O-f-: 1253
  • M+O-f-: 1309

Total gene count does closely match sum of all transcript counts but not the exon counts. There are 8 variations of exon counts which is explained below.

This is the script used for featureCounts. I have single-end reads.

module load subread/2.0.0

featureCounts \
-a "genes.gtf" \
-o "counts-exon.txt" \
-F "GTF" \
-t "exon" \
-g "exon_id" \
-s 0 \
-T 4 \
-f \
-M \
-O \
*.bam

I included or excluded the parameters -M, -O and -f

-O  Assign reads to all their overlapping meta-features 
-M  Multi-mapping reads will also be counted.
-f  Perform read counting at feature level 

which gives rise to the 8 variations of exon counts. But none of them matches the expected count of roughly 11500.

• I understand that gene count and sum of transcript counts are very close since they come from the same salmon pipeline. I expect the featurecounts based exon counts to be further away but not wildly different. Is there an explanation for this mismatch?
• Is there a parameter in featureCounts that is set incorrectly?
• Is there a way to generate exon counts using the rnaseq pipeline or salmon, so I don't need to use featureCounts?

Ideas, suggestions and solutions are appreciated.