How to Annotate Your Own RNA-Seq Data

Oak
rna-seq annotation gencode single-cell bioinformatics
Glowing jellyfish-like creatures in distinct colors drift through a dark forest, each ringed by a hovering band of light.

By the end of this you will have a gene-level count/TPM matrix where every row carries a stable identifier, a gene symbol, a biotype, a chromosome, and a MANE Select transcript where one exists, plus a manifest file recording exactly which annotation release produced it. You will also have a functional layer (GO, Reactome, MSigDB Hallmark) attached by identifier rather than by symbol, and, if you have single-cell data, per-cell labels with confidence scores. You need FASTQs or a quantification directory, roughly 50 GB of free disk for a human index and genome, 16 GB of RAM for salmon (32 GB if you build a decoy-aware index), R ≥ 4.3 with tximport, GenomicFeatures, rtracklayer, and AnnotationDbi, and a shell with wget, samtools, and gffread. Everything below is human-specific and assumes GRCh38.

Annotation in this context means two different things that people run together. The first is structural: which intervals of the genome are genes and transcripts, and what are their identifiers. That choice changes your numbers before you interpret anything. The second is functional: what those genes are known to do. The first is a versioning problem. The second is a mapping and statistics problem.

1. Pick one annotation release and pin it

Three annotation sources dominate human work: GENCODE (identical to Ensembl’s human gene set for the basic/comprehensive gene models), RefSeq, and UCSC’s derived tables. They do not agree. GENCODE’s comprehensive set for a recent release contains roughly 62,000 gene entries and about 250,000 transcripts, of which only about 20,000 genes are protein-coding. RefSeq’s curated set is far smaller and more conservative. If you quantify the same library against both, the count for a given gene will differ, sometimes by a lot, because the transcript models differ in UTR length, in the number of retained-intron isoforms, and in whether an overlapping lncRNA exists to absorb reads.

We use GENCODE comprehensive (CHR, primary assembly only) for personal data, for three reasons. It includes lncRNA and pseudogene models, which matter because reads from those loci otherwise get misassigned to neighboring protein-coding genes. It carries tag "MANE_Select" on the transcripts that RefSeq and Ensembl agree on, which gives you a clean one-transcript-per-gene view when you want one. And its versioned identifiers (ENSG00000141510.17) make provenance unambiguous.

REL=47
mkdir -p ref && cd ref
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.primary_assembly.annotation.gtf.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.transcripts.fa.gz
sha256sum *.gz > annotation.sha256

Write the release number into a manifest now, not later. If you re-sequence in a year and quantify against release 51, the two matrices are not directly comparable at the transcript level and only approximately comparable at the gene level. Longitudinal comparison means re-quantifying the old FASTQs against the new annotation, which costs an hour of CPU and saves you from a fake effect.

2. Build a decoy-aware index

Selective alignment against the transcriptome plus a genome decoy avoids the main failure mode of naive transcriptome pseudoalignment: intronic and intergenic reads getting forced onto the nearest transcript.

zcat GRCh38.primary_assembly.genome.fa.gz | grep '^>' | cut -d ' ' -f1 | sed 's/>//g' > decoys.txt
cat gencode.v47.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz
salmon index -t gentrome.fa.gz -d decoys.txt -k 31 -p 16 --gencode -i salmon_idx_v47

--gencode splits the pipe-delimited GENCODE FASTA headers and keeps only the transcript ID. Skip it and every downstream join breaks in a confusing way, because your transcript IDs will look like ENST00000456328.2|ENSG00000290825.1|-|OTTHUMT.... Use -k 31 for reads ≥ 75 bp; drop to -k 23 or so for 50 bp reads or you will lose mapping rate.

3. Quantify with bias correction and posterior samples

salmon quant -i ref/salmon_idx_v47 -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  --seqBias --gcBias --posBias \
  --numGibbsSamples 20 \
  --validateMappings -p 16 \
  -o quant/sample01

Read quant/sample01/logs/salmon_quant.log before anything else. A poly-A selected library from whole blood should map at 80–90% against a decoy-aware index. Below 70%, suspect adapter contamination, rRNA carryover, or a species/build mismatch. The 20 Gibbs samples cost little and give you per-transcript inferential variance, which is the only honest way to know whether a transcript-level difference is real or just multi-mapping ambiguity. (For samples where two isoforms share nearly all their exons, the point estimate is close to meaningless without that variance.)

4. Build tx2gene from the GTF you used

Do not download a tx2gene table from somewhere else. Derive it from the same GTF file whose checksum you recorded.

library(rtracklayer); library(tximport); library(data.table)

gtf <- import("ref/gencode.v47.primary_assembly.annotation.gtf.gz")
tx  <- as.data.table(mcols(gtf[gtf$type == "transcript"]))

tx2gene <- tx[, .(TXNAME = transcript_id, GENEID = gene_id)]
files <- file.path("quant", list.files("quant"), "quant.sf")
names(files) <- list.files("quant")

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

lengthScaledTPM gives counts that are safe to hand to DESeq2 or limma-voom while still correcting for isoform-usage changes in effective gene length. If you want raw counts plus an offset instead, use the default "no" and pass txi directly to DESeqDataSetFromTximport, which handles the length matrix for you.

Version suffixes are the usual first crash. If your quant files carry ENST00000456328.2 and your tx2gene carries ENST00000456328, tximport errors out on missing transcripts. Either keep versions on both sides (preferred, since it proves they came from the same release) or pass ignoreTxVersion = TRUE. Keep versions.

5. Attach gene-level metadata from the same GTF

Symbols are labels, not keys. Build your table keyed on versionless Ensembl gene IDs and treat gene_name as display text.

gn <- as.data.table(mcols(gtf[gtf$type == "gene"]))
genes <- gn[, .(gene_id, gene_name, gene_type, seqnames = as.character(seqnames(gtf[gtf$type == "gene"])))]
genes[, ensg := sub("\\..*$", "", gene_id)]

# MANE Select transcript per gene, where one exists
mane <- tx[sapply(tx$tag, function(x) "MANE_Select" %in% x),
           .(ensg = sub("\\..*$", "", gene_id), mane_tx = transcript_id)]
genes <- merge(genes, mane, by = "ensg", all.x = TRUE)

Three things bite here. PAR_Y genes appear twice with _PAR_Y appended to the ID and the same symbol, so a naive symbol join silently duplicates about 45 rows on chrX/chrY. Symbols change between releases (the MARCH/MARCHF and SEPT/SEPTIN renamings broke a lot of pipelines, largely because spreadsheet software had been mangling them into dates for years). And roughly 1–2% of symbols in GENCODE map to more than one gene ID. Deduplicate explicitly, by summing counts or by keeping the highest-expressed ID, and write down which you did.

Filter after annotation, not before. A reasonable default for bulk tissue: keep genes with gene_type in protein_coding, lncRNA, and the IG/TR segment types, then apply edgeR::filterByExpr. Dropping pseudogenes before quantification would have pushed their reads onto real genes. Dropping them after is free.

6. Run the annotation-dependent QC checks

These checks exist to catch a wrong annotation or a mislabeled sample, and each one is a single line once the biotype column is attached.

Mitochondrial fraction: sum TPM over seqnames == "chrM". In whole blood this is typically a few percent. Very high values point to degraded RNA or a 3’-biased library.

Globin fraction: HBB, HBA1, HBA2, HBD. In whole blood without globin depletion these transcripts can consume well over half the library, which crushes effective depth for everything else. If you see that, treat low-expression genes in that sample as unreliable rather than as biologically low.

rRNA and MT-rRNA: gene_type == "rRNA" plus MT-RNR1/MT-RNR2. Poly-A selection should leave this low. Ribo-depletion libraries behave differently.

Sex check: XIST (ENSG00000229807) against RPS4Y1, UTY, DDX3Y. This catches sample swaps faster than any genotype concordance check and needs nothing but the expression matrix. If the RNA and the genome from the same person disagree here, stop and resolve it before interpreting anything.

Biotype composition: the fraction of the library in protein_coding should be stable across your own samples over time. A sudden shift usually means a library-prep change, not biology.

7. Map to functional sets by identifier

Over-representation tests on a thresholded gene list throw away most of the data and are sensitive to your cutoff. We prefer a ranked test on the full list, using shrunken log-fold-changes or a signed Wald statistic.

library(fgsea); library(msigdbr)

h <- msigdbr(species = "Homo sapiens", collection = "H")
pathways <- split(h$ensembl_gene, h$gs_name)

stats <- setNames(res$stat, sub("\\..*$", "", rownames(res)))
stats <- stats[!is.na(stats)]
fg <- fgseaMultilevel(pathways, stats, minSize = 15, maxSize = 500)

Use the Ensembl ID column from msigdbr, not the symbol column. Symbol-based joins lose 5–10% of genes silently and the loss is biased toward recently renamed and poorly characterized genes, which is exactly the set you would want to notice.

Pathway databases can also be built rather than borrowed. Projects have assembled organism-specific metabolic pathway databases directly from RNA-seq assemblies plus enzyme-function prediction, which is the route to take when no curated resource exists for your system1. For human personal data, curated resources exist and building your own is not worth the time. Curated multi-omics resources that tie transcript changes to protein and clinical layers are a better source of prior structure than any single-gene lookup2.

One level up from pathway enrichment is upstream regulator inference: asking which transcription factors or signaling nodes would explain the observed pattern. Integrating expression with other molecular layers is what makes those inferences hold up, as opposed to reading a TF’s own mRNA level, which is frequently uninformative3. Multi-omics integration with machine learning is the general direction personal profiling is heading, and expression alone is the weakest of the layers to reason from in isolation4.

8. Annotate the noncoding and isoform layers deliberately

Gene-level counts collapse a lot. Isoform switching, alternative polyadenylation, circular RNAs, and small RNA classes each carry signal that gene-level summarization erases, and the annotation you chose determines which of those you can see at all5. If you care about isoform usage, keep the transcript-level quant.sf files and the Gibbs samples, and use a method that models inferential replicates rather than treating transcript TPMs as observed data. If you care about lncRNAs, check whether the specific loci you are interested in are even present in your release before concluding they are not expressed.

9. Cell-type annotation for single-cell data

If your data is bulk, most of the variance between timepoints in blood is cell-composition change, not per-cell regulation. Deconvolution (CIBERSORTx with LM22, or BayesPrism with a matched single-cell reference) gives you fraction estimates you can regress out or model explicitly. Do this before you interpret any gene as “up”.

If you have single-cell data, the annotation problem is assigning a label to each cell after clustering. Three approaches, in increasing order of automation:

Marker-based. Score canonical marker sets per cluster (AddModuleScore in Seurat, sc.tl.score_genes in Scanpy) and label manually. Slow, fully auditable, and still the ground truth you check the others against.

Reference-based. SingleR against a labeled bulk or single-cell reference, or CellTypist with a pretrained model such as Immune_All_Low.pkl for PBMC data. Fast, and the per-cell probability output tells you where the reference did not cover your data.

import celltypist
pred = celltypist.annotate(adata, model="Immune_All_Low.pkl",
                           majority_voting=True, over_clustering="leiden")
adata = pred.to_adata(insert_prob=True)

Model-based. Transformer architectures trained across many datasets handle large atlases and generalize across batches better than nearest-centroid reference mapping; CIForm treats each cell’s expression vector as a sequence of sub-vectors and scales to datasets where reference-mapping methods slow down6. Cross-omics approaches using vector-quantized autoencoders extend the same idea to labeling cells in one modality from labels in another, which matters if you have scATAC or CITE-seq alongside RNA7.

Whatever you use, keep a confidence measure and do not accept labels for cells that sit below it. Methods that read out deep network training dynamics, rather than only final predictions, can flag which cells and which regions the model is uncertain about, and that uncertainty map is usually more informative than the labels themselves for finding a population the reference did not contain8. Rare cell types, doublets, and stressed cells all present as confidently mislabeled if you only look at the argmax.

Common problems

Chromosome naming mismatch. UCSC uses chr1, Ensembl uses 1. Mixing a UCSC genome with an Ensembl GTF produces zero counts everywhere and an alignment step that appears to succeed. Check with grep '^>' genome.fa | head against cut -f1 annotation.gtf | head.

Missing --gencode at index time. Transcript IDs come through as pipe-delimited monsters and tximport reports that none of your transcripts are in tx2gene. Rebuild the index rather than patching strings downstream.

Duplicated symbols after a symbol-keyed join. Always key on versionless ENSG. If a downstream tool demands symbols, convert at the last possible step and check sum(duplicated(symbols)) before you do.

Annotation drift across timepoints. Comparing a sample quantified against v43 with one against v47 will produce differences in lncRNA and pseudogene counts that look like biology. Re-quantify everything against one release whenever you change releases.

Low mapping rate with high duplication. Usually low-complexity library or over-amplification, not an annotation problem. Check salmon mapping rate, the duplication rate from fastp, and the insert-size distribution before touching the annotation.

Over-reading single genes. A 1.5-fold change in one transcript between two timepoints in one person, with no replicates, is not an observation you can act on. Personal RNA-seq is most useful for composition, for large coordinated pathway shifts, and for tracking the same features across many timepoints. If a result looks clinically meaningful (an immune signature, a marker associated with a disease process, anything you would want to act on), that is the point to bring the data and the methods to a physician or a clinical geneticist. Research-grade RNA-seq is not a diagnostic and the annotation choices above change the numbers enough that treating them as one would be a mistake.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Alex Van Moerkercke, Michele Fabris, Jacob Pollier, et al. CathaCyc, a Metabolic Pathway Database Built from Catharanthus roseus RNA-Seq Data. Plant and Cell Physiology, 2013. https://doi.org/10.1093/pcp/pct039 ↩

  2. Hao Luo, Yunhao Yang, Zhipeng Gong, et al. LOSTdb: a manually curated multi-omics database for lung cancer research. BMC Bioinformatics, 2025. https://doi.org/10.1186/s12859-025-06319-6 ↩

  3. J. Pei, M. Schuldt, E. Nagyova, et al. Multi-omics integration identifies key upstream regulators of pathomechanisms in hypertrophic cardiomyopathy due to truncating MYBPC3 mutations. Clinical Epigenetics, 2021. https://doi.org/10.1186/s13148-021-01043-3 ↩

  4. Hany E. Marei, Carlo Cenciarelli, David Vagni. Integration of artificial intelligence and multi-omics for precision medicine. Functional & Integrative Genomics, 2026. https://doi.org/10.1007/s10142-026-02024-6 ↩

  5. Wenbin Dao, Simeng Zhang, Tao Zhang, et al. Advances in Poultry RNA-Omics Research: Technologies, RNA Information Layers, and Applications in Complex Traits. Animals, 2026. https://doi.org/10.3390/ani16172700 ↩

  6. Jing Xu, Aidi Zhang, Fang Liu, et al. CIForm as a Transformer-based model for cell-type annotation of large-scale single-cell RNA-seq data. Briefings in Bioinformatics, 2023. https://doi.org/10.1093/bib/bbad195 ↩

  7. Han Peng, Wuchao Liu, Yifang Cai, et al. Effective and Robust Single-cell Cross-omics Annotation via Vector-Quantized Autoencoders. Proceedings of the 32nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining V.2, 2026. https://doi.org/10.1145/3770855.3818995 ↩

  8. Jonathan Karin, Reshef Mintz, Barak Raveh, et al. Interpreting single-cell and spatial omics data using deep neural network training dynamics. Nature Computational Science, 2024. https://doi.org/10.1038/s43588-024-00721-5 ↩