📚 课程目录第3章设计公式设计

3. 设计公式设计：因子、交互项与批次效应处理

到了这一章，我们终于要直面 DESeq2 最核心的环节——设计公式。说实话，很多用户在跑差异分析时结果不对劲，十有八九是公式写错了。我自己就踩过这个坑，而且不止一次。

设计公式说白了就是告诉 DESeq2：“你要怎么解释我的表达量数据？”它是你实验设计的数学映射。因子、交互项、批次效应，这三块要是没处理好，后续的一切分析都可能是泡沫。

  核心观点： 设计公式决定了你的差异比较能否反映真实生物学差异，还是被混淆变量污染。

3.1 因子类型与编码

DESeq2 里的因子分三种：字符型、数值型、有序型。但我真正关心的只有两种——离散因子和连续因子。离散因子比如“处理组 vs 对照组”，连续因子比如“年龄”、“测序深度”。

R 语言默认会把字符型自动转为因子，但顺序是按字母排的。你想想看，如果对照组是“Control”，处理组是“Treatment”，字母顺序是 Control > Treatment，那对照组就成了基线。这不一定是你想要的。

因子类型	R 中存储方式	DESeq2 默认处理	我个人的建议
离散因子（如分组）	`factor()`	取第一个水平为基线	用 `relevel()` 手动指定基线
连续因子（如时间）	`numeric`	作为协变量，不做对比	确保没有缺失值，且量纲合理
有序因子（如剂量低/中/高）	`ordered()`	多项式对比，慎用	除非明确知道趋势，否则别乱用

一个小技巧： 在构建 dds 对象前，先用 colData(dds)$condition <- relevel(colData(dds)$condition, ref = "control") 明确基线。这个习惯能让你避免很多莫名其妙的比较结果。

3.2 交互项：什么时候加，什么时候不加

交互项处理的是“一个因子的效应是否依赖于另一个因子”。比如，药物 A 在女性中有效，在男性中无效，这就是交互效应。如果你只写主效应模型 ~ sex + treatment，那你会错过这个关键信息。

但交互项不是免费的午餐。每加一个交互项，参数数量就翻倍，自由度消耗也翻倍。样本量小的项目，交互项会让模型不稳。

我记得有一次帮一个合作伙伴分析数据，他们有 3 个时间点 × 2 种处理，每组只有 3 个生物重复。他们非要把交互项 ~ time * treatment 写进去，结果 Wald 检验出来的 p 值几乎全是 1。为什么？因为交互项的参数太多，模型几乎饱和了。后来我建议他们换用 LRT 检验，只比较 ~ time + treatment 和 ~ time * treatment，才找到几个真正的交互基因。

# 主效应模型
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData = coldata,
                               design = ~ time + treatment)

# 交互项模型
dds_interact <- DESeqDataSetFromMatrix(countData = counts,
                                        colData = coldata,
                                        design = ~ time + treatment + time:treatment)
# 或者简写为
dds_interact <- DESeqDataSetFromMatrix(countData = counts,
                                        colData = coldata,
                                        design = ~ time * treatment)

避坑指南： 我曾经见过一位用户在交互项模型里做 resultsNames()，发现有一堆 "time2:treatmentB" 之类的名称，然后不知道怎么提取结果。记住：交互项的系数名是 因子1水平:因子2水平，用 contrast 参数指定时一定要匹配这些名字。

3.3 批次效应：你的数据里可能藏着一个“幽灵因子”

批次效应是 RNA-seq 数据里最烦人的问题，没有之一。它可能是测序仪不同、试剂批次不同、甚至操作人员不同造成的系统性偏差。如果你不把它写进设计公式，它就会混杂到生物学差异里。

我个人习惯：只要样本不是同一天、同一批处理完的，一定把批次信息放到设计公式里。哪怕你觉得“应该没影响”，也放进去。因为如果你放对了，它只消耗 1 个自由度；如果你没放，它可能虚假产生几千个差异基因。

比如你的实验有两个批次，设计公式应该是：

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

这样 DESeq2 会先把批次效应从数据中“扣除”，然后再比较 condition 之间的差异。注意，是扣除，不是消除。批次效应仍然存在，但它不会干扰你的组间比较。

  一个重要判断： 如果你的批次和实验条件是完美共线的——比如批次1全部是对照组，批次2全部是处理组——那你就不能把批次放进设计公式，因为模型无法区分“批次差异”和“处理效应”。这种情况，要么重新设计实验，要么用 sva 或 RUVseq 先估计隐藏因子。

举个极端例子：我有一次帮人看数据，批次1有3个对照+3个处理，批次2有3个对照+3个处理。看起来完美平衡，对吧？但批次1的测序深度比批次2高20%。如果不放批次因子，差异基因里会出现大量核糖体蛋白基因——因为它们的表达量和测序深度强相关。放了批次因子后，这些假阳性就消失了。

3.4 设计公式的“可辨识性”陷阱

这是新手最容易踩的坑：你把所有可能的因子都放进公式，然后 DESeq2 报错说“design matrix is not full rank”。什么意思？就是你给的因子之间有线性依赖关系，比如你同时放了 ~ condition + batch + condition:batch，但条件 × 批次的交互项和两个主效应是线性相关的。

DESeq2 的解决方式是自动丢弃那些冗余的列，但你能看到的系数名就变得莫名其妙。我教你一个检查方法：

# 构建模型矩阵
model.matrix(~ condition + batch, data = colData(dds))

# 查看秩
library(Matrix)
rankMatrix(model.matrix(~ condition + batch, data = colData(dds)))

如果秩小于列数，说明有共线性。一般我的做法是：从最简单的主效应模型开始，确认能跑通，再逐步加交互项和批次因子。每加一个，检查一下 resultsNames 是不是可解释的。

我的个人经验： 当设计公式复杂时，优先使用 LRT 检验而不是 Wald 检验。因为 LRT 比较的是全模型和简化模型，它对系数之间的相关性更鲁棒。而且 LRT 可以测试整个交互项的显著性，而不仅仅是某个系数。

最后，我用一张图总结设计公式的决策流。记住核心：因子要完整、交互要谨慎、批次要纳入、模型要可辨识。

看到这个流程图了吗？核心逻辑就是一层层过滤。从实验因子出发，先处理批次，再决定要不要交互项。千万别一上来就写个大公式，除非你确定数据能撑住那么多参数。

好了，关于设计公式的因子、交互项和批次处理，我把自己踩过的坑和总结的经验都摊出来了。下次写设计公式时，多想想这几点，你的差异结果会干净很多。