02 DESeq2 差异表达分析

DESeq2 是 bulk RNA-seq 差异分析里最常用的工具。它假设 counts 服从负二项分布，用每个基因的表达均值和离散度估计显著性，加上 Benjamini-Hochberg 校正给出最终的 padj。流程非常规整，稳到几乎是默认选择。

这一章用 DESeq2 官方教程数据 airway 走一遍完整流程。数据是 8 个人类气道上皮细胞样本 —— 4 个对照 + 4 个糖皮质激素（地塞米松，dex）处理，来自 Himes et al. 2014（PLoS ONE 9(6): e99625）。

核心流程是三步

DESeq2 的调用界面非常简洁。真正的 API 就三步：

dds <- DESeqDataSet(se, design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "trt", "untrt"))

1. 构建 DESeqDataSet 对象，声明设计公式（design 最后一列是要检验的变量）
2. DESeq() 一个调用里做完归一化、离散度估计、Wald 检验
3. results() 提取指定对比的差异结果

真实项目的工作重点不在这三行代码上，而在围绕它的 QC 和解释：样本是否正常、模型是否合适、结果是否可信。下面按这个顺序走。

加载 airway 并构建 DESeq 对象

library(DESeq2)
library(airway)

data(airway)
se <- airway
# 把 untrt 设为参考水平，这样 log2FC 是 "trt 相对 untrt"
se$dex <- relevel(se$dex, ref = "untrt")

# design: 控制细胞系差异，检验 dex 处理效应
dds <- DESeqDataSet(se, design = ~ cell + dex)

# 预过滤：counts 总和 < 10 的基因直接丢掉
# 不影响统计正确性，但能加速并减少内存
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

airway 的 colData 里 cell 是细胞系（4 个不同来源），dex 是处理变量。design = ~ cell + dex 的意思是"在控制细胞系差异之后看 dex 的影响"。设计公式里变量的顺序无所谓，但要检验的变量一般放最后，results() 默认会拿最后一个变量的第一个对比。

跑 DESeq 并看总体结果

dds <- DESeq(dds)

res <- results(dds, contrast = c("dex", "trt", "untrt"))
summary(res)

summary(res) 会告诉你上下调基因各有多少个（默认阈值 padj < 0.1）。对 airway 这份数据，通常会看到几千个显著差异基因 —— 糖皮质激素是强效应通路，差异很明显。

真实示例：airway 上的完整 DESeq2 分析

配套脚本 bulk02_deseq2_sci.R 把上面的流程完整跑了一遍，输出 6 张真实数据图和一份 DE 基因表。装好 DESeq2 + airway + apeglm 就能直接跑，不需要下载 FASTQ：

Rscript scripts/bulk02_deseq2_sci.R

脚本里做的事：预过滤 → DESeq() → results() → lfcShrink(apeglm) → 6 张图 → 输出 DE 表。下面逐张看。

图 1：每个样本的库大小

横轴是样本 ID，纵轴是预过滤之后每个样本的总 reads（百万）。这张图是最基础的 QC：库大小应该在一个量级，如果有某个样本只有其他样本的 1/5，就要回去看测序报告。airway 里 8 个样本大小相当，没问题。

图 2：VST 归一化后的 PCA

vst(dds, blind = FALSE) 做方差稳定变换 —— 让高表达和低表达基因的方差都差不多，这样 PCA 才不被少数高表达基因主导。颜色是 dex 处理，形状是细胞系。

这张图回答两个问题：

• PC1 是什么：理想情况下 PC1 直接把处理组和对照组分开。airway 里 PC1 正是 dex 轴。
• 细胞系差异大不大：同样颜色但不同形状的点会偏开，说明细胞系间有个体差异。这就是 design 里加 cell 的原因。

图 3：样本间欧式距离

样本两两之间的欧式距离（在 VST 空间）画成热图，再聚一次类。同一条件 + 同一细胞系的样本应该靠在一起。如果某个样本跟同组的样本距离反而大过别组，要小心标签弄错或者该样本有问题。

图 4：离散度估计

DESeq2 的核心是"用所有基因一起拟合 mean-dispersion 关系"，低表达基因数据少、方差估计不稳，就借用趋势回到那条红色的拟合曲线上。灰色是每个基因的原始离散度，红线是拟合趋势。

看这张图主要是确认：趋势是不是随均值单调下降。如果红线走势很乱，可能说明样本里有严重的技术偏差。正常的曲线应该像这张一样，左下右上缓慢下降。

图 5：p 值分布与 MA 图

左边是原始 p 值直方图。一张健康的直方图长这样：00.05 一段突出（真正显著的基因），0.051 平坦分布（无效假设下的均匀分布）。如果直方图两端高中间低，说明离散度估计有问题；如果整体偏左（p 值都很小），可能是 design 漏了重要的协变量。

右边是 MA plot，横轴 log10(均值表达量)，纵轴 log2 fold change。红点是显著 DE 基因（|LFC| > 1 且 padj < 0.05）。这里用的是 apeglm shrunken LFC，低表达基因（左侧）的 fold change 被往 0 拉，避免"表达量低 + 噪声大 → fold change 虚高"的假象。

图 6：Top DE 基因的归一化 counts

按 padj 排最前面的 6 个基因的 normalized counts，每条细胞系画成一个颜色。y 轴是 log10。每个基因都能看到 4 个 untrt 和 4 个 trt 点，处理组和对照组在 y 轴上分得很开，且不同细胞系的点虽然绝对值有差异，但处理效应方向一致 —— 这是 DESeq2 在控制了 cell 之后仍然能拿到强显著性的原因。

提取结果表

脚本末尾把 shrunken log2FC + 原始 padj 写成一张 TSV：

de_tbl <- as.data.frame(res_shrink) |>
  tibble::rownames_to_column("gene_id") |>
  mutate(padj_raw = results(dds, contrast = c("dex", "trt", "untrt"))$padj) |>
  arrange(padj_raw)

write.table(de_tbl, "de_genes_airway.tsv", sep = "\t", quote = FALSE, row.names = FALSE)

这张表的列：gene_id、baseMean、log2FoldChange（shrunken）、lfcSE、pvalue、padj_raw。进 clusterProfiler 或 fgsea 做富集分析就用这张表。

套到自己数据上

脚本的前 15 行把 airway 换成自己的数据就行。常见两种输入：

# 输入 1：tximport 的结果
library(tximport)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + condition)

# 输入 2：featureCounts 的 counts 矩阵
count_mat <- read.delim("counts.tsv", row.names = 1)
samples <- read.delim("metadata/samples.tsv", row.names = 1)
dds <- DESeqDataSetFromMatrix(count_mat, colData = samples, design = ~ batch + condition)

其余的 DESeq() / results() / lfcShrink() 调用不用改。真实项目里最容易踩的坑：

• 样本表里 rownames 要和 counts 矩阵的列名严格一致（包括顺序）
• 如果样本有批次（不同测序 run、不同试剂盒批次），一定要加进 design
• 小于 3 对 3 的实验 DESeq2 能跑，但统计力很弱，不建议做复杂 design

下载资源

bulk02_deseq2_sci.R12 KB

↓

下载 airway DESeq2 完整脚本

下一步

• 03 功能富集：GO / KEGG / GSEA
• 04 火山图、热图与富集可视化

参考资源

• DESeq2 官方教程
• RNA-seq workflow (Bioconductor)
• apeglm 方法
• Himes et al. 2014, airway 数据源