Hi all, 

I am having a weird issue (weird to me being a new R user) using featurecounts to generate counts for genes but also for individual exons. I have successfully mapped gene counts to their respective genes and the total reads mapped and annotated to genes match those generated by a bioinformatics department so its a good quality control. 

I also want to generate counts for individual exons so that I can test the knockdown of one of our genes at specific exons in a mouse model. The issue is that the exon counts do not add up or match I am seeing with total read mapped to genes. The code I am using is as follows:

For annotating genes:

fc <- featureCounts(bamDir, isGTFAnnotationFile = TRUE, annot.ext = "Mus_musculus.GRCm38.94.chr.gtf", isPairedEnd = FALSE)

Count result for gene of interest = 

This matches results obtained from bioinformatcis department  for the gene above as well as all other genes (with +/- of 5 counts, I am guessing for parameters they have used in their script which I am not )...

However, when I attempt to annotate exons specifically, I am running in to trouble for two reasons. The total counts generated and also the amount of exons being incorrect (possible duplicates?)

The code I am using is as follows:

fc <- featureCounts(bamDir, isGTFAnnotationFile = TRUE, annot.ext = "Mus_musculus.GRCm38.94.chr.gtf", isPairedEnd = FALSE, useMetaFeatures = FALSE) OR

fc <- featureCounts(bamDir, isGTFAnnotationFile = TRUE, annot.ext = "Mus_musculus.GRCm38.94.chr.gtf", isPairedEnd = FALSE, useMetaFeatures = "exons")

These both produce the following results: 

This is icnredible confusing as there are only 5 exons for this gene, and clearly the counts are far too low to reflect the total counts annotated to the gene as a whole.. Can anybody advice what I might be doing wrong here and why this result is being produced?

And finally, is there an intuitive way of plotting counts per exons in an easy to follow way, maybe thorugh ggplo2 just plotting counts on y axis by exon number on x axis? or a way of passing this festurecount result to ggbio?

Thank you for your help!

Yes as you suspect there should be duplicates in the annotated exons for this gene. You can check this by looking at the exon annotation data featureCounts extracted from the GTF file (the extracted annotation is saved in the component fc$annotation).

When you count reads at exon level, by default a read overlapping more than one exon is not counted. Therefore all duplicate exons will receive zero count. Also junction reads will not be counted either because they span more than one exon. However if you set allowMultiOverlap=TRUE you will get all these reads counted, but the total number of reads counted to exons will be greater than the total number of reads counted to the gene.

For gene level counting, reads overlapping duplicate exons or multiple different exons are always counted once only (the allowMultiOverlap parameter controls how to count those reads overlapping more than one gene instead of more than one exon).

You can use the flattenGTF function in Rsubread to flatten your annotation before you run featureCounts.