Mowgli方法细节与可解释性实验设置

tzc_fly

已于 2024-09-23 12:06:29 修改

阅读量755

点赞数 29

分类专栏：生物计算工具单细胞多组学分析文章标签：数据挖掘深度学习

于 2024-09-23 11:25:46 首次发布

本文链接：https://blog.csdn.net/qq_40943760/article/details/142442920

版权

生物计算工具同时被 2 个专栏收录

52 篇文章 69 订阅

订阅专栏

单细胞多组学分析

50 篇文章 43 订阅

订阅专栏

Mowgli方法细节

在Mowgli中，结合了NMF和OT，将多个组学矩阵分解为组学特定的字典和共享的嵌入。OT被用作无监督重建的损失函数，Mowgli的目标是无监督学习，并且提高聚类性能。

在Mowgli之前，OT已经被用于提高聚类性能，主要是单模态数据集下的聚类：Optimal transport improves cell–cell similarity inference in single-cell omics data（动机：细胞异质性可以用无监督聚类来实现，而无监督聚类依赖于相似性的度量，所以选择用OT作为细胞相似性的度量）
fig0

OT在scRNA-seq分析中提高聚类性能。从输入预处理数据到性能评估所采用的流程分别总结为（A）基线metric 和（B）OT metric。在scanpy中，相当于是sc.pp.neighbors中添加了OT以优化细胞-细胞图。

符号

考虑 $n$ 个细胞，跨模态测量，每个模态 $p$ 的特征数量是 $m_{p}$ （比如基因的数量）。记 $\textbf{A}^{(p)}\in\R_{+}^{m_{p}\times n}$ 为模态 $p$ 的数据集。此外，要求 $\textbf{A}^{(p)}$ 的每一列之和为1，即离散概率分布（每个细胞的基因表达量总和为1）。

最优传输

最佳传输 (OT) 由 Monge 和 Kantorovich 定义，旨在通过计算将一个分布传输到另一个分布的最小成本来比较两个概率分布 A 和 B。在离散的情况下，经典的OT距离，也称为Wasserstein distance，在 $m$ 维的直方图 $\textbf{a}=(a_{1},...,a_{m})$ 和 $\textbf{b}=(b_{1},...,b_{m})$ 中被定义为：
eq1
其中：
eq1-1
在离散场景下，coupling $\textbf{P}\in\prod(\textbf{a},\textbf{b})$ 是一个矩阵，表示离散概率分布 $\textbf{a}$ 中的mass（资源量）如何从一个bin（例如基因）移动到另一个bin，以将 $\textbf{a}$ 转换为 $\textbf{b}$ 。换句话说， $P_{k,l}$ 是在将细胞 $\textbf{a}$ 的基因表达谱运输到细胞 $\textbf{b}$ 的基因表达谱时，基因 $k$ 和基因 $l$ 之间运输的基因表达量。

ground cost $\textbf{C}\in\R_{+}^{m\times m}$ 是一个pairwise distance 矩阵，对将mass（资源量）从一个特征（比如基因）转移到另一个特征的惩罚进行编码。因此， $\textbf{C}$ 的选择应该使相似的bins（gene） $k$ 和 $l$ 具有较低的cost $C_{k,l}$ 。事实上，这将有利于相似基因之间的基因表达传输。对于一个具体的组学 $p$ ，我们以数据驱动的方式将 $\textbf{C}$ 定义为特征之间的成对余弦距离矩阵，特征即数据集 $\textbf{A}^{(p)}$ 中的行。换言之， $\textbf{u}_{k}\in\R_{+}^{n}$ 和 $\textbf{u}_{l}\in\R_{+}^{n}$ 为数据集中的两行：
eq2
其中 $<\textbf{x},\textbf{y}>$ 为dot product。以这个方法计算ground cost已经在Optimal transport improves cell–cell similarity inference in single-cell omics data中得到验证。

由于单细胞数据的高维性，作者使用经典 OT 的近似值。下面公式中描述的 OT 熵正则化是使用快速且支持 GPU 的 Sinkhorn 算法计算的。
eq3
其中，熵 $\textbf{E}$ 被定义为：
eq4
如果 $\epsilon$ 被设为0，OT退化为经典OT。如果值越大，则coupling $\textbf{P}$ 越分散。Optimal transport improves cell–cell similarity inference in single-cell omics data已经证明，这种正则化的OT可以提高细胞-细胞相似性。熵正则化有望控制由于技术误差和单细胞水平基因表达的随机性而引起的系统噪声。此外，更分散的耦合增加了特征之间的mass交换。这使 OT 能够利用特征（例如基因）之间的关系，从而进一步激发其在单细胞数据中的应用。

Mowgli

作者旨在把矩阵 $\textbf{A}^{(p)}$ 分解为矩阵 $\textbf{H}^{(p)}\in\R_{+}^{m_{p}\times k}$ （模态特定字典）和 $\textbf{W}\in\R_{+}^{k\times n}$ （跨模态共享的cell embeddings）。整数 $k$ 是latent space的维数。作者使用正则化的OT作为重建损失用于比较 $\textbf{H}^{(p)}\textbf{W}$ 和 $\textbf{A}^{(p)}$ 。

最终损失为：
eq5
其中， $j$ 代表细胞 $j$ ，参数 $\rho_{p}$ 和 $\mu$ 控制 $\textbf{H}^{(p)}$ 和 $\textbf{W}$ 的稀疏性。参数的具体取值参考mowgli原文的methods部分。

其次，在Liu数据集上，因子数（隐空间维数）为5，其余数据集为50。

作者交替最小化 $\textbf{H}^{(p)}$ 和 $\textbf{W}$ 。可以证明， $\textbf{H}^{(p)}$ 和 $\textbf{W}$ 上的这些光滑最小化问题等价于新对偶变量 $\textbf{G}^{(p)}$ 上的以下光滑最小化问题。这些问题可以使用标准优化方法解决，首选的方法是 L-BFGS，一种有限内存拟牛顿方法。

数据预处理方法

RNA预处理

细胞的质量控制过滤是根据线粒体基因表达的比例、表达基因的数量和计数总数进行的（使用 Muon 的 filter_obs）。基因的质量控制过滤是根据表达基因的细胞数量进行的（使用 Muon 的 filter_var）。将细胞标准化为总和为 10000（使用 Scanpy 的 normalize_total），然后进行对数转换（使用 Scanpy 的 log1p）。选择变异性最大的前 2500 个基因（对于 Liu 数据集为 1500 个）进行下游分析（使用 Scanpy 的 highly_variable_genes，其中 flavor=‘seurat’）。

ATAC预处理

对开放峰（open peaks）的数量和计数总数进行细胞质量控制过滤（使用 Muon 的 filter_obs）。对峰开放的细胞数量进行峰质量控制过滤（使用 Muon 的 filter_var）。在 Liu、TEA 和 10X PBMC 中，将细胞标准化为总和为 10000（使用 Scanpy 的 normalize_total），然后进行对数转换（使用 Scanpy 的 log1p）。在 OP Multiome 中，使用 TF-IDF（使用 Muon 的 tfidf）对细胞进行标准化，以遵循其作者选择的预处理。选择变化最大的峰进行下游分析（使用 Scanpy 的 highly_variable_genes，其中 flavor=‘seurat’）。由于数据在数据集中的分布存在差异，我们选择在 Liu 中保留 1500 个峰，在 PBMC 中保留 5,000 个峰，在 OP Multiome 和 TEA 中保留 15,000 个峰。

ADT预处理

由于蛋白质数量较少，数据噪声比 RNA 或 ATAC 低，因此未进行质量控制或特征选择。数据通过 Center Log Ratio 进行归一化（使用 Muon 的 clr）。

数据分析

基因集富集分析-GSEA

gProfiler API 通过 Scanpy 的 rich 功能使用。自定义源 GO:CC、GO: MF、GO: BP、Azimuth 和 ImmuneSigDB 均从 Enrichr 网站检索。选择调整后的 p 值低于 0.05（经过 Bonferroni 校正）的基因集进行进一步分析。为了使基因具有可比性，将矩阵 $H^{(rna)}$ 的行标准化为 1。然后，每个因子的前 150 个基因被用作 gProfiler 的无序输入。

Motif富集分析

Signac 用于执行 Motif 富集分析，使用 JASPAR2022 基序数据库。为了使峰值具有可比性，我们将矩阵 $H^{(atac)}$ 的行标准化为 1。每个因子的 100 个最高峰值被用作 Signac 的 FindMotifs 的输入。跨因子最高峰值的并集构成背景。

然而，如何搜索重要因子，比如细胞类型特异性因子，然后，如何找到因子中的top特征？

可视化

为了可视化 MOFA+、integrative NMF 和 Mowgli 模型中细胞的表征，作者计算了 kNN 图 (k = 20)，并使用细胞低维嵌入之间的欧几里得距离（使用 Scanpy 的邻居）。使用这些图计算 2-D UMAP（使用 Scanpy 的 umap）。对于 Seurat v4，使用 Seurat v4 的 RunUMAP 函数执行基于 WNN 图的 2-D UMAP 投影。

聚类

对于 Mowgli、integrative NMF 和 MOFA+，作者使用具有不同分辨率的 Leiden 算法对数据集进行聚类（使用 Scanpy 的 leiden）。与 UMAP 可视化类似，Leiden 算法的输入是之前计算的 kNN 图。对于 Seurat v4，使用 Seurat v4 的 FindClusters 函数执行 Leiden 聚类。

TEA-seq上的注释：手动注释

MOFA+ 使用了 15 个因子，integrative NMF 使用了 30 个因子。Mowgli 使用了 50 个因子。对于这三种方法，都应用了粗 Leiden 聚类（使用 Scanpy 的 leiden）。对于这三种方法，每个簇都根据典型marker 基因和蛋白的表达分配一种细胞类型。为了确认此注释，在数据集的 RNA 信号上运行了 Azimuth（使用 Azimuth 网络工具和 PBMC reference）。在 Sankey 图中确认了四个独立注释的一致性。由于聚类粗糙，手动注释中没有树突状细胞。同样，ADT 信号告诉我们，所有四个注释中都遗漏了一个 CD4-CD8-T 细胞群。

可解释实验1：TEA-seq上的细胞类型特异性factors

为了测试 Mowgli、NMF 和 MOFA+ 的生物学可解释性，作者评估了它们的因子与已注释的免疫细胞类型之间的关联。在此做出的基本假设是，可解释的方法应该提供并非在所有细胞中广泛活跃的因子，而是选择性地与细胞类型相关的因子（或者其他表型，比如临床标签）。事实上，表征由多种因子组合产生的细胞类型是一项艰巨的任务。相反，拥有细胞类型特异性因子使相关细胞类型的生物学解释变得简单。

为了评估这种特异性，对于每种细胞类型，作者绘制了 Mowgli、NMF 和 MOFA+ 的因子如何根据它们在细胞类型内的平均权重和它们在细胞类型外的平均权重分布（图 4C）。特定于细胞类型的因子在细胞类型内的平均权重应较高，在细胞类型外的平均权重应较低，因此落在图的左上角。由于 MOFA+ 的因子不限于正向，其正向和负向部分可能与不同的生物学信息相关，作者将每个因子分为两部分，就像 MOFA+ 的解释工具中所做的那样。此外，作者用特异性分数量化了每个因子的表现，该分数也在图 4C 中以粗体显示，并在方法部分中定义（Methods-Specificity小节）。
fig4c

如图 4C 所示，MOFA+ 和 NMF 倾向于将多个因子与同一种细胞类型相关联，Mowgli 定义的因子与细胞类型之间存在明确的一对一关联。此外，这些因子在 Mowgli 中的特异性得分高于在 MOFA+ 和整合性 NMF 中的特异性得分。这在 NK 细胞、CD8 T 细胞和 CD4 T 细胞中尤其明显，其中 MOFA+ 和 NMF 似乎都聚合了来自许多因子的信息，而 Mowgli 更具选择性。

气泡图的绘制可参考mowgli_reproducibility/visualize/visualize_tea.py文件下的plot_bubbles。

当选择指定的因子后，我们可以获得该因子下的 top 特征，输入为指定模态的特征embedding，比如 $\textbf{H}^{(rna)}\in\R_{+}^{m_{rna}\times k}$ ，对应参数H_mowgli，dim为指定的因子位置（是一个整数，范围 $0$ 到 $k - 1$ ），比如15：

def top_mowgli(dim, n, H_mowgli):
    """
    Get the top n feature for a given dimension.
    """
    H_scaled = H_mowgli / H_mowgli.sum(axis=1, keepdims=True)
    return H_scaled[:, dim].argsort()[::-1][:n]

当选出top genes或者peaks后，Mowgli的教程提供了富集分析R示例：https://mowgli.readthedocs.io/en/latest/vignettes/TF%20enrichment.html

富集结果可以解释与特定细胞类型相关的通路，Mowgli用这个功能演示了免疫细胞亚群相关的生物信息发现。作者重点关注四种免疫细胞亚群相关的因子：effector memory CD8 T-cells (因子 49)， naive B cells (因子 33)， memory B cells (因子 44)， CD56dim NK cells (因子 2)。对于每个因子，考虑了三个模态的embeddings，其中，对 $\textbf{H}^{(rna)}$ 进行基因富集分析，对 $\textbf{H}^{(atac)}$ 进行motifs富集分析，图5C展示了结果。

下面以effector memory CD8 T-cells (因子 49)为例说明Mowgli发现的生物学信息。对应于因子 49 的效应记忆 CD8 T 细胞 (CD8 TEM 细胞)，Mowgli 提取出两个top基因 (CRTAM 和 KLRK1)，已知这两个基因对 CD8+T 细胞介导的细胞毒性至关重要，两个top蛋白质 (CD45RO、TCR-a/b) 分别是已知的记忆 T 细胞marker和 T 细胞受体。更有趣的是，这里还确定了该亚群的几种转录因子 (TF) 候选调节子，其中包括 EOMES 和 TBX21 (又名 T-bet)，已知它们对 CD8 TEM 发育很重要。此外，五种top候选 TF 调节子 (TBR1、TBX21、TBX4、TBX5 和 MGA) 以同一因子下的三个top基因为target genes (CCL5、CRTAM 和 IL21R)，从而表明调控过程可能对 CD8 TEM 细胞很重要。

fig5c