打开近期几篇高分期刊的组学文章，你会发现一个有意思的现象：哪怕内容千差万别，图表部分里总有一张结构相似的火山图。横轴是 log2 差异倍数，纵轴是 -log10 校正后 P 值，左上角和右上角散落着蓝点红点，中间铺着一片灰。这张图的统计逻辑不复杂，代码也不超过五十行，但几乎每隔一段时间就会有人来问我同一个问题：多分组的数据，到底该怎么画火山图？
这句话背后藏着一个非常实际的痛点。常规火山图教程几乎全部围绕"两组比较"展开——处理组对对照组、肿瘤对正常、突变对野生型。一旦遇到多分组设计，比如三组、四组甚至带时间序列的分组，很多人立刻就懵了：有人把几组两两组合拆成好几张图，有人直接拿 ANOVA 的 P 值套火山图模板，还有人把所有比较的结果全部叠在一张图上。结果做出来的图，要么信息冗余，要么统计不严谨，投出去直接被审稿人质疑。这篇文章围绕"多分组差异分析火山图"把整条链路讲透，包括多分组差异分析的合理路线、聚合指标的计算逻辑、用 R 实现发表级火山图的完整代码，以及一批常规文档里不会写的踩坑经验。面向正在做转录组、蛋白组或代谢组多组比较的科研人员，目标是解决"如何用一张图讲清楚多组差异"这个问题。
1. 先理解火山图的坐标轴：它到底在表达什么
1.1 横轴与纵轴的统计含义
火山图之所以叫火山图，是因为把差异分析结果画成散点后，显著上调的基因在右上角聚成一个"喷发口"，显著下调的在左上角形成另一个"喷发口"，中间大量无差异的基因像平地一样铺在底部，整体轮廓像一座正在喷发的火山。
横轴的 log2FC 解决的是"变化幅度"的问题。FC（Fold Change）本质上是处理组表达量除以对照组表达量得到的比值，取 log2 之后，上调两倍对应 +1，下调两倍对应 -1，这样让对称的上下调变化在坐标轴上具有相同的视觉距离，避免 2 倍和 0.5 倍在普通线性坐标上不对称的尴尬。纵轴的 -log10(P) 解决的是"变化可信度"的问题。原始 P 值越小，-log10 转换后的数值越大，点在图上的位置越高。一个 P 值为 0.05 的点对应到纵轴约 1.3，一个 P 值为 1e-10 的点对应 10，两者在图上差异极其明显。
实际操作中我经常用一句大白话跟学生解释：横轴量的是"变化有多大"，纵轴量的是"这个变化有多靠谱"。只有幅度大而且可靠的点，才配被点亮成红色或蓝色；幅度大但不可靠的点，只能安静地待在灰色区域里。
1.2 被忽略的细节：纵轴到底放哪个 P 值
这里有一个值得停下来细看的细节：真正决定火山图格局的，往往不是阈值本身，而是你放在横轴和纵轴上的统计量是哪一种。纵轴用原始 P 值还是校正后的 P 值（padj），用 BH 法还是 Bonferroni 法，横轴用普通 log2FC 还是 shrinkage 收缩后的 log2FC，这些选择直接决定你最终能圈出多少基因，也决定整张图的可复现性。
高分期刊普遍要求纵轴使用校正后的 P 值。原因很简单，组学数据动辄对上万个基因做检验，如果不做多重假设检验校正，假阳性数量会非常可观。举个例子，两万个基因全部没有真实差异，用 0.05 作为显著性阈值做两万次检验，平均也会有一千个左右被误判为显著。这批假阳性放到图上就是一堆没有生物学意义的红点，审稿人一眼就能看出问题。所以我的默认配置是：RNA-seq 用 DESeq2/edgeR 自己输出的 padj，芯片或蛋白组学数据用 limma 的 adj.P.Val，作图时再看一眼分布，确认没有出现整体偏移。
2. 多分组差异分析的三种主流路线
2.1 路线一：两两比较，各画各的图
这是最老实、也最容易通过审稿的做法。假设你有三组样本（Control、Treat_A、Treat_B），就跑两次差异分析：Control vs Treat_A、Control vs Treat_B，各自得到一套 log2FC 和 padj，分别画两张火山图。
这种做法统计上挑不出毛病，思路清晰，实验记录也好写，但它有两大实际问题。第一，分组数量一多，图的张数随之爆炸——四组就要六张图，六组就要十五张图，期刊版面根本放不下，整页图全被火山图占满，其他信息无从展示。第二，两两比较之间没有统一基准，读者很难快速抓住"哪些基因在所有比较里都显著"，跨图对比只能靠肉眼来回扫，信息提取效率极低。
所以我的判断是：两两比较适合分组数量少（两组或三组）且关注点集中在一两个特定对比的实验。一旦分组超过三组，就必须考虑下面两种更集约的方案。
2.2 路线二：聚合指标，一张图讲完所有比较
为了在一张图里呈现多组比较的全局信息，很多高分论文采用"聚合策略"：对每个基因，取所有两两比较中绝对 log2FC 的最大值，以及所有比较中校正后 P 值的最小值，然后用这两个聚合指标画火山图。
这个策略背后的逻辑非常直白：只要有一个比较里该基因达到了显著且变化幅度足够大，这个基因就值得被标记出来。聚合图的优点是信息密度高，适合回答"哪个基因在整体上最值得关注"这类问题；缺点是你必须接受信息压缩——图上的点不再对应单一比较，而是多个比较的"上包络"。实际使用中，我会在论文方法部分明确写清楚聚合规则，并在图注里注明"FC 取最大绝对值，P 值取最显著值"，审稿人对这种做法的接受度很高。我见过不少发表在顶级期刊上的文章，方法部分就是一两句话带过，审稿人并不会因此为难作者。
2.3 路线三：全局检验，先筛再定位
还有一种思路在许多资深生信工程师那里很受欢迎：先做全局检验，把组间存在总体差异的基因筛出来，再对筛出来的基因做两两比较，用两两比较中的最大 log2FC 画火山图。
这种做法的好处是逻辑链条完整——先问"这个基因在任意两组之间是否存在显著差异"，再问"具体差多少、方向如何"，相当于给火山图加了一道统计关卡。RNA-seq 的 count 数据通常用 edgeR 的 F 检验或 DESeq2 的 LRT（似然比检验）来做全局筛选，而不是对标准化后的表达量直接跑普通 ANOVA，这一点很多教程没有讲清楚。原因在于 count 数据服从离散分布且方差随均值变化，直接套用适用于正态连续数据的 ANOVA，会严重高估统计功效。芯片和蛋白组学数据相对接近连续分布，用常规 ANOVA 或者其非参数版本 Kruskal-Wallis 是可以接受的，但样本量很小时要特别谨慎，方差齐性假设很容易被打破。
下面这张表汇总了三条路线的适用场景，方便你快速对号入座。
| 路线 | 统计逻辑 | 适合场景 | 图表形式 | 审稿接受度 |
|---|---|---|---|---|
| 两两比较 | 每个比较独立检验 | 分组少（2-3组）、关注具体对比 | 每个比较一张火山图 | 高 |
| 聚合指标 | 取最大绝对FC与最显著P | 分组多、版面有限、找核心基因 | 单张聚合火山图 | 高（需注明规则） |
| 全局检验后两两 | 先全局筛选再定位差异 | 分组复杂、需要全局扫描 | 单张或多张火山图 | 中高 |
3. 实操准备：从表达矩阵到可靠的差异结果表
3.1 数据格式与分组信息的常见坑
先聊数据准备。无论用哪种差异分析工具，你手里至少要有两样东西：一是表达矩阵，基因在行、样本在列；二是分组信息表，至少包含样本名和组别两列。表达矩阵的来源不同，预处理要求也不一样。RNA-seq 的 count 矩阵可以直接交给 DESeq2 或 edgeR 处理，芯片数据或蛋白组学定量数据则要先做背景校正、归一化，再交给 limma，否则后续的 log2FC 和 P 值都不可靠。
分组信息表最常见的坑是样本顺序和表达矩阵列名不一致。我处理过不止一次这种情况：Excel 里手工整理分组表时，样本顺序跟矩阵列名错位，merge 之后组别标签张冠李戴，整个差异分析的结果全部作废。写了多次脚本之后，我的固定做法是：读入矩阵后先把列名排序，再用 merge 按样本名关联分组信息，并且用 factor 显式指定组别水平的顺序。R 里 factor 默认按字母序排列，如果不显式指定，Control 组可能被排在 Treatment 后面，后面所有比较基准都会错位，图上标签也跟着乱。
3.2 差异分析的执行要点与结果表结构
差异分析算完，结果表里通常包含这些列：gene_id、logFC、AveExpr、t 统计量（或 z 统计量）、PValue、FDR/padj。多分组两两比较时，我会把每一组比较的 logFC 和 padj 单独命名，再合并成一张"宽表"，方便后面做聚合和画图：
merge_res <- Reduce(function(x, y) merge(x, y, by = "gene_id"), list(res_c_vs_t1, res_c_vs_t2, res_c_vs_t3))merge 的时候要特别注意基因 ID 的去重。如果原始矩阵里有重复的基因符号，merge 会把行数悄悄放大，画图时同一个基因出现多个点，乍看没什么异常，实际上统计的样本数据和基因数目对不上，投稿后被要求提供源数据时会非常被动。建议在差异分析之前就做好基因注释去重，保留表达量最高的转录本或取均值，把重复问题解决在最前面。
另外，DESeq2 的 results() 函数在多分组设计里默认只输出一个比较方向的结果，很多人在这里栽过跟头。正确做法是先用 resultsNames() 查看所有可用的比较名称，再用 contrast 参数显式指定你要提取哪一组对比。养成差分析完第一件事先打印 resultsNames() 的习惯，能避开一大批隐性错误。
4. 画图全流程：从聚合指标到发表级火山图
4.1 聚合指标怎么算才不出错
拿到宽表之后，第一步是计算两个聚合量：max_abs_logFC 和 min_padj。这里有一个非常容易踩的细节：取 min_padj 时不能直接对包含 NA 的向量用 pmin()，因为只要有一个比较里 padj 是 NA（比如某个基因在某一组样本里完全不表达），pmin 会直接返回 NA，整行基因都会被丢弃，导致图上的基因数量明显偏少。
稳妥的做法是先把所有比较的 padj 列里的 NA 替换成 1（1 在 -log10 转换后等于 0，视觉上等同于不显著），再取最小值。同样，max_abs_logFC 那边也要带上 na.rm = TRUE。完整代码如下：
library(dplyr) final_res <- merge_res %>% mutate( max_abs_logFC = pmax(abs(logFC_c_vs_t1), abs(logFC_c_vs_t2), abs(logFC_c_vs_t3), na.rm = TRUE), min_padj = pmin(replace(padj_c_vs_t1, is.na(padj_c_vs_t1), 1), replace(padj_c_vs_t2, is.na(padj_c_vs_t2), 1), replace(padj_c_vs_t3, is.na(padj_c_vs_t3), 1)), direction = case_when( max_abs_logFC >= 1 & min_padj < 0.05 ~ "up", max_abs_logFC <= -1 & min_padj < 0.05 ~ "down", TRUE ~ "ns" ) )这里把 |log2FC| ≥ 1（即差异倍数 ≥ 2）和 padj < 0.05 当作默认阈值。阈值不是死规矩，完全可以根据实验性质调整：药物处理实验差异通常很大，阈值可以收紧到 |log2FC| ≥ 2；临床样本异质性大，可以放宽到 |log2FC| ≥ 0.58（即 1.5 倍差异）。关键是阈值必须在方法部分交代清楚，不要画完图再倒推一个看起来好看的数字，那是数据造假的前奏。
4.2 ggplot2 火山图主体代码
ggplot2 画火山图是绝对主流，配上 ggrepel 做基因标签，整套流程非常成熟。核心代码并不复杂：
library(ggplot2) library(ggrepel) p <- ggplot(final_res, aes(x = max_abs_logFC, y = -log10(min_padj))) + geom_point(aes(color = direction), size = 1.8, alpha = 0.75) + scale_color_manual( values = c("up" = "#D32F2F", "down" = "#1976D2", "ns" = "#BDBDBD"), name = "Significance" ) + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") + geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey40") + labs(x = "Max |log2(Fold Change)|", y = "-log10(adjusted P-value)") + theme_classic(base_size = 14) + theme(legend.position = "top")几个细节值得展开说。第一，geom_point 的 size 和 alpha 要配合点总数来调。基因数在两万左右时，size 1.5 到 2、alpha 0.6 到 0.8 是比较稳的区间；点数超过三万，建议先对 ns 类别做下采样，否则中间灰色区域会密集成一片黑，根本看不出点的疏密变化。第二，配色不要直接用 Python matplotlib 的默认亮色，期刊打印出来容易失真。我长期用 #D32F2F 红色、#1976D2 蓝色、#BDBDBD 灰色这套 Material Design 配色，在白色背景上对比度足够，也相对色盲友好。第三，坐标轴标签一定写清楚是 Max |log2(FC)|，而不是笼统的 log2FC，否则读者会误以为这是单一比较的结果，研究方法部分前后对不上。
4.3 基因标注与标签防重叠
一张发表级火山图，通常只会在显著的基因里挑一部分标上基因名。全部标注会糊成一团，我的做法是：优先标注 |log2FC| 最大或 padj 最小的前 10 到 20 个基因，如果论文有关注的特定基因（比如通路核心成员、前期验证过的候选基因），再单独用 ggrepel 强制标注。
top_genes <- final_res %>% filter(direction != "ns") %>% arrange(min_padj) %>% head(15) p <- p + geom_point(data = top_genes, aes(x = max_abs_logFC, y = -log10(min_padj)), color = "black", size = 2.2, shape = 21, stroke = 0.5) p_labeled <- p + geom_text_repel(data = top_genes, aes(label = gene_id), size = 3.2, max.overlaps = 20, segment.color = "grey50", segment.size = 0.3)要注意的是，ggrepel 的 max.overlaps 参数在不同 ggplot2 版本里的默认值不一致，不显式设置时可能莫名其妙丢标签，尤其当你把数据过滤后重新绘图，标签消失得悄无声息。segment 连接线颜色用灰色、线宽 0.3 左右最自然，太粗会抢散点的视觉权重。标签字号在最终导出 180 mm 宽的单栏图里对应 5 到 7 pt 比较合适，字体大了显得业余，太小在打印稿里看不清。
5. 常见翻车现场与排查方法
5.1 显著基因过多，整张图一片红海
这种情况十有八九是纵轴用了原始 P 值，或者没有过滤低表达基因。低表达基因的 count 数很低，差异倍数波动极不稳定，会产生大量虚假显著。RNA-seq 分析前用 edgeR 的 filterByExpr 或手动保留在至少一组样本中 CPM 大于 1 的基因，能消掉一大部分噪音。如果过滤后还是红点多，可以考虑把阈值从 padj 0.05 收紧到 0.01，同时给纵轴设上限（cap），避免个别 P 值小到 1e-300 的点把纵轴拉到失真。
5.2 图边缘有基因"飞出去"
当某个基因的 log2FC 特别大（比如基因敲除后完全不表达），点会直接冲出绘图区域。审稿人不会喜欢这种图。我的做法是定义坐标轴范围，同时对超出范围的基因做截断标记。用 scale_x_continuous(limits = c(-8, 8)) 之前，一定先看看数据的真实分布，硬截断会导致图内点的数量与统计结果不一致，最好在图注里注明截断范围。有些人会用 coord_cartesian 来缩放，这个函数不会删点，只是改变显示区域，比直接 limits 更安全。
5.3 聚合时方向信息丢失
聚合指标里藏着一个隐患：取绝对值最大 logFC 会让"上调"和"下调"信息变成单一方向。比如基因 X 在比较 1 中上调 3 倍，在比较 2 中下调 2.5 倍，max_abs_logFC 是 3，按上面的分类逻辑会被标成 up，但这个基因在不同比较里的方向其实并不一致。这在生物学上恰恰是值得注意的现象，聚合图却把它掩盖了。我的处理方式是在聚合表里额外生成一列 direction_consistency，如果所有显著比较的方向一致才标 up/down，方向冲突的标为 conflict，用第三种颜色（比如紫色）在图中单独标出。这样既保留了聚合图的简洁，又不丢失多组比较特有的矛盾信息。
5.4 分组顺序导致比较基准错乱
前面提过 factor 顺序的问题，这里再补充一个具体例子。三组样本名称分别是 ctrl、treatment、recovery，字母序排列是 ctrl、recovery、treatment。如果直接用默认排序，你的"recovery vs ctrl"会被当成"recovery vs treatment"来解读，结果完全错位。DESeq2 的 results() 函数默认只输出一个比较方向，多组时要用 contrast 参数显式指定。养成跑完差异分析先打印 resultsNames() 的习惯，能省掉大量返工时间。
5.5 输出格式与分辨率的最后一道坎
高分期刊对图片格式有硬性要求：位图至少 300 dpi，线图和散点图优先矢量格式。ggplot2 的 ggsave 可以一次性满足两个要求：
ggsave("volcano_multi_group.pdf", p_labeled, width = 180, height = 150, units = "mm") ggsave("volcano_multi_group.tiff", p_labeled, width = 180, height = 150, units = "mm", dpi = 300, compression = "lzw")PDF 是矢量格式，无论放大到多少都保持清晰；TIFF 用于投稿系统强制要求位图的场景，300 dpi 是标配。宽度 180 mm 对应单栏或 1.5 栏宽度，高度 150 mm 保持比例协调。注意导出前把图表里的中文字体统一替换成英文——ggplot2 默认主题在 PDF 里嵌入中文字体会出现字体警告，严重的甚至导致生成的 PDF 文件打不开，这个问题在 Windows 系统上尤其常见。
6. 从图表到高影响力期刊：投稿前的自查清单
最后聊一点务虚但很关键的内容。IF=33.2 这个级别的期刊，类型多半是 Nature 系列大子刊或 Lancet 系列，这类期刊对图表的要求不是"华丽"，而是"信息准确、自解释、经得起审稿人逐项挑刺"。
一张火山图能不能达到这个标准，我会按下述清单逐项自查：
- 横轴纵轴的统计量名称是否准确，log2FC 的底数、P 值的校正方法是否在方法部分有交代；
- 阈值线是否与正文方法部分一致，红色和蓝色点的数量是否与文中报告的差异基因数目完全吻合；
- 图注里是否写清楚多组聚合规则（取最大绝对 FC、最显著 P），是否说明图中每个点代表一个基因；
- 配色是否照顾色盲读者，打印成灰度图后 up/down/ns 三类点是否仍然可区分；
- 分辨率是否满足 300 dpi 或矢量格式，图片在双栏排版缩到 80 mm 宽度后标签是否依然可读。
其中最容易忽略的是第二点。很多人在图上用颜色标了 up 和 down，但全文从头到尾没有明确说"我们定义 |log2FC| ≥ 1 且 padj < 0.05 为显著差异基因"，审稿人只能自己猜，猜错了就是一轮 major revision。建议把阈值定义直接写进图注，或者至少在图注里引导读者去看方法部分的对应段落，这是投入产出比最高的一个动作。
多分组聚合火山图本质上是一个信息压缩与信息保真的权衡。压缩做得太好，图很漂亮但丢了细节；保真做得太多，图就乱成一片。我个人的经验是先画一个"全信息版"（每个比较单独出图）自己核对基因方向和显著状态，再画"聚合版"给读者看。两个版本在分析报告里都保留，投稿时根据期刊偏好选择，图自然经得起推敲。每次做完一张图，顺手把聚合规则和阈值记录在 R 脚本的注释里，三个月后返修时你一定会感谢当时的自己。
