4. DESeq2简介：为什么选择DESeq2做批次校正？

当你面对一堆RNA-seq数据，发现不同批次、不同测序时间、不同操作人……结果差异比生物差异还大时，你的第一反应是什么？

我当年第一次遇到这种情况，差点把数据直接扔了。后来才知道，批次效应是RNA-seq里最坑人的陷阱之一。而DESeq2，就是专门用来解决这个问题的利器。

4.1 什么是DESeq2？

DESeq2是一个基于负二项分布的差异表达分析R包，由Love、Huber和Anders在2014年推出。它不仅能做差异分析，还能内建批次校正。

说白了，它把你的表达矩阵、分组信息和批次信息一起喂进去，模型自动把批次效应“剥离”掉，然后告诉你哪些基因在实验条件下真的变了。

核心公式长这样（你不需要背，但理解原理有用）：

K_ij ~ NB(mu_ij, alpha_i)
mu_ij = s_j * q_ij
log2(q_ij) = beta_i0 + beta_i_condition + beta_i_batch + ...

K_ij是观测到的counts，NB是负二项分布，s_j是文库大小因子，q_ij是归一化后的表达水平。最关键的是最后一行：感兴趣的实验条件（condition）和批次（batch）都作为线性项加在log2尺度上。这样一来，模型就能把batch效应单独估计出来，不影响condition的差异比较。

4.2 为什么非得用DESeq2做批次校正？

说实话，市面上能做批次校正的工具不少，比如ComBat、limma的removeBatchEffect、甚至是简单的线性回归。但我个人更偏爱DESeq2，原因有三：

1. 模型原生支持批次变量：你不用先校正再分析，所有步骤在一个模型里完成。标准误、p值、统计推断都考虑了批次不确定性。
2. 处理低表达基因更稳健：RNA-seq有很多低表达基因，counts离散程度高。负二项模型比正态模型更贴合实际。我在项目中遇到过用ComBat校正后低表达基因方差被严重低估的情况，换了DESeq2就好了。
3. 自带“shrinkage”（压缩估计）：对于小样本、高方差基因，DESeq2会把极端log2 fold change往0方向拉一下，避免假阳性。这个在批次校正后的差异分析中特别有用。

4.3 与其他方法的对比

你可能想问：ComBat不是专门做批次校正的吗？为什么还要用DESeq2？我给你列个表，一看就明白。

方法	统计模型	是否保留原始counts	适用场景	我的经验
ComBat	基于经验贝叶斯的正态模型	否（输出校正后表达值）	芯片数据、表达谱校正	对RNA-seq低表达基因效果不稳定
limma + removeBatchEffect	线性模型（正态假设）	否（输出校正后log值）	标准化后的表达数据	常用于芯片，counts需先vst变换
DESeq2	负二项GLM	是（模型内部处理）	RNA-seq count数据	我最常选，一步到位
edgeR	负二项似然比检验	是	类似DESeq2	支持批次但不如DESeq2灵活

看到了吗？DESeq2是唯一一个在counts层面、用负二项模型、批次作为协变量同时分析的。其他工具要么要你先变换数据，要么统计假设不匹配RNA-seq的离散特性。

4.4 批次校正的流程概览

咱们用DESeq2做批次校正的完整流程，我习惯分成5步。下面这张SVG图帮你把逻辑理清楚，后面每一章都会对应到这些步骤。

4.5 什么时候不该用DESeq2做批次校正？

当然，DESeq2不是万能的。我遇到过几种情况，它并不是最佳选择：

• 批次和条件完全混杂：比如所有对照组是第一批次，所有处理组是第二批次。模型无法区分到底是批次效应还是处理效应。这时候你得重新设计实验，没有任何统计方法能救你。
• 只有2个样本但批次不同：DESeq2在估计批次效应时需要每个批次有至少2个观测值。如果某个批次只有1个样本，模型会报错或给出不可靠估计。
• 需要可视化校正后的“干净”表达值：DESeq2的批次校正在模型内部。如果你想拿到去除批次效应后的表达矩阵做聚类或热图，用DESeq2的vst()或rlog()配合blind=FALSE会更合适。

4.6 一个最简单的例子先热热身

别急，咱们先不深入，先看最小可运行代码是什么样子。假设你已经有了count矩阵和样本信息表：

library(DESeq2)

# 假设你的数据
counts <- as.matrix(read.csv("counts.csv", row.names=1))
colData <- read.csv("metadata.csv", row.names=1)

# 检查批次和条件都是factor
colData$batch <- as.factor(colData$batch)
colData$condition <- as.factor(colData$condition)

# 建对象，design里batch在前
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = colData,
                              design = ~ batch + condition)

# 核心分析（这一步搞定了所有估计和检验）
dds <- DESeq(dds)

# 提取处理组vs对照组的差异结果
res <- results(dds, contrast = c("condition", "treated", "control"))
resOrdered <- res[order(res$padj),]
head(resOrdered)

看到了吗？就这几行。DESeq()函数内部自动估计了离散度、拟合了包含批次项的GLM、计算了p值。你不需要手动做任何“批次校正”步骤。