1. DESeq2简介：RNA-seq差异表达分析背景与DESeq2的设计哲学
大家好，我是这次课的主讲。做RNA-seq数据分析这些年，我接触过不少工具。但说实话，DESeq2一直是我工具箱里的首选。为什么？因为它稳。嗯，咱们先从背景聊起。
1.1 RNA-seq差异表达分析，到底在解决什么问题？
你想想看，一个RNA-seq实验，最终拿到的是一堆测序读数。经过比对和定量，我们得到一个计数矩阵——行是基因，列是样本，每个格子是一个整数，代表某个基因在某个样本里检测到的片段数。
问题是：这些计数值，能直接拿来比大小吗？
不能。原因有三个：
- 测序深度不一样——有的样本总reads多，有的少。直接比计数，深度高的样本所有基因看起来都上调。
- 基因长度不一样——长基因更容易被测到，但这不代表它表达量更高。
- 技术噪声大——生物学重复之间的波动，比普通实验大得多。低表达基因尤其明显。
所以，差异表达分析的核心任务就三个：
- 标准化——消除测序深度和基因长度带来的偏差
- 建模——用合适的统计模型描述计数数据的分布
- 假设检验——判断某个基因在处理组和对照组之间，表达量变化是否显著
我刚开始做RNA-seq时，用过老式的倍数变化直接筛选。结果呢？假阳性率高得吓人。后来才明白，统计建模这一步省不了。
1.2 DESeq2的设计哲学：稳健、严谨、可解释
DESeq2是Love等人2014年在Genome Biology上发表的。它脱胎于早期的DESeq，但做了大量改进。我个人习惯把它的设计哲学归纳为三点：
| 设计原则 | 具体做法 | 为什么重要 |
|---|---|---|
| 负二项分布建模 | 用负二项分布而非泊松分布描述计数数据 | 泊松要求均值等于方差，但RNA-seq数据方差远大于均值。负二项多了一个离散度参数，能更好拟合真实情况。 |
| 收缩估计 | 利用所有基因的信息，对每个基因的离散度进行收缩 | 低表达基因的方差估计极不稳定。收缩相当于借力——用全局趋势稳定局部估计。我在项目中遇到过一个小样本实验，没有收缩的话，差异基因列表几乎全是噪声。 |
| 对数倍数变化的收缩 | 对低表达、高变异的基因，强行拉低其倍数变化 | 这步很多人不理解。确实，它会降低部分基因的fold change。但换来的是排名更可靠，假阳性更少。说白了，我们宁可得不到一个巨大的倍数变化，也要确保这个变化是真的。 |
核心一句话：DESeq2不做简单的标准化，它用统计模型来"理解"数据的结构，然后把真实信号从噪声中剥离出来。
1.3 从计数矩阵到差异结果：一条清晰的流水线
下面这张图，是我每次讲DESeq2必画的。它把整个分析流程串起来了。你先看，我再解释。
我解释一下这张图：
- 第一步：输入原始计数矩阵。注意，DESeq2要求整数计数，不接受FPKM或TPM。我曾经被人问过能不能用RPKM跑DESeq2，答案是不行——那个标准化方式不一样。
- 第二步：构建DESeqDataSet对象。这一步你提供counts、实验设计公式和样本信息。实验设计这一步很关键，写错了公式，后面所有结果都跑偏。
- 第三步：核心建模。估计大小因子做标准化，然后估计每个基因的离散度，最后做收缩。这里有几个参数可以调，比如
fitType、minReplicatesForReplace，会影响后面的结果。 - 第四步：拟合负二项模型，做Wald检验，给出p值和校正后的padj。
- 第五步：输出结果表。你用它来做火山图、MA图、热图，都是基于这里产出的log2FoldChange和padj。
1.4 为什么参数调整会影响可视化？
你做一张火山图，x轴是log2FoldChange，y轴是-log10(padj)。那问题来了：
- 如果收缩强度不同，同一个基因的log2FoldChange可能差0.5甚至1倍。它在火山图上的位置就变了——可能从显著下调变成不显著。
- 如果离散度估计的方法不同，p值分布会整体偏移，显著基因的数量可能翻倍或腰斩。
- 如果独立过滤阈值不同，padj校正后的基因列表可以差几十到上百个。
所以，参数调整不是单纯的数学游戏。它直接决定了你看到的图长什么样，以及你最后汇报的差异基因列表是哪1000个还是哪500个。
我个人习惯：做正式分析前，先跑一遍默认参数，看看数据的整体质量。然后再根据实验设计调一两处关键参数，比较输出表的变化。千万不要一次调七八个参数，到时候变了都不知道哪个参数引起的。
1.5 一段最简的R代码
按惯例，看一段核心代码。你感受一下DESeq2的接口有多简洁：
library(DESeq2)
# 构建数据对象
dds <- DESeqDataSetFromMatrix(
countData = counts_matrix,
colData = sample_info,
design = ~ condition
)
# 核心分析
dds <- DESeq(dds)
# 提取结果
res <- results(dds, contrast = c("condition", "treated", "control"))
如果你跑过这段代码，你会发现它几乎没有选项。那是因为DESeq把大部分默认参数设得很合理。但恰恰是那些没写出来的参数，比如betaPrior、independentFiltering、cooksCutoff，才是我们后面要重点讲解的。
注意：千万不要以为跑通默认参数就万事大吉了。对于小样本、高噪声、低表达基因占比大的数据，默认参数不一定最优。我见过有人拿3个生物学重复跑，出来2000个差异基因，结果复查时发现大部分是假阳性——离散度没有收缩好。
1.6 本章小结
到这里，我对DESeq2的背景和设计哲学做了一个梳理。你至少要知道：它为什么用负二项分布，哪些参数可以调，以及这些调整最终会怎么影响你画的图。后面的章节，我们会逐个参数拆解，看它们到底在做什么，以及影响有多大。
好了，这一章就到这里。记住一句话：用好DESeq2，不是会跑代码，而是理解每一步在干什么。