📚 课程目录第2章数据输入与DESeqDataSet构建

2. 数据输入与 DESeqDataSet 构建：countData、colData 与 design 公式

讲差异分析，第一步就是喂数据。

很多人觉得这步简单，不就是读个表格吗？
我在项目里见过太多次，跑完才发现样本名对不上、分组写错了、design 公式漏了关键变量——全得重来。

今天我们就掰开揉碎，把 countData、colData、design 公式讲透。

核心一句话：DESeqDataSet 是 DESeq2 一切计算的起点。它把计数矩阵、样本元数据和统计模型绑在一起。这三样东西，一个都不能错。

2.1 三大输入要素

参数	数据类型	作用
`countData`	整数矩阵 / data.frame	存放每个基因在每个样本中的原始 read count
`colData`	data.frame	描述每个样本的实验条件、批次、协变量等
`design`	公式 (formula)	告诉 DESeq2 你的实验设计，比如 ~ condition 或 ~ batch + condition

说白了，countData 是「谁有多少」，colData 是「谁是谁」，design 是「我要比什么」。

2.2 countData：你的计数矩阵

countData 长什么样？我给你一个最小示例：

# 假设我们只有 4 个基因、3 个样本
countData <- matrix(
  c(100, 85, 92,
    250, 233, 245,
    45,  52,  48,
    780, 755, 801),
  nrow = 4, byrow = TRUE
)
rownames(countData) <- c("GeneA", "GeneB", "GeneC", "GeneD")
colnames(countData) <- c("sample1", "sample2", "sample3")
countData

几点硬规则：

必须是整数——不是 TPM、FPKM，也不是 log 转换后的值。原始 count，一颗老鼠屎都不能少。
行是基因，列是样本。行列名必须和 colData 的样本一一对应。
不允许有 NA、负值。遇到就报错，我没见过例外。

警告：千万不要把标准化后的表达量扔进 DESeq2。标准化是 DESeq2 内部自己算的，你提前归一化会影响它内置的 size factor 估计。我见过有人直接把 TPM 丢进去——出来的 p-value 全不对，根本没法用。

2.3 colData：样本的「身份证」

colData 告诉你每个样本属于哪个组、哪个批次、什么条件。

colData <- data.frame(
  row.names = c("sample1", "sample2", "sample3"),
  condition = c("control", "treated", "treated"),
  batch     = c("A", "A", "B")
)
colData

这里的关键：行名必须和 countData 的列名完全一致。顺序可以不一样，DESeq2 会自动对齐。但名字不能写错也别多写少写。

我有个经验：colData 的第一列尽量用 factor，不要用 character。你想想看，character 是按字母序的，control 和 treated 在 R 里会被当成两种字符串，factor 才能指定谁是基线水平。

关于基线水平：DESeq2 默认按字母序选对照组。如果你的对照组是 "control"，处理组是 "treatment"，那没问题。但如果是 "treated" 和 "untreated"——untreated 比 treated 字母更靠前，DESeq2 会拿 untreated 当基线。不放心的话手动设置：

colData$condition <- relevel(colData$condition, ref = "control")

小技巧：我喜欢在项目一开始就把 colData 写好，写成一个 CSV 文件。每次跑差异分析前读进来，改改 condition 列就完事。这样不容易出错，也方便别人复现。

2.4 design 公式：告诉模型你要比什么

design 是 R 的公式对象，格式固定：

design <- ~ condition
design <- ~ batch + condition
design <- ~ condition + time + condition:time   # 交互项

最常用的就是单因素（~ condition）和带批次校正的（~ batch + condition）。

为什么会需要批次？举个例子：你周一跑的样本和周五跑的样本，批次效应很常见。如果不加到 design 里，DESeq2 会把它当成生物学差异——结果就是一堆假阳性。

关于交互项：如果你研究的是「处理在不同时间点的效果」，就需要 ~ condition + time + condition:time。复杂设计我建议你画个表再写公式，别靠脑子硬想。

注意：design 公式里的变量必须在 colData 中存在。如果你写了 ~ batch + condition，但 colData 里没有 batch 列——DESeq2 会直接报错。我年轻时犯过这错，找 bug 找了半小时，最后发现是 variable 名称写错了。

2.5 构建 DESeqDataSet

前面东西准备好了，一行代码搞定：

library(DESeq2)

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

如果你之前用 tximport 导入的数据，那就用另一个函数：

dds <- DESeqDataSetFromTximport(
  txi = txi,
  colData = colData,
  design = ~ condition
)

构建完的 dds 是一个 SummarizedExperiment 对象。你可以用 dim(dds)、colData(dds)、assay(dds) 查看内容。

我个人习惯在构建后先做一步筛选：

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

去掉那些在所有样本中 read count 总和不到 10 的基因。这算是个经验值，能减少后续计算量，还能过滤掉很多噪音。

关于 pre-filtering：有些人喜欢直接对所有基因跑，觉得「万一淹没的基因有意义呢？」我做过实验对比——先过滤一下，结果基本一致。而且跑得更快，p-value 的分布也更稳健。你试试就知道了。

2.6 常见错误与避坑

行名或列名不匹配：countData 和 colData 的样本名不一致，DESeq2 会报 "sample names in colData do not match"。改名字就好。
countData 不是整数：DESeq2 要求整数。如果你输入了小数，它会尝试取整。我建议你提前检查：all(countData == floor(countData))。
design 公式写反了：~ condition 没问题。但如果你写成 ~ 1 + condition，多了一个截距项——DESeq2 可以用，但含义变了。保持简单。
colData 里有 NA：DESeq2 会自动删掉含有 NA 的列。但这可能会打乱你设计。最好提前清理干净。

我曾经帮同事调试一个项目，他的 dds 构建完一直报错。我一看——colData 里把 "control" 写成了 "contorl"。就一个字母，跑了一整天都是错的。你说冤不冤。

好了，dds 构建完成。之后就可以直接跑 DESeq()，然后拿结果。记住：数据输入这步花再多时间检查都不为过。花 10 分钟清理数据，剩下 3 小时的调试时间——这笔账不亏。