3. 数据准备：表达矩阵、分组信息与基因注释

做差异分析，八成时间花在数据准备上。这话不夸张。我刚开始做 RNA-seq 时，经常是矩阵格式错了、分组顺序反了、基因 ID 没对上版本号，最后跑出来一堆奇怪的 p-value。这一章咱们就把这三个基础构件讲透，后面用 DESeq2 和 edgeR 时才能一路顺畅。

3.1 表达矩阵构建

说白了，表达矩阵就是一个二维表格：行是基因，列是样本，单元格是每个基因在每个样本中的 read count（原始整数计数）。注意，是整数，别用 FPKM/RPKM 之类的归一化值，DESeq2 和 edgeR 内部有自己的归一化方法。

我习惯从 featureCounts 或 HTSeq-count 的输出开始。每个样本一个文件，第一列是基因 ID，第二列是计数。你需要把它们合并成一个完整的矩阵。下面是我常用的方式——用 R 的 read.table 循环读入，再用 do.call(cbind) 拼起来：

# 假设你有一个向量 sample_files，里面是各个样本的计数文件路径
count_list <- lapply(sample_files, function(f) {
  df <- read.table(f, header = FALSE, row.names = 1, skip = 2)
  colnames(df) <- gsub("_counts\\.txt", "", basename(f))
  return(df)
})
count_matrix <- do.call(cbind, count_list)
# 看看结果
head(count_matrix[, 1:4])

注意，有些工具会在文件头写统计摘要，比如 featureCounts 默认前两行是额外信息，所以 skip = 2 跳过。千万记得检查列名是否跟样本名一一对应——我遇到过文件名跟实际样本弄混的坑，那是某个项目的脑洞命名规则导致的。

3.2 样本分组信息

光有矩阵还不够，我们需要告诉工具：哪个样本是对照，哪个是处理。分组信息就是一个数据框，至少有两列：sample 和 condition。condition 必须是个因子（factor），而且因子的水平顺序决定了谁是谁的对照。

举个例子：

group_info <- data.frame(
  sample = c("WT_1", "WT_2", "KO_1", "KO_2"),
  condition = c("WT", "WT", "KO", "KO"),
  row.names = "sample"
)
group_info$condition <- factor(group_info$condition, levels = c("WT", "KO"))
# 这样 WT 是基线，KO 相对于 WT 比较

嗯，这里要注意——levels 的顺序很重要！DESeq2 默认将第一个水平作为分母。如果 R 自动按字母序排序，你可能发现“KO”变成了对照，结果完全反转。我曾经在某个合作项目中，因为没显式设置因子水平，给生物学意义解释带来了麻烦……从那以后，我每次都在 factor() 里明确写 levels，不敢偷懒。

如果你的实验有批次效应，还可以加一列 batch，后面在 DESeq2 的 ~ batch + condition 公式里用上。

3.3 基因注释文件

差异分析结果里只有基因 ID（比如 ENSG...）是不够的。你要看懂结果，最好有基因符号、染色体位置、基因类型等注释。常见做法是从 GTF 文件解析，或者用 biomaRt 在线查询。

我自己的套路是：先建一个注释数据框，至少包含 gene_id（与表达矩阵的行名一致）和 gene_name。如果要用 GTF 文件，可以这样：

# 假设你有一个 GTF 文件，用 rtracklayer 读取
library(rtracklayer)
gtf <- import("genes.gtf")
anno <- as.data.frame(gtf)
# 提取基因级别的记录，保留 gene_id 和 gene_name
gene_anno <- unique(anno[anno$type == "gene", c("gene_id", "gene_name")])

但如果 GTF 很大（比如人类的 GENCODE 版本），我更喜欢用 org.Hs.eg.db 这类包做 ID 转换：

library(org.Hs.eg.db)
ens_ids <- rownames(count_matrix)
gene_symbols <- mapIds(org.Hs.eg.db, 
                       keys = ens_ids, 
                       column = "SYMBOL", 
                       keytype = "ENSEMBL",
                       multiVals = "first")
# 检查有多少没匹配上
sum(is.na(gene_symbols))

为什么用 multiVals = "first"？因为有些 Ensembl ID 映射到多个符号，这里取第一个避免异常。当然，你可以事后检查那些多对一的情况。

3.4 三样东西怎么搭在一起？

我画了一张图，把三个数据的准备流程串起来。你看到最后，表达矩阵、分组信息、基因注释文件必须保持样本顺序一致，基因 ID 一致，这样才能喂给差异分析工具。

3.5 常见数据结构的对比

我把三种数据的格式要求和注意事项整理成了一张表，方便你对照检查：