差异表达分析的统计基础：假设检验与多重检验校正

做差异表达分析，说白了就是在问一个问题：某个基因在两组间的表达差异，到底是真实的生物学差异，还是随机波动造成的？

这个问题，统计学家早就帮我们搭好了框架。就是假设检验。我刚开始接触RNA-seq时，觉得这步就是个"走流程"，直到有一次自己踩了坑——没做多重检验校正，结果搞出来一堆假阳性基因，验证实验全废了……从那以后，我再也不敢小看这一步。

这一章，我们就来拆解DESeq2背后的统计逻辑。搞懂了这些，你才能明白为什么DESeq2的输出里有 pvalue 和 padj 两列，以及到底该信哪一个。

差异表达分析统计流程 基因表达定量 Count 矩阵 统计假设检验 Wald 检验 原始 p 值 每个基因一个 多重检验校正 BH / FDR 差异 基因 ⚠ 关键认知 不做多重检验校正 → 假阳性激增（2 万个基因同时检验，纯随机也能出 1000 个 "显著"） 零假设 H₀：基因在两组间无差异 备择假设 H₁：基因在两组间有差异

1. 假设检验：统计决策的基本框架

假设检验的思路其实很直观。你先设立一个"零假设"（H₀），也就是"这个基因在两组间没有差异"。然后拿数据去检验，看看能不能推翻它。

如果能推翻，就接受备择假设（H₁）——基因有差异。推翻的依据是什么？就是 p 值。p 值表示的是：在零假设成立的前提下，观察到当前结果（或更极端结果）的概率。

通俗理解：

  • p 值小（比如 < 0.05） → 零假设不太可能成立 → 拒绝 H₀ → 基因显著差异
  • p 值大 → 没有足够证据拒绝 H₀ → 无法判断显著

但注意，p 值不是"基因差异的概率"。这是新手最容易误解的地方。我见过不少人在论文里写 "p=0.03 说明基因有 97% 的概率是差异表达的"——这是错的。p 值衡量的是数据与零假设的兼容性，不是效应概率。

💡 我的个人习惯：在汇报差异基因数量时，我一般会把 p 值和 padj 都列出来，让审稿人看到校正前后的对比。这既是透明的做法，也能体现你对统计方法的把控。

2. DESeq2 采用的检验策略

DESeq2 没有用传统的 t 检验。为什么？因为 RNA-seq 数据是计数数据（整数），而且方差往往比均值大——存在"过度离散"问题。你想想看，同一种处理条件下，几个生物学重复的 count 数可能波动很大，直接用 t 检验会低估方差，导致假阳性。

DESeq2 用的是 负二项分布模型（Negative Binomial），也叫 Gamma-Poisson 分布。这个模型有两个参数：均值和离散度。离散度参数专门用来处理过度离散的情况。

# 负二项模型示意（非实际代码）
# 对于基因 i，样本 j：
# K_ij ~ NB(μ_ij, α_i)
# 其中 μ_ij 是期望表达量，α_i 是该基因的离散度参数

具体到检验环节，DESeq2 默认使用 Wald 检验。简单说，就是用基因的 log2 倍数变化（log2FC）除以它的标准误，得到一个 z 统计量，再算出 p 值。

# Wald 检验统计量
# W = β / SE(β)
# β 是 log2FC 估计值，SE 是标准误
# p 值来自标准正态分布的双侧检验

我记得有一次项目里，用户自己用 t 检验跑了同样数据，结果差异基因数量是 DESeq2 的 3 倍多。后来越挖越深才发现，那些多出来的基因基本都是假阳性——尤其是低表达基因，方差被严重低估了。

3. 多重检验校正：为什么不能只看 p 值

到现在为止，我们拿到的还是"原始 p 值"。但问题来了——

一个典型的 RNA-seq 实验会同时检测大约 20,000 个基因。如果每个基因都用 p < 0.05 来判断显著性，即使两组的表达完全一样（零假设全部成立），你也预计会发现约 1,000 个"显著"基因。

这就是 多重检验问题。检验次数越多，假阳性累积越严重。

⚠ 我曾经踩过的坑：早期做项目时，我直接用 p < 0.01 筛基因，筛出 800 多个差异基因。组学验证做了 20 个，只有 3 个能重复出来……从那以后，我每次都会先看一眼 padj 列。

校正的核心思路很简单：控制"假阳性发现"的比例。主要有两个方向：

  • 控制 I 类错误率（FWER）：如 Bonferroni 校正，非常严格
  • 控制假发现率（FDR）：如 BH 方法，更常用，也是 DESeq2 的默认方法

4. 常用校正方法对比

我整理了一张对比表，方便你看清差异：

方法 控制目标 严格程度 适用场景
Bonferroni FWER（至少一个假阳性） 极高 候选基因极少、要求零假阳性
Benjamini-Hochberg (BH) FDR（假阳性比例） 中等 转录组差异分析（DESeq2 默认）
BY (Benjamini-Yekutieli) FDR（依赖相关时） 较高 基因间存在强相关结构
q-value (Storey) FDR（基于估计） 略宽松 希望最大化检验效力时

在 DESeq2 中，默认采用 BH 方法。它的做法是：将所有基因按 p 值从小到大排序，然后逐个计算调整后的 p 值（padj），使得在 padj < 0.05 的基因集合中，假阳性比例不超过 5%。

# R 中的 p 值调整示例（DESeq2 内部逻辑）
pvals <- c(0.001, 0.005, 0.01, 0.03, 0.05, 0.08, 0.20, 0.50)

# Bonferroni
padj_bonf <- p.adjust(pvals, method = "bonferroni")

# BH
padj_bh   <- p.adjust(pvals, method = "BH")

# 对比结果
cbind(pvalue = pvals, bonferroni = padj_bonf, BH = padj_bh)

关键结论：Bonferroni 在基因数量大时过于保守，容易漏掉真实差异基因。BH 方法在敏感性和假阳性控制之间取得了平衡，是差异分析的首选。

5. 在 DESeq2 中选择合适的阈值

DESeq2 的结果表中，你至少要看三列：log2FoldChange、pvalue 和 padj。筛选差异基因时，我建议以 padj 为准。

常用的筛选标准：

  • padj < 0.05 （默认，适合大多数场景）
  • padj < 0.01 （更严格，适合后续实验验证成本高的情况）
  • |log2FC| > 1 且 padj < 0.05 （加上倍数变化门槛，减少生物无关的微弱差异）
# DESeq2 结果提取示例
res <- results(dds, alpha = 0.05)
summary(res)

# 筛出显著基因
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig_genes)

⚠ 注意：不要先按 p 值筛选，再挑 log2FC 大的，那样会引入选择偏倚。正确的做法是将 padj 和 log2FC 的阈值同时用在 results() 或 subset() 中。

我个人习惯是，先把 padj < 0.05 作为初筛，然后根据实验需求决定是否再收紧。如果是要做下游功能富集分析，稍微宽松一点（比如 padj < 0.1）也能接受，因为富集分析本身有额外的稳健性。

← 上一章 📖 返回目录 下一章 →