Ya aprendí DESeq2, ¿por qué debo aprender otra herramienta?

En S19, creamos un gráfico de volcanes con DESeq2. Sin embargo, al enviar un artículo, los revisores a menudo comentan: "Verifique si los resultados se reproducen con edgeR" o, si se tiene una cohorte grande con decenas de muestras y DESeq2 es demasiado lento, lo que obliga a trabajar toda la noche.

Las tres herramientas — DESeq2, edgeR y limma-voom — no son competidoras, sino tres filosofías estadísticas diferentes para el mismo problema. Al finalizar esta sección, sabrá "cuál usar en cada situación". Además, verá que todas comparten la misma idea central: Bayes empírico. Lo vimos una vez en S19 y lo veremos dos veces más aquí. Si se repite tres veces, no es coincidencia, sino un principio.

edgeR: el hermano de DESeq2

Al igual que DESeq2, edgeR modela los conteos con una distribución binomial negativa. Las diferencias radican en dos aspectos:

Primero, el método de normalización. DESeq2 usa el factor de tamaño de la mediana de las razones, mientras que edgeR usa TMM (media recortada de los valores M). La idea de TMM es calcular la diferencia logarítmica de las expresiones (valor M) y la expresión media (valor A) entre dos muestras, luego recortar los extremos y usar el promedio ponderado del centro restante para ajustar la escala de la biblioteca.

donde $M_g = \log_2(y_{gk}/N_k) - \log_2(y_{gr}/N_r)$ y $G^*$ es el conjunto de genes después de eliminar los valores extremos. El supuesto clave es que "la mayoría de los genes no cambian su expresión". Este recorte evita el sesgo de composición, donde unos pocos genes altamente expresados (por ejemplo, la globina) dominan la biblioteca y hacen que otros genes parezcan artificialmente menos expresados.

Segundo, la estimación de la dispersión. edgeR reduce $\alpha$ en tres etapas: común → con tendencia → específica del gen. En la etapa específica del gen, la dispersión de cada gen se acerca a la tendencia común; esto es exactamente lo mismo que el ajuste de DESeq2 y corresponde al Bayes empírico. Solo cambia el nombre.

limma-voom: un enfoque completamente diferente

limma fue originalmente una herramienta para la era de los microarrays (valores continuos). Es un motor estadístico refinado durante 20 años bajo la suposición de una distribución normal y modelos lineales. Dado que los conteos de RNA-seq no siguen una distribución normal, voom cierra esta brecha.

La idea de voom es sorprendentemente práctica. Al convertir los recuentos en log-CPM, los genes con una media baja tienden a tener una mayor varianza (tendencia media-varianza). Voom ajusta esta tendencia mediante una curva suave y asigna un peso de precisión a cada observación.

Es decir, "reduce la confianza en las observaciones cuya varianza se predice que será alta". Las observaciones inestables de genes poco expresados contribuyen menos al ajuste de regresión. Al ajustar un modelo lineal mediante mínimos cuadrados ponderados con este peso, podemos llevar directamente los datos de recuento a las herramientas potentes del mundo de la distribución normal.

t moderado: Bayes empírico de limma

Aquí aparece el tercer Bayes empírico. La prueba t estándar utiliza $s_g$ (el error estándar por gen) para cada gen individualmente. Con solo tres réplicas, $s_g$ fluctúa mucho. Limma construye la estadística t utilizando un valor para la varianza por gen $s_g^2$ que ha sido atraído hacia la varianza previa global $s_0^2$.

$\tilde{s}_g^2$ es un promedio ponderado entre la varianza previa y la varianza observada. Cuando hay pocas réplicas y $d_g$ es pequeño, se atrae fuertemente hacia $s_0^2$. Es exactamente el mismo esquema que el ajuste de dispersión de S19. La única diferencia es si el objetivo es la dispersión o la varianza residual.

Para comprenderlo intuitivamente

Supongamos que los grados de libertad previos son $d_0 = 4$ y la varianza previa es $s_0^2 = 0,10$, y que la varianza observada de un gen determinado es $s_g^2 = 0,90$ (con 3 réplicas, por lo que $d_g = 2$).

La varianza observada de 0,90 se ha reducido a 0,367. Es probable que este gen haya tenido una "alta dispersión por casualidad", y la varianza moderada atenúa esa exageración. Este es el principio detrás de la reducción de falsos positivos.

¿Cuándo usar qué herramienta?

Los criterios prácticos para elegir entre las tres metodologías se resumen a continuación.

Recordemos una regla empírica: si hay pocos repetidos, DESeq2/edgeR; si hay muchas muestras, limma-voom. Y si las tres herramientas difieren significativamente en los genes más significativos con expresión diferencial, eso no es un problema de la herramienta, sino una señal de que hay efectos de lote u valores atípicos ocultos en los datos.

Comparación directa de los tres métodos con R en Colab

Ahora, la práctica. Ejecutemos las tres herramientas en el mismo entorno con el conjunto de datos airway (tratamiento farmacológico frente a control en asma, 8 muestras).

Paso 1: Datos y paquetes

if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("airway", "edgeR", "limma", "DESeq2"))

library(airway); data(airway)
counts <- assay(airway)                 # gene × sample raw count
group  <- airway$dex                    # trt vs untrt
counts <- counts[rowSums(counts) > 10, ] # remove low expression

Paso 2: edgeR

library(edgeR)
y <- DGEList(counts = counts, group = group)
y <- calcNormFactors(y)                 # TMM normalization
design <- model.matrix(~ group)
y <- estimateDisp(y, design)            # common→trended→tagwise
fit <- glmQLFit(y, design)              # quasi-likelihood
qlf <- glmQLFTest(fit, coef = 2)
edger_res <- topTags(qlf, n = Inf)$table

Paso 3: limma-voom

library(limma)
v <- voom(y, design, plot = TRUE)       # calculate precision weight (check curve)
vfit <- lmFit(v, design)
vfit <- eBayes(vfit)                    # moderated t
voom_res <- topTable(vfit, coef = 2, number = Inf)

voom(plot = TRUE) dibuja la curva de media-varianza, que debe verificarse visualmente. Si esta curva es plana, significa que la normalización ya se realizó correctamente; si está muy curvada, indica que los pesos de Voom tienen un papel importante.

Paso 4: Superponer tres conjuntos de genes significativos

sig <- function(tab, p, lfc) rownames(tab)[tab[[p]] < 0.05 & abs(tab[[lfc]]) > 1]
e <- sig(edger_res, "FDR", "logFC")
v_ <- sig(voom_res, "adj.P.Val", "logFC")
length(intersect(e, v_)) / length(union(e, v_))  # Jaccard similarity

Por lo general, se obtienen valores superiores a 0.8. Las tres herramientas suelen llegar a conclusiones similares, pero difieren únicamente en los genes cercanos al límite. Estos genes límite son la prioridad para la validación posterior (qPCR).

Mapeo de CS

• Empirical Bayes = regularización: Atraer las estimaciones individuales hacia un objetivo común constituye el esqueleto de la regularización de CS (L2, ridge). La dispersión de S19 y la varianza residual —solo difieren en el objeto de estudio, pero son exactamente las mismas matemáticas—.
• precision weight = mínimos cuadrados ponderados (WLS): voom introduce la confianza por observación como un peso inverso a la varianza. Es la técnica estándar en estadística e interpretación numérica para manejar la heterocedasticidad.
• TMM trimming = estadísticas robustas: Recortar los valores extremos y utilizar solo la parte central es el principio de las estadísticas robustas, como la media recortada o la mediana.

Defectos comunes

• Realizar pruebas sin aplicar factores de normalización → El sesgo de composición distorsiona completamente los resultados. calcNormFactors() no es opcional.
• Introducir valores ya log/CPM en limma-voom → voom espera datos en bruto (DGEList). La doble transformación arruina la curva.
• Insistir en limma-voom con solo 3 muestras → La aproximación normal se debilita. Para pocas réplicas, DESeq2/edgeR son más seguros.
• Los resultados de las tres herramientas difieren mucho y se culpa a la "herramienta" → Generalmente se trata de un efecto de lote oculto. Corrijámoslo primero con design = ~ batch + group.

Para profundizar más

El texto es una narración reconstruida directamente por BPD. Para profundizar, utilicemos los siguientes materiales canónicos en orden.

• Harvard STAT115 — La clase de Xiaole Shirley Liu Differential Expression: edgeR & limma (subtítulos completos, buena traducción automática). Resume estadísticamente las diferencias entre los modelos de las tres herramientas.
• Artículo original (voom): Law, Chen, Shi, Smyth (2014), voom: precision weights unlock linear model analysis tools for RNA-seq read counts, Genome Biology 15:R29. Es la fuente original de la derivación del peso de precisión.
• Artículo original (edgeR): Robinson, McCarthy, Smyth (2010), edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26:139.
• Artículo original (t moderado): Smyth (2004), Linear models and empirical Bayes methods for assessing differential expression in microarray experiments, SAGMB 3:Artículo 3. Es un clásico de los métodos de Bayes empírico.
• Documentación oficial: Guía del usuario de limma (bioconductor.org/packages/limma) — El capítulo de RNA-seq es la referencia para las recetas prácticas.

Ahora que tenemos estas tres herramientas, en la siguiente entrega, S21, introduciremos esta lista de genes significativos en un análisis GSEA para ampliar el enfoque de "genes individuales" a "vías biológicas". En ese momento, se completará la última parte del proceso "FASTQ to Paper".