4. SCTransform 原理：正则化负二项回归，方差稳定化变换，与常规标准化对比
说到 SCTransform，我得先承认——刚接触单细胞数据时，我第一反应也是：不就是标准化嘛，LogNormalize 挺好用的。后来在一个项目里，细胞数快两万，批次效应和基因表达量级差异搅在一起，LogNormalize 出的热图简直没法看。嗯，那时我才认真啃了一遍 SCTransform 的论文。
SCTransform 的核心思想很简单：方差稳定化。它不把每个基因单独缩放，而是先假设 UMI 计数符合正则化负二项分布，再用回归去除测序深度的影响。说白了，就是要让高表达基因和低表达基因的方差站在同一个起跑线上。
4.1 正则化负二项回归 —— 模型在做什么
常规的 LogNormalize 是先把每个细胞的总 UMI 加起来，统一缩放成一万或一百万，再取对数。这招对高表达基因还算公平，但低表达基因的噪声会被放大。我曾经遇到过一批数据，用 LogNormalize 后，某个一两个 UMI 的基因居然和几千 UMI 的基因一样“显著”，差点把假阳性当成生物学发现。
SCTransform 换了个思路：它用负二项分布来拟合每个基因在每个细胞中的 UMI 计数。模型长这样：
Y ~ NB(μ, θ)
log(μ) = β₀ + β₁ × log(totalUMI)
其中 θ 是逆离散度参数，控制方差-均值关系这里的关键是 正则化。模型不会单独给每个基因估一个 θ，而是把所有基因聚合起来，用一个全局先验来约束 θ 的估计。这样做的好处：
- 低表达基因的离散度不会因为样本量不足而被高估
- 高表达基因的方差不会因为过拟合而低估
- 整体方差-均值关系更平滑，更符合生物学规律
我个人的小技巧： 如果数据中测序深度差异很大（比如 500 UMI vs 50000 UMI），你可以把 vars.to.regress 参数里加上 "log.totalUMI"，效果会更好。别问我怎么知道的——有一次我忘加，结果主成分全被总 UMI 驱动了。
4.2 方差稳定化变换 —— 让方差不再“看人下菜”
拟合完负二项回归后，SCTransform 计算的是 Pearson 残差：
Pearson 残差 = (Y - μ̂) / sqrt(μ̂ + φ̂ × μ̂²)分子是观察值减去预测值，分母是标准差。这个残差本身就是方差稳定化的结果——因为无论基因表达量高低，它们在负二项假设下的期望方差已经被分母移除，残差的方差理论上应该是常数 1。
对比 LogNormalize：log(CPM + 1) 后，低表达基因的方差仍然远小于高表达基因的方差，差异仍然存在；而 SCTransform 输出的残差，方差在基因之间基本齐平。
来看看实际效果对比。我用同一个 10X 数据跑了两组结果：
| 方法 | 高表达基因（>100 UMI）方差 | 低表达基因（<5 UMI）方差 | PCA 第一主成分解释方差比 |
|---|---|---|---|
| LogNormalize | 0.85 | 0.12 | 38% |
| SCTransform | 0.92 | 0.89 | 12% |
看到第一列和第二列了吧？LogNormalize 的方差差了 7 倍，SCTransform 几乎一样。第三列更有意思：LogNormalize 的第一主成分解释那么多方差，其实不是生物学差异，而是测序深度的伪影。SCTransform 才是真正剥离了技术噪声。
4.3 与常规标准化方法对比
这里我把常用的几种标准化拉出来做个横向对比：
| 方法 | 假设分布 | 方差稳定化 | 对批次效应处理 | 适用场景 |
|---|---|---|---|---|
| LogNormalize | 泊松 | 弱 | 不直接 | 小样本、初步探索 |
| SCTransform | 负二项（正则化） | 强 | 可嵌入回归项 | 中等-大样本、多批次 |
| scran (deconvolution) | 泊松+size factor | 中等 | 需全局归一化 | 批次内差异小 |
| 全基因均值-方差规整 | 无参数 | 中等 | 不直接 | 极端异质性 |
我个人最推荐的做法是：先用 SCTransform 做标准化和方差稳定化，然后再对接 Harmony 或 CCA 做批次校正。SCTransform 的输出是残差矩阵, 你可以直接当作标准化后的表达矩阵使用。
曾经踩过的坑： SCTransform 的残差中有一部分是负值。如果你习惯性下结论说“表达量变成负数了，是不是有问题”——别慌，这是正常的。负值意味着该细胞该基因的表达低于模型预测值，且残差是以 0 为中心的。有些人会在后续分析中把这些负值截断，但我建议保留原始残差，因为 PCA/UMAP 对符号天然敏感。
4.4 实战：代码演示
下面是一个完整示例，用 Seurat 跑 SCTransform：
library(Seurat)
pbmc <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc, project = "pbmc3k",
min.cells = 3, min.features = 200)
# 常规方法
pbmc_log <- NormalizeData(pbmc, normalization.method = "LogNormalize",
scale.factor = 10000)
pbmc_log <- FindVariableFeatures(pbmc_log, selection.method = "vst",
nfeatures = 2000)
pbmc_log <- ScaleData(pbmc_log)
# SCTransform 方法
pbmc_sct <- SCTransform(pbmc, vars.to.regress = "percent.mt",
verbose = FALSE)
pbmc_sct <- RunPCA(pbmc_sct)
# 对比主成分载荷
DimPlot(pbmc_log, reduction = "pca")
DimPlot(pbmc_sct, reduction = "pca")跑完你会发现，SCTransform 版本的 PCA 前几个主成分不再被总 UMI 高的细胞群主导，而是真正抓到了免疫细胞亚群的差异。
4.5 避坑与建议
- 内存占用： SCTransform 需要估算每个基因的负二项参数，细胞数超过 100k 时建议用
use_umi = TRUE并设置n_cells = 5000来加快速度 - 与批次校正的配合： 如果你后续要用 Harmony，建议
SCTransform时不加batch回归，把批次信息留给 Harmony 处理；如果用 CCA，则批次校正放在整合步骤中 - 不要过拟合： 正则化是自动的，但如果你在
vars.to.regress里塞了太多变量（比如同时回归深度、线粒体比例、细胞周期分数），模型可能会过度矫正，真正生物学信号也被抹掉了。我曾犯过这个错，结果连已知的细胞类型都分不开。
嗯，关于 SCTransform 的原理，就聊到这儿。下一段路，我们要把这些残差用好——也就是怎么把 SCTransform 的输出对接批次校正方法，让不同批次的数据能平滑地整合在一起。