# How to Retrieve Gene Annotations from UCSC, Ensembl, and GENCODE in R

UCSC, Ensembl, and GENCODE solve different annotation problems. GENCODE supplies the canonical, manually curated gene models for human and mouse, Ensembl supplies a versioned database with stable identifiers and rich cross-references, and UCSC supplies browser tracks, conservation data, and coordinate conversion tools. The rule of thumb is straightforward. Start with GENCODE for human or mouse transcript-level work, switch to Ensembl for attribute-heavy queries or species outside the vertebrate mainstream, and use UCSC when browser tracks, repeats, or genome-build conversion enter the picture.

This article doubles as a compact Bioconductor tutorial for the annotation stack in R for bioinformatics. It shows retrieval code for all three sources, compares them on the criteria that decide real projects, and walks through four worked scenarios from RNA-seq to GWAS.

## What Gene Annotation Means in an R Workflow

Gene annotation is the mapping between identifiers and genomic features. A complete annotation tells you where a gene sits (chromosome, start, end, strand), how it is transcribed (transcripts, exons, coding sequence), what it is called (stable ID, symbol, synonyms), and what it is (biotype, meaning the classification into protein_coding, lncRNA, miRNA, and similar categories, plus cross-references to GO terms, PubMed, and orthologs, which are genes in other species descended from a common ancestor). In R these records arrive as GRanges objects from the GenomicRanges package, as TxDb or EnsDb containers, as OrgDb organism databases, or as plain data frames parsed from GTF and GFF3 files. GTF (Gene Transfer Format) and GFF3 (General Feature Format) are tab-delimited text standards with one genomic feature per row and an attribute column holding key-value pairs.

The source matters because identifiers are messy. A survey of protein-coding gene nomenclature pulled from Ensembl release 113 found that symbols and synonyms range from 16,600 in pig to 64,580 in mouse, and that symbols mapping to more than one gene complicate the identification of roughly 18 percent of human protein-coding genes and about 10 percent of rat and mouse genes [1]. The authors recommend stable gene IDs issued by Ensembl, NCBI, HGNC, RGD, or VGNC, with genomic information retrieved directly from the authoritative database [1].

The source also matters because annotations carry errors and change between versions. Plant mitochondrial genomes are a clear case. A curated Rhodiola reference resolved 30 mitochondrial protein-coding genes with corrected exon-intron boundaries and validated 5-prime terminal variants in ccmC, ccmFn, and nad9, and the associated tools exist specifically to handle trans-splicing, RNA editing, and non-canonical start and stop codons that generic pipelines mis-annotate [2]. Whenever your annotation comes from an automated pipeline rather than a curated reference, treat the coordinates as provisional and record the release you used.

## Gene Annotation Databases in R: The Three Options at a Glance

![Decision path asking whether the species is human or mouse, whether attributes or browser tracks are needed, leading to GENCODE, Ensembl or UCSC.](/article-images/illustrations/gene-annotation-databases-r-annotation-source-triage.webp)

| Criterion | UCSC | Ensembl | GENCODE |
|--|--|--|--|
| Core product | Genome browser track data and liftOver chain files | Versioned gene database with BioMart and REST APIs | Curated GTF and GFF3 for human and mouse |
| Human builds | hg19, hg38, and hundreds of others | GRCh38 with patch releases, GRCh37 archive | GRCh38, plus a GRCh37 archive |
| Main R access | rtracklayer, makeTxDbFromUCSC, prebuilt TxDb packages, AnnotationHub | biomaRt, ensembldb, EnsDb via AnnotationHub, REST | download.file plus rtracklayer or makeTxDbFromGFF |
| Release model | Rolling track updates | Numbered releases, for example release 113 [1] | Numbered releases tied to Ensembl for human and mouse |
| Identifier style | Track-specific IDs (knownGene, refGene) | Stable ENSG and ENST IDs with version suffixes | Same Ensembl IDs plus tags such as basic and MANE_Select |
| Transcript set | Depends on the track | GENCODE gene set for human and mouse | Comprehensive and basic subsets |
| Species coverage | Very broad assembly list | Ensembl Genomes divisions for plants, bacteria, fungi, metazoa, protists | Human and mouse only |
| Best for | Browser tracks, conservation, repeats, liftover | Attribute-rich queries, orthologs, version pinning | Canonical gene models for counting and reporting |
| Main watch-out | Track contents differ from Ensembl's own release | Live service changes with each release | PAR_Y duplicates and versioned IDs |

Three relationships in that table deserve emphasis. GENCODE is the annotation that Ensembl ships for human and mouse, repackaged with extra tags and file flavors. UCSC tracks are assembled by the UCSC team from several inputs, so gene counts on the same build can differ from Ensembl's own release. None of the three is a static file, they are moving targets that you must pin down by release number.

## Retrieving UCSC Annotations in R

UCSC hosts assemblies across many clades and layers annotation tracks on top of each one. For human hg38, gene-level tracks include knownGene (UCSC's own gene set), refGene (mirroring NCBI RefSeq), ensGene (mirroring Ensembl), and GENCODE mirrors, alongside regulatory, repeat, and conservation tracks. Each track carries its own identifier namespace, so pick the track that matches your question rather than grabbing whichever one loads first.

### Track tables with rtracklayer

The rtracklayer package queries the UCSC Table Browser from R.

```r
library(rtracklayer)
session <- browserSession()
genome(session) <- "hg38"
q  <- ucscTableQuery(session, "knownGene")
kg <- getTable(q)
head(kg)
```

getTable returns a data frame of the requested table. Track identifiers occasionally change between assemblies, so confirm the exact name in the Table Browser before scripting. The Table Browser can also export BED, GTF, and other formats directly, and the UCSC REST API at api.genome.ucsc.edu handles quick lookups from any language. For the knownGene track specifically, the companion kgXref table links track IDs to gene symbols and descriptions, which is useful when you need readable labels.

### TxDb objects and AnnotationHub

The GenomicFeatures package turns UCSC tracks into TxDb objects, SQLite-backed containers that answer coordinate questions quickly.

```r
library(GenomicFeatures)
supportedUCSCtables()
txdb <- makeTxDbFromUCSC(genome = "hg38", tablename = "knownGene")
tx_by_gene <- transcriptsBy(txdb, by = "gene")
```

Prebuilt packages such as TxDb.Hsapiens.UCSC.hg38.knownGene and TxDb.Hsapiens.UCSC.hg38.refGene install the same content without a live connection at analysis time. AnnotationHub distributes these resources as well and caches whatever you download, so later sessions work offline.

One structural note. TxDb objects store coordinates and identifiers, not richer attributes such as complete symbol lists or biotypes, so keep a separate key table if you need those fields.

### LiftOver for build conversion

UCSC publishes chain files that map intervals between genome builds, a process called liftover. rtracklayer reads and applies them.

```r
chain   <- import.chain("hg38ToHg19.over.chain")
hg19_gr <- unlist(liftOver(gr_hg38, chain))
```

Chain files live under the goldenPath area of hgdownload.soe.ucsc.edu, for example hg38ToHg19.over.chain.gz for converting hg38 coordinates back to hg19. Variant-level annotation platforms maintain parallel hg38 and hg19 support for exactly this reason, because both builds remain in active circulation [3].

## Retrieving Ensembl Annotations in R

### biomaRt for attribute-rich queries

biomaRt queries the live Ensembl BioMart web service. Its strength is attributes, meaning you can pull GO terms, PubMed IDs, orthologs, variant consequence predictions, and UniProt cross-references in one call.

```r
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
annot <- getBM(
  attributes = c("ensembl_gene_id", "hgnc_symbol", "chromosome_name",
                 "start_position", "end_position", "gene_biotype"),
  filters    = "hgnc_symbol",
  values     = c("TP53", "BRCA1", "EGFR"),
  mart       = ensembl
)
```

Attribute and filter names differ between datasets and releases. listAttributes(ensembl) and listFilters(ensembl) print what is currently available. One-to-many attributes such as go_id return one row per relationship, so deduplicate after merging. Because the service reflects the current release, results can change between runs unless you pin a version.

### ensembldb and EnsDb for GenomicRanges workflows

The ensembldb package provides EnsDb objects, self-contained SQLite databases exposing Ensembl gene models through GenomicRanges-style accessors. AnnotationHub is the easiest source.

```r
library(AnnotationHub)
ah   <- AnnotationHub()
hits <- query(ah, c("EnsDb", "Homo sapiens"))
hits
edb  <- hits[[length(hits)]]    # the most recent EnsDb your cache offers

library(ensembldb)
g <- genes(edb, columns = c("gene_id", "gene_name", "gene_biotype"))
```

For keyed lookups, select() behaves like the classic AnnotationDbi interface.

```r
sel <- select(edb,
              keys    = c("ENSG00000141510", "ENSG00000012048"),
              keytype = "GENEID",
              columns = c("SYMBOL", "GENEBIOTYPE", "SEQSTRAND"))
```

EnsDb objects work offline once cached, return GRanges natively, and support the AnnotationFilter query grammar. Vertebrate releases also carry protein annotations where the underlying Ensembl data includes them. For transcript-level work, exonsBy(edb, by = "tx") gives per-transcript exon coordinates suitable for quantification pipelines.

### Archives and Ensembl Genomes

Reproducible scripts should pin a release. biomaRt accepts a version argument that routes to the matching archive server.

```r
ensembl_v106 <- useEnsembl(biomart = "genes", version = 106)
```

For species outside vertebrates, useEnsemblGenomes reaches the plants, metazoa, fungi, bacteria, and protists divisions.

```r
plants <- useEnsemblGenomes(biomart = "plants_mart",
                            dataset = "athaliana_eg_gene")
```

Release pinning is standard practice in identifier research. The nomenclature survey discussed above drew every protein-coding identifier for human, mouse, rat, and pig from a single Ensembl release, 113, precisely so the counts would be comparable [1].

## Retrieving GENCODE Annotations in R

### What GENCODE adds over plain Ensembl files

GENCODE merges manual annotation from the Havana team with Ensembl automatic annotation for human and mouse. Files come in comprehensive and basic flavors. The comprehensive set includes every plausible transcript, while the basic set restricts to a high-confidence subset suited to counting. Each feature carries tags such as basic, CCDS (the consensus coding sequence project), and MANE_Select, plus a level field grading evidence from 1 (verified loci) down to 3 (automatically annotated). MANE_Select marks the single transcript per gene agreed between NCBI RefSeq and Ensembl, which makes it a sensible default for reporting gene-level results.

Release numbering differs by species. Human releases carry plain numbers (release 44, for example), mouse releases carry an M prefix, and human GENCODE numbers historically offset the Ensembl release they accompany, so GENCODE 44 shipped alongside Ensembl 109. Record both numbers in your methods.

### Downloading and parsing the GTF

GENCODE files are hosted on EMBL-EBI infrastructure. Pick the primary assembly GTF when your reference FASTA is a plain primary assembly, since the full-genome files add scaffolds your alignment index may lack.

```r
url  <- paste0(
  "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/",
  "gencode.v44.annotation.gtf.gz"
)
dest <- file.path(tempdir(), "gencode.v44.annotation.gtf.gz")
download.file(url, dest, mode = "wb")
```

Substitute the current release number, and verify checksums for large downloads. Parsing takes one call with rtracklayer.

```r
library(rtracklayer)
gencode       <- import(dest, format = "gtf")
genes_gencode <- gencode[mcols(gencode)$type == "gene"]
table(mcols(genes_gencode)$gene_type)
```

Building a TxDb from the same file works through GenomicFeatures.

```r
library(GenomicFeatures)
txdb_gencode <- makeTxDbFromGFF(dest, format = "gtf")
```

Because the TxDb keeps coordinates and IDs but discards most attributes, park the metadata in a side table keyed by gene_id.

```r
gene_meta <- as.data.frame(mcols(genes_gencode))[
  , c("gene_id", "gene_name", "gene_type", "level")
]
saveRDS(gene_meta, "gencode_v44_gene_metadata.rds")
```

### PAR regions and versioned identifiers

Two GENCODE quirks break naive pipelines. First, pseudoautosomal region (PAR) genes, which sit in the segments of the sex chromosomes that recombine, appear twice, once on chrX and once on chrY with a _PAR_Y suffix appended to the Ensembl ID. Parsers that demand unique transcript IDs can fail on these duplicates, so filter them deliberately.

```r
par_y       <- grep("_PAR_Y", genes_gencode$gene_id, value = TRUE)
genes_main  <- genes_gencode[!grepl("_PAR_Y", genes_gencode$gene_id)]
```

Second, GENCODE writes versioned identifiers such as ENSG00000141510.18, where the suffix after the dot is the annotation version. Downstream resources often use unversioned IDs, so strip versions consistently on both sides of any join.

```r
genes_main$gene_id_unversioned <- sub("\\.\\d+$", "", genes_main$gene_id)
```

## Pros and Cons of the Three Options

Which source is the best gene annotation database in R depends on the job, and the criteria below decide most real choices.

### Coordinate systems and genome builds

UCSC names chromosomes with a chr prefix (chr1), Ensembl and GENCODE without (1). hg19 equals GRCh37 and hg38 equals GRCh38, while Ensembl adds periodic patch releases to GRCh38. Pick one build per project, record it, and convert deliberately rather than mixing sources and hoping intervals line up. The persistence of both major human builds is why variant annotation resources such as FAVOR support hg38 and hg19 side by side [3].

### Identifier stability

Symbols are convenient labels and poor keys. The release 113 survey quantified how often symbols collide across human, rat, mouse, and pig, and concluded that stable IDs from Ensembl, NCBI, HGNC, RGD, or VGNC are the dependable choice [1]. In practice, keep Ensembl gene IDs as the primary key in every table, attach symbols only for display, and convert through one pinned release.

### Transcript models

For human and mouse, Ensembl and GENCODE describe the same models, so the choice is about packaging rather than biology. UCSC tracks are a different matter. knownGene and refGene each follow their own model-building logic and will disagree with each other and with Ensembl on transcript boundaries. For counting reads, the GENCODE basic set is the common default, while the comprehensive set suits transcript-level exploration.

### Species coverage

UCSC covers the widest range of assemblies, Ensembl Genomes adds plants, fungi, bacteria, metazoa, and protists, and GENCODE covers human and mouse only. Outside those, the realistic paths are building a custom annotation package from NCBI data [4], using community platforms such as ShEnrich for penaeid shrimp [5], or taking the GFF3 shipped alongside a newly published genome [6].

### Reproducibility

Live services change under running scripts. Pin Ensembl releases with the version argument or archive hosts, snapshot GENCODE GTF files locally, and rely on AnnotationHub caching for offline work. Save parsed annotation objects as RDS files next to your count matrix, and record the release in every methods paragraph and figure caption.

## Worked Scenarios

### Scenario 1: Attach symbols and biotypes to an RNA-seq result table

Typical setup. Counts were generated against a GENCODE GTF, differential expression ran in limma, DESeq2, or edgeR, and the output table needs symbols and biotypes. Strip versions, map through one EnsDb, and merge.

```r
library(ensembldb)

res_df <- as.data.frame(de_results)
res_df$gene_id <- sub("\\.\\d+$", "", rownames(res_df))

mapped <- select(edb,
                 keys    = res_df$gene_id,
                 keytype = "GENEID",
                 columns = c("SYMBOL", "GENEBIOTYPE"))

res_annot <- merge(mapped, res_df,
                   by.x = "GENEID", by.y = "gene_id", all.y = TRUE)
```

select() can return multiple rows per key when identifiers map to several entries, so deduplicate before merging. Downstream enrichment inherits whatever identifier discipline you applied here. A published atherosclerosis analysis retrieved four GEO expression datasets (GSE100927, GSE43292, GSE28829, and GSE20129), ran limma for differential expression, then applied GO and KEGG annotation and enrichment to the resulting gene lists [7]. Clean mapping upstream is what makes that enrichment interpretable.

### Scenario 2: Annotate ChIP-seq peaks with ChIPseeker

ChIP-seq (chromatin immunoprecipitation sequencing) produces peaks, genomic intervals with enriched signal, and ChIPseeker annotates them relative to the nearest TSS (transcription start site), classifies peak locations, and compares profiles across datasets [8]. It accepts any TxDb, so both the UCSC and GENCODE routes above feed it.

```r
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

anno <- annotatePeak(
  "peaks.bed",
  tssRegion = c(-2000, 2000),
  TxDb      = TxDb.Hsapiens.UCSC.hg38.knownGene,
  annoDb    = "org.Hs.eg.db"
)
plotAnnoPie(anno)
as.data.frame(anno)
```

The package also indexes around 15,000 bed files from GEO, NCBI's Gene Expression Omnibus, letting you test your peaks for significant overlap with published datasets to infer co-regulation or transcription factor complexes [8]. ChIPseeker is released under the Artistic-2.0 license through Bioconductor [8].

### Scenario 3: Move GWAS lead variants to genes

Post-GWAS (genome-wide association study) pipelines lean on annotation constantly. A lead variant is the most strongly associated variant at a locus, and assigning it to a gene requires coordinates, gene models, and functional context. One cardiac conduction block study combined GenomicSEM, fine mapping with SuSiE and FINEMAP, and used the FUMA platform to annotate lead SNPs and candidate genes drawn from the GWAS Catalog and FinnGen summary statistics [9]. At the variant level, FAVOR provides whole-genome functional annotations through an API and interactive browser [3]. Functional priors also feed scoring methods directly. In one ischemic stroke comparison, SBayesRC, which integrates functional annotation priors, showed the largest effect size and liability-scale R2 among European ancestry participants in the REGARDS holdout validation cohort, with an odds ratio per standard deviation of 1.24 [10].

The R-side lesson is coordinate hygiene. Lead variants arrive as rsIDs or chr-position strings on a stated build. Confirm the build, lift intervals if needed, and intersect with gene bodies or promoter windows from a pinned annotation.

### Scenario 4: Non-model organisms

When no Ensembl division covers your species, build the annotation layer yourself. One documented approach constructs custom annotation R packages from NCBI genomic data, imports them as ordinary libraries, and keeps them compatible with clusterProfiler, demonstrated for Mycobacterium tuberculosis H37Rv with a rifampicin-response RNA-seq dataset [4]. Community databases fill other gaps. ShEnrich integrates 61,287 annotated proteins representing 34,169 genes across five penaeid shrimp species, with 442 KEGG pathways and 12,897 GO terms served through clusterProfiler-based enrichment [5]. For freshly sequenced species, the genome paper's own GFF3 is often the best annotation available, as with the Rosa bracteata assembly of 540.63 Mb carrying 42,789 predicted protein-coding genes, 38,394 of them functionally annotated [6].

## Common Mistakes and Limitations

- **Mixing genome builds.** hg19 coordinates intersected with hg38 gene bodies produce silent nonsense. Confirm the build of every input file, then convert once with UCSC chain files and verify a few known loci.
- **Chromosome naming mismatches.** chr1 versus 1 breaks range operations. Fix with seqlevelsStyle(gr) <- "UCSC" from the GenomeInfoDb package before intersecting.
- **Versioned IDs on one side of a join only.** ENSG00000141510.18 will not match ENSG00000141510. Strip versions on both sides with the same regular expression.
- **Symbol-first joins.** Symbols collide across genes, which the Ensembl release 113 survey documented at scale [1]. Join on stable IDs and add symbols at the end for readability.
- **PAR_Y duplicates inflating chrY counts.** Deduplicate pseudoautosomal entries before counting or plotting sex chromosome coverage.
- **Assuming UCSC knownGene equals the Ensembl gene set.** The tracks have different building rules and different ID namespaces. Match the track to the question.
- **Unpinned releases.** Annotation drifts between releases, and a result that cannot name its release cannot be reproduced. Record Ensembl, GENCODE, and UCSC track versions explicitly.
- **TxDb attribute loss.** Building a TxDb discards most GTF attributes. Keep the metadata side table you saved as RDS.
- **Treating annotation as ground truth.** Curated resources exist because automated pipelines mis-annotate difficult biology, including the trans-spliced and RNA-edited transcripts that motivated dedicated mitochondrial annotation tools [2].

A final limitation belongs in the open. Annotation quality varies by species, by track, and by release, coordinates shift with patches, and any claim headed for publication or clinical interpretation deserves review by the person accountable for the analysis. Checking the release notes of the annotation you ship is part of the analysis, not an optional extra.

## The AnnotationHub Ecosystem: Querying, Caching, and Offline Work

AnnotationHub is the distribution layer that most Bioconductor annotation travels through, and treating it as more than a download button pays off. Every resource in the hub carries metadata you can inspect before committing to a download: species, genome build, source database, release, and the date the record was added. That metadata is the fastest way to answer the question that sinks reproducibility most often, namely which release of an annotation you actually have.

```r
library(AnnotationHub)
ah <- AnnotationHub()

# Inspect what you are about to download
hits <- query(ah, c("EnsDb", "Homo sapiens"))
mcols(hits)[, c("title", "dataprovider", "species",
                "rdataclass", "sourceversion", "pretryview")]

# Record the snapshot date alongside your project notes
snapshotDate(ah)
```

The `snapshotDate` matters because the hub is versioned itself. When a hub snapshot updates, new resources appear and old ones can be reorganized, so two labs running the same script a month apart can retrieve different EnsDb builds. Writing `snapshotDate(ah)` into a project log, alongside the hub accession of the resource you took, turns a vague "downloaded from AnnotationHub" methods statement into something a reviewer can trace.

Once a resource is downloaded, the hub caches it locally. Later sessions load from the cache, which means your analysis keeps working when the network does not. Two habits make this reliable:

1. **Do not index the cache by position alone.** `hits[]` works today and breaks when the query result reorders. Filter on a stable field such as the title or genome build, then subset.
2. **Copy cached resources into your project directory.** The cache is a shared user-level store that can be pruned; an RDS or SQLite copy inside the project folder survives cache resets.

AnnotationHub also distributes resources beyond EnsDb, including UCSC-derived TxDb packages, OrgDb objects, and sequence data, so a single query interface covers most of the annotation stack. When you cannot find exactly the resource you want, for example a GENCODE-flavored EnsDb at your preferred release, you can build one locally rather than waiting for the hub to ship it.

## Building Your Own Containers from GENCODE Files

The article's main walkthrough shows `makeTxDbFromGFF` on a GENCODE GTF. Two refinements make that step production-grade rather than exploratory.

First, supply chromosome metadata explicitly. `makeTxDbFromGFF` inspects the sequences it finds and will warn about, and in some cases drop, contigs it cannot classify, which produces a TxDb with silently missing sequences. Passing a `chrominfo` object built from the same assembly removes the guesswork:

```r
library(GenomeInfoDb)
chrominfo <- getChromInfoFromUCSC("hg38")
# Align the naming style to the GTF (Ensembl-style, no chr prefix)
chrominfo$circles <- chrominfo$circRNA  # keep UCSC circular- chromosome flags
txdb_gencode <- makeTxDbFromGFF(
  dest,
  format    = "gtf",
  chrominfo = chrominfo
)
```

If you use this route, verify with `seqlevels(txdb_gencode)` that the sequence names match your alignment index. A TxDb that carries sequences absent from your BAM files will generate false-positive "unannotated" regions during counting.

Second, know the limits of `ensDbFromGtf`. The ensembldb package can construct EnsDb objects directly from Ensembl-formatted GTF files, which is convenient for building versioned, offline containers. GENCODE GTF files differ in attribute layout from Ensembl's own exports, and that difference can trip the parser. When the goal is a GENCODE-faithful container, the safer pattern is the TxDb route plus the metadata side table described earlier. When the goal is an Ensembl-release-faithful container, download the Ensembl GTF for that release and build the EnsDb from it directly. Matching the file flavor to the container builder avoids an entire class of parsing bugs.

A third pattern worth knowing: you do not always need a container at all. For many projects, a filtered GRanges plus a metadata data frame is sufficient, faster to build, and easier to version:

```r
gencode        <- import(dest, format = "gtf")
exons_gencode  <- gencode[mcols(gencode)$type == "exon"]
tx_gencode     <- gencode[mcols(gencode)$type == "transcript"]
saveRDS(list(genes = genes_gencode, exons = exons_gencode,
             tx = tx_gencode, meta = gene_meta),
        "gencode_v44_annotation.rds")
```

The trade-off is that GRanges objects live in memory, and a comprehensive human GTF is large. The performance section below covers ways to keep that manageable.

## Choosing Between TxDb, EnsDb, OrgDb, and Raw GTF Objects

Newcomers often treat these four containers as interchangeable, then discover mid-project that the one they picked cannot answer their question. Each stores a different slice of annotation:

- **TxDb** stores coordinates and identifiers: genes, transcripts, exons, CDS, keyed by internal IDs. It does not store symbols, biotypes, or cross-references. It is the right container when the question is spatial: which peaks fall in promoters, which exons overlap my reads.
- **EnsDb** stores coordinates plus Ensembl attributes: gene symbols, biotypes, and, where the source release includes them, protein-level mappings. It answers spatial questions and identifier questions in one object, which is why it is the default for human and mouse RNA-seq annotation.
- **OrgDb** stores identifier relationships and functional category membership, not coordinates. It is the engine behind GO and KEGG enrichment and symbol-to-Entrez conversion. Asking an OrgDb for a gene's chromosomal location returns nothing useful.
- **Raw GTF as GRanges** stores everything the file contains, at the cost of memory and the absence of indexed queries. It is the most faithful representation and the only one that preserves GENCODE tags such as `basic` and `MANE_Select` without an extra step.

You can interrogate any of them before committing:

```r
keytypes(edb)          # what can I use as a lookup key
columns(edb)           # what can I pull back
head(keys(edb, keytype = "GENEBIOTYPE"))  # sample key values
```

A useful discipline is to declare, at the top of your analysis script, which container answers which sub-question. For example: coordinate overlap uses the TxDb built from GENCODE 44, identifier mapping uses EnsDb at the matching release, and enrichment uses org.Hs.eg.db keyed through Entrez IDs obtained from the EnsDb. Writing this down exposes mismatches early, such as the subtle failure mode where you map symbols to Entrez via one release and coordinates via another.

## Identifier Conversion in Depth

Identifier conversion looks trivial and quietly destroys more results than any other annotation step. The release 113 nomenclature survey gives the reason: symbols collide, synonyms proliferate, and the authors' recommendation to key on stable identifiers from Ensembl, NCBI, HGNC, RGD, or VGNC exists precisely because display names cannot serve as primary keys [1]. Here is how to convert identifiers in R without losing genes.

### The three conversion tools and when each applies

`select()` from AnnotationDbi-style interfaces, including ensembldb, is the workhorse for keyed lookups. `mapIds()` is the same engine with a one-value-per-key contract, which is often what you want for building a clean lookup vector:

```r
sym <- mapIds(edb,
              keys      = res_df$gene_id,
              keytype   = "GENEID",
              column    = "SYMBOL",
              multiVals = "first")
```

`clusterProfiler::bitr()` wraps the same idea with explicit reporting of unmapped keys, and it plugs directly into the enrichment machinery that downstream analyses typically use [7]:

```r
library(clusterProfiler)
converted <- bitr(res_df$gene_id,
                  fromType = "ENSEMBL",
                  toType   = c("SYMBOL", "ENTREZID"),
                  OrgDb    = org.Hs.eg.db)
n_input <- length(unique(res_df$gene_id))
n_map   <- length(unique(converted$ENSEMBL))
message(sprintf("Mapped %d of %d genes (%.1f%%)",
                n_map, n_input, 100 * n_map / n_input))
```

That mapping percentage belongs in your results. If 15 percent of genes fail to map, every enrichment analysis you run afterward is missing 15 percent of its input, and reviewers who ask will want the number [7].

### Why keys fail, and what each failure means

Failures cluster into four causes, each with a different remedy:

- **Version suffixes.** Covered earlier; strip on both sides with one shared helper function so the regex cannot drift between scripts.
- **Retired identifiers.** Genes merge, split, and get renamed between releases. An ID absent from your pinned release may exist in a later one. Check the current release before assuming a dead key.
- **Out-of-scope IDs.** A table keyed partly on pseudogene IDs, miRNA IDs, or scaffold contigs will not map through a container filtered to a subset of features. Count these separately from true failures.
- **Symbol-based input.** If your table arrives keyed by symbols rather than IDs, collisions mean some rows are ambiguous at the outset. The survey found symbols mapping to more than one gene complicate identification for roughly 18 percent of human protein-coding genes [1], so flag rather than silently pick the first match.

### The alias trap

OrgDb objects expose an `ALIAS` keytype, and it is tempting to use it as a symbol resolver. Be careful: aliases include outdated names, names of other things, and names that have migrated between genes over the years. When a symbol-based conversion returns a surprising gene, cross-check against a second source before believing it. The disciplined pattern is: convert symbols to stable IDs once, at project start, with a human eyeballing the low-confidence rows, and never let raw symbols near a merge again.

## Sanity Checks Before You Trust an Annotation

Treat an annotation like any other dataset: it gets validation before use. The following checks take minutes and catch the errors that otherwise surface as inexplicable results weeks later.

### Count and biotype profile

```r
gene_tbl <- as.data.frame(table(mcols(genes_gencode)$gene_type))
print(gene_tbl)
```

You are not checking against a memorized number; you are checking that the profile looks like what you expect for your species and release, and recording it. A GTF that accidentally contains both comprehensive and basic subsets, or a file from the wrong release, shows up here as an implausible biotype distribution.

### Identifier overlap with your count matrix

```r
count_ids <- sub("\\.\\d+$", "", rownames(counts))
annot_ids <- sub("\\.\\d+$", "", genes_gencode$gene_id)
cat("Counts in annotation:", mean(count_ids %in% annot_ids), "\n")
```

If the overlap is below what you expect for your alignment reference, the likeliest cause is that counting ran against a different GTF than the one you just loaded. This single check prevents the most embarrassing category of annotation bug: annotating a count matrix with gene models from a different file than the one used for alignment.

### Known-locus spot check

Pick five to ten genes you know well, from different chromosomes and biotypes, and confirm they land on the expected chromosomes and strands, sit inside plausible ranges, and carry the expected symbols. This is the genomic equivalent of a plot-your-data-first habit.

### LiftOver reconciliation

![Worked calculation of the percentage of ranges kept after LiftOver, using before and after range counts.](/article-images/illustrations/gene-annotation-databases-r-liftover-retention.webp)

After conversion, quantify what survived:

```r
before <- length(gr_hg38)
after  <- length(hg19_gr)
cat(sprintf("LiftOver retained %d of %d ranges (%.1f%%)\n",
            after, before, 100 * after / before))
```

LiftOver is not lossless. Regions without orthologous sequence in the target build are dropped, and regions that moved between chromosomes in the interval between builds can land in unexpected places. Know what was discarded before interpreting an analysis on the converted coordinates. This matters practically whenever older variant-level resources on hg19 must be joined to newer annotation on hg38, a situation that persists because both builds remain in active circulation [3].

### Sex chromosome and PAR accounting

```r
table(seqnames(genes_gencode) == "chrY",
      grepl("_PAR_Y", genes_gencode$gene_id))
```

If your pipeline deduplicated PAR entries earlier, confirm the counts moved in the direction you intended.

## Comparing Gene Models Across Sources in R

Because UCSC tracks, Ensembl releases, and GENCODE files describe the same genome with different model-building rules, a direct comparison in R is the honest way to see how much your results depend on the source. The comparison also makes the abstract claim that UCSC tracks follow their own logic concrete and quantifiable.

Start by pulling two sources for the same build:

```r
library(GenomicFeatures)
library(rtracklayer)

txdb_ucsc <- makeTxDbFromUCSC(genome = "hg38", tablename = "refGene")
gencode   <- import(dest, format = "gtf")
genes_gc  <- gencode[mcols(gencode)$type == "gene"]

genes_ucsc <- genes(txdb_ucsc)
```

Then compare three quantities that matter for downstream work:

```r
# 1. How many genes does each source call on this build?
length(genes_ucsc)
length(genes_gc)

# 2. How many exons per transcript, distributionally?
exons_ucsc <- exonsBy(txdb_ucsc, by = "tx")
exons_gc   <- exonsBy(gencode[mcols(gencode)$type == "exon"],
                      by = "gene")  # illustrative grouping
summary(width(exons_ucsc))

# 3. How far apart are the two sources' TSS calls for the same gene?
#    (requires a shared key; versioned IDs need stripping first)
```

The third comparison is the most informative. Build per-gene transcript start positions from both sources, match on the unversioned Ensembl ID where the source provides one, and compute distances with `GenomicRanges::distanceToNearest()`. A histogram of those distances shows you, rather than tells you, that transcript boundaries differ between sources, and the heavy tail of large distances identifies exactly the genes where your counting or peak-annotation results will be sensitive to the source choice.

For a targeted comparison around a specific question, restrict to a promoter window and ask how many transcripts each source places there:

```r
prom <- promoters(genes_gc, upstream = 2000, downstream = 200)
ov   <- countOverlaps(prom, exons_ucsc)
summary(ov)
```

Genes where sources disagree strongly are candidates for a per-gene curation decision, and the MANE_Select tag in GENCODE exists precisely to give you a defensible single transcript per gene when such decisions are needed. For projects headed toward publication, running the core analysis under two annotation sources and reporting where conclusions agree is inexpensive insurance, and it converts the "which database did you use" reviewer question from a vulnerability into a documented robustness check.

## Performance, Memory, and Network Resilience

Annotation steps are where long analyses stall, and the fixes are mostly structural.

**Parse big GTF files selectively.** rtracklayer's `import()` accepts a `which` argument taking a GRanges, so you can pull only the region you need rather than the whole genome:

```r
region <- GRanges("chr17", IRanges(1, 8e6))  # example: one chromosome arm
gtf_sub <- import(dest, format = "gtf", which = region)
```

For whole-file parsing where you only want certain feature types, reading the compressed file through a filtered connection and subsetting lines before full parsing keeps memory bounded. The gene and transcript rows are a small fraction of a comprehensive GTF; exon rows dominate.

**Let SQLite do the work.** TxDb and EnsDb objects are backed by SQLite, so range queries and key lookups run against disk-backed indexes rather than in-memory data frames. This is why a 30-gene lookup against an EnsDb completes instantly while the equivalent operation on a parsed GTF data frame can crawl. For repeated per-gene queries in a loop, batch them into a single `select()` call rather than looping.

**Chunk large biomaRt queries.** The BioMart web service performs best with focused queries. Pulling dozens of attributes for tens of thousands of genes in one call invites timeouts, and one-to-many attributes such as GO terms multiply the result rows. Split the work: coordinates and symbols in one query, GO terms in a second, orthologs in a third, then join on the stable ID. Keep the join keyed on identifiers, never on row order.

**Expect and plan for service interruptions.** Live services, whether BioMart or the UCSC Table Browser, fail occasionally and change without notice. The structural defense is the one already emphasized: cache everything, snapshot parsed objects as RDS, and write your pipeline so the network-touching annotation step runs once, at project setup, not inside every analysis rerun. If a track name disappears between assemblies, `supportedUCSCtables()` and the Table Browser itself are the first places to look for the replacement.

## A Troubleshooting Field Guide

These are the failure patterns that recur across projects, with the first thing to try for each.

**biomaRt returns an empty data frame.** Nine times out of ten the attribute or filter name does not exist in the current release, and the service returns an empty result rather than an error. Run `listAttributes()` and `listFilters()` against your mart object and confirm every name character for character. The remaining one time in ten, the values themselves are wrong: versioned IDs passed to a filter expecting unversioned IDs silently match nothing.

**makeTxDbFromUCSC fails or produces an odd TxDb.** UCSC track schemas evolve, and the helper sometimes lags. Check the exact table name in the Table Browser, try `supportedUCSCtables()` for the current mapping, and consider whether a prebuilt TxDb package or a GENCODE GTF covers the need instead. A failing live fetch is often a signal that the pinned-local-file approach fits your project better anyway.

**LiftOver silently drops intervals.** Not an error, a property. Compare lengths before and after as shown above, and inspect a sample of dropped intervals. If the drops concentrate in one chromosome or region, suspect an assembly gap or a sequence renamed between builds rather than a bug in your code.

**GRanges operations throw seqlevels mismatch errors.** Two objects with different chromosome naming or different seqlevel sets cannot be intersected directly. Harmonize with `GenomeInfoDb` tools: `seqlevelsStyle()` to convert naming, then `keepSeqlevels()` or `dropSeqlevels()` to make the two objects' sequence sets agree before the operation.

**EnsDb keys return NA rows.** Check the keytype. `select()` with `GENEID` expects Ensembl gene IDs, and passing symbols with `keytype = "SYMBOL"` is a different query that may legitimately return nothing for IDs that were valid gene IDs. Also confirm you did not strip versions on only one side of the workflow.

**AnnotationHub query returns nothing.** Loosen the query. Multi-term queries are conjunctive, so `query(ah, c("EnsDb", "Homo sapiens", "103"))` fails if the release string is formatted differently in the metadata. Search with fewer terms, inspect `mcols(hits)`, then narrow.

**GENCODE ID joins match at near zero percent.** Almost always the version-suffix problem in one direction, or a PAR_Y-suffixed ID set that never matches. Apply the shared version-stripping helper and recheck.

## Common Misconceptions, Corrected

**"The GENCODE release number matches the Ensembl release number."** It does not, and the offset is historical rather than accidental. GENCODE 44 shipped alongside Ensembl 109, so recording only one number leaves the other ambiguous. Record both, every time.

**"refGene on UCSC is just Ensembl data with another name."** NCBI RefSeq and Ensembl maintain separate gene model sets with separate curation processes. The tracks will disagree on transcript boundaries, and treating them as interchangeable produces results that depend on which track happened to be loaded.

**"Liftover is a faithful coordinate translation."** It is a best-effort mapping with real losses. Build-level conversions exist because both hg19 and hg38 remain in active use across the variant annotation ecosystem [3], but every converted interval deserves the reconciliation check described above.

**"OrgDb packages contain gene coordinates."** They do not. OrgDb objects hold identifier cross-references and functional category membership. Coordinates come from TxDb, EnsDb, or GTF-derived GRanges. Confusing the two leads to futile searches through enrichment databases for chromosomal locations.

**"Everything in the comprehensive GENCODE set is a confirmed transcript."** The comprehensive set deliberately includes every plausible model, and the `level` field grades evidence down to level 3, annotated automatically. The basic subset, and tags such as `MANE_Select` and `CCDS`, exist so you can restrict to higher-confidence models when the analysis calls for it.

**"AnnotationHub will always give me the latest release."** The hub serves what its snapshot contains. That behavior is a feature for reproducibility, but only if you record the snapshot date and the resource you retrieved, as described earlier.

**"Curated annotation means correct annotation."** Curation improves annotation; it does not make it final. The plant mitochondrial literature makes the point vividly: dedicated tools and curated references exist because trans-splicing, RNA editing, and non-canonical start and stop codons defeat generic pipelines, and a curated Rhodiola reference was needed to correct exon-intron boundaries and validate transcript termini that automated annotation had wrong [2]. Difficulty is not confined to plant mitochondria; it lives wherever genome biology is unusual, and unusual is a relative term.

## A Project Layout and Methods Paragraph for Reproducible Annotation

Annotation reproducibility is mostly file management plus prose discipline. A layout that has aged well:

```text
project/
  annotation/
    gencode.v44.annotation.gtf.gz
    gencode.v44.annotation.gtf.gz.md5
    gencode_v44_gene_metadata.rds
    txdb_gencode_v44.sqlite
    annotation_log.txt          # snapshotDate, hub IDs, fetch dates
  results/
    de_results_annotated.rds
  scripts/
    01_build_annotation.R
    02_annotate_results.R
  sessionInfo.txt
```

The annotation log is a plain text file with one line per resource: what it is, which release, where it came from, and when. Two minutes of bookkeeping at fetch time replaces an afternoon of archaeology when a collaborator asks which GTF produced the count matrix.

The companion script, `01_build_annotation.R`, performs every network operation: downloads the GTF, verifies the checksum, builds containers, runs the sanity checks, and saves the RDS and SQLite artifacts. Everything downstream loads from those artifacts and never touches the network. Rerunning the analysis, on a train or on a compute cluster without internet access, then works identically.

For the methods paragraph, a template that covers what reviewers actually ask:

```text
Gene annotations were obtained from GENCODE release 44 (human, GRCh38,
primary assembly GTF). Gene models were parsed with rtracklayer X.Y.Z
and stored as a TxDb via GenomicFeatures A.B. Identifier mapping used
EnsDb built from Ensembl release [matching release] via ensembldb C.D.
Symbols were attached for display only; all joins used unversioned
Ensembl gene IDs. The mapping rate for the differential expression
table was NN percent. LiftOver conversions, where applied, used
UCSC chain files and retained NN percent of intervals.
```

Every bracketed value comes from your annotation log. If a value cannot be filled in, the pipeline is missing a record, which is exactly the failure this layout exists to prevent.

## Applications Beyond the RNA-seq Default

The scenarios in the main article cover RNA-seq, ChIP-seq, GWAS, and non-model genomes. Three further applications show how the same annotation stack flexes.

**Promoter panel design.** Designing a targeted panel, whether for methylation assays, capture sequencing, or qPCR arrays, is an annotation retrieval problem: take a gene list, extract promoter windows from a pinned annotation, add strand awareness so the window tracks the correct TSS, and export as BED. The `promoters()` accessor on an EnsDb or TxDb does the heavy lifting, and the sanity-check habits above apply directly, since a panel built against the wrong build is expensive to discover at the bench.

**Cross-species work.** Orthology questions route through biomaRt's ortholog attributes or through comparative resources, and the identifier discipline matters doubly because symbols rarely correspond one-to-one across species. The release 113 survey that quantified symbol ambiguity in human, rat, mouse, and pig [1] is a reminder that cross-species symbol matching, a pattern still seen in older pipelines, should be replaced by stable-ID orthology mapping with a pinned release on each side.

**Community-scale annotation for emerging organisms.** When a species sits outside the major databases, the realistic workflow combines whatever exists. The penaeid shrimp resource ShEnrich integrates tens of thousands of annotated proteins and hundreds of KEGG pathways served through clusterProfiler-based enrichment [5], and custom NCBI-derived annotation packages keep results compatible with the same enrichment machinery [4]. For a freshly assembled genome, the publication's own GFF3 is often the best annotation available, as with the Rosa bracteata assembly carrying 42,789 predicted protein-coding genes, of which 38,394 were functionally annotated [6]. In each case the R-side pattern is identical to the model-organism workflow: parse once, containerize, sanity-check, pin, and record.

**Post-GWAS gene prioritization.** The cardiac conduction block study cited earlier combined fine mapping with platform-based annotation of lead SNPs and candidate genes [9], and functional priors flow directly into modern scoring methods, as the ischemic stroke comparison of SBayesRC against alternative approaches demonstrates [10]. The R contribution to such pipelines is the coordinate-hygiene layer: confirming the build of every summary statistic file, lifting where required, and intersecting lead variants with gene bodies and promoter windows from one pinned annotation so that fine-mapping credibility intervals and gene models speak the same coordinate language.

Across all of these, the through-line is the same. Choose the source that matches the question, containerize it once, validate it before use, and write down exactly what you used. The code in this article is the mechanical part; the habits are what make results hold up.

## Frequently Asked Questions

### Which database should I use for a human RNA-seq project?

Use a GENCODE GTF for counting and an Ensembl EnsDb or biomaRt query for identifier mapping, pinned to the same release. GENCODE defines the gene models, Ensembl supplies the cross-references, and matching releases keep coordinates and IDs consistent end to end.

### How do I convert Ensembl gene IDs to gene symbols in R?

Strip version suffixes first, then use ensembldb's select() with keytype "GENEID", biomaRt's getBM with the ensembl_gene_id filter, or clusterProfiler's bitr with an OrgDb. Always report how many keys failed to map so nothing disappears silently.

### Are GENCODE and Ensembl annotations the same thing?

For human and mouse, Ensembl ships GENCODE models, so the underlying biology matches. The differences are packaging, since GENCODE adds tags such as basic and MANE_Select, an evidence level field, and comprehensive versus basic file flavors.

### How do I handle hg19 versus hg38 coordinate mismatches?

Confirm the build of every input, choose one build for the project, and convert with UCSC chain files through rtracklayer's liftOver. Spot-check several well-known loci after conversion before trusting the rest.

### Can I get these annotations for non-model organisms?

UCSC and the Ensembl Genomes divisions cover many species, while GENCODE covers human and mouse only. For uncovered species, build a custom annotation package from NCBI data or adopt a community resource, and expect to curate manually.

### How do I make my annotation step reproducible?

Record release numbers everywhere, use archived Ensembl marts or local AnnotationHub caching, save parsed annotation objects as RDS files beside your results, and capture sessionInfo output with the analysis.

<script type="application/ld+json">
{
  "@context": "https://schema.org",
  "@type": "FAQPage",
  "mainEntity": [
    {
      "@type": "Question",
      "name": "Which database should I use for a human RNA-seq project?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Use a GENCODE GTF for counting and an Ensembl EnsDb or biomaRt query for identifier mapping, pinned to the same release. GENCODE defines the gene models, Ensembl supplies the cross-references, and matching releases keep coordinates and IDs consistent end to end."
      }
    },
    {
      "@type": "Question",
      "name": "How do I convert Ensembl gene IDs to gene symbols in R?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Strip version suffixes first, then use ensembldb's select() with keytype GENEID, biomaRt's getBM with the ensembl_gene_id filter, or clusterProfiler's bitr with an OrgDb. Always report how many keys failed to map so nothing disappears silently."
      }
    },
    {
      "@type": "Question",
      "name": "Are GENCODE and Ensembl annotations the same thing?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "For human and mouse, Ensembl ships GENCODE models, so the underlying biology matches. The differences are packaging, since GENCODE adds tags such as basic and MANE_Select, an evidence level field, and comprehensive versus basic file flavors."
      }
    },
    {
      "@type": "Question",
      "name": "How do I handle hg19 versus hg38 coordinate mismatches?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Confirm the build of every input, choose one build for the project, and convert with UCSC chain files through rtracklayer's liftOver. Spot-check several well-known loci after conversion before trusting the rest."
      }
    },
    {
      "@type": "Question",
      "name": "Can I get these annotations for non-model organisms?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "UCSC and the Ensembl Genomes divisions cover many species, while GENCODE covers human and mouse only. For uncovered species, build a custom annotation package from NCBI data or adopt a community resource, and expect to curate manually."
      }
    },
    {
      "@type": "Question",
      "name": "How do I make my annotation step reproducible?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Record release numbers everywhere, use archived Ensembl marts or local AnnotationHub caching, save parsed annotation objects as RDS files beside your results, and capture sessionInfo output with the analysis."
      }
    }
  ]
}
</script>

## Sources

1. [Protein-coding genes in humans and model mammals (mouse, rat and pig): gene identifiers and disambiguation of gene nomenclature retrieved from the Ensembl genome browser.](https://pubmed.ncbi.nlm.nih.gov/41408139/)
2. [RhoMitoAnnotator and Polypods, Bioinformatics Tools for the Rhodiola Mitochondrial Gene Assembly, Annotation and Phylogenetic Analysis.](https://pubmed.ncbi.nlm.nih.gov/42196420/)
3. [FAVOR 2.0: A reengineered functional annotation of variants online resource for interpreting genomic variation.](https://pubmed.ncbi.nlm.nih.gov/41335103/)
4. [Enhanced RNA-Seq Expression Profiling and Functional Enrichment in Non-model Organisms Using Custom Annotations.](https://pubmed.ncbi.nlm.nih.gov/42367215/)
5. [ShEnrich: a functional annotation database and enrichment platform for penaeid shrimp.](https://pubmed.ncbi.nlm.nih.gov/42963765/)
6. [The chromosome-scale genome assembly and annotation of Rosa bracteata (Macartney Rose).](https://pubmed.ncbi.nlm.nih.gov/41794947/)
7. [CASP1 as a key autophagy-related gene in atherosclerosis: Identification by bioinformatics analysis and experimental validation.](https://pubmed.ncbi.nlm.nih.gov/42616740/)
8. [ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization.](https://pubmed.ncbi.nlm.nih.gov/25765347/)
9. [Genomic Structural Equation Modeling Combined With Post-GWAS Analysis Identifies Two Risk Gene Loci and Functionally Sensitive Genes Associated With Cardiac Conduction Block.](https://pubmed.ncbi.nlm.nih.gov/41541644/)
10. [Evaluation of functional annotation-informed and ancestry-specific polygenic risk scores for ischemic stroke.](https://pubmed.ncbi.nlm.nih.gov/42787420/)