edgeR/limma 必做的 5 组验证图和 2 个判断原则引言跑完差异表达分析拿到一张火山图看到红红绿绿的差异基因第一反应通常是不错有结果了。但能出火山图 ≠ 结果可靠。我在审稿和帮同学看数据时反复遇到一个场景火山图漂漂亮亮但翻到后面的样本聚类图发现对照组和处理组的样本完全混在一起。这意味着差异基因可能不是处理导致的而是其他混杂因素比如建库批次、样本采集时间、RNA 质量差异造成的。这篇文章不讲怎么做差异分析而是讲拿到结果之后必须做的质量控制。我会给出 5 组验证图 2 个判断原则帮你快速判断你的差异基因到底能不能信一、先问自己一个问题你的火山图里那些显著差异基因是真实的生物学差异还是技术噪音的放大下面这个表列出了常见的假阳性陷阱问题典型表现后果低表达基因 P 值劫持counts 10 的基因大量出现在显著列表功能注释全是错的这些基因根本没被检测到批次效应混杂处理组和对照组分别在不同批次建库你以为的差异其实是批次差异离群样本某个样本 reads 数是其他的 2 倍整个分析被一个异常样本拉偏样本分组不成立聚类图/PCA 图里两组完全重叠没有分组差异 所有差异基因都是噪音P 值分布异常QQ 图严重偏离对角线统计模型假设被违反P 值不可信个人经验这五个问题我在真实项目中每个都踩过至少一次。最惨的一次是帮师弟调代码调了两周最后发现差异基因全是批次效应。先做 QC再相信结果——这句话写在每个生信流程的教科书里但真正做到的人不多。二、验证图 1MA 图 —— 看低表达区域的假阳性MA 图Mean-Abundance Plot是差异分析最基础的诊断图。x 轴是基因的平均表达量y 轴是 log2 Fold Change。它的核心功能检查差异是否集中在低表达基因。怎么看这张图正常情况差异基因均匀分布在各个表达水平红蓝点散布在全图。有问题差异基因集中在左侧低表达区域橙色圈标注的可疑假阳性大量出现。这些基因 counts 极低离散度估计不稳定P 值不可信。个人经验edgeR 默认不会过滤低表达基因。我一般先做一步过滤——只保留 CPM 1 且至少在一半样本中检测到的基因。这一步能砍掉 30%-50% 的低表达假阳性。三、验证图 2样本聚类热图 —— 看分组是否成立如果你的实验设计是 A vs B那在样本聚类热图里A 组的样本应该聚在一起B 组的样本应该聚在一起。如果两组样本完全混在一起说明组间差异远小于个体差异差异分析的结果基本不可信。怎么看这张图正常情况列方向样本方向的聚类树明显分成两支Control 和 Treatment 各占一支。有问题两组样本交错排列没有分组趋势。说明实验处理的效果太弱或者混杂因素太强。个人经验样本量太少每组 3时聚类热图的结论要谨慎。如果经费允许每组至少 4-5 个重复。四、验证图 3QQ 图 —— 看 P 值分布是否合理QQ 图把每个基因的 P 值跟均匀分布的期望值做比较——如果你的模型完全正确、没有任何差异信号所有点应该落在对角线上。怎么看这张图正常情况有真实差异信号大部分点在对角线上尾部右上角明显偏离对角线向上翘。翘得越厉害差异信号越强。有问题P 值膨胀从早期左下角就偏离对角线意味着大量基因的 P 值比期望值小——可能是模型假设不对、有未控制的混杂因素。更严重P 值收缩所有点都在对角线下方通常是样本量太小或过度离散造成的。个人经验QQ 图的早期偏离是一个经常被忽略的信号。如果左下角就偏离了对角线先别急着看差异基因列表回去检查是不是有未校正的批次效应。五、验证图 4样本 QC —— 看有没有离群样本差异分析对离群样本极其敏感。一个 reads 数异常高的样本能把整个分析的 P 值分布搅乱。怎么看这两张图左图箱线图所有样本总 reads 应在合理范围内变异系数 20%。某个样本异常偏高需要检查建库浓度或是否有污染。右图密度图所有样本表达分布应大致重叠。某个样本分布明显偏移说明归一化可能无法完全校正。个人经验我见过一个真实案例——某课题组 12 个 RNA-seq 样本其中一个总 reads 是其他的 3 倍后来发现是测序中心跑了两条 lane 合并了。不做这一步 QC整个分析被拉偏。六、验证图 5PCA 图 —— 看批次效应PCA 图是检查批次效应最直观的工具。同一个实验分组的样本应该靠近如果按建库批次、测序日期分开你看到的差异可能不是生物学差异。怎么看这张图正常情况两组的椭圆明显分开说明处理效应大于个体差异。有问题两组椭圆完全重叠说明处理效应太弱。如果上面看起来不理想换个角度着色——按批次着色个人经验批次效应是差异分析里最常见的陷阱没有之一。我的建议是拿到数据先画 PCA如果发现批次效应立刻在 design formula 里加 Batch 因子比如design ~ Batch Group。不要先跑完差异分析再回头修。七、判断原则 1生物学重复的一致性统计检验的效力取决于重复数和组内变异。每组至少 3 个重复这是底线。2 个重复的差异分析只比较了 2 个样本。同一组内样本的皮尔逊相关系数应该 0.9。如果组内有样本相关系数 0.8优先排查。样本量越大检测到的差异基因越可信但不意味着越多越好。个人经验我常用的快速判断标准——组内中位数相关 0.95组间中位数相关 0.9。满足这两个条件分组合理差异分析结果基本可信。八、判断原则 2P 值与效应量双重阈值只看 P 值是差异分析里最普遍的错误。P 值只告诉你有统计差异不告诉你差异有多大。样本量越大越小的 fold change 也能显著。推荐标准阈值参数推荐值说明P 值 / FDRFDR 0.05BH 校正后的 FDR 比 raw P 值更严格log2 Fold Change|log2FC| 12 倍差异可根据领域调整到 1.5 或 0.58表达量过滤CPM 1至少一半样本中排除低表达噪音个人经验很多审稿人专门挑 log2FC 在 0.5-1 之间的弱差异基因质疑结论。如果论文核心论点依赖这些基因最好准备 qPCR 或 Western Blot 验证。九、完整代码5 图 1 次出下面是从原始 counts 矩阵到 5 张验证图的一体化 R 脚本可直接复制替换自己的数据十、写在最后这篇文章的核心观点只有一句话先验证再相信。差异分析工具DESeq2、edgeR、limma不会告诉你结果对不对——它们只负责做数学运算。质量控制是你的责任不是算法的责任。我的建议流程是过滤低表达基因filterByExpr或 CPM 1画 5 张 QC 图MA 图、热图、QQ 图、样本 QC、PCA检查两个判断标准组内相关 0.95、FDR logFC 双阈值全部通过 → 相信结果继续下游分析有异常 → 排查问题修完再重跑这 5 张图 2 个标准大概 10 分钟就能跑完。比起发现结果不可靠后返工两周这 10 分钟是性价比最高的投入。如果你正在做差异分析建议把这段 QC 代码加到你的流程开头。有问题欢迎在公众号后台留言交流。