FeatureCounts output does not contain all the genes in the annotation file

Question

[devbt15]

Dear all, I use Subread in Bash for mapping and FeatureCounts for count quantification for my RNASeq data. I have used the same pipeline for other RNASeq datasets in the past and FeatureCounts gives count output for all the genes in the annotation file.

Recently I performed RNASeq 150 bp paired-ends reads in rice (Illumina). However, this time, after running featureCounts the count output file contained only 24000 genes out of total 42000 genes in the annotation file. I cannot understand why is it so? Does it suggest that mapping failed for the 18000 ignored genes? Does it suggest PCR amplification-based duplicated reads?

Can someone explain why is it so?

Regards.

[viDub]

Hi, Did you finally find out why it is so? I am running into the same problem with my own data and I can't find where this comes from...

[lieven.sterck]

do you mean that those genes are completely missing from the output or that they have 0 count ?

[devbt15]

Thank you for the response. They are completely missing from the output. I also ran HTSeq-counts and it also did not provide the count data for all the genes.

[ATpoint]

Please provide code.

[devbt15]

Here it is:

Subread index: subread-buildindex -F -o subreadindex ricegen.fa

Annotation file: gffread -E riceann.gff3 -T -o riceann.gtf

Subread mapping: subread-align -a riceann.gtf -F GTF -T 8 --sortReadsByCoordinates --multiMapping -i subreadindex -t 0 -r Q1read1.fq.gz -R Q1read2.fq.gz -o Q1.bam

FeatureCounts: featureCounts -p -T 8 -G ricegen.fa -t exon -g gene_id -a riceann.gtf -F GTF -M -O --fraction -o STRcounts.txt Q1.bam Q2.bam Q3.bam V1.bam V2.bam V3.bam

HTSeq-counts: htseq-count -t exon -i gene_id -f bam -m intersection-nonempty Q1.bam Q2.bam Q3.bam V1.bam V2.bam V3.bam riceann.gtf > result_file.txt

ricegen.fa and riceann.gtf are the genome and annotation files used.

[cpad0112]

Did you check the chromosome names in reference vs GTF?

[devbt15]

Yes, the names are the same. Anyways it would have been highlighted while running subread if it was not the same.

[cpad0112]

is RNAseq library prep strand specific?

[devbt15]

Yes, it is strand-specific.

[cpad0112]

How is read coverage for missing genes per exon (for eg take one and check)? Please also post featurecounts summary, if possible.

[devbt15]

Thanks, I will check the bigwig on the genome browser and let you know.

[devbt15]

Here is the FeatureCounts output screenshot image for the 6 bam files as input in the code above: https://ibb.co/PFf9WxQ

I also checked the read coverage for gene IDs missing from the FeatureCounts output vs annotation file. Here is the screenshot image: https://ibb.co/SmD3cK2

[Hyper_Odin]

Can you also post the snapshot of the gtf file?

[devbt15]

Please find it here

https://ibb.co/KK1kwb3

[devbt15]

Hi viDub, Unfortunately not. Since I have little time left for this work, I shifted my analysis from SubRead to Salmon based reference transcriptome analysis where I do not need the FeatureCounts or HTSeqCount. I am sorry, but I cannot help you anymore at this point in time. Regards.

Answer 1

[el.dansch]

Did you check if the missing genes have a CDS annotated, or do they maybe have missing gene_ids (or whatever you specified to count) in the GTF file?

Answer 2

[Beáta]

I ran into the same problem. Apparently the missing genes are all single-exon genes, represented in the gtf file by a single row, as "gene" only, no exons. Since featurecounts in default mode specifies "exon" as the feature for counting, these genes were left out of the counts file. If the feature is specified as "gene", these "no-exon-specified" genes are also counted (metafeature can be left as "gene_id"). I would recommend counting in two steps: first, by default setting (exon as feature), second, by "gene" for feature. Extract the counts for the "no-exon-specified" genes from the second counts file, and concatenate with the first counts file.

[lieven.sterck]

Not immediately convinced this is a valid approach.

in my book, counting on 'gene' will also include reads that map in introns for instance (because the gene is from start gene to stop gene). As a consequence you will over-count, as typically you do not want to include reads in introns as expression signal.

In any way it's clear this is a GFF/GTF issue as even those single exon genes should have an exon (and CDS) annotated & included in the GFF file

Perhaps try to run something like AGAT on it , it may be able fix this.

[ATpoint]

That is a poor GTF file then...