Generador de volcano plots

Qué genes cambiaron entre dos condiciones, cuánto y con qué fiabilidad, todo a la vez. Cada punto es un gen: su posición de izquierda a derecha es el log₂ fold change y su altura el −log₁₀ del valor p. Cuanto más alto y más hacia fuera está un punto, mayor y más seguro es el cambio; la base gris cerca de cero son los genes que no cambiaron.

Lee en horizontal el tamaño y la dirección del cambio (la derecha es al alza, la izquierda a la baja, ±1 es una duplicación o una reducción a la mitad) y luego en vertical la evidencia (1.3 es p = 0.05, y cada unidad por encima es diez veces menor). Los resultados quedan en las dos esquinas superiores, fuera de las líneas discontinuas. El recorrido de arriba lee las seis regiones en un ejemplo de 20 genes.

|log₂FC| ≥ 1, un cambio de 2 veces, es el más común; 0.58 (1.5 veces) se usa donde se esperan efectos pequeños, como en proteómica, microarrays o datos de célula única. El umbral es una decisión de presentación, así que indícalo en la leyenda de la figura. En el generador de volcano plots, escríbelo en la instrucción («|log2FC| > 0.58») o como un cambio de N veces («1.5 veces»), que se convierte a log₂.

El que se usó para el umbral de significación, para que la línea horizontal signifique lo que dice la leyenda. Con miles de genes analizados, la mayoría de los artículos aplican el umbral y representan el p ajustado por Benjamini–Hochberg (padj, FDR, q). El generador de volcano plots titula el eje −log₁₀(valor p ajustado) cuando los valores están ajustados y −log₁₀(valor p) en caso contrario.

Porque subieron más genes de los que bajaron, o al revés: un volcano plot es tan simétrico como la biología. Que toda la base gris se desplace lejos de x = 0 es otra cosa. Suele indicar un problema de normalización y conviene revisarlo en un MA plot antes de leer los resultados. El generador de volcano plots nunca refleja la nube para que parezca equilibrada.

Sus valores p eran demasiado pequeños para almacenarse y se escribieron como 0, así que −log₁₀ p es infinito y el programa de gráficos los fija a una misma altura. Es frecuente con el p_val_adj de Seurat y con RNA-seq muy profundo. El generador de volcano plots dibuja esos genes en el menor p distinto de cero de la tabla y lo indica en la figura; darlos como p < 10⁻³⁰⁰ es honesto, inventar un valor exacto no.

Con ggplot2, añade una columna de clase (al alza, a la baja, no significativo) y luego ggplot(res, aes(log2FoldChange, -log10(padj), colour = class)) + geom_point() + geom_vline(xintercept = c(-1, 1), linetype = 2) + geom_hline(yintercept = -log10(0.05), linetype = 2), con ggrepel para las etiquetas. EnhancedVolcano reúne los mismos pasos en sus argumentos pCutoff y FCcutoff. Para saltarte el código, adjunta al generador de volcano plots la tabla que guardó write.csv().

Prism dibuja un volcano plot automáticamente cuando analizas una tabla agrupada con múltiples pruebas t (Prism 8 y posteriores): los resultados incluyen un gráfico de la diferencia frente al −log₁₀ p. Con una tabla de resultados hecha en otro sitio, introduce el log₂ fold change como X y el −log₁₀ p como Y en una tabla XY y añade las líneas de umbral como conjuntos de datos adicionales, o adjunta esa tabla al generador de volcano plots.

En Excel, añade una columna =-LOG10(C2) junto a los valores p, haz una dispersión XY del fold change frente a ella y añade cada línea de umbral como su propia serie de dos puntos; colorear las clases exige una serie por clase. En Python, plt.scatter(df.log2FoldChange, -np.log10(df.padj), c=colours) con plt.axvline y plt.axhline hace lo mismo, y adjustText coloca las etiquetas. El generador de volcano plots acepta la misma hoja en .xlsx o .csv sin columna de fórmula.

Sí: el generador de volcano plots lee los resultados de DESeq2, edgeR, limma, Seurat y Perseus por los encabezados de sus columnas, en .csv, .tsv o .xlsx, hasta 60 000 filas. Las filas con NA se omiten, se usa el p ajustado cuando el archivo lo tiene y una columna Fold change simple se convierte a log₂.

No: el generador de volcano plots lee los valores p de tu tabla y nunca ejecuta la estadística. Haz la prueba en DESeq2, edgeR, limma o tu programa de proteómica; a partir de esos números el generador de volcano plots calcula exactamente el −log₁₀ p, la clase de cada gen, los recuentos, los ejes y las líneas de umbral.

Por defecto el generador de volcano plots etiqueta los diez resultados más significativos, repartidos entre el lado sobreexpresado y el infraexpresado. Escribe «etiqueta los 20 primeros», «sin etiquetas» o nombra los genes («resalta NANOG, SOX2 y MYC»); los genes nombrados se etiquetan aunque estén en gris.

El generador de volcano plots admite hasta 60 000 filas, suficiente para un transcriptoma completo. Hasta 200 genes, cada punto se coloca en sus coordenadas exactas. Por encima, los genes más significativos (hasta 300) se colocan exactamente y el resto se dibuja a partir de sus recuentos y su dispersión reales, franja a franja, mientras la leyenda conserva los recuentos exactos de sobreexpresados, infraexpresados y no significativos.