RNA-Seq-volcano plot: Grenzwerte, Code und Fallstricke

Ein volcano plot bringt jedes Gen aus einem Test auf differentielle Expression in einem einzigen Panel zusammen. Die x-Achse ist der log2-Fold-Change. Die y-Achse ist −log10 des adjustierten p-Werts. Gene, die weit vom Zentrum entfernt und weit oben liegen, sind Ihre Kandidaten.

Der volcano plot entsteht mit zehn Codezeilen. Ihn korrekt zu lesen, erfordert mehr Sorgfalt. Die beiden gestrichelten Linien, die die meisten kopieren, |log2FC| ≥ 1 und padj < 0.05, sind eine Darstellungsgewohnheit, keine Regel — und das edgeR-Benutzerhandbuch sagt, dass die Verwendung dieses Paars zur Auswahl von Genen „destroy the control of FDR in general“ werde (Abschnitt 2.14).

Dieser Beitrag behandelt, was jede Achse bedeutet, welchen Grenzwert Sie einzeichnen sollten und warum, wie Sie jede Region der Punktwolke lesen und wie funktionierender R- und Python-Code aussieht. Wenn Sie nur das Diagramm benötigen, nimmt unser volcano plot generator eine Ergebnistabelle entgegen und gibt ein beschriftetes Diagramm zurück.

Volcano plot aus simulierten Daten (12,000 Gene, 5% mit einem echten Effekt). Erstellt mit matplotlib; der Code zur Datengenerierung steht im Python-Abschnitt unten, und das vollständige Skript zum Erstellen der Grafik liegt im Repository neben diesem Bild.

Die beiden Achsen und warum beide logarithmiert sind

Beide Achsen sind aus demselben Grund logarithmiert: Die Rohskalen verbergen den für Sie wichtigen Bereich.

Fold Change auf einer Rohskala ist asymmetrisch. Ein Gen, das sich verdoppelt, erhält 2.0, und ein Gen, das sich halbiert, erhält 0.5; nach oben reicht die Skala bis unendlich, während sie nach unten unter 1 zusammengedrängt wird. log2 behebt das: Verdoppeln wird zu +1, Halbieren zu −1, und die beiden Seiten des Diagramms sind Spiegelbilder.

P-Werte stellen das umgekehrte Problem dar. Sie häufen sich nahe null, wo alle interessanten Werte liegen. −log10 verteilt diese Ansammlung, und das Minuszeichen setzt die stärkste Evidenz nach oben.

Aus dieser Tabelle folgen zwei Fakten. Ein Punkt bei y = 4 hat einen hundertmal kleineren padj-Wert als ein Punkt bei y = 2. Die y-Achse hat außerdem keine Obergrenze, sodass ein Gen mit winzigem p-Wert das gesamte Diagramm in die Länge ziehen kann.

padj auf die y-Achse setzen, nicht den Roh-p-Wert

Roh-p-Werte auf der y-Achse liefern Ihnen Hunderte falscher Treffer. Bulk-RNA-seq testet jedes Gen, daher liefert ein Lauf mit 15,000 getesteten Genen und keinem echten Effekt trotzdem etwa 750 Gene bei p < 0.05.

Die Lösung ist eine Korrektur der False Discovery Rate. Die Benjamini-Hochberg-Methode skaliert p-Werte neu, sodass ein Grenzwert von 0.05 „etwa 5 % der Gene, die ich behalte, sind falsch“ bedeutet, nicht „5 % aller Tests“. R setzt dies in p.adjust mit method = "BH" um, nach Benjamini und Hochberg (1995).

Beide wichtigsten Tools liefern die korrigierte Spalte bereits:

• DESeq2 (Vignette): log2FoldChange, pvalue, padj. Die padj-Spalte verwendet standardmäßig BH, da results() mit pAdjustMethod = "BH" ausgeführt wird.
• edgeR (Benutzerhandbuch): topTags() gibt logFC, logCPM, eine Teststatistik-Spalte (F oder LR), PValue und FDR zurück. FDR ist die gewünschte Spalte.

Zwei DESeq2-Details führen hier häufig zu Fehlern. Erstens verwendet results() standardmäßig alpha = 0.1, und sein unabhängiger Filterungsschritt ist auf diesen Wert abgestimmt. Wenn Ihre gestrichelte Linie bei 0.05 liegt, rufen Sie results(dds, alpha = 0.05) auf, damit die Filterung zur gezeichneten Linie passt.

Zweitens kommen manche Gene mit padj = NA zurück. Das ist unabhängige Filterung und Ausreißerentfernung bei der Arbeit, wie im Vignette-Abschnitt „Warum werden einige p-Werte auf NA gesetzt?“ beschrieben. Lassen Sie diese Zeilen vor dem Erstellen des Diagramms weg. Wandeln Sie NA niemals in 1 oder 0 um.

Wo die Grenzwertlinien eingezeichnet werden sollten

Es gibt kein standardmäßiges Grenzwertpaar, nur gängige Paare. Hier sehen Sie, woher die verbreiteten Zahlen tatsächlich kommen:

Wählen Sie anhand der Signalstärke, nicht aus Gewohnheit. Wenn 6,000 Gene durchkommen, ist die Liste unbrauchbar, und ein höherer Fold-Change-Grenzwert hilft. Wenn neun Gene durchkommen, lassen Sie die Fold-Change-Linie weg und sortieren Sie stattdessen nach padj.

Der edgeR-Leitfaden enthält die beste Ein-Satz-Regel für den x-Grenzwert. Lesen Sie ihn als „den Fold Change, unter dem wir definitiv nicht an dem Gen interessiert sind“, nicht als die Änderung, oberhalb derer ein Gen interessant wird.

Der wichtige Grenzwert-Fallstrick

Das Filtern nach p und anschließend nach Fold Change ist kein gültiger Test. Der edgeR-Leitfaden formuliert es unmissverständlich: Solche Kombinationen „sind beide ad hoc und liefern keine aussagekräftigen p-Werte für Tests auf differentielle Expression relativ zu einem Fold-Change-Schwellenwert“. Sie „begünstigen gering exprimierte, aber stark variable Gene“ und „zerstören im Allgemeinen die Kontrolle der FDR“.

Wenn die Fold-Change-Linie für Ihre Aussage wichtig ist, testen Sie gegen sie, statt nachträglich zu filtern:

# 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)

Jetzt berücksichtigt padj den Größengrenzwert bereits, und die vertikalen Linien beschreiben den durchgeführten Test. Beide Optionen sind dokumentiert: siehe Tests des log2-Fold-Changes über einem Grenzwert in DESeq2 und die TREAT-Publikation hinter glmTreat.

Den volcano plot Region für Region lesen

Lesen Sie die Punktwolke, bevor Sie die Beschriftungen lesen. Die Form sagt mehr über das Experiment aus als jedes einzelne Gen.

Oben in der Mitte ist die Region, die am häufigsten falsch gelesen wird. Diese Gene sind meist Gene mit hohen Counts, bei denen der Standardfehler klein ist, sodass selbst eine Änderung von 15% die p-Schwelle überschreitet. Die Änderung ist real, aber ob eine Änderung von 15% relevant ist, ist eine biologische Frage, keine statistische. Ein Warnsignal: Wenn Housekeeping- oder ribosomale Gene dort oben liegen, vermuten Sie ein Problem der Normalisierung oder Bibliothekszusammensetzung, nicht der Biologie.

Die unteren Ränder sind das Spiegelbild. Ein Gen mit 4 Counts in einer Gruppe und 30 in der anderen ergibt einen riesigen log2-Fold-Change und kaum Evidenz. Diese Punkte machen die x-Achse breit und die Geschichte spannend, replizieren sich aber selten. Eine Shrinkage der Fold Changes korrigiert die Darstellung: lfcShrink(dds, coef = 2, type = "apeglm") zieht schlecht gemessene Gene in Richtung null und lässt gut gemessene unverändert.

Zwei Formprüfungen runden die Lektüre ab:

• Symmetrie. Eine asymmetrische Punktwolke, etwa 900 nach oben und 40 nach unten, kann eine echte Aktivierung sein. Sie kann aber auch von einer Verschiebung der Zusammensetzung stammen, die durch die Normalisierung nicht aufgefangen wurde. Sehen Sie sich das MA-Diagramm und die Size Factors an, bevor Sie den Satz schreiben.
• Breite des V. Die Lücke um x = 0, mit Punkten, die an beiden Flanken nach oben steigen, ist zu erwarten, weil größere Änderungen die p-Schwelle leichter überschreiten. Ein durchgehender vertikaler Ausschlag bei x = 0 bedeutet dagegen meist, dass Sie Roh-p-Werte dargestellt haben.

Den volcano plot in R aus DESeq2-Ergebnissen erstellen

Dieses Beispiel läuft mit einem DESeqDataSet namens dds und benötigt ggplot2, ggrepel und 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()

Wenn Sie von edgeR kommen, benennen Sie zwei Spalten um; der Rest des Skripts funktioniert:

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

Für eine vollständig gestaltete Version mit Verbindungslinien und schattierten Grenzwertfeldern lohnt sich ein Blick in die EnhancedVolcano-Vignette.

Den volcano plot in Python aus einem pandas-DataFrame erstellen

Exportieren Sie die DESeq2-Tabelle als CSV und führen Sie anschließend diesen Code aus. Er verwendet ausschließlich pandas, numpy und 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")

Wenn Sie in Python mit edgeR-Ausgaben arbeiten, benennen Sie zuerst logFC in log2FoldChange und FDR in padj um.

Der Code hinter dem obigen Diagramm

Das Beispieldiagramm verwendet simulierte Daten, keinen echten Datensatz. Das Modell ist bewusst einfach: Basismittelwerte sind über vier Größenordnungen verteilt, 5% der Gene tragen einen echten Effekt, und der Standardfehler wächst, wenn der Basismittelwert sinkt. Dieser letzte Teil erzeugt die Punkte mit schwacher Evidenz an den äußeren Rändern.

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

Fünf Fehler, die einen volcano plot irreführend machen

• Roh-p-Werte auf der y-Achse. Hunderte Gene überschreiten die Linie durch Zufall. Verwenden Sie padj oder FDR und nennen Sie in der Achsenbeschriftung, welches Sie verwenden.
• Die Fold-Change-Linie als biologische Bedeutung behandeln. Eine 2-fache Änderung in einem Transkriptionsfaktor und in einem Strukturprotein sind keine vergleichbaren Aussagen. Die Linie ist ein Filter, keine Evidenz.
• Ausreißer die Achsen strecken lassen. Ein Gen bei log2FC = 12 macht alles andere flach. Begrenzen Sie den Ausschnitt mit coord_cartesian() oder set_xlim() und markieren Sie die abgeschnittenen Punkte mit Dreiecken am Rand, statt sie wegzulassen.
• Nach dem Fold Change allein einfärben. Rot für jedes Gen jenseits von x = 1 kennzeichnet verrauschte Gene mit niedrigen Counts als Treffer. Die Farbe braucht beide Bedingungen: jenseits der Fold-Change-Linie und jenseits der padj-Linie.
• Ganz auf gestrichelte Linien verzichten. Ohne sie kann der Leser den Grenzwert nicht aus dem Diagramm erkennen. Zeichnen Sie beide und geben Sie die Zahlen in der Bildunterschrift oder Legende an.

FAQ

Ist ein volcano plot dasselbe wie ein MA-Diagramm?

Nein, sie beantworten unterschiedliche Fragen. Ein MA-Diagramm setzt die mittlere Expression auf die x-Achse und den log2-Fold-Change auf die y-Achse; es zeigt, ob Ihre großen Änderungen von Genen mit niedriger Count-Zahl stammen. Ein volcano plot blendet das Expressionsniveau aus und zeigt stattdessen die Evidenz. Erstellen Sie während der Analyse beide, veröffentlichen Sie den volcano plot.

Was mache ich, wenn padj genau 0 ist?

Das ist ein Gleitkomma-Unterlauf, kein p-Wert von null. -log10(0) ergibt Unendlich, und der Punkt verschwindet oder beschädigt die Achse. Begrenzen Sie padj vor dem Logarithmieren auf den kleinsten positiven Double-Wert, wie im Python-Code oben, und geben Sie das Gen im Text als padj < 1e-300 an.

Wie viele Gene sollte ich beschriften?

Zehn bis zwanzig sind lesbar, mehr werden unübersichtlich. Sortieren Sie nach padj unter den Genen, die auch die Fold-Change-Linie passieren, und beschriften Sie eine feste Anzahl. Verwenden Sie ggrepel in R oder verschobene annotate-Aufrufe in Python, damit der Text nicht auf den Punkten liegt. Beschriften Sie jedes Gen, das Sie im Abstract nennen, auch wenn es nicht unter den ersten zehn ist.

Kann ein volcano plot direkt in eine Publikation?

Ja, wenn er als Vektorgrafik oder in der vom Journal geforderten Auflösung exportiert wird und die zugrunde liegende Analyse reproduzierbar ist. Verlage behandeln Daten-Diagramme anders als Illustrationen, und die Regeln sind vor der Einreichung wichtig — siehe unseren Leitfaden zu Richtlinien von Zeitschriften für KI-generierte Diagramme. Für das Übersichts-Panel einer Publikation deckt ein graphical abstract maker den schematischen Teil ab, während der volcano plot ein Diagramm Ihrer eigenen Zahlen bleibt.