Volcano plot de RNA-Seq: puntos de corte, código y trampas

Un volcano plot pone cada gen de una prueba de expresión diferencial en un solo panel. El eje x es el cambio de pliegue log2. El eje y es −log10 del p-value ajustado. Los genes lejos del centro y en la parte alta son tus candidatos.

El gráfico se hace con diez líneas de código. Leerlo correctamente requiere más cuidado. Las dos líneas discontinuas que la mayoría copia, |log2FC| ≥ 1 y padj < 0.05, son un hábito de visualización, no una regla — y la guía de usuario de edgeR dice que usar ese par para seleccionar genes "destruirá el control del FDR en general" (sección 2.14).

Esta entrada cubre qué significa cada eje, qué punto de corte dibujar y por qué, cómo leer cada región de la nube y código R y Python funcional. Si solo necesitas la figura, nuestro Generador de volcano plot toma una tabla de resultados y devuelve un gráfico etiquetado.

Volcano plot a partir de datos simulados (12 000 genes, 5 % con un efecto real). Hecho con matplotlib; el código que generó los datos está en la sección de Python más abajo, y el script de trazado completo se encuentra junto a esta imagen en el repositorio.

Los dos ejes, y por qué ambos están en logaritmo

Ambos ejes están en logaritmo por la misma razón: las escalas crudas ocultan la parte que te importa.

El cambio de pliegue en una escala cruda es asimétrico. Un gen que se duplica obtiene 2.0 y un gen que se reduce a la mitad obtiene 0.5, así que lo que sube tiene espacio hasta el infinito mientras que lo que baja queda comprimido bajo 1. El log2 lo arregla: duplicarse se convierte en +1, reducirse a la mitad en −1, y los dos lados del gráfico son imágenes especulares.

Los p-values son el problema opuesto. Se acumulan cerca de cero, donde viven todos los interesantes. Aplicar −log10 extiende esa acumulación, y el signo menos coloca la evidencia más fuerte en la parte superior.

De esa tabla se desprenden dos hechos. Un punto en y = 4 tiene un padj cien veces menor que un punto en y = 2. Y el eje y no tiene límite superior, así que un solo gen con un p-value diminuto puede estirar todo el gráfico.

Pon padj en el eje y, no el p-value crudo

Los p-values crudos en el eje y te darán cientos de falsos positivos. El RNA-seq en masa analiza todos los genes, así que un experimento con 15 000 genes analizados y sin efecto real aún da unos 750 genes con p < 0.05.

La solución es una corrección de la tasa de descubrimientos falsos. El método de Benjamini-Hochberg reescala los p-values para que un punto de corte de 0.05 signifique "alrededor del 5 % de los genes que conservo son falsos", no "el 5 % de todas las pruebas". R lo hace en p.adjust con method = "BH", según Benjamini y Hochberg (1995).

Ambas herramientas principales ya te dan la columna corregida:

• DESeq2 (viñeta): log2FoldChange, pvalue, padj. La columna padj es BH por defecto, ya que results() se ejecuta con pAdjustMethod = "BH".
• edgeR (guía de usuario): topTags() devuelve logFC, logCPM, una columna de estadístico de prueba (F o LR), PValue y FDR. FDR es la columna que quieres.

Dos detalles de DESeq2 causan problemas aquí. Primero, results() usa alpha = 0.1 por defecto, y su paso de filtrado independiente está ajustado a ese valor. Si tu línea discontinua está en 0.05, llama a results(dds, alpha = 0.05) para que el filtrado coincida con la línea que dibujaste.

Segundo, algunos genes devuelven padj = NA. Eso es el filtrado independiente y la eliminación de valores atípicos haciendo su trabajo, según la sección "¿Por qué algunos p-values se establecen en NA?" de la viñeta. Elimina esas filas antes de graficar. Nunca conviertas NA en 1 o 0.

Dónde dibujar las líneas de corte

No existe un par de puntos de corte estándar, solo los habituales. Aquí está de dónde vienen realmente los números populares:

Elige según la cantidad de señal que tengas, no por costumbre. Si 6 000 genes pasan, la lista es inútil y una barra de cambio de pliegue más alta ayuda. Si nueve genes pasan, elimina la línea de cambio de pliegue y ordena por padj en su lugar.

La guía de edgeR tiene la mejor regla de una línea para el punto de corte en x. Lee la como "el cambio de pliegue por debajo del cual definitivamente no nos interesa el gen", no el cambio por encima del cual un gen se vuelve interesante.

La trampa del punto de corte que conviene conocer

Filtrar por p y luego por cambio de pliegue no es una prueba válida. La guía de edgeR es contundente al respecto: tales combinaciones "son ad hoc y no dan p-values con sentido para probar expresiones diferenciales relativas a un umbral de cambio de pliegue". "favorecen a genes poco expresados pero muy variables" y "destruyen el control del FDR en general".

Si la barra de cambio de pliegue importa para tu afirmación, prueba contra ella en lugar de filtrar a posteriori:

# DESeq2: test whether |log2FC| is above 1, not just above 0
res <- results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs", alpha = 0.05)

# edgeR: same idea, via the TREAT method
fit <- glmQLFit(y, design)
tr  <- glmTreat(fit, coef = 2, lfc = 1)

Ahora padj ya tiene en cuenta el punto de corte de tamaño, y las líneas verticales describen la prueba que ejecutaste. Ambas opciones están documentadas: consulta pruebas de cambio de pliegue log2 por encima de un umbral en DESeq2 y el artículo de TREAT detrás de glmTreat.

Cómo leer el gráfico, región por región

Lee la nube antes de leer las etiquetas. La forma dice más sobre el experimento que cualquier gen individual.

La parte superior central es la región que la gente lee mal. Esos genes suelen ser genes con recuentos altos, donde el error estándar es pequeño, así que incluso un cambio del 15 % supera el punto de corte de p. El cambio es real, pero si un cambio del 15 % importa es una cuestión de biología, no de estadística. Una señal de advertencia: si los genes de mantenimiento o ribosómicos están ahí arriba, sospecha de un problema de normalización o de composición de la librería, no de biología.

Los bordes inferiores son la imagen especular. Un gen con 4 recuentos en un grupo y 30 en el otro da un cambio de pliegue log2 enorme y casi ninguna evidencia. Estos puntos ensanchan el eje x y hacen la historia emocionante, y rara vez se replican. Los cambios de pliegue reducidos arreglan la visualización: lfcShrink(dds, coef = 2, type = "apeglm") acerca los genes mal medidos a cero y deja intactos los bien medidos.

Dos comprobaciones de forma cierran la lectura:

• Simetría. Una nube desequilibrada, por ejemplo 900 sobreexpresados y 40 reprimidos, puede ser una activación real. También puede venir de un cambio de composición que la normalización no absorbió. Mira el gráfico MA y los factores de tamaño antes de escribir la frase.
• Anchura de la V. El hueco alrededor de x = 0 con puntos que ascienden por ambos flancos es esperable, ya que los cambios más grandes superan el punto de corte de p más fácilmente. En cambio, un pico vertical sólido en x = 0 suele significar que trazaste los p-values crudos.

Gráficalo en R a partir de resultados de DESeq2

Esto se ejecuta sobre un DESeqDataSet llamado dds y necesita ggplot2, ggrepel y dplyr.

library(DESeq2)
library(ggplot2)
library(ggrepel)
library(dplyr)

lfc_cut  <- 1
padj_cut <- 0.05

# keep alpha equal to the line you draw, or filtering optimises for 0.1
res <- results(dds, alpha = padj_cut)

df <- as.data.frame(res)
df$gene <- rownames(df)
df <- df[!is.na(df$padj), ]
df$padj <- pmax(df$padj, .Machine$double.xmin)  # padj can underflow to 0

df$class <- "Not significant"
df$class[df$padj < padj_cut & df$log2FoldChange >=  lfc_cut] <- "Up"
df$class[df$padj < padj_cut & df$log2FoldChange <= -lfc_cut] <- "Down"

top <- df %>%
  filter(class != "Not significant") %>%
  arrange(padj) %>%
  slice_head(n = 10)

ggplot(df, aes(x = log2FoldChange, y = -log10(padj), colour = class)) +
  geom_point(size = 1, alpha = 0.7) +
  geom_hline(yintercept = -log10(padj_cut), linetype = "dashed") +
  geom_vline(xintercept = c(-lfc_cut, lfc_cut), linetype = "dashed") +
  geom_text_repel(
    data = top, aes(label = gene),
    colour = "black", size = 3, max.overlaps = Inf, show.legend = FALSE
  ) +
  scale_colour_manual(values = c(
    "Up" = "#c23b3b", "Down" = "#2f6fb2", "Not significant" = "#b8bec9"
  )) +
  coord_cartesian(xlim = c(-5, 5)) +
  labs(x = "log2 fold change", y = "-log10 padj", colour = NULL) +
  theme_classic()

Si vienes de edgeR, renombra dos columnas y el resto del script funciona:

df <- topTags(qlf, n = Inf)$table
df$log2FoldChange <- df$logFC
df$padj <- df$FDR

Para una versión completamente estilizada con conectores y cajas de corte sombreadas, la viñeta de EnhancedVolcano merece la pena leerla.

Gráficalo en Python desde un DataFrame de pandas

Exporta la tabla de DESeq2 a CSV y luego ejecuta esto. Solo usa pandas, numpy y matplotlib.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# first column is the gene name; needs DESeq2's log2FoldChange and padj
df = pd.read_csv("deseq2_results.csv", index_col=0)
df = df.dropna(subset=["padj", "log2FoldChange"]).copy()

lfc_cut, padj_cut = 1.0, 0.05

# padj can underflow to 0, which makes -log10 infinite; clip it first
df["y"] = -np.log10(df["padj"].clip(lower=np.finfo(float).tiny))
df["cls"] = "not significant"
df.loc[(df["padj"] < padj_cut) & (df["log2FoldChange"] >= lfc_cut), "cls"] = "up"
df.loc[(df["padj"] < padj_cut) & (df["log2FoldChange"] <= -lfc_cut), "cls"] = "down"

colours = {"not significant": "#b8bec9", "down": "#2f6fb2", "up": "#c23b3b"}
fig, ax = plt.subplots(figsize=(7, 5), dpi=200)
for name, grp in df.groupby("cls"):
    ax.scatter(grp["log2FoldChange"], grp["y"], s=6, c=colours[name],
               label=name, alpha=0.7, linewidths=0)

ax.axhline(-np.log10(padj_cut), ls="--", lw=0.9, c="#444444")
for v in (-lfc_cut, lfc_cut):
    ax.axvline(v, ls="--", lw=0.9, c="#444444")

top = df[df["cls"] != "not significant"].nsmallest(10, "padj")
for gene, row in top.iterrows():
    ax.annotate(gene, (row["log2FoldChange"], row["y"]), xytext=(4, 4),
                textcoords="offset points", fontsize=7)

ax.set_xlabel("log2 fold change")
ax.set_ylabel("-log10 padj")
ax.set_xlim(-5, 5)
ax.legend(frameon=False, markerscale=2)
fig.tight_layout()
fig.savefig("volcano.png")

Si trabajas con la salida de edgeR en Python, renombra primero logFC a log2FoldChange y FDR a padj.

El código detrás de la figura anterior

La figura de ejemplo usa datos simulados, no un conjunto de datos real. El modelo es deliberadamente simple: las medias base se extienden a lo largo de cuatro órdenes de magnitud, el 5 % de los genes tienen un efecto real, y el error estándar crece a medida que la media base disminuye. Esa última parte es lo que crea los puntos de baja evidencia en los bordes lejanos.

import numpy as np
from math import erfc

rng = np.random.default_rng(11)
n = 12000

base_mean = 10 ** rng.uniform(0, 4.2, n)
true_lfc = np.where(rng.random(n) < 0.05, rng.normal(0, 1.1, n), 0.0)
se = np.sqrt(0.035 + 9.0 / base_mean)        # low-count genes get a larger standard error
lfc = true_lfc + rng.normal(0, se)
pval = np.array([erfc(abs(s) / np.sqrt(2)) for s in lfc / se])

# Benjamini-Hochberg, same as R's p.adjust(method = "BH")
order = np.argsort(pval)
ranked = pval[order] * n / np.arange(1, n + 1)
padj = np.empty(n)
padj[order] = np.clip(np.minimum.accumulate(ranked[::-1])[::-1], 0, 1)
padj[base_mean < 1.5] = np.nan                # stand-in for independent filtering

Cinco errores que hacen que un volcano plot sea engañoso

• P-values crudos en el eje y. Cientos de genes cruzan la línea por azar. Usa padj o FDR, y di cuál en la etiqueta del eje.
• Tratar la línea de cambio de pliegue como significado biológico. Un cambio de 2 veces en un factor de transcripción y en una proteína estructural no son afirmaciones comparables. La línea es un filtro, no evidencia.
• Dejar que los valores atípicos estiren los ejes. Un gen con log2FC = 12 aplana todo lo demás. Recorta la vista con coord_cartesian() o set_xlim(), y marca los puntos recortados con triángulos en el borde en lugar de eliminarlos.
• Colorear solo por el cambio de pliegue. Rojo para cada gen más allá de x = 1 marca como resultados a genes ruidosos de bajo recuento. El color necesita ambas condiciones: más allá de la línea de cambio de pliegue y más allá de la línea de padj.
• Ninguna línea discontinua. Sin ellas, el lector no puede distinguir tu punto de corte en el gráfico. Dibuja ambas, y pon los números en el pie de figura o en la leyenda.

Preguntas frecuentes

¿Es un volcano plot lo mismo que un gráfico MA?

No, y responden a preguntas distintas. Un gráfico MA pone la expresión media en el eje x y el cambio de pliegue log2 en el eje y, así que muestra si tus cambios grandes provienen de genes de bajo recuento. Un volcano plot deja de lado el nivel de expresión y muestra evidencia en su lugar. Haz ambos mientras trabajas, publica el volcano.

¿Qué hago cuando padj es exactamente 0?

Eso es subdesbordamiento de coma flotante, no un p-value de cero. -log10(0) da infinito y el punto desaparece o rompe el eje. Recorta padj al double positivo más pequeño antes del logaritmo, como en el código de Python anterior, y reporta el gen como padj < 1e-300 en el texto.

¿Cuántos genes debería etiquetar?

De diez a veinte es legible, más es un desastre. Ordena por padj entre los genes que también pasan la línea de cambio de pliegue, y luego etiqueta un número fijo. Usa ggrepel en R o llamadas annotate desplazadas en Python para que el texto no se asiente sobre los puntos. Etiqueta cualquier gen que menciones en el resumen, aunque no esté entre los diez primeros.

¿Puede un volcano plot ir directamente a un artículo?

Sí, si se exporta como arte vectorial o a la resolución que pide la revista, y si el análisis subyacente es reproducible. Las editoriales tratan las figuras de datos de forma distinta a las ilustraciones, y las reglas importan antes de que envíes — consulta nuestra guía sobre políticas de las revistas sobre figuras generadas por IA. Para el panel resumen de un artículo, un Generador de Resumen Gráfico cubre la parte esquemática, mientras que el volcano sigue siendo un gráfico de tus propios números.