1. 差异基因分析到底在解决什么问题做转录组的人绕来绕去最后都要落到同一个问题上哪些基因在实验组和对照组之间变了。差异基因分析就是回答这个问题的标准动作而火山图是把结果一次性摊开给你看的图。标题里说一秒学会我理解这种说法背后是很多新手被一堆报错和参数劝退的无奈——真正上手过的人都知道画图那一步确实快难的是前面把数据和对齐逻辑准备利索。我这次就按先把原理讲透再给能直接抄的代码这个路数来写。适合谁看如果你手头有一张表达矩阵counts 表或者 TPM/FPKM 表有分组信息想拿到差异基因列表并且想要一张能放进论文或者汇报 PPT 的火山图那这篇基本能覆盖你从零到出图的全流程。哪怕你之前没怎么用过 R只要照着步骤走也能跑通。先说清楚这两件事的定位。差异基因分析本质是一次带统计约束的对比它不是简单地算个比值就完事而是要考虑基因表达本身的波动、测序深度差异、生物学重复的离散程度。火山图则是把每个基因的两个核心指标——变化倍数和显著性——同时画在二维平面上。横轴放变化幅度纵轴放显著性越往两边越靠上的点就是越值得你关注的基因。理解了这一点后面所有的参数选择、代码调整都有了落脚点不会变成盲人摸象。2. 分析前的整体设计与关键参数选型2.1 从一张表达矩阵说起一切从表达矩阵开始。行是基因列是样本格子里的数字是该基因在该样本里的表达量。这个矩阵通常来自上游的比对和定量流程比如用比对工具得到每个基因的 read 计数。这里有个新手最容易忽略的点你必须明确手里的数值到底是原始 count 还是已经归一化过的表达量。这决定了你后面用哪套差异分析工具。如果拿 FPKM 去喂 DESeq2结果会出问题因为 DESeq2 的统计模型假设输入是整数计数。我个人的判断标准很简单打开文件看一眼数据如果里面是带小数点的连续值大概率是归一化后的表达量如果全是整数那基本就是 count。别嫌这一步土我见过太多人因为没确认这一点跑出来的差异基因列表全是假阳性或者全空。样本分组信息同样关键。你需要一个表标明每个样本属于哪一组。分组变量将来要进设计公式决定模型怎么算。分组命名尽量用英文、别有空格和特殊符号比如 control 和 treat别写成对照组和处理组否则某些函数处理起来会给你添麻烦。2.2 三大主流工具怎么选这个领域里最常被拿来讨论的就是 DESeq2、edgeR 和 limma 这三家。它们不是随便挑的背后对应的数据情况和统计假设不一样。工具适用输入核心思路什么时候优先用DESeq2原始 count负二项分布中位数比值归一化样本量少、有生物学重复、想省心edgeR原始 count负二项分布TMM 归一化样本量小、需要更灵活的自定义limma-voomcount 或连续值线性模型 voom 权重样本量大、设计复杂多因素我一般给新手的建议是只要你有正经的生物学重复直接用 DESeq2。它的封装程度高一行DESeq()就把归一化、离散估计、统计检验全做了出错概率低。edgeR 灵活但在参数上要你操心的地方多。limma 更适合样本量较大、需要处理复杂实验设计比如多时间点、多处理组合的场景。选工具的本质是选一个和你数据结构匹配的统计假设。生物学重复少比如每组 3 个数据波动大负二项模型更贴合测序计数的过离散特性所以 DESeq2 和 edgeR 更稳。样本量上去了数据接近正态线性模型那套就够用了。2.3 阈值设定的学问跑完统计你会得到每个基因的变化倍数和显著性。接下来要划两条线一条是变化倍数的门槛一条是显著性的门槛。变化倍数通常用 log2FoldChange 表示。为什么要取 log2因为表达上调 2 倍和下调到 1/2 在生物学上是对称的取 log2 后上调 2 倍是 1下调到一半是 -1绝对值相同画图和筛选都公平。常见阈值是 |log2FC| 1对应表达变化 2 倍以上。显著性门槛这边我强烈建议用校正后的 P 值padj / FDR而不是原始的 P 值。原因在于你一次要检验上万个基因假阳性会累积。原始 P 值小于 0.05 的基因里可能一大半是噪声。校正后 P 值控制的是错误发现率更靠谱。常用阈值是 padj 0.05严格一点的用 0.01。这两条线怎么定没有绝对标准。我的经验是先看数据整体分布再定阈值。如果差异基因少得可怜可以考虑把 log2FC 门槛降到 0.58即 1.5 倍但要在文章里说清楚理由。硬套一个 2 倍阈值可能把真正有生物学意义的小幅变化基因全滤掉了。3. 核心概念拆解与实操要点3.1 log2FoldChange 到底怎么算很多新手以为 log2FC 就是两组均值一比再取对数其实不完全对。在 DESeq2 里它用的是经过归一化和收缩估计后的值。归一化是为了消除测序深度差异——同样一个基因测序深的那次读数天然就多不校正的话会误判成上调。这里还有一个容易被忽视的机制log2FC 收缩shrinkage。当某个基因在对照组里几乎不表达处理组里表达了一点原始 log2FC 可能算出个很大的值但这种大变化往往是噪声。收缩估计会把这类不稳定的大倍数往中间拉让结果更稳。我的实操心得是对于下游要做富集分析或者筛选重点基因的场景用收缩后的 log2FC 更放心但如果只是想快速看个趋势原始值也够用。计算过程大致是先估计每个基因的离散度再拟合模型得到处理组相对对照组的系数这个系数就是 log2 尺度的变化量。之所以强调理解过程是因为当你看到某个基因 log2FC 大得离谱时你要能判断这是真信号还是低表达基因带来的假象而不是盲目相信数字。3.2 P 值与校正后 P 值的区别原始 P 值回答的是如果两组其实没有差异我观察到当前这么大差异的概率有多大。P 越小说明这个差异越不可能是偶然。问题在于你同时做了上万个这样的检验哪怕全是噪声按 5% 的显著性水平也会有一批基因碰巧显著。校正后 P 值就是来解决这个问题的。它把多重检验的因素考虑进去控制的是我判定为差异显著的基因里有多少比例其实是假的。所以你会看到同一个基因原始 P 值 0.001校正后可能变成 0.1直接从显著变不显著。这不是 bug而是必要的严谨。我的建议是报告结果时优先看 padj。如果审稿人或者导师想看原始 P 值你可以两个都给但结论一定要基于校正后的值下。3.3 数据标准化为什么绕不开标准化这件事说穿了就是让不同样本之间可比。测序深度不同、文库组成不同都会让计数产生系统性偏差。DESeq2 默认用的是中位数比值法edgeR 用的是 TMMlimma 有自己的归一化流程。它们的目标是一致的把技术因素造成的差异压下去把生物学差异留下来。有个实操细节值得提醒如果你做的是配对样本或者有时间序列结构标准化的处理需要相应调整不能套用最基础的流程。配对样本要在设计公式里加上配对因子否则组内个体差异会被误当成组间差异。这个坑我踩过当时差异基因列表里全是和个体相关的基因跟处理因素根本不沾边查了半天才发现是设计公式写错了。4. 一步步做差异分析从表达矩阵到结果表4.1 数据准备与环境搭建先把 R 环境和需要的包装好。DESeq2 通过 Bioconductor 安装别用普通 CRAN 的方式装装不上。# 安装 Bioconductor 管理器如果还没装 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) # 安装核心包 BiocManager::install(DESeq2) install.packages(ggplot2) install.packages(ggrepel)数据这块你需要两个文件一个是 count 矩阵一个是样本分组表。count 矩阵建议存成 CSV第一列是基因 ID后面每列一个样本。分组表两列一列样本名一列分组。我习惯在读入后立刻检查一下数据对不对countData - read.csv(counts.csv, row.names 1, check.names FALSE) colData - read.csv(metadata.csv, row.names 1, check.names FALSE) # 检查列名是否对得上 all(colnames(countData) rownames(colData)) # 检查分组情况 table(colData$condition)那个all(...)返回 TRUE说明样本名完全对齐这是后面不出错的前提。如果返回 FALSE说明两个表的样本顺序或者命名不一致得先手动调整。注意count 矩阵里不能有 NA也不建议有全是 0 的行。全 0 的基因没有任何信息会干扰离散度估计建议先过滤掉。4.2 DESeq2 完整流程下面是核心流程我把它写成一段可以直接跑的形式每步加了注释。library(DESeq2) # 构建 DESeq 数据集对象 dds - DESeqDataSetFromMatrix( countData countData, colData colData, design ~ condition ) # 过滤低表达基因至少在一半样本中有一定表达 keep - rowSums(counts(dds) 10) ceiling(ncol(dds) / 2) dds - dds[keep, ] # 设定参考水平对照组作为基准 dds$condition - relevel(dds$condition, ref control) # 跑差异分析 dds - DESeq(dds) # 提取结果 res - results(dds, contrast c(condition, treat, control)) resdesign ~ condition这一句是告诉模型我要比较的是 condition 这个变量。relevel那步是设定谁当基准这样算出来的 log2FC 方向才符合你的预期——处理组相对对照组上调就是正值。DESeq()这一步是最耗时的它会完成归一化、离散度估计和统计检验。样本不多的话几分钟就出来了。跑完之后results()提取表格里面有几列baseMean、log2FoldChange、lfcSE、stat、pvalue、padj。如果想让结果更稳可以加一步收缩library(apeglm) res_shrunk - lfcShrink(dds, coef condition_treat_vs_control, type apeglm)4.3 结果提取与初步筛选拿到结果表之后先整理成能用的格式再做筛选。res_df - as.data.frame(res) res_df - res_df[order(res_df$padj), ] res_df$gene - rownames(res_df) # 按阈值筛选差异基因 res_df$group - not_sig res_df$group[res_df$log2FoldChange 1 res_df$padj 0.05] - up res_df$group[res_df$log2FoldChange -1 res_df$padj 0.05] - down # 看看各类有多少 table(res_df$group) # 导出 write.csv(res_df, deg_results.csv, row.names FALSE)跑完这一步你就能拿到一个带分组标签的结果表。我通常先看table(res_df$group)的结果上调和下调各有多少个。如果上调下调数量严重失衡或者总数少得异常就得回头查问题别急着往下走。4.4 实操记录与参数说明这里记几个我实际跑的时候的观察。第一过滤低表达基因那一步很重要不做的话一堆全 0 或者只有一两个 read 的基因会拉低离散度估计的稳定性结果图会很难看。第二relevel不可省略否则 log2FC 的正负方向可能和你认知相反导致火山图上下调标反。第三padj 那一列如果有 NA通常是该基因在独立过滤中被剔除了这是正常现象不用慌。关于过滤阈值counts 10这是个经验值。定得太严一些低表达但有功能的基因会被丢掉定得太松噪声基因一堆。我的做法是结合实验本身的深度来调如果整体测序深度高可以适当把阈值调到 20 甚至更高。5. 一分钟画出火山图ggplot2 实操5.1 数据整理成火山图需要的格式火山图的核心是三个映射横轴 log2FoldChange纵轴 -log10(padj)颜色分组。整理一下数据plot_df - res_df[!is.na(res_df$padj), ] # 纵轴用校正后 P 值取负对数 plot_df$logP - -log10(plot_df$padj) # 避免极端值把图压扁可以设个上限 plot_df$logP[plot_df$logP 300] - 300为什么要取 -log10因为原来是 P 越小越显著取负对数之后P 越小值越大点就越高。这样越显著越靠上符合直觉。设上限这一步是为了防止个别极小 P 值把整个纵轴拉得太长其他点全挤在底部看不清。5.2 基础火山图代码library(ggplot2) library(ggrepel) ggplot(plot_df, aes(x log2FoldChange, y logP, color group)) geom_point(alpha 0.6, size 1.2) scale_color_manual(values c(up #E64B35, down #4DBBD5, not_sig grey70)) geom_vline(xintercept c(-1, 1), linetype dashed, color grey40) geom_hline(yintercept -log10(0.05), linetype dashed, color grey40) labs(x log2 Fold Change, y -log10(FDR), color Group) theme_bw()这段代码跑出来就是标准的火山图。两条竖虚线是 log2FC 的门槛横虚线是显著性门槛四个区域里右上和左上就是显著上调和下调的基因。颜色上我习惯上调用暖色、下调用冷色灰色留给不显著的视觉上一眼能分清。5.3 配色、标注与美化基础图出来之后为了能放进报告或者论文通常还要标注几个重点基因。# 挑出最显著的若干个基因做标注 top_genes - head(plot_df[order(plot_df$padj), ], 10) ggplot(plot_df, aes(x log2FoldChange, y logP, color group)) geom_point(alpha 0.6, size 1.2) scale_color_manual(values c(up #E64B35, down #4DBBD5, not_sig grey70)) geom_vline(xintercept c(-1, 1), linetype dashed, color grey40) geom_hline(yintercept -log10(0.05), linetype dashed, color grey40) geom_text_repel(data top_genes, aes(label gene), size 3, max.overlaps 20) labs(x log2 Fold Change, y -log10(FDR), color Group) theme_bw(base_size 13) theme(panel.grid.minor element_blank())geom_text_repel会自动帮你调整标签位置避免文字和点重叠比手动调坐标省心太多。max.overlaps控制最多允许重叠的数量调大一点可以标更多基因。配色上给两个方向参考一种是像上面这样的红蓝对比另一种是蓝黄对比。别用太刺眼的纯色也别搞一堆渐变色火山图上点的颜色是分类用的不是连续量用离散配色就对了。还有个小细节alpha控制点透明度点多了会糊成一片调低透明度能看出一层层的密度。6. 常见问题与排查技巧实录6.1 报错与异常排查跑差异分析的过程中报错是家常便饭但大部分集中在几个地方。DESeqDataSetFromMatrix报错八成是 countData 和 colData 的样本名对不上或者 countData 里有非整数。这时候回头看all(colnames(...) rownames(...))是不是 TRUE。DESeq()报错可能是设计公式里的变量名和 colData 里的列名不一致检查一下colData的列名到底是什么。还有一个隐蔽的问题分组变量被 R 当成数值处理了。如果你写的是condition 1和condition 2R 可能把它当连续变量那你算出来的就不是两组比较而是回归。解决办法是把分组转成因子或者直接用字母命名。6.2 结果不合理的排查思路结果出来之后先别急着高兴或者沮丧几个关键点要查。差异基因数量是不是合理如果两组差异明显通常会有几百到几千个差异基因如果只有个位数可能阈值太严或者分组有问题如果几乎所有基因都显著那更可能是数据或者归一化出了问题。火山图的形状也能提供信息。正常情况是一个漏斗形或者对称的 V 形中间密集、两边稀疏。如果图形很奇怪比如一边倒、或者点分布完全不对称首先怀疑 log2FC 的方向是不是搞反了其次看是不是某一组样本质量特别差。我踩过的一个坑是忘了过滤低表达基因结果图上出现一条横在顶部的点带全是那些在对照组为 0、处理组有一点表达的基因log2FC 大得吓人。过滤之后就正常了。所以过滤这一步不是可选项是必须项。6.3 常见问题速查表现象可能原因解决方向样本名对不上报错两个表顺序或命名不一致检查并统一样本名差异基因特别少阈值太严 / 分组错误放宽阈值 / 核对分组几乎所有基因显著归一化异常 / 批次效应检查标准化 / 加批次因子火山图一边倒log2FC 方向反了检查 relevel 设置顶部一横排异常点低表达基因未过滤增加过滤步骤padj 大量为 NA独立过滤剔除正常现象可降低过滤强度这张表我建议存下来遇到问题先对照着查能省下大量瞎试的时间。差异分析这东西代码本身不难难的是数据准备和结果判读。把这两头把握住中间那段流程跑起来其实很快。最后分享一个我自己的习惯每次跑完差异分析除了火山图我都会顺手存一份结果表并记下这次用的参数阈值、过滤条件、工具版本。因为过段时间回头看或者换一批数据复现没有这些记录就得从头猜。参数记清楚了下一次复现就是几分钟的事这才是真正的一秒学会背后该有的底气。 
长河君 · 2026/9/17 8:14:03 · 阅读 1.2k RNA-seq差异基因分析与火山图绘制：DESeq2实战
1. 差异基因分析到底在解决什么问题做转录组的人绕来绕去最后都要落到同一个问题上哪些基因在实验组和对照组之间变了。差异基因分析就是回答这个问题的标准动作而火山图是把结果一次性摊开给你看的图。标题里说一秒学会我理解这种说法背后是很多新手被一堆报错和参数劝退的无奈——真正上手过的人都知道画图那一步确实快难的是前面把数据和对齐逻辑准备利索。我这次就按先把原理讲透再给能直接抄的代码这个路数来写。适合谁看如果你手头有一张表达矩阵counts 表或者 TPM/FPKM 表有分组信息想拿到差异基因列表并且想要一张能放进论文或者汇报 PPT 的火山图那这篇基本能覆盖你从零到出图的全流程。哪怕你之前没怎么用过 R只要照着步骤走也能跑通。先说清楚这两件事的定位。差异基因分析本质是一次带统计约束的对比它不是简单地算个比值就完事而是要考虑基因表达本身的波动、测序深度差异、生物学重复的离散程度。火山图则是把每个基因的两个核心指标——变化倍数和显著性——同时画在二维平面上。横轴放变化幅度纵轴放显著性越往两边越靠上的点就是越值得你关注的基因。理解了这一点后面所有的参数选择、代码调整都有了落脚点不会变成盲人摸象。2. 分析前的整体设计与关键参数选型2.1 从一张表达矩阵说起一切从表达矩阵开始。行是基因列是样本格子里的数字是该基因在该样本里的表达量。这个矩阵通常来自上游的比对和定量流程比如用比对工具得到每个基因的 read 计数。这里有个新手最容易忽略的点你必须明确手里的数值到底是原始 count 还是已经归一化过的表达量。这决定了你后面用哪套差异分析工具。如果拿 FPKM 去喂 DESeq2结果会出问题因为 DESeq2 的统计模型假设输入是整数计数。我个人的判断标准很简单打开文件看一眼数据如果里面是带小数点的连续值大概率是归一化后的表达量如果全是整数那基本就是 count。别嫌这一步土我见过太多人因为没确认这一点跑出来的差异基因列表全是假阳性或者全空。样本分组信息同样关键。你需要一个表标明每个样本属于哪一组。分组变量将来要进设计公式决定模型怎么算。分组命名尽量用英文、别有空格和特殊符号比如 control 和 treat别写成对照组和处理组否则某些函数处理起来会给你添麻烦。2.2 三大主流工具怎么选这个领域里最常被拿来讨论的就是 DESeq2、edgeR 和 limma 这三家。它们不是随便挑的背后对应的数据情况和统计假设不一样。工具适用输入核心思路什么时候优先用DESeq2原始 count负二项分布中位数比值归一化样本量少、有生物学重复、想省心edgeR原始 count负二项分布TMM 归一化样本量小、需要更灵活的自定义limma-voomcount 或连续值线性模型 voom 权重样本量大、设计复杂多因素我一般给新手的建议是只要你有正经的生物学重复直接用 DESeq2。它的封装程度高一行DESeq()就把归一化、离散估计、统计检验全做了出错概率低。edgeR 灵活但在参数上要你操心的地方多。limma 更适合样本量较大、需要处理复杂实验设计比如多时间点、多处理组合的场景。选工具的本质是选一个和你数据结构匹配的统计假设。生物学重复少比如每组 3 个数据波动大负二项模型更贴合测序计数的过离散特性所以 DESeq2 和 edgeR 更稳。样本量上去了数据接近正态线性模型那套就够用了。2.3 阈值设定的学问跑完统计你会得到每个基因的变化倍数和显著性。接下来要划两条线一条是变化倍数的门槛一条是显著性的门槛。变化倍数通常用 log2FoldChange 表示。为什么要取 log2因为表达上调 2 倍和下调到 1/2 在生物学上是对称的取 log2 后上调 2 倍是 1下调到一半是 -1绝对值相同画图和筛选都公平。常见阈值是 |log2FC| 1对应表达变化 2 倍以上。显著性门槛这边我强烈建议用校正后的 P 值padj / FDR而不是原始的 P 值。原因在于你一次要检验上万个基因假阳性会累积。原始 P 值小于 0.05 的基因里可能一大半是噪声。校正后 P 值控制的是错误发现率更靠谱。常用阈值是 padj 0.05严格一点的用 0.01。这两条线怎么定没有绝对标准。我的经验是先看数据整体分布再定阈值。如果差异基因少得可怜可以考虑把 log2FC 门槛降到 0.58即 1.5 倍但要在文章里说清楚理由。硬套一个 2 倍阈值可能把真正有生物学意义的小幅变化基因全滤掉了。3. 核心概念拆解与实操要点3.1 log2FoldChange 到底怎么算很多新手以为 log2FC 就是两组均值一比再取对数其实不完全对。在 DESeq2 里它用的是经过归一化和收缩估计后的值。归一化是为了消除测序深度差异——同样一个基因测序深的那次读数天然就多不校正的话会误判成上调。这里还有一个容易被忽视的机制log2FC 收缩shrinkage。当某个基因在对照组里几乎不表达处理组里表达了一点原始 log2FC 可能算出个很大的值但这种大变化往往是噪声。收缩估计会把这类不稳定的大倍数往中间拉让结果更稳。我的实操心得是对于下游要做富集分析或者筛选重点基因的场景用收缩后的 log2FC 更放心但如果只是想快速看个趋势原始值也够用。计算过程大致是先估计每个基因的离散度再拟合模型得到处理组相对对照组的系数这个系数就是 log2 尺度的变化量。之所以强调理解过程是因为当你看到某个基因 log2FC 大得离谱时你要能判断这是真信号还是低表达基因带来的假象而不是盲目相信数字。3.2 P 值与校正后 P 值的区别原始 P 值回答的是如果两组其实没有差异我观察到当前这么大差异的概率有多大。P 越小说明这个差异越不可能是偶然。问题在于你同时做了上万个这样的检验哪怕全是噪声按 5% 的显著性水平也会有一批基因碰巧显著。校正后 P 值就是来解决这个问题的。它把多重检验的因素考虑进去控制的是我判定为差异显著的基因里有多少比例其实是假的。所以你会看到同一个基因原始 P 值 0.001校正后可能变成 0.1直接从显著变不显著。这不是 bug而是必要的严谨。我的建议是报告结果时优先看 padj。如果审稿人或者导师想看原始 P 值你可以两个都给但结论一定要基于校正后的值下。3.3 数据标准化为什么绕不开标准化这件事说穿了就是让不同样本之间可比。测序深度不同、文库组成不同都会让计数产生系统性偏差。DESeq2 默认用的是中位数比值法edgeR 用的是 TMMlimma 有自己的归一化流程。它们的目标是一致的把技术因素造成的差异压下去把生物学差异留下来。有个实操细节值得提醒如果你做的是配对样本或者有时间序列结构标准化的处理需要相应调整不能套用最基础的流程。配对样本要在设计公式里加上配对因子否则组内个体差异会被误当成组间差异。这个坑我踩过当时差异基因列表里全是和个体相关的基因跟处理因素根本不沾边查了半天才发现是设计公式写错了。4. 一步步做差异分析从表达矩阵到结果表4.1 数据准备与环境搭建先把 R 环境和需要的包装好。DESeq2 通过 Bioconductor 安装别用普通 CRAN 的方式装装不上。# 安装 Bioconductor 管理器如果还没装 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) # 安装核心包 BiocManager::install(DESeq2) install.packages(ggplot2) install.packages(ggrepel)数据这块你需要两个文件一个是 count 矩阵一个是样本分组表。count 矩阵建议存成 CSV第一列是基因 ID后面每列一个样本。分组表两列一列样本名一列分组。我习惯在读入后立刻检查一下数据对不对countData - read.csv(counts.csv, row.names 1, check.names FALSE) colData - read.csv(metadata.csv, row.names 1, check.names FALSE) # 检查列名是否对得上 all(colnames(countData) rownames(colData)) # 检查分组情况 table(colData$condition)那个all(...)返回 TRUE说明样本名完全对齐这是后面不出错的前提。如果返回 FALSE说明两个表的样本顺序或者命名不一致得先手动调整。注意count 矩阵里不能有 NA也不建议有全是 0 的行。全 0 的基因没有任何信息会干扰离散度估计建议先过滤掉。4.2 DESeq2 完整流程下面是核心流程我把它写成一段可以直接跑的形式每步加了注释。library(DESeq2) # 构建 DESeq 数据集对象 dds - DESeqDataSetFromMatrix( countData countData, colData colData, design ~ condition ) # 过滤低表达基因至少在一半样本中有一定表达 keep - rowSums(counts(dds) 10) ceiling(ncol(dds) / 2) dds - dds[keep, ] # 设定参考水平对照组作为基准 dds$condition - relevel(dds$condition, ref control) # 跑差异分析 dds - DESeq(dds) # 提取结果 res - results(dds, contrast c(condition, treat, control)) resdesign ~ condition这一句是告诉模型我要比较的是 condition 这个变量。relevel那步是设定谁当基准这样算出来的 log2FC 方向才符合你的预期——处理组相对对照组上调就是正值。DESeq()这一步是最耗时的它会完成归一化、离散度估计和统计检验。样本不多的话几分钟就出来了。跑完之后results()提取表格里面有几列baseMean、log2FoldChange、lfcSE、stat、pvalue、padj。如果想让结果更稳可以加一步收缩library(apeglm) res_shrunk - lfcShrink(dds, coef condition_treat_vs_control, type apeglm)4.3 结果提取与初步筛选拿到结果表之后先整理成能用的格式再做筛选。res_df - as.data.frame(res) res_df - res_df[order(res_df$padj), ] res_df$gene - rownames(res_df) # 按阈值筛选差异基因 res_df$group - not_sig res_df$group[res_df$log2FoldChange 1 res_df$padj 0.05] - up res_df$group[res_df$log2FoldChange -1 res_df$padj 0.05] - down # 看看各类有多少 table(res_df$group) # 导出 write.csv(res_df, deg_results.csv, row.names FALSE)跑完这一步你就能拿到一个带分组标签的结果表。我通常先看table(res_df$group)的结果上调和下调各有多少个。如果上调下调数量严重失衡或者总数少得异常就得回头查问题别急着往下走。4.4 实操记录与参数说明这里记几个我实际跑的时候的观察。第一过滤低表达基因那一步很重要不做的话一堆全 0 或者只有一两个 read 的基因会拉低离散度估计的稳定性结果图会很难看。第二relevel不可省略否则 log2FC 的正负方向可能和你认知相反导致火山图上下调标反。第三padj 那一列如果有 NA通常是该基因在独立过滤中被剔除了这是正常现象不用慌。关于过滤阈值counts 10这是个经验值。定得太严一些低表达但有功能的基因会被丢掉定得太松噪声基因一堆。我的做法是结合实验本身的深度来调如果整体测序深度高可以适当把阈值调到 20 甚至更高。5. 一分钟画出火山图ggplot2 实操5.1 数据整理成火山图需要的格式火山图的核心是三个映射横轴 log2FoldChange纵轴 -log10(padj)颜色分组。整理一下数据plot_df - res_df[!is.na(res_df$padj), ] # 纵轴用校正后 P 值取负对数 plot_df$logP - -log10(plot_df$padj) # 避免极端值把图压扁可以设个上限 plot_df$logP[plot_df$logP 300] - 300为什么要取 -log10因为原来是 P 越小越显著取负对数之后P 越小值越大点就越高。这样越显著越靠上符合直觉。设上限这一步是为了防止个别极小 P 值把整个纵轴拉得太长其他点全挤在底部看不清。5.2 基础火山图代码library(ggplot2) library(ggrepel) ggplot(plot_df, aes(x log2FoldChange, y logP, color group)) geom_point(alpha 0.6, size 1.2) scale_color_manual(values c(up #E64B35, down #4DBBD5, not_sig grey70)) geom_vline(xintercept c(-1, 1), linetype dashed, color grey40) geom_hline(yintercept -log10(0.05), linetype dashed, color grey40) labs(x log2 Fold Change, y -log10(FDR), color Group) theme_bw()这段代码跑出来就是标准的火山图。两条竖虚线是 log2FC 的门槛横虚线是显著性门槛四个区域里右上和左上就是显著上调和下调的基因。颜色上我习惯上调用暖色、下调用冷色灰色留给不显著的视觉上一眼能分清。5.3 配色、标注与美化基础图出来之后为了能放进报告或者论文通常还要标注几个重点基因。# 挑出最显著的若干个基因做标注 top_genes - head(plot_df[order(plot_df$padj), ], 10) ggplot(plot_df, aes(x log2FoldChange, y logP, color group)) geom_point(alpha 0.6, size 1.2) scale_color_manual(values c(up #E64B35, down #4DBBD5, not_sig grey70)) geom_vline(xintercept c(-1, 1), linetype dashed, color grey40) geom_hline(yintercept -log10(0.05), linetype dashed, color grey40) geom_text_repel(data top_genes, aes(label gene), size 3, max.overlaps 20) labs(x log2 Fold Change, y -log10(FDR), color Group) theme_bw(base_size 13) theme(panel.grid.minor element_blank())geom_text_repel会自动帮你调整标签位置避免文字和点重叠比手动调坐标省心太多。max.overlaps控制最多允许重叠的数量调大一点可以标更多基因。配色上给两个方向参考一种是像上面这样的红蓝对比另一种是蓝黄对比。别用太刺眼的纯色也别搞一堆渐变色火山图上点的颜色是分类用的不是连续量用离散配色就对了。还有个小细节alpha控制点透明度点多了会糊成一片调低透明度能看出一层层的密度。6. 常见问题与排查技巧实录6.1 报错与异常排查跑差异分析的过程中报错是家常便饭但大部分集中在几个地方。DESeqDataSetFromMatrix报错八成是 countData 和 colData 的样本名对不上或者 countData 里有非整数。这时候回头看all(colnames(...) rownames(...))是不是 TRUE。DESeq()报错可能是设计公式里的变量名和 colData 里的列名不一致检查一下colData的列名到底是什么。还有一个隐蔽的问题分组变量被 R 当成数值处理了。如果你写的是condition 1和condition 2R 可能把它当连续变量那你算出来的就不是两组比较而是回归。解决办法是把分组转成因子或者直接用字母命名。6.2 结果不合理的排查思路结果出来之后先别急着高兴或者沮丧几个关键点要查。差异基因数量是不是合理如果两组差异明显通常会有几百到几千个差异基因如果只有个位数可能阈值太严或者分组有问题如果几乎所有基因都显著那更可能是数据或者归一化出了问题。火山图的形状也能提供信息。正常情况是一个漏斗形或者对称的 V 形中间密集、两边稀疏。如果图形很奇怪比如一边倒、或者点分布完全不对称首先怀疑 log2FC 的方向是不是搞反了其次看是不是某一组样本质量特别差。我踩过的一个坑是忘了过滤低表达基因结果图上出现一条横在顶部的点带全是那些在对照组为 0、处理组有一点表达的基因log2FC 大得吓人。过滤之后就正常了。所以过滤这一步不是可选项是必须项。6.3 常见问题速查表现象可能原因解决方向样本名对不上报错两个表顺序或命名不一致检查并统一样本名差异基因特别少阈值太严 / 分组错误放宽阈值 / 核对分组几乎所有基因显著归一化异常 / 批次效应检查标准化 / 加批次因子火山图一边倒log2FC 方向反了检查 relevel 设置顶部一横排异常点低表达基因未过滤增加过滤步骤padj 大量为 NA独立过滤剔除正常现象可降低过滤强度这张表我建议存下来遇到问题先对照着查能省下大量瞎试的时间。差异分析这东西代码本身不难难的是数据准备和结果判读。把这两头把握住中间那段流程跑起来其实很快。最后分享一个我自己的习惯每次跑完差异分析除了火山图我都会顺手存一份结果表并记下这次用的参数阈值、过滤条件、工具版本。因为过段时间回头看或者换一批数据复现没有这些记录就得从头猜。参数记清楚了下一次复现就是几分钟的事这才是真正的一秒学会背后该有的底气。