使用totalVI进行单细胞多组学数据联合概率建模.

• paper：Joint probabilistic modeling of single-cell multi-omic data with totalVI

0. TL; DR

TotalVI（Total Variational Inference）是一种端到端的深度概率生成模型，能够联合分析单细胞的 RNA 转录组和表面蛋白质组数据。该方法通过变分自编码器框架，显式建模 RNA 的技术噪声（如测序深度变异、有限灵敏度）和蛋白质数据的技术偏差（如背景抗体信号、非特异性结合）。

TotalVI 能够解耦蛋白质的前景与背景信号，整合不同批次或不同抗体面板的实验数据，插补缺失的蛋白质测量，并进行差异表达分析。作者在小鼠脾脏和淋巴结的 CITE-seq 数据集上验证了方法的有效性，展示了其在细胞亚型鉴定、批效应校正和跨模态相关性分析中的优越性能。

1. 背景介绍

1.1 单细胞多组学技术的发展

近年来，单细胞 RNA 测序（scRNA-seq）技术极大地推动了细胞异质性的研究，使得研究人员能够在全转录组水平上解析细胞类型和状态。然而，转录组水平并不总是与蛋白质丰度直接对应，而蛋白质才是细胞功能的直接执行者。流式细胞术（flow cytometry）长期以来一直是分析细胞表面蛋白质的金标准，但受限于荧光检测通道，只能同时检测少量标记蛋白。

CITE-seq 技术的出现弥合了这一差距。该技术通过使用结合有寡核苷酸条形码的抗体（barcoded antibodies），在单细胞水平上同时测量 RNA 表达和数百种表面蛋白质的丰度。这种方法不仅显著增加了可同时检测的蛋白质数量，还提供了与 scRNA-seq 兼容的高通量分析流程，为理解基因表达与细胞表型之间的关系提供了前所未有的机会。

1.2 数据分析的挑战

尽管 CITE-seq 提供了丰富的多模态数据，但其联合分析面临多重技术挑战：

• 模态特异性的技术噪声：RNA 数据受到测序深度限制、检测灵敏度不足（dropout 现象）和过离散（over-dispersion）的影响；而蛋白质数据虽然稀疏性较低，但存在显著的背景噪声，包括未结合的游离抗体（ambient antibodies）和非特异性抗体结合（non-specific binding），导致几乎所有细胞都检测到非零的蛋白质计数。
• 批次效应与数据整合：不同实验、不同实验室或不同抗体面板（antibody panels）引入的批次效应（batch effects）会掩盖真实的生物学信号。特别是当整合使用不同抗体面板的实验时，许多算法要求特征（features）完全匹配，这会导致蛋白质信息的损失。
• 缺乏统一的统计框架：现有方法通常将 RNA 和蛋白质分开处理，或仅用其中一种模态进行聚类，另一种用于验证。缺乏能够同时建模两种模态、量化不确定性并支持下游任务（如差异表达、去噪和缺失值插补）的端到端解决方案。

2. totalVI 方法

TotalVI 是一个基于深度生成模型的统计框架，使用变分推断（variational inference）联合建模 RNA 和蛋白质计数数据。

2.1 模型架构与生成过程

TotalVI 的输入包括：RNA 计数矩阵 $X \in \mathbb{N}^{N \times G}$（$N$ 个细胞，$G$ 个基因），蛋白质计数矩阵 $Y \in \mathbb{N}^{N \times T}$（$T$ 种蛋白质），以及批次指示变量 $s_n \in {0,1}^B$（$B$ 个批次）。模型假设每个细胞的观测数据由以下潜在变量生成：

• $z_n$：细胞的低维生物状态（latent cell representation），服从 Logistic Normal 分布 $z_n \sim \text{LogisticNormal}(0, I)$，确保 $z_n$ 位于概率单纯形上（非负且和为1），便于后续的原型分析（archetypal analysis）。
• $\ell_n$：RNA 大小因子（size factor），捕获测序深度和细胞大小变异，服从对数正态分布 $\ell_n \mid s_n \sim \text{LogNormal}(\ell_\mu^\top s_n, \ell_\sigma^{2\top} s_n)$。
• $\beta_{nt}$：蛋白质特异性背景强度（background intensity），服从 $\beta_{nt} \mid s_n \sim \text{LogNormal}(c_t^\top s_n, d_t^\top s_n)$。

⚪ RNA 似然函数

对于每个基因 $g$ 和细胞 $n$，TotalVI 假设观测计数 $x_{ng}$ 服从负二项分布（Negative Binomial）：

其中 $\rho_{ng} = f_\rho(z_n, s_n)$ 是通过解码器神经网络生成的归一化基因频率，$\theta_g$ 是基因特异性逆离散参数（inverse dispersion）。该建模策略与 scVI 类似，通过 $\ell_n$ 吸收测序深度变异，$z_n$ 捕获生物学变异。

⚪ 蛋白质似然函数

为建模蛋白质数据中的背景与前景信号，TotalVI 使用负二项混合分布（Negative Binomial Mixture）：

其中 $\pi_{nt} = h_\pi(z_n, s_n)$ 是混合权重，表示观测计数来自背景成分的概率；$\alpha_{nt} = g_\alpha(z_n, s_n)$ 是前景偏移参数，确保前景成分的均值 $\beta_{nt}\alpha_{nt}$ 始终大于背景成分 $\beta_{nt}$；通过引入隐变量 $v_{nt} \sim \text{Bernoulli}(\pi_{nt})$ 指示计数来源（$v_{nt}=1$ 为背景，$v_{nt}=0$ 为前景）。

该混合模型允许每个细胞-蛋白质对具有不同的背景概率，从而自适应地解耦技术噪声与真实生物学信号。

2.2 变分推断与优化

由于后验分布难以直接计算，TotalVI 使用变分自编码器（VAE）框架进行近似推断。编码器（encoder）神经网络将输入 $(x_n, y_n, s_n)$ 映射到潜在变量的近似后验分布参数：

模型通过最大化证据下界（ELBO, Evidence Lower Bound）来联合优化生成参数 $\nu$ 和变分参数 $\eta$：

训练采用随机梯度下降和 Adam 优化器，并使用 KL 退火（warm-up）策略避免陷入局部最优。

2.3 缺失数据与数据整合

TotalVI 通过 totalVI-union 模式处理不同抗体面板的数据整合。对于缺失的蛋白质测量，模型在编码时使用零填充，在解码时仅计算观测特征的 ELBO 项。这使得模型能够在保留所有蛋白质信息的同时，利用共享的潜在空间 $z_n$ 进行跨批次整合。此外，模型可以通过对缺失蛋白质的后验预测分布进行采样来实现插补（imputation）：

2.4 下游分析任务

• 去噪（Denoising）：通过将背景成分设为零，TotalVI 可以生成去噪后的蛋白质表达估计 $(1-\pi_{nt})\beta_{nt}\alpha_{nt}$，以及控制测序深度后的 RNA 表达 $\rho_{ng}$。

• 差异表达分析：TotalVI 使用 Bayes Factor 框架检测差异表达特征。对于 RNA，计算对数 fold change（LFC）$\lambda_{a,b}^g = \log_2 \rho_{a} - \log_2 \rho_{b}$；对于蛋白质，计算去噪后的 LFC $\gamma_{a,b}^t$。通过比较后验样本中 $|\lambda| \geq \delta$ 与 $|\lambda| \leq \delta$（$\delta=0.2$）的比例计算 Bayes Factor：

• 原型分析（Archetypal Analysis）：由于 $z_n$ 位于单纯形上，每个维度可解释为对特定”原型”（archetype）的权重。通过解码单位向量（one-hot vectors），可以将潜在维度与特定基因和蛋白质表达程序关联，增强模型可解释性。

3. 实验分析

作者使用小鼠脾脏和淋巴结（spleen and lymph node, SLN）的 CITE-seq 数据以及公共数据集（PBMC, MALT）评估 TotalVI 的性能。小鼠数据集包含4个样本：SLN111-D1、SLN111-D2（111种抗体，两天重复）、SLN208-D1、SLN208-D2（208种抗体，其中111种为子集）。

3.1 蛋白质背景解耦

作者首先验证了 TotalVI 区分蛋白质前景与背景的能力：

• a：散点图显示 TotalVI 前景概率与 CD20 蛋白质对数计数的关系。高斯混合模型（GMM, Gaussian Mixture Model）基于全局分布设定的阈值（垂直虚线）导致大量假阴性。TotalVI 能识别出低于 GMM 阈值但具有高前景概率的 B cells（蓝色点）。
• b：UMAP 可视化显示，高前景概率细胞（>0.8，蓝色）聚集在 B cell 区域，而低前景概率细胞（<0.2，绿色）聚集在 T cell 区域。
• c：对应细胞的 Ms4a1（CD20 编码基因）RNA 表达分布验证了上述分类的准确性。
• d：对于 T cell 标记 CD28，GMM 设定的低阈值导致大量假阳性。TotalVI 能区分高前景概率的 T cells（橙色）和本应属于背景的低前景概率 B cells（绿色）。
• e-f：UMAP 投影和 RNA 表达分布确认了 CD28 在 T cells 中的特异性表达。
• g-h：原始 CD20 蛋白质分布在 B 和 T cells 间有重叠，而去噪后分布完全分离。
• i-j：类似地，CD28 去噪后仅在 T cells 中检测到信号。

3.2 数据整合与批次校正

作者评估了 TotalVI 的整合性能：

• a：未经校正的 PCA 显示 SLN111-D1（蓝色）和 SLN208-D2（黄色）批次间存在显著混杂。
• b：TotalVI-intersect（使用蛋白质交集）整合后，批次效应被消除，细胞按类型（CD8 T、B cells、Myeloid 等）聚类。
• c：TotalVI-union（使用蛋白质并集，处理缺失值）同样实现了有效的批次混合。
• d：在潜在混合指标（latent mixing metric，越低越好）和特征保留指标（feature retention metric，越高越好）上，TotalVI 优于 Scanorama 和 Seurat v3。
• e：在高维测量混合指标（measurement mixing metric）上，TotalVI 表现与竞品相当或更优。
• f：TotalVI 成功将 SLN111-D1（有蛋白质）与 SLN111-D2（蛋白质 held-out）整合，潜在空间中细胞按类型聚类而非按批次。
• g：对 CD4、CD8a、CD23、CD357 等关键标记物，插补的蛋白质表达模式（上排）与观察值（下排）高度一致。
• h：与 Seurat v3 相比，TotalVI 的插补蛋白质与观察值的 Pearson 相关性更高，尤其对高表达蛋白质。

3.3 差异表达分析

作者展示了 DE 分析结果：

• a：TotalVI 整合的 SLN-all 数据集 UMAP，包含32个精细注释的免疫亚群（如 pDCs、NKT、Ifit3-high B cells）。
• b-c：热图显示 RNA 和蛋白质标记物在各细胞类型中的差异表达模式，对角线结构验证标记物的特异性。
• d：火山图比较 TotalVI、t-test 和 Wilcoxon 检验。传统方法错误地将 I-A/I-E、IgD、CD19 等 putative negatives 鉴定为上调（假阳性），而 TotalVI 正确识别了 CD73、CD357、CD122 等 putative positives。
• e：热图显示 TotalVI 鉴定的上调蛋白质在 ICOS-high Tregs 中特异性高表达。
• f：比较 TotalVI 估计的 LFC 与原始数据计算的 LFC。对于混杂背景噪声的蛋白质（如 CD20、CD28），TotalVI 通过去除背景获得了更大的对比度（绝对 LFC 值更大）。

3.4 B细胞异质性分析

作者深入分析了 B cell 亚群：

• a：B cell 亚群的 UMAP，包括过渡性 B（Transitional B）、成熟 B（Mature B）、B1 B、边缘区 B（MZ B）和 Ifit3-high B。
• b-c：六种经典表面标记（CD93、CD24、CD21、CD43、IgD、CD23）的蛋白质和 RNA 表达空间分布高度一致。
• d-e：脾脏和淋巴结的 B cell 组成差异，B1 和 MZ B 主要分布于脾脏。
• f：蛋白质 DE 热图显示各亚群特异性标记（如 CD93 和 CD24 在过渡性 B 细胞中高表达）。
• g：RNA DE 热图揭示了非表面蛋白编码基因（如转录因子 Bhlhe41 标记 B1 B 细胞）的调控模式。
• h：过渡性 B 细胞中 RNA 和蛋白质的相关矩阵，显示两个反相关模块分别对应过渡性特征（如 CD93）和成熟特征（如 IgD）。
• i：潜在维度 $Z_{16}$ 捕获了从过渡性到成熟 B 细胞的连续分化轨迹。
• j：随着 $1-Z_{16}$ 增加（沿分化轴），过渡性模块特征表达下降，成熟模块特征表达上升，展示了协调的基因-蛋白质共调控程序。