第3章 表达数据预处理：让数据说话前先“洗个澡”
我做了这么多年RNA-seq分析，最怕什么？不是模型调参，而是拿到原始表达矩阵后直接往里塞。你想想看，原始counts里有多少垃圾：那些背景噪声、批次偏移、基因长度带来的偏差……不处理干净，模型学到的全是假信号。
这一章我就带着你过一遍我平时必做的三步预处理。说白了，就是给数据“洗澡”：先剃掉没价值的低表达基因，再换成统一度量衡（TPM/FPKM/RPKM），最后抹掉批次之间的差异。
3.1 过滤低表达基因：先把垃圾清出去
为什么一定要过滤？我举个实际例子。当年我跑一个肿瘤样本对比项目，原始矩阵有6万多个基因，但半数以上在超过一半样本里是零。这些基因不仅没用，还让后续的归一化产生很大偏差——它们比真实信号更容易被离群值左右。
常用的方法就是设定一个CPM（counts per million）阈值。比如要求基因在至少两个样本中CPM > 1。EDGER的filterByExpr函数就能自动搞定，但我喜欢自己写，心里踏实。
# 假设 count_matrix 是 raw counts (基因×样本)
library(edgeR)
# 手动过滤：每个基因至少在一半样本中 CPM > 1
keep <- rowSums(cpm(count_matrix) > 1) >= ncol(count_matrix) / 2
filtered_counts <- count_matrix[keep, ]
还有种做法是先过滤掉在少于最小文库组大小样本中表达量为零的基因。EDGER的filterByExpr默认就采用这个逻辑——它会根据实验设计自动计算最小文库组的大小。
# 使用 edgeR 自动过滤
y <- DGEList(counts = count_matrix, group = group)
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes = FALSE]
keep.lib.sizes参数留为TRUE，结果发现归一化因子计算错误——因为你删掉了基因，但每个样本的总文库大小没变，导致后续归一化因子的分母用错。所以，一定要显式设置keep.lib.sizes = FALSE，让edgeR重新计算文库大小。
3.2 TPM / FPKM / RPKM 转换：统一度量衡
原始counts受两个因素影响：基因长度和文库测序深度。长基因更容易被检测到，深度的样本counts整体更高。所以，直接拿counts跨基因或跨样本比较，是不公平的。
三种转换的区别，我直接列个表你一看就懂：
| 度量 | 分子 | 分母 | 跨样本可比性 | 跨基因可比性 |
|---|---|---|---|---|
| RPKM | reads per gene | 基因长度(kb) × 总reads(百万) | 弱 | 有 |
| FPKM | fragments per gene | 基因长度(kb) × 总fragments(百万) | 弱 | 有 |
| TPM | RPK (per gene) | ∑RPK per million | 强 | 有 |
我个人习惯：一律用TPM。原因很简单——TPM在样本间才是真正可比的。你看RPKM/FPKM，每个样本除的是自己的总reads，所以同一个基因在不同样本里的TPM之和总是相同的，这保证了一个基因在样本A里的表达量高，就一定是因为它相对其他基因表达高了，而不是因为A样本深度大。
转换的代码也很直接：
# 从 counts 计算 TPM
# gene_length 是每个基因的长度向量 (bp)
counts_to_tpm <- function(counts, gene_length) {
# 计算 RPK (reads per kilobase)
rpk <- counts / (gene_length / 1000)
# 每个样本的 RPK 求和
rpk_sum <- colSums(rpk)
# TPM = RPK / sum(RPK) * 1e6
tpm <- t(t(rpk) / rpk_sum) * 1e6
return(tpm)
}
tpm_matrix <- counts_to_tpm(filtered_counts, gene_length)
tximeta或tximport从转录本定量直接得到转录本级别长度，再聚合到基因。另外，别在过滤低表达之前做TPM，因为过滤后的计数总和更小，TPM值会整体放大，产生偏差。
3.3 批次UQ归一化：抹平“实验手抖”
批次效应是RNA-seq里最头疼的问题之一。我接过一个项目，两个实验室测的同一批细胞，结果热图一画，样本直接按测序批次聚类——而不是按细胞类型。这哪行？
归一化方法很多，TMM常用，但上四分位数（UQ）归一化是我在跨平台数据整合时优先选用的方法。它的思路很简单：假设大多数基因在样本间不是差异表达的，那么每个样本的75%分位数（即上四分位数）应该差不多。我们把这个共同的75%分位数作为一个“锚”，算出每个样本的缩放因子。
EDGER已经内置了这个方法：
# 在过滤后的 DGEList 上计算 UQ 归一化因子
y <- calcNormFactors(y, method = "upperquartile", p = 0.75)
# 查看归一化因子
y$samples$norm.factors
拿到归一化因子之后，EDGER会自动用它缩放文库大小，然后cpm(y, normalized.lib.sizes=TRUE)就能得到UQ归一化的CPM。我自己也写过纯手动的版本，方便查错：
# 手动计算 UQ 归一化因子
uq_values <- apply(y$counts, 2, function(x) {
quantile(x[x > 0], probs = 0.75)
})
uq_ref <- mean(uq_values)
norm_factors <- uq_values / uq_ref
预处理完成，检查一下效果
做完这三步，我一般会做个PCA或热图，看看样本是否按实验分组聚类，而不是按批次。如果效果不理想，就返回调整过滤阈值或者换归一化方法（比如用TMM试试）。