📚 课程目录第4章差异分析结果的解读

4 差异分析结果的解读归纳总结

差异分析跑完了——拿到一长串基因，p值、padj、log2FC，眼花缭乱。然后呢？

说实话，很多人卡在这一步。我见过的项目里，至少有三分之一的分析报告，在这里翻了车。要么把噪声当成信号，要么错过了真正的生物学故事。

说白了，你得知道哪些数字是“真家伙”，哪些是“偶然的幸运”。今天我就把解读差异结果的核心思路掰开揉碎，结合我踩过的坑，给你一套可以落地的方法。

4.1 差异基因的筛选标准——别死搬教条

最常用的筛选条件是padj ＜ 0.05 且 |log2FC| ＞ 1。但你想想，这个阈值真的适合所有数据吗？

阈值组合	适用场景	我踩过的坑
padj＜0.05 + \|log2FC\|＞1	常规RNA-seq，样本数≥3	遇到高变异性数据，很多真基因被漏掉
padj＜0.1 + \|log2FC\|＞0.58	单细胞RNA-seq或样本量少的实验	曾经用了标准阈值，结果什么都没找到……
padj＜0.05 + \|log2FC\|＞2	需要高度保守的候选基因列表	适合做后续实验验证，但会错过很多中度变化

我建议你：先跑一个默认阈值看看火山图，如果上下调基因太少或太多，再调整。我个人习惯会同时保留两套阈值的结果——一套严格用于后续富集，一套宽松用于探索性分析。

⚠ 千万不要只看p值不看padj！
原始p值在多重检验下毫无意义。我见过有人拿p＜0.01直接说是差异基因——结果一半是假阳性。

4.2 火山图：一眼看穿真假信号

火山图就是差异分析的“第一印象”。横轴是log2FC，纵轴是 -log10(padj)。左上方和右上方那些点，就是我们要找的差异基因。

这里有几个经验：

看对称性：如果上下调数量差不多，说明结果可靠。如果全部偏向一侧，要么是标准化出了问题，要么是实验设计有偏向。
留意纵轴极端值：有些基因padj小到10⁻³⁰，但log2FC只有0.5。这种基因虽然统计显著，但可能是由于极低的变异造成的，生物学意义有限。
注意横轴边界：如果log2FC的绝对值超过10或更多，多半是表达量为0或极低导致的假象。我曾在项目中看到log2FC=20的基因，一查原来是某个样本完全没表达——这事得排除。

💡 避坑小技巧
绘制火山图时，我总是用 geom_text_repel 给top10标志性基因加上标签。这样一眼就能看到是不是已知的marker基因被富集出来，快速验证分析的合理性。

# R 代码：绘制基础火山图（使用DESeq2结果的示例）
library(ggplot2)
library(ggrepel)

# 假设 res 是 DESeq2 结果，已转为数据框
res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df)
res_df$sig <- "no"
res_df$sig[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1] <- "up"
res_df$sig[res_df$padj < 0.05 & res_df$log2FoldChange < -1] <- "down"

ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1.5) +
  scale_color_manual(values = c("down" = "#3498db", "no" = "grey60", "up" = "#e74c3c")) +
  theme_minimal() +
  labs(title = "Volcano Plot", x = "log2 Fold Change", y = "-log10(adjusted p-value)")

4.3 MA图：揪出表达量依赖的偏差

MA图的横轴是基因的平均表达量（log2 counts per million），纵轴是log2FC。正常情况下的点应该对称分布在y=0两侧，且无明显的喇叭口形状。

为什么要看MA图？因为差异分析的结果常常受到低表达基因的影响。低表达基因的log2FC波动大，容易产生假阳。如果MA图上看到低表达区域（靠左侧）的log2FC变异非常大，就需要考虑过滤掉低表达基因再分析。

  📊 我常用的评估标准

  一个好的MA图：

  红色（显著上调）和蓝色（显著下调）基本对称
点分布呈水平的带状，没有明显的上下弯曲
低表达区域的点不能比高表达区域明显更多或更分散

4.4 热图与聚类：把抽象数字变成视觉模式

筛选出差异基因后，最直观的就是画热图。行是基因，列是样本，颜色表示表达量高低。

我个人习惯用 pheatmap 包，配色用蓝-白-红（低-中-高）。注意几点：

一定要做行聚类和列聚类，看看样本是否按实验分组自动分开。如果对照组和治疗组混在一起，要么差异基因选错了，要么样本有问题。
不要全部基因都画进去。我会先按padj排序取前50个差异最大的基因，或者用方差筛选前500个高变基因。画太多基因画面太拥挤，看不出模式。
归一化方式：建议每行做Z-score，这样上下调一目了然。我踩过的坑是直接用了raw counts，结果上下调被表达量高的基因主导。

# R 代码：绘制热图（使用pheatmap）
library(pheatmap)

# 假设 dds 是 DESeq2 对象，已进行 varianceStabilizingTransformation
vsd <- vst(dds, blind = FALSE)
top_genes <- head(order(res$padj), 50)
mat <- assay(vsd)[top_genes, ]
# 按行做Z-score
mat_z <- t(scale(t(mat)))

# 注释信息
annotation_col <- data.frame(
  condition = colData(dds)$condition,
  row.names = colnames(dds)
)

pheatmap(mat_z,
         annotation_col = annotation_col,
         color = colorRampPalette(c("#3498db", "white", "#e74c3c"))(100),
         main = "Top 50 DEGs Heatmap",
         scale = "none",  # 已手动缩放
         show_rownames = TRUE,
         show_colnames = TRUE)

4.5 功能富集分析：给基因穿上生物学意义的外衣

差异基因列表本身只是名字。要理解它们参与了什么通路、影响了哪些生物学过程，必须做功能富集分析。

常用工具：clusterProfiler（R包）、DAVID、WebGestalt。KEGG和GO是最常见的数据库。

这里有一个我反复强调的要点：富集结果要警惕“假阳性”。如果你丢进去500个差异基因，得到几十个显著富集的通路，这很正常。但如果每个通路都显著（p值全在0.001以下），反而有问题——可能是你的基因列表太长，或者背景数据库太小。我曾在项目中遇到这种情况，最后发现是没过滤掉非编码RNA导致的……

💡 我的筛选习惯

先用 padj ＜ 0.05 做严格筛选
富集时使用 qvalue ＜ 0.05 作为显著性阈值
只看 包含6～200个基因 的通路，太小或太大的通路往往不靠谱
手动检查1-2个通路的表达模式，确保不是由少数基因驱动

# R 代码：使用 clusterProfiler 做GO富集
library(clusterProfiler)
library(org.Hs.eg.db)

# 将基因ID转换为ENTREZID
gene_list <- rownames(res_df[res_df$sig %in% c("up", "down"), ])
entrez <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

ego <- enrichGO(gene          = entrez$ENTREZID,
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.05)

head(ego)

4.6 常见陷阱与避坑指南

做了这么多年差异分析，我总结出下面几个最容易出问题的地方——希望你一个也别踩到。

误区一：只关注padj最小的基因
很多时候padj最小的基因其实log2FC很小。我建议你先看log2FC绝对值超过1的基因，再在其中挑padj小的。
误区二：不同软件的结果直接对比
DESeq2、edgeR、limma用的模型不同，log2FC数值会有差异。不要直接对比两个软件的log2FC值，要对比方向和显著性。
误区三：忽略批次效应
如果热图显示样本按批次聚类而非按分组聚类，你的差异基因里混入了大量批次相关基因。需要用ComBat等工具校正后再做差异分析。
误区四：过度依赖自动化注释
富集结果里的通路名字有时候会误导人。我记得有一次富集出来“嗅觉转导通路”，以为是神经相关的，后来才发现是因为基因列表里混入了嗅觉受体基因……所以一定要回看基因列表的真实情况。

⚠ 最终建议
差异分析结果解读不是一道计算题，它是一道综合题。你需要把统计显著性、表达变化倍数、生物学背景、数据质量——几个维度融合在一起，才能下结论。不要迷信任何一个数值。

好了，这一套流程下来，你已经可以自信地面对任何差异分析结果了。记住：可视化是手段，不是目的。最终你要回答的问题是——这些基因的变化，能解释我的实验假设吗？