1. 认识DESeq2：差异表达分析的核心思想与参数概览

差异表达分析，说白了就是要找出两组样本里哪些基因的表达量真的变了。RNA-seq 数据一出来就是几十个G的测序片段，但真正做分析时，我们面对的是一个计数矩阵——每个基因在每个样本里有多少个读数。你想想看，这个数字本身并不直接告诉我们哪个基因上调了、哪个下调了。它受很多因素影响：测序深度、基因长度、文库复杂度……甚至同一批样本里，同一个基因的计数也能差好几倍。

DESeq2 是这几年最主流的差异表达工具之一。我个人习惯用它来跑所有 RNA-seq 项目，因为它对低表达基因的处理特别稳。我记得最早我刚入行时，用过 edgeR 和 limma，但后来发现 DESeq2 的方差模型更贴合真实数据——它从根本上解决了“低表达基因为什么看起来差异那么大”这个问题。好，我们直接进正题。

1.1 核心思想：方差均值关系与负二项分布

为什么不用泊松分布？很多入门教程会告诉你“RNA-seq 计数符合泊松”，但实际项目中你会发现，基因表达不是那么乖的。泊松分布假设方差等于均值——可现实中低表达基因的方差往往远大于均值，这叫「过离散」（overdispersion）。

我曾在一次肿瘤 vs 正常组织的项目里，看到某个基因在正常组平均只有5个读数，但肿瘤组里有20个——你看着觉得很显著吧？结果 DESeq2 把它的 p 值调得很高。为什么？因为它发现这个基因在组内的波动太大了，5 和 20 不过是随机波动的错觉。这就是负二项分布出场的原因：它多了一个离散参数，允许方差比均值大。

DESeq2 的核心假设是：

• 每个基因的计数服从 负二项分布（NB）
• 均值由实验条件和测序深度决定
• 方差 = 均值 + 离散参数 × 均值²

而那个离散参数，DESeq2 会借用所有基因的信息一起估计。这种做法在统计里叫「信息共享」——对于低表达基因，单靠它自己那几个计数根本算不準离散参数，但借了高表达基因的力量就能稳住。

1.2 核心参数一：设计公式 (design formula)

设计公式是 DESeq2 里第一个必须认真对待的参数。它告诉模型“哪些因素会影响表达量”。最简单的写法：

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = colData,
                              design = ~ condition)

其中 condition 是你要比较的分组变量。如果你有批次效应，就得加上：

design = ~ batch + condition

这里要注意顺序：DESeq2 会用设计公式最后一列做对比。所以如果你写 ~ condition + batch，默认比较的是 condition 的水平差异，batch 只是协变量。我见过有人把顺序搞反了，结果比较的是批次差异，哭笑不得。

另外，如果你的因素有多个水平，比如三种处理（Ctrl, TreatA, TreatB），你可以用 results 的 contrast 参数指定要比较哪两个。我不会把所有对比都设成公式里的交互项，那样参数太多，小样本下容易过拟合。

1.3 核心参数二：差异检验方法

DESeq2 提供两种检验：Wald test 和 LRT（似然比检验）。Wald 是默认的，速度快，适合简单的两组比较。但如果你要比较多个条件，或者要评估某个因素的整体显著性（比如“批次这个因素到底有没有影响”），LRT 更合适。

# Wald test (默认)
dds <- DESeq(dds)
res <- results(dds)

# LRT：比较全模型 vs 简化模型
dds <- DESeq(dds, test = "LRT", reduced = ~ 1)
res <- results(dds)

除了检验方法，还有两个过滤参数我也经常调：

• cooksCutoff：用来标记异常样本。默认是 TRUE，DESeq2 会自动计算 Cook's distance 并过滤掉那些单个样本影响力过大的基因。如果你样本数太少（比如每个组只有2个），自动过滤可能会把好基因也杀掉。我一般样本数少于4时会把 cooksCutoff 改成 FALSE，然后手动检查。
• independentFiltering：默认开启，它会根据基因的平均表达量自动选择一个阈值，过滤掉那些表达太低、没法做差异检验的基因。这能减少多重检验负担。我个人觉得这个功能很好用，一般保持默认即可。

1.4 核心参数三：倍数变化收缩 (shrinkage)

这是 DESeq2 最让我喜欢的地方。原始的 log2 倍数变化（LFC）对于低表达基因常常高估得离谱。比如一个基因在对照组计数为2，处理组为10，直接算 LFC = log2(10/2) ≈ 2.32。但如果这个基因本身波动很大，这个 2.32 的效应很可能只是噪声。

DESeq2 通过 lfcShrink 函数做收缩（shrinkage），把那些方差大、证据弱的倍数变化往0缩。这个收缩不是随便减一个数，而是利用了先验分布——说白了就是“借力”，借所有基因的信息来告诉你，这个基因的效应不可靠，得缩回去。

res_shrink <- lfcShrink(dds, contrast = c("condition", "treated", "control"),
                        type = "apeglm")

常用的收缩方法有 apeglm 和 ashr。我个人更推荐 apeglm，它在大多数情况下性能稳定，而且收缩后的 LFC 仍然可以用标准阈值（比如 |LFC| > 1）来筛选。

1.5 核心流程与数据要求

整个 DESeq2 分析流程可以概括为三步：

1. 构建对象：DESeqDataSetFromMatrix
2. 差异分析：DESeq
3. 提取结果：results 和 lfcShrink

输入数据也很明确：

我通常从 featureCounts 或 HTSeq 得到计数矩阵，然后直接用 R 读入。注意计数矩阵必须全是整数，不要用 RPKM/FPKM 之类的值——DESeq2 的设计就是基于原始计数的，它自己会做归一化。

1.6 参数概览速查表

下表列出了本章涉及的核心参数，方便你快速回顾：

这份表格可以贴在墙上参考，但实际使用时，还是要根据你的数据特点来做选择。比如样本量很小，就不要轻易用自动的 cooksCutoff；如果是为了做功能富集，收缩后的 LFC 比原始 LFC 更可靠。