批次效应矫正：那些默认值背后的秘密

· 18 · Tutorials

批次效应矫正：那些默认值背后的秘密

你知道 Harmony 的 θ 默认值是 2，scVI 的 n_latent 默认是 10——但你知道为什么吗？知道什么时候调它们会让结果变好，什么时候会让一切崩掉吗？

这篇文章不讲"Harmony vs scVI vs Seurat MNN 谁更好"——那类文章已经很多了。这篇讲的是：当你拿到一个新数据集，那些关键参数在数学上意味着什么，它们的默认值为何在大多数情况下管用，以及哪些边界条件会让它们失效。

先建立直觉：批次效应是什么在数学上在说什么

批次效应的本质是：同一细胞类型，在不同实验条件下测出来的基因表达向量，在高维空间里系统性偏移。

矫正方法大体分两类思路：

  • 对齐型（Harmony、Seurat MNN）：把不同批次的细胞嵌入拉到一起
  • 生成模型型（scVI、scANVI、totalVI）：显式建模"批次是噪声变量"，学一个无批次的潜在表示

两者参数逻辑完全不同。下面分开讲。

一、scVI / scANVI 参数深潜

scVI 的核心是一个变分自编码器（VAE），目标函数是：

$$\mathcal{L} = \underbrace{\mathbb{E}{q(z|x,s)}[\log p(x|z,s)]}{\text{重建损失}} - \underbrace{\beta \cdot D_{KL}(q(z|x,s) | p(z))}_{\text{正则化损失}}$$

其中 $x$ 是观测到的基因计数，$s$ 是批次标签，$z$ 是潜在空间向量。理解这个公式，你就理解了 scVI 的所有核心参数。

1.1 n_latent：潜在空间维度

默认值：10

$z \in \mathbb{R}^{n_latent}$，是编码器学到的"无批次"细胞表示。

为什么默认 10 够用？

一般 scRNA-seq 数据集，细胞的生物学变异主要由有限数量的细胞类型/状态驱动。10 个维度足以区分 PBMCs 的 8-12 个主要细胞类型，以及它们之间的状态连续体。研究表明，对于 10k–50k 的典型数据集，前 10 个 PCA 成分能解释 > 70% 的方差——n_latent=10 的选择与这个经验规律对应。

什么时候应该调大？

场景 建议
数据集 > 100k 细胞、细胞类型极多（如全器官图谱） 调到 20–30
包含罕见细胞亚群（< 0.5%），需要精细区分 调到 15–20
发育轨迹数据（连续细胞状态） 调到 15–20
仅做批次矫正、不做下游分型 保持 10，不必大

什么时候会崩？

  • n_latent 过大（如 50+）：模型参数过多，在小数据集（< 5k 细胞）上容易过拟合——训练损失漂亮，但 UMAP 中批次点仍然分开（模型把批次信息"藏"进了多余的维度）
  • n_latent 过小（如 5）：信息瓶颈太窄，不同细胞类型的表示开始重叠，下游聚类精度下降

快速判断法： 训练完后跑 model.get_latent_representation()，然后做 PCA，看前 n_latent 个成分的解释方差曲线。如果曲线很快趋于平坦（比如第 6 个成分后基本为 0），说明你的 n_latent 可以更小，模型在浪费容量。

1.2 n_hidden：隐藏层神经元数

默认值：128

编码器/解码器每个隐藏层的神经元数。这控制的是非线性映射的容量。

为什么 128 够用： 对于一般的 2000–5000 个高变基因输入，128 个神经元的隐藏层能捕捉足够的非线性特征。增大它在大多数数据集上没有显著提升，只会增加训练时间。

例外： 输入基因数非常大（如全基因组转录组，10k+ 基因），建议同步调大 n_hidden=256，否则第一层压缩步幅过大，容易丢细节。

1.3 dispersion：负二项分布的离散参数

选项："gene" | "gene-batch" | "gene-label" | "gene-cell"
默认："gene"

scVI 用负二项（NB）或零膨胀负二项（ZINB）建模基因计数：

$$p(x_{ng} | \mu_{ng}, \theta_g) = \text{NB}(\mu_{ng}, \theta_g)$$

$\theta_g$（有时记为 dispersion）控制每个基因计数分布的"散布程度"：$\theta_g$ 越小，分布越分散（尾部越重）；$\theta_g \to \infty$ 时趋近于泊松分布。

选项 含义 适用场景
"gene" 每个基因一个 $\theta_g$，跨所有细胞/批次共享 绝大多数场景，默认首选
"gene-batch" 每个（基因 × 批次）一个 $\theta$，允许批次间离散度不同 批次间测序深度差异巨大（如 10x vs Smart-seq2 混合）
"gene-label" 每个（基因 × 细胞类型）一个 $\theta$ 有强先验知识，知道不同细胞类型对同一基因离散度差异大
"gene-cell" 每个（基因 × 细胞）各一个 $\theta$ 最灵活，但参数最多，需要大数据集（> 50k 细胞）

为什么 "gene" 默认会崩的场景：
当你混合了单端测序（10x）和全长测序（Smart-seq2）时，同一基因的计数分布差异极大——10x 产生稀疏的计数，Smart-seq2 产生更连续的高计数。此时用同一个 $\theta_g$ 建模两个批次，等于强迫一个参数拟合两个完全不同的分布。改用 "gene-batch" 能显著改善这类数据集的 ELBO。

1.4 gene_likelihood：基因计数的概率分布

选项："zinb" | "nb" | "poisson"
默认："zinb"（零膨胀负二项）

ZINB 在 NB 基础上加了一个伯努利门控，用来建模技术性零值：

$$p(x_{ng} = 0) = \pi_{ng} + (1 - \pi_{ng}) \cdot \text{NB}(0 | \mu_{ng}, \theta_g)$$

什么时候用 "nb" 而非 "zinb"：

  • 数据集较新、测序深度较高（> 3000 UMI/cell），零值较少
  • 数据已经过质控过滤，低质量细胞和 doublets 已去除
  • 用 "zinb" 时 UMAP 出现奇怪的细胞簇分裂——这通常是零膨胀参数在过拟合噪声

1.5 β（kl_weight）：KL 散度的权重

默认值：1.0

这是最容易被忽视、影响最大的参数之一。

回看损失函数：

$$\mathcal{L} = \text{重建损失} - \beta \cdot D_{KL}(q | p)$$

  • KL 项的作用：把后验 $q(z|x,s)$ 拉向标准正态先验 $p(z) = \mathcal{N}(0, I)$，强迫潜在空间平滑、连续、可插值。正是这个约束让不同批次的细胞在潜在空间里"混"在一起。
  • 重建项的作用：确保解码器能从 $z$ 还原出原始计数，保留生物学信息。

两者是拮抗关系：

β 值 效果
β = 0 纯自编码器，没有正则化，批次效应完全保留在 z 中
β < 1 轻正则化，重建好，批次矫正弱，适合批次效应很小的数据
β = 1（默认） 标准 VAE，平衡点，大多数场景适用
β > 1（β-VAE） 强正则化，批次混合好，但生物学细节可能被压缩丢失
β >> 10 潜在空间坍缩（posterior collapse），所有细胞 z 趋同，聚类失去意义

实际调参建议：

  • 默认 kl_weight=1.0，先跑，看 scVI 的 UMAP
  • 如果批次分离仍然明显：试 kl_weight=1.5–2.0
  • 如果细胞类型之间边界模糊：适当降低到 kl_weight=0.5–0.8
  • scVI 有 kl_weight warm-up 机制（前几个 epoch 从 0 线性增加到目标值），这个 warm-up 本身很重要，不要禁用

注意：scVI 的 Python API 在 train() 里通过 kl_weight 或 plan_kwargs={"kl_weight": ...} 传入，而不是模型初始化时。不同版本 API 位置略有不同，查文档确认。

1.6 alpha（scANVI 特有）：监督信号强度

scANVI 在 scVI 基础上加入了细胞类型标签的监督：

$$\mathcal{L}{scANVI} = \mathcal{L}{scVI} + \alpha \cdot \mathcal{L}_{classifier}$$

alpha 控制分类器损失的权重。默认 alpha=1。

  • 调大 alpha：分类器更强，标注细胞的表示更规整，但无标注细胞更依赖分类器预测（过拟合风险）
  • 调小 alpha：更多依赖无监督部分，适合标注细胞比例 < 10% 的场景

二、Harmony2 参数深潜

Harmony 的核心思路是迭代软聚类 + 质心移动：

  1. 在 PCA 空间中，对所有细胞做 K-means 软聚类，得到每个细胞属于每个簇的隶属度
  2. 计算每个簇内各批次细胞的质心偏差
  3. 将每个细胞的嵌入向各自批次质心的反方向平移（消除批次偏差）
  4. 重复直到收敛

Harmony2（最新版本）引入了软余弦矫正，允许矫正方向在空间中连续变化，而非对每个簇做硬性平移。

2.1 theta（θ）：批次多样性惩罚

默认值：2（可以是一个向量，每个批次协变量一个值）

$\theta$ 控制 Harmony 在计算软聚类质心时，对批次不均匀性的惩罚强度：

$$\text{簇内惩罚} \propto \theta \cdot \sum_b \frac{n_b}{n} \log \frac{n_b / n}{R_{kb}}$$

其中 $R_{kb}$ 是批次 $b$ 的细胞对簇 $k$ 的隶属度，$n_b/n$ 是批次 $b$ 的细胞比例。

直觉理解：$\theta$ 越大，Harmony 越努力让每个簇内各批次细胞的比例接近全局比例，即强制混合。

θ 值 效果
θ = 0 不惩罚，等于不做批次矫正
θ = 2（默认） 适中矫正，适合批次效应中等的数据
θ = 4–6 强矫正，批次混合好，但可能将不同细胞类型强行压在一起
θ > 8 过矫正风险——生物学差异开始被当成批次效应去掉

什么时候调高 θ：

  • 批次差异来源非常强（不同测序平台）
  • 数据集有经验证的同种细胞类型，UMAP 里却分成了两团

什么时候降低 θ：

  • 你的"批次"其实是真实的生物学差异（如疾病 vs 正常，时间点不同）
  • 矫正后发现一些本来应该分开的细胞类型被合并了

多个协变量时的 θ 向量：

harmony_kwargs = {
    "theta": [2, 1],  # 第一个协变量（测序批次）矫正强，第二个（供体）矫正弱
}

2.2 sigma（σ）：软聚类带宽

默认值：0.1

控制 K-means 软聚类的"模糊度"。$\sigma$ 越小，每个细胞越确定地只属于一个簇（趋向硬聚类）；$\sigma$ 越大，隶属度分布越均匀。

为什么 0.1 通常够： 这个值使得 PCA 空间中距离最近簇的细胞隶属度 >> 其他簇，同时又不完全变成硬分配。

崩的场景： 数据分布非常不均匀（细胞类型比例差异极大），小 $\sigma$ 时稀有细胞类型可能被划入相邻大簇，导致稀有群体在矫正后被"同化"。这时适当增大 $\sigma$（如 0.2）给稀有群体更多"缓冲"。

2.3 nclust：软聚类数目

默认值：min(100, N_cells / 30)（大约）

Harmony 内部进行 K-means 软聚类的簇数。这个参数经常被忽视，但影响显著。

  • nclust 太小：每个"超簇"混入了太多不同细胞类型，矫正时把真实异质性平均掉了
  • nclust 太大：超过实际生物学群体数，收敛慢，且每簇细胞数太少时质心估计不稳定

建议： 至少是你预期细胞类型数的 5–10 倍。50k 细胞的数据集，nclust=100–200 通常是合理的。

2.4 block_size：每次更新使用的细胞比例

默认值：0.05（即每次用 5% 的细胞）

Harmony 每次迭代并不用所有细胞更新参数，而是采样一个 block。这是一个随机梯度思路，作用是：

  • 降低内存占用
  • 为收敛过程引入随机性，避免局部极小值

什么时候调：

  • 数据集非常小（< 5k 细胞）：block_size=0.05 每次只用 250 个细胞，随机性过大，收敛不稳定。调到 0.2–0.5。
  • 数据集非常大（> 500k 细胞）：默认 5% 仍然很多（25k 细胞每次），可适当降到 0.01–0.02 加速每步。

2.5 max_iter_harmony 和 eps_harmony：收敛控制

默认：max_iter_harmony=10，eps_harmony=1e-4

Harmony 运行到目标函数变化 < eps_harmony 或达到最大迭代次数停止。

10 次迭代通常够用，但有时不够： 批次效应复杂、协变量多时，观察 Harmony 的 objective function 曲线（harmony_object.objective_harmony），如果到第 10 次还没收敛（曲线还在下降），调大 max_iter_harmony=20。

三、totalVI：多模态的额外参数

totalVI 同时处理 RNA 和蛋白质（CITE-seq / ASAP-seq）数据，在 scVI 基础上加了：

3.1 蛋白质背景先验

CITE-seq 的蛋白质信号里有非特异性背景噪声，totalVI 用混合模型建模：

$$p(\text{protein}) = \pi \cdot \mathcal{N}(\mu_{bg}, \sigma_{bg}) + (1-\pi) \cdot \mathcal{N}(\mu_{fg}, \sigma_{fg})$$

  • protein_background_prior_mean（默认 -10，log-normal 空间）：背景信号均值的先验
  • protein_background_prior_scale（默认 0.5）：背景先验的不确定性

什么时候这两个参数会崩：
如果你的蛋白质 panel 含有很多不在目标细胞上表达的蛋白（阴性对照大量），背景均值会系统性偏高。此时 protein_background_prior_mean 应该调得更高（如 -5），否则模型会把真实的低表达蛋白也判为背景。

3.2 protein_dispersion

类似 scVI 的 dispersion，控制蛋白质计数的离散度建模。通常保持默认（"protein"，每个蛋白一个参数），但当你有多个测量批次时，"protein-batch" 更合适。

四、通用调参流程

不管用哪个工具，建议按这个顺序：

1. 先跑默认参数
   → 看 UMAP：批次分离？细胞类型边界清晰？
   → 有批次分离 → 考虑加强矫正参数（θ↑ 或 kl_weight↑）
   → 细胞类型边界模糊 → 考虑降低矫正强度

2. 调一个参数，固定其他
   → 对比 LISI 分数（iLISI 衡量批次混合，cLISI 衡量细胞类型保留）
   → 两个分数要同时看，不能只追求高 iLISI

3. 检查生物学合理性
   → marker gene 是否在正确的簇里高表达？
   → 稀有细胞类型有没有被"同化"消失？

速查表

问题 可能原因 调参方向
UMAP 里批次仍然分开 矫正强度不足 Harmony: θ↑；scVI: kl_weight↑
不同细胞类型被合并 过矫正 Harmony: θ↓；scVI: kl_weight↓
稀有细胞类型消失 nclust 太小 / σ 太小 Harmony: nclust↑, sigma↑
训练不收敛 学习率 / warm-up 问题 检查 kl_weight 曲线，加 warm-up
大数据集 OOM block_size 太大 Harmony: block_size↓；scVI: batch_size↓
ZINB 出现奇怪小簇 零膨胀参数过拟合 换 gene_likelihood="nb"
scANVI 标注迁移差 alpha 太小 alpha↑，或增加有标注细胞比例

小结

批次效应矫正参数的选择，本质上是在两个对立目标之间找平衡：

  • 批次混合（integration） vs 生物学保留（biology preservation）
  • 模型容量（capacity） vs 泛化能力（generalization）

默认值是经过大量数据集验证的合理起点，但它们不是万能的。理解每个参数在数学上控制的是什么，你才能在默认值失效时知道往哪里走——而不是盲目网格搜索。

相关工具文档：

scRNA-seq 批次效应 scVI Harmony 参数调优 单细胞