Diseño experimental y contrastes

El alma estadística de DESeq2. Cómo se interpretan las fórmulas, qué cambia con el orden, cómo ajustar por batch, qué significa una interacción y cómo se especifican contrastes arbitrarios con contrast=c(...).

El núcleo: la fórmula

En DESeq2, la fórmula del design no es decorativa, define qué se modela y, por consecuencia, qué contrastes son posibles y qué se ajusta. Es la decisión estadística más importante del análisis.

dds <- DESeqDataSet(se, design = ~ condicion)              # solo condición
dds <- DESeqDataSet(se, design = ~ batch + condicion)      # condición ajustada por batch
dds <- DESeqDataSet(se, design = ~ sex + age + condicion)  # varias covariables
dds <- DESeqDataSet(se, design = ~ condicion * tiempo)     # interacción

Sintaxis idéntica a lm o glm en R base. Convención DESeq2: la variable de interés va al final. Las anteriores se tratan como covariables a ajustar (su efecto se controla, pero no es el foco del contraste).

Diseños comunes

Una variable, dos grupos

design(dds) <- ~ condicion

El caso más simple. results() por defecto compara los dos niveles de condicion.

levels(dds$condicion)
# [1] "ctrl" "tto"
res <- results(dds)
# log2 fold change (MLE): condicion tto vs ctrl

LFC positivo = sube en tto (porque ctrl es la referencia, primer nivel del factor).

Una variable, tres o más grupos

levels(dds$tratamiento)
# [1] "ctrl" "dosis1" "dosis2"

res_d1 <- results(dds, contrast = c("tratamiento", "dosis1", "ctrl"))
res_d2 <- results(dds, contrast = c("tratamiento", "dosis2", "ctrl"))
res_d12 <- results(dds, contrast = c("tratamiento", "dosis2", "dosis1"))

contrast = c("factor", "numerador", "denominador"), sintaxis explícita. LFC positivo significa subida del numerador frente al denominador. Esta forma es preferible a name = cuando hay tres o más niveles, porque controlas el contraste exacto sin pensar en cómo se llaman los coeficientes internos.

Ajustar por batch técnico

design(dds) <- ~ batch + condicion

Si las muestras se procesaron en lotes técnicos distintos (días de secuenciación, técnicos diferentes, kits distintos), el batch es una fuente de varianza no biológica que hay que ajustar para no confundirlo con efecto del tratamiento.

~ batch + condicion significa: “modela el efecto medio del batch, y sobre lo que queda, mide el efecto de la condición”. Es lo que en estadística llaman adjustment for confounding.

El batch debe ser una variable observada, no puedes ajustar si no sabes qué muestra fue de qué lote. Y debe ser no totalmente confundido con la condición: si todas las muestras de tto fueron al batch A y todas las de ctrl al batch B, no se puede separar estadísticamente.

Varias covariables

design(dds) <- ~ sex + age + batch + condicion

Cada + añade un término. El orden importa menos matemáticamente que en regresión clásica (DESeq2 hace inferencia tipo III), pero por convención: covariables primero, variable de interés al final.

age puede ser numérica continua, no tiene por qué ser factor.

Interacciones

Cuando el efecto de A depende del nivel de B:

design(dds) <- ~ condicion + tiempo + condicion:tiempo
# o equivalente
design(dds) <- ~ condicion * tiempo

Si la respuesta a tto cambia entre tiempo_1 y tiempo_2, hay interacción. Detectarla requiere modelarla.

Para extraer el contraste de la interacción:

resultsNames(dds)
# [1] "Intercept" "condicion_tto_vs_ctrl" "tiempo_t2_vs_t1" "condicionetto.tiempo_t2"

res_int <- results(dds, name = "condicionetto.tiempo_t2")

LFC positivo significa que el efecto de tto es mayor en t2 que en t1. Las interacciones son sutiles, si tu diseño no las necesita, no las metas (sube el coste estadístico).

Recodificar referencia de un factor

A menudo el nivel de referencia que R toma alfabéticamente no es el que quieres:

dds$condicion <- relevel(dds$condicion, ref = "ctrl")

Después de relevel, hay que volver a llamar a DESeq() (porque el modelo cambia). Si lo haces antes del primer DESeq(), ya estás bien.

Contrastes explícitos: la forma definitiva

contrast = c("factor", "numerador", "denominador") es la sintaxis que recomiendo aprender bien. Funciona siempre, controla el sentido del LFC, es legible.

results(dds, contrast = c("condicion", "tto", "ctrl"))     # LFC positivo = sube tto
results(dds, contrast = c("condicion", "ctrl", "tto"))     # LFC positivo = sube ctrl

Para diseños complejos, también admite una lista de coeficientes (notación de combinaciones lineales):

results(dds, contrast = list(c("condicion_tto_vs_ctrl"),
                             c("condicion_doble_tto_vs_ctrl")))
# diferencia entre dos coeficientes ya estimados

Lo usarás raramente, para análisis estándar, la forma c("factor", "A", "B") cubre el 95 %.

Cambiar el diseño después

Puedes cambiar el design sin reconstruir el objeto:

design(dds) <- ~ batch + condicion
dds <- DESeq(dds)

Conserva los datos, vuelve a ajustar el modelo. Útil si después de PCA descubres un batch effect que merece ajustar.

Diseño sin réplicas

DESeq2 técnicamente funciona sin réplicas (1 muestra por grupo) pero no se puede estimar dispersión correctamente. Los resultados saldrán pero con confianza estadística sospechosa. No publiques análisis sin al menos 3 réplicas por grupo.

Si solo tienes 2 grupos × 2 réplicas, sí funciona, pero el poder estadístico es bajo, y conviene ser conservador con la interpretación.

Diseños balanceados vs no balanceados

Si tu batch está confundido con la condición (todos los tto en batch A, todos los ctrl en batch B), el modelo ~ batch + condicion falla o da resultados raros, DESeq2 no puede distinguir el efecto.

Cuando se detecta este caso:

table(dds$batch, dds$condicion)
#     ctrl  tto
# A      3    0
# B      0    3

Solo se puede ajustar por batch si hay algo de mezcla entre niveles. Diseñar bien al inicio del experimento (asignación aleatoria de batches) es la única solución real. Si ya está hecho mal, hay técnicas como sva o RUV para inferir factores latentes, pero salen del scope de esta ruta.

Diagnosticar el modelo

Antes de confiar en resultados, dos comprobaciones rápidas:

Confluencia del fit

plotDispEsts(dds)

Curva de dispersión razonable (cubierto en tutorial 5).

Visualizar los efectos

plotCounts(dds, gene = "ENSG00000123456", intgroup = "condicion")

Para un gen sospechoso, ver los counts crudos por grupo. Si el gen tiene padj < 0.01 pero los counts no se ven claramente distintos, revisar.

El flujo idiomático con diseño cuidado

# Definir factores con referencia explícita
dds$condicion <- factor(dds$condicion, levels = c("ctrl", "tto"))
dds$batch     <- factor(dds$batch)

# Diseño
design(dds) <- ~ batch + condicion

# Verificar que el modelo es identificable
table(dds$batch, dds$condicion)   # debe haber mezcla

# Ajustar
dds <- DESeq(dds)

# Contraste explícito
res <- results(dds, contrast = c("condicion", "tto", "ctrl"), alpha = 0.05)
summary(res)

Esto es lo que distingue un análisis profesional de uno “rápido”.

Trampas habituales

• Olvidar definir el nivel de referencia. relevel(factor, ref = "control") te ahorra encontrar que LFC positivo “significa” lo contrario.
• Confundir confounding. Si batch y condición están totalmente confundidos, no se puede ajustar. Mucha gente añade ~ batch + condicion sin mirar table(batch, condicion) y obtiene resultados sin sentido.
• Meter interacciones sin necesidad. Cada interacción consume grados de libertad y reduce el poder. Solo si hay razón biológica para esperarla.
• contrast = list(...) para algo que puede hacerse con c("factor", "A", "B"). Más simple es mejor.
• Ajustar dos veces el mismo factor. Por ejemplo, ~ batch + sex + age + condicion con sex y age correlacionados, los efectos se confunden entre ambos. Comprueba colinealidad antes con with(colData(dds), cor(...)).

En la siguiente entrega

Tienes el modelo bien definido. La siguiente entrega completa el ciclo: shrinkage de log-fold-changes con lfcShrink, qué hace, cuándo importa, y cómo cambia (o no) los volcano plots. Es lo que da resultados publicables. Lo siguiente.