Comprensión completa de GSEA y fgsea: desde la puntuación de enriquecimiento hasta la práctica con MSigDB.

GSEA y fgsea: Cómo demostrar que no es un solo gen, sino toda una vía metabólica la que está cambiando.

Cómo transformar una lista de genes con expresión diferencial en una descripción de vías. Analice paso a paso cómo se calcula la puntuación de enriquecimiento acumulado (running enrichment score) en GSEA y aplique las vías Hallmark de MSigDB con fgsea.

Intermedio

20min

Verificado (2026-07-24)

GSEApathway enrichmentMSigDBfgsea

Se han trazado 500 puntos en el gráfico de dispersión. ¿Qué ha ocurrido?

Al analizar los datos de S19 y S20, obtuvimos una lista de genes significativos. Sin embargo, los revisores, los colaboradores e incluso ustedes se preguntan: "¿Qué información aportan estos 500 genes?". No existe ningún artículo que simplemente enumere 500 nombres de genes. Los artículos suelen indicar que "la vía de la proliferación celular está activada y la vía de la respuesta inmune está inhibida en el tumor".

Cambiar la perspectiva de los genes individuales a las vías (pathways) es el último paso en el proceso de transformación de datos "de FASTQ a artículo" y el objetivo de esta sección. La clave está en diferenciar entre dos enfoques estadísticos: ORA y GSEA.

Primera idea: simplemente contemos (ORA)

El método más intuitivo es el análisis de sobre-representación (Over-Representation Analysis). Se utiliza la distribución hipergeométrica para responder a la pregunta: "¿Cuántos de mis 500 genes significativos pertenecen a la vía de la proliferación celular? ¿Es un número mayor del esperado por azar?".

P(X \ge k) = \sum_{i=k}^{n} \frac{\binom{K}{i}\binom{M-K}{n-i}}{\binom{M}{n}}

Donde $M$ es el número total de genes, $K$ es el número de genes de la vía de la proliferación celular, $n$ es el número de genes en mi lista significativa y $k$ es el número de genes de la vía de la proliferación celular dentro de esa lista. Es un cálculo similar a la prueba exacta de Fisher.

Aunque ORA es sencillo, tiene una limitación importante: solo considera los genes que superan un umbral (padj < 0.05) y descarta el resto. Si todas las vías muestran un aumento ligero pero consistente (por ejemplo, un 1,3 veces), pero ningún gen individual supera el umbral, ORA no detectará esa vía. En biología, la señal más interesante suele ser precisamente este "cambio pequeño pero consistente".

Segunda idea: analicemos todo el rango (GSEA)

El análisis de enriquecimiento de conjuntos de genes (Gene Set Enrichment Analysis) elimina los umbrales. Ordena los 20.000 genes en una sola lista según su cambio de expresión (o estadística con signo) y analiza si los genes de una vía específica están concentrados en la parte superior de la lista.

La herramienta clave es la puntuación de enriquecimiento acumulada (running enrichment score). Al recorrer la lista ordenada de arriba a abajo:

Si se encuentra un gen perteneciente a la vía, se aumenta la puntuación (el peso es proporcional a la magnitud de la estadística de ese gen),
Si se encuentra un gen que no pertenece a la vía, se disminuye la puntuación.

ES = \max_i \left( \sum_{g_j \in S, j \le i} \frac{|r_j|^p}{N_R} - \sum_{g_j \notin S, j \le i} \frac{1}{N - N_S} \right)

Si los genes de la ruta están agrupados en la parte superior de la lista, la suma acumulada aumenta drásticamente al principio y el ES (su valor máximo) es grande. Si están dispersos aleatoriamente, los aumentos y las disminuciones se compensan, y el ES se mantiene cerca de 0. Esta curva tiene exactamente la misma estructura que la estadística de Kolmogorov-Smirnov: muestra la desviación máxima entre dos distribuciones.

Seguimiento manual de la puntuación acumulada

Supongamos que hemos ordenado 8 genes por rango y que el conjunto de genes de interés es $S = \{g_1, g_2, g_5\}$ (para simplificar los pesos, asumimos que un acierto es +0.33 y un fallo es −0.20).

Rango	Gen	¿Elemento de la ruta?	Cambio	Suma acumulada
1	g1	● Acierto	+0.33	0.33
2	g2	● Acierto	+0.33	0.66
3	g3	Fallo	−0.20	0.46
4	g4	Fallo	−0.20	0.26
5	g5	● Acierto	+0.33	0.59
6	g6	Fallo	−0.20	0.39
7	g7	Fallo	−0.20	0.19
8	g8	Fallo	−0.20	−0.01

El valor máximo de la suma acumulada es ES = 0.66 (en el rango 2). Dado que los genes de la ruta están agrupados en la parte superior, se forma un pico al principio. La altura de este pico es la estadística, y la ubicación del pico indica "hasta qué gen hay contribuyentes clave" (borde principal).

Significación: permutación

Para determinar si el ES ocurre por casualidad, se necesita una distribución nula. GSEA mezcla aleatoriamente las etiquetas de los genes (o las etiquetas fenotípicas), recalcula el ES miles de veces y obtiene el valor p según la posición del ES real dentro de esa distribución. Dado que el tamaño varía entre las rutas, se utiliza el NES (ES normalizado) para comparar entre las rutas.

Aquí radica el problema de las versiones antiguas de GSEA: miles de permutaciones × miles de rutas = varias horas. fgsea (GSEA rápido) acelera este cálculo cientos de veces mediante un método de Monte Carlo multinivel adaptativo, estimando con precisión valores p extremadamente pequeños. Es el estándar actual en la práctica.

Ejecutar MSigDB Hallmark con fgsea

Ahora, hagamos la práctica. Tomaremos directamente los resultados de DESeq2 de S19 y continuaremos hasta el análisis de rutas.

Paso 1: Crear el vector de rangos

if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("fgsea", "msigdbr"))
library(fgsea); library(msigdbr)

# res = resultado DESeq2 de S19 (incluye columnas log2FoldChange, stat)
# Generar vector de rangos con estadística firmada (eliminar NA, nombre de gen = SYMBOL)
ranks <- res$stat
names(ranks) <- rownames(res)
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)

Es más estable utilizar stat (estadístico de Wald) basado en el orden que utilizar solo el cambio de magnitud, ya que refleja tanto la magnitud como la fiabilidad.

Paso 2: Cargar y ejecutar conjuntos de vías

# 50 vías Hallmark de MSigDB (humano)
h <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(h$gene_symbol, h$gs_name)

set.seed(42)
fgseaRes <- fgsea(pathways = pathways, stats = ranks,
                  minSize = 15, maxSize = 500)
head(fgseaRes[order(padj)], 10)

Paso 3: Diagrama de la ruta principal

topUp <- fgseaRes[NES > 0][order(padj)][1, pathway]
plotEnrichment(pathways[[topUp]], ranks) +
  ggplot2::labs(title = topUp)

La curva que dibuja plotEnrichment es precisamente la puntuación acumulada (running score) que hemos seguido manualmente. Si el pico se encuentra a la izquierda (lado de alta expresión), esa vía está activada en el tumor. En los datos de cáncer de mama, en la parte superior aparecerán vías como las del ciclo celular HALLMARK_E2F_TARGETS y HALLMARK_G2M_CHECKPOINT. BRCA1 y MKI67, que vimos en S19, son miembros de esa vía. La explicación queda completa.

Mapeo CS

Estadístico de Kolmogorov-Smirnov: ES es la desviación máxima entre la distribución acumulada de los genes de la vía y la distribución uniforme. Es la aplicación biológica de la prueba de KS.
Suma acumulada / máximo subconjunto parcial: La curva de puntuación de enriquecimiento se obtiene buscando el valor máximo de la suma acumulada, lo que la convierte en un problema similar al de encontrar el subarreglo de máximo subconjunto parcial, relacionado con el algoritmo de Kadane.
Prueba de permutación: Crear la distribución mediante el intercambio de etiquetas en lugar de una distribución nula analítica es la esencia de la estadística no paramétrica; el enfoque multinivel de fgsea es una aceleración de Monte Carlo para la estimación de eventos raros.

Errores comunes

Incompatibilidad de identificadores de genes: si el vector de rangos utiliza SYMBOL y el conjunto de vías utiliza ENTREZ, la coincidencia será 0. No ignore la advertencia minSize.
Clasificación basada únicamente en el cambio de expresión (fold change): los valores extremos de LFC de los genes con baja expresión contaminan la parte superior. Se recomienda utilizar stat (LFC/SE ajustado).
Confusión entre ORA y GSEA: clusterProfiler enrichKEGG, que "incluye solo genes significativos", es ORA. El análisis de rango completo sin umbral corresponde a GSEA/fgsea. Elija según el objetivo.
Falta de corrección por pruebas múltiples: al realizar pruebas con miles de vías, se decide mediante padj (BH). Si se informa con pval, se obtendrán muchas vías falsas positivas.

Para profundizar

El texto es una narración reconstruida directamente por BPD. Para un estudio más profundo, utilice los siguientes materiales de referencia.

Harvard STAT115: Pathway and Gene Set Enrichment Analysis del profesor Xiaole Shirley Liu (con subtítulos completos y una excelente traducción automática). La distinción estadística entre ORA y GSEA es clara.
Artículo original (GSEA): Subramanian et al. (2005), Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles, PNAS 102:15545. Es la fuente de la puntuación acumulada.
Artículo original (fgsea): Korotkevich, Sukhov, Sergushichev (2019), Análisis rápido de enriquecimiento de conjuntos de genes, bioRxiv 060012. Implementa una aceleración multinivel.
Base de datos: MSigDB (gsea-msigdb.org) — Ofrece gratuitamente las colecciones Hallmark (H), C2 (vías canónicas de KEGG/Reactome) y C5 (GO). La licencia es del tipo CC para uso en investigación, por lo que es necesario verificar las condiciones para su uso comercial.

De este modo, se completa el proceso "de FASTQ a artículo", en el que los datos FASTQ sin procesar se transforman en figuras y descripciones de vías. A partir del siguiente capítulo, S22, pasaremos a un escenario completamente diferente: el ensamblaje de novo, que consiste en ensamblar secuencias a partir de cero sin un genoma de referencia.

BioPlayground