Resultados y shrinkage con lfcShrink

r
bioconductor
rnaseq
Por qué los log2FoldChange crudos son ruidosos en genes poco expresados, cómo lfcShrink encoge los LFCs hacia cero según la incertidumbre, qué algoritmo elegir (apeglm, ashr, normal) y por qué los volcano plots publicables casi siempre usan LFCs encogidos.

El problema con LFCs crudos

results() devuelve un LFC máxima verosimilitud (MLE). Para un gen con muchos counts, MLE es buen estimador, ruido bajo, poca incertidumbre. Para un gen con counts bajos, MLE es muy ruidoso: pueden salir LFCs de ±8 simplemente por una muestra con 3 counts vs otra con 0.

Esto da dos patologías visibles:

  1. Top de LFC dominado por basura: si ordenas por abs(log2FoldChange), los primeros suelen ser genes con expresión casi nula y LFCs extremos pero no fiables.
  2. MA plot con cola de puntos extremos a la izquierda (counts bajos), visualmente molesta y sustantivamente engañosa.

La solución estadística: shrinkage. Encoger los LFCs hacia cero proporcionalmente a su incertidumbre. Un gen con counts bajos y LFC nominal alto pero error grande → se encoge mucho hacia cero. Un gen con counts altos y LFC moderado pero error pequeño → no se mueve apenas.

lfcShrink(): el patrón estándar

library(DESeq2)

dds <- DESeq(dds)
resultsNames(dds)
# [1] "Intercept" "condicion_tto_vs_ctrl"

res_shr <- lfcShrink(dds, coef = "condicion_tto_vs_ctrl", type = "apeglm")

Argumentos:

  • coef: qué coeficiente encoger (el nombre del resultsNames()).
  • type: qué algoritmo de shrinkage usar.

Output: igual estructura que results(), pero con LFCs encogidos. Los p-valores no se tocan, siguen siendo los del test Wald sobre los LFCs crudos.

Los tres algoritmos: apeglm, ashr, normal

Tres métodos disponibles, en orden de preferencia para uso general:

apeglm (recomendado por defecto)

res_shr <- lfcShrink(dds, coef = "condicion_tto_vs_ctrl", type = "apeglm")

Implementación del prior adaptativo (Zhu, Ibrahim, Love 2018). Es el default recomendado por la vignette de DESeq2 desde 2018. Buen balance: encoge LFCs ruidosos sin perjudicar genes con buena señal.

Limitación: solo funciona con coef = ..., no con contrast = .... Para contrastes arbitrarios (entre dos niveles no-referencia), usa ashr.

ashr

res_shr <- lfcShrink(dds, contrast = c("tratamiento", "dosis2", "dosis1"),
                     type = "ashr")

Adaptive Shrinkage (Stephens 2017). Más flexible, admite contrast = c(...). Buena alternativa cuando apeglm no es viable.

Requiere el paquete ashr instalado: install.packages("ashr").

normal (legacy)

res_shr <- lfcShrink(dds, coef = "condicion_tto_vs_ctrl", type = "normal")

El método clásico de DESeq2 (anterior a 2018). Aún funciona pero apeglm lo supera en casi todo. Mencionado solo porque verás tutoriales antiguos que lo usan.

Método Cuándo Notas
apeglm Por defecto Solo con coef
ashr Si necesitas contrast=c(...) arbitrario Requiere paquete
normal Solo por compatibilidad con análisis viejos Legacy

Cuándo aplicar shrinkage

Tres casos donde sí:

  • Ranking de genes por LFC (top up/down regulated).
  • Volcano plot publicable: el eje X queda mucho más limpio.
  • Heatmap de top genes: si seleccionas por LFC, mejor con LFCs encogidos.

Un caso donde no importa:

  • Filtrado por padj: usa la res sin encoger. Los p-valores son los mismos. Es solo si después vas a ordenar por LFC o pintar volcano cuando el shrinkage cambia algo.

Recomendación pragmática: calcula ambos, usa los crudos para inspección y los encogidos para visualización publicable.

res_crudo <- results(dds, alpha = 0.05)
res_shr   <- lfcShrink(dds, coef = "condicion_tto_vs_ctrl", type = "apeglm")

Comparar visualmente: MA plots

par(mfrow = c(1, 2))
plotMA(res_crudo, main = "Sin shrinkage", ylim = c(-5, 5))
plotMA(res_shr,   main = "apeglm",       ylim = c(-5, 5))

plotMA (que viene con DESeq2): media de counts en X, LFC en Y, puntos rojos = padj < alpha. En el plot sin shrinkage verás “explosión” de LFCs a la izquierda (counts bajos). En el plot con apeglm, la cola se calma y solo quedan los LFCs robustos.

Diferencias prácticas

Con apeglm:

  • LFCs de genes muy expresados con señal clara → casi no cambian.
  • LFCs de genes con pocos counts o señal débil → se encogen hacia cero, a veces mucho.
  • p-valores: idénticos a los de results(). Shrinkage no es un nuevo test.
  • padj: idénticos. Mismo orden de genes significativos.

Lo que cambia es el LFC, lo que afecta:

  • Ranking por log2FoldChange: distinto.
  • Filtro abs(LFC) > 1: más estricto con LFCs encogidos (filtra basura).
  • Volcano y heatmap visuales: más limpios.

Inspeccionar los cambios

plot(res_crudo$log2FoldChange, res_shr$log2FoldChange,
     pch = 16, cex = 0.4,
     xlab = "LFC sin shrinkage", ylab = "LFC apeglm")
abline(0, 1, col = "red", lty = 2)

Lo esperable: nube de puntos cerca de la diagonal para LFCs cercanos a 0, y desviación creciente para LFCs extremos (esos son los que se encogen). El plot del LFC encogido está visiblemente menos disperso en la cola.

El flujo idiomático actualizado

Lo que se incluye en un análisis profesional moderno:

library(DESeq2)

# Filtro mínimo
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

# Modelo
dds <- DESeq(dds)

# Resultados crudos para tests (p-valores)
res_crudo <- results(dds, alpha = 0.05)
summary(res_crudo)

# Resultados con shrinkage para visualización y ranking
res <- lfcShrink(dds, coef = "condicion_tto_vs_ctrl", type = "apeglm")

# Diagnóstico
plotDispEsts(dds)
plotMA(res)

# Filtros y export
res_sig <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]
nrow(res_sig)

library(tibble)
as.data.frame(res) |>
    rownames_to_column("gen_id") |>
    write.csv("DE_results.csv", row.names = FALSE)

Las decisiones que tomas:

  • alpha = 0.05 (en lugar del default 0.1), más estricto pero estándar en publicación.
  • apeglm para shrinkage, default moderno.
  • Filtro abs(LFC) > 1: al menos 2× cambio. Razonable para biología.
  • Export con IDs como columna: para análisis downstream y anotación.

Trampas habituales

  • Usar LFCs crudos para gráficos publicables. Los volcano plots con LFCs sin encoger son estéticamente los típicos de la vignette antigua. Hoy se prefieren con shrinkage.
  • Confundir “shrinkage” con “test”. El shrinkage no cambia p-valores. Si esperabas que apareciesen menos significativos con shrinkage, no es así, el filtrado de significancia es independiente.
  • apeglm con contrast = c(...). No funciona. Usa type = "ashr" si necesitas contrastes que no son coeficientes simples.
  • No comparar plotMA antes y después. El antes y después de shrinkage es la mejor herramienta de diagnóstico. Compara siempre.
  • Reportar LFCs sin aclarar si son encogidos. En un paper o informe, especifica: “log2 fold-changes shrunken using apeglm”. Las dos versiones (encogida y no) dan números distintos para los mismos genes.

En la siguiente entrega

Cerramos el bloque de análisis. La siguiente parte es calidad y visualización: PCA con expectativa, heatmaps de muestras, MA y volcano publicables, y la lectura honesta de un PCA. Lo siguiente.