PAPER KEY: 8FKE6UJ8
TITLE: SCOT: Single-Cell Multi-Omics Alignment with Optimal Transport
AUTHORS: Demetci, Pinar; Santorella, Rebecca; Sandstede, Björn; Noble, William Stafford; Singh, Ritambhara

Research Articles
SCOT: Single-Cell Multi-Omics Alignment with Optimal
Transport
PINAR DEMETCI,1,2,*,i REBECCA SANTORELLA,3 BJO ̈ RN SANDSTEDE,3,* WILLIAM STAFFORD NOBLE,4,5 and RITAMBHARA SINGH1,2,ii
ABSTRACT
Recent advances in sequencing technologies have allowed us to capture various aspects of the genome at single-cell resolution. However, with the exception of a few of co-assaying technologies, it is not possible to simultaneously apply different sequencing assays on the same single cell. In this scenario, computational integration of multi-omic measurements is crucial to enable joint analyses. This integration task is particularly challenging due to the lack of sample-wise or feature-wise correspondences. We present single-cell alignment with optimal transport (SCOT), an unsupervised algorithm that uses the Gromov–Wasserstein optimal transport to align single-cell multi-omics data sets. SCOT performs on par with the current state-of-the-art unsupervised alignment methods, is faster, and requires tuning of fewer hyperparameters. More importantly, SCOT uses a self-tuning heuristic to guide hyperparameter selection based on the Gromov–Wasserstein distance. Thus, in the fully unsupervised setting, SCOT aligns single-cell data sets better than the existing methods without requiring any orthogonal correspondence information.
Keywords: data integration, manifold alignment, multi-omics, optimal transport, single-cell genomics.
1. INTRODUCTION
T
he growing variety of single-cell assays allows us to measure the heterogeneous landscape of cell state in a sample, revealing distinct subpopulations and their developmental and regulatory trajectories across time. Different technologies can interrogate different molecular aspects of the cell, such as gene expression, protein synthesis, chromatin accessibility, DNA methylation, histone modifications, and chromatin three-dimensional (3D) confirmation. Combining data generated by these single-cell assays can
1Center for Computational Molecular Biology, Brown University, Providence, Rhode Island, USA. 2Department of Computer Science, Brown University, Providence, Rhode Island, USA. 3Division of Applied Mathematics, Brown University, Providence, Rhode Island, USA. 4Department of Genome Sciences, University of Washington, Seattle, Washington, USA. 5Paul G. Allen School of Computer Science and Engineering, University of Washington, Seattle, Washington, USA. *These authors contributed equally to this work. iORCID ID (https://orcid.org/0000-0002-5644-0326). iiORCID ID (https://orcid.org/0000-0002-7523-160X).
JOURNAL OF COMPUTATIONAL BIOLOGY Volume 29, Number 1, 2022 # Mary Ann Liebert, Inc. Pp. 3–18 DOI: 10.1089/cmb.2021.0446
3
Downloaded by COLUMBIA UNIVERSITY LIBRARY from www.liebertpub.com at 02/16/25. For personal use only.


provide novel insights into the interactions between these molecular views and their joint regulatory mechanisms. Hence, learning this combined information is critical to our understanding of complex biological processes and heterogeneous diseases. Despite its importance, combining single-cell multi-omics data is a challenging task. Aside from a few recent co-assay procedures that simultaneously isolate separate molecular material for each measurement, applying multiple assays on the same single cell is impossible. Sometimes, sequencing assays need access to the same molecular material, such as with chromatin accessibility and 3D chromatin conformation capture assays. In such cases, the measurements are taken by dividing a cell population into subpopulations and assaying them separately, losing the potential for 1–1 correspondence of cells that is required for easy data integration. Moreover, in cases where we can take measurements in the same cell and preserve the 1–1 correspondences, the choice of the experimental method for processing cells and isolating molecular materials of interest can introduce additional challenges and noise in the co-assayed data (Hu et al., 2018). For example, for simultaneous isolation of DNA and RNA, there are two general approaches: physical separation of DNA and RNA followed by separate amplification, or simultaneous preamplification followed by physical separation of the two materials. For the first approach, separation techniques such as centrifugation and micropipetting are not high-throughput; however, high-throughput approaches (Macaulay et al., 2015; Angermueller et al., 2016) have been found to introduce variability in coverage and sequencing depth of various genomic regions in the isolated DNA (Hu et al., 2018). In recent years, computational methods have been developed to solve the single-cell data integration problem. Many of these methods combine different experiments from a single modality such as RNA sequencing for correcting batch effects (Welch et al., 2017, 2019; Amodio and Krishnaswamy, 2018; Barkas et al., 2019; Stuart et al., 2019). However, integrating data from multiple modalities such as gene expression and DNA methylation presents unique challenges. For example, when we measure different properties of a cell, we cannot a priori identify correspondences between features in the two domains. Accordingly, integrating two or more single-cell data modalities requires methods that rely on neither common cells nor features across the data types. This aspect prevents the application of some existing single-cell alignment methods to unsupervised settings because they require some correspondence information to perform alignment (Welch et al., 2017, 2019; Amodio and Krishnaswamy, 2018; Barkas et al., 2019; Stuart et al., 2019). Earlier versions of the popular batch integration method Seurat required correspondence information in the form of cells from a similar biological state that are shared across the two data sets (known as ‘‘anchor points’’). While a more recent version automatically selects these anchor points, it still requires features from one domain to be mapped to the other domain to perform the single-cell alignment (Stuart et al., 2019). This mapping might be possible for experiments such as gene expression and chromatin accessibility, where one can map the chromatin region read counts to the corresponding gene regions. However, it can be difficult to perform for other sequencing assay combinations. Furthermore, Cao et al. (2020) have shown that such methods do not yield quality alignments in unsupervised settings. Multiple approaches have tried to align data sets in an entirely unsupervised manner. One of the earliest attempts, the joint Laplacian manifold alignment algorithm, constructs eigenvector projections based on k-nearest neighbor (k-NN) graph Laplacians of the data (Wang and Mahadevan, 2009). The generalized unsupervised manifold alignment (GUMA) (Cui et al., 2014) algorithm seeks a 1–1 correspondence between two data sets based on optimization of a local geometry matching term. Liu et al. (2019) showed that these methods do not perform well on the single-cell alignment task and proposed a manifold alignment (MA) algorithm based on the maximum mean discrepancy (MMD) measure, called MMD-MA. Another method, UnionCom (Cao et al., 2020), extends GUMA to perform unsupervised topological alignment and makes it more suitable for single-cell multi-omics integration. While MMD-MA aims to match the global distributions of the data sets in a shared latent space, UnionCom emphasizes learning both local and global alignments between the two distributions. Neither method requires any correspondence information, either among samples or features, to perform an alignment. The respective articles demonstrate state-of-the-art performance on simulated and real data sets. Although these results are encouraging, MMD-MA and UnionCom require that the user specify three and four hyperparameters, respectively. Hyperparameter selection can significantly affect the quality of alignments. Therefore, in an unsupervised real-world setting with no validation data on correspondences, hyperparameter tuning can be difficult to perform and can lead to subpar alignments.
4 DEMETCI ET AL.
Downloaded by COLUMBIA UNIVERSITY LIBRARY from www.liebertpub.com at 02/16/25. For personal use only.


In this article, we propose an unsupervised alignment method based on optimal transport theory. Optimal transport finds the most cost-effective way to move data points from one domain to another. One way to think about it is as the problem of moving a pile of sand to fill in a hole through the least amount of work. Traditionally, optimal transport problems have been difficult to compute, especially for large-scale data sets. However, subsequent relaxations (Kantorovich, 1942; Peyre ́ et al., 2019) modify the original optimal transport problem, making it more applicable and easier to compute. Recently, several regularization procedures (Peyre ́ et al., 2016) have further improved the computational scalability of optimal transport. In biology, an emerging number of applications are using optimal transport to learn a mapping between data distributions (Alvarez-Melis and Jaakkola, 2018; Yang et al., 2018; Schiebinger et al., 2019; Yang and Uhler, 2019; Cang and Nie, 2020). Schiebinger et al. (2019) used it to study temporal changes in gene expression by using regularized unbalanced optimal transport to compute expression differences between time points. SpaOTsc (Cang and Nie, 2020) maps cells with high ligand expression onto cells with high receptor expression to recover cell signaling relationships in spatially resolved single-cell RNA-seq data sets. ImageAEOT (Yang et al., 2018) maps single-cell images to a common latent space through an autoencoder and then uses optimal transport to track cell trajectories. In related work, the same authors used autoencoders and optimal transport to learn transport maps among multiple domains (Yang and Uhler, 2019). However, the application of their method to single-cell data sets requires some form of supervision, such as class labels, to be used during transport. The classic optimal transport problem requires data sets from the same metric space. Me ́moli (2011) generalized optimal transport to the Gromov–Wasserstein distance, which compares metric spaces directly instead of comparing samples across spaces, making optimal transport suitable for multimodal alignment. In natural language processing, Alvarez-Melis and Jaakkola (2018) used this approach to measure similarities between pairs of words across languages to compute the similarity between languages. As far as we are aware, the only biological application of the Gromov–Wasserstein optimal transport comes from the study by Nitzan et al. (2019), which uses it to reconstruct the spatial organization of cells from transcriptional profiles. We present single-cell alignment with optimal transport (SCOT), an unsupervised algorithm that uses the Gromov–Wasserstein-based optimal transport to align single-cell multi-omics data sets (presented schematically in Fig. 1). Like UnionCom, SCOT aims to preserve local geometry when aligning single-cell data. SCOT achieves this by constructing a k-NN graph for each data set (or domain) and then computing graph distance matrices for each k-NN graph to capture the intra-domain distances. SCOT then finds a probabilistic coupling matrix that minimizes the discrepancy between the intra-domain distance matrices. Finally, it uses the coupling matrix to project one single-cell data set onto another through barycentric projection, thus aligning them. Unlike MMD-MA and UnionCom, SCOT requires tuning only two hyperparameters and is robust to the choice of one. We compare the alignment performance of SCOT with MMD-MA and UnionCom on four simulated and two real-world data sets. SCOT aligns data sets as well as the state-of-the-art methods and scales well with increasing numbers of samples. Moreover, we demonstrate that the Gromov–Wasserstein distance can guide SCOTs hyperparameter tuning in a fully unsupervised setting when no orthogonal alignment information is available. Thus, unlike other methods, SCOT provides a heuristic for hyperparameter selection without validation data. The source code for SCOT is publicly available at http:// rsinghlab.github.io/SCOT.
FIG. 1. Schematic of SCOT alignment of single-cell multi-omics data. A population of cells is aliquoted for different single-cell sequencing assays. SCOT constructs k-NN graphs based on sample-wise correlations and finds a probabilistic coupling between the samples of each domain that minimizes the distance between the two intra-domain graph distance matrices. Barycentric projection projects one domain onto another based on this coupling matrix. SCOT, single-cell alignment with optimal transport.
SINGLE-CELL ALIGNMENT WITH OPTIMAL TRANSPORT 5
Downloaded by COLUMBIA UNIVERSITY LIBRARY from www.liebertpub.com at 02/16/25. For personal use only.


2. METHODS
SCOT relies on the Gromov–Wasserstein optimal transport to move data points from one domain to another while preserving the original local geometry. The goal of the transport problem at the core of SCOT is to find an ideal ‘‘coupling’’ (also called ‘‘correspondence’’) matrix that describes the probability of alignment between each point across domains. In this section, we first introduce optimal transport theory, followed by its extension to the Gromov–Wasserstein distance. Then, we present the details of our algorithm. We have two data sets representing two domains, X = (x1‚ x2‚ . . . ‚ xnx ) from X and Y = (y1‚ y2‚ . . . ‚ yny ) from Y. The data sets have nx and ny points, respectively. We do not require any correspondence information or assume that there is any ground truth for 1–1 correspondence between samples or features, but we do assume that there is some underlying shared biology (e.g., cells across the data sets sharing a lineage or belonging to shared cell types), so that the data sets can be meaningfully aligned.
2.1. Optimal transport
The Kantorovich optimal transport problem seeks to find a minimal cost mapping between two probability distributions or discrete measures (Peyre ́ 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 l and  as such
l=
n Xx
i=1
pidxi and  =
n Xy
j=1
qjdyj ‚ (1)
where Pnx
i = 1 pi = 1 = Pny
j = 1 qj‚ pi  0‚ qj  0 and dxi is the Dirac measure. This optimal transport problem finds a minimal coupling p that attains
min
p2P(‚ l)
n Xx
i=1
n Xy
j=1
c(i‚ j)p(i‚ j) (2)
subject to : p(i‚ j)  0‚
n Xx
i=1
p(i‚ j) = qj‚
n Xy
j=1
p(i‚ j) = pi
where c(i‚ j) is a cost function defined over the samples from the two data sets and P(l‚ ) is the set of couplings of l and  given by
P(l‚ ) = fp 2 Rnx · ny
+ : p1ny = l‚ pT 1nx = g: (3)
Intuitively, the cost function says how many resources it will take to move point xi in the first data set to point yj in the second data set, and the coupling p relates the two discrete measures l and  by correspondence probabilities. Each row pi tells us how to split the mass of data point xi onto the points yj for j = 1‚ . . . ‚ ny, and the condition p1ny = p requires that the sum of each row pi is equal to pi, the probability of sample xi. The discrete optimal transport problem finds a coupling matrix, G, that minimizes the cost of moving samples through the linear program:
min
G2P(l‚ ) ÆG‚ Cæ: (4)
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
min
G2P(l‚ ) ÆG‚ Cæ -  H(G)‚ (5)
where  > 0 and H(G) is the Shannon entropy (Pnx
i=1
Pny
j = 1 Gij log Gij).
6 DEMETCI ET AL.
Downloaded by COLUMBIA UNIVERSITY LIBRARY from www.liebertpub.com at 02/16/25. For personal use only.


Equation (5) is a strictly convex optimization problem, and for some unknown vectors u 2 Rnx and v 2 Rny , the solution has the form G = diag(u)K diag(v)‚ with K = exp - C
e
 , element-wise. This solution can be obtained efficiently via Sinkhorn’s algorithm, which iteratively computes
u)l%Kv and v)%KT u‚ (6)
where % denotes element-wise division. This derivation immediately follows from solving the corresponding dual problem for Equation (5) (Peyre ́ 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 l and  as above, the cost function in the formulation of the optimal transport problem will now be defined over sample-wise pairwise distances dx(i‚ k) and dy(j‚ l) in the X and Y data sets, respectively:
GW(l‚ ) : = min
p2P(l‚ )
n Xx
i‚ k
n Xy
j‚ l
L(dx(i‚ k)‚ dy(j‚ l)) p (i‚ j) p(k‚ l): (7)
where L indicates the cost function. The main change from basic optimal transport [Eq. (2)] to the GromovWasserstein optimal transport [Eq. (7)] is that we consider the effect of transporting pairs of samples rather than single samples. Intuitively, L(dx(i‚ k)‚ dy(j‚ l)) captures how transporting xi to yj and xk to yl would distort the original distances between i and k and between xj and xl. This change ensures that the optimal transport plan p will preserve some local geometry. For solving the Gromov–Wasserstein optimal transport formulation, we compute pairwise distance matrices Dx and Dy for the two domains separately, as well as the fourth-order tensor L 2 Rnx · nx · ny · ny , where Lijkl = L(Dx
ik‚ Dy
jl). Then, the discrete Gromov–Wasserstein problem can also be expressed as the inner product
GW(l‚ ) = min
G2 P(l‚ ) ÆL(Dx‚ Dy) G‚ Gæ (8)
Equation (8) is now both nonlinear and nonconvex and involves operations on a fourth-order tensor, including the O(n2
x n2
y) operation tensor product L(Dx‚ Dy) G for a naive implementation. Peyre ́ et al.
(2016) showed that for some choices of loss function this product can be computed in O(n2
x ny + nx n2
y) cost. In particular, for the case L = L2, the inner product can be computed by
L(Dx‚ Dy) G = (Dx)2l1T
ny + 1nx T ((Dy)2)T - DxG(Dy)T : (9)
As in the classic optimal transport case, the coupling matrix can be efficiently computed for an entropically regularized optimization problem:
GW(l‚ ) = min
G2P(l‚ ) ÆL(Dx‚ Dy) G‚ Gæ - H(G): (10)
Larger values of  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  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). Peyre ́ 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 = L2:
SINGLE-CELL ALIGNMENT WITH OPTIMAL TRANSPORT 7
Downloaded by COLUMBIA UNIVERSITY LIBRARY from www.liebertpub.com at 02