第三章 数据准备与表达矩阵格式
做批次校正，最怕什么？数据还没跑就翻车。我见过不少新手，把精力全放在调参上，结果第一步的表达矩阵格式就写错了，RUVseq 直接报错。嗯，今天咱们专门聊聊怎么准备好这份“入场券”。
3.1 表达矩阵长什么样？
RUVseq 要求的输入，说白了就是一个整数计数矩阵。行为基因或转录本，列为样本，每个格子放的是比对上的 read 数。别用 FPKM 或 TPM，那些小数会被当成浮点数，RUVseq 的负二项模型根本没法处理。
我个人习惯用 DESeq2 的 DESeqDataSetFromMatrix 函数来构建，但前提是你得有一个干净的计数矩阵。来看个表格：
| 特征 | 要求 | 常见错误 |
|---|---|---|
| 数据类型 | 整数 | 用了 FPKM / RPKM 的小数 |
| 行名 | 唯一基因标识符 (ENSEMBL / SYMBOL) | 有重复或缺失 |
| 列名 | 唯一样本 ID | 包含空格或特殊符号 |
| 缺失值 | 不允许有 NA | 用了 NA 或用 0 填充但未检查 |
我曾经碰过一个项目，合作方给的矩阵列名带上了样本处理日期，R 读取后自动把空格转成点，结果 colnames 跟元数据表对不上。折腾了两个小时才发现问题。所以，列名里别放空格、横线、括号，统一用字母、数字和下划线最稳。
3.2 数据准备四步走
我习惯把这个流程拆成四个步骤，你可以直接对着做。
- 原始计数：从 featureCounts、HTSeq 或 Salmon 的量化结果读进来。
- 过滤低表达基因：至少要在几个样本里有一定计数，否则 RUVseq 的估计会不稳定。
- 归一化：一般用 upper quartile 或 TMM，但 RUVseq 内部也做归一化，所以这一步可跳过，不过我个人会先算一下看看有没有离群样本。
- 构建目标矩阵：确保是所有样本的并集基因，不要因为某个基因在一个样本里是 0 就删掉它——0 也是有效信息。
⚠️ 强调一点：过滤时别太激进。我遇到过有人把平均 count < 10 的基因全删了，结果把一些只在某批次特异表达的基因也丢掉了，RUVseq 的因子估计严重失准。
3.3 一个常见的错误源头——低表达基因的处理
为什么低表达基因会搞事情？RUVseq 里的负二项模型，对低计数非常敏感。如果某个基因在多数样本里都是 0，少数样本里突然蹦出 10 个 reads，那它会被当作一个很大的波动，直接影响你估计的 unwanted factors。你想想，这不就等于让噪声去主导校正方向吗？
所以我的建议是：先用 edgeR::filterByExpr 或自己设一个最小阈值。比如：
# 示例：保留至少在5个样本中 counts >= 10 的基因
keep <- rowSums(counts >= 10) >= 5
counts_filtered <- counts[keep, ]这里 5 和 10 是经验值，你可以根据你的测序深度调。我一般先跑一下 summary(rowSums(counts))，看看中位数在哪里，再决定阈值。
💡 小技巧：如果你用的是 UMI 数据，可以适当把阈值降低，比如 counts >= 3。因为 UMI deduplication 已经去掉了 PCR 重复，每个计数更可信。
3.4 表达矩阵格式的代码演示
假设你手里有一个 counts 矩阵（data.frame 或 matrix），还有一个包含批次信息的元数据表。我一般把两部分准备好，然后丢给 newSeqExpressionSet：
library(RUVSeq)
# 假设 counts 是 data.frame，行名是基因，列名是样本
counts_matrix <- as.matrix(counts)
# 元数据至少要有两列：样本 ID 和批次/条件
pheno <- data.frame(
sample = colnames(counts_matrix),
batch = c(rep("batch1", 3), rep("batch2", 3)),
condition = c("treat","treat","control",
"treat","treat","control"),
row.names = 1
)
# 构建 RUVSeq 的输入对象
set <- newSeqExpressionSet(
counts = counts_matrix,
phenoData = pheno
)
set # 查看对象结构注意：phenoData 必须是一个 AnnotatedDataFrame 或 data.frame，行名与 colnames(counts_matrix) 完全一致。我早期犯过一个错：copy-paste 时多了一个空格，R 没报错但样本名字符串对不上，后面的分析全乱了。从那以后，我每次都用 identical() 比对一下：
stopifnot(identical(rownames(pheno), colnames(counts_matrix)))3.5 要不要先做“基因过滤”？
这个问题我经常被问到。答案是：一定要，但方法要对。RUVseq 选择 negative control genes 的时候（比如用不上 empirical 方法），它会根据 counts 矩阵自己算。如果里面有太多噪声基因，筛选出来的 control set 也会不靠谱。
下面这张图展示了数据准备的核心逻辑：
⚠️ 警告：有些工具会输出“稀疏矩阵”格式（比如 Matrix 包的 dgCMatrix），RUVseq 不支持这种格式。你要先转成常规的 matrix，用 as.matrix() 就行。我见过有人直接拿 dgCMatrix 喂进去，RUVseq 没报错但结果全是 NA。
3.6 元数据——另一个陷阱
表达矩阵搞定后，元数据（phenoData）也别忽视。RUVseq 的 newSeqExpressionSet 要求你传一个 phenoData 参数，里面至少包含批次信息（或你关心的生物学变量）。注意，行名必须与表达矩阵的列名完全一致，顺序可以不一样（函数内部会匹配）。
我自己习惯把所有样本元信息放在一个 CSV 里，用 read.csv 读取后，用 merge 或直接提取需要的列：
# 从 CSV 读取元数据
meta <- read.csv("sample_info.csv", row.names = 1)
# 筛选需要的列：批次、条件
pheno <- meta[, c("batch", "condition")]
# 再次验证顺序
stopifnot(all(rownames(pheno) %in% colnames(counts_matrix)))
pheno <- pheno[colnames(counts_matrix), , drop = FALSE]
# 此时 pheno 的行名顺序与 counts_matrix 列名一致这段代码最后一行做了重新排序，让元数据和表达矩阵严格对齐。嗯，这一步看起来多余，但能让你避免后面莫名其妙的 “error in checkSlotAssignment” 错误。我以前偷懒跳过这个 pheno <- pheno[colnames(counts_matrix), , drop = FALSE]，结果跑了三个小时的模型全错了——白白浪费一个下午。
3.7 常见问题速查表
我把平时踩过的坑整理成一个小表格，你准备数据时对着 check：
| 问题 | 表现 | 解决 |
|---|---|---|
| counts 矩阵不是整数 | RUVseq 结果全是 NA 或负数 | 检查是否用了 FPKM，换回原始 count |
| 列名有空格/特殊字符 | 匹配元数据时失败 | 用 make.names() 清洗 |
| 行名重复 | RUVseq 报错 “duplicate rownames” | 合并重复基因（取平均或求和） |
| 存在全零基因行 | 影响方差估计，但模型仍能跑 | 过滤掉全零行，不然会浪费自由度 |
| 元数据列名不匹配 | 报错 “no columns found” | 严格用 colnames(counts) 来命名 |
💡 最后一个小提醒：如果你是用 Salmon 或 kallisto 产生的转录本定量，别直接拿来用。那个是 transcript-level 的计数（甚至可能已经是 TPM 了），RUVseq 一般建议用 gene-level 的 counts。用 tximport 聚合成基因级别，再丢给 newSeqExpressionSet，这才是正规操作。