📚 课程目录第3章批次效应检测方法

数据预处理：归一化、缺失值处理与异常值过滤

数据预处理这一步，看着像杂活，实际上整条流水线的成败一大半都压在这。我踩过不少坑，后来才明白——原始表达量数据直接往下游扔，结果基本没法看。今天聊聊归一化、缺失值和异常值这三个核心动作，帮你把底子打扎实。

3.1 为什么必须归一化

说白了，不同样本的测序深度、文库大小、基因长度都不一样。你直接拿原始counts去比较，好比拿不同秤称出来的重量比大小——没意义。归一化的目标就一个：让样本间的表达量能公平对比。

不同场景要选不同的归一化策略。我习惯先想清楚一个问题：我的下游分析需要什么尺度的数据？

方法	适用场景	核心思想	常用工具
TMM (Trimmed Mean of M-values)	差异表达分析（edgeR, limma）	去除高表达基因的干扰，取修剪后的均值比	edgeR::calcNormFactors
RPKM/FPKM	基因长度不同时比较基因间表达	每百万reads每千碱基的reads数	stringTie, cufflinks
TPM	跨样本、跨基因比较	先归一化基因长度，再归一化文库大小	RSEM, Salmon, kallisto
CPM (Counts Per Million)	快速粗筛，简单直接	每个count除以总counts乘百万	edgeR::cpm
Quantile Normalization	芯片数据或多组学整合	强制让每个样本的分位数分布一致	preprocessCore

我的习惯：如果做差异表达，优先TMM。如果做基因表达丰度比较，用TPM。如果只是可视化，CPM够用。别一上来就Quantile——它抹掉了很多生物变异，我吃过这个亏。

3.2 TMM归一化实战

TMM是差异表达分析里的老大哥。edgeR的calcNormFactors默认就用它。原理不复杂：选一个参考样本，对每个样本计算M值（log2 fold change）和A值（平均log2表达量），砍掉两端的极端值，取加权平均。

上代码看看：

library(edgeR)

# 假设counts是基因×样本的表达矩阵
counts <- matrix(rpois(1000, lambda = 10), nrow = 100, ncol = 10)
rownames(counts) <- paste0("gene", 1:100)
colnames(counts) <- paste0("sample", 1:10)

# 创建DGEList对象
y <- DGEList(counts = counts)

# 计算TMM归一化因子
y <- calcNormFactors(y, method = "TMM")

# 查看归一化因子
y$samples$norm.factors

# 获取归一化后的CPM
cpm_norm <- cpm(y, normalized.lib.sizes = TRUE)

我曾经翻车的地方：calcNormFactors要求输入是整数counts。如果你手上有FPKM或TPM，千万别直接往里喂。老老实实回到原始counts矩阵。

3.3 RPKM/FPKM与TPM的区别

RPKM和FPKM当年很火，但现在我基本只用TPM了。为什么？RPKM的分子和分母用的都是总reads数，但总reads数本身受高表达基因影响很大。TPM先归一化基因长度，再归一化文库大小，这样跨样本比较的时候更稳定。

换算公式也好记：

# RPKM计算（手动实现）
rpkm <- function(counts, gene_lengths, total_reads) {
  rpm <- counts / (total_reads / 1e6)
  rpkm <- rpm / (gene_lengths / 1000)
  return(rpkm)
}

# TPM计算
tpm <- function(counts, gene_lengths) {
  rate <- counts / gene_lengths
  tpm <- rate / sum(rate) * 1e6
  return(tpm)
}

用Salmon或kallisto直接就能得到TPM，省心。我最近一个项目里处理10个公共数据集，全部用kallisto跑了一遍quant，然后直接读TPM，再转成log2(TPM+1)做后续分析，效果不错。

3.4 缺失值处理

表达量矩阵里出现0值是很正常的，但如果是技术性缺失（比如某个基因在某个样本里就是没测到），就需要处理。我个人把缺失值分成两类：

技术性缺失：文库构建或测序失败导致。需要补值或删掉。
生物学缺失：基因确实不表达。保留0值，别乱补。

区分这两类的方法：看同组内的其他样本。如果其他样本都有表达，就这个样本是0，大概率技术问题；如果整组都是0，那就是生物学沉默。

# 简单的缺失值标记
# 假设counts矩阵里0值占比过高可能是技术缺失
missing_ratio <- apply(counts, 2, function(x) sum(x == 0) / length(x))
# 如果某个样本缺失率超过80%，建议剔除该样本
samples_to_keep <- names(which(missing_ratio < 0.8))
counts_filtered <- counts[, samples_to_keep]

# 对个别基因的缺失值用KNN插补（谨慎使用）
library(impute)
# 注意：KNN插补前需要先做log转换或者至少保证数据不是counts
# 这里以log2(CPM+1)为例
log_cpm <- log2(cpm(counts_filtered) + 1)
imputed <- impute.knn(log_cpm, k = 10)
log_cpm_imputed <- imputed$data

注意：我一般不推荐对counts数据直接插补，更不要用均值填充。零膨胀的counts分布用均值填充会引入严重偏差。如果非要做，至少先用CPM或log转换。

3.5 异常值过滤

异常值分两种：样本级别的异常和基因级别的异常。样本级别的异常通常用PCA、聚类树、或者样本相关性热图来发现。基因级别的异常则根据表达量分布来检测。

我常用的样本异常过滤步骤：

计算每个样本的总reads数，看是否有数量级差异
做PCA，标记离群样本
计算样本间spearman相关系数，相关系数低于0.8的样本标记为可疑
检查基因检出率（表达量>0的基因数）

# R代码示例：基于PCA剔除离群样本
pca_res <- prcomp(t(log_cpm), scale. = TRUE)
# 计算每个样本在PC1和PC2上的马氏距离
mahal_dist <- mahalanobis(pca_res$x[, 1:2], 
                           center = colMeans(pca_res$x[, 1:2]), 
                           cov = cov(pca_res$x[, 1:2]))
# 以0.01的p值阈值筛选异常样本
outliers <- which(p.adjust(pchisq(mahal_dist, df = 2, lower.tail = FALSE), 
                           method = "fdr") < 0.01)
if (length(outliers) > 0) {
  message("发现异常样本: ", paste(colnames(log_cpm)[outliers], collapse = ", "))
  log_cpm_filtered <- log_cpm[, -outliers]
}

3.6 本章知识体系总览

下面这张图把预处理的关键流程串起来了，你看一眼就明白：

3.7 常见问题与避坑指南

避坑1：不要对counts做log转换后再归一化

log转换前必须已经完成文库大小归一化，不然低表达基因的log操作会放大噪音。我见过有人直接把counts取log2(counts+1)然后做TMM——结果归一化因子全乱了。

避坑2：缺失值补全后记得检查分布

补完记得画个密度图，看看插补后的数值是否破坏了原始分布形态。如果出现突兀的峰值，说明插补参数要调。

避坑3：异常值不一定是坏样本

有时候离群样本恰恰是生物变异的体现。比如疾病组里一个极端表型的样本。我建议先看临床信息，再决定是保留还是剔除。

我个人的预处理checklist：

输入数据确认是原始counts矩阵（整数，无负值）
运行calcNormFactors，检查norm.factors是否在0.7~1.3之间，超出太多说明某样本质量差
用boxplot看归一化后的logCPM分布是否一致
用PCA看是否有明显离群点，结合实验设计判断
导出clean表达矩阵，同时记录下所有被剔除的样本ID，放进论文方法部分

数据预处理这块看起来琐碎，但每一步都影响下游的结果。你想想看，归一化做错了，差异表达分析出来的差异基因可能有一大半是文库大小差异造成的，根本不是真正的生物学差异。所以别急着跑统计模型，先把数据整干净了再动手。