📚 课程目录第4章 DESeq2运行与结果提取

第4章：DESeq2运行与结果提取——从对象创建到差异表

上一章我们把表达矩阵和样本信息都准备好了。嗯，现在终于进入正题——跑DESeq2。

我个人觉得，DESeq2之所以受欢迎，不仅仅是算法好，关键是它的流程设计得特别“人性化”。你只需要创建好数据对象，剩下的事情一个函数搞定。但别高兴太早，创建对象这一步有个坑，我见过不少人在这里翻车。

4.1 创建DESeqDataSet对象

DESeq2的所有分析都围绕着一个核心对象——DESeqDataSet。说白了，它就是把三样东西绑一起：

原始计数矩阵（count matrix）
样本信息表（colData）
分析设计公式（design formula）

创建代码其实很简单，但细节决定成败：

library(DESeq2)

count_matrix <- read.table("counts.txt", header = TRUE, row.names = 1)
col_data <- read.table("metadata.txt", header = TRUE, row.names = 1)

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

这里面有几个地方我要特别提醒：

要注意的点：

count_matrix 必须是整数，不能是归一化后的值。DESeq2内部自己会做标准化，你传FPKM进去反而会出问题。
col_data的行名必须与count_matrix的列名一一对应。顺序可以不一样，但名字要对得上。
design公式里的变量名要和col_data里的列名完全一致。大小写也要注意。

我曾经帮一个同事看代码，他的一直报错，查了半天发现是列名里有个空格没去掉。嗯，这种事情真的很烦人，所以建议你们在预处理阶段就把列名统一成纯字母加下划线。

还有一个细节——因子水平顺序。DESeq2默认会把第一个水平当作对照。如果你想手动指定，可以这样做：

col_data$condition <- factor(col_data$condition, 
                             levels = c("control", "treatment"))
dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = col_data,
    design = ~ condition
)

我的习惯：我会在创建对象之前，先用colData(dds)确认一下样本信息对不对。多看一眼，省得后面浪费时间。

4.2 差异分析流程——一个函数搞定

对象创建好之后，核心步骤就是一个函数：

dds <- DESeq(dds)

没错，就这一行。DESeq()函数内部会自动执行以下几个步骤：

估计每个基因的size factor（相当于标准化因子）
估计基因的离散度（dispersion）
拟合负二项分布模型
进行Wald检验或LRT检验

你想想看，如果手动做这些，光离散度估计那块就得写不少代码。DESeq2把这些打包在一起，确实省了很多事。

不过要注意，这一步可能会跑几分钟，特别是基因数多的时候。我记得有一次跑了五万个基因，等了快十分钟。别着急，也别中途中断。

4.3 结果表提取

差异分析跑完后，提取结果用的是results()函数：

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

这里contrast参数指定你要比较哪两组。顺序是：先处理组，后对照组。这样算出来的log2FC就是正数表示处理组上调。

结果表里的列长这样：

列名	含义	备注
`baseMean`	所有样本的平均表达量	归一化后的值
`log2FoldChange`	处理组 vs 对照组的log2倍数变化	正值表示上调
`lfcSE`	log2FC的标准误差	用于计算p值
`stat`	Wald检验统计量
`pvalue`	原始p值	未经多重假设校正
`padj`	校正后的p值（BH方法）	我们通常用这个筛选差异基因

你会注意到，结果表里没有直接显示基因名。因为DESeq2默认把基因名放在行名里。如果你想把行名变成一列，方便后续分析，可以用：

res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df)

日志示例：运行results()后，控制台会输出类似这样的信息：

log2 fold change (MLE): condition treatment vs control 
Wald test p-value: condition treatment vs control 
DataFrame with 58438 rows and 6 columns

仔细看这些输出信息，它会告诉你比较的是哪两组，用了什么检验方法。

4.4 筛选差异表达基因

结果表拿到了，接下来就是筛选差异表达基因。常规标准是：

padj < 0.05：显著差异的阈值
|log2FC| > 1：倍数变化大于2倍（或小于0.5倍）

代码实现很直接：

deg <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(deg)  # 看看筛选出多少个基因

不过我个人建议，先别急着用太严格的阈值。你可以先看看padj小于0.05的有多少，再结合log2FC调整。如果差异基因太少，说不定是实验设计的问题。

4.5 结果表的导出与可视化准备

分析完了，结果要保存下来。常用的做法是：

write.csv(as.data.frame(res), file = "DESeq2_results.csv")
write.csv(as.data.frame(deg), file = "DESeq2_deg.csv")

这两个文件后面用来做火山图、热图什么的。嗯，这一步虽然简单，但我建议大家养成好习惯——保留完整结果表，别只保存筛选后的。万一后面想调整阈值，还得重新跑analysis。

我的做法：我会把完整结果表、差异基因表、以及归一化后的表达矩阵都导出来。这些后面做可视化的时候都用得上。

归一化表达矩阵可以用这个提取：

norm_counts <- counts(dds, normalized = TRUE)
write.csv(norm_counts, file = "normalized_counts.csv")

4.6 本章知识结构图

我把这一整套流程画成了图，你看了应该能理清脉络：

4.7 避坑指南

最后分享几个我踩过的坑：

我曾经遇到的问题：

运行DESeq()后报"error in estimateSizeFactors"——后来发现是矩阵里有NA值。DESeq2不接受缺失值。
结果表里有大量NA的p值——那是因为有些基因表达量太低，离散度估计不了。
padj值全是1——检查一下是不是对照组和处理组的样本名搞反了。

嗯，这一章的内容就到这里。你只要把这三个核心函数用好——DESeqDataSetFromMatrix、DESeq、results，整个分析流程就算拿下了。