Unsupervised Integration of Single-Cell Multi-omics Datasets with Disproportionate Cell-Type Representation

Abstract

Integrated analysis of multi-omics data allows the study of how different molecular views in the genome interact to regulate cellular processes; however, with a few exceptions, applying multiple sequencing assays on the same single cell is not possible. While recent unsupervised algorithms align single-cell multi-omic datasets, these methods have been primarily benchmarked on co-assay experiments rather than the more common single-cell experiments taken from separately sampled cell populations. Therefore, most existing methods perform subpar alignments on such datasets. Here, we improve our previous work Single Cell alignment using Optimal Transport (SCOT) by using unbalanced optimal transport to handle disproportionate cell-type representation and differing sample sizes across single-cell measurements. We show that our proposed method, SCOTv2, consistently yields quality alignments on five real-world single-cell datasets with varying cell-type proportions and is computationally tractable. Additionally, we extend SCOTv2 to integrate multiple ($M\ge 2$) single-cell measurements and present a self-tuning heuristic process to select hyperparameters in the absence of any orthogonal correspondence information.

Available at: http://rsinghlab.github.io/SCOT.

Appendix

1.1 Embedding Method Details

The full details of t-SNE can be found in [13]. For each domain m, we compute $P^{m}$, an $n_m \times n_m$ cell-to-cell transition matrix; each entry $P^{m}_{j|i}$ is the conditional probability that a data point $x_i^m$ would pick $x_j^m$ as its neighbor when chosen according a Gaussian distribution centered at $x_i^m$:

$$\begin{aligned} P_{j|i}^m =\frac{\exp (- ||x_i^m - x_j^m||^2 / 2\sigma _i^2)}{\sum _{k \ne i} \exp (-||x_i^m - x_k^m||^2/2\sigma _i^2)}. \end{aligned}$$

(9)

The bandwidth $\sigma _i$ is chosen according to the density of the data points through a binary search for the value of $\sigma _i$ that achieves the user-supplied perplexity value. $P^m$ is computed by averaging $P^m_{i|j}$ and $P^m_{j|i}$ to give more weight to outlier points:

$$\begin{aligned} P^m_{ij} = \frac{P_{i|j}^m + P_{j|i}^m}{2 n_m} \end{aligned}$$

(10)

Then, to jointly embed all domains through the anchor domain $X^1$, the optimization problem is:

$$\begin{aligned} \min _{X^{1'}, \dots , X^{M'}} \sum _{m=1}^M\text {KL} (P^{m} || Q^{m'}) + \beta \sum _{m=2}^M ||X^{1'}-X^{m'} (\varGamma ^m)^T ||^2_F, \end{aligned}$$

(11)

where $X^{m'}$ is the lower dimensional embedding of $X^m$, $P^m$ is defined as in Eq. 9, and $\varGamma ^m$ is the coupling matrix from solving Eq. 6 for $m=1,2, \dots , M$, $X^{m'}$. The probability matrix $Q^m$ is computed through a Student-t distribution with one degree of freedom:

$$\begin{aligned} Q^{m'}_{ij} = \frac{ (1 + ||x_i^{m'} - x_j^{m'}||)^{-1} }{ \sum _{k \ne l} 1 + ( || x_k^{m'} - x_l^{m'} || )^{-1} }. \end{aligned}$$

(12)

The intuition behind the cost $\text {KL} (P^{m} || Q^{m'})$ is very similar to that of GW; if two points have a high transition probability in the original space, then they should also have a high transition probability in the latent space.

1.2 Hyperparameter Tuning Procedure Details

For each alignment method, we define a grid of hyperparameters and choose the best performing combination for each experiment. If methods share similar hyperparameters in their formulation, we keep the range defined for these consistent across all algorithms. We refer to the publication and the code repository for each method to choose a hyperparameter ranges whenever possible.

For Pamona, we search the number of neighbors in the cell neighborhood graphs, $k \in \{20, 30, \dots , 150\}$, the entropic regularization coefficient, $\epsilon \in \{5e{-}4, 3e{-}4, 1e{-}4, 7e{-}3, 5e{-}3, \dots , 1e{-}2 \}$, geometry preservation trade-off coefficient, $\lambda \in \{0.1, 0.5, 1, 5, 10\}$, and lastly, embedding dimensionality, $p \in \{3, 4, 5, 10, 30, 32\}$, the output dimension for embedding. For UnionCom, we search the trade-off parameter $\beta \in \{0.1, 1, 5, 10, 15, 20\}$, the regularization coefficient $\rho \in \{0, 0.1, 1, 5, 10, 15, 20\}$, the maximum neighborhood size permitted in the neighborhood graphs, $k_{max} \in \{40, 100, 150\}$, and embedding dimensionality $p \in \{3, 4, 5, 10, 30, 32\}$. For MMD-MA:, we tune the weights $\lambda _1$ and $\lambda _2$ $\in \{1e{-}2, 5e{-}3, 1e{-}3, 5e{-}4, \dots , 1e{-}9\}$, and the embedding dimensionality, $p \in \{3,4,5,10,30,32\}$. For bindSC, we choose the coefficient that assigns weight to the initial gene activity matrix $\alpha \in \{0, 0.1, 0.2, \dots 0.9\}$, the coefficient that assigns weight factor to multi-objective function $\lambda \in \{0.1, 0.2, \dots , 0.9\}$, and the number of canonical vectors for the embdedding space $K \in \{3, 4, 5, 10, 30, 32\}$. Lastly, for Seuratv4, we tune the number of neighbors to consider when finding anchors, $k \in \{5, 10, 15, 20\}$, co-embedding dimensionality, $p \in \{3, 4, 5, 10, 30, 32\}$ and the choice of the reference and anchor domains when finding anchors.