RNA-Seq 火山图：阈值、代码与常见坑

火山图把一次差异表达检验里的每个基因都画在同一张图上。x 轴是 log2 fold change，y 轴是校正后 p 值的 −log10。离中心越远、位置越高的基因，就是你的候选。

画图只要十行代码，读懂它要费些心。多数人照抄的那两条虚线——|log2FC| ≥ 1 与 padj < 0.05——是展示惯例，不是规定。edgeR 用户手册甚至写明，用这一对条件挑基因"总体上会破坏 FDR 控制"（第 2.14 节）。

这篇讲清四件事：两个轴各自的含义、该画哪条阈值线以及为什么、点云各个区域怎么读、以及可直接运行的 R 与 Python 代码。只想要图的话，我们的火山图生成器接受一份结果表，直接输出带标注的图。

基于模拟数据（12,000 个基因，其中 5% 存在真实效应）的火山图。用 matplotlib 绘制，生成代码见下文 Python 一节。

两个轴，以及为什么都取对数

两个轴都取对数，原因相同：原始刻度会把你真正关心的部分藏起来。

原始刻度上的 fold change 是不对称的。翻倍是 2.0，减半是 0.5——向上可以一直到无穷，向下却被压在 1 以下。取 log2 就解决了：翻倍变成 +1，减半变成 −1，图的左右两侧成为镜像。

p 值的问题正好相反。它们全挤在 0 附近，而有意思的基因都在那儿。取 −log10 把这一堆值摊开，负号则让证据最强的点落在最上方。

这张表带出两个事实。y = 4 的点，padj 是 y = 2 那个点的百分之一。而 y 轴没有上界，所以一个 p 值极小的基因就能把整张图拉长。

y 轴要放 padj，不要放原始 p 值

y 轴放原始 p 值，会白送你几百个假阳性。bulk RNA-seq 对每个基因都做检验，一次检验了 15,000 个基因、实际毫无效应的实验，按 p < 0.05 依然会挑出约 750 个基因。

解决办法是做假发现率校正。Benjamini-Hochberg 方法把 p 值重新缩放，使 0.05 这条线的含义变成"我留下的基因里约有 5% 是假的"，而不是"全部检验的 5%"。R 里用 p.adjust 加 method = "BH" 即可，方法出自 Benjamini 与 Hochberg（1995）。

两个主流工具都已经给出了校正列：

• DESeq2（vignette）：log2FoldChange、pvalue、padj。padj 默认就是 BH 校正，因为 results() 的默认参数是 pAdjustMethod = "BH"。
• edgeR（用户手册）：topTags() 返回 logFC、logCPM、一列检验统计量（F 或 LR）、PValue 和 FDR。要用的是 FDR 这一列。

DESeq2 有两个细节最容易踩坑。第一，results() 默认 alpha = 0.1，它的独立过滤步骤是按这个值优化的。如果你的虚线画在 0.05，就该调用 results(dds, alpha = 0.05)，让过滤与你画的线一致。

第二，有些基因返回的 padj 是 NA。按 vignette 里"Why are some p values set to NA?"一节的说明，这是独立过滤和离群值剔除在起作用。画图前把这些行去掉，绝不要把 NA 改成 1 或 0。

阈值线该画在哪

不存在标准阈值，只有常见阈值。那些流行数字的实际出处如下：

按手上信号的强弱来选，别按习惯选。如果 6,000 个基因过线，这份名单没有用，把 fold change 的门槛抬高才有意义。如果只有 9 个过线，那就干脆去掉 fold change 线，直接按 padj 排序。

x 轴阈值的最好一句解释来自 edgeR 手册：把它读成"低于这个 fold change 的基因我们肯定不感兴趣"，而不是"高于这个 fold change 的基因就值得关注"。

一个值得知道的阈值陷阱

先按 p 值筛、再按 fold change 筛，并不构成一次合法的检验。edgeR 手册说得很直接：这类组合"都是权宜之计，就相对某个 fold change 阈值检验差异表达而言给不出有意义的 p 值"。它们"偏向低表达但高变异的基因"，并且"总体上破坏了 FDR 控制"。

如果 fold change 门槛对你的结论很关键，那就直接针对它做检验，而不是事后过滤：

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

这样 padj 本身就已经把幅度门槛算进去了，两条竖线描述的正是你做过的检验。两种做法都有文档：DESeq2 见 tests of log2 fold change above a threshold，glmTreat 背后的方法见 TREAT 论文。

分区读图

先读点云，再读标注。图的形状比任何单个基因都更能说明这次实验的情况。

正上方是最容易被误读的区域。这些通常是高表达基因，标准误很小，所以哪怕只变化 15% 也能过 p 值线。变化是真的。但 15% 的变化有没有意义，是生物学问题，不是统计问题。一个警告信号：如果内参基因或核糖体基因出现在这里，该怀疑归一化或文库组成有问题，而不是生物学发现。

底部两侧是镜像的情况。一个基因在一组里 4 个 count、在另一组里 30 个，log2 fold change 会很大，而证据几乎没有。这些点会把 x 轴撑宽、把故事讲得很动人，却极少能重复出来。用收缩后的 fold change 能修正展示效果：lfcShrink(dds, coef = 2, type = "apeglm") 会把测得不准的基因拉向 0，同时几乎不动测得准的基因。

读图收尾时再做两项形状检查：

• 对称性。 点云一边倒，比如 900 个上调、40 个下调，可能是真实的激活，也可能是归一化没吸收掉的组成偏移。下笔之前先看 MA 图和 size factor。
• V 形的宽度。 x = 0 附近留一个缺口、两侧的点往上爬，是正常形态，因为变化越大越容易过 p 值线。如果 x = 0 处反而立着一根实心尖峰，通常说明你画的是原始 p 值。

用 R 画：从 DESeq2 结果出发

下面这段基于名为 dds 的 DESeqDataSet，需要 ggplot2、ggrepel 和 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()

如果结果来自 edgeR，改两个列名，脚本其余部分照用：

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

想要连线标注、阴影阈值框这类完整样式，可以读一读 EnhancedVolcano vignette。

用 Python 画：从 pandas DataFrame 出发

把 DESeq2 的结果表导成 CSV，然后跑下面这段。只用到 pandas、numpy 和 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")

如果在 Python 里用 edgeR 的输出，先把 logFC 改名为 log2FoldChange、FDR 改名为 padj。

上面那张图的生成代码

示例图用的是模拟数据，不是真实数据集。模型故意做得很简单：base mean 跨四个数量级，5% 的基因带真实效应，标准误随 base mean 下降而增大。正是最后这一条，造出了图上远端那些证据很弱的点。

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

让火山图产生误导的五个错误

• y 轴放原始 p 值。 会有几百个基因纯靠碰运气过线。请用 padj 或 FDR，并在轴标签里写明用的是哪个。
• 把 fold change 线当成生物学意义。 转录因子变化 2 倍和结构蛋白变化 2 倍，不是同一种结论。这条线是筛子，不是证据。
• 让离群点把坐标轴撑爆。 一个 log2FC = 12 的基因就能把其余全部压平。用 coord_cartesian() 或 set_xlim() 裁掉视野，并把被裁掉的点用三角形标在边缘，而不是直接删掉。
• 只按 fold change 上色。 只要 x > 1 就标红，等于把低表达的噪声基因认成命中。上色必须同时满足两个条件：过 fold change 线并且过 padj 线。
• 完全不画虚线。 没有虚线，读者无法从图上看出你的阈值。两条都画上，并把具体数值写进图注或图例。

常见问题

火山图和 MA 图是一回事吗？

不是，两者回答的问题不同。MA 图的 x 轴是平均表达量、y 轴是 log2 fold change，因此能看出大幅变化是否来自低表达基因。火山图舍弃表达量，换成展示证据强度。做分析时两张都画，投稿放火山图。

padj 正好是 0 怎么办？

那是浮点下溢，不是 p 值真等于 0。-log10(0) 得到无穷，点会消失或把坐标轴弄坏。取对数前把 padj 夹到最小正双精度数，做法见上面的 Python 代码，正文里则写成 padj < 1e-300。

该标注多少个基因？

10 到 20 个还能读，再多就乱了。在同时过 fold change 线的基因里按 padj 排序，然后标固定个数。R 里用 ggrepel，Python 里用带偏移的 annotate，别让文字压在点上。摘要里点名的基因一定要标出来，即使它不在前十。

火山图能直接放进论文吗？

可以，前提是导成矢量图或达到期刊要求的分辨率，而且背后的分析可复现。出版社对数据图和示意图的要求不同，这些规则要在投稿前就搞清楚——见我们这篇各出版社对 AI 生成图表的规定。论文的总览图部分可以交给图摘生成器，火山图则始终是你自己数据画出来的图。