SCOT: Single-Cell Multi-Omics Alignment with Optimal Transport

2.1. Optimal transport

The Kantorovich optimal transport problem seeks to find a minimal cost mapping between two probability distributions or discrete measures (Peyré et al., 2019). Referring back to the problem of moving a sand pile to fill in a hole, the Kantorovich optimal transport allows us to split the mass of a grain of sand instead of moving the whole grain; therefore, the mappings need not be 1–1. Consider discrete measures $\mu$ and $\nu$ as such

where ${\sum }_{i=1}^{n_x}{p}_i=1={\sum }_{j=1}^{n_y}{q}_j,p_i\ge 0,q_j\ge 0$ and ${\delta }_{x_i}$ is the Dirac measure. This optimal transport problem finds a minimal coupling $\pi$ that attains

where $c\left(i,j\right)$ is a cost function defined over the samples from the two data sets and $\Pi \left(\mu ,\nu \right)$ is the set of couplings of $\mu$ and $\nu$ given by

Intuitively, the cost function says how many resources it will take to move point x_i in the first data set to point y_j in the second data set, and the coupling $\pi$ relates the two discrete measures $\mu$ and $\nu$ by correspondence probabilities. Each row ${\pi }_i$ tells us how to split the mass of data point x_i onto the points y_j for $j=1,\dots ,n_y$, and the condition $\pi 1_{n_y}=p$ requires that the sum of each row ${\pi }_i$ is equal to p_i, the probability of sample x_i. The discrete optimal transport problem finds a coupling matrix, $\Gamma$, that minimizes the cost of moving samples through the linear program:

Although this problem can be solved with minimum cost flow solvers, it is usually regularized with entropy for more efficient optimization and empirically better results (Cuturi, 2013). Entropy diffuses the optimal coupling, meaning that more masses will be split. Thus, the numerical optimal transport problem is

where $\in >0$ and $H\left(\Gamma \right)$ is the Shannon entropy (${\sum }_{i=1}^{n_x}{\sum }_{j=1}^{n_y}{\Gamma }_{ij}log{\Gamma }_{ij}$).

Equation (5) is a strictly convex optimization problem, and for some unknown vectors $u\in {ℛ}^{n_x}$ and $v\in {ℛ}^{n_y}$, the solution has the form ${\Gamma }^{\ast }=diag\left(u\right)Kdiag\left(v\right),$ with $K=exp\left(-\frac{C}{\delta }\right)$, element-wise. This solution can be obtained efficiently via Sinkhorn's algorithm, which iteratively computes

where $\%$ denotes element-wise division. This derivation immediately follows from solving the corresponding dual problem for Equation (5) (Peyré et al., 2019).

2.2. The Gromov–Wasserstein optimal transport

While the classic optimal transport formulation requires us to define a cost function across domains [Eq. (2)], this is difficult to do when working with data from different metric spaces. This is because we cannot directly compare data points with different modalities, such as in the case of multi-omic alignment. The Gromov–Wasserstein distance extends optimal transport by comparing distances between data points rather than directly comparing the data points themselves (Alvarez-Melis and Jaakkola, 2018) and allows us to work with data from different modalities. Consider the same discrete measures $\mu$ and $\nu$ as above, the cost function in the formulation of the optimal transport problem will now be defined over sample-wise pairwise distances $d_x\left(i,k\right)$ and $d_y\left(j,l\right)$ in the X and Y data sets, respectively:

where L indicates the cost function. The main change from basic optimal transport [Eq. (2)] to the Gromov–Wasserstein optimal transport [Eq. (7)] is that we consider the effect of transporting pairs of samples rather than single samples. Intuitively, captures how transporting x_i to y_j and x_k to y_l would distort the original distances between i and k and between x_j and x_l. This change ensures that the optimal transport plan $\pi$ will preserve some local geometry.

For solving the Gromov–Wasserstein optimal transport formulation, we compute pairwise distance matrices D^x and D^y for the two domains separately, as well as the fourth-order tensor $L\in {ℛ}^{n_x\times n_x\times n_y\times n_y}$, where $L_{ijkl}=L\left(D_{ik}^x,D_{jl}^y\right)$. Then, the discrete Gromov–Wasserstein problem can also be expressed as the inner product

Equation (8) is now both nonlinear and nonconvex and involves operations on a fourth-order tensor, including the $O\left(n_x^2{n}_y^2\right)$ operation tensor product $L\left(D^x,D^y\right)\otimes \Gamma$ for a naive implementation. Peyré et al. (2016) showed that for some choices of loss function this product can be computed in $O\left(n_x^2{n}_y+n_x{n}_y^2\right)$ cost. In particular, for the case $L=L_2$, the inner product can be computed by

As in the classic optimal transport case, the coupling matrix can be efficiently computed for an entropically regularized optimization problem:

Larger values of $\in$ lead to not only an easier optimization problem but also a denser coupling matrix, meaning that solutions will indicate significant correspondences between more data points. Smaller values of $\in$ lead to sparser solutions, meaning that the coupling matrix is more likely to find the correct one-to-one correspondences for data sets where there are one-to-one correspondences. However, it also yields a harder (more nonconvex) optimization problem (Alvarez-Melis and Jaakkola, 2018).

Peyré et al. (2016) proposed using a projected gradient descent approach for optimization, where both the projection and the gradient are taken with respect to the Kullback–Leibler divergence. These projections are computed via the Sinkhorn iterations. Algorithm 1 in the Supplementary Materials presents the algorithm for $L=L_2.$

Algorithm 1: Unsupervised hyperparameter search procedure
Input: Data sets $X,Y$.
$n\leftarrow min\left(n_x,n_y\right)$, $k_1\leftarrow min\left(0.2n,50\right)$
// Fix k1 and vary $\in$
// Fix ${\in }_1$ and vary k
if $n>250$ then
$k_2\leftarrow arg{min}_{k\in \left[20,100\right]}SCOT\left(X,Y,k,{\in }_1\right)$
end
else
$k_2\leftarrow arg{min}_{k\in \left[0.05n,0.2n\right]}SCOT\left(X,Y,k,{\in }_1\right)$
end
// Do a more refined search around k2 and ${\in }_1$
Return: $k_{best},{\in }_{best}$

Open in a new tab

2.3. Single-cell alignment with optimal transport

Our method, SCOT, works as follows. First, we compute the pairwise distances on our data in a way similar to Nitzan et al. (2019). To do this, we use the correlations between data points within each data set to construct k-NN connectivity graphs. We find that connectivity graphs, which connect nodes with binary edges, empirically work better than weighted edges. This could be because connectivity graphs potentially denoise the data. Next, we compute the shortest path distance on the graph between each pair of nodes via Dijkstra's algorithm.

We set the distance of any unconnected nodes to be the maximum finite distance in the graph and normalize the matrix by dividing the elements by this maximum distance. If k is the number of samples, then the k-NN graph is the complete graph, so the corresponding distance matrix is a matrix of all ones. In this case, the distance matrix does not provide information about the local geometry, so we recommend keeping k small relative to the number of samples to avoid this scenario. We find that our approach is robust to the choice of k (Fig. 5).

Since we do not know the true distribution of the original data sets, we follow the study by Alvarez-Melis and Jaakkola (2018) and empirically set $\mu$ and $\nu$ to be the uniform distributions on the data points. Then, we solve for the optimal coupling $\Gamma$, which minimizes Equation (10). To implement this method, we use the Python Optimal Transport toolbox (https://pot.readthedocs.io/en/stable) (Flamary and Courty, 2017).

One of the advantages of using optimal transport is the probabilistic interpretation of the resulting coupling matrix $\Gamma$, where the entries of the normalized row $\frac{1}{p_i}{\Gamma }_i$ are the probabilities that the fixed data point x_i corresponds to each y_j. However, to use the evaluation metrics previously used in the field and to visualize alignment, we need to project the two data sets into the same space. The Procrustes approach proposed in the study by Alvarez-Melis and Jaakkola (2018) does not generalize to data sets with different feature and sample dimensions, so we use a barycentric projection:

2.4. Alternative unsupervised alignment procedure

In the description of SCOT, the number k for nearest neighbors and the entropy weight $\in$ are hyperparameters. One way to set these hyperparameters for optimal alignment is to use some orthogonal correspondence information to select the best alignment either directly (Liu et al., 2019; Cao et al., 2020) or by performing cross-validation (Singh et al., 2020). This selection strategy is problematic for truly unsupervised setting, where no correspondence information is available a priori upon sequencing separate cell cultures.

As a solution, we provide an alternative procedure to learn reasonable alignments based on tracking the Gromov–Wasserstein distance [Eq. (8)]. This procedure is based on our observation that the Gromov–Wasserstein distance serves as a proxy for measuring alignment quality (Fig. 4A). In this procedure, we alternate between optimizing $\in$ and k to minimize the Gromov–Wasserstein distance between the domains (detailed in Algorithm 1). Although the lowest Gromov–Wasserstein distance is not always the best alignment, it consistently appears to be one of the better alignments.

3.1. Simulated data sets

We follow the study by Liu et al. (2019) and benchmark SCOT on three different simulations (https://noble.gs.washington.edu/proj/mmd-ma). All three simulations contain two domains with 300 samples that have been nonlinearly projected to 1000- and 2000-dimensional feature spaces, respectively. The three simulations are a bifurcation, a Swiss roll, and a circular frustum (Fig. 2) with points belonging to three different groups. In addition to these three previously existing simulations, we use Splatter (Zappia et al., 2017) to create simulated single-cell RNA sequencing count data, which we call synthetic RNA-seq. We generate 5000 cells with 1000 genes from 3 cell groups and reduce the count matrix to the 5 genes with the highest variances. This count matrix is mapped into two new domains with dimensions $p_1=50$ and $p_2=500$ by multiplying it with two randomly generated matrices, resulting in data with dimensions $5000\times 50$ and $5000\times 500$.

All four data sets were simulated with 1–1 sample-wise correspondences, which are solely used for evaluating model performance. Each domain is projected to a different dimension, so there is no feature-wise correspondence either. In all simulations, we Z-score normalize the features before running the alignment algorithms as in the study by Liu et al. (2019).

3.2. Single-cell multi-omics data sets

We use two sets of single-cell multi-omics data to demonstrate the applicability of our model to real data sets. Both data sets are generated by co-assays; thus, we have known cell-level correspondence information for benchmarking. The first data set is generated using the scGEM assay (Cheow et al., 2016), which simultaneously profiles gene expression and DNA methylation. The data set (Sequence Read Archive accession SRP077853) is derived from human somatic cell samples undergoing conversion to induced pluripotent stem cells and shows a continuous trajectory. This data set was also used by Cao et al. (2020) to demonstrate the performance of their UnionCom algorithm. We preprocessed the data as described in the original publications (Cheow et al., 2016; Cao et al., 2020) and ended up with dimensions 177 × 34 for the gene expression data and 177 × 27 for the chromatin accessibility data.

The second data set is generated by the SNAREseq assay (Chen et al., 2019), which links chromatin accessibility with gene expression. The data (Gene Expression Omnibus accession GSE126074) is derived from a mixture of human cell lines: BJ, H1, K562, and GM12878 and show distinct cell type clusters. We preprocess the data sets following the study by Chen et al. (2019). The resulting data matrices for the SNARE-seq data set were of size 1047 × 19 and 1047 × 10 for ATAC-seq and RNA-seq, respectively. We unit normalize all real data sets as done in the study by Singh et al. (2020). Both data sets that we work with have been previously published and made publicly available by the original authors. Therefore, they do not require IRB approval.

3.3. Evaluation metrics

We compare SCOT with the two state-of-the-art unsupervised single-cell alignment methods MMD-MA (Liu et al., 2019) and UnionCom (Cao et al., 2020). None of these methods use any correspondence information for aligning the data sets. However, all data sets have 1–1 sample-level correspondence information, which we use to quantify the alignment performance through the “fraction of samples closer than the true match” (FOSCTTM) metric introduced by Liu et al. (2019). For each domain, we compute the Euclidean distances between a fixed sample point and all the data points in the other domain. Next, we use these distances to compute the fraction of samples that are closer to the fixed sample than its true match. Finally, we average these values for all the samples in both domains. For perfect alignment, all samples would be closest to their true match, yielding an average FOSCTTM of zero. Therefore, a lower average FOSCTTM corresponds to better alignment performance.

Since all the data sets have group-specific (simulations) or cell-type-specific (real experiments) labels, we also adopt the metric used by Cao et al. (2020) called “label transfer accuracy” (LTA) to assess the quality of the cell label assignment and to allow for a more direct comparison with their results. This metric measures the ability to correctly transfer sample labels from one domain to another based on their neighborhood in the aligned domain. As described in the study by Cao et al. (2020), we train a k-NN classifier on one of the domains and predict the sample labels in the other domain. The LTA is the proportion of correctly predicted labels, so it ranges from 0 to 1, and higher values indicate good performance. We apply this metric to alignments selected by the FOSCTTM measure.

3.4. Hyperparameter tuning

We run each method over a grid of hyperparameters and select the setting that yields the lowest average FOSCTTM. For SCOT, the grid covers the regularization weight $\in \in \left\{0.0001,0.0005,0.001,0.005,\dots ,0.1\right\}$ and number of neighbors $k\in \left\{10,15,20,25,30,35,\dots 100,\frac{1}{6}{n}_x\right\}$. We observe empirically that going above $\frac{1}{6}n$ for k does not yield any improvement in alignment.

We pick the hyperparameters for MMD-MA and UnionCom based on the default values and recommended ranges. MMD-MA has three hyperparameters: weights ${\lambda }_1,{\lambda }_2\in \left\{10^{-3},10^{-4},10^{-5},10^{-6},10^{-7}\right\}$ for the terms in the optimization problem and the dimensionality $p\in \left\{4,5,6,16,32,64\right\}$ of the embedding space. UnionCom requires the user to specify four hyperparameters: the number $kmax\in \left\{40,100\right\}$ of maximum number of neighbors in the graph, the dimensionality $p\in \left\{4,5,6,16,32,64\right\}$ of the embedding space, the trade-off parameter $\beta \in \left\{0.1,1,10,15,20\right\}$ for the embedding, and a regularization coefficient $\rho \in \left\{0,5,10,15,20\right\}$. We select the embedding dimension $p\in \left\{16,32,64\right\}$ around the default value of 32 set by UnionCom but also add $p\in \left\{4,5,6\right\}$ to match the recommended values for MMD-MA. We keep the hyperparameter search space size approximately consistent across the three methods. For each data set, we present alignment and runtime results for the best performing hyperparameters.

Furthermore, we consider the scenario where correspondence information is unavailable to pick the optimal hyperparameters. For SCOT, we apply the alternative unsupervised alignment algorithm (Algorithm 2 in the Supplementary Materials S1) to align all the data sets. Since MMD-MA and UnionCom do not provide a hyperparameter selection strategy, we rely on the default hyperparameters; we use UnionCom's provided default parameters of $kmax=40,p=32,\rho =10,$ and $\beta =1$, and the center values of MMD-MA's recommended range: $p=5,{\lambda }_1=10^{-5},$ and ${\lambda }_2=10^{-5}.$ We also present the alignment results for all three methods in this fully unsupervised setting.

4.1. SCOT successfully aligns the simulated data sets

In this experiment, we align the three simulation data sets from the study by Liu et al. (2019), as well as the synthetic single-cell RNA-seq count data generated with Splatter (Zappia et al., 2017). Before alignment, we first select the best performing hyperparameters for each method using the ground-truth correspondence information, as described in Section 3.4.

In Figure 2, we visualize the original domains, as well as the alignment performed by SCOT. We color the samples by their domain and cell-type identity. We observe that the global structure is matched, and cells cluster correctly based on cell-type identity. We then sort and plot the FOSCTTM score for each sample in Figure 2C. The mean FOSCTTM values are summarized in Table 1. We also report the LTA values in Table 2 when the first domain is used to train a classifier to predict the labels in the second domain. Overall, we observe that SCOT consistently achieves one of the lowest average FOSCTTM scores, thereby demonstrating its ability to recover the correct correspondences. SCOT also consistently yields high LTA scores indicating that samples are correctly mapped to their assigned groups.

4.2. SCOT gives state-of-the-art performance for single-cell multi-omics alignment

Next, we apply our method to real single-cell sequencing data and visualize the alignments in Figure 3. To have ground-truth information on cell–cell correspondences solely for benchmarking purposes, we use data sets generated by co-assaying technology. Overall, SCOT gives the lowest average FOSCTTM measure in comparison to MMD-MA and UnionCom (Table 1) and recovers accurate 1–1 correspondences in single-cell data sets. For the scGEM data, we report LTA using the DNA methylation domain for predicting the cell-type labels in the gene expression domain.

For the SNARE-seq data set, we use the gene expression domain for predicting cell labels in the chromatin accessibility domain (Table 2). SCOT yields the best LTA result on SNAREseq data set and performs comparably to the other methods for scGEM. All methods have higher LTA performance on SNAREseq data set compared with scGEM data set because SNAREseq data set contains a mixture of different cell types that cluster separately, whereas scGEM data set contains cells going through a continuous differentiation.

While MMD-MA and UnionCom project both data sets to a shared low-dimensional space, SCOT projects one data set onto the other. We find that the direction of projection makes no significant difference in performance (Supplementary Table S1).

4.3. SCOTs alternative unsupervised hyperparameter tuning procedure achieves quality alignments

We compare the alignment performances in fully unsupervised settings, when we have no validation data on correspondences to use for hyperparameter tuning, as described in Section 3.4. We present the alignment performances, measured by average FOSCTTM measures, in Table 3, and by LTA in Table 4, when using SCOTs alternative self-tuning procedure. In this procedure, hyperparameter choice is guided by the Gromov–Wasserstein distance, as we have observed a correlation between the Gromov–Wasserstein distances between the aligned data sets and alignment quality (Fig. 4A).

In this unsupervised setting, we use MMD-MA's and UnionCom's default parameters since they lack self-tuning capability. SCOT returns nearly the same alignments for simulated data and only marginally worse alignments for real data. In contrast, MMD-MA and UnionCom show inconsistent alignment performance and fail to align some of the simulated and all real data sets with the default parameter values. Therefore, the proposed procedure could guide a user to an alignment close to the optimal result when no orthogonal information is available.

4.4. SCOTs computation speed scales well with the sample size

We compare SCOTs running times with the baseline methods for the best performing hyperparameters on the synthetic RNA-seq data set by varying the number of cells to demonstrate how each algorithm scales to larger data sets. While SCOT is implemented for CPU, both MMD-MA and UnionCom algorithms provide GPU versions, which run faster. Therefore, we use them for benchmarking. We run CPU computations on an Intel Xeon e5-2670 with 16 GB memory and GPU computations on a single NVIDIA GTX 1080ti with VRAM of 11 GB. SCOTs running time scales similar to that of MMD-MA, even though SCOT runs on a CPU and MMD-MA runs on a GPU (Fig. 4B). Both methods scale better than the GPU-based UnionCom implementation.

4.5. Investigating algorithmic choices and hyperparameters of SCOT

To better understand our method, we investigated the effects of different algorithmic choices and hyperparameter combinations on the alignment performance of the real-world data sets. Figure 5 shows the range of average FOSCTTM values we receive for alignments with different combinations of k (number of neighbors in k-NN graphs) and $\in$ (entropic regularization coefficient) values for the two real-world sequencing data sets. Overall, we observe that the choice of $\in$ tends to make a larger impact on the alignment performance than k. Next, we consider the effect of different algorithmic choices on the alignment performance of SCOT.

We compare the final SCOT model with (1) no entropic regularization, (2) using Euclidean distances for intra-domain distance matrices, and (3) using correlation-based intra-domain distance matrices in lieu of graph distances. For each of these settings, we run alignments for the same combinations of hyperparameters as described in Section 3.4 and record the average FOSCTTM measure we receive for each alignment. In Figure 6, we compare these in violin plots for scGEM and SNARE-seq data sets. This experiment shows that both entropic regularization and modeling the single-cell data sets as graphs for intra-domain distance computations yield lower FOSCTTM measures, corresponding to higher quality alignments.