用于单细胞多组学整合的无监督拓扑对齐方法

摘要

动机：单细胞多组学数据提供了细胞的全面分子视图（molecular view of cells）。然而，单细胞多组学数据集通常由未配对的细胞组成，这些细胞在不同模态（组学）下具有明显的不匹配特征，这使得数据整合具有挑战性（两大困难：unpaired cells和unmatched features）。

结果：在本研究中，提出了一种称为UnionCom的新算法，用于单细胞多组学数据整合。UnionCom不需要任何监督信息，无论是细胞之间的配对标记还是特征之间的调控关系。它不需要跨数据集的细胞之间的一对一对应关系，并且可以保留特定于数据集细胞类型的样本。UnionCom在模拟和真实单细胞多组学数据集上都优于最先进的方法。UnionCom对参数选择以及特征的子集采样具有鲁棒性。

1.单细胞多组学数据整合背景

单细胞测序的出现提供了细胞水平的高分辨率组学图谱，例如基因组学的单细胞DNA测序、转录组的单细胞RNA测序（sRNA-seq）、染色质可及性的单细胞测序（ATAC），为揭示与基本生物学问题相关的分子机制提供了机会，例如细胞的分化。虽然我们已经对单细胞基因组、转录组和表观遗传学数据的计算分析进行了广泛的研究，但大多数方法仅限于处理某个特定组学的单细胞数据，从而造成了一系列关于生物过程研究的分隔。最近，从同一生物样本或组织中采样的细胞中提取的跨模态单细胞多组学数据正在出现。整合单细胞多组学数据将建立多模态的联系，提供更全面的内在生物过程的分子视图。

单细胞数据集的整合正引起机器学习和数据科学领域的关注。已经开发了许多单细胞数据整合方法，用于在一种模态下对单细胞数据集进行批次效应校正（batch-effect correction）。然而，与单细胞批次效应校正问题相比，跨模态整合单细胞多组学数据在两方面提出了新的挑战。首先，单细胞多组学数据集由未配对细胞组成，这些细胞在不同模态下具有不匹配的特征。由于大多数单细胞测序分析对细胞具有破坏性，相同模态或不同模态的单细胞数据集通常具有未配对的细胞。此外，单细胞多组学数据集在不同模态下没有共同特征，因为每个模态都旨在从特定方面描述细胞。不同模态的不同特征存在于不同的数学空间中。因此，一般来说，单细胞多组学数据在样本（细胞）之间或特征之间没有任何的直接对应关系。在实践中，在应用最先进的单细胞数据整合方法之前，通常会根据先验知识将不同的特征在不同的模态之间进行经验上的预匹配（比如指定调控关系），再形成一个公共空间。其次，单细胞多组学数据是通过不同的测量分析方法生成的，具有不同的基础数据生成机制。因此，跨模态的数据集上的复杂非线性几何畸变可以引入固有的低维结构。这种复杂的畸变使得单细胞多组学数据的整合比一种模态中单细胞数据集的批次效应校正困难得多。因此，对于单模态的批次效应校正，线性运算（如平移和缩放）通常工作良好，但是不能对齐具有复杂非线性几何畸变（geometrical distortions）的单细胞多组学数据集。

相关工作

批次效应校正

目前已经开发了用于单细胞数据整合的计算方法。然而，许多现有方法是为批次效应校正而设计的。单细胞多组学整合的挑战是由于不配对的细胞，不匹配的特征以及跨模态的复杂非线性几何畸变而产生的，并且这些问题仍然没有解决。开创性的单细胞数据整合方法，相互最近邻法（MNN，mutual nearest neighbor），旨在通过对齐MNN消除数据集之间的批次效应，MNN是通过公共特征空间中跨数据集的细胞之间的欧几里德距离计算的。然而，MNN仅限于批次效应变化比不同细胞类型之间的生物效应变化小得多的情况。已建立的单细胞数据分析管道Seurat使用典型相关分析（CCA，canonical correlation analysis）将数据集上的特征空间投影到公共子空间中，该公共子空间最大化了数据集间的相关性，然后采用类似于MNN的策略，通过基于从公共子空间计算的MNN在数据集之间找到anchor cell（即对应点）来对齐数据集之间的细胞。然而，由于CCA是一种线性降维方法，它无法捕捉跨模态的单细胞多组学数据之间的非线性相互关联关系。Scanorama建立在基于MNN的策略之上，用于高效地整合大规模scRNA-seq数据集。scAlign开发了一种基于自动编码器的方法，用于整合scRNA-seq数据集（公共特征空间的非线性降维），跨数据集对齐细胞。Harmony使用主成分分析（PCA）将细胞投影到共享embedding空间中。

单细胞多组学整合方法

为跨模态的单细胞多组学整合开发的方法正在出现。Liger（Welch等人2019年提出）研究了具有跨模态预匹配公共特征空间的数据集，采用非负矩阵分解找到公共特征的共享低维因子，以匹配单细胞多组学数据集。MATCHER（Welch等人2017年提出）取消了基于流形对齐的单细胞多组学数据集的特征对应信息要求。然而，MATCHER仅限于一维轨迹的对齐，它无法对齐复杂的轨迹，如分支树结构。MAGAN（Amodio和Krishnaswamy在2018年提出）通过生成对抗网络方法对齐生物流形，整合了单细胞多组学数据集。然而，当跨数据集的样本之间的对应信息不可用时，它在无监督任务中缺乏效果。Harmonic（Stanley等人2018年提出）基于扩散图方法在数据集上对齐低维几何结构，但它需要部分特征的匹配信息。最大平均差异流形对齐（MMD-MA；Liu等人2019年提出），一种用于单细胞多组学数据集的无监督流形对齐算法，将潜在的生物低维结构嵌入到再生核希尔伯特空间（RKHSs）中，并通过最小化不同模态间的MMD，找到用于流形对齐的公共子空间。

2.Materials and Methods

2.1.UnionCom algorithm

UnionCom是一种用于单细胞多组学数据整合的无监督拓扑对齐算法（见图1）。我们对UnionCom的详细信息描述如下。这里，我们为两个数据集的情况制定了我们的方法。然而，对于任何的单细胞多组学数据集，它都可以很容易地应用。假设我们有两个单细胞多组学数据集$\text{X}=[x_1,...,x_{n_x}]\in R^{d_x\times n_x}$和$\text{Y}=[y_1,...,y_{n_y}]\in R^{d_y\times n_y}$。对于$\text{X}$，$d_x$是特征的维数，$n_x$是细胞的数量，特征可以是不同组学下的（比如gene expression，DNA methylation）。在不丧失一般性的情况下，我们假设$n_x\le n_y$，我们假设$\text{X}$和$\text{Y}$的固有低维流形是$M_x$和$M_y$，并且它们具有相似的拓扑结构。给定$\text{X}$和$\text{Y}$，UnionCom由以下三个主要步骤组成（A1-A3）。

• 图1：UnionCom的概况图。
• a：给定两个单细胞不同组学的数据集（Dataset1和Dataset2），两个数据集具有相似的embedding拓扑结构。
• b：UnionCom将每个单细胞数据集转换为细胞之间的距离矩阵。
• c：通过全局缩放参数$\alpha$（global scaling parameter）重新缩放跨数据集的拓扑结构上的全局失真（global distortions）。
• d：通过矩阵优化方法匹配距离矩阵，跨数据集对齐细胞；
• e：最后，将这些不匹配特征跨模态投影到公共的embedding空间中，达到细胞的特征可比性。
• UnionCom不需要在跨数据集的细胞之间标注对应关系，并且可以保留具有数据集特定的样本（比如Dataset2中的黑色分支）

A1.从单细胞数据集到距离矩阵

UnionCom通过使用距离矩阵将每个单细胞数据集的固有低维结构embedding到同一度量空间，距离矩阵被定义为$[\text{K}{]}_{ij}=d(x_i,x_j)$，其中，$d$是流形上细胞之间的测地距离（geodesic distance）。因此，${\text{K}}_x$和${\text{K}}_y$代表$M_x$和$M_y$的测地距离矩阵。为了计算测地距离，UnionCom首先基于欧几里得距离为每个数据集构建一个$k$最近邻graph（$k$为使得k最近邻graph连通的最小数），然后使用Floyd-Warshall算法计算图上每对节点（每对细胞）之间的最短距离，将这个最短距离近似为测地距离（Floyd算法解决的是图中任意两点的最短路径问题）。

---

测地距离geodesic distance：比如在三维Mesh中，Geodesic Distance 就是两顶点沿网格表面的最短路径的距离。

---

A2.匹配距离矩阵从而对齐细胞

通过扩展无监督流形对齐算法GUMA（Nips2014），UnionCom提出全局缩放参数$\alpha$，用于重新缩放数据集拓扑的全局结构，另外，UnionCom放宽了GUMA要求的跨数据集细胞与细胞对应的限制。并且UnionCom能够保留特定于数据集的细胞样本。

具体而言，UnionCom将无监督拓扑对齐问题转为矩阵优化问题：$$mi{n}_{\text{F},\alpha }||\alpha 1_{n_x\times n_y}\odot {\text{K}}_x-\text{F}{\text{K}}_y{\text{F}}^T|{|}_F^2$$其中，$\text{F}=[{\text{F}}_{ij}{]}_{n_x\times n_y}$是点匹配矩阵（point matching matrix）用于满足某些约束需要，$1_n\in R^n$代表全1向量，$1_{m\times n}\in R^{m\times n}$代表全1矩阵，$\odot$代表哈达姆乘积。$||X|{|}_F^2=tr(X^{T}X)$代表矩阵的Frobenius范数。

在GUMA中， F ∈ { 0 , 1 } n x × n y \textbf{F}\in\left\{0,1\right\}^{n_{x}\times n_{y}} F∈{0,1}nx​×ny​是硬编码的匹配矩阵，用于标记$\text{X}$和$\text{Y}$的对应关系，其中，$[{\text{F}}_{ij}]=1$表示$x_i$和$y_j$是对等的或者匹配的（counterpart）。因此，其实GUMA隐式地假设两个数据集之间需要存在一些细胞具有强烈一对一的匹配关系，这一现象通常与实际数据不一致。GUMA在所有可能的0-1整数矩阵集合（该集合记为$Ⅱ$）上找合适的$\text{F}$。 Ⅱ = { F ∣ F ∈ { 0 , 1 } n x × n y , F 1 n y = 1 n x , F T 1 n x ≤ 1 n y , n x ≤ n y } Ⅱ=\left\{\textbf{F}|\textbf{F}\in\left\{0,1\right\}^{n_{x}\times n_{y}},\textbf{F}1_{n_{y}}=1_{n_{x}},\textbf{F}^{T}1_{n_{x}}\leq 1_{n_{y}},n_{x}\leq n_{y}\right\} Ⅱ={F∣F∈{0,1}nx​×ny​,F1ny​​=1nx​​,FT1nx​​≤1ny​​,nx​≤ny​}UnionCom寻求数据集之间细胞的soft匹配，它放宽了之前$\text{F}$的约束，转为：$${Ⅱ}^{\prime }=\left\{\text{F}|\text{F}\ge 0,\text{F}{1}_{n_y}=1_{n_x},{\text{F}}^T{1}_{n_x}\le 1_{n_y},n_x\le n_y\right\}$$其中，$\text{F}\ge 0$表示矩阵中每个元素都是非零的，并且每行的总和为1。此时$\text{F}\in {Ⅱ}^{\prime }$提供了两个数据集之间样本匹配的概率解释：$[\text{F}{]}_{ij}$表示$x_i$和$y_j$的匹配概率。

因此，UnionCom通过解决约束集上的矩阵优化问题来调整对齐两个数据集之间的拓扑结构，并且我们可以发现，当数据集之间的细胞匹配信息可用时，我们可以手动将$\text{F}$的对应元素设置为1和0，UnionCom就能扩展到有监督或者半监督数据整合。

A3.将多组学数据投影到公共embedding空间

为了同时保留特定数据集固有的局部分支，UnionCom建立在t-SNE方法之上。具体的，UnionCom旨在embed $\text{X}\in R^{d_x\times n_x}$和$\text{Y}\in R^{d_y\times n_y}$到一个公共的embedding空间（$p$维，$p\le min(d_x,d_y)$），降维数据集表示为${\text{X}}^{\prime }\in R^{p\times n_x}$和${\text{Y}}^{\prime }\in R^{p\times n_y}$。UnionCom通过最小化$L({\text{X}}^{\prime },{\text{Y}}^{\prime })$得到降维数据集：$$mi{n}_{{\text{X}}^{\prime },{\text{Y}}^{\prime }}L({\text{X}}^{\prime },{\text{Y}}^{\prime })=KL({\text{P}}^{\text{X}}||{\text{Q}}^{{\text{X}}^{\prime }})+KL({\text{P}}^{\text{Y}}||{\text{Q}}^{{\text{Y}}^{\prime }})+\beta ||{\text{X}}^{\prime }-{\text{Y}}^{\prime }{\text{F}}^T|{|}_F^2$$其中，$\text{F}$是A2获得的点匹配矩阵，${\text{P}}^{\text{X}}=[{\text{P}}_{ij}^{\text{X}}{]}_{n_x\times n_x}$是$\text{X}$原始空间中定义的细胞到细胞的转移概率矩阵，${\text{P}}^{\text{Y}}=[{\text{P}}_{ij}^{\text{Y}}{]}_{n_y\times n_y}$是$\text{Y}$原始空间中定义的细胞到细胞的转移概率矩阵。${\text{Q}}^{{\text{X}}^{\prime }}=[{\text{Q}}_{ij}^{{\text{X}}^{\prime }}{]}_{n_x\times n_x}$和${\text{Q}}^{{\text{Y}}^{\prime }}=[{\text{Q}}_{ij}^{{\text{Y}}^{\prime }}{]}_{n_y\times n_y}$则是在$\text{X}$和$\text{Y}$降维后的公共空间中定义的细胞到细胞的转移概率矩阵。KL散度定义为：$$KL(\text{P}||\text{Q})=\sum _i\sum _j{\text{P}}_{ij}log\frac{{\text{P}}_{ij}}{{\text{Q}}_{ij}}$$损失函数$L({\text{X}}^{\prime },{\text{Y}}^{\prime })$有三项：两个KL散度，和一个用于测量公共embedding空间中匹配数据之间距离的耦合项。$\beta$是一个数值平衡参数。我们用梯度下降优化该问题，并且在原始空间特征数量大于100时，将$p=32$作为默认值。

2.2.Data

2.2.1.模拟数据

我们模拟了两组单细胞多组学数据集，如下所示：

• Simulation 1 包含两个数据集，它们共享一个类似的低维拓扑结构（图2a）
• Simulation 2 包含两个数据集，Dataset 1具有embedding在2D空间的分叉树（图2b左上图），为了测试提出的全局缩放因子$\alpha$的能力，Dataset 2被非线性扭曲并embedding到3D空间，此外还添加了一个分支，该分支具有与其他分支不同的统计信息，以评估UnionCom识别特定于数据集细胞类型的能力。

因此，Simulation 2 的 Dataset 2具有embedding在3D空间中的三叉树，其中一个分支作为数据集特定的细胞类型（图2b，左下图，绿色分支为数据集特定的细胞类型）。

同一Simulation中两个数据集之间的细胞对应信息从Simulation中已知。我们将两类模拟数据集投影到特征空间中（对于每个Simulation，Dataset 1和Dataset 2的维数分别为1000和500）。

2.2.2.真实单细胞多组学数据

我们使用genetype、expression和methylation（sc-GEM）数据（Cheow等人在2016年提出）和single-cell nucleosome、methylation和transcription（sc-NMT）数据（Clark等人在2018年提出）的两组真实单细胞多组学数据。由sc-GEM生成的第一个实数集包含两个单细胞组学数据集，分别是人类细胞重新编程为诱导多能干细胞（iPS）的样本的gene expression和DNA methylation。通过sc-NMT生成的第二个真实数据集包含三个单细胞组学数据集，分别是在三个时间阶段（胚胎第5.5天（E5.5）、E6.5和E7.5）收集的小鼠原肠胚化样本的gene expression、DNA methylation（DNA甲基化）和chromatin accessibility（染色质可及性）。

sc-GEM同时测定同一细胞的基因表达和DNA甲基化，sc-NMT同时测量同一细胞的基因表达、DNA甲基化和染色质可及性。因此，对于每个真实数据集，单细胞多组学数据集中的细胞对应信息都是可用的。

2.3.评价方法

我们使用两个指标来评估单细胞多组学整合方法：Neighborhood Overlap和Label Transfer Accuracy。两个指标都基于整合数据集的公共嵌入空间计算。

当多组学数据集之间的细胞-细胞对应信息可用时，Stanley等人在2018提出的邻域重叠（Neighborhood Overlap）用于测量两个数据集之间细胞一对一匹配关系的能力，邻域重叠百分比定义为一个数据集可以从另一个数据集中找到对应细胞（即Neighborhood）的百分比，我们使用两个数据集的平均Neighborhood overlap作为最终的Neighborhood Overlap。邻域重叠的平均百分比在0%到100%之间，较高的百分比表示两个数据集之间的细胞到细胞匹配关系得到更好的复原。

当细胞标签信息（例如细胞类型）可用时，Label Transfer Accuracy标签转移准确性（已在Transfer Learning中广泛使用，并被Johansen和Quon在2019提出的scAlign采用）用于测量细胞标签在公共嵌入空间中从一个数据集转移到另一个数据集中的能力。假设Dataset 2比Dataset 1有更多细胞，我们使用Dataset 2作为training set，Dataset 1作为testing set。我们构建$k^{acc}$-$nn$分类器并在Dataset 2上训练，Label Transfer Accuracy是在testing set上的细胞标签预测准确性。标签转移准确性的值，范围从0%到100%，更高的百分比表示两个对齐数据集的转移标签的性能更好。我们设置$k^{acc}=5$，但实际上，不同大小的$k^{acc}$对标签转移准确性的影响是稳定的。

3.实验结果

• 图2：模拟数据集的整合；
• a和b：Simulation 1和Simulation 2，数据集1（左上图）和数据集2（左下图）的可视化，在对齐前后分别使用t-SNE；通过MMD-MA（中上面板：点根据其对应的数据集着色；中下面板：点按照其对应的分支着色）和UnionCom可视化两个对齐数据集的公共嵌入空间（右上面板：点根据其对应的数据集着色；右下面板：点按照其相应的分支着色）。模拟2的数据集2的绿色分支（图b的左下面板）是数据集特定的细胞类型，具有特定于数据集2唯一的拓扑结构。
• c：不同邻域大小下邻域重叠的平均百分比（左面板：Simulation 1；右面板：Simulation 2）。

• 图3：在2项模拟数据实验中，通过8种方法确定了标签转移精度。UnionCom1：UnionCom使用欧几里德距离代替测地距离；UnionCom2：UnionCom使用固定的$\alpha =\sqrt{d_y/d_x}$。

• 图4：sc-GEM的多组学数据整合（gene expression 和 DNA methylation）；
• a：在比对前先使用t-SNE对gene expression 和 DNA methylation数据集进行可视化；
• 分别由Seurat（b）、MMD-MA（c）、Harmony（d）和UnionCom（e）对两个数据集整合的公共嵌入空间进行可视化，左面板：点（代表细胞）根据其相应的数据集着色；右面板：点（代表细胞）根据相应的细胞类型着色。
• f：不同邻域大小下邻域重叠的平均百分比；
• g：不同$k^{acc}$的标签转移准确性；

• 图5：sc-NMT的多组学数据整合（gene expression，DNA methylation 和 chromatin accessibility）；
• a：在比对前分别使用UMAP对基因表达、DNA甲基化和染色质可及性数据集进行可视化；
• b：分别通过Seurat（左面板）、MMD-MA（中面板）和Harmony（右面板）可视化DNA甲基化和染色质可及性的两个数据集整合的公共嵌入空间；
• c：UnionCom对DNA甲基化和染色质可及性的两个数据集整合（左图）和DNA甲基化、染色质可及性和基因表达的三个数据集中的公共嵌入空间进行可视化（右图）；
• d：在DNA甲基化和染色质可及性两个数据集的整合中，邻域重叠的平均百分比（上图）和标签转移准确性（下图）；

• 图6：UnionCom在DNA甲基化和染色质可及性两个数据集（sc-NMT）上的数据子采样特征和参数选择方面的稳健性。
• a：在整合之前，分别从每个DNA甲基化和染色质可及性数据集中随机抽样特征子集，然后再整合，再记录UnionCom、Seurat、MMD-MA和Harmony的标签转移准确性。
• b：当在UnionCom的步骤A1中选择k-nn graph的不同$k$时，UnionCom的标签转移精度。
• c和d和e：$\rho$是在做矩阵优化时，构造拉格朗日乘子引入的参数（这部分属于附加材料内容），$\beta$为损失函数$L$中的数值平衡参数；
• f：在UnionCom的步骤A3中，将两个数据集嵌入不同维度$p$的公共空间时，UnionCom的标签传输精度。