Volcano y MA plots publicables

r
bioconductor
rnaseq
Cómo se construyen los dos gráficos de resultados más usados en RNA-seq. Volcano con ggplot2 desde cero, MA plot de DESeq2 vs implementación manual, y los detalles que distinguen un plot exploratorio de uno que va en un paper.

Los dos plots de resultados

Cuando un paper publica un análisis de RNA-seq, casi siempre incluye volcano plot o MA plot (a menudo los dos). Cumplen funciones distintas:

  • Volcano: log2FC en X, −log10(padj) en Y. Combina magnitud + significación. Es el plot resumen.
  • MA plot: media de expresión en X, log2FC en Y. Muestra cómo se distribuyen los cambios sobre el rango de expresión. Diagnóstico clave.

Ambos sobre 20.000+ genes a la vez. Con shrinkage aplicado (lfcShrink), ambos quedan limpios.

Volcano plot desde cero con ggplot2

DESeq2 no trae función nativa de volcano. La idiomática es construirlo con ggplot2, tienes control total.

library(ggplot2)
library(dplyr)
library(tibble)

# Convertir DESeqResults a tibble
df <- as.data.frame(res) |>
    rownames_to_column("gen_id") |>
    filter(!is.na(padj))            # quitar NAs del independent filtering

# Categoría para colorear: up / down / no significativo
df <- df |>
    mutate(estado = case_when(
        padj < 0.05 & log2FoldChange >  1 ~ "Up",
        padj < 0.05 & log2FoldChange < -1 ~ "Down",
        TRUE                              ~ "NS"
    ))

table(df$estado)
# Down    NS    Up
#  340 18500   420

Plot básico:

ggplot(df, aes(log2FoldChange, -log10(padj), color = estado)) +
    geom_point(size = 1.2, alpha = 0.6) +
    scale_color_manual(values = c(Up = "#A24B4B", Down = "#3D6FA8", NS = "grey70")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
    labs(x = "log2 fold-change", y = "-log10 p-valor ajustado",
         color = NULL) +
    theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          legend.position = "top")

Elementos clave del plot:

  • Líneas verticales en ±1: umbral de magnitud (al menos doble cambio).
  • Línea horizontal en -log10(0.05): umbral de significación.
  • Colores discretos: rojo (up), azul (down), gris (NS). Evita rainbow.
  • Alpha < 1: con miles de puntos, transparencia revela la densidad.

Etiquetar los top genes

Lo que añade carácter al plot: etiquetar los genes más destacados.

library(ggrepel)

top_genes <- df |>
    filter(padj < 0.05, abs(log2FoldChange) > 1) |>
    slice_min(padj, n = 12)

ggplot(df, aes(log2FoldChange, -log10(padj), color = estado)) +
    geom_point(size = 1.2, alpha = 0.6) +
    geom_text_repel(data = top_genes,
                    aes(label = gen_id),
                    size = 3.2,
                    box.padding = 0.5,
                    max.overlaps = 20,
                    seed = 42) +
    scale_color_manual(values = c(Up = "#A24B4B", Down = "#3D6FA8", NS = "grey70")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
    labs(x = "log2 fold-change", y = "-log10 p-valor ajustado",
         color = NULL) +
    theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          legend.position = "top")

geom_text_repel (del paquete ggrepel) evita superposición de etiquetas. seed hace el plot reproducible, sin él, las etiquetas se mueven entre renders.

Cuántos etiquetar: 8-15 es el número típico. Más es ilegible. Suele seleccionarse por mejor padj + LFC alto, o por relevancia biológica conocida (los del paper anterior, candidatos previos).

EnhancedVolcano: la alternativa pre-hecha

Si no quieres construir desde cero, el paquete EnhancedVolcano (Bioconductor) tiene defaults razonables:

library(EnhancedVolcano)

EnhancedVolcano(res,
    lab = rownames(res),
    x = "log2FoldChange",
    y = "padj",
    pCutoff = 0.05,
    FCcutoff = 1,
    pointSize = 1.5,
    labSize = 3,
    col = c("grey70", "#5F8575", "#3D6FA8", "#A24B4B"))

Pros: una línea, defaults razonables. Contras: estética opinable, menos control fino.

Para publicación, construyo casi siempre con ggplot2. Para slides exploratorios o reports rápidos, EnhancedVolcano está bien.

MA plot: la versión de DESeq2

Lo más rápido:

plotMA(res, ylim = c(-5, 5))

plotMA viene con DESeq2. Recibe directamente el DESeqResults. Por defecto:

  • Eje X: media normalizada (baseMean).
  • Eje Y: log2FC.
  • Puntos rojos: significativos (padj < alpha).
  • Eje X en escala log10 automática.

Es el diagnóstico canónico del análisis. Lo importante de leer:

  • Sin shrinkage: cola de puntos extremos a la izquierda (baja expresión).
  • Con shrinkage: cola contenida, “embudo” característico.
  • Si el plot es asimétrico (muchos más up que down, o al revés), puede haber sesgo de normalización.
  • Si los puntos rojos están solo a la derecha (alta expresión), puede haber problema de potencia en genes poco expresados.

MA plot publicable con ggplot2

Para publicación con control fino:

ggplot(df, aes(baseMean, log2FoldChange, color = estado)) +
    geom_point(size = 0.8, alpha = 0.5) +
    scale_x_log10(labels = scales::comma) +
    scale_color_manual(values = c(Up = "#A24B4B", Down = "#3D6FA8", NS = "grey70")) +
    geom_hline(yintercept = c(-1, 0, 1), linetype = c("dashed", "solid", "dashed"),
               color = c("grey50", "black", "grey50")) +
    labs(x = "Mean expression (normalized)",
         y = "log2 fold-change",
         color = NULL) +
    theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          legend.position = "top")

Detalles:

  • scale_x_log10: esencial para que los counts ocupen todo el eje.
  • Línea horizontal en 0 (referencia) + en ±1 (umbral). Tres líneas, dos discontinuas y una continua.
  • scales::comma para labels legibles (1,000 en lugar de 1e3).

Combinar varios contrastes

Cuando tu análisis tiene varios contrastes (dosis1_vs_ctrl, dosis2_vs_ctrl…), un patrón útil es un volcano por contraste en grid:

library(patchwork)

p1 <- volcano_plot(res_d1, title = "Dosis 1 vs ctrl")
p2 <- volcano_plot(res_d2, title = "Dosis 2 vs ctrl")
p3 <- volcano_plot(res_d12, title = "Dosis 2 vs Dosis 1")

(p1 | p2 | p3) + plot_layout(guides = "collect")

patchwork combina ggplots fácilmente. plot_layout(guides = "collect") une las leyendas en una sola.

Guardar para publicación

ggsave("volcano_publicable.pdf", plot = last_plot(),
       width = 6, height = 5, units = "in",
       device = cairo_pdf)

Por qué cairo_pdf: gestiona bien transparencia (alpha < 1) en PDF. El device por defecto a veces aplasta puntos transparentes a opacos.

Para PNG: device = "png", dpi = 300.

Tamaño: un volcano típico cabe bien en 6×5 inches. Más grande pierde densidad visual, más pequeño pierde legibilidad.

Trampas habituales

  • Rainbow palette. scale_color_brewer(palette = "Set1") o similar. Suficiente para 3-4 categorías. Para volcano, dos colores (up/down) + gris (NS) basta, más es ruido.
  • Etiquetar 50 genes. Ilegible. Selecciona 8-15 con slice_min(padj, n = 12).
  • No fijar seed en ggrepel. Las etiquetas se mueven en cada render. Pon seed = 42 o similar.
  • Volcano con LFC crudo. La cola fea de la izquierda no aporta. Siempre con lfcShrink.
  • MA plot sin ylim razonable. Si un gen tiene LFC = 12, el plot pierde resolución en la mayoría. ylim = c(-5, 5) recorta y se ve mejor, anota los recortados como triángulos si te importa.

En la siguiente entrega

Los plots están listos. La siguiente entrega es anotar genes: hasta ahora has trabajado con IDs (ENSG00000123456). Para reportar, necesitas convertir a símbolos (TP53, BRCA1). Lo siguiente.