📚 课程目录第4章数据预处理

4. 数据预处理

数据预处理是整个单细胞分析里最容易被轻视、也最容易翻车的环节。我做过的第一个单细胞项目，上来就跳过了仔细的标准化，直接用原始count跑了PCA，结果前三个主成分全是批次效应……嗯，从那以后我再也不敢偷懒了。

咱们这一章要聊的，说白了就是把原始的UMI计数矩阵，洗成适合后续CCA、MNN校正的样子。你想想看，原始数据里每个细胞的有效分子数差好几倍，高表达基因和低表达基因不在一个量级，不洗一洗，后面的批次校正根本没法做。

4.1 标准化：让细胞之间能公平比较

标准化干的事很简单：消除不同细胞间测序深度差异。我最常用的方法是“library-size normalization + log1p转换”。

标准化的核心公式：

normalized_counts = log1p( (raw_counts / total_counts_per_cell) * scale_factor )

常见的scale_factor取10000，也就是CPM（counts per million）的对数。但在实际项目中，我遇到过一些极端的异质性样本——有的细胞文库特别小（只有几百UMI），这时候直接除scale_factor会让稀疏矩阵变得极其不平衡。我的习惯是先过滤掉低质量细胞（UMI小于1000的），再标准化。

另外，SCTransform也是一个不错的选择。它用正则化负二项回归拟合方差-均值关系，能保留更多生物学变异。不过计算代价大，我一般做10万以内的细胞用SCT，更大的数据集还是老老实实走LogNormalize。

我的小技巧：如果后续要做CCA，建议统一用LogNormalize。因为CCA对非负值更友好，SCT出来的残差可能有负值，不做额外处理会影响CCA的交互相似度计算。

下面是一个典型的标准化代码（R语言Seurat生态）：

library(Seurat)
# 假设已经创建了Seurat对象，原始数据在RNA assay的counts里
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize",
                      scale.factor = 10000)

4.2 高变基因选择：只挑信息量最大的基因

标准化的数据里有两万多个基因，但大部分基因在细胞间表达差异很小，贡献了噪声。高变基因选择就是把这“大部分”丢掉，只留下变异大的那几百到几千个基因。我习惯选2000个，够用。

Seurat提供了三种方法：vst（方差稳定变换）、mean.var.plot（均值-方差分箱）和dispersion（离散度）。我个人最推荐vst，原因是它对稀疏数据更鲁棒，而且速度很快。

我曾经踩过的坑：有一次我选了10000个HVG去做CCA，结果很多批次特异的基因被混了进来。校正后反而破坏了生物信号。后来我规规矩矩用2000个，效果就对了。数量不是越多越好。

选高变基因的代码：

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst",
                             nfeatures = 2000)

选完以后，记得看一眼“平均表达量 vs 方差”的散点图，通常高变基因会呈“火山”状分布。如果发现所有基因都挤在左下角，说明标准化可能没做对。

另外，下面这个表格把三种方法的特点对比一下，方便你根据数据量选择：

方法	原理	适用场景	我个人的偏好
vst	拟合loess曲线，标准化方差	稀疏数据、中等大小数据集	★★★★★ 日常首选
mean.var.plot	按表达量分箱，选离群基因	表达量差异大的数据（比如混合细胞类型）	★★★☆☆ 特殊场景
dispersion	直接取变异系数	数据噪声较小时可用	★★☆☆☆ 不推荐

4.3 降维准备：给PCA和CCA铺好路

降维准备这一步，说白了就是“把数据整理好，让PCA能直截了当地跑”。很多人以为直接拿标准化后的矩阵就能PCA了，其实不对。

我习惯把这一步拆成两个动作：ScaleData（中心化+缩放）和RunPCA前的检查。

ScaleData：为什么必须做？

高变基因的表达量范围可能差很多——有的基因平均表达10，有的只有0.1。PCA对量级敏感，如果不缩放到单位方差，高表达基因会主导前几个主成分。ScaleData做的就是Z-score变换：z = (x - mean) / sd。

实操注意：ScaleData默认会回归掉细胞周期、线粒体比例等无关变异。如果你做批次校正，建议先不要回归“批次”这一列，留到后面MNN或CCA去处理。我曾经因为提前回归了批次，导致校正效果变差。

pbmc <- ScaleData(pbmc, vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"))

降维前的最后检查

跑PCA之前，我习惯把数据再做个可视化排查：

检查scale后的表达矩阵是否有NaN或Inf（常见于某些细胞所有高变基因表达为0）
检查细胞和基因的总和是否合理
如果用了多个样本，先画个合并后的基因-基因热图，看看有无极端批次效应

这些检查看起来很笨，但救过我很多次。上次有个同事的脚本跑出来PCA第一主成分全是线粒体基因，就是因为忘关了线粒体基因的回归——你看，预处理出问题，后面全白搭。

# 检查NaN
any(is.na(pbmc[["RNA"]]@scale.data))

# 跑PCA
pbmc <- RunPCA(pbmc, npcs = 50, verbose = FALSE)

PCA做完了，就可以拿前20~50个主成分去喂给CCA或MNN了。我个人经验是，如果数据比较简单（比如同一个组织不同技术批次），30个PC足够；如果混合了多种细胞状态，可以取到50个。

整个数据预处理的逻辑流，我画成了一张图，方便你回忆：

看到这张图你就明白了：标准化是根基，高变基因选择是“去噪”，降维准备是“统一量纲”。三步缺一不可。下一个章节咱们就拿着这份预处理好的数据，正式进入CCA实战。