DESeq2: el flujo básico

r
bioconductor
rnaseq
El ciclo central de DESeq2 en tres llamadas: DESeq(), results() y resultsNames(). Qué hace cada paso por dentro, size factors, dispersiones, modelo NB, y por qué la interpretación de un log2FoldChange exige saber qué se comparó.

Tres líneas

El núcleo de DESeq2 son tres líneas. El resto son matices.

library(DESeq2)

dds <- DESeq(dds)
res <- results(dds)
summary(res)

DESeq() ajusta el modelo. results() extrae el contraste por defecto. summary() da el resumen.

Para entender lo que hay debajo, conviene ver cada paso por separado.

Qué hace DESeq() por dentro

DESeq() es un atajo que ejecuta tres pasos:

dds <- estimateSizeFactors(dds)        # 1. tamaño de librería ajustado
dds <- estimateDispersions(dds)        # 2. dispersión gen a gen
dds <- nbinomWaldTest(dds)             # 3. test Wald sobre el modelo NB

1. estimateSizeFactors(): normalización entre muestras

Cada muestra se secuenció con profundidad distinta, una con 30M reads, otra con 50M. Comparar counts crudos da diferencias artificiales. DESeq2 calcula un size factor por muestra (un escalar) y divide los counts por él internamente.

El método (la mediana de ratios de Anders & Huber 2010): es robusto a unos pocos genes con expresión muy alta. Mucho mejor que dividir simplemente por el total (colSums).

sizeFactors(dds)
# untrt1    trt1   untrt2    trt2   ...
# 1.034   0.892   1.121   0.987    ...

Valores cercanos a 1 = profundidades similares. Si una muestra tiene size factor 0.3, está claramente subsecuenciada, probable problema de QC, conviene investigar.

2. estimateDispersions(): la dispersión por gen

Para counts de RNA-seq, la distribución estadística adecuada es la binomial negativa, Poisson con sobre-dispersión. Cada gen tiene su propia dispersión, pero con N pequeño (3-6 muestras por grupo) estimar la dispersión por separado para cada gen es ruidoso.

DESeq2 resuelve esto con shrinkage: estima la dispersión cruda por gen, ajusta una curva global como referencia, y encoge las estimaciones individuales hacia esa curva. Resultado: dispersiones más estables, especialmente para genes con counts bajos.

plotDispEsts(dds)

El gráfico estándar de diagnóstico: muestra puntos crudos (rojo), ajuste global (negro) y estimaciones encogidas (azul). Lo esperable: la curva baja monotónicamente con la media, más expresión → menos dispersión relativa. Si no se ve así, algo va mal (batch fuerte sin modelar, contaminación, datos mezclados).

3. nbinomWaldTest(): el test estadístico

Para cada gen y cada coeficiente del modelo, se calcula:

  • log2FoldChange (LFC): la diferencia estimada en log2.
  • lfcSE: error estándar del LFC.
  • stat: el estadístico Wald = LFC / lfcSE.
  • pvalue: p-valor del test.

Después, padj = p-valor ajustado por múltiples comparaciones (Benjamini-Hochberg).

results(): extraer el contraste

res <- results(dds)
head(res)
# log2 fold change (MLE): condicion tto vs ctrl
# Wald test p-value: condicion tto vs ctrl
#
#         baseMean log2FoldChange     lfcSE      stat   pvalue     padj
# ENSG...  142.32          0.512    0.123     4.163  3.1e-05  0.00034
# ENSG...   89.14         -0.230    0.198    -1.161  0.2456   0.4912

Columnas:

  • baseMean: media (normalizada) de counts del gen entre todas las muestras. Solo informativo, no del contraste.
  • log2FoldChange: el cambio estimado en log2. Un LFC de 1 = doble. LFC de -1 = mitad.
  • lfcSE: error estándar del LFC. Útil para evaluar incertidumbre.
  • stat: el estadístico Wald.
  • pvalue: p-valor crudo.
  • padj: p-valor ajustado (FDR). Es el que usas para filtrar significantes.

La cabecera (log2 fold change (MLE): condicion tto vs ctrl) te dice qué se comparó. Es la información crítica. Si dice “tto vs ctrl”, un LFC positivo significa subida en tto respecto a ctrl.

resultsNames(): qué coeficientes hay

resultsNames(dds)
# [1] "Intercept"             "condicion_tto_vs_ctrl"

El intercepto + un coeficiente por condición (relativo al nivel de referencia). Para extraer un contraste específico:

res <- results(dds, name = "condicion_tto_vs_ctrl")

Por defecto, results() devuelve el último coeficiente del último factor de la fórmula. Si tu fórmula es ~ batch + condicion, el default es la comparación de condición. Es razonable, pero sé explícito con name = o contrast = si tu diseño es complejo.

Resumen del resultado

summary(res)
# out of 25000 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up)      : 1250, 5.0%
# LFC < 0 (down)    : 1180, 4.7%
# outliers          : 0,    0%
# low counts        : 4500, 18%
# (mean count < 5)

Cuántos genes pasan el umbral (FDR < 0.1 por defecto). El independent filtering quita genes con baja media de counts antes del ajuste de p-valor, sube la potencia para los que sí tienen señal.

Para usar FDR < 0.05 (más estricto):

res <- results(dds, alpha = 0.05)
summary(res)

Filtrar significantes

El patrón típico:

res_sig <- res[!is.na(res$padj) & res$padj < 0.05, ]
nrow(res_sig)
# 1280

!is.na(res$padj) es importante: el independent filtering deja NA en padj para genes filtrados. Sin la comprobación, padj < 0.05 con un NA da NA y se cuela.

Filtrar también por magnitud del LFC, no solo por significación:

res_fuerte <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]

LFC > 1 = al menos doble. Para análisis biológico, suele ser más útil que solo padj, un cambio del 5 % puede ser “significativo” con N alto y a la vez irrelevante en biología.

Ordenar y exportar

res_ord <- res[order(res$padj), ]    # más significativos arriba

library(tibble)
res_df <- as.data.frame(res_ord) |>
    rownames_to_column("gen_id")

write.csv(res_df, "resultados_DE.csv", row.names = FALSE)

results() devuelve un DESeqResults (extensión de DataFrame). Conviértelo a data.frame para write.csv, dplyr, etc.

El flujo idiomático completo

library(DESeq2)

# Cargar (ya cubierto en tutorial 4)
dds <- DESeqDataSetFromMatrix(counts, meta, design = ~ batch + condicion)

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

# Análisis
dds <- DESeq(dds)

# Resultados del contraste principal
res <- results(dds, alpha = 0.05)
summary(res)

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

# Significantes
sig <- res[!is.na(res$padj) & res$padj < 0.05, ]

Diez líneas para un análisis estándar. El resto son refinamientos (shrinkage, anotación, visualización publicable) que veremos en los próximos tutoriales.

Trampas habituales

  • Confundir LFC con fold-change normal. LFC = log2(FC). Para “veces de cambio” reales, 2^LFC (subida) o 2^|LFC| (en valor absoluto). Un LFC = 2 significa 4× más, no 2×.
  • Confiar en pvalue en vez de padj. Con 25.000 genes testeados, sin corrección saldrán ~1.250 falsos positivos por azar. Usa siempre padj.
  • Filtrar por padj < 0.05 sin manejar NAs. Los NA se filtran como FALSE en < 0.05, lo que es deseado, pero conviene escribirlo explícito: !is.na(padj) & padj < 0.05.
  • Ignorar plotDispEsts. Si la curva de dispersión es rara (puntos crudos muy por encima de la curva, o sin tendencia con la media), hay problema. Revisar antes de seguir.
  • Asumir el orden de la comparación. La cabecera de results() te dice “tto vs ctrl” o “ctrl vs tto” según el orden de factor. Léelo siempre antes de interpretar signos del LFC.

En la siguiente entrega

Ya tienes resultados. La siguiente entrega es transformaciones, vst y rlog, para visualización, PCA y heatmaps. Counts crudos no sirven para esos análisis exploratorios. Los logs ingenuos tampoco. DESeq2 tiene dos transformaciones diseñadas precisamente para esto. Lo siguiente.