差异表达：统计模型与假设检验

所属板块 分析方向与案例

把基础对象与算法方法重新放回真实分析任务与工作流。

适合谁读 建议在以下阶段阅读

已理解 RNA-seq 定量流程，希望掌握差异表达统计建模原理的读者

建议起点 推荐阅读路径

从 count 数据的统计分布特性出发，理解负二项分布模型、离散度估计与多重检验校正的算法基础。

快速概览

差异表达分析的核心是构建统计模型来判断观察到的表达差异是否超出随机波动范围。本节介绍基于负二项分布的 count 数据建模、离散度估计方法与多重检验校正算法。

• 统计模型：为什么 RNA-seq 差异分析使用负二项分布而非正态分布
• 离散度估计：技术噪声与生物学变异的分离
• 多重检验：FDR 控制的算法原理（Benjamini-Hochberg 方法）

前置知识

问题定义

差异表达检测问题：给定两个或多个生物学条件下的 RNA-seq count 数据，识别表达水平发生系统性变化的基因。

• 输入：

  • 基因表达计数矩阵 $C = [c_{ij}]$，$c_{ij}$ 为基因 $i$ 在样本 $j$ 中的 counts
  • 实验设计矩阵 $X$，描述样本分组、批次、协变量
  • 基因长度信息（可选，用于某些归一化方法）

• 输出：

  • 每个基因的 log fold change 估计值 $\hat{\beta}_i$
  • 统计显著性度量（p-value、调整后的 p-value / FDR）
  • 差异表达基因列表（满足 FDR 阈值的基因集合）

• 核心难点：

  • Count 数据的离散分布特性（方差 > 均值）
  • 样本间技术偏差与生物学变异的分离
  • 多重假设检验中的假阳性控制

为什么从 Counts 出发：统计分布基础

RNA-seq 数据的统计特性

RNA-seq counts 服从离散概率分布，具有以下特征：

1. 非负整数：counts $\in \{0, 1, 2, ...\}$
2. 过度离散（Overdispersion）：方差显著大于均值
3. 均值-方差依赖：高表达基因的绝对波动更大

泊松分布的局限

简单模型假设 counts 服从泊松分布：

$c_{ij} \sim \text{Poisson}(\lambda_{ij})$

其中 $\text{Var}(c_{ij}) = \lambda_{ij} = \mathbb{E}[c_{ij}]$。

问题：泊松分布假设方差等于均值，但 RNA-seq 数据中观察到：

$\text{Var}(c_{ij}) \gg \mathbb{E}[c_{ij}]$

这种过度离散来源于生物学重复间的真实变异。

负二项分布模型

更合适的模型是负二项分布（Negative Binomial, NB）：

$c_{ij} \sim \text{NB}(\mu_{ij}, \phi_i)$

其中：

• $\mu_{ij}$：基因 $i$ 在样本 $j$ 中的期望 counts
• $\phi_i$：离散度参数（dispersion），控制方差大小

负二项分布的方差为：

$\text{Var}(c_{ij}) = \mu_{ij} + \phi_i \cdot \mu_{ij}^2$

• 当 $\phi_i \rightarrow 0$：退化为泊松分布
• 当 $\phi_i > 0$：允许方差大于均值（过度离散）

广义线性模型（GLM）框架

模型构建

差异表达分析使用广义线性模型来关联实验条件与表达水平：

$\log(\mu_{ij}) = \log(s_j) + x_j^T \beta_i$

其中：

• $s_j$：样本 $j$ 的大小因子（size factor，归一化因子）
• $x_j$：样本 $j$ 的设计矩阵行向量
• $\beta_i$：基因 $i$ 的系数向量（包含 log fold change）

对数线性关系

模型假设条件对数期望的线性关系：

$\log \mathbb{E}[c_{ij} | x_j] = \alpha_i + \log s_j + \sum_k x_{jk} \beta_{ik}$

参数解释：

• $\alpha_i$：基因 $i$ 的基础表达水平（截距）
• $\beta_{ik}$：条件 $k$ 相对于对照的 log fold change
• $s_j$：样本特异性归一化因子（DESeq2 使用 median ratio 方法估计）

假设检验

对于两个条件（control vs treatment），检验：

$H_0: \beta_i = 0 \quad \text{vs} \quad H_1: \beta_i \neq 0$

使用 Wald 检验 或 似然比检验（LRT）计算 p-value。

Wald 统计量：

$W_i = \frac{\hat{\beta}_i}{\text{SE}(\hat{\beta}_i)} \sim \mathcal{N}(0, 1) \text{ (under } H_0)$

离散度估计：算法核心

为什么需要估计离散度

离散度参数 $\phi_i$ 决定了统计检验的灵敏度：

• 低估离散度 $\rightarrow$ 假阳性增加（把噪声当信号）
• 高估离散度 $\rightarrow$ 假阴性增加（错过真实差异）

估计策略

方法 1：基因特异性估计（Gene-wise）

使用每个基因的数据独立估计：

$\hat{\phi}_i^{\text{gene}} = \frac{s_i^2 - \bar{c}_i}{\bar{c}_i^2}$

其中 $s_i^2$ 和 $\bar{c}_i$ 为样本方差和均值。

问题：低表达基因样本数少，估计不稳定。

方法 2：平均离散度（Mean-based）

假设具有相似平均表达水平的基因具有相似离散度，将基因按平均表达量分 bin，在每个 bin 内共享离散度估计。

方法 3：经验贝叶斯收缩（DESeq2/edgeR 采用）

算法核心：将基因特异性估计向全局趋势收缩：

$\hat{\phi}_i^{\text{shrunk}} = \frac{n_i \cdot \hat{\phi}_i^{\text{gene}} + m \cdot \hat{\phi}^{\text{trend}}}{n_i + m}$

其中：

• $n_i$：基因 $i$ 的有效样本量
• $m$：收缩强度参数
• $\hat{\phi}^{\text{trend}}$：平均表达-离散度趋势线估计

算法复杂度：

• 时间：$O(G \times n)$，$G$ 为基因数，$n$ 为样本数
• 空间：$O(G)$

多重检验校正：FDR 控制算法

问题背景

当同时对 $G \approx 20{,}000$ 个基因进行假设检验时，即使所有基因都无真实差异，在 $\alpha = 0.05$ 水平下仍预期有：

$G \times \alpha = 1{,}000 \text{ 个假阳性}$

Benjamini-Hochberg FDR 控制

算法步骤：

1. 将所有基因按 p-value 从小到大排序：$p_{(1)} \leq p_{(2)} \leq ... \leq p_{(G)}$

2. 找到最大的 $k$ 使得：

$p_{(k)} \leq \frac{k}{G} \cdot \alpha$

3. 拒绝所有 $i \leq k$ 的假设（即标记为显著差异）

输出：每个基因的调整 p-value（FDR）

算法复杂度：$O(G \log G)$，主要由排序步骤决定。

独立筛选（Independent Filtering）

策略：在 FDR 校正前，过滤掉低表达基因以提高检验效能。

算法原理：

如果基因在所有样本中的 counts 都低于阈值 $c_{\min}$，则认为其缺乏足够信息支持差异检测。

过滤后的基因数 $G' < G$，FDR 控制更严格：

$p_{(k)} \leq \frac{k}{G'} \cdot \alpha$

算法实例：假设检验过程

假设有两组 RNA-seq 样本：control 与 treatment。某个基因在 treatment 组 counts 明显更高，统计建模需要回答：

1. 测序深度归一化：两组样本的 $\log s_j$ 是否一致？
2. 离散度估计：该基因在组内重复之间的变异 $\phi_i$ 有多大？
3. 假设检验：在 NB-GLM 模型下，$\beta_i$ 是否显著偏离 0？
4. 多重检验校正：在同时检验 20,000 基因的背景下，这个差异是否还能在 FDR 控制后保留？

只有这些问题都被纳入模型，这个”更高”才更接近真正可报告的差异表达信号。

算法实例：DESeq2 工作流程

标准分析流程

Raw counts matrix
    ↓
1. 归一化因子估计（median ratio method）
    ↓
2. 离散度估计（empirical Bayes shrinkage）
    ↓
3. GLM 拟合（negative binomial）
    ↓
4. Wald 检验 / LRT
    ↓
5. Independent filtering
    ↓
6. Benjamini-Hochberg FDR correction
    ↓
Differential expression results

关键参数设置

参数	作用	默认值
alpha	FDR 控制阈值	0.1
lfcThreshold	最小 log2 fold change	0
altHypothesis	双侧或单侧检验	”greaterAbs”
minReplicatesForReplace	异常值替换所需重复数	7

检验效能与实验设计

样本量计算

差异表达分析的统计效能取决于：

1. 效应大小：log fold change 越大，越容易检测
2. 离散度：生物学变异越小，效能越高
3. 样本量：重复数 $n$ 增加，效能提升

经验法则：

• 对于中等效应（LFC = 1）和中等变异，每组 $n = 3$ 可检测约 50-60% 的真阳性
• 每组 $n = 6$ 可提升至 80-90%
• 稀有转录本或高变异基因需要 $n \geq 10$

批次效应控制

设计矩阵扩展：

$\log \mu_{ij} = \log s_j + \text{condition}_j \cdot \beta_i + \text{batch}_j \cdot \gamma_i$

在模型中显式包含批次变量，可分离技术变异与生物学信号。

与上下游分析的连接

差异表达分析在 RNA-seq 工作流中的位置：

• 上游依赖：稳定的定量结果、正确的注释版本、一致的 gene/transcript 层级
• 下游应用：功能富集分析、通路分析、可视化验证

一个”差异表达结果不合理”的问题，根源可能来自：

• 前处理或 QC 不足
• counts 构建方式不合适
• 样本设计存在混杂因素
• transcript 层与 gene 层目标混用

常见误区

• 混淆统计显著性与生物学重要性：
• p-value 小 ≠ 生物学意义大
• 小 fold change（LFC < 0.5）即使统计显著，可能无生物学功能
• 建议同时考虑效应大小和显著性
• 忽略模型假设检验：
• DESeq2/edgeR 假设：

相关算法与扩展阅读

• RNA-seq 工作流概览：差异分析在完整流程中的位置
• 基因层与转录本层表达定量：DE 分析的输入层级选择
• TPM、FPKM、CPM 与有效长度：为什么 DE 从 counts 而非 TPM 出发
• 数据库与资源：差异基因的功能注释与通路分析

相关页面

流程入口

RNA-seq 工作流概览

先看整个 RNA-seq 主线，再把差异表达放回 workflow 中理解。

延伸阅读

单位与归一化

TPM、FPKM、CPM 与有效长度

先区分 counts 与 TPM 的角色，再理解为什么 DE 常从 counts 出发。

延伸阅读

分析层级

基因层与转录本层表达定量

差异表达前要先明确你的分析单位是 gene 还是 transcript。

延伸阅读

解释层

常用数据库与资源

差异基因列表的下游解释通常要回到数据库、通路与公开资源。

延伸阅读