Anotación de genes: AnnotationDbi, org.*.eg.db

r
bioconductor
rnaseq
Cómo pasar de IDs (ENSG00000123456) a símbolos (TP53) y descripciones. Los paquetes org.*.eg.db, AnnotationDbi::mapIds, las trampas de mapeos uno-a-muchos y una pasada rápida a clusterProfiler para enriquecimiento funcional.

El problema: IDs vs símbolos

Tras DESeq2, tu tabla de resultados tiene IDs como:

ENSG00000141510   (Ensembl)
NM_000546.6       (RefSeq)
7157              (Entrez)

Para reportar, necesitas el nombre humano: TP53. Y a veces más: descripción, función, longitud, cromosoma.

Bioconductor lo resuelve con los paquetes org.*.eg.db, bases de datos de anotación por organismo, y la API genérica de AnnotationDbi para consultarlas.

Los paquetes org.*.eg.db

Cada organismo tiene su paquete:

Paquete Organismo
org.Hs.eg.db Homo sapiens
org.Mm.eg.db Mus musculus
org.Rn.eg.db Rattus norvegicus
org.Dm.eg.db Drosophila melanogaster
org.Sc.sgd.db Saccharomyces cerevisiae
org.At.tair.db Arabidopsis thaliana

Instalación estándar:

BiocManager::install("org.Hs.eg.db")

Cada uno contiene mapeos entre muchos tipos de identificador y enriquecimientos (GO terms, KEGG paths). Suelen pesar 100-200 MB cada uno. Pesados al descargar, ligeros en uso.

mapIds: la función central

library(AnnotationDbi)
library(org.Hs.eg.db)

ids <- c("ENSG00000141510", "ENSG00000012048", "ENSG00000139618")

simbolos <- mapIds(org.Hs.eg.db,
                   keys      = ids,
                   column    = "SYMBOL",
                   keytype   = "ENSEMBL",
                   multiVals = "first")

simbolos
# ENSG00000141510 ENSG00000012048 ENSG00000139618
#         "TP53"        "BRCA1"        "BRCA2"

Argumentos:

  • keys: los IDs que tienes.
  • column: qué quieres recuperar (SYMBOL, GENENAME, ENTREZID, …).
  • keytype: qué son tus IDs (ENSEMBL, REFSEQ, ENTREZID, UNIPROT, …).
  • multiVals: cómo manejar mapeos uno-a-muchos: "first", "list", "CharacterList".

Para saber qué columnas tiene un paquete:

columns(org.Hs.eg.db)
# [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID"
#     "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL" "IPI" "MAP"
#     "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ"
#     "SYMBOL" "UCSCKG" "UNIGENE" "UNIPROT"

Y los keytypes válidos:

keytypes(org.Hs.eg.db)

Suelen ser los mismos de columns, pero no todos son válidos como punto de partida.

Añadir columnas al DESeqResults

El patrón idiomático:

library(dplyr)
library(tibble)

res_anotado <- as.data.frame(res) |>
    rownames_to_column("ensembl_id") |>
    mutate(
        simbolo = mapIds(org.Hs.eg.db, keys = ensembl_id,
                         column = "SYMBOL", keytype = "ENSEMBL",
                         multiVals = "first"),
        descripcion = mapIds(org.Hs.eg.db, keys = ensembl_id,
                             column = "GENENAME", keytype = "ENSEMBL",
                             multiVals = "first"),
        entrez_id = mapIds(org.Hs.eg.db, keys = ensembl_id,
                           column = "ENTREZID", keytype = "ENSEMBL",
                           multiVals = "first")
    ) |>
    arrange(padj)

head(res_anotado, 5) |> select(ensembl_id, simbolo, descripcion, log2FoldChange, padj)

Tres mapeos en un solo chain. Símbolos para reporte, descripción para readability, ENTREZID para enriquecimiento funcional.

La trampa: uno-a-muchos

Un ID de Ensembl puede mapear a varios símbolos (versiones distintas del gen, sinónimos legales). El argumento multiVals controla qué hacer:

  • "first" (default), devuelve el primero. Pierdes info pero el output es un vector simple.
  • "asNA": si hay más de uno, devuelve NA.
  • "list": devuelve lista por ID. Tienes que decidir qué hacer después.

Para análisis exploratorios, "first" está bien. Para reportes formales, conviene avisar al lector: “Cuando un ID Ensembl mapeaba a múltiples símbolos, se reportó el primero alfabéticamente.”

Anotación desde rowData

Si tu SummarizedExperiment vino de tximeta, parte de la anotación ya está en rowData(dds):

rowData(dds)
# DataFrame con columnas como gene_id, gene_name, gene_biotype...

En ese caso, en vez de mapIds, puedes usar directamente la metadata:

res$simbolo <- rowData(dds)$gene_name[match(rownames(res), rowData(dds)$gene_id)]

Más directo y consistente con tu transcriptoma exacto. Si dispones de rowData bien poblado, prefiérelo a mapIds, evita desajustes entre versiones de Ensembl.

Filtrar por biotipo

Una operación frecuente: quedarse solo con genes que codifican proteínas:

keep <- rowData(dds)$gene_biotype == "protein_coding"
dds_pc <- dds[keep, ]

Quita pseudogenes, lncRNA, microRNA, etc. Reduce ruido en análisis funcionales. Decide al inicio si quieres restringir o no, afecta a los p-valores ajustados (menos tests = menos pena por multiplicidad).

Enriquecimiento funcional con clusterProfiler

Una vez tienes símbolos y ENTREZID, el paso natural es preguntarse: ¿qué procesos biológicos están enriquecidos en mis genes diferenciales?

clusterProfiler es el paquete estándar. Una vista rápida del uso:

library(clusterProfiler)
library(org.Hs.eg.db)

# Lista de genes diferencialmente expresados
genes_up <- res_anotado |>
    filter(padj < 0.05, log2FoldChange > 1) |>
    pull(entrez_id) |>
    na.omit() |>
    unique()

# Lista universo (todos los testados)
universo <- res_anotado |>
    pull(entrez_id) |>
    na.omit() |>
    unique()

# Enriquecimiento GO (Biological Process)
go_bp <- enrichGO(gene          = genes_up,
                  universe      = universo,
                  OrgDb         = org.Hs.eg.db,
                  ont           = "BP",
                  pAdjustMethod = "BH",
                  pvalueCutoff  = 0.05,
                  qvalueCutoff  = 0.10)

head(as.data.frame(go_bp))

Argumentos clave:

  • gene: tus DE up (o down, o todos).
  • universe: los genes testeados, no todo el genoma. Crítico para no inflar enriquecimientos.
  • ont: ontología: "BP" (procesos), "MF" (función molecular), "CC" (componentes).

Visualización rápida:

dotplot(go_bp, showCategory = 15)
barplot(go_bp, showCategory = 15)

KEGG pathways:

kegg <- enrichKEGG(gene = genes_up,
                   universe = universo,
                   organism = "hsa",            # human; mmu, rno, dme...
                   pvalueCutoff = 0.05)

Esto es una pasada rápida, clusterProfiler tiene mucha más profundidad (GSEA, comparación de listas, networks). Para análisis serio, consulta su vignette.

La filosofía: universo bien elegido

El error más común en enriquecimiento: usar todo el genoma como universo. No.

El universo debe ser el conjunto de genes testeados, los que entraron a DESeq2 después del filtrado de baja expresión. Si testeaste 18.000 genes y usas 25.000 (todo el genoma) como universo, los enriquecimientos se inflan: los términos enriquecidos en los 7.000 genes filtrados aparecen como falsa señal.

Regla: universo = unique(na.omit(res$entrez_id)) tras mapIds sobre todos los genes del análisis.

Trampas habituales

  • Mezclar IDs de versiones distintas. Ensembl actualiza IDs (a veces). Si tu transcriptoma es de Ensembl 105 y org.Hs.eg.db se actualizó con Ensembl 110, parte del mapeo falla con NAs.
  • Versiones con sufijo. ENSG00000141510.18 (con .18) no mapea, org.Hs.eg.db espera ENSG00000141510. Quita el sufijo: sub("\\..*$", "", id).
  • Universo mal elegido en enriquecimiento. Probablemente el error más caro en bioinformática aplicada. Universo = testeados.
  • multiVals = "first" sin saberlo. Pierdes genes ambiguos. Para validación, comprueba cuántos NAs: sum(is.na(simbolos)).
  • Usar org.Hs.eg.db con datos de ratón. Suena obvio, pero pasa. El paquete del organismo debe coincidir con la muestra.

En la siguiente entrega

Última entrega de la ruta: un caso completo de RNA-seq, desde el airway cargado hasta el informe, pasando por filtrado, modelo, shrinkage, QC visual, volcano publicable, anotación y enriquecimiento. Todo unido, en el orden que pasa en un análisis real. Lo siguiente.