跟着论文复现一张火山图

跟完这篇，您能挑一篇公开了数据的转录组论文，把它的原始数据取回来、按它的口径重跑一遍差异分析，最后画出一张和原文趋势一致的火山图，并在对不上的时候知道去哪里找原因。

复现别人的一张图，是练生信最扎实的方式。它逼您把「论文里一句话带过的方法」拆成一步步能照做的动作，也让您第一次直面「为什么我的结果和人家的不完全一样」。这篇就带您完整走一遍，用一张火山图当靶子。

整件事拆成五步，跟着往下走即可：

读懂论文的数据来源与阈值口径；
取回对应的公共数据；
按论文口径重跑差异分析；
照着原图复刻火山图；
对不上时，逐项排查原因。

先读懂论文用了什么数据、什么阈值口径

动手之前，先把论文的方法部分（Methods）读透。复现能不能对上，八成取决于您是否抄准了下面这几项。建议边读边填一张表：

要素	在论文哪里找	为什么关键
数据编号	Data availability / 附录，形如 `GSExxxxx`、`PRJNAxxxxx`	决定您去哪儿取数据
样本与分组	Methods、图注、附表	哪些样本算处理组、哪些算对照，直接决定差异方向
比较方向	Methods、图注（如 treated vs. control）	火山图的上下调是相对谁而言
差异分析工具	Methods，如 DESeq2 / edgeR / limma-voom	不同工具的 log2FC 和 p 值口径不同
显著性阈值	Methods 或图注（如 `padj < 0.05`）	火山图那条水平虚线画在哪
倍数变化阈值	Methods 或图注（如 `log2FC` 绝对值 > 1）	那两条竖直虚线画在哪
参考基因组 / 注释版本	Methods（如 GRCh38、Ensembl release 104）	影响基因 ID 与表达量，是复现对不上的常见根源

几个读论文时容易漏、却很要命的细节：

阈值到底是 padj 还是原始 pvalue。 有的论文写 p < 0.05 其实指校正后的值，有的真的是原始 p 值，口径差一个字，显著基因数能差出好几倍。拿不准时以「校正后」为准更稳妥，同时在复现说明里注明您的选择。
有没有对 log2FC 做收缩（shrinkage）。 很多 DESeq2 流程画火山图前会用 lfcShrink()，横轴会因此明显「收窄」。论文没写不代表没做，这是后面对不上的高发点。
样本有没有被剔除。 论文常在正文或附表里提一句「排除了 N 个质控不合格样本」，漏掉这句，您的样本数就和人家不一样了。

复现的目标是趋势一致、结论可重现，不是像素级一模一样。只要上下调的主要基因、显著数量的量级、火山图的整体形状对得上，就算复现成功。追求每个点分毫不差，多数时候既做不到、也没必要。

取到对应公共数据

拿到数据编号后，去公共数据库把数据取回来。转录组数据最常见的落脚点是 GEO（GSExxxxx）和它背后的 SRA（原始测序 reads）。这里有一条省时的岔路，先判断走哪条：

论文只给了原始 reads（fastq）：您需要自己从头比对、定量，工作量大，也最容易引入版本差异。
论文在 GEO 附件里放了处理好的表达矩阵（Supplementary file，常是一个 counts 或 TPM 表）：直接下这个，能跳过比对定量，复现难度骤降。优先走这条路。

判断和下载的思路：

打开 GEO 上该 GSE 的页面，先看页面底部的 Supplementary file。若有一个明显是表达矩阵的文件（名字里带 counts、matrix、expression），先取它。
若只有原始数据，才走 SRA，把 SRR 编号对应的 fastq 取下来，再自己比对定量。检索与下载的完整做法见在 GEO 中检索并获取公共数据。
数据取回后，务必核对样本数、分组标签是否和论文一致。GEO 的样本描述（characteristics）里通常写了每个样本属于哪一组，照着它把分组表建好。表达矩阵怎么整理成规范形态、分组表怎么建，见准备一份规范的表达矩阵。

# 若论文提供了处理好的 counts 矩阵，通常就是一个可直接读入的压缩表
# 具体文件名以 GEO 页面为准，这里只示意形态
gunzip GSExxxxx_raw_counts.tsv.gz
head GSExxxxx_raw_counts.tsv    # 先看一眼：行是基因、列是样本？

一个务实的提醒：能用论文自带的 counts 矩阵，就别从 fastq 重头跑。 从原始 reads 复现要对齐比对软件、参考基因组、注释版本一大串参数，任何一处不同都会让结果漂移。用现成矩阵，您复现的是「从 counts 到火山图」这一段，反而更干净、更可控。

按论文口径跑差异

数据和分组就位后，按论文写明的工具和口径跑差异分析。如果论文用的是 DESeq2，就照它的比较方向和阈值来，别自作主张换工具或换阈值——复现阶段，忠实还原比「跑得更好」更重要。

library(DESeq2)

# counts：原始整数矩阵（行基因、列样本）；colData：样本分组表，行名与 counts 列名对齐
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = colData,
  design    = ~ condition
)

# 参考水平设成论文里的对照组，log2FC 方向才和原文一致
dds$condition <- relevel(dds$condition, ref = "control")

dds <- DESeq(dds)

# contrast 写清「谁比谁」，方向和论文图注保持一致
res <- results(dds, contrast = c("condition", "treated", "control"))

# 若论文画火山图前做了收缩，这里补上，横轴范围才对得上
# apeglm 是 DESeq2 常用的收缩方法之一，按论文所用的方法选
# res <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")

跑之前，先对着上一节那张表把三件事核对一遍：

比较方向：contrast 里处理组和对照组的先后，决定上调下调的正负。写反了，火山图会左右镜像。
是否收缩：论文若画火山图前做了 lfcShrink()，您也要做，否则横轴范围对不上。
阈值：把论文的 padj 线和 log2FC 线记下来，下一步画图直接用。

DESeq2 建对象、跑流程、读结果表、定阈值的完整讲解，见用 DESeq2 做差异表达分析。若论文用的是 edgeR 或 limma-voom，思路一样，只是换成对应工具的标准流程。

阈值 padj < 0.05、log2FC 绝对值 > 1 是文献里最常见的默认，属于常用惯例、按数据调，不是铁律。但在复现语境下，您不需要纠结哪套阈值「更好」——照抄论文的那套就对了。阈值本身合不合理，是原作者的判断，不是复现要解决的问题。

复刻它的火山图配色与标注

差异表结果一出，就能照着原图把火山图搭起来。复刻一张火山图，核心是对齐四样东西：两个轴、阈值虚线、上下调配色、标注的基因。

library(ggplot2)
library(dplyr)

p_cut   <- 0.05   # 抄论文的显著性阈值
lfc_cut <- 1      # 抄论文的 log2FC 阈值

res <- as.data.frame(res) %>%
  mutate(change = case_when(
    padj < p_cut & log2FoldChange >  lfc_cut ~ "Up",
    padj < p_cut & log2FoldChange < -lfc_cut ~ "Down",
    TRUE                                      ~ "NS"
  ))

ggplot(res, aes(log2FoldChange, -log10(padj), color = change)) +
  geom_point(alpha = 0.6, size = 1.2) +
  scale_color_manual(values = c(Up = "#E64B35", Down = "#3182BD", NS = "grey70")) +
  geom_hline(yintercept = -log10(p_cut), linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = c(-lfc_cut, lfc_cut), linetype = "dashed", color = "grey40") +
  theme_bw()

对着原图逐项校准：

配色：多数论文上调用红、下调用蓝、不显著用灰。若原图用了别的配色（比如上调橙、下调绿），把 scale_color_manual 里的颜色换成对应值即可。
标注哪些基因：原图标了名字的那几个，就在您的图上标同样的几个，别自己另选。用 ggrepel 的 geom_text_repel() 把它们错开。
轴范围：若原图横轴收得比您窄，多半是原文做了 lfcShrink()，回上一步补上收缩即可。

火山图两个轴各放什么、阈值线怎么画、配色和标注的完整讲究，见火山图：差异基因可视化的标准做法。

对不上怎么排查

第一次复现，几乎不会一次对齐——这很正常。别急着怀疑自己，按下面几类原因从大到小逐个排除。它们对结果的影响，通常也是从上往下递减：

可能原因	典型表现	怎么查、怎么办
分组或方向错了	上下调整体反了，或显著基因几乎不重合	回头核对 `contrast` 方向、每个样本的分组标签是否和论文一致
样本没对齐	样本数和论文不同	查论文有没有剔除样本、有没有多个批次只用了其中一批
注释 / 基因组版本不同	基因数量、部分基因 ID 对不上	看论文用的参考版本（如 GRCh38 与 GRCh37、Ensembl release 几），尽量对齐
定量软件不同	从 fastq 复现时 counts 系统性偏移	若论文用 featureCounts、您用了 Salmon，counts 会不同；能用论文的矩阵就别自己重定量
收缩与否	显著基因大致对得上，但横轴范围 / 极端点不同	确认论文画图前是否做 `lfcShrink()`，补上或去掉以对齐
软件版本差异	结果非常接近但仍有细微出入	记下您的 DESeq2 版本；不同版本默认行为（如独立过滤、收缩方法）可能微调
阈值口径	显著基因数量级差很多	复查论文的线是 `padj` 还是原始 `pvalue`、倍数阈值是 1 还是 1.5

排查时的两条实用心法：

从影响最大的往下查。 分组、方向、样本这类「结构性」错误，会让结果面目全非；软件小版本差异只带来毫厘出入。先排大的，再抠小的。
记录您的环境。 跑完 sessionInfo() 把 R 和各个包的版本存下来。既方便自己回溯，也是将来别人复现您结果时的依据。

sessionInfo()   # 记下 R 与 DESeq2、ggplot2 等的版本号

诚实地说：从原始 fastq 出发，想和几年前的论文分毫不差，往往做不到——当年的软件版本、参考注释可能已经很难完整还原。这不是您做错了。复现追求的是「主要差异基因和结论能重现」，而不是复刻每一个小数点。把差异查清、能解释清楚，本身就是有价值的产出。

在平台上跑一遍

上面这套「取数据 → 跑差异 → 画火山图」，在百沐一下不写代码也能走完：上传论文提供的表达矩阵和分组表，用一句话说清「按 treated vs. control 跑 DESeq2、阈值照论文用 padj < 0.05 且 log2FC 绝对值 > 1、画火山图并标出这几个基因」，平台会零代码替您跑完差异分析、画出火山图，并连同可复现的 R 代码与用到的包版本一起交付，方便您和原文逐项比对。完整流程见组学分析助手使用指南。