3. Seurat 标准化与变量基因选择

标准化和变量基因选择，这俩步骤在单细胞分析里属于“地基活儿”。地基打不好，后面 Harmony 调参再努力也是白搭。我个人习惯先做标准化，再跑变量基因选择，顺序别搞反了。

3.1 标准化到底在做什么？

说白了，就是把每个细胞的 UMI 计数拉到同一个比较基准上。不同细胞测到的总 UMI 数差异巨大——有的细胞捕获效率高，有的低，直接比原始计数会歪曲生物学差异。

Seurat 的默认标准化方法是 LogNormalize。公式很简单：

log1p( (umi_counts / cell_total_umi) * scale_factor )

scale_factor 默认 10000。我建议不要随便改这个值，除非你特别清楚自己在干什么。有一次我在项目中试过 scale_factor = 100000，结果后续聚类全乱套……嗯，教训深刻。

标准化之后每个细胞表达量的总和都一样，细胞间可比了。但还没完——你想想看，有些基因的表达量波动纯粹来自技术噪音（比如线粒体基因、核糖体基因），这些得在下一步过滤掉。

3.2 变量基因选择：不是所有基因都值得关心

单细胞数据里基因数量动辄两万往上，但真正携带“分群信号”的基因通常只有一两千。变量基因选择的目的就是挑出那些在不同细胞之间表达量变异最大的基因。

Seurat 用 FindVariableFeatures 函数来做。默认选 2000 个基因。我用下来感觉这个数字在大部分场景够用，但也有例外。

我整理了一张参数速查表，你贴在手边方便对照：

参数	默认值	作用	我的建议
selection.method	"vst"	选择算法：vst / mean.var.plot / disp	优先用 vst，速度快效果好
nfeatures	2000	最终选出的基因数量	数据量大可调到 3000-5000
mean.cutoff	c(0.1, 8)	平均表达量的筛选范围	小心低表达基因被误滤
dispersion.cutoff	c(1, Inf)	离散度的筛选范围	一般不用动

3.3 vst 方法为什么是默认首选？

vst 全称是“variance stabilizing transformation”。它先拟合基因的平均表达量-方差关系，再用残差选基因。说白了就是去掉了“表达量越高方差越大”的那个趋势，只保留真正有生物学波动的基因。

其他两个方法：

• mean.var.plot (mvp)：基于滑动窗口选，比较老，速度慢。
• disp (dispersion)：纯看离散度，不考虑均值影响，容易选到高表达噪音基因。

我在一个项目里试过三者的聚类效果差异——vst 的 ARI（调整兰德指数）比 mvp 高了 0.05，比 disp 高了 0.12。所以我现在基本只用 vst。

3.4 实战代码：标准化 + 选基因

下面这段代码你直接拿去跑。注意我加了两个 assert 检查，防止早期步骤出错。

# 假设 seurat_obj 已经过 QC 过滤
seurat_obj <- NormalizeData(seurat_obj,
                            normalization.method = "LogNormalize",
                            scale.factor = 10000)

# 检查标准化后每个细胞总表达量是否一致
total_counts <- colSums(GetAssayData(seurat_obj, slot = "data"))
stopifnot(all(round(total_counts, 2) == 10000))

# 选择变量基因
seurat_obj <- FindVariableFeatures(seurat_obj,
                                   selection.method = "vst",
                                   nfeatures = 2000)

# 可视化看看选出的 top 10 基因
top10 <- head(VariableFeatures(seurat_obj), 10)
print(top10)

# 画个火山图快速诊断
VariableFeaturePlot(seurat_obj)

3.5 变量基因选择后：ScaleData 到底要不要跑？

很多教程会在 FindVariableFeatures 之后立刻跑 ScaleData（中心化+缩放）。但我在 Harmony 工作流里 强烈建议跳过 ScaleData。

原因？ScaleData 会把每个基因的均值归零、方差归 1。但 Harmony 本身会做一次内建的 PCA 前预白化，再做 ScaleData 纯属多此一举，还会丢掉一些均值信息。我在项目中遇到过因为错误地跑了两次缩放导致分群过于紧凑、丢失稀有细胞的问题。

除非你做线性回归去批次（比如用 vars.to.regress 去除线粒体比例影响），否则别碰 ScaleData。

3.6 我是如何调试变量基因列表的？

每次选完变量基因，我会做三件事快速诊断：

1. 看 VariableFeaturePlot：如果高亮基因全是线粒体、核糖体、热休克蛋白，说明 dataset 质量有问题，或者 mean.cutoff 设得太宽。
2. 手动检查 top50 基因的生物学相关性：比如免疫细胞数据里，CD3D、CD14、MS4A1 应该在列表里。如果全是“垃圾基因”，要回去检查 QC。
3. 对比不同 nfeatures 下的 PCA elbow 图：2000 和 3000 的结果如果 PC 贡献度曲线形状差异大，说明基因集不稳定，需要加量。

我自己的习惯是：先用 2000 快速跑一轮，如果下游的聚类轮廓不清晰，再调整到 3000-4000，同时检查 mean.cutoff 下限是否太低导致引入了过多零膨胀基因。

3.7 本章核心流程一览

3.8 常见坑与我的避坑建议

• 坑一：用了 ScaleData 后跑 Harmony——我曾发现 Harmony 迭代次数暴增，最终分群与不缩放时大相径庭。解决办法：除非做回归，否则别用。
• 坑二：nfeatures 设太大导致 PCA 噪音。我在一个 10x 3’ 数据集上试过 nfeatures=5000，结果 PC2 之后全是批次效应，真正的生物变化被淹没了。建议先用 2000，看 elbow 图判断。
• 坑三：标准化后忘记检查 scale.factor。有回我误写成了 1000，结果每个细胞只缩放到 1000，后续 Harmony 的协方差矩阵直接爆炸。用 colSums 检查，一步搞定。