第3章：DESeq2基本流程——构建dds对象、执行差异分析

说实话，我第一次用DESeq2的时候，被它的一长串函数名搞得有点懵。但后来你会发现，核心流程其实就三步：装数据、跑一步、看结果。嗯，今天我们就来拆解这三步里最关键的——构建dds对象和执行差异分析。我手头有个水稻转录组项目，正好拿那套数据来举例。

3.1 准备工作：数据长什么样

在调用任何DESeq2函数之前，你得确保手里的表达矩阵是“干净”的。说白了，就是原始整数counts矩阵，不能是RPKM/FPKM，不能是小数，更不能有负值——我见过有人把log2标准化后的数据塞进去，结果DESeq2直接报错。

矩阵长这样：

• 行：基因ID或symbol（唯一标识）
• 列：样本名（要与样品信息表的行名对应）
• 每个格子：非负整数，代表比对到该基因的read数

同时你还需要一个样品信息表（colData），至少包含一列分组信息（比如condition），以及可能的技术重复批次等。我个人习惯在构建dds之前，先把colData的行名设成样本名，并确保顺序和counts矩阵的列顺序一致——这一步踩过坑，顺序不对跑出的结果全是NaN。

# 检查顺序是否一致
all(colnames(counts_matrix) == rownames(colData))  
# 返回TRUE才放心

3.2 构建dds对象：那句核心代码

准备好了数据，就该干正活了。DESeq2用一个叫DESeqDataSet的对象来装所有数据和分析结果。构建它的函数是DESeqDataSetFromMatrix。你要是看过别的教程可能会用DESeqDataSetFromHTSeq，那只是读取格式不同，核心逻辑一样。

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

参数你别写反了：countData是表达式矩阵，colData是样品信息。design公式里的condition是colData里的列名。为什么要用波浪号？这是R公式语法，告诉DESeq2你要比较哪些组。

这里有个小细节——因子的水平顺序。DESeq2默认把第一个水平作为对照组。比如你的condition有“control”和“treatment”，我建议用factor把control放在前面：

colData$condition <- factor(colData$condition, levels = c("control","treatment"))

不这样做的话，比较方向可能反过来，log2FoldChange的符号就错了。我记得有次帮同事排查，折腾了一小时才发现是因子顺序的问题。

3.3 执行差异分析：一行命令就够

对象建好了，跑差异分析只需要一行：

dds <- DESeq(dds)

别小看这行代码，它后台干了四件事：

1. 估计size factors（校正测序深度差异）
2. 估计基因的离散度（dispersion）
3. 拟合广义线性模型（负二项分布）
4. 进行Wald检验（或LRT，默认是Wald）

你会看到终端刷出一堆信息，比如“estimating size factors”、“estimating dispersions”等。如果数据里有outlier，它还会提示“gene-wise dispersion estimates”。这些都是正常的。如果没看到任何警告……嗯，那可能是数据太完美了，现实中不常见。

3.4 提取结果：results() 的几种用法

跑完DESeq()，差异分析结果其实已经存在dds里了。你需要用results()把它拿出来。最简单的是：

res <- results(dds, contrast = c("condition","treatment","control"))

contrast参数指定要比较哪两组：第三个元素是实验组，第四个是对照组。你也可以用name参数，但掌握contrast就够用了。

返回的结果是一个DataFrame，关键列有：

通常会根据padj < 0.05 和 |log2FoldChange| > 1 来筛选差异基因。我习惯加上一个“baseMean”的过滤，比如 res$baseMean > 10，避免那些低表达基因带来的噪音。

3.5 流程总览：一张图看懂

下面这张图把刚才讲的几个核心步骤串了起来。你对照着图，再看看代码，应该就清楚整个管线了。

3.6 避坑指南——我当年踩过的雷

避坑1：非整数counts
我曾经从某个比对工具直接导出RPKM矩阵，习惯性地用round取整后塞进DESeqDataSetFromMatrix，结果报“negative values in counts matrix”的错。后来才发现round会把零变成负数？不，其实是因为RPKM里有过小的值，取整后有些变成负数（负的极小值）。正确的做法是——从源头就保存raw counts。

避坑2：design公式写错
如果你有多个批次，design公式里要包含批次效应，比如~ batch + condition。但别把批次的级别写成数字（比如1,2,3），否则R会当成连续变量。我那时候batch设成1-4的数字，结果模型给你拟合出一条斜线……

避坑3：过滤参数
results()默认帮我们做了独立过滤（independent filtering），目的是减少低表达基因的假阳性。但如果你设了alpha=0.01，过滤阈值也会跟着变。我建议先不改alpha，看一遍结果，再根据padj分布调整。

3.7 关键函数速查表