# 第 07 章 差异表达分析

## 本章目标

学完本章，你将能够：

1. 用一句话说清差异表达分析（differential expression analysis, DEA）回答什么问题，以及它在整个课题流程中的位置；
2. 整理“表达矩阵 + 分组信息”两类输入，并用 `model.matrix()` 构建设计矩阵；
3. 用 limma 包跑通 log2 转换 → lmFit → eBayes → topTable 的完整流程，并读懂结果表的每一列；
4. 用 `|log2FC| > 1 且 adj.P < 0.05` 的标准筛选差异基因，能解释为什么不能只看 p 值；
5. 把差异基因与焦亡基因集取交集，导出 CSV 结果表，为后续富集分析与预后分析提供输入。

## 正文

### 1. 差异表达分析是什么

在肝细胞癌（hepatocellular carcinoma, HCC）课题中，我们最想知道的一件事是：与癌旁正常组织相比，肿瘤组织里“谁变了”。差异表达分析（differential expression analysis, DEA）就是系统地回答这个问题的统计方法：对每一个基因，比较它在两组样本（如 Normal vs Tumor）中的平均表达水平，检验差异是否显著，并估计变化幅度。

一句话版本：**DEA 给每个基因算两个数——变化幅度（fold change, FC）和显著性（p 值/校正后 p 值），然后按统一标准筛选出“变化大且可信”的基因。**

打个比方，表达谱像一张“基因点名册”，DEA 就是逐个点名、看谁在两组之间“冒头”或“缩头”。差异基因列表是整个课题的“原材料”，后面的功能分析（第 09 章）、网络分析（第 10 章）和生存分析（第 11 章）都建立在它之上。

### 2. 数据准备：表达矩阵与分组信息

DEA 需要两类输入：

- **表达矩阵**：行是基因，列是样本，值是表达量；
- **分组信息**：每个样本属于哪个组（Normal 还是 Tumor）。

本教程课题使用 TCGA-LIHC（肝细胞癌）的 RNA-seq 表达谱，约 371 个肿瘤样本 + 50 个癌旁正常样本（数据下载与整理见第 06 章）。注意：TCGA 的 RNA-seq 原始数据是整数 counts，公开渠道常以 TPM/FPKM 等归一化形式提供。**本章教学演示用 log2 转换后的表达量跑通流程；若手头是 counts，建议用 limma 的 voom() 或 limma-trend 处理**——新手常踩的坑。

在 R 中，分组信息是一个和表达矩阵列顺序一一对应的因子向量：

```r
group <- factor(c(rep("Normal", 50), rep("Tumor", 371)))
```

limma 需要把分组信息转成**设计矩阵**（design matrix），它告诉模型“每个样本属于哪个组、我们要比较哪两组”：

```r
design <- model.matrix(~0 + group)      # 不设截距，每组一个系数
colnames(design) <- levels(group)       # 列名改为 Normal / Tumor
```

新手先记住：`model.matrix()` 就是把分组因子翻译成一张 0/1 矩阵——行是样本，列是组，某样本属于某组则对应位置为 1。真正要比较的两组之差，再用 `makeContrasts()` 显式声明：

```r
contrast.matrix <- makeContrasts(Tumor - Normal, levels = design)
```

### 3. limma 流程：三件套与关键参数

limma（linear models for microarray data）是 Bioconductor 上最常用的差异分析包。它的核心思想是：对每个基因拟合一个线性模型，再用**经验贝叶斯（empirical Bayes）**方法借用全基因组信息来稳定单个基因的方差估计。标准流程是“三件套”：

```r
library(limma)
expr_log2 <- log2(expr + 1)                  # TPM 表达量加 1 再取 log2，避免 log(0)

fit  <- lmFit(expr_log2, design)             # 第一步：对每个基因拟合线性模型
fit2 <- contrasts.fit(fit, contrast.matrix)  # 第二步：取出我们关心的对比 Tumor - Normal
fit2 <- eBayes(fit2)                         # 第三步：经验贝叶斯方差收缩，得到统计量
res  <- topTable(fit2, coef = 1, number = Inf, sort.by = "none")
```

逐步解释：

- `lmFit()`：逐基因做线性回归，得到每个基因在每个组系数下的均值估计与标准误；
- `contrasts.fit()`：把系数组合成我们关心的对比（Tumor 减去 Normal），得到每个基因的 logFC 与标准误；
- `eBayes()`：把全基因组的信息“借”给每个基因，让方差估计更稳定——这是 limma 相对普通 t 检验的核心优势；
- `topTable()`：汇总成结果表。**`number = Inf` 表示输出全部基因**（默认只输出前 10 个，新手常在这里“丢基因”）；`sort.by = "none"` 保持原始基因顺序，方便和表达矩阵对应。

`topTable` 输出的结果表包含以下关键列：

| 列名 | 含义 |
| --- | --- |
| logFC | log2 变化倍数：肿瘤相对癌旁的表达倍数（log2 尺度） |
| AveExpr | 该基因在所有样本中的平均表达量 |
| t | 该对比的 t 统计量 |
| P.Value | 未校正的 p 值 |
| adj.P.Val | BH 校正后的 p 值（即 FDR 估计值） |
| B | 该基因“确实有差异”的对数后验优势（log-odds） |

其中 `logFC` 和 `adj.P.Val` 是我们筛选差异基因最关心的两列。

### 4. 筛选标准：为什么不能只用 p 值

本教程采用的筛选标准（也是生信论文最常用的标准）：

> **|log2FC| > 1 且 adj.P < 0.05**

- `|log2FC| > 1`：变化倍数超过 2 倍（FC > 2 或 FC < 0.5）。这保证差异**有生物学意义**——一个基因只变了 10% 可能不值得关注；
- `adj.P < 0.05`：保证差异**统计上可信**。

为什么不能只看 p 值？因为一次要比较成千上万个基因：即使两组完全没差异，按 p < 0.05 的标准，也会有约 5% 的基因“纯靠运气”显著。假设检测 20,000 个基因，就会随机出现约 1,000 个假阳性。这就是**多重检验问题**（multiple testing）。

解决方法是对 p 值做校正。本教程用 **Benjamini–Hochberg（BH）校正**：把 p 值从小到大排序、逐级调整，得到 adj.P.Val，它控制的是**假发现率**（false discovery rate, FDR）。直观理解：adj.P < 0.05 意味着“被判定为显著的这批基因里，预期约 5% 是假阳性”。R 里的等价写法是 `p.adjust(p, method = "BH")`，topTable 默认已做这一步。

筛选代码：

```r
deg_idx <- abs(res$logFC) > 1 & res$adj.P.Val < 0.05
deg     <- res[deg_idx, ]
table(res$adj.P.Val < 0.05, abs(res$logFC) > 1)   # 交叉计数，看两个条件各自筛掉多少
```

### 5. 结果解读：上调/下调与焦亡基因集交集

拿到差异基因后，先看整体格局：

```r
up   <- deg[deg$logFC > 0, ]   # 上调：肿瘤中表达升高
down <- deg[deg$logFC < 0, ]   # 下调：肿瘤中表达降低
nrow(up); nrow(down)           # 分别计数
```

注意方向约定：`logFC > 0` 表示在 Tumor − Normal 中升高，即**肿瘤中上调**，反之则下调。

下一步回到课题主题——细胞焦亡。把差异基因与“焦亡基因集”取交集，得到“在 HCC 中异常表达的焦亡基因”，这正是后续生存分析与免疫分析的核心对象：

```r
pyroptosis_genes <- c("GSDMD", "GSDME", "CASP1", "IL1B", "IL18",
                      "NLRP3", "AIM2", "PYCARD")   # 示例列表
intersect(rownames(deg), pyroptosis_genes)
```

**重要说明**：焦亡基因集的“权威版本”需要从文献（焦亡机制综述或相关生信论文的补充材料）中获取，不同论文收录的基因从十几个到几十个不等。本教程的列表仅为演示用示例；正式课题中请注明基因集来源并核对基因名拼写（如是否收录 GSDMB、GSDMC、CASP4/5 等）。

### 6. 输出结果：保存 CSV

把结果保存下来，供第 08、09 章继续使用：

```r
dir.create("results", showWarnings = FALSE)
write.csv(res, "results/DEA_TCGA_LIHC_all.csv")    # 全部基因的结果表
write.csv(deg, "results/DEA_TCGA_LIHC_DEGs.csv")   # 筛选后的差异基因
write.csv(intersect(rownames(deg), pyroptosis_genes),
          "results/pyroptosis_DEGs.csv")           # 焦亡相关的差异基因
```

**保存时务必保留“未筛选”的全基因结果表**：第 09 章做 GSEA 时需要对所有基因按 logFC 排序，只用差异基因子集会丢失大量信息。

**完整可运行示例**：用模拟数据跑通上述全部流程（`set.seed` 保证可复现；真实分析时把 `expr`、`group` 换成第 06 章得到的真实数据即可）：

```r
set.seed(2024)
n_gene <- 2000; n_normal <- 30; n_tumor <- 30
expr <- matrix(rnbinom(n_gene * (n_normal + n_tumor), size = 10, mu = 100),
               nrow = n_gene)
rownames(expr) <- paste0("Gene", 1:n_gene)
expr[1:200, (n_normal + 1):(n_normal + n_tumor)] <-   # 模拟 200 个上调基因（约 4 倍）
  expr[1:200, (n_normal + 1):(n_normal + n_tumor)] * 4
group <- factor(c(rep("Normal", n_normal), rep("Tumor", n_tumor)))
# ……随后依次执行第 3 节的 design / limma / 筛选代码即可看到结果……
```

![图06 差异表达火山图（示例，模拟数据）](../figures/fig06_volcano.png)

上图为火山图示例（第 08 章会教怎么画、怎么读）：每个点是一个基因，横轴 log2FC，纵轴 −log10(adj.P.Val)，越靠左右两侧高处，越是“变化大且显著”的差异基因。下图 PCA 图则从样本层面展示 Normal 与 Tumor 的整体分离，说明分组确实能解释表达差异（两张图均为模拟数据示意，读者用真实数据跑脚本后得到自己的图）。

![图08 PCA 图（示例，模拟数据）](../figures/fig08_pca.png)

## 本章小结

- DEA 逐基因比较两组表达，输出变化幅度（logFC）与显著性（adj.P.Val），回答“肿瘤 vs 癌旁哪些基因变了”。
- 输入是表达矩阵 + 分组信息；分组用 `model.matrix()` 转成设计矩阵，用 `makeContrasts()` 声明比较。
- limma 三件套：`lmFit` → `contrasts.fit` → `eBayes`，最后 `topTable` 出表；`number = Inf` 别忘。
- 筛选标准 `|log2FC| > 1 且 adj.P < 0.05`：前者保证幅度，后者（BH 校正）控制多重检验下的假阳性，不能只看 p 值。
- 与焦亡基因集取交集得到课题核心基因；基因集要注明来源。
- 全基因结果表也要保存，第 09 章 GSEA 要用。

## 思考题

1. `topTable` 为什么默认只返回 10 行？`number = Inf` 和 `sort.by = "none"` 各解决什么问题？
2. 如果两组样本量都只有 3 个，limma 的经验贝叶斯校正还可靠吗？为什么？
3. 假设一个基因 p 值极小（1e-10）但 log2FC = 0.2，你会把它算作差异基因吗？为什么？
4. 焦亡基因集来自不同论文时结果可能不同，你认为论文里应该如何报告基因集的选取与版本？
