1 El escenario

Imagina que acabas de incorporarte a un laboratorio que estudia la enfermedad de Alzheimer. Tu asesora te da una lista de 20 genes de riesgo identificados por estudios GWAS y genética familiar, y te dice: “Necesito que descargues lecturas crudas de un estudio de RNA-seq en cerebro, anotes estos genes, explores su expresión, hagas un análisis de enriquecimiento y busques sus asociaciones genéticas. Para el viernes.”

Si hicieras esto manualmente — descargar archivos del SRA, abrir Ensembl, buscar cada gen, copiar la información, ir a GEO, descargar el dataset, abrir Reactome… — te tomaría días. Y si la lista crece a 200 genes, el enfoque manual simplemente no escala.

La respuesta: acceso programático. Hoy vamos a recorrer el camino completo — desde las lecturas crudas en la terminal hasta el análisis funcional en R — con un solo ejercicio de principio a fin usando los genes de riesgo de Alzheimer.

Nota conceptual — Manual vs. programático: El acceso manual (navegador web) es ideal para exploración y para entender qué datos existen. El acceso programático (terminal, R, Python) es indispensable cuando necesitas escala, reproducibilidad y automatización. No son enfoques que compitan — se complementan. Primero exploras a mano para entender (sesiones 07–08), luego automatizas para producir (esta sesión).

Al finalizar esta sesión, serás capaz de:

  1. Descargar e inspeccionar lecturas crudas desde el SRA usando la terminal.
  2. Acceder programáticamente a bases de datos genómicas desde R usando Bioconductor.
  3. Anotar genes masivamente con biomaRt, descargar datos de GEO con GEOquery, y hacer enriquecimiento con clusterProfiler.
  4. Crear un flujo de trabajo reproducible que integre datos de múltiples fuentes — todo con los mismos genes de Alzheimer.

2 El cambio de paradigma: de clics a código

Antes de escribir una sola línea de código, veamos exactamente qué ganamos:

Aspecto Manual (web) Programático (terminal + R)
Reproducibilidad Baja (depende de clics) Alta (código documentado)
Escalabilidad Un gen a la vez Miles de genes simultáneamente
Automatización No Sí (scripts, pipelines)
Documentación Capturas de pantalla Código fuente versionado
Integración Copiar y pegar Objetos en memoria

En el campo del Alzheimer, donde los estudios GWAS generan listas de cientos de loci de riesgo, el acceso programático no es un lujo — es una necesidad.

Nota de buenas prácticas: Cada paso de nuestro análisis quedará documentado en código, con versiones de las bases de datos registradas, y con resultados que cualquier persona puede regenerar ejecutando el mismo script.

3 Prerequisitos: instalación de paquetes

Antes de empezar, necesitas tener los paquetes instalados. Ejecuta este bloque una sola vez (no es necesario repetirlo cada sesión):

# Instalar BiocManager si no lo tienes
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Paquetes de Bioconductor
BiocManager::install(c(
    "biomaRt", "GEOquery", "org.Hs.eg.db",
    "clusterProfiler", "ReactomePA", "enrichplot",
    "limma", "AnnotationHub"
))

# Paquetes de CRAN
install.packages(c("tidyverse", "pheatmap", "gwasrapidd", "httr", "jsonlite"))

Nota: La instalación de Bioconductor puede tardar 10-15 minutos la primera vez. Si algún paquete falla, revisa los mensajes de error — frecuentemente se trata de dependencias del sistema que faltan.

4 FASE 1: La terminal — Obtener datos crudos del SRA

Nuestro primer paso como bioinformáticos: ir a la fuente. Antes de analizar datos procesados, vamos a descargar lecturas crudas de un estudio de RNA-seq en corteza cerebral humana — el tipo de dato que alimenta los análisis de expresión que después exploraremos en R.

Nota conceptual — ¿Por qué la terminal? Aunque es posible llamar herramientas de línea de comandos desde R con system(), para descargas y procesamiento de datos crudos la terminal es el entorno natural. En bioinformática, muchos flujos de trabajo combinan la terminal (para descargas y procesamiento de datos crudos) con R/RStudio (para análisis estadístico y visualización).

4.1 Paso 1 — Preparar el entorno

Abre una terminal en el servidor y carga las herramientas necesarias:

# Cargar el módulo de SRA Toolkit
module load sra/3.1.1
# module load: activa software preinstalado en el servidor (sistema de módulos).
# Sin este paso, fastq-dump no estaría disponible en tu PATH.

# Verificar que se cargó correctamente
which fastq-dump
# which: muestra la ruta completa del ejecutable. Si no devuelve nada, el módulo no se cargó.

# Crear una carpeta de trabajo
mkdir -p /mnt/atgc-d1/bioinfoII/$USER/alzheimer_ejercicio
# mkdir -p: crea el directorio y todos los intermedios necesarios (-p = parents).
# Si ya existe, no da error.

cd /mnt/atgc-d1/bioinfoII/$USER/alzheimer_ejercicio

Nota: Usamos /mnt/atgc-d1/bioinfoII/$USER/ porque es el espacio de trabajo asignado con suficiente almacenamiento. La variable $USER se reemplaza automáticamente por tu nombre de usuario.

Referencia: SRA Toolkit es desarrollado por NCBI. Documentación oficial: SRA Toolkit Wiki | SRA Toolkit GitHub.

4.2 Paso 2 — Descargar lecturas de RNA-seq de cerebro humano

Vamos a descargar un subconjunto de lecturas de un estudio de RNA-seq en corteza cerebral humana:

# Descargar las primeras 1,000 lecturas (subconjunto para no saturar el servidor)
fastq-dump --split-files -X 1000 SRR1382242
  • --split-files: separa las lecturas paired-end en dos archivos (_1.fastq y _2.fastq)
  • -X 1000: descarga solo las primeras 1,000 lecturas
  • SRR1382242: run de RNA-seq de corteza cerebral humana

La descarga puede tardar 1-3 minutos. Verás un mensaje como Read 1000 spots... cuando termine.

4.3 Paso 3 — Inspeccionar los archivos FASTQ

# Verificar los archivos descargados
ls -lh SRR1382242_*.fastq
# ls -lh: lista archivos con detalles (-l) y tamaños legibles (-h, ej. "5.2M").
# El asterisco (*) es un comodín que expande a todos los archivos que coincidan.

# Ver las primeras 2 lecturas (8 líneas = 2 lecturas × 4 líneas cada una)
head -8 SRR1382242_1.fastq
# head -N: muestra las primeras N líneas de un archivo.

# Contar el número total de lecturas
echo "Forward reads:" && grep -c "^@SRR" SRR1382242_1.fastq
echo "Reverse reads:" && grep -c "^@SRR" SRR1382242_2.fastq
# grep -c "^@SRR": cuenta (-c) las líneas que comienzan (^) con "@SRR" (identificador de lectura).

# Verificar la longitud de las lecturas
head -2 SRR1382242_1.fastq | tail -1 | wc -c
# Pipeline de 3 comandos: head -2 toma las primeras 2 líneas,
# tail -1 se queda con la última (la secuencia), wc -c cuenta caracteres.

Resultado esperado: Cada lectura ocupa 4 líneas en formato FASTQ (identificador, secuencia, separador +, calidades). Deberías ver 1,000 lecturas en cada archivo.

Reflexión:

  • ¿Por qué cada lectura ocupa exactamente 4 líneas en FASTQ?
  • Si un archivo FASTQ tiene 4,000 líneas, ¿cuántas lecturas contiene?
  • ¿Por qué descargamos solo 1,000 lecturas y no todas? (El run completo puede tener millones de lecturas y ocupar decenas de GB.)

Conexión con el ejercicio: Estos datos crudos son la materia prima. Después de alinear y cuantificar (pasos que aprenderás en otro módulo), se generan las matrices de expresión que descargamos ya procesadas con GEOquery en la Fase 2.

4.4 Paso 3.b — Práctica de formatos: de la teoría a la terminal

En la sesión 06 vimos los formatos fundamentales de bioinformática de manera teórica. Ahora vamos a tocarlos con las manos: descargar archivos reales y explorarlos en la terminal.

4.4.1 Tabla de referencia

Formato Extensión Contenido Coordenadas Uso principal
FASTA .fa, .fasta Secuencias N/A Genomas de referencia, secuencias proteicas
FASTQ .fq, .fastq Secuencias + calidad N/A Lecturas crudas de secuenciación
GFF3 .gff3, .gff Anotación genómica 1-based, closed Bases de datos (Ensembl, NCBI)
GTF .gtf Anotación genómica 1-based, closed RNA-seq (HTSeq, featureCounts)
BED .bed Regiones genómicas 0-based, half-open Intervalos, picos ChIP-seq
VCF .vcf Variantes 1-based SNPs, indels, variantes clínicas
SAM/BAM .sam, .bam Alineamientos 1-based Reads alineados al genoma

4.4.2 Ejercicio práctico: explorar formatos reales

Ya tienes archivos FASTQ del paso anterior. Ahora vamos a crear y explorar otros formatos usando datos del cromosoma 21 (donde se encuentra APP, uno de nuestros genes de Alzheimer).

# --- 1. FASTA: descargar la secuencia de APP desde Ensembl ---
# Descargamos solo la secuencia de la proteína APP (acceso rápido vía REST API)
curl -s "https://rest.ensembl.org/sequence/id/ENSP00000284981?type=protein" \
  -H "Content-type:text/x-fasta" > APP_protein.fasta
# curl: herramienta de línea de comandos para transferir datos desde/hacia un servidor.
#   -s: modo silencioso (no muestra barra de progreso).
#   -H: agrega un encabezado HTTP (aquí pedimos formato FASTA).
#   > : redirige la salida al archivo APP_protein.fasta.
# Documentación de la REST API de Ensembl: https://rest.ensembl.org/

# Inspeccionar: ¿cuántos aminoácidos tiene?
head -1 APP_protein.fasta                   # Línea de encabezado (>)
grep -v "^>" APP_protein.fasta | wc -c      # Contar caracteres (≈ aminoácidos)
# grep -v "^>": excluye (-v) las líneas de encabezado (que comienzan con >).
# wc -c: cuenta caracteres.

# --- 2. GFF3: descargar anotación del gen APP ---
curl -s "https://rest.ensembl.org/overlap/id/ENSG00000142192?feature=gene;feature=exon" \
  -H "Content-type:text/x-gff3" > APP_annotation.gff3
# El endpoint /overlap/id/ devuelve todos los features que se solapan con un ID dado.
# feature=gene;feature=exon pide tanto la anotación del gen como de cada exón.

# Inspeccionar: ¿cuántos exones tiene?
head -5 APP_annotation.gff3
grep -c "exon" APP_annotation.gff3
# grep -c: cuenta las líneas que contienen "exon" (en lugar de mostrarlas).

# --- 3. BED: crear un archivo BED a partir de nuestras coordenadas ---
# Recordatorio: BED usa coordenadas 0-based, half-open
# APP está en chr21:25,880,550-26,170,754 (1-based, Ensembl)
# En BED (0-based): chr21  25880549  26170754
echo -e "chr21\t25880549\t26170754\tAPP\t0\t+" > alzheimer_genes.bed
echo -e "chr21\t25891744\t26055015\tAPP_exon1-17\t0\t+" >> alzheimer_genes.bed
# echo -e: imprime texto interpretando secuencias de escape (\t = tabulador).
# > : crea el archivo (sobreescribe si existe).
# >>: agrega al final del archivo sin sobreescribirlo.
# BED tiene 6 columnas mínimas: cromosoma, inicio, fin, nombre, score, hebra.

# Inspeccionar
cat alzheimer_genes.bed
# cat: muestra el contenido completo del archivo.
# Pregunta: ¿por qué la posición de inicio es 25880549 y no 25880550?

# --- 4. VCF: descargar variantes de APP desde Ensembl ---
curl -s "https://rest.ensembl.org/overlap/id/ENSG00000142192?feature=variation" \
  -H "Content-type:application/json" | head -c 500
# | head -c 500: toma solo los primeros 500 caracteres de la salida (el JSON puede ser enorme).
# Nota: la API REST devuelve JSON, no VCF directamente.
# En ClinVar y dbSNP puedes descargar archivos VCF reales.

Resultado esperado: Después de este ejercicio deberías poder distinguir visualmente cada formato, entender sus columnas y la diferencia entre coordenadas 0-based (BED) y 1-based (GFF3, VCF).

Regla de oro: BED es el único formato común que usa coordenadas 0-based. Todos los demás (GFF3, GTF, VCF, SAM) usan 1-based. Si conviertes entre formatos, siempre ajusta las coordenadas o usa herramientas que lo hagan automáticamente (como bedtools o gffread).

Conexión con S06: En la sesión 06 aprendiste la estructura teórica de estos formatos. Ahora los has visto en la terminal con datos reales de nuestro caso de estudio (Alzheimer). En la Fase 2 verás cómo R puede leer y manipular estos datos programáticamente.

5 FASE 2: R — De genes a conocimiento

Ahora pasamos a R. Con las mismas preguntas sobre Alzheimer que motivaron la descarga de datos crudos, vamos a anotar los genes de riesgo, explorar su expresión, descubrir las rutas donde convergen y buscar sus asociaciones genéticas. Cada paso usa la salida del anterior.

5.1 Paso 4 — Anotar los genes con biomaRt (Ensembl)

En sesiones anteriores, buscar un solo gen en Ensembl tomaba varios minutos de clics. Ahora vamos a anotar 20 genes en segundos.

5.1.1 ¿Qué es BioMart?

biomaRt (Durinck et al., 2009, DOI: 10.1038/nprot.2009.97) es la interfaz de R para BioMart, el sistema de consultas masivas de Ensembl. La documentación completa y tutoriales están en la vignette oficial de biomaRt en Bioconductor.

Nota técnica — Vocabulario de BioMart: Necesitas conocer cuatro términos:

  • Mart: La base de datos a la que te conectas (ej. "genes" para el mart de genes de Ensembl).
  • Dataset: La especie dentro del mart (ej. "hsapiens_gene_ensembl" para humano).
  • Filters: Los criterios de búsqueda — lo que defines para filtrar (ej. gene name, chromosome name).
  • Attributes: La información que quieres obtener (ej. gene ID, descripción, GO terms).

Piensa en ello como una consulta SQL: los filters son el WHERE y los attributes son el SELECT.

5.1.2 Conectar a Ensembl y anotar un gen

# Si no están instalados, descomentar las siguientes líneas:
# BiocManager::install("biomaRt")
# install.packages(c("tidyverse", "httr"))

library(biomaRt)
library(tidyverse)

# Configurar timeout extendido (Ensembl puede tardar en responder)
options(timeout = 300)
httr::set_config(httr::timeout(300))

# Conectar al mart de genes de Ensembl para humano
ensembl <- useMart(
  biomart = "ENSEMBL_MART_ENSEMBL",
  dataset = "hsapiens_gene_ensembl",
  host = "https://www.ensembl.org"
)

# --- Nota sobre el código ---
# useMart() se conecta directamente a los servidores de Ensembl.
# Alternativa: useEnsembl() ofrece selección automática de mirrors,
# pero puede dar timeout en servidores con conexiones lentas.
# Si useMart() falla, intenta con un mirror específico:
# ensembl <- useMart(
#   biomart = "ENSEMBL_MART_ENSEMBL",
#   dataset = "hsapiens_gene_ensembl",
#   host = "https://useast.ensembl.org"
# )

# Empecemos con APOE: el gen de riesgo más importante para Alzheimer tardío
apoe_info <- getBM(
  attributes = c(
    "hgnc_symbol",           # Símbolo del gen (ej. APOE)
    "ensembl_gene_id",       # ID de Ensembl (ej. ENSG00000130203)
    "chromosome_name",       # Cromosoma
    "start_position",        # Posición de inicio
    "end_position",          # Posición de fin
    "strand",                # Hebra (+1 o -1)
    "gene_biotype",          # Tipo de gen (protein_coding, etc.)
    "description"            # Descripción funcional
  ),
  filters = "hgnc_symbol",
  values = "APOE",
  mart = ensembl
)

# --- Nota sobre el código ---
# getBM() es LA función central de biomaRt.
# "BM" = BioMart. Recibe: qué quieres (attributes),
# cómo filtrar (filters + values), y contra qué base (mart).

print(apoe_info)
##   hgnc_symbol ensembl_gene_id chromosome_name start_position end_position
## 1        APOE ENSG00000130203              19       44903787     44909396
##   strand   gene_biotype                                        description
## 1      1 protein_coding apolipoprotein E [Source:HGNC Symbol;Acc:HGNC:613]

Nota de buenas prácticas: Validar resultados programáticos contra datos conocidos es una práctica fundamental. Si sabes que APOE está en el cromosoma 19, verifica que el resultado lo confirme.

5.1.3 Escalar: de 1 gen a 20

Lo que manualmente tomaría una hora de copiar y pegar, aquí toma una línea más:

# Genes de riesgo de Alzheimer: genética familiar + GWAS
genes_alzheimer <- c(
  "APP", "PSEN1", "PSEN2", "APOE", "TREM2",
  "CLU", "BIN1", "ABCA7", "PICALM", "CR1",
  "CD33", "MS4A6A", "EPHA1", "CD2AP", "SPI1",
  "SORL1", "ADAM10", "MAPT", "GRN", "TOMM40"
)

# Convertir entre sistemas de identificadores — el cross-reference programático
conversion_ids <- getBM(
  attributes = c(
    "hgnc_symbol",       # Para humanos, el nombre estándar
    "ensembl_gene_id",   # Llave para Ensembl y AnnotationHub
    "entrezgene_id",     # Llave para NCBI, GEO, PubMed
    "uniprotswissprot"   # Llave para UniProt (solo curadas)
  ),
  filters = "hgnc_symbol",
  values = genes_alzheimer,
  mart = ensembl
)

print(conversion_ids)
##    hgnc_symbol ensembl_gene_id entrezgene_id uniprotswissprot
## 1        EPHA1 ENSG00000284816          2041           P21709
## 2         MAPT ENSG00000277956          4137                 
## 3         MAPT ENSG00000277956          4137           P10636
## 4         MAPT ENSG00000276155          4137           P10636
## 5         MAPT ENSG00000276155          4137                 
## 6          APP ENSG00000142192           351           P05067
## 7          APP ENSG00000142192           351                 
## 8        ABCA7 ENSG00000064687         10347           Q8IZY2
## 9        ABCA7 ENSG00000064687         10347                 
## 10        APOE ENSG00000130203           348                 
## 11        APOE ENSG00000130203           348           P02649
## 12        SPI1 ENSG00000066336          6688                 
## 13        SPI1 ENSG00000066336          6688           P17947
## 14       TREM2 ENSG00000095970         54209           Q9NZC2
## 15       TREM2 ENSG00000095970         54209                 
## 16       SORL1 ENSG00000137642          6653           Q92673
## 17       SORL1 ENSG00000137642          6653                 
## 18      ADAM10 ENSG00000137845           102           O14672
## 19      ADAM10 ENSG00000137845           102                 
## 20       EPHA1 ENSG00000146904          2041           P21709
## 21       EPHA1 ENSG00000146904          2041                 
## 22       PSEN1 ENSG00000080815          5663                 
## 23       PSEN1 ENSG00000080815          5663           P49768
## 24        CD33 ENSG00000105383           945                 
## 25        CD33 ENSG00000105383           945           P20138
## 26      TOMM40 ENSG00000130204         10452                 
## 27      TOMM40 ENSG00000130204         10452           O96008
## 28         CR1 ENSG00000203710          1378           P17927
## 29         CR1 ENSG00000203710          1378                 
## 30         CLU ENSG00000120885          1191           P10909
## 31         CLU ENSG00000120885          1191                 
## 32         GRN ENSG00000030582          2896                 
## 33         GRN ENSG00000030582          2896           P28799
## 34       CD2AP ENSG00000198087         23607           Q9Y5K6
## 35       CD2AP ENSG00000198087         23607                 
## 36        BIN1 ENSG00000136717           274           O00499
## 37        BIN1 ENSG00000136717           274                 
## 38      PICALM ENSG00000073921          8301           Q13492
## 39      PICALM ENSG00000073921          8301                 
## 40      MS4A6A ENSG00000110077         64231           Q9H2W1
## 41      MS4A6A ENSG00000110077         64231                 
## 42        MAPT ENSG00000186868          4137                 
## 43        MAPT ENSG00000186868          4137           P10636
## 44       PSEN2 ENSG00000143801          5664           P49810
## 45       PSEN2 ENSG00000143801          5664
# Guardar para uso en pasos posteriores
write_csv(conversion_ids, "tabla_conversion_IDs_alzheimer.csv")

Nota técnica — Mapeo de identificadores: Un gen puede tener múltiples IDs (transcritos, isoformas). No te alarmes si tu tabla tiene más filas que genes — es normal. La clave es entender la relación uno-a-muchos.

5.1.4 Ortólogos en ratón

Los modelos murinos son esenciales en la investigación de Alzheimer (ratones transgénicos APP/PS1, 5xFAD, etc.). biomaRt puede buscar ortólogos:

# Consulta 1: IDs de Ensembl (página feature)
genes_ids <- tryCatch(
  getBM(
    attributes = c("hgnc_symbol", "ensembl_gene_id"),
    filters = "hgnc_symbol",
    values = genes_alzheimer,
    mart = ensembl
  ),
  error = function(e) {
    message("Error al consultar Ensembl: ", conditionMessage(e))
    data.frame(hgnc_symbol = character(), ensembl_gene_id = character())
  }
)

# Consulta 2: Ortólogos en ratón (página homologs)
if (nrow(genes_ids) > 0) {
  ortologos_raw <- tryCatch(
    getBM(
      attributes = c(
        "ensembl_gene_id",
        "mmusculus_homolog_ensembl_gene",
        "mmusculus_homolog_associated_gene_name",
        "mmusculus_homolog_perc_id"
      ),
      filters = "ensembl_gene_id",
      values = genes_ids$ensembl_gene_id,
      mart = ensembl
    ),
    error = function(e) {
      message("Error al consultar ortólogos: ", conditionMessage(e))
      data.frame(ensembl_gene_id = character(),
                 mmusculus_homolog_ensembl_gene = character(),
                 mmusculus_homolog_associated_gene_name = character(),
                 mmusculus_homolog_perc_id = numeric())
    }
  )

  # --- Nota sobre el código ---
  # biomaRt no permite mezclar atributos de distintas "páginas" en una sola
  # consulta. hgnc_symbol está en "feature" y los homólogos en "homologs".
  # Solución: dos consultas y unirlas con merge().
  ortologos_raton <- merge(genes_ids, ortologos_raw, by = "ensembl_gene_id")

  ortologos_raton %>%
    dplyr::filter(mmusculus_homolog_ensembl_gene != "") %>%
    print()
}
##    ensembl_gene_id hgnc_symbol mmusculus_homolog_ensembl_gene
## 1  ENSG00000030582         GRN             ENSMUSG00000034708
## 2  ENSG00000064687       ABCA7             ENSMUSG00000035722
## 3  ENSG00000066336        SPI1             ENSMUSG00000002111
## 4  ENSG00000073921      PICALM             ENSMUSG00000039361
## 5  ENSG00000080815       PSEN1             ENSMUSG00000019969
## 6  ENSG00000095970       TREM2             ENSMUSG00000023992
## 7  ENSG00000110077      MS4A6A             ENSMUSG00000024677
## 8  ENSG00000110077      MS4A6A             ENSMUSG00000079419
## 9  ENSG00000110077      MS4A6A             ENSMUSG00000024679
## 10 ENSG00000120885         CLU             ENSMUSG00000022037
## 11 ENSG00000130203        APOE             ENSMUSG00000002985
## 12 ENSG00000130204      TOMM40             ENSMUSG00000002984
## 13 ENSG00000136717        BIN1             ENSMUSG00000024381
## 14 ENSG00000137642       SORL1             ENSMUSG00000049313
## 15 ENSG00000137845      ADAM10             ENSMUSG00000054693
## 16 ENSG00000142192         APP             ENSMUSG00000022892
## 17 ENSG00000143801       PSEN2             ENSMUSG00000010609
## 18 ENSG00000146904       EPHA1             ENSMUSG00000029859
## 19 ENSG00000186868        MAPT             ENSMUSG00000018411
## 20 ENSG00000198087       CD2AP             ENSMUSG00000061665
## 21 ENSG00000203710         CR1             ENSMUSG00000016481
##    mmusculus_homolog_associated_gene_name mmusculus_homolog_perc_id
## 1                                     Grn                   74.8735
## 2                                   Abca7                   77.7726
## 3                                    Spi1                   88.5185
## 4                                  Picalm                   97.6994
## 5                                   Psen1                   92.7195
## 6                                   Trem2                   68.2609
## 7                                  Ms4a6b                   51.5556
## 8                                  Ms4a6c                   47.1111
## 9                                  Ms4a6d                   53.3333
## 10                                    Clu                   76.3920
## 11                                   Apoe                   70.6625
## 12                                 Tomm40                   93.0748
## 13                                   Bin1                   77.2344
## 14                                  Sorl1                   93.2249
## 15                                 Adam10                   96.1230
## 16                                    App                   94.8052
## 17                                  Psen2                   95.9821
## 18                                  Epha1                   88.2172
## 19                                   Mapt                   65.4262
## 20                                  Cd2ap                   86.2285
## 21                                   Cr1l                    9.6826

Conexión con el ejercicio: Ahora tenemos los identificadores de nuestros 20 genes. Con estos IDs podemos buscar su expresión en GEO (Paso 5) y hacer enriquecimiento funcional (Paso 6). El cross-reference programático en acción.

5.2 Paso 5 — Explorar la expresión con GEOquery

¿Los genes de riesgo de Alzheimer están diferencialmente expresados en cerebros de pacientes? Vamos a descargar un dataset completo y explorarlo.

Conexión: En la Fase 1 descargamos lecturas crudas de RNA-seq. Esas lecturas, después de procesarse (alineamiento + cuantificación), generan exactamente el tipo de matriz que GEOquery descarga ya lista.

5.2.1 Descargar el dataset

GEOquery (Davis & Meltzer, 2007) permite descargar datasets completos de Gene Expression Omnibus directamente a R. Documentación: vignette oficial de GEOquery.

# Si no están instalados, descomentar las siguientes líneas:
# BiocManager::install("GEOquery")
# install.packages("pheatmap")

library(GEOquery)
library(pheatmap)

# Descargar GSE33000: expresión en corteza prefrontal
# Narayanan et al. (2014) — 310 Alzheimer + 157 controles
gse <- getGEO("GSE33000", GSEMatrix = TRUE, AnnotGPL = TRUE)
# getGEO(): descarga un dataset completo de GEO.
#   GSEMatrix = TRUE: descarga la matriz de expresión ya procesada.
#   AnnotGPL = TRUE: incluye la anotación de la plataforma (qué gen mide cada sonda).
# Devuelve una lista de objetos ExpressionSet (uno por plataforma).

gse <- gse[[1]]
# [[1]]: tomamos el primer (y normalmente único) ExpressionSet de la lista.

# Registrar la descarga — buena práctica para reproducibilidad
cat("Dataset:", "GSE33000", "\n")
## Dataset: GSE33000
cat("Título:", experimentData(gse)@title, "\n")
## Título: Gene expression profiles of human prefrontal cortex brain tissues
# experimentData(): extrae los metadatos generales del experimento (título, autores, etc.).
cat("Muestras:", ncol(gse), "| Sondas:", nrow(gse), "\n")
## Muestras: 624 | Sondas: 39280
cat("Fecha de descarga:", as.character(Sys.Date()), "\n")
## Fecha de descarga: 2026-03-05

5.2.2 Extraer las piezas del dataset

Un objeto ExpressionSet tiene tres componentes principales que se acceden con funciones específicas:

# 1. Matriz de expresión: las mediciones numéricas
expr_matrix <- exprs(gse)
# exprs(): extrae la matriz de expresión (genes × muestras) del ExpressionSet.
# Cada fila es una sonda/gen, cada columna es una muestra.
cat("Dimensiones:", nrow(expr_matrix), "sondas x", ncol(expr_matrix), "muestras\n")
## Dimensiones: 39280 sondas x 624 muestras
# 2. Metadatos de las muestras: ¿quién es quién?
pheno_data <- pData(gse)
# pData(): extrae los datos fenotípicos de las muestras (phenotype data).
# Contiene la información clínica: condición, tejido, edad, etc.
pheno_data %>%
  dplyr::select(title, geo_accession, source_name_ch1) %>%
  head(10)
##             title geo_accession  source_name_ch1
## GSM1423780  PFC_1    GSM1423780  HBTRC_PF_Pool_1
## GSM1423781  PFC_2    GSM1423781  HBTRC_PF_Pool_2
## GSM1423782  PFC_3    GSM1423782  HBTRC_PF_Pool_3
## GSM1423783  PFC_4    GSM1423783  HBTRC_PF_Pool_4
## GSM1423784  PFC_5    GSM1423784  HBTRC_PF_Pool_5
## GSM1423785  PFC_6    GSM1423785  HBTRC_PF_Pool_6
## GSM1423786  PFC_7    GSM1423786  HBTRC_PF_Pool_7
## GSM1423787  PFC_8    GSM1423787  HBTRC_PF_Pool_8
## GSM1423788  PFC_9    GSM1423788  HBTRC_PF_Pool_9
## GSM1423789 PFC_10    GSM1423789 HBTRC_PF_Pool_10
# 3. Anotación de las sondas: ¿qué gen mide cada sonda?
feat_data <- fData(gse)
# fData(): extrae los datos de las features (feature data).
# Mapea cada sonda a su gen, cromosoma, GO terms, etc.

5.2.3 Heatmap de genes de Alzheimer

Conectamos los datos de GEO con nuestra lista de genes (Paso 4):

# Buscar las sondas que corresponden a nuestros genes
gene_col <- grep("Gene.symbol|Gene symbol|GENE_SYMBOL|gene_assignment",
                 colnames(feat_data), value = TRUE, ignore.case = TRUE)[1]

if (!is.na(gene_col)) {
  sondas_interes <- feat_data %>%
    dplyr::filter(!!sym(gene_col) %in% genes_alzheimer)
  cat("Sondas encontradas:", nrow(sondas_interes), "para",
      length(unique(sondas_interes[[gene_col]])), "genes\n")
} else {
  gene_col_alt <- colnames(feat_data)[
    sapply(feat_data, function(x) any(grepl("APOE", x)))
  ][1]
  if (!is.na(gene_col_alt)) {
    sondas_interes <- feat_data %>%
      dplyr::filter(!!sym(gene_col_alt) %in% genes_alzheimer)
    gene_col <- gene_col_alt
  }
}
## Sondas encontradas: 32 para 19 genes
# Generar el heatmap
if (exists("sondas_interes") && nrow(sondas_interes) > 0) {
  sondas_validas <- intersect(rownames(sondas_interes), rownames(expr_matrix))
  if (length(sondas_validas) == 0) {
    sondas_validas <- intersect(sondas_interes$ID, rownames(expr_matrix))
  }

  if (length(sondas_validas) > 0) {
    expr_interes <- expr_matrix[sondas_validas, , drop = FALSE]

    row_labels <- sondas_interes[[gene_col]][
      match(rownames(expr_interes), rownames(sondas_interes))
    ]
    if (any(is.na(row_labels))) {
      row_labels <- sondas_interes[[gene_col]][
        match(rownames(expr_interes), sondas_interes$ID)
      ]
    }
    rownames(expr_interes) <- make.unique(row_labels)

    # pheatmap(): crea un heatmap con clustering jerárquico.
    # Documentación: https://cran.r-project.org/package=pheatmap
    #   scale = "row": normaliza cada fila (gen) a media 0 y SD 1, para
    #     comparar patrones de expresión entre muestras independientemente de la magnitud.
    #   clustering_distance_rows = "correlation": agrupa genes por similitud de perfil.
    #   show_colnames = FALSE: oculta nombres de muestras (son demasiadas para leer).
    #   colorRampPalette(): crea un gradiente de color suave de 100 valores.
    pheatmap(
      expr_interes,
      scale = "row",
      clustering_distance_rows = "correlation",
      show_colnames = FALSE,
      main = "Expresión de genes de riesgo de Alzheimer (corteza prefrontal)",
      fontsize_row = 10,
      color = colorRampPalette(c("#003366", "white", "#AC1911"))(100)
    )
  }
}
*Heatmap* de expresión de genes de riesgo de Alzheimer en GSE33000 (corteza prefrontal). Colores: azul = sub-expresión, blanco = nivel medio, carmesí = sobre-expresión.

Heatmap de expresión de genes de riesgo de Alzheimer en GSE33000 (corteza prefrontal). Colores: azul = sub-expresión, blanco = nivel medio, carmesí = sobre-expresión.

Nota conceptual — Interpretando el heatmap: El clustering agrupa genes con patrones de expresión similares. En Alzheimer, esperamos ver que genes involucrados en procesos similares (ej. procesamiento de amiloide: APP, PSEN1, ADAM10; o neuroinflamación: TREM2, CD33, CR1) se agrupen juntos.

5.2.4 PCA: ¿las muestras se separan?

# PCA (Análisis de Componentes Principales) reduce miles de dimensiones a unas pocas.
expr_clean <- expr_matrix[complete.cases(expr_matrix), ]
# complete.cases(): elimina filas (sondas) con valores NA — prcomp no tolera NAs.

pca_result <- prcomp(t(expr_clean), scale. = TRUE)
# prcomp(): calcula el PCA. Es la función base de R para esto.
#   t(): transpone la matriz porque prcomp espera muestras en filas y genes en columnas.
#   scale. = TRUE: estandariza cada gen (media=0, sd=1) antes del PCA, para que
#     genes altamente expresados no dominen el resultado.

pca_df <- as.data.frame(pca_result$x[, 1:2])
# pca_result$x: contiene las coordenadas de cada muestra en el nuevo espacio de componentes.
# [, 1:2]: tomamos solo PC1 y PC2 (los dos ejes de mayor varianza).

# Buscar la columna de condición (Alzheimer vs. Control)
condition_col <- NULL
for (col in colnames(pheno_data)) {
  vals <- unique(pheno_data[[col]])
  if (length(vals) >= 2 && length(vals) <= 10) {
    vals_lower <- tolower(as.character(vals))
    if (any(grepl("alzheimer|control|normal|disease|dementia|ad$|^ad ", vals_lower))) {
      condition_col <- col
      break
    }
  }
}

if (is.null(condition_col)) {
  candidates <- grep("disease|status|diagnosis|condition|group|type",
                     colnames(pheno_data), value = TRUE, ignore.case = TRUE)
  for (col in candidates) {
    if (length(unique(pheno_data[[col]])) >= 2 &&
        length(unique(pheno_data[[col]])) <= 10) {
      condition_col <- col
      break
    }
  }
}

if (!is.null(condition_col)) {
  pca_df$condition <- pheno_data[[condition_col]][
    match(rownames(pca_df), rownames(pheno_data))
  ]
}

var_explained <- round(100 * pca_result$sdev^2 / sum(pca_result$sdev^2), 1)

n_conditions <- length(unique(na.omit(pca_df$condition)))
if (!is.null(condition_col) && n_conditions >= 2 && n_conditions <= 10) {
  ggplot(pca_df, aes(x = PC1, y = PC2, color = condition)) +
    geom_point(size = 3, alpha = 0.7) +
    labs(
      title = "PCA — Corteza prefrontal (GSE33000)",
      x = paste0("PC1 (", var_explained[1], "% varianza)"),
      y = paste0("PC2 (", var_explained[2], "% varianza)"),
      color = "Condición"
    ) +
    theme_minimal(base_size = 14)
} else {
  ggplot(pca_df, aes(x = PC1, y = PC2)) +
    geom_point(size = 3, color = "#003366") +
    labs(
      title = "PCA — Corteza prefrontal (GSE33000)",
      x = paste0("PC1 (", var_explained[1], "% varianza)"),
      y = paste0("PC2 (", var_explained[2], "% varianza)")
    ) +
    theme_minimal(base_size = 14)
}
PCA de las muestras de corteza prefrontal (GSE33000). Cada punto es una muestra.

PCA de las muestras de corteza prefrontal (GSE33000). Cada punto es una muestra.

Nota conceptual — PCA en neurogenómica: En enfermedades neurodegenerativas, la separación entre pacientes y controles en el PCA suele ser menos clara que en cáncer, porque los cambios de expresión son más sutiles. No ver una separación nítida no significa que no haya diferencias. Buena práctica: Siempre revisa un PCA como primer paso de control de calidad.

Conexión con el ejercicio: El heatmap nos muestra qué genes cambian. Ahora la pregunta es: ¿convergen esos genes en rutas biológicas coherentes? Para eso necesitamos el análisis de enriquecimiento (Paso 6).

5.3 Paso 6 — Enriquecimiento funcional con clusterProfiler y ReactomePA

Tenemos una lista de genes de riesgo de Alzheimer. La pregunta no es “¿cuáles son?” sino “¿qué rutas biológicas convergen?”. Los genes individuales son piezas de un rompecabezas — el enriquecimiento nos muestra la imagen completa.

5.3.1 ORA vs. GSEA: dos filosofías

Nota conceptual — Dos enfoques para enriquecimiento:

  • ORA (Over-Representation Analysis): tienes una lista de genes de interés y preguntas “¿hay rutas biológicas sobre-representadas?”. Usa la prueba hipergeométrica.
  • GSEA (Gene Set Enrichment Analysis): tienes todos los genes ranqueados por una métrica (ej. log2 fold change) y preguntas “¿hay rutas cuyos genes tienden a estar arriba o abajo del ranking?”. No requiere umbral.
Método Input Ventaja Función en R
ORA Lista de genes Simple, intuitivo enrichGO(), enrichKEGG(), enrichPathway()
GSEA Ranking completo No pierde información gseGO(), gseKEGG(), gsePathway()

En la práctica, es bueno hacer ambos y comparar.

5.3.2 Preparar los genes y convertir IDs

Usaremos cuatro paquetes complementarios:

# Si no están instalados, descomentar las siguientes líneas:
# BiocManager::install(c("clusterProfiler", "ReactomePA", "org.Hs.eg.db", "enrichplot"))

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

# Lista expandida de genes asociados a Alzheimer
genes_symbols <- c(
  # Genética familiar (inicio temprano)
  "APP", "PSEN1", "PSEN2",
  # Principal factor de riesgo (inicio tardío)
  "APOE", "TOMM40",
  # Inmunidad innata y microglía
  "TREM2", "CD33", "CR1", "SPI1", "INPP5D", "PLCG2",
  # Endocitosis y tráfico vesicular
  "BIN1", "PICALM", "CD2AP", "SORL1", "RIN3",
  # Metabolismo lipídico
  "CLU", "ABCA7", "ABCA1",
  # Señalización y sinapsis
  "EPHA1", "PTK2B", "ADAM10",
  # Citoesqueleto y tau
  "MAPT", "CASS4",
  # Otros loci GWAS
  "MS4A6A", "GRN", "MEF2C", "ZCWPW1",
  # Procesamiento de amiloide
  "BACE1", "BACE2", "NCSTN", "APH1A"
)

# bitr() = "Biological Id TRanslator" — convierte entre tipos de identificadores.
# clusterProfiler necesita Entrez IDs, no gene symbols.
gene_conversion <- bitr(
  genes_symbols,
  fromType = "SYMBOL",      # Tipo de ID de entrada
  toType = c("ENTREZID", "ENSEMBL"),  # Tipos de ID de salida
  OrgDb = org.Hs.eg.db      # Base de datos de anotación (humano)
)
# bitr() usa org.Hs.eg.db como diccionario de traducción entre IDs.
# Si un símbolo no tiene equivalente Entrez, se omite con un warning.

entrez_ids <- unique(gene_conversion$ENTREZID)
cat("Genes convertidos exitosamente:", length(entrez_ids), "IDs únicos de",
    length(genes_symbols), "símbolos\n")
## Genes convertidos exitosamente: 32 IDs únicos de 32 símbolos

Nota de buenas prácticas: bitr() puede fallar silenciosamente para algunos genes. Siempre verifica cuántos genes se convirtieron vs. cuántos enviaste.

5.3.3 Enriquecimiento Gene Ontology (ORA)

# enrichGO(): realiza Over-Representation Analysis (ORA) contra Gene Ontology.
# Pregunta: ¿hay categorías GO con más genes de mi lista de lo esperado por azar?
ego_bp <- enrichGO(
  gene = entrez_ids,        # Nuestros genes de interés (Entrez IDs)
  OrgDb = org.Hs.eg.db,    # Base de datos de anotación (mapeo gen → GO)
  ont = "BP",               # Ontología: BP (Biological Process), MF o CC
  pAdjustMethod = "BH",    # Corrección de Benjamini-Hochberg para tests múltiples
  pvalueCutoff = 0.05,     # Umbral de significancia (p ajustado)
  readable = TRUE           # Convierte Entrez IDs a gene symbols en el resultado
)

# dotplot(): visualiza las categorías enriquecidas.
# Tamaño del punto = número de genes; color = p-value ajustado.
dotplot(ego_bp, showCategory = 15) +
  ggtitle("GO Biological Process — Genes de riesgo de Alzheimer")
*Dotplot* de Gene Ontology (Biological Process) — top 15 categorías enriquecidas en genes de riesgo de Alzheimer.

Dotplot de Gene Ontology (Biological Process) — top 15 categorías enriquecidas en genes de riesgo de Alzheimer.

Nota conceptual — ¿Qué esperamos ver? Procesamiento de APP (APP, PSEN1, ADAM10, BACE1), respuesta inmune innata (TREM2, CD33, CR1), endocitosis (BIN1, PICALM, SORL1), y metabolismo lipídico (APOE, CLU, ABCA7).

5.3.4 Enriquecimiento Reactome

# enrichPathway(): ORA contra la base de datos Reactome (rutas curadas manualmente).
# Similar a enrichGO() pero usando rutas biológicas de Reactome en lugar de GO.
ereactome <- enrichPathway(
  gene = entrez_ids,        # Genes de interés (Entrez IDs)
  organism = "human",       # Organismo
  pvalueCutoff = 0.05,     # Umbral de significancia
  readable = TRUE           # Convierte IDs a gene symbols
)

dotplot(ereactome, showCategory = 15) +
  ggtitle("Reactome Pathways — Genes de riesgo de Alzheimer")
*Dotplot* de rutas de Reactome enriquecidas en genes de riesgo de Alzheimer.

Dotplot de rutas de Reactome enriquecidas en genes de riesgo de Alzheimer.

5.3.5 Visualizaciones avanzadas

# cnetplot(): red bipartita gen–término. Muestra qué genes participan en qué categorías.
# Nodos grandes = categorías GO, nodos pequeños = genes.
# Permite ver genes "hub" que conectan múltiples procesos biológicos.
cnetplot(ego_bp, showCategory = 5) +
  ggtitle("Red Gen-Término GO (top 5 categorías)")
Red gen-término GO (*cnetplot*): muestra qué genes participan en qué categorías.

Red gen-término GO (cnetplot): muestra qué genes participan en qué categorías. 

# pairwise_termsim(): calcula la similitud semántica entre categorías GO.
# Necesario antes de emapplot() para saber qué categorías agrupar.
ego_bp_pairwise <- pairwise_termsim(ego_bp)

# emapplot(): mapa de enriquecimiento donde categorías similares se agrupan cerca.
# Útil para identificar "clusters temáticos" (ej. inmunidad, metabolismo, señalización).
emapplot(ego_bp_pairwise, showCategory = 20) +
  ggtitle("Mapa de enriquecimiento GO BP")
Red gen-término GO (*cnetplot*): muestra qué genes participan en qué categorías.

Red gen-término GO (cnetplot): muestra qué genes participan en qué categorías.

5.3.6 GSEA: usando el ranking completo

# Para GSEA necesitas un vector NOMBRADO con un ranking numérico.
# En un análisis real, el ranking sería el log2 fold change de un análisis diferencial.
# Aquí simulamos un ranking como ejemplo didáctico.
set.seed(42)  # Fijar semilla para reproducibilidad
all_entrez <- keys(org.Hs.eg.db, "ENTREZID")
# keys(): extrae todos los IDs de un tipo dado de la base de datos de anotación.

background_ids <- sample(setdiff(all_entrez, entrez_ids), 500)
# setdiff(): genes del background que NO están en nuestra lista.
# sample(): selecciona 500 genes aleatorios como background.

gene_list <- setNames(
  c(rnorm(length(entrez_ids), mean = 2, sd = 0.5),
    rnorm(500, mean = 0, sd = 0.3)),
  c(entrez_ids, background_ids)
)
# setNames(): asigna nombres (IDs) a cada valor del vector.
# Nuestros genes de Alzheimer reciben valores altos (simulando sobreexpresión).

gene_list <- sort(gene_list, decreasing = TRUE)
# GSEA requiere el vector ORDENADO de mayor a menor.

# gseGO(): realiza Gene Set Enrichment Analysis contra Gene Ontology.
# A diferencia de enrichGO() (ORA), no requiere un umbral — usa TODOS los genes ranqueados.
gsea_go <- gseGO(
  geneList = gene_list,    # Vector nombrado y ordenado
  OrgDb = org.Hs.eg.db,   # Base de datos de anotación
  ont = "BP",              # Ontología: Biological Process
  pvalueCutoff = 0.05,    # Umbral de significancia
  verbose = FALSE          # No imprimir mensajes de progreso
)

if (nrow(gsea_go@result %>% dplyr::filter(p.adjust < 0.05)) > 0) {
  # gseaplot2(): genera el gráfico clásico de GSEA con tres paneles:
  #   1) Enrichment score a lo largo del ranking.
  #   2) Posición de los genes del gene set en el ranking (barcode plot).
  #   3) Señal de ranking (métrica de cada gen).
  gseaplot2(gsea_go, geneSetID = 1,
            title = gsea_go@result$Description[1])
}
Gráfico GSEA clásico mostrando el enrichment score.

Gráfico GSEA clásico mostrando el enrichment score.

Nota conceptual — El gráfico GSEA: La curva verde muestra el “enrichment score” a lo largo del ranking. Si los genes de una ruta están concentrados en la parte superior, la curva sube — el pico es el enrichment score.

Conexión con el ejercicio: Las rutas enriquecidas nos dicen por qué convergen estos genes. Ahora veamos la evidencia genética directa: ¿qué dice el GWAS Catalog? (Paso 7).

5.4 Paso 7 — Asociaciones GWAS con gwasrapidd

Los estudios GWAS han identificado docenas de loci de riesgo para Alzheimer. Con gwasrapidd (Magno & Couto, 2020, DOI: 10.1093/bioinformatics/btz605), exploramos el GWAS Catalog programáticamente para APOE y TREM2. Documentación: vignette oficial de gwasrapidd.

# Si no está instalado, descomentar la siguiente línea:
# install.packages("gwasrapidd")

library(gwasrapidd)

# === APOE: el gen de riesgo más fuerte ===
# get_variants(): consulta la REST API del GWAS Catalog y devuelve un objeto S4
# con las variantes (SNPs) asociadas a un gen. Contiene: variant_id (rsID),
# chromosome, position, functional_class, etc.
apoe_variants <- tryCatch(
  get_variants(gene_name = "APOE"),
  error = function(e) {
    message("No se pudo conectar al GWAS Catalog para APOE: ", conditionMessage(e))
    NULL
  }
)

if (!is.null(apoe_variants) && nrow(apoe_variants@variants) > 0) {
  cat("Variantes GWAS encontradas para APOE:",
      nrow(apoe_variants@variants), "\n")

  variant_ids_apoe <- apoe_variants@variants$variant_id
  # get_associations(): obtiene los estudios y métricas de asociación para cada variante.
  # Devuelve: pvalue, odds ratio (or_per_copy_number), beta, tamaño de efecto.
  apoe_gwas <- tryCatch(
    get_associations(variant_id = variant_ids_apoe),
    error = function(e) {
      message("No se pudieron obtener asociaciones: ", conditionMessage(e))
      NULL
    }
  )

  if (!is.null(apoe_gwas)) {
    cat("Asociaciones GWAS para APOE:",
        nrow(apoe_gwas@associations), "\n")
    apoe_gwas@associations %>%
      dplyr::select(association_id, pvalue, or_per_copy_number, beta_number) %>%
      dplyr::arrange(pvalue) %>%
      head(10)
  }
} else {
  cat("GWAS Catalog no disponible en este momento.\n")
  cat("Consulta manualmente: https://www.ebi.ac.uk/gwas/genes/APOE\n")
}
## Variantes GWAS encontradas para APOE: 620 
## Asociaciones GWAS para APOE: 9504
## # A tibble: 10 × 4
##    association_id pvalue or_per_copy_number beta_number
##    <chr>           <dbl>              <dbl>       <dbl>
##  1 87304303            0                 NA       0.206
##  2 87309478            0                 NA       0.253
##  3 93377798            0                 NA       0.206
##  4 87301391            0                 NA       0.179
##  5 96179228            0                 NA       0.371
##  6 96181295            0                 NA       0.383
##  7 96227020            0                 NA       0.361
##  8 96184445            0                 NA       0.390
##  9 96226246            0                 NA       0.372
## 10 96170842            0                 NA       0.347
# === TREM2: inmunidad innata y microglía ===
trem2_variants <- tryCatch(
  get_variants(gene_name = "TREM2"),
  error = function(e) {
    message("No se pudo conectar al GWAS Catalog para TREM2: ", conditionMessage(e))
    NULL
  }
)

if (!is.null(trem2_variants) && nrow(trem2_variants@variants) > 0) {
  trem2_gwas <- tryCatch(
    get_associations(variant_id = trem2_variants@variants$variant_id),
    error = function(e) NULL
  )
  if (!is.null(trem2_gwas)) {
    cat("\nAsociaciones GWAS para TREM2:",
        nrow(trem2_gwas@associations), "\n")
  }
}
## 
## Asociaciones GWAS para TREM2: 109

Nota conceptual — APOE y el riesgo de Alzheimer: APOE tiene tres alelos principales: ε2 (protector), ε3 (neutro), y ε4 (riesgo). Un portador heterocigoto de ε4 tiene ~3x más riesgo; un homocigoto ε4/ε4 tiene ~12x más riesgo. Las variantes rs429358 y rs7412 definen estos alelos — son de las asociaciones genéticas más fuertes y reproducibles en toda la medicina.

Nota conceptual — TREM2 y la hipótesis neuroinflamatoria: El descubrimiento de variantes de riesgo en TREM2 (un receptor de microglía) demostró que la neuroinflamación no es solo una consecuencia del Alzheimer, sino una causa activa. Esto conecta directamente con las rutas de inmunidad innata que encontramos enriquecidas en el Paso 6.

5.5 Paso 8 — AnnotationHub: recursos centralizados

¿Qué pasa cuando necesitas un recurso de anotación que no está en biomaRt? AnnotationHub es un portal centralizado con miles de recursos de anotación genómica. Documentación: vignette oficial de AnnotationHub.

# Si no está instalado, descomentar la siguiente línea:
# BiocManager::install("AnnotationHub")

library(AnnotationHub)

# AnnotationHub(): crea una conexión al repositorio central de recursos.
# La primera vez descarga los metadatos del catálogo (puede tardar unos segundos).
ah <- AnnotationHub()
cat("Total de recursos:", length(ah), "\n")
## Total de recursos: 70865
# query(): filtra los recursos por palabras clave.
# Busca en título, especie, proveedor, tipo, etc.
ah_human <- query(ah, "Homo sapiens")
cat("Recursos para humano:", length(ah_human), "\n")
## Recursos para humano: 26722
# Combinar criterios: Ensembl databases para humano
ah_ensembl <- query(ah, c("Homo sapiens", "EnsDb"))
ah_ensembl
## AnnotationHub with 28 records
## # snapshotDate(): 2025-10-29
## # $dataprovider: Ensembl
## # $species: Homo sapiens
## # $rdataclass: EnsDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH53211"]]' 
## 
##              title                             
##   AH53211  | Ensembl 87 EnsDb for Homo Sapiens 
##   AH53715  | Ensembl 88 EnsDb for Homo Sapiens 
##   AH56681  | Ensembl 89 EnsDb for Homo Sapiens 
##   AH57757  | Ensembl 90 EnsDb for Homo Sapiens 
##   AH60773  | Ensembl 91 EnsDb for Homo Sapiens 
##   ...        ...                               
##   AH109606 | Ensembl 109 EnsDb for Homo sapiens
##   AH113665 | Ensembl 110 EnsDb for Homo sapiens
##   AH116291 | Ensembl 111 EnsDb for Homo sapiens
##   AH116860 | Ensembl 112 EnsDb for Homo sapiens
##   AH119325 | Ensembl 113 EnsDb for Homo sapiens

Nota técnica — AnnotationHub vs. biomaRt: biomaRt consulta Ensembl en tiempo real (necesita internet). AnnotationHub descarga recursos y los cachea localmente. Para análisis reproducibles a largo plazo, AnnotationHub es preferible porque puedes fijar la versión exacta del recurso.

# OrgDb: bases de datos de anotación de organismos (mapeo gen → GO, KEGG, etc.)
# Son la fuente que usa org.Hs.eg.db internamente.
query(ah, c("OrgDb", "Homo sapiens"))
## AnnotationHub with 1 record
## # snapshotDate(): 2025-10-29
## # names(): AH121953
## # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
## # $species: Homo sapiens
## # $rdataclass: OrgDb
## # $rdatadateadded: 2025-10-29
## # $title: org.Hs.eg.db.sqlite
## # $description: NCBI gene ID based annotations about Homo sapiens
## # $taxonomyid: 9606
## # $genome: NCBI genomes
## # $sourcetype: NCBI/ensembl
## # $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.ensembl.org/p...
## # $sourcesize: NA
## # $tags: c("NCBI", "Gene", "Annotation") 
## # retrieve record with 'object[["AH121953"]]'
# Chain files para liftOver: permiten convertir coordenadas entre versiones del genoma.
# Esto es crucial porque muchos datasets de Alzheimer más antiguos usan hg19 (GRCh37),
# pero los actuales usan hg38 (GRCh38).
query(ah, c("chainfile", "hg19", "hg38"))
## AnnotationHub with 4 records
## # snapshotDate(): 2025-10-29
## # $dataprovider: UCSC, NCBI
## # $species: Homo sapiens
## # $rdataclass: ChainFile
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH14108"]]' 
## 
##             title                                        
##   AH14108 | hg38ToHg19.over.chain.gz                     
##   AH14150 | hg19ToHg38.over.chain.gz                     
##   AH78915 | Chain file for Homo sapiens rRNA hg19 to hg38
##   AH78916 | Chain file for Homo sapiens rRNA hg38 to hg19
# Para usar un recurso: ah[["AH14150"]]  (reemplaza con el ID del recurso deseado).
# El recurso se descarga y se cachea localmente la primera vez.

6 El flujo completo: visión de conjunto

Hemos recorrido un camino continuo desde las lecturas crudas hasta las asociaciones genéticas — todo con los genes de riesgo de Alzheimer:

Terminal (SRA)          R (biomaRt)           R (GEOquery)
Lecturas crudas    →    Anotación de genes →  Expresión en cerebro
FASTQ inspection        20 genes, 4 IDs       Heatmap + PCA
     │                  + ortólogos ratón           │
     │                       │                      │
     └───────────────────────┼──────────────────────┘
                             │
                             ▼
                    R (clusterProfiler)        R (gwasrapidd)
                    Enriquecimiento GO    →    Asociaciones GWAS
                    Enriquecimiento Reactome   APOE ε4, TREM2 R47H
                    Dotplots + cnetplots       Evidencia genética
Paso Herramienta ¿Qué obtuvimos? Base de datos
1-3 Terminal + SRA Toolkit Lecturas crudas de RNA-seq cerebral SRA (NCBI)
4 biomaRt Anotación, IDs cruzados, ortólogos Ensembl
5 GEOquery Expresión diferencial, heatmap, PCA GEO (NCBI)
6 clusterProfiler + ReactomePA Rutas biológicas convergentes GO, Reactome
7 gwasrapidd Variantes de riesgo y asociaciones GWAS Catalog (EBI)
8 AnnotationHub Recursos centralizados y cacheados Bioconductor

La lección: Ninguna herramienta sola te da la imagen completa. La fuerza está en la integración — cada paso alimenta al siguiente y agrega una capa de evidencia. Esto es exactamente lo que hicimos manualmente en las sesiones 07–08 con TP53 y EGFR, pero ahora de forma programática, escalable y reproducible.

7 Ejercicio 9.1: Tu turno — el pipeline integrador

Crea un script de R (o un RMarkdown) que automatice el siguiente flujo de trabajo para los genes de riesgo de Alzheimer:

  1. Anotación (biomaRt): obtener la anotación completa de los 20 genes (IDs, coordenadas, descripción, ortólogos en ratón).
  2. Expresión (GEOquery): descargar GSE33000, extraer la expresión de los genes de interés, generar un heatmap.
  3. Enriquecimiento (clusterProfiler + ReactomePA): realizar ORA para Gene Ontology (BP) y Reactome. Generar dotplots.
  4. Visualización: producir al menos tres figuras: heatmap, dotplot de GO, y red gen-término (cnetplot).
  5. Interpretación: incluir un párrafo que contextualice los genes en las rutas del Alzheimer (procesamiento de amiloide, neuroinflamación, endocitosis, metabolismo lipídico).

Nota de buenas prácticas — Estructura del script: Organiza tu código en secciones claras con comentarios. Usa sessionInfo() al final. Guarda resultados intermedios (write_csv()) para no tener que re-descargar datos si algo falla.

8 Scripts complementarios

Todo el código de esta sesión está organizado en scripts independientes en la carpeta scripts/:

Script Contenido Bases de datos
script_S09_biomaRt.R Consultas a Ensembl con biomaRt Ensembl
script_S09_GEOquery.R Descarga y exploración de datos de GEO GEO/NCBI
script_S09_enrichment.R Análisis de enriquecimiento GO y Reactome GO, Reactome
script_S09_integracion.R Flujo completo que integra los tres anteriores Todas

9 Preguntas de reflexión

  1. ¿Qué ventajas tiene usar biomaRt sobre descargar archivos manualmente de Ensembl? ¿En qué situación preferirías el acceso manual?
  2. Si el servidor de Ensembl está caído, ¿qué alternativas tienes? (Pista: mirrors, AnnotationHub, archivos locales).
  3. ¿Cuál es la diferencia fundamental entre ORA y GSEA? Si tuvieras un experimento con solo 5 genes diferencialmente expresados, ¿cuál usarías y por qué?
  4. ¿Por qué es importante registrar la versión de las bases de datos en tus scripts de R?
  5. Los genes de riesgo de Alzheimer convergen en unas pocas rutas biológicas. ¿Qué nos dice esta convergencia sobre la arquitectura genética de una enfermedad compleja?

10 Cierre

En esta sesión recorrimos el camino completo desde las lecturas crudas en la terminal hasta el análisis funcional en R — todo con un hilo conductor: los genes de riesgo de Alzheimer. Las mismas bases de datos que exploramos con clics en sesiones anteriores ahora son accesibles con código. Esto no solo ahorra tiempo: hace que nuestros análisis sean reproducibles, escalables y documentados.

El proyecto integrador (disponible como documento aparte en otros/proyecto_integrador.docx) te permitirá llevar esta integración un paso más allá: una plantilla tipo checklist que te guiará por todas las bases de datos que exploramos en las sesiones 07-09, te pedirá escribir código, construir un diagrama de tu proceso de búsqueda y responder preguntas de reflexión de todo el módulo.

11 Información de la sesión

sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.3.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Mexico_City
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] gwasrapidd_0.99.18     AnnotationHub_4.0.0    BiocFileCache_3.0.0   
##  [4] dbplyr_2.5.1           enrichplot_1.30.4      org.Hs.eg.db_3.22.0   
##  [7] AnnotationDbi_1.72.0   IRanges_2.44.0         S4Vectors_0.48.0      
## [10] ReactomePA_1.54.0      clusterProfiler_4.18.4 pheatmap_1.0.13       
## [13] GEOquery_2.78.0        Biobase_2.70.0         BiocGenerics_0.56.0   
## [16] generics_0.1.4         lubridate_1.9.5        forcats_1.0.1         
## [19] stringr_1.6.0          dplyr_1.2.0            purrr_1.2.1           
## [22] readr_2.1.6            tidyr_1.3.2            tibble_3.3.1          
## [25] ggplot2_4.0.2          tidyverse_2.0.0        biomaRt_2.66.0        
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.5.1               urltools_1.7.3.1           
##   [3] ggplotify_0.1.3             filelock_1.0.3             
##   [5] R.oo_1.27.1                 triebeard_0.4.1            
##   [7] polyclip_1.10-7             graph_1.88.1               
##   [9] XML_3.99-0.22               lifecycle_1.0.5            
##  [11] httr2_1.2.2                 lattice_0.22-9             
##  [13] MASS_7.3-65                 magrittr_2.0.4             
##  [15] limma_3.66.0                sass_0.4.10                
##  [17] rmarkdown_2.30              jquerylib_0.1.4            
##  [19] yaml_2.3.12                 otel_0.2.0                 
##  [21] ggtangle_0.1.1              cowplot_1.2.0              
##  [23] DBI_1.2.3                   RColorBrewer_1.1-3         
##  [25] abind_1.4-8                 GenomicRanges_1.62.1       
##  [27] R.utils_2.13.0              ggraph_2.2.2               
##  [29] yulab.utils_0.2.4           tweenr_2.0.3               
##  [31] rappdirs_0.3.4              gdtools_0.5.0              
##  [33] ggrepel_0.9.6               tidytree_0.4.7             
##  [35] rentrez_1.2.4               reactome.db_1.95.0         
##  [37] codetools_0.2-20            DelayedArray_0.36.0        
##  [39] DOSE_4.4.0                  xml2_1.5.2                 
##  [41] ggforce_0.5.0               tidyselect_1.2.1           
##  [43] aplot_0.2.9                 farver_2.1.2               
##  [45] viridis_0.6.5               matrixStats_1.5.0          
##  [47] Seqinfo_1.0.0               jsonlite_2.0.0             
##  [49] tidygraph_1.3.1             systemfonts_1.3.1          
##  [51] tools_4.5.1                 ggnewscale_0.5.2           
##  [53] progress_1.2.3              treeio_1.34.0              
##  [55] Rcpp_1.1.1                  glue_1.8.0                 
##  [57] gridExtra_2.3               SparseArray_1.10.8         
##  [59] xfun_0.56                   qvalue_2.42.0              
##  [61] MatrixGenerics_1.22.0       withr_3.0.2                
##  [63] BiocManager_1.30.27         fastmap_1.2.0              
##  [65] digest_0.6.39               timechange_0.4.0           
##  [67] R6_2.6.1                    gridGraphics_0.5-1         
##  [69] GO.db_3.22.0                RSQLite_2.4.6              
##  [71] R.methodsS3_1.8.2           utf8_1.2.6                 
##  [73] fontLiberation_0.1.0        data.table_1.18.2.1        
##  [75] prettyunits_1.2.0           graphlayouts_1.2.3         
##  [77] httr_1.4.7                  htmlwidgets_1.6.4          
##  [79] S4Arrays_1.10.1             scatterpie_0.2.6           
##  [81] graphite_1.56.0             pkgconfig_2.0.3            
##  [83] gtable_0.3.6                blob_1.3.0                 
##  [85] S7_0.2.1                    XVector_0.50.0             
##  [87] htmltools_0.5.9             fontBitstreamVera_0.1.1    
##  [89] fgsea_1.36.2                scales_1.4.0               
##  [91] png_0.1-8                   ggfun_0.2.0                
##  [93] knitr_1.51                  rstudioapi_0.18.0          
##  [95] tzdb_0.5.0                  reshape2_1.4.5             
##  [97] nlme_3.1-168                curl_7.0.0                 
##  [99] cachem_1.1.0                BiocVersion_3.22.0         
## [101] parallel_4.5.1              pillar_1.11.1              
## [103] grid_4.5.1                  vctrs_0.7.1                
## [105] tidydr_0.0.6                cluster_2.1.8.2            
## [107] evaluate_1.0.5              cli_3.6.5                  
## [109] compiler_4.5.1              rlang_1.1.7                
## [111] crayon_1.5.3                plyr_1.8.9                 
## [113] fs_1.6.6                    ggiraph_0.9.4              
## [115] stringi_1.8.7               viridisLite_0.4.3          
## [117] BiocParallel_1.44.0         assertthat_0.2.1           
## [119] Biostrings_2.78.0           lazyeval_0.2.2             
## [121] GOSemSim_2.36.0             fontquiver_0.2.1           
## [123] Matrix_1.7-4                hms_1.1.4                  
## [125] patchwork_1.3.2             bit64_4.6.0-1              
## [127] KEGGREST_1.50.0             statmod_1.5.1              
## [129] SummarizedExperiment_1.40.0 igraph_2.2.1               
## [131] memoise_2.0.1               bslib_0.10.0               
## [133] ggtree_4.0.4                fastmatch_1.1-8            
## [135] bit_4.6.0                   ape_5.8-1                  
## [137] gson_0.1.0

12 Referencias

12.1 Artículos

  • Davis, S. & Meltzer, P.S. (2007). “GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor.” Bioinformatics, 23(14), 1846–1847. DOI: 10.1093/bioinformatics/btm254
  • Durinck, S. et al. (2009). “Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.” Nature Protocols, 4(8), 1184–1191. DOI: 10.1038/nprot.2009.97
  • Magno, R. & Couto, A. (2020). “gwasrapidd: an R package to query, download and wrangle GWAS Catalog data.” Bioinformatics, 36(2), 649–650. DOI: 10.1093/bioinformatics/btz605
  • Morgan, M. (2023). “AnnotationHub: Client to access AnnotationHub resources.” Bioconductor. Página del paquete
  • Narayanan, M. et al. (2014). “Common dysregulation network in the human prefrontal cortex underlies two neurodegenerative diseases.” Molecular Systems Biology, 10, 743. DOI: 10.15252/msb.20145304
  • Wu, T. et al. (2021). “clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.” The Innovation, 2(3), 100141. DOI: 10.1016/j.xinn.2021.100141
  • Yu, G. & He, Q.Y. (2016). “ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization.” Molecular BioSystems, 12(2), 477–479. DOI: 10.1039/C5MB00663E

12.2 Documentación oficial y tutoriales de los paquetes

Paquete Repositorio Documentación / Vignette
biomaRt Bioconductor Accessing Ensembl (vignette)
GEOquery Bioconductor GEOquery vignette
clusterProfiler Bioconductor Vignette · Libro del autor
ReactomePA Bioconductor ReactomePA vignette
enrichplot Bioconductor enrichplot vignette
org.Hs.eg.db Bioconductor Anotación del genoma humano (Entrez Gene)
gwasrapidd CRAN gwasrapidd vignette
AnnotationHub Bioconductor AnnotationHub vignette
pheatmap CRAN Documentación en CRAN
SRA Toolkit GitHub (NCBI) SRA Toolkit Wiki
Ensembl REST API rest.ensembl.org Documentación de endpoints

13 Errores y sugerencias

Este material está en desarrollo continuo. Si encuentras errores, enlaces rotos o tienes sugerencias para mejorarlo, por favor repórtalos a:

Dra. Yalbi I. Balderas Martínez — Licenciatura en Ciencias Genómicas, UNAM — Marzo 2026
CC BY-NC-SA 4.0 Este material está licenciado bajo CC BY-NC-SA 4.0