QC visual: PCA y heatmap de muestras

r
bioconductor
rnaseq
Cómo se diagnostica un análisis RNA-seq antes de mirar resultados. PCA de muestras leído con criterio, heatmap de distancias entre muestras y qué significa que las réplicas no se agrupen donde tocan.

Antes de mirar resultados, mira muestras

Primera regla del análisis serio: el primer plot que ves no es el volcano. Es el PCA y el heatmap de muestras. Si esos dos están mal, los resultados no se interpretan, se ignoran.

Lo que diagnostica el QC visual:

  • ¿Las réplicas se agrupan?
  • ¿Los grupos experimentales se separan en PC1 o PC2?
  • ¿Hay alguna muestra claramente outlier?
  • ¿Hay un batch effect dominando la variación?
  • ¿Las etiquetas están bien?

Cinco preguntas antes de pasar al análisis estadístico.

PCA: el plot fundamental

library(DESeq2)
vsd <- vst(dds, blind = FALSE)

plotPCA(vsd, intgroup = "condicion")

plotPCA viene con DESeq2. Por defecto:

  • Usa los 500 genes con más varianza tras vst.
  • Centra y escala.
  • Devuelve un ggplot listo.

Lectura inicial:

  • Eje X (PC1) captura la mayor fuente de variación. Lo ideal: que separe los grupos biológicos.
  • Eje Y (PC2) captura la segunda fuente. A veces es la condición (si el efecto biológico está en PC2), a veces es batch, a veces es covariable.

El % de varianza explicada sale en los ejes. Si PC1 captura el 80 %, hay una fuente dominante. Si PC1 + PC2 < 25 %, la variación está repartida, análisis posiblemente difícil.

Versión publicable

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, alpha = 0.85) +
    xlab(paste0("PC1: ", percent_var[1], "% varianza")) +
    ylab(paste0("PC2: ", percent_var[2], "% varianza")) +
    scale_color_manual(values = c("#5F8575", "#A24B4B")) +
    coord_fixed() +
    theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank())

coord_fixed() mantiene la relación 1:1 entre ejes, importante porque las distancias en PC1 y PC2 son comparables. Sin coord_fixed, el plot puede engañar.

shape = batch y color = condicion permiten ver dos factores a la vez. Si los grupos se separan por color (condicion) y los lotes se mezclan dentro de cada color (formas mezcladas), el diseño está limpio. Si los lotes se agrupan dentro del color, hay batch effect.

Casos típicos en el PCA

Caso ideal

PC1 separa por condicion (color), PC2 sin patrón claro

Análisis straightforward. Sigue adelante.

Batch effect ligero

PC1 separa por condicion, PC2 separa por batch

Razonable. Modelar batch como covariable (~ batch + condicion) lo controla. No es problema.

Batch effect dominante

PC1 separa por batch, condicion solo se ve en PC3 o no se ve

Problema serio. Hay más varianza técnica que biológica. Si batch y condicion están balanceados (no confundidos), ~ batch + condicion rescata el análisis. Si están confundidos, no.

Muestra outlier

Una muestra sola en una esquina, el resto agrupado

Investigar. Comprueba: - Library size (colSums(counts)). - Tasa de mapeo (del MultiQC). - Si es real (problema biológico inesperado, error de etiquetado).

A veces hay que excluir la muestra y rehacer. Documenta la decisión.

Etiquetas intercambiadas

Algunas muestras de "tto" aparecen en el grupo de "ctrl" en PC1

Bandera roja inmediata. Antes de “explicar biológicamente”, verifica la hoja de muestras. Es la causa #1 de bugs en bioinformática.

Heatmap de distancias entre muestras

Una segunda visualización del mismo principio, más informativa para experimentos pequeños:

library(pheatmap)
library(RColorBrewer)

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

colores <- colorRampPalette(rev(brewer.pal(9, "Greens")))(255)

pheatmap(sample_mat,
         clustering_distance_rows = sample_dist,
         clustering_distance_cols = sample_dist,
         col = colores,
         border_color = NA,
         fontsize = 10)

Pasos clave:

  • dist(t(assay(vsd))): distancia euclídea entre muestras (transpongo porque dist opera sobre filas).
  • clustering_distance_*: usa la misma matriz para el dendrograma.
  • paste(condicion, batch): labels que muestran la información clave.

Cómo leerlo:

  • Bloques diagonales = grupos que se parecen entre sí. Lo deseable: bloques que coincidan con la condición.
  • Una muestra que se aleja de todas las demás = posible outlier.
  • Dos muestras con distancia casi cero = sospecha de duplicación técnica o etiqueta mal.

Combinar PCA y heatmap

PCA y heatmap muestran la misma información desde dos ángulos. Reportar ambos suele ser excesivo, pero mirar los dos durante el análisis es buena práctica:

  • Para datasets pequeños (n < 20), heatmap de distancias es más informativo, cada muestra etiquetada.
  • Para datasets grandes (n > 50), PCA es esencial, heatmap se vuelve ilegible.

En un informe profesional, PCA en figura principal, heatmap en suplemento (o al revés).

Sample-sample correlation: variante útil

En vez de distancia euclídea, correlación:

cor_mat <- cor(assay(vsd))
pheatmap(cor_mat, ...)

Valores cerca de 1 = muestras muy parecidas. Más interpretable que distancia para gente no técnica. La distribución típica:

  • Correlación intra-grupo > 0.95, normal.
  • Correlación entre grupos 0.85-0.95, normal con efecto biológico moderado.
  • Correlación < 0.8 entre réplicas técnicas, sospechoso.

Qué hacer si el QC va mal

Cinco salidas posibles ordenadas de menor a mayor coste:

  1. Modelar batch (~ batch + condicion), gratis si tienes la info.
  2. Excluir outliers documentados: pierdes potencia, ganas claridad. Justifica en el informe.
  3. Inferir factores latentes (sva, RUVSeq), útil cuando hay batch sin observar. Requiere tutorial aparte.
  4. Repetir secuenciación de muestras problemáticas, si presupuesto y tiempo lo permiten.
  5. Rehacer el experimento con mejor diseño: la solución correcta cuando lo demás falla.

La opción 5 nadie quiere oír, pero a veces es la honesta. Mejor reconocer un experimento mal diseñado que publicar resultados que no se replican.

Trampas habituales

  • Ignorar el porcentaje de varianza explicada. Un PCA con PC1 al 80 % es muy distinto de uno con PC1 al 12 %. Reporta siempre los %.
  • coord_fixed() olvidado. Sin él, dos ejes con escalas distintas engañan visualmente, parece que los grupos se separan en X cuando solo se ven separados por la escala.
  • PCA con counts crudos. Dominado por housekeeping de alta expresión. Siempre vsd o rld antes.
  • PCA con todos los genes. Mejor con los 500-2000 más variables. plotPCA lo hace por defecto.
  • Heatmap con escala distinta a dist. Es habitual usar correlación para el cluster y distancia para el color, dando dendrogramas inconsistentes con el heatmap. Mantén consistencia.

En la siguiente entrega

Tienes el QC. Lo siguiente son los gráficos del resultado del análisis: volcano plots y MA plots publicables, con los detalles que separan un plot exploratorio de uno presentable. Lo siguiente.