Normalización y transformaciones (vst, rlog)

r
bioconductor
rnaseq
Por qué los counts crudos y el log ingenuo no sirven para visualizar, qué hacen vst y rlog, cuándo usar cada uno y cómo encajan en el flujo de PCA y heatmap.

El problema con counts crudos

Los counts de RNA-seq tienen tres propiedades incómodas para visualización:

  1. Profundidades distintas entre muestras (cubierto con size factors).
  2. Distribución muy sesgada: la mayoría de genes tienen pocos counts, unos pocos tienen miles.
  3. Heterocedasticidad: la varianza crece con la media. Genes muy expresados muestran fluctuaciones mucho mayores que genes poco expresados.

Para PCA, heatmap, clustering, herramientas que asumen escalas comparables y varianzas controladas, los counts crudos llevan a que los pocos genes muy expresados dominen todo. El resultado: un PCA donde solo se ve la variación de 10 genes housekeeping y nada del biológicamente interesante.

¿Por qué no log(counts + 1)?

La intuición inmediata: log2(counts + 1) aplasta la escala. Funciona… a medias.

El problema es el pseudocount. Con +1, un gen con 0 counts queda en log2(1) = 0. Uno con 1 count, en log2(2) = 1. Esa diferencia de 1 unidad en log para genes apenas detectados es ruido puro, no señal. Y se ve gigante al lado de una diferencia de 1 unidad en log entre 1000 y 2000 counts (que sí significa algo).

Resultado: el log ingenuo invierte el problema, ahora los genes con counts bajos dominan con su ruido en lugar de los muy expresados con su escala.

vst() y rlog(): la solución de DESeq2

DESeq2 trae dos transformaciones diseñadas para estabilizar la varianza sin sobre-amplificar el ruido de counts bajos:

  • vst(): Variance Stabilizing Transformation. Aproximación analítica.
  • rlog(): regularized log. Encogimiento bayesiano hacia el promedio.

Ambas producen datos en escala parecida a log2, pero con varianza aproximadamente constante entre genes con distintas medias.

library(DESeq2)

vsd <- vst(dds, blind = FALSE)
rld <- rlog(dds, blind = FALSE)

El objeto resultante es un DESeqTransform, un SummarizedExperiment con la matriz transformada en assay().

vst vs rlog: cuándo cada uno

vst() rlog()
Velocidad Rápida, segundos para 30k×100 Lenta, minutos para 30k×100
Datasets pequeños (n < 30) Buena Mejor
Datasets grandes (n > 30) Buena Innecesariamente cara
Resultado Parecido a rlog para grandes Más cuidadoso con counts bajos

Regla práctica (la oficial de la vignette):

  • vst() por defecto. Especialmente con n > 30 o cuando el tiempo importa.
  • rlog() con datasets pequeños y heterogéneos donde rlog suaviza mejor.

Para 99 % de los análisis, vst() está bien.

El argumento blind

vsd <- vst(dds, blind = FALSE)
  • blind = TRUE: transformación ignora el diseño. Útil para QC inicial donde quieres ver agrupamientos naturales sin que el modelo “informe” la transformación.
  • blind = FALSE: usa el diseño para calcular dispersiones. Recomendado para análisis downstream (PCA con expectativa biológica, heatmap de top DE).

La vignette de DESeq2 explica: para QC antes del análisis (¿se mezclan tratamientos?), blind = TRUE. Para visualización después del análisis, blind = FALSE. La diferencia visual suele ser pequeña.

Uso en PCA

Esta es la combinación canónica:

plotPCA(vsd, intgroup = c("condicion", "batch"))

plotPCA viene con DESeq2. Por defecto usa los 500 genes con más varianza tras la transformación. Eso es lo correcto, sin la transformación, “más varianza” significaría “más expresado”.

Argumentos útiles:

plotPCA(vsd, intgroup = "condicion", ntop = 1000)
plotPCA(vsd, intgroup = "condicion", returnData = TRUE)  # devuelve los puntos para ggplot

Si quieres control total con ggplot2:

library(ggplot2)
pca_data <- plotPCA(vsd, intgroup = c("condicion", "batch"), returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = condicion, shape = batch)) +
    geom_point(size = 4) +
    xlab(paste0("PC1: ", percent_var[1], "% variance")) +
    ylab(paste0("PC2: ", percent_var[2], "% variance")) +
    coord_fixed()

Lo veremos a fondo en el tutorial 9 (QC visual).

Uso en heatmap de muestras

Otro uso clásico: matriz de distancias entre muestras para detectar agrupamientos esperados:

library(pheatmap)
sample_dist <- dist(t(assay(vsd)))
sample_mat <- as.matrix(sample_dist)
rownames(sample_mat) <- paste(vsd$condicion, vsd$batch, sep = "_")
colnames(sample_mat) <- NULL

pheatmap(sample_mat,
         clustering_distance_rows = sample_dist,
         clustering_distance_cols = sample_dist)

dist(t(assay(vsd))), calcula distancia euclídea entre muestras tras transformar y transponer (porque dist opera sobre filas y queremos distancias entre muestras, que son columnas).

Si las muestras del mismo grupo no se agrupan, hay un problema: o batch effect fuerte, o muestras mal etiquetadas, o el efecto biológico es muy débil.

Uso en heatmap de top genes

Mostrar los N genes más diferencialmente expresados:

top_genes <- rownames(res)[order(res$padj)][1:30]

mat <- assay(vsd)[top_genes, ]
mat <- mat - rowMeans(mat)        # centrar por gen

annot <- as.data.frame(colData(vsd)[, c("condicion", "batch")])

pheatmap(mat, annotation_col = annot,
         show_rownames = TRUE,
         scale = "none",            # ya centramos a mano
         color = colorRampPalette(c("#3F6354", "white", "#A24B4B"))(50))

Centrar por gen (rowMeans) es necesario: sin centrar, el heatmap se domina por la magnitud de cada gen, no por la diferencia entre muestras.

El flujo idiomático

Lo que se ejecuta típicamente justo después de DESeq():

dds <- DESeq(dds)
vsd <- vst(dds, blind = FALSE)

# QC visual
plotPCA(vsd, intgroup = "condicion")

# Análisis diferencial
res <- results(dds)

DESeq() para inferencia (counts crudos por dentro), vst() para visualización. Dos objetos paralelos sobre el mismo dataset.

Trampas habituales

  • Usar la matriz transformada para DESeq(). DESeq necesita counts crudos enteros. El vsd es solo para visualización y exploración.
  • rlog con datasets grandes. Tarda mucho y no mejora vst materialmente. Si N > 30 y tienes prisa, vst.
  • blind = TRUE cuando ya hay diseño que considerar. Para visualizar resultados del modelo, blind = FALSE da una vista más coherente.
  • Heatmap sin centrar por gen. Domina la magnitud absoluta de expresión, no las diferencias. Centrar (- rowMeans()) o usar scale = "row" en pheatmap.
  • PCA con counts crudos o log(counts+1). Producen PCAs que parecen razonables pero están dominados por housekeeping de alta expresión. Es la causa #1 de PCAs “que no separan nada interesante”.

En la siguiente entrega

Has visto el flujo con un factor simple (~ condicion). El siguiente paso son diseños reales: ajustar por batch, modelar interacciones, definir contrastes con contrast = c(...). Es donde más se nota el rigor estadístico y donde más errores aparecen. Lo siguiente.