第3讲：差异表达分析

从计数到差异基因

本讲概述

学习目标： - 理解RNA-seq数据的统计特性（负二项分布） - 掌握DESeq2分析流程 - 理解多重检验校正原理 - 掌握结果筛选和可视化

重点难点： - 🔑 重点：DESeq2分析步骤、结果筛选标准 - ⚠️ 难点：离散度估计、log2FC收缩原理

1. RNA-seq数据统计特性

1.1 计数数据特点

特性	说明
离散型	非负整数
过度离散	方差 > 均值
零膨胀	大量基因在某些样本中为0

1.2 为什么不能用t检验？

正态分布：均值 = 方差
RNA-seq：  方差 >> 均值（过度离散）

1.3 负二项分布

均值 = μ
方差 = μ + α×μ²  (α为离散度参数)

离散度	分布	方差 vs 均值
α = 0	泊松分布	方差 = 均值
α > 0	负二项分布	方差 > 均值

实际应用： - 低表达基因：离散度高（采样误差大） - 高表达基因：离散度低（采样误差小）

2. DESeq2工作流程 🔑

2.1 核心原理

负二项广义线性模型（GLM）：

log2(q_ij) = β_0 + β_1 × condition + ε

β_0：截距（对照组基线）
β_1：处理效应（log2 Fold Change）

2.2 分析步骤

原始counts矩阵
      ↓
  计算Size Factors（标准化）
      ↓
  估计基因离散度（趋势拟合）
      ↓
  拟合负二项GLM
      ↓
  Wald检验/LRT
      ↓
  多重检验校正（BH）
      ↓
  差异基因筛选

2.3 DESeq2代码

# 1. 创建DESeq2对象
dds <- DESeqDataSetFromMatrix(
  countData = counts_matrix,
  colData = sample_info,
  design = ~ condition
)

# 2. 过滤低表达基因
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

# 3. 运行DESeq2（包含所有步骤）
dds <- DESeq(dds)

# 4. 提取结果
res <- results(dds, contrast = c("condition", "treat", "ctrl"))

2.4 Size Factor标准化

中位数比值法：

size factor_j = median_i (counts_ij / geometric_mean_i)

假设：大多数基因在样本间表达稳定

2.5 离散度估计 ⚠️

DESeq2采用趋势拟合： - 先计算每个基因的离散度 - 再拟合离散度-均值趋势曲线 - 最后将各基因离散度向趋势收缩

为什么需要收缩？ - 小样本时，单个基因的离散度估计不稳定 - 收缩可提高估计的稳健性

3. edgeR简介

特性	DESeq2	edgeR
离散度估计	趋势拟合	经验贝叶斯
收缩性	更强	较弱
小样本稳健性	更好	中等
速度	较慢	较快

推荐： - 样本少（n=3）：DESeq2 - 复杂设计：两者均可

4. 多重检验校正 🔑

4.1 为什么需要校正？

检测20,000个基因（α=0.05）：

假阳性期望数 = 20,000 × 0.05 = 1,000个！

4.2 Benjamini-Hochberg (BH) 方法

控制假发现率（FDR）：

FDR = E[假阳性数 / 拒绝数]

找到最大的k，使得 pₖ ≤ (k/m) × α

4.3 p值 vs FDR

指标	符号	使用场景
p值	pvalue	不推荐直接使用
调整p值	padj	差异基因筛选

警告

必须使用padj筛选差异基因，不能用原始pvalue！

5. 结果解读 🔑

5.1 DESeq2结果表

列名	含义	说明
baseMean	平均表达量	标准化后的平均reads
log2FoldChange	log2倍数变化	处理/对照
lfcSE	log2FC标准误	估计不确定性
stat	Wald统计量	检验统计量
pvalue	原始p值	未校正
padj	校正p值	BH校正后

5.2 差异基因筛选标准

推荐阈值：
- |log2FC| > 1  （倍数变化 > 2）
- padj < 0.05  （FDR < 5%）

分类矩阵：

           padj < 0.05
            Yes       No
       ┌─────────┬─────────┐
log2FC │  上调   │ nominally│
  > 1  │  (Up)   │   上调   │
       ├─────────┼─────────┤
-1 to 1│ 无差异  │  无差异  │
       ├─────────┼─────────┤
log2FC │  下调   │ nominally│
 < -1  │ (Down)  │   下调   │
       └─────────┴─────────┘

5.3 log2FC收缩 ⚠️

问题：低表达基因的高log2FC不可靠

DESeq2解决方案：

resLFC <- lfcShrink(dds, coef = "condition_treat_vs_ctrl")

效果： - 低表达基因：log2FC向0收缩 - 高表达基因：保持相对稳定

6. 结果可视化

6.1 MA图

X轴：平均表达量（log10）
Y轴：log2 Fold Change
特点：红点为差异基因，低表达基因log2FC波动大

plotMA(res, ylim = c(-5, 5))

6.2 火山图

X轴：log2 Fold Change
Y轴：-log10(padj)
特点：同时展示变化倍数和统计显著性

6.3 热图

展示Top差异基因的表达模式，需使用标准化后的数据

7. 常见问题

问题	可能原因	解决方案
无差异基因	分组错误	检查colData
	批次效应	添加批次到设计公式
差异基因过多	生物学差异大	检查样本质量

8. 本讲小结

要点	内容
统计基础	RNA-seq数据符合负二项分布，存在过度离散
DESeq2核心	负二项GLM + 中位数比值标准化 + Wald检验
多重校正	使用BH方法控制FDR
筛选标准	\|log2FC\| > 1 + padj < 0.05
结果解读	关注表达量、变化倍数、统计显著性

思考题

为什么RNA-seq数据不能使用t检验？
Size factor标准化与TPM标准化的区别？
如何解释padj > 0.05但pvalue < 0.05的基因？
火山图和MA图各自的优势是什么？

参考文献

Love MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014
Robinson MD, et al. edgeR: a Bioconductor package for differential expression analysis. Bioinformatics. 2010

--- title: "第3讲：差异表达分析" subtitle: "从计数到差异基因" --- ## 本讲概述 **学习目标**： - 理解RNA-seq数据的统计特性（负二项分布） - 掌握DESeq2分析流程 - 理解多重检验校正原理 - 掌握结果筛选和可视化 **重点难点**： - 🔑 **重点**：DESeq2分析步骤、结果筛选标准 - ⚠️ **难点**：离散度估计、log2FC收缩原理 --- ## 1. RNA-seq数据统计特性 ### 1.1 计数数据特点 | 特性 | 说明 | |------|------| | **离散型** | 非负整数 | | **过度离散** | 方差 > 均值 | | **零膨胀** | 大量基因在某些样本中为0 | ### 1.2 为什么不能用t检验？ ``` 正态分布：均值 = 方差 RNA-seq：方差 >> 均值（过度离散） ``` ### 1.3 负二项分布 ``` 均值 = μ 方差 = μ + α×μ² (α为离散度参数) ``` | 离散度 | 分布 | 方差 vs 均值 | |--------|------|--------------| | α = 0 | 泊松分布 | 方差 = 均值 | | α > 0 | 负二项分布 | 方差 > 均值 | **实际应用**： - 低表达基因：离散度高（采样误差大） - 高表达基因：离散度低（采样误差小） --- ## 2. DESeq2工作流程 🔑 ### 2.1 核心原理 **负二项广义线性模型（GLM）**： ``` log2(q_ij) = β_0 + β_1 × condition + ε β_0：截距（对照组基线） β_1：处理效应（log2 Fold Change） ``` ### 2.2 分析步骤 ``` 原始counts矩阵 ↓ 计算Size Factors（标准化） ↓ 估计基因离散度（趋势拟合） ↓ 拟合负二项GLM ↓ Wald检验/LRT ↓ 多重检验校正（BH） ↓ 差异基因筛选 ``` ### 2.3 DESeq2代码 ```r # 1. 创建DESeq2对象 dds <- DESeqDataSetFromMatrix( countData = counts_matrix, colData = sample_info, design = ~ condition ) # 2. 过滤低表达基因 keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep, ] # 3. 运行DESeq2（包含所有步骤） dds <- DESeq(dds) # 4. 提取结果 res <- results(dds, contrast = c("condition", "treat", "ctrl")) ``` ### 2.4 Size Factor标准化 **中位数比值法**： ``` size factor_j = median_i (counts_ij / geometric_mean_i) ``` 假设：大多数基因在样本间表达稳定 ### 2.5 离散度估计 ⚠️ DESeq2采用**趋势拟合**： - 先计算每个基因的离散度 - 再拟合离散度-均值趋势曲线 - 最后将各基因离散度向趋势收缩 **为什么需要收缩？** - 小样本时，单个基因的离散度估计不稳定 - 收缩可提高估计的稳健性 --- ## 3. edgeR简介 | 特性 | DESeq2 | edgeR | |------|--------|-------| | 离散度估计 | 趋势拟合 | 经验贝叶斯 | | 收缩性 | 更强 | 较弱 | | 小样本稳健性 | **更好** | 中等 | | 速度 | 较慢 | 较快 | **推荐**： - 样本少（n=3）：**DESeq2** - 复杂设计：两者均可 --- ## 4. 多重检验校正 🔑 ### 4.1 为什么需要校正？检测20,000个基因（α=0.05）： ``` 假阳性期望数 = 20,000 × 0.05 = 1,000个！ ``` ### 4.2 Benjamini-Hochberg (BH) 方法 **控制假发现率（FDR）**： ``` FDR = E[假阳性数 / 拒绝数] ``` 找到最大的k，使得 pₖ ≤ (k/m) × α ### 4.3 p值 vs FDR | 指标 | 符号 | 使用场景 | |------|------|----------| | p值 | pvalue | **不推荐直接使用** | | 调整p值 | **padj** | **差异基因筛选** | ::: {.callout-warning} **必须使用padj筛选差异基因，不能用原始pvalue！** ::: --- ## 5. 结果解读 🔑 ### 5.1 DESeq2结果表 | 列名 | 含义 | 说明 | |------|------|------| | baseMean | 平均表达量 | 标准化后的平均reads | | **log2FoldChange** | log2倍数变化 | 处理/对照 | | lfcSE | log2FC标准误 | 估计不确定性 | | stat | Wald统计量 | 检验统计量 | | pvalue | 原始p值 | 未校正 | | **padj** | 校正p值 | **BH校正后** | ### 5.2 差异基因筛选标准 ``` 推荐阈值： - |log2FC| > 1 （倍数变化 > 2） - padj < 0.05 （FDR < 5%） ``` **分类矩阵**： ``` padj < 0.05 Yes No ┌─────────┬─────────┐ log2FC │ 上调 │ nominally│ > 1 │ (Up) │ 上调 │ ├─────────┼─────────┤ -1 to 1│ 无差异 │ 无差异 │ ├─────────┼─────────┤ log2FC │ 下调 │ nominally│ < -1 │ (Down) │ 下调 │ └─────────┴─────────┘ ``` ### 5.3 log2FC收缩 ⚠️ **问题**：低表达基因的高log2FC不可靠 **DESeq2解决方案**： ```r resLFC <- lfcShrink(dds, coef = "condition_treat_vs_ctrl") ``` **效果**： - 低表达基因：log2FC向0收缩 - 高表达基因：保持相对稳定 --- ## 6. 结果可视化 ### 6.1 MA图 - **X轴**：平均表达量（log10） - **Y轴**：log2 Fold Change - **特点**：红点为差异基因，低表达基因log2FC波动大 ```r plotMA(res, ylim = c(-5, 5)) ``` ### 6.2 火山图 - **X轴**：log2 Fold Change - **Y轴**：-log10(padj) - **特点**：同时展示变化倍数和统计显著性 ### 6.3 热图展示Top差异基因的表达模式，需使用标准化后的数据 --- ## 7. 常见问题 | 问题 | 可能原因 | 解决方案 | |------|----------|----------| | 无差异基因 | 分组错误 | 检查colData | | | 批次效应 | 添加批次到设计公式 | | 差异基因过多 | 生物学差异大 | 检查样本质量 | --- ## 8. 本讲小结 | 要点 | 内容 | |------|------| | 统计基础 | RNA-seq数据符合**负二项分布**，存在过度离散 | | DESeq2核心 | 负二项GLM + **中位数比值标准化** + Wald检验 | | 多重校正 | 使用**BH方法控制FDR** | | 筛选标准 | **\|log2FC\| > 1 + padj < 0.05** | | 结果解读 | 关注表达量、变化倍数、统计显著性 | --- ## 思考题 1. 为什么RNA-seq数据不能使用t检验？ 2. Size factor标准化与TPM标准化的区别？ 3. 如何解释padj > 0.05但pvalue < 0.05的基因？ 4. 火山图和MA图各自的优势是什么？ --- ## 参考文献 1. Love MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014 2. Robinson MD, et al. edgeR: a Bioconductor package for differential expression analysis. Bioinformatics. 2010