单细胞多组学主题嵌入揭示了细胞类型特异性和与 COVID-19 严重程度相关的免疫特征.

0. TL; DR

moETM（multi-omics embedded topic model）是一种可解释的深度生成模型，用于整合单细胞多组学数据。该方法基于变分自编码器（variational autoencoder, VAE）框架，采用高斯乘积（product-of-Gaussians, PoG）作为编码器融合多模态信息，并使用线性解码器学习低维主题（topic）嵌入与特征嵌入。通过引入批次效应校正参数，moETM能够有效消除技术批次差异，同时保留生物学变异。

作者在7个公开数据集上与6种 state-of-the-art 方法（SMILE、scMM、Cobolt、MultiVI/TotalVI、MOFA+、Seurat V4）进行了系统比较，证明moETM在细胞聚类、批次校正和跨模态插补任务中均表现优异。特别地，moETM应用于 COVID-19 CITE-seq 数据集时，成功识别出与疾病严重程度相关的免疫细胞亚群特征（如B cell恶性病变相关主题），并发现IL6和CD5等潜在治疗靶点。

1. 背景介绍

近年来，单细胞测序技术已从单一模态发展到多模态并行测量。CITE-seq（Cellular Indexing of Transcriptomes and Epitopes by Sequencing）技术能够同时捕获单个细胞的转录组（transcriptome）和表面蛋白质组（surface proteome）；而 scRNA-seq 与 scATAC-seq（assay for transposase-accessible chromatin using sequencing）的联合分析则提供了基因表达与染色质可及性（chromatin accessibility）的关联视图。这些多组学数据为解析细胞异质性、基因调控程序和疾病机制提供了前所未有的机会。

然而，单细胞多组学数据分析面临以下挑战：

目前已发展出多种单细胞多组学整合方法。这些方法在特定任务上表现良好，但往往需要牺牲可扩展性、可解释性或灵活性。尤其当使用神经网络编码高维多组学数据时，其黑盒特性使得识别细胞类型特异性标志物和调控程序变得困难。因此，迫切需要一种既能有效整合多模态数据、校正批次效应，又能提供生物学可解释性的统一框架。

2. moETM 方法

moETM 将单细胞多组学数据类比为多语言文档：每个细胞是一个文档，每种组学模态（基因、蛋白质、染色质峰）是一种语言，测序读数（reads）是文档中的词符（tokens）。模型假设每个细胞的多模态观测数据是由一组潜在的主题（topics）混合生成的，这些主题代表细胞状态或生物学程序。

2.1 数据生成过程

对于包含 $M$ 种模态的数据，细胞 $n$ 的观测为计数向量集合 \(\{\mathbf{x}_n^{(m)}\}_{m=1}^M\) ，其中每种模态 $m$ 有 $V^{(m)}$ 个特征（如基因、蛋白质或峰）。

每个细胞 $n$ 的主题混合比例 $\boldsymbol{\theta}_n \in \mathbb{R}^K$（$K$ 为主题数）服从Logistic Normal分布：

\[\boldsymbol{\delta}_n \sim \mathcal{N}(0, \mathbf{I}), \quad \boldsymbol{\theta}_n = \text{softmax}(\boldsymbol{\delta}_n)\]

对于模态 $m$ 中的第 $i$ 个读数，其对应的特征索引 $v^{(m)}$ 服从分类分布：

\[\mathbf{w}_{n,i}^{(m)} \sim \text{Cat}(\mathbf{r}_{n}^{(m)})\]

其中 $\mathbf{r}_{n}^{(m)} \in \mathbb{R}^{V^{(m)}}$ 是观测速率向量，参数化为：

\[\hat{r}_{n,v}^{(m)} = \boldsymbol{\theta}_n \mathbf{\alpha} \mathbf{p}_{:,v}^{(m)} + \lambda_{s(n),v}^{(m)}\]

其中$\mathbf{\alpha} \in \mathbb{R}^{K \times L}$ 是主题嵌入矩阵（$L$ 为嵌入维度）；$\mathbf{p}^{(m)} \in \mathbb{R}^{L \times V^{(m)}}$ 是模态 $m$ 的特征嵌入矩阵;$\lambda_{s(n),v}^{(m)}$ 是依赖于批次 $s(n)$ 的批次效应校正参数。

归一化后的期望观测速率为：

\[r_{n,v}^{(m)} = \frac{\exp(\hat{r}_{n,v}^{(m)})}{\sum_{v'=1}^{V^{(m)}} \exp(\hat{r}_{n,v'}^{(m)})}\]

给定观测计数 \(\mathbf{x}_n^{(m)}\) （特征 $v$ 的计数为 $x_{n,v}^{(m)}$ ），似然服从多项分布（Multinomial）：

\[p\left(\{\mathbf{x}_n^{(m)}\}_{m=1}^M \mid \mathbf{r}_n\right) \propto \prod_{m=1}^M \prod_{v=1}^{V^{(m)}} \left(r_{n,v}^{(m)}\right)^{x_{n,v}^{(m)}}\]

2.2 编码器：高斯乘积（Product-of-Gaussians）

为整合多模态信息，moETM 采用 PoG 框架替代传统的混合专家模型（MoE）。对于每种模态 $m$，编码器神经网络（两层全连接层）输出该模态特有的高斯参数：

\[[\boldsymbol{\mu}_n^{(m)}; \log \boldsymbol{\sigma}_n^{(m)}] = \text{Encoder}(\tilde{\mathbf{x}}_n^{(m)}; \mathbf{W})\]

其中 $\tilde{\mathbf{x}}_n^{(m)}$ 是模态 $m$ 的归一化计数（每个特征计数除以该模态总计数）。

PoG 融合各模态的高斯分布，得到联合后验分布 \(q(\boldsymbol{\delta}_n) = \mathcal{N}(\boldsymbol{\mu}^*, \boldsymbol{\sigma}^{*2})\) ，其参数为：

\[\boldsymbol{\mu}^* = \frac{\sum_{m=1}^M \boldsymbol{\mu}_n^{(m)} / \boldsymbol{\sigma}_n^{(m)2}}{1 + \sum_{m=1}^M 1/\boldsymbol{\sigma}_n^{(m)2}}, \quad \boldsymbol{\sigma}^{*2} = \frac{1}{1 + \sum_{m=1}^M 1/\boldsymbol{\sigma}_n^{(m)2}}\]

与 MoE 方法（对每个模态采样后平均）相比，PoG 只需从联合高斯分布采样一次，减少了蒙特卡洛近似带来的采样噪声，提高了变异推断的稳定性。

2.3 解码器与批次校正

解码器采用线性矩阵分解结构，重建各模态的观测：

\[\hat{\mathbf{x}}_n^{(m)} \propto \boldsymbol{\theta}_n \mathbf{\alpha} \mathbf{p}^{(m)} + \boldsymbol{\lambda}_{s(n)}^{(m)}\]

这种设计将细胞-主题矩阵 $\boldsymbol{\Theta}$、主题嵌入 $\mathbf{\alpha}$ 和特征嵌入 $\mathbf{p}^{(m)}$ 分离，使得不同模态共享相同的细胞-主题表示，确保跨模态一致性；各模态拥有独立的特征嵌入，捕捉模态特异性；批次效应通过加性参数 $\boldsymbol{\lambda}_{s}^{(m)}$ 进行校正，允许模型学习跨批次的共享生物学主题。

2.4 变分推断与模型训练

由于边缘似然难以直接计算，moETM 通过最大化证据下界（ELBO）进行优化：

\[\mathcal{L}_n = \mathbb{E}_{q(\boldsymbol{\delta}_n)}[\log p(\{\mathbf{x}_n^{(m)}\}_{m=1}^M \mid \boldsymbol{\delta}_n, \boldsymbol{\Theta})] - \text{KL}[q(\boldsymbol{\delta}_n) \| p(\boldsymbol{\delta}_n)]\]

其中 KL 散度有闭式解。通过重参数化技巧（reparameterization trick）采样 $\tilde{\boldsymbol{\delta}}_n = \boldsymbol{\mu}^* + \boldsymbol{\sigma}^* \odot \boldsymbol{\epsilon}$（$\boldsymbol{\epsilon} \sim \mathcal{N}(0, \mathbf{I})$），模型参数（编码器权重 $\mathbf{W}$、解码器权重 $\mathbf{\alpha}$、$\mathbf{p}^{(m)}$ 和批次参数 $\boldsymbol{\lambda}^{(m)}$）可通过反向传播端到端优化。

2.5 跨模态插补与下游分析

训练后的模型可利用已观测模态推断缺失模态。例如，利用 RNA 编码器生成主题比例，再通过蛋白质解码器重建蛋白质表达（RNA2protein），反之亦然（protein2RNA）。

模型推断的细胞主题比例可用于聚类和可视化。每个主题的特征权重 $\mathbf{\alpha}\mathbf{p}^{(m)}$ 可识别该主题的标记基因、蛋白质或调控峰。结合 GSEA（Gene Set Enrichment Analysis）和 motif 富集分析，可揭示细胞类型特异性通路和转录因子调控网络。

3. 实验分析

作者在7个公开数据集上评估了 moETM，包括4个 gene+peak 数据集（BMMC1、MSLAC、MKC、MBC）和3个 gene+protein 数据集（BMMC2、HWBC、HBIC）。对比方法包括 SMILE、scMM、Cobolt、MultiVI/TotalVI、MOFA+ 和 Seurat V4。评估指标涵盖生物学保守性（ARI、NMI）和批次校正（kBET、GC）。

3.1 多组学整合性能评估

作者展示了 moETM 与对比方法在多种实验设置下的性能（60/40随机分割、整数据集训练、跨批次训练）：

3.2 细胞聚类与批次校正可视化

作者通过 UMAP 可视化验证聚类质量：

3.3 跨模态插补准确性

作者评估了 moETM 的跨模态插补能力：

3.4 细胞类型特征识别与调控网络

作者展示了 BMMC2 数据的主题分析：

作者分析了 BMMC1（gene+peak）数据的调控关系：

3.5 COVID-19疾病严重程度分析

作者展示了 moETM 在 COVID-19 CITE-seq 数据集（HBIC，130个患者，781,123个细胞）上的应用：

这些结果表明 moETM 不仅能识别已知的细胞类型标记，还能揭示与疾病表型相关的稀有细胞亚群（如恶性 B cell）和潜在治疗靶点（IL6、CD5）。