Hello there,
I my RNA-seq analysis in Galaxy, I aligned my reads using STAR with the UCSC hg38 reference genome (hg38.knownGene.gtf.gz ). After quantification with featurecounts and running DESeq2, My IDs were in UniProt format, so I converted them to gene symbols.
However, I noticed a large number of genes still had NA values after conversion. When I manually checked some of these IDs, I found that they included transcripts, UniParc entries, and some merged identifiers.
Also my DESeq2 results contained more than 200,000 rows, which made me suspect that featureCounts may not have been working properly.
So, I repeated the entire pipeline, this time using the GENCODE v50 annotation file (gencode.v50.annotation.gtf.gz) in STAR and selecting the built-in reference in featureCounts. This time, DESeq2 produced ~30,000 rows, which seems much more reasonable. The gene identifiers were Entrez IDs, which I converted to gene symbols, but I still ended up with ~700 NA values.
My question is:
What was likely wrong in my initial workflow that led to the inflated gene count and annotation issues? Could the main issue have been not using the built-in annotation/reference in featureCounts?
Also why do I still get NA values?
In previous analyses, I ran featureCounts using the same reference I used for STAR and I don’t recall encountering this problem.
Many things can go wrong, and it is hard to say what went wrong in your situation without looking at the data.
Annotations of human genes are comprehensive, and, unless you are after chimeric transcript in cancer or something similar, where novel genes/transcripts are expected, I recommend any canonical gene annotation. Do not create/update the gene annotation data. It is OK to use the same annotation file in mapping and counting step. You can also do not use the annotation file in the mapping step.
DESeq2 returns number of features in the read count files. featureCounts output depends on input data and settings. If you updated the annotation file, e.g., RNA_STAR followed by STringTie followed StringTie_merge, the new annotation may have many new gene annotations, depending on data.
Read counts works as following: reads are counted against features defined in 3rd column of the annotation file, for example, exon, and aggregated against features described in the last column (attributes), for example, gene_id. If you select another attribute, like transcript_id (if it is present in the annotation file), you’ll get counts for transcripts (more lines in output). If you disable aggregation on attributes, you’ll get counts on exons. Changes in featureCounts job setup have very strong effect on output.
Consider using relevant GTN tutorials as a guide, for example: