📚 课程目录第4章输入数据格式

第 4 章：输入数据格式 — count 矩阵与样本信息表

说到 DESeq2 的输入数据，说白了就两个东西：count 矩阵和样本信息表。我刚开始用的时候，总觉得这步简单，拿来就用，结果没少在数据格式上栽跟头。今天咱们把这两个家伙彻底搞清楚。

4.1 count 矩阵长什么样？

count 矩阵，就是每个基因在每一个样本中的原始 read 计数。行是基因，列是样本。注意，这里必须是整数，不能是 RPKM、FPKM 或者 TPM 这类归一化后的值。为什么？DESeq2 内部用的是负二项分布，它依赖原始计数来估计离散度。你给个浮点数进去，算法直接懵掉。

一个标准的 count 矩阵大概像下面这样：

gene_id   sample1 sample2 sample3 sample4
AT1G01010 1234    987     1567    2110
AT1G01020 56      34      89      102
AT1G01030 0       2       1       0
...

小提示：我习惯把基因 ID 放在第一列，并且保证它没有重复。如果有多行重复的基因名，DESeq2 会报错，因为内部要用行名做索引。

4.2 样本信息表的作用

样本信息表告诉我们每个样本对应的实验条件。比如处理组和对照组，时间点，批次等。DESeq2 根据这些信息来构建模型公式。

它的格式很简单：第一列是样本名字，必须和 count 矩阵的列名 完全一致。后面几列是实验变量。举个例子：

sample   condition   batch
sample1  control     batch1
sample2  control     batch1
sample3  treatment   batch2
sample4  treatment   batch2

你看，condition 列里只有两个水平：control 和 treatment。这就是 DESeq2 要对比的变量。batch 列如果也放在模型里，就可以校正批次效应。

千万注意：我曾经遇到一个项目，样本信息表的行名写成了 “Sample_1”，但 count 矩阵的列名是 “sample1”。结果下划线匹配不上，直接报错。检查行名对齐是数据预处理中最多数的小坑。

4.3 两个文件怎么对齐？

最简单的方法：让样本信息表的行名（第一列）和 count 矩阵的列名完全同顺序、同拼写。很多教程说 DESeq2 会自动对齐，但我不建议完全依赖它。手动排序更稳妥。

看一下我对齐的步骤：

先读取 count 矩阵，把第一列设为行名。
读取样本信息表，把第一列设为行名。
用 all(colnames(countData) == rownames(colData)) 检查是否完全一致。
如果不一致，用 match() 或者直接重排列。

代码示例：

countData <- read.csv("counts.csv", row.names = 1)
colData <- read.csv("metadata.csv", row.names = 1)

# 确保样本顺序一致
colData <- colData[colnames(countData), , drop = FALSE]
# 确认
all(colnames(countData) == rownames(colData))  # 应该返回 TRUE

4.4 常见格式陷阱

嗯，这里我整理了几个我实际碰到的坑：

问题	后果	解决办法
count 矩阵中有缺失值 (NA)	DESeq2 报错，模型无法拟合	检查上游定量流程，保证每个基因每个样本都有整数计数
基因名有重复	DESeq2 合并行或报错	去重，保留表达量最高的转录本，或使用 ensembl ID 去重
样本信息表与 count 矩阵样本数不一致	部分样本丢失或报错	用 intersect 取交集，确保两者样本完全对应
列名包含特殊字符（如空格、连字符）	R 中将列名自动转换为非法名称，后续引用困难	建议统一用小写字母和下划线，如 `sample_01`

重点提醒：你最好在读取数据后，先用 str() 和 head() 看一眼。别嫌麻烦，这一步能省掉后面半小时的排错时间。

4.5 为什么必须是原始整数计数？

底层原因：DESeq2 假设每个基因的计数服从负二项分布，其方差依赖于均值。归一化后的值（比如 RPKM）破坏了原始的均值-方差关系，导致离散度估计完全失真。说白了，标准化之后再跑差异分析，就像拿烤过的面包去量面粉的质量——意义不大。

经验谈：我建议所有样本名用 `_` 分隔，不用空格、连字符、括号。比如 `WT_day0_rep1`。这样在 R 里直接用 `make.names()` 也基本不会变形。

4.6 构建 DESeqDataSet 对象

数据准备好了，一行代码就把它们拼起来：

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

注意这里的 design 参数。它是一个公式，告诉 DESeq2 我们的实验设计。如果还有批次，就写成 ~ batch + condition。顺序可以调，但核心比较变量放在最后。

避坑：我曾经把 batch 放在 condition 后面，结果 DESeq2 默认用最后一个变量做差异分析，导致我比较的不是处理 vs 对照，而是 batch1 vs batch2。所以，如果你用多因子模型，记得把 results() 中指定 contrast。

4.7 数据清洗：过滤低表达基因

DESeq2 内部其实会过滤，但提前过滤掉那些在所有样本里都只有几个计数的基因，能降低多重检验负担，也让结果更稳健。我通用的过滤标准：

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

这样保留至少在 3 个样本中计数 ≥10 的基因。当然你可以根据测序深度调整参数。

最后提醒一句：输入数据这事，看起来琐碎，但 80% 的运行报错都出在这里。多花 10 分钟检查格式，能省下后面几小时的 debug 时间。