Bulk RNA-seq 分析入门（二）：DESeq2 设计公式、对比与差异表达

上一篇把 FASTQ 处理到来源明确的计数或定量结果。这一篇往前走一步：如何从一个计数矩阵得到“比较了谁、考虑了什么因素、结果有多不确定”都能够说清楚的差异表达表？

本文只讨论 bulk RNA-seq，不是单细胞分析。本篇练习使用 Bioconductor airway 包附带的公开计数对象，不需要先下载原始 FASTQ。它包含四个气道平滑肌细胞来源的处理与未处理样本，适合演示配对关系。数据背景来自 airway 官方说明。这里学习统计流程，不将演示结果当成新的实验发现或医疗建议。

代码按 2026-10-06 核验的官方 API 组织，DESeq2 官方教程标注的包版本为 1.52.0。读者实际运行要记录自己的版本，不能只写“用了 DESeq2”。本机未执行下列 R 分析，因此本文不提供杜撰的差异基因数量、p 值或结果图片。

1. 先确认输入，而不是先调用 DESeq

最容易让整个分析失去意义的错误，是把样本标签与计数列错配。一个矩阵有六列，元数据也有六行，这只是尺寸相同，不是身份相同。必须先检查 ID 唯一、双方集合相同，再按 ID 重排，最后检查顺序完全一致。

第二个检查是输入单位。普通的计数矩阵入口要求非负、有限、整数的原始计数，不接受 TPM、FPKM、VST、rlog 或已经做过库大小缩放的数值。发现小数时应追问生成方法，不能为了通过构造函数而直接四舍五入。Salmon 估计计数的小数可以走 tximport 的专门入口，和 TPM 伪装计数不是同一回事。DESeq2 计数输入说明

第三个检查是元数据与实验设计。没有足够独立生物学重复、关键标签存疑或明显失败文库时，应暂停检验。已经知道的配对关系不能因为“不加这个变量也能跑”就忽略；batch 也不能因为 PCA 看起来有影响就随意追加，而应先追查真实实验记录以及它是否可与处理因素分别估计。

2. 设计公式是一张解释承诺，不是装饰

考虑三个原创的假设问题：独立来源的处理/对照比较可写 ~ condition；两组都分布在各批次且批次效应可估计时可写 ~ batch + condition；同一个来源各自提供处理与对照样本时可写 ~ subject + condition。这些公式回答的问题不一样，不能只以哪张火山图更漂亮为选择标准。

本篇的 airway 练习采用 ~ subject + condition，将来源之间的稳定差别和我们感兴趣的处理差别区分开。这个加性模型估计一个共同的处理系数，不估计来源与处理的交互；它不能证明各来源的真实效应相同，也不能涵盖所有重复测量问题。复杂时间序列、多层随机效应或处理交互需要另外设计模型；不要把两个条件的配对模板复制到所有实验上。DESeq2 官方配对与多因素分析入口

若所有对照都来自 B1、所有处理都来自 B2，那么 batch 与 condition 携带同一份分组信息。~ batch + condition 无法告诉我们变化究竟来自处理还是批次，这不是把 batch 变量删掉就解决了。设计矩阵的满秩检查可以发现不可分别估计的参数，但“满秩”也不保证混杂完全消失或样本量足够，只是必要的数学检查。

3. 完整练习：公开计数、配对模型与结果保存

运行前需要可用的 R 与相容的 Bioconductor 环境，并已经具备 airway、SummarizedExperiment、DESeq2、apeglm。本段检查依赖是否存在，不会自动安装。依赖缺失时先根据 Bioconductor 官方安装说明准备独立环境，确认安装权限、联网条件与版本兼容，再运行练习。

将下面代码保存为 airway_deseq2_learning.R，在自己的学习工作目录运行 Rscript airway_deseq2_learning.R。它只读取包自带的公开数据，在当前工作目录下创建一个带唯一名称的全新结果目录，不覆盖已有目录，也不会因 Rscript 退出而随会话临时目录消失。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
required <- c("airway", "SummarizedExperiment", "DESeq2", "apeglm")
available <- vapply(required, requireNamespace, logical(1), quietly = TRUE)
if (!all(available)) {
  stop("Prepare compatible R/Bioconductor packages first: ",
       paste(required[!available], collapse = ", "))
}
suppressPackageStartupMessages({
  library(airway)
  library(SummarizedExperiment)
  library(DESeq2)
})
data("airway", package = "airway")
cts <- assay(airway, "counts")
meta <- as.data.frame(colData(airway))

stopifnot(is.matrix(cts), is.numeric(cts), all(is.finite(cts)),
          all(cts >= 0), all(cts == floor(cts)),
          all(cts <= .Machine$integer.max),
          !anyDuplicated(rownames(cts)), !anyDuplicated(colnames(cts)),
          !anyDuplicated(rownames(meta)),
          setequal(colnames(cts), rownames(meta)))
meta <- meta[colnames(cts), , drop = FALSE]
stopifnot(identical(colnames(cts), rownames(meta)))
meta$subject <- factor(meta$cell)
meta$condition <- factor(meta$dex,
                         levels = c("untrt", "trt"),
                         labels = c("control", "treated"))
stopifnot(!anyNA(meta$subject), !anyNA(meta$condition),
          nlevels(meta$subject) >= 3L,
          all(table(meta$subject, meta$condition) == 1L))

design_formula <- ~ subject + condition
design_matrix <- model.matrix(design_formula, data = meta)
if (qr(design_matrix)$rank != ncol(design_matrix)) {
  stop("Design is not full rank: review confounding and unused levels")
}
if (nrow(design_matrix) <= ncol(design_matrix)) {
  stop("No residual degrees of freedom for this teaching model")
}
storage.mode(cts) <- "integer"
dds <- DESeqDataSetFromMatrix(countData = cts, colData = meta,
                             design = design_formula)

# A declared, condition-independent teaching prefilter.
smallest_group <- min(table(meta$condition))
keep <- rowSums(counts(dds) >= 10L) >= smallest_group
stopifnot(any(keep))
dds <- dds[keep, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05)
coef_name <- "condition_treated_vs_control"
stopifnot(coef_name %in% resultsNames(dds))
res_shrunk <- lfcShrink(dds, coef = coef_name, type = "apeglm")

tab <- data.frame(gene_id = rownames(res), as.data.frame(res),
                  row.names = NULL, check.names = FALSE)
tab$log2FC_shrunk <- res_shrunk$log2FoldChange[
  match(tab$gene_id, rownames(res_shrunk))
]
tab$has_test_pvalue <- is.finite(tab$pvalue)
tab$has_adjusted_pvalue <- is.finite(tab$padj)

out <- tempfile("perry-airway-de-", tmpdir = getwd())
stopifnot(dir.create(out))
saveRDS(cts, file.path(out, "input_raw_counts.rds"))
write.table(meta, file.path(out, "sample_metadata.tsv"),
             sep = "\t", quote = FALSE, col.names = NA)
write.table(tab[order(tab$padj, na.last = TRUE), ],
             file.path(out, "treated_vs_control.tsv"),
             sep = "\t", row.names = FALSE, quote = FALSE)
saveRDS(list(dds = dds, res = res, res_shrunk = res_shrunk,
             metadata = meta, design = design_formula,
             contrast = c("condition", "treated", "control"),
             input = "Public airway teaching dataset"),
         file.path(out, "airway_bundle.rds"))
capture.output(sessionInfo(), file = file.path(out, "sessionInfo.txt"))
capture.output(resultsNames(dds), file = file.path(out, "coefficients.txt"))
message("Teaching outputs saved to: ", normalizePath(out))
message("Keep the printed bundle directory path for lesson 3.")

这里的预过滤规则是在最小条件组样本数那么多的样本中，计数至少为 10。它是为了教学声明的一种规则，不是所有物种、测序深度和实验的唯一正确选择，也不是先看对比结果再把不喜欢的基因删掉。保存的原始计数没有因过滤而消失，结果表则对应进入模型的基因。

运行结束后，结果目录保留在你启动 Rscript 时的当前学习工作目录里。保存打印出的完整目录路径，其中的 airway_bundle.rds 将作为下一篇的输入；不要移动或重命名后还使用旧路径。这里 tempfile() 只是生成一个未被占用的目录名称，实际位置明确指定为 getwd()，并非 R 会话自动清理的临时目录。

4. 如何读一行结果：效应与证据要分开

baseMean 描述归一化计数的平均水平；log2FoldChange 是该对比的效应估计；lfcSE 表达其标准误；stat、pvalue 和 padj 是相应统计检验与多重比较处理的产物。这里的方向固定为 treated 相对 control，正值表示模型下处理组更高，负值表示更低。将分子分母互换，效应符号随之改变，所以导出的表必须附有对比说明。DESeq2 结果字段手册

如果一个基因估计 log2FC 为 1，表示该模型中的相对倍数约为 2；这只是算术含义，不意味着这个基因一定显著，更不意味着它是整个表里最重要的基因。另一个基因可以有很小的效应，却因为信息更充分而得到较小 p 值。解释时应同时看效应、精度、表达水平和研究问题，而不是把单个数字当作综合评分。

padj < 0.05 是在特定模型、数据和假设下控制多重检验错误的一种判定规则，不是“这个基因有 95% 的概率是真的”，也不是任意一个选中基因有恰好 5% 假阳性概率。阈值之外的基因也不自动等于无变化。BH 等多重比较方法的定义可参阅 R 官方 p.adjust 文档。

5. shrink、筛选与效应阈值检验不是同一件事

当一个基因信息少或变异较大时，未经收缩的倍数估计可能很极端。lfcShrink 提供更稳定的效应估计，常用于可视化和效应排序。这不等于重新定义了本篇普通 Wald 检验的原假设，也不是把所有小效应变成“没有生物学意义”。本练习明确保留原检验的 pvalue/padj，另列 log2FC_shrunk，避免读者不知道某一列来自哪一步。收缩的模型依据可阅读 DESeq2 原始方法论文。

“普通检验 padj 小于阈值，并且收缩后绝对 log2FC 大于 1”可以作为声明过的候选展示规则，却不等价于已经检验“真实绝对效应大于 1”。若问题确实是超过一个预先指定的最小效应，应使用相应的效应阈值检验，例如：

1
2
3
4
5
6
7
8
9
# 在完整练习后继续运行：检验 |log2FC| 是否超过预定阈值1。
# 它产生新的检验结果，不能混称为上面的普通差异检验。
res_min_effect <- results(
  dds,
  contrast = c("condition", "treated", "control"),
  lfcThreshold = 1,
  altHypothesis = "greaterAbs",
  alpha = 0.05
)

默认 DESeq() 足够用于这个常规 bulk 学习入口。不要看到某篇单细胞性能教程，就盲目切换到 glmGamPoi 或同时改多项参数。任何拟合方式、离群值处理和阈值变更，都应说明对应的数据特征与验证依据，而不只是追求运行更快或基因更多。

6. Salmon 输入的替代入口：不把 TPM 当成计数

如果上一篇实际采用 Salmon 定量，建议从每样本 quant.sf 导入，而不是手动抽出 TPM 交给 DESeqDataSetFromMatrix()。下面是有真实文件之后的替代入口模板，不需要和 airway 练习同时运行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# metadata.tsv: first column sample_id; other columns include condition.
# tx2gene.tsv: two named columns, transcript_id and gene_id, from the same reference.
library(tximport)
library(DESeq2)
meta <- read.delim("metadata.tsv", row.names = 1, check.names = FALSE)
stopifnot(!anyDuplicated(rownames(meta)))
meta$condition <- factor(meta$condition, levels = c("control", "treated"))
stopifnot(!anyNA(meta$condition))
files <- setNames(file.path("salmon", rownames(meta), "quant.sf"), rownames(meta))
stopifnot(all(file.exists(files)))
tx2gene <- read.delim("tx2gene.tsv", check.names = FALSE)
stopifnot(identical(names(tx2gene), c("transcript_id", "gene_id")),
          !anyNA(tx2gene), !anyDuplicated(tx2gene$transcript_id))
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                 countsFromAbundance = "no")
stopifnot(identical(colnames(txi$counts), rownames(meta)))
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ condition)

这条路线通过完整对象向 DESeq2 传递估计计数与长度相关信息；不要又手动对矩阵做一次长度缩放，造成重复校正。countsFromAbundance 的不同选项有不同下游约定，本文只展示 no 加专门构造器的一条路线。转录本映射缺失、版本号不匹配、一个转录本映到多个基因等必须调查并记录；模板的结构校验不能替代参考来源审计。tximport 官方方法说明

~ condition 仅用于已经确认的简单独立设计。真实配对或批次实验应携带相应字段、检查设计并替换公式，不能因为换了导入入口就忘了配对关系。3′ 标签型 RNA-seq 的定量和长度处理也可能不同，应查相应协议，而不是把本模板视为所有 RNA 测序的标准答案。

7. 看到 NA 和报错，先判断原因

结果出现 NA，不应直接全部改为 1 或删除。全零、异常计数或独立过滤可能造成不同字段不可用，含义不一样。需要分别记录是否拥有检验 p 值和调整后 p 值，再按定义决定候选集或富集背景。报出非满秩时，优先检查完全混杂、没有样本的因子水平和冗余协变量，而不是不断删变量直到成功。

如果两组没有明确分开，也不要立即重做模型。样本差异可能主要由来源、细胞组成、时间或其他真实因素驱动。先回到原始 QC 与元数据，再解释探索图。计数模型加入 batch 与在展示图上去除 batch 是不同的事情；后一种视觉处理不能创造前一种设计里缺失的信息。

8. 本篇的自检与交付

在教学副本中，读入元数据后、按 ID 对齐之前插入 meta <- meta[rev(rownames(meta)), , drop = FALSE]。保留随后按 ID 对齐的那行时，顺序检查应通过，因为程序已经正确恢复匹配；再注释掉 meta <- meta[colnames(cts), , drop = FALSE]，顺序检查才应阻止继续。这是在验证对齐逻辑，不是在修改真实分组。把 condition 的参考水平反过来，检查 resultsNames() 与你指定的对比方向有何区别。构造全部对照在 B1、处理在 B2 的元数据，检查设计矩阵秩，练习解释为什么这是设计问题。

正式交付至少应有输入来源、过滤规则、模型公式、参考水平、对比清单、完整结果表、NA 的处理说明、实际会话版本与 QC 记录。保留整张表，而不只保留显著基因；未来改变阈值或重查某个候选时，才不必依赖截图猜测。

上一章：实验设计、FASTQ 与计数；下一章：QC 图、功能富集与结果解释。

官方与原始方法参考