跳到正文
百沐一下百沐一下
登录
本节目录

火山图：把差异分析结果画成一张经得起审稿的图

跟完这篇，您能把一张差异分析结果表，一步步画成读得懂、经得起审稿的火山图：两个轴各放什么、阈值怎么定、上下调怎么配色、关键基因怎么标，最后导出成投稿能用的图片文件——中间常见的几个坑也一并避开。

火山图（volcano plot）是差异表达分析最常见的一张图。它把成千上万个基因压进一张散点图里，一眼就能看出「哪些基因变了、变了多少、有多可信」。工具不复杂，但要画对，得先想清楚两个轴。

火山图的两个轴

火山图是一张散点图，每个点是一个基因。两个轴分别回答两个问题：

轴变量回答的问题越往哪边
横轴 Xlog2FoldChange（log2 倍数变化）变了多少、往哪个方向变越往右上调越强，越往左下调越强
纵轴 Y-log10(padj)（校正 p 值取负对数）这个变化有多可信越往上越显著

两点需要说清楚：

  1. 横轴用 log2 而不是原始倍数，是为了让「上调 4 倍」和「下调 4 倍」在图上对称（分别落在 +2 和 −2），零点表示没有变化。
  2. 纵轴取 -log10，是因为 p 值越小越显著，直接画 p 值会把最显著的点压到最底下。取负对数之后，padj = 0.01 变成 2、padj = 0.001 变成 3，越显著的点越高，符合直觉。

纵轴用的是校正后的 p 值（padj，也叫 FDR、adjusted p-value），不是原始 p 值。成千上万个基因同时检验会带来大量假阳性，必须做多重检验校正（DESeq2 默认用 Benjamini-Hochberg 法）。画火山图请优先用 padj。

准备一张差异结果表

火山图的输入就是您跑完差异分析得到的那张结果表，至少要有两列：倍数变化和校正 p 值。不同工具的列名不一样，先对照一下：

工具倍数变化列校正 p 值列
DESeq2（results()）log2FoldChangepadj
edgeR（topTags()）logFCFDR
limma（topTable()）logFCadj.P.Val

本文按 DESeq2 的列名往下写。若您用的是 edgeR 或 limma，把代码里的列名换成对应的即可。先把结果表读进来，并把基因名单独放一列（后面标注要用）：

# 若结果已存成 CSV（DESeq2 / edgeR / limma 导出的都行）
res <- read.csv("deg_results.csv", row.names = 1)
res$gene <- rownames(res)          # 把基因名放到单独一列

# 若刚跑完 DESeq2，也可以直接转成数据框：
# res <- as.data.frame(results(dds))
# res$gene <- rownames(res)

数据从哪来、结果表各列什么意思，见 用 DESeq2 做差异表达分析。

给每个基因贴上下调标签

判断一个基因「值得关注」，通常同时看两件事：显著（padj 足够小）且变化幅度足够大。所以先定两个阈值，再按阈值把每个基因分成三类。

library(dplyr)

p_cut   <- 0.05   # 显著性阈值：常用惯例，按数据调，不是铁律
lfc_cut <- 1      # log2FC 阈值：常用惯例，按研究问题松紧

res <- res %>%
  mutate(change = case_when(
    padj < p_cut & log2FoldChange >  lfc_cut ~ "Up",     # 显著上调
    padj < p_cut & log2FoldChange < -lfc_cut ~ "Down",   # 显著下调
    TRUE                                     ~ "NS"       # 其余算不显著
  ))

table(res$change)   # 看看三类各有多少个

这两个阈值也决定了待会儿三条虚线画在哪：

  • 一条水平线在显著性阈值处：y = -log10(0.05) ≈ 1.30；
  • 两条竖直线在倍数变化阈值处：x = -1 和 x = 1（即上调或下调 2 倍以上）。

阈值是惯例，不是定律。 padj < 0.05 与 |log2FC| > 1 是文献里最常见的默认，但都不是硬标准：样本量大、信噪比高时可以更严（如 padj < 0.01）；探索阶段、差异本就微弱时可以放宽 log2FC。关键是在图注或正文里写清您用的阈值，让读者知道线画在哪。padj 为 NA 的基因（DESeq2 独立过滤后置空的）在这里会自然落进 NS，正是我们想要的。

画出带配色和阈值线的火山图

配色的作用是把三类点一眼分开：显著上调、显著下调、其余不显著。通行惯例是红上、蓝下、灰不显著。

类别含义常用颜色
Up显著上调红色
Down显著下调蓝色
NS未达阈值灰色
library(ggplot2)

cols <- c(Up = "#E64B35", Down = "#3182BD", NS = "grey80")

p <- ggplot(res, aes(x = log2FoldChange, y = -log10(padj), color = change)) +
  geom_point(alpha = 0.6, size = 1.2) +
  scale_color_manual(values = cols) +
  geom_hline(yintercept = -log10(p_cut), linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = c(-lfc_cut, lfc_cut), linetype = "dashed", color = "grey40") +
  labs(x = "log2 fold change", y = "-log10(padj)", color = NULL) +
  theme_bw()

p

几个让图更干净的小处理：

  • 让灰点退到后面做背景，红蓝的信号点才跳得出来：给 NS 更低的透明度或更小的点即可。
  • 点很多时降低透明度（alpha = 0.4 左右），避免中间黑成一团、看不出密度。
  • 红蓝别用最刺眼的纯红纯蓝，稍微降一点饱和度更耐看，也更接近期刊常见配色。
  • 画图时若看到 Removed N rows containing missing values 的提示，多半是 padj = NA 的基因被丢掉了——它们本就归入不显著，不影响主图。

标注要讲的关键基因

火山图上不需要、也不应该标出每一个基因的名字——只标少数几个真正要讲的。标注用 ggrepel 包，它会自动把文字错开、避免重叠。挑选方式常见两种，按需选一种：

  1. 按生物学意义手动选：您研究里关心的那几个基因，直接列名字最稳妥。
  2. 按统计量自动选 Top N：例如取「显著且 log2FC 最大」的前若干个。
library(ggrepel)

# 方式一：手动指定要标的基因
genes_to_label <- c("TP53", "MYC", "CDKN1A")
lab <- dplyr::filter(res, gene %in% genes_to_label)

# 方式二：自动取显著变化里最靠边的 10 个
# lab <- res %>% filter(change != "NS") %>% slice_max(abs(log2FoldChange), n = 10)

p_lab <- p +
  geom_text_repel(
    data = lab,
    aes(label = gene),
    size = 3, max.overlaps = 20,
    box.padding = 0.4, min.segment.length = 0,
    show.legend = FALSE            # 避免图例里多出一个「a」
  )
p_lab

标注的分寸：宁少勿多。一张图标 5～15 个名字通常就够讲清故事；标几十上百个，文字会糊成一片，反而什么都看不清。

导出成出版级图片

最后把图存成文件。投稿优先存矢量图（PDF / SVG，放大不糊）；预览或贴进文档用位图（PNG，记得给足分辨率）。

ggsave("volcano.pdf", p_lab, width = 6, height = 5)              # 矢量图，投稿首选
ggsave("volcano.png", p_lab, width = 6, height = 5, dpi = 300)   # 位图，预览或贴进文档

width / height 单位默认是英寸，按目标版面调；位图的 dpi 一般给到 300 以上才够清晰。

想更省事，可以用 Bioconductor 的 EnhancedVolcano 包：一行 EnhancedVolcano(res, lab = res$gene, x = "log2FoldChange", y = "padj") 就能出一张自带阈值线和标注的火山图。代价是细节不如手写 ggplot 灵活，微调配色、字号时仍要翻它的参数。

常见坑与自查

出图前对照这张表快速自查一遍：

错误后果怎么避免
两个轴放反（把 log2FC 放纵轴）图彻底读不对，方向和显著性混在一起记牢：横轴 log2FC、纵轴 −log10(padj)
忘了取负对数（纵轴直接放 padj）最显著的点被压到最底部，图上下颠倒纵轴永远是 -log10(padj)
用原始 p 值而非 padj假阳性没被控制，显著数被高估优先用校正后的 padj / FDR
不画、不标阈值线读者不知道「显著」的界在哪画出虚线，并在图注写清阈值
一次标几十个基因文字重叠成团，喧宾夺主只标要讲的几个，5～15 个为宜
没处理 padj = NA 或 = 0点消失或飞到无穷高NA 归入不显著或先剔除；padj = 0 时可设一个下限再取对数

关于 padj = 0：个别极显著的基因 padj 会被记为 0，-log10(0) 是无穷大，点会飞出画面。稳妥做法是给它设个下限，例如 res$padj[res$padj == 0] <- min(res$padj[res$padj > 0]) / 10，再取对数。关于 NA：DESeq2 会对表达量过低的基因做独立过滤，把它们的 padj 置为 NA，画图时应剔除或归入不显著，不要当成缺失去插补。

在平台上跑一遍

不想写 R 代码也能得到同样的图：在 百沐一下 传入一张差异分析结果表（含 log2FoldChange 与 padj 两列），用一句话说清「画火山图、阈值 padj<0.05 且 |log2FC|>1、标出这几个基因」，就会给您一张出版级火山图，连同可复现的 R 代码一起交付，配色和阈值都能对话式再调。操作细节见 智能绘图使用指南。

相关