# Joint Analysis of Single-Cell Data across Cohorts with Missing Modalities

Marianne Arriola  Department of Computer Science, Cornell Tech, Cornell
UniversityNew YorkNYUSA10065<EMAIL_ADDRESS>,  Weishen Pan  Department of
Population Health Sciences, Weill Cornell Medicine, Cornell UniversityNew
YorkNYUSA10065<EMAIL_ADDRESS>,  Manqi Zhou  Department of
Computational Biology, Cornell UniversityIthacaNYUSA14853<EMAIL_ADDRESS>,
Qiannan Zhang  Department of Population Health Sciences, Weill Cornell
Medicine, Cornell UniversityNew YorkNYUSA10065<EMAIL_ADDRESS>,
Chang Su  Department of Population Health Sciences, Weill Cornell Medicine,
Cornell UniversityNew YorkNYUSA10065<EMAIL_ADDRESS>and  Fei Wang
Department of Population Health Sciences, Weill Cornell Medicine, Cornell
UniversityNew YorkNYUSA10065<EMAIL_ADDRESS>
###### Abstract.

Joint analysis of multi-omic single-cell data across cohorts has significantly
enhanced the comprehensive analysis of cellular processes. However, most of
the existing approaches for this purpose require access to samples with
complete modality availability, which is impractical in many real-world
scenarios. In this paper, we propose $SC^{5}$ (Single-Cell Cross-Cohort Cross-
Category) integration, a novel framework that learns unified cell
representations under domain shift without requiring full-modality reference
samples. Our generative approach learns rich cross-modal and cross-domain
relationships that enable imputation of these missing modalities. Through
experiments on real-world multi-omic datasets, we demonstrate that $SC^{5}$
offers a robust solution to single-cell tasks such as cell type clustering,
cell type classification, and feature imputation.

single-cell, multi-omics

††conference: ; ;

##  1\. Introduction

Advancements in multi-omics single-cell high-throughput sequencing
technologies have introduced novel approaches to integrate rich information
across modalities. Techniques such as single-cell RNA sequencing (scRNA-seq)
combined with the assay for transposase-accessible chromatin using sequencing
(ATAC-seq) allow for the concurrent analysis of the transcriptome and
chromatin accessibility within the same cell(Buenrostro et al., 2015).
Additionally, Cellular Indexing of Transcriptomes and Epitopes by Sequencing
(CITE-seq) enables the measurement of surface protein and transcriptome data
within individual cells using oligonucleotide-tagged antibodies(Stoeckius et
al., 2017). By integrating data from diverse omics, including the
transcriptome, proteome, and epigenome, we can achieve a more comprehensive
understanding of genome regulation from various perspectives(Xu et al., 2022;
Gayoso et al., 2021; Zhou et al., 2023; Ashuach et al., 2023; Cao et al.,
2022; Demetci et al., 2022; Baysoy et al., 2023).

Despite the great efforts on single-cell integration, integrating multi-omic
data across cohorts from different datasets remains a relatively underexplored
problem. To make it consistent and more generalizable, we treat each cohort
from each dataset as a specific domain and each omic as a specific modality.
This problem becomes even more challenging due to the incompleteness of omics
in some cohorts, as shown in Figure 1. This setting is also referred to as
cross-cohort cross-category ($C^{4}$) learning (Rajendran et al., 2024). Most
current single-cell integration methods either lack mechanisms to learn from
incomplete multi-omic datasets or they assume the existence of reference
samples with complete modalities(Zhou et al., 2023; Braams, 1991; Gong et al.,
2021). When some modalities are totally missing in some/all datasets, the lack
of mechanisms to address this gap overlooks the opportunity to leverage the
rich cross-modality interactions present in paired samples. Similarly, when
real-world datasets lack reference samples with full modality
coverage(Rajendran et al., 2023; Lotfollahi et al., 2022) (Case 1 in Figure
1), the presumption of a complete reference dataset becomes untenable.
Moreover, the discrepancy of the feature distributions across datasets, which
is also called batch effect in single-cell settings, is a critical issue for
single-cell integration to be addressed (Luecken et al., 2022). In this work,
we consider the setting of integrating multi-omics data under batch effect
across cohorts where modalities may be entirely missing within a domain,
departing from previous notions of modality scarcity.

Figure 1. Problem overview. Left: Different cases of modality availability.
Right: Different cases of modality imputation.

We propose $SC^{5}$ (Single-Cell Cross-Cohort Cross-Category) to jointly
analyze single-cell data across different domains under missing modality
settings (Cases 2 and 3 in Figure 1). At a high level, our framework models
latent topics underlying single-cell multi-modal features in each domain. We
implement a variational autoencoder that learns modality-dependent and domain-
dependent topic feature distributions in an unsupervised manner. In
particular, we model the data as being generated by a combination of domain-
dependent and modality-dependent latent factors. These domain-dependent
factors describe shared information across modalities within each domain. On
the other hand, the modality-dependent factors describe shared information
across domains with common modalities. This information sharing across
modalities and domains enables feature learning that is robust in the setting
where full-modality samples are unavailable. $SC^{5}$ can transfer rich
multimodal features between domains and impute missing modalities using
observed modalities. The contributions of this proposed framework are:

  1. (1)

Modality-invariant feature learning without requiring full modality
availability.

  2. (2)

Adaptation of disparate domain-specific, multi-modal feature distributions
under missing modality settings.

  3. (3)

Imputation of modality features that are entirely unseen in a domain.

Through experiments on real-world datasets, we show that $SC^{5}$ can
effectively learn cross-modality relationships and thus infer more
representative single-cell embeddings even under missing modality settings.
With both qualitative analysis and quantitive metrics, we demonstrate that our
approach outperforms the baselines in different tasks. Code is available at
https://github.com/anonsc5kdd/sc5.

##  2\. Related Works

###  2.1. Integration with Missing Modalities

Recent methods in multi-modal integration under missing modality settings
perform imputation as a pre-processing step. SMIL (Ma et al., 2021) addresses
bias from modality scarcity by performing missing modality imputation before
learning joint embeddings. MultiModN (Swamy et al., 2023) proposes a
multimodal fusion architecture that maintains robustness to missing modalities
by sequentially encoding each modality, skipping over missing samples. Ding
et. al (Ding et al., 2015) use knowledge transfer to produce embeddings for a
missing modality in a target domain by using an available common modality to
learn low-rank factors. Konwer et. al (Konwer et al., 2023) propose an
adversarial architecture that learns modality-agnostic embeddings by using a
discriminator to predict the absence/presence of modalities. However, their
approach requires their discriminator to be trained over $2^{M}-1$ possible
missing modality combinations for $M$ modalities.

Another line of recent methods apply multi-modal integration for single-cell
data under missing modality settings but require full-modality reference
samples. Multigrate (Lotfollahi et al., 2022) uses a product-of-experts
encoder to integrate available modalities and a condition vector to
incorporate batch information. However, the conditional autoencoder has
limited interpretability due to the neural network decoder. Cobolt (Gong et
al., 2021) is a hierarchical Bayesian generative model that handles missing
modalities by estimating posteriors for each unique modality availability
setting. scMoMAT (Zhang et al., 2023) pre-computes missing features before
using matrix factorization in the integration step. moETM (Zhou et al., 2023)
deals with missing modalities by training a model to reconstruct one modality
from another. Most of these methods require samples with full modality
availability at training. UINMF (Kriebel and Welch, 2022) can integrate
features from single-cell datasets containing shared and unshared features
using non-negative matrix factorization. However, it is not designed to
jointly leverage multiple modalities from the same cell, integrate across
cohorts, or perform missing modality imputation.

###  2.2. Integration under Batch Effect

There are some existing works on addressing batch effect for single-cell
integration (Haghverdi et al., 2018; Li et al., 2020). Li et al. propose an
unsupervised deep embedding algorithm to remove batch effect by clustering (Li
et al., 2020). Haghverdi et al. propose to use matching mutual nearest
neighbors for batch effect removal. However, these techniques to address batch
effect are mostly specific to the structure of integration model. Furthermore,
they assume access to full-modality data.

##  3\. Problem Settings

###  3.1. Notations

We consider the single-cell data from $D$ different domains (i.e. cohorts,
hospitals) and $M$ modalities which might be measured in at least one of these
domains. For domain $d$, there are $N^{(d)}$ cells measured on available
modalities $\mathbf{m}^{(d)}$ ($\mathbf{m}^{(d)}\subseteq\\{1,...,M\\}$ ). We
denote $M_{d}:=|\mathbf{m}^{(d)}|$ as the number of available modalities in
domain $d$. We assume all samples in each domain $d$ have corresponding
measurements for each available modality in $\mathbf{m}^{(d)}$. For a modality
which is not in $\mathbf{m}^{(d)}$, it is completely missing in domain $d$
(there exist no partial samples). The activity of the cell is independently
measured by each modality $m$ ($m\in\mathbf{m}^{(d)}$) by a $V^{(m)}$-dim
feature vector $\mathbf{X}^{(d,m)}$, where each feature corresponds to unique
modality-specific measurements. $\mathbf{x}^{(d,m)}_{n}$ is the value of
$\mathbf{X}^{(d,m)}$ for the $n$-th sample in domain $d$. To simplify our
notation, we will write $\\{\mathbf{x}^{(d,m)}_{n}\\}_{m\in\mathbf{m}_{d}}$ as
$\mathbf{x}^{(d,\cdot)}_{n}$ without causing any confusion.

Table 1. Variable Definitions Variable |  Definition
---|---
$D\in\mathbb{Z}^{+}$ |  # domains (e.g. donors/cohorts)
$M\in\mathbb{Z}^{+}$ |  # modalities in all domains
$\mathbf{m}^{(d)}$ |  The set of modalities in domain $d$
$M^{(d)}$ |  Number of modalities in domain $d$
$N^{(d)}\in\mathbb{Z}^{+}$ |  # cells in domain $d$
$V^{(m)}\in\mathbb{Z}^{+}$ |  # modality-specific features
$L\in\mathbb{Z}^{+}$ |  Hidden feature dimension
$K\in\mathbb{Z}^{+}$ |  # cross-domain, cross-modality topics
$\mathbf{X}^{(d,m)}$ |  Single-cell features
$\boldsymbol{\delta}^{(d,\cdot)}$ |  Cell topic embedding
$\theta^{(d)}$ |  Cell-topic proportion
$\boldsymbol{\alpha}\in\mathbb{R}^{K\times L}$ |  Global topic embedding
$\boldsymbol{\beta}^{(d,\cdot)}\in\mathbb{R}^{K\times L}$ |  Domain-specific topic variation
$\boldsymbol{\omega}^{(d,\cdot)}\in\mathbb{R}^{K\times L}$ |  Domain-specific topic embedding
$\boldsymbol{\rho}^{(\cdot,m)}\in\mathbb{R}^{L\times V^{(m)}}$ |  Cross-domain & modality-specific feature embedding
$\boldsymbol{\lambda}^{(d,m)}\in\mathbb{R}^{1\times V^{(m)}}$ |  Sample-specific noise parameter

###  3.2. Task

####  3.2.1. Learning Integrated Representations

For a given cell $\mathbf{x}^{(d,\cdot)}_{n}$, we aim to integrate features
across the available modalities by learning modality-invariant features
$\mathbf{z}^{(d,\cdot)}_{n}$ describing its common cell states.

####  3.2.2. Missing modality imputation

After obtaining the latent variable $\mathbf{Z}^{(d,\cdot)}$ from the observed
modalities $\\{\mathbf{X}^{(d,m)}\\}_{m\in\mathbf{m}^{(d)}}$, we further aim
to impute features $\hat{\mathbf{X}}^{(i,\bar{m})}$ for a missing modality
$\bar{m}$ in domain $d\in[1,\cdots,D]$.

##  4\. Method

We propose $SC^{5}$, a novel framework to learn integrated cell embeddings
under data heterogeneity by performing information sharing across modalities
and domains. We represent cellular processes as being generated by ‘topics’ as
in (Zhou et al., 2023) and implement variational inference to learn the
generative parameters and integrated features. To model the discrepancy across
different domains, we enable the framework to learn domain-specific topics.
Moreover, our design supports imputation of modalities that are entirely
missing from a domain. In this section, we first formulate the data generative
process from a topic modeling perspective in Section 4.1. Next, we describe
the generative model for integration using a variational approach in Section
4.2 and how to train this model in Section 4.3. In Section 4.4, we incorporate
cell-cell similarities as an auxiliary signal during multi-modal integration
to preserve modality-specific feature relationships in the embedding space.
Finally, we describe how $SC^{5}$ computes integrated representations and
performs imputation for missing modalities on testing samples in Section 4.5.

###  4.1. Data Generative Process

Figure 2. Graphical illustration of generative model. Figure 3. Model
overview. Top: The encoder integrates the available modalities for each domain
via the product-of-experts (PoE). The decoder reconstructs the modality-
specific features by topic and feature embeddings which capture global,
domain-dependent, and modality-dependent variation. Bottom: Integrated feature
representations are obtained by using global and domain-dependent embeddings.

The molecular activities in each cell $n$ can be measured with $M$ omics (such
as gene expression, protein expression or chromatin accessibility) in a domain
$d\in D$. Profiling these omics leads to $M^{(d)}$ feature vectors
$\\{\mathbf{x}^{(d,m)}_{n}\\}_{m\in\mathbf{m}^{(d)}}$, each with modality-
specific feature dimension $V^{(m)}$. We extend topic modeling, similar to
(Zhou et al., 2023), to perform information-sharing between both modalities
and domains in multi-omic feature integration. We consider each cell as a
”document” written in $M$ ”languages” (modalities) and by $D$ ”authors”
(domains). Each feature from the modality $m\in\mathbf{m}^{(d)}$ is a ”word”
from the $m$-th vocabulary whose vocabulary size is $V^{(m)}$. Each sequencing
read is a ”token” in the document and the abundance of reads mapped to the
same feature is the ”word count” in the ”document”. That is, the value of the
$v$-th feature of $\mathbf{x}^{(d,m)}_{n}$ (denoted as
$\mathbf{x}^{(d,m)}_{n,v}$) is count of the word $v^{(m)}$ in the ”document”.

This generative process is described by Figure 2. We use $K$ latent topics to
describe cells across modalities and domains. We assume these $K$ topics are
shared across domains while their proportion distributions might be different
across domains. For each cell (indexed by $n$) in domain $d$, its topic
proportion $\boldsymbol{\theta}^{(d,\cdot)}_{n}$ is sampled from a logistic
normal distribution:

(1)
$\boldsymbol{\delta}^{(d,\cdot)}_{n}\sim\mathcal{N}(\mathbf{M}^{{(d,\cdot)}},\mathbf{C}^{(d,\cdot)}),\boldsymbol{\theta}^{(d,\cdot)}_{n}=\text{softmax}(\boldsymbol{\delta}^{(d,\cdot)}_{n})$

where $\mathbf{M^{(d,\cdot)}},\mathbf{C^{(d,\cdot)}}\in\mathbb{R}^{K}$ are
learnable domain-specific parameters corresponding to the mean vector and
diagonal covariance matrix.

After obtaining $\boldsymbol{\theta}^{(d,\cdot)}_{n}$, the expected rate of
observing different ”words” in cell $n$ is parameterized by the following
equation:

(2)
$\boldsymbol{r}^{(d,m)}_{n}=\text{softmax}\left(\boldsymbol{\theta}^{(d,\cdot)}_{n}\boldsymbol{\omega}^{(d,\cdot)}\boldsymbol{\rho}^{(\cdot,m)}+\boldsymbol{\lambda}^{(d,m)}\right)$

where $\boldsymbol{r}^{(d,m)}_{n}$ is $V^{(m)}$-dim vector, whose $v$-th
element, $\boldsymbol{r}^{(d,m)}_{n,v}$,is the expected rate of observing
feature $v^{(m)}$ in cell $n$.
$\boldsymbol{\omega}^{(d,\cdot)}=\boldsymbol{\alpha}+\boldsymbol{\beta}^{(d,\cdot)}$
is composed of a shared term $\boldsymbol{\alpha}$ and a domain-specific term
$\boldsymbol{\beta}^{(d,\cdot)}$. The shared term $\boldsymbol{\alpha}$ is a
unified topic embedding to describe patterns across modalities and domains,
while the additive term $\boldsymbol{\beta}^{(d,\cdot)}$ flexibly allows
domain-specific variation in the topic features.
$\boldsymbol{\rho}^{(\cdot,m)}\in\mathbb{R}^{L\times V^{(m)}}$ is the cross-
domain feature embedding while
$\boldsymbol{\lambda}^{(d,m)}\in\mathbb{R}^{1\times V^{(m)}}$ is the scalar
noise.

Finally, for each ”read” of the $n$-th cell (”document”) from the $m$-th
modality and the $d$-th domain, a feature index is drawn from a categorical
distribution represented by $r^{(d,m)}_{n}$. And $\mathbf{x}^{(d,m)}_{n,v}$ is
the count of the word $v^{(m)}$. And express the data likelihood in terms of
the read count :

(3)
$p(\boldsymbol{x}^{(d,m)}_{n}|\boldsymbol{r}^{(d,m)}_{n})=\prod_{v=1}^{V^{(m)}}[\boldsymbol{r}^{(d,m)}_{n,v}]^{\boldsymbol{x}^{(d,m)}_{n,v}}$

Based on the above generative process, in order to obtain modality-invariant
representation $\mathbf{z}^{(d,\cdot)}_{n}$ for the cell
$\mathbf{x}^{(d,\cdot)}_{n}$ , we use topic embedding
$\boldsymbol{\omega}^{(d,\cdot)}$ and topic assignments
$\boldsymbol{\theta}^{(d,\cdot)}_{n}$:

(4)
$\mathbf{z}^{(d,\cdot)}_{n}=\boldsymbol{\theta}^{(d,\cdot)}_{n}\boldsymbol{\omega}^{(d,\cdot)}$

###  4.2. Model Design

Based on the data generative process in Section 4.1, we design a variational
model as illustrated in Figure 3. We denote the parameter set of the model as
$\Theta$. We aim to optimize $\Theta$ by maximizing the data log-likelihood
over all cells across domains:
$\mathcal{L}=\sum_{d=1}^{D}\sum_{n=1}^{N^{(d)}}\mathcal{L}^{d}_{n}$, where
$\mathcal{L}^{d}_{n}$ is:

(5) $\displaystyle\mathcal{L}_{n}^{d}$ $\displaystyle=\log
p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\Theta)$
$\displaystyle=\log\int_{\boldsymbol{\delta}}p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta},\Theta)p^{(d,\cdot)}(\boldsymbol{\delta}|\Theta)d\mathbf{\boldsymbol{\delta}}$

where $p^{(d,\cdot)}\left(\boldsymbol{\delta}|\Theta\right)$ and
$p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta},\Theta)$ are
parameterized by $\Theta$. In the following section, we will neglect $\Theta$
and write them as $p^{(d,\cdot)}(\boldsymbol{\delta})$ and
$p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})$ for
simplification. Since $\log p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n})$ is not
tractable, we propose to maximize the evidence lower bound (ELBO) with the
variational posterior
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$:

(6) $\displaystyle\log p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n})$
$\displaystyle=\log\int_{\boldsymbol{\delta}}p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})p^{(d,\cdot)}(\boldsymbol{\delta})d\boldsymbol{\delta}$
$\displaystyle\geq\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\log
p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})$
$\displaystyle\quad\quad-\textsc{KL}(q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})\parallel
p^{(d,\cdot)}(\boldsymbol{\delta}))$
$\displaystyle=\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\sum_{m\in\mathbf{m}^{(d)}}\log
p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$
$\displaystyle\quad\quad-\textsc{KL}(q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})\parallel
p^{(d,\cdot)}(\boldsymbol{\delta}))$

where we name such ELBO term for $\log
p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n})$ as $ELBO^{(d,\cdot)}_{n}$. The full
ELBO formulation is given in Appendix B.

For $p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$ in this
equation, we follow the data generative process in Section 4.1 and
parameterize $p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$ by
$\\{\boldsymbol{\alpha},\boldsymbol{\beta}^{(d,\cdot)},\boldsymbol{\rho}^{(\cdot,m)},\boldsymbol{\lambda}^{(d,m)}\\}_{1\leq
d\leq D,m\in\mathbf{m}^{(d)}}$. This is corresponding to the ”decoder” module
in Figure 3.

For $q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$, we
estimate this posterior distributions as a normal distribution. To infer
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$ from multi-
modality $\mathbf{x}^{(d,\cdot)}_{n}$, we model it as a product of normal
distributions (Wu and Goodman, 2018) across modality-specific topic assignment
distributions within the domain. Particularly, we first model
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,m)}_{n})$ as a normal
distribution for each $m$ and compute its mean and covariance as follows:

(7)
$[\boldsymbol{\mu}^{(d,m)}_{n},\boldsymbol{\Sigma}^{(d,m)}_{n}]=\text{NNET}^{(\cdot,m)}\left(\tilde{\mathbf{x}}^{(d,m)}_{n};\mathbf{W}^{(\cdot,m)}\right)$

where $\tilde{\mathbf{x}}^{(d,m)}_{n}$ is normalized counts for each feature
as the raw count of the feature divided by the total counts of modality $m$ in
cell $n$. $\text{NNET}^{(\cdot,m)}$ is a modality-specific neural network
shared across domains with parameters denoted as $\mathbf{W}^{(\cdot,m)}$.

After obtaining
$[\boldsymbol{\mu}^{(d,m)}_{n},\boldsymbol{\Sigma}^{(d,m)}_{n}]$ for each
available modality,
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$ is
estimated as the product of these normal distributions. The mean and
covariance of this joint normal for domain $d$ are computed as:

(8)
${\boldsymbol{\mu}^{*}}^{(d,\cdot)}_{n}=\frac{\sum_{m\in\mathbf{m}^{(d)}}\boldsymbol{\mu}^{(d,m)}_{n}\boldsymbol{\Sigma}^{(d,m)}_{n}}{1+\sum_{m\in\mathbf{m}^{(d)}}\boldsymbol{\Sigma}^{(d,m)}_{n}},{\boldsymbol{\Sigma}^{*}}^{(d,\cdot)}_{n}=\frac{\prod_{m\in\mathbf{m}^{(d)}}\Sigma^{(d,m)}_{n}}{1+\sum_{m\in\mathbf{m}^{(d)}}\Sigma^{(d,m)}_{n}}$

The neural networks and the product-of-experts operation correspond to the
”encoder” module in Figure 3. Notably, our framework can flexibly handle
variable modality availability in computing the topic assignment distribution,
as the product is taken over available modalities $\mathbf{m}^{(d)}$ in domain
$d$.

###  4.3. Training Procedure

With the model introduced in Section 4.2, the evidence lower bound (ELBO) in
Equation 6 can be calculated as follows:

We resample $\boldsymbol{\delta}$ from
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$ and then
compute the term
$\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\sum_{m\in\mathbf{m}^{(d)}}\log
p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$ by
reparameterization trick (Kingma and Welling, 2013), where $\log
p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$ can be
calculated via the procedure and equations in Section 4.2.

The KL-divergence in Equation (6) between the univariate Gaussian posterior
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})$ and
univariate Gaussian learnable priors $p^{(d,\cdot)}(\boldsymbol{\delta})$ for
a domain $d$, which we express as
$\textsc{KL}\left(q(\boldsymbol{\delta})\parallel
p_{\theta}(\boldsymbol{\delta})\right)$ for simplicity, has a closed-form
solution:

(9) $\displaystyle\textsc{KL}\left(q(\boldsymbol{\delta})\parallel
p(\boldsymbol{\delta})\right)=\log\left(\frac{\sqrt{\boldsymbol{\Sigma}_{p}}}{\sqrt{\boldsymbol{\Sigma}_{q}}}\right)+\frac{\boldsymbol{\Sigma}_{q}+(\boldsymbol{\mu}_{q}-\boldsymbol{\mu}_{p})^{2}}{2\boldsymbol{\Sigma}_{p}}-\frac{1}{2}$

where the full derivation is given in Appendix B.

All the parameters $\Theta$ of the model, including the decoder weights
$\\{\boldsymbol{\alpha},\boldsymbol{\beta}^{(d,\cdot)},\boldsymbol{\rho}^{(\cdot,m)},\boldsymbol{\lambda}^{(d,m)}\\}_{1\leq
d\leq D,m\in\mathbf{m}^{(d)}}$ and the encoder weights
$\\{\mathbf{W}^{(\cdot,m)}\\}_{1\leq m\leq M}$ are jointly optimized by
minimizing loss $\mathcal{L}$ across $D$ domains:

(10)
$\Theta\leftarrow\operatorname*{arg\,min}_{\Theta}\sum_{d=1}^{D}\mathcal{L}^{(d,\cdot)}$

where $\mathcal{L}^{(d,\cdot)}$ is computed as follows:

(11)
$\mathcal{L}^{(d,\cdot)}=-\frac{1}{N^{(d,\cdot)}}\sum_{n=1}^{N^{(d,\cdot)}}ELBO^{(d,\cdot)}_{n}+\lambda_{\beta}\|\boldsymbol{\beta}^{(d,\cdot)}\|_{2}$

Here we also constrain the norm of $\boldsymbol{\beta}$ to ensure that
$\alpha$ is representative of global topics and $\lambda_{\beta}$ is the
weight.

###  4.4. Neighborhood Contrastive Loss

One issue that may arise in the missing modality setting is the bias to
certain modalities. For example, the embedding may rely on features from
modalities that are available across the majority of the domains, leading to
poor performance in some domains where this modality is unavailable.
Therefore, we ensure that available modalities in each domain $d$ are equally
represented by regularizing the learned embedding
$\boldsymbol{\delta}^{(d,\cdot)}$ to preserve cell-cell relations for all
cells in domain $d$.

For two different cells $\mathbf{x}^{(d,m)}_{n},\mathbf{x}^{(d,m)}_{i}$ from
the same domain $d$ and modality $m\in\mathbf{m}^{(d)}$, we measure their
distance as the Euclidean distance between $\mathbf{x}^{(d,m)}_{n}$ and
$\mathbf{x}^{(d,m)}_{n}$. Then the k-nearest neighbor samples for a given
sample $\mathbf{x}^{(d,\cdot)}_{n}$ on modality $m$ can be computed with this
distance metric, which is denoted as $\mathbf{nn}^{(d,m)}_{k}(n)$. Note that
$\mathbf{nn}^{(d,m)}_{k}(n)$ might be different across modalities.

The neighborhood contrastive loss (NCL) objective maximizes the similarity
between cell topic features from cells and their nearest neighboring cells in
domain $d$ and modality $m$, while minimizing the embedding similarity of
cells from separate neighborhoods:

(12)
$\mathcal{L}^{(d,\cdot)}_{cont}(n)=-\textstyle\sum_{m\in\mathbf{m}^{(d)}}\sum_{i\in\mathbf{nn}^{(d,m)}_{k}(n)}\log\left(\frac{\sigma(\boldsymbol{\delta}_{n}^{(d,m)},\boldsymbol{\delta}_{i}^{(d,m)})}{\sum_{j\neq
n}\sigma(\boldsymbol{\delta}_{n}^{(d,m)},\boldsymbol{\delta}_{j}^{(d,m)})}\right)$

where embedding similarity
$\sigma(\boldsymbol{\delta}_{n}^{(d,\cdot)},\boldsymbol{\delta}_{i}^{(d,\cdot)})$
between cells $n$ and $i$ is computed as:

(13)
$\sigma(\boldsymbol{\delta}_{n}^{(d,m)},\boldsymbol{\delta}_{i}^{(d,m)})=\text{exp}\left(\frac{\langle\boldsymbol{\delta}_{n}^{(d,m)},\boldsymbol{\delta}_{i}^{(d,m)}\rangle}{\kappa\|\boldsymbol{\delta}_{n}^{(d,m)}\|\cdot\|\boldsymbol{\delta}_{i}^{(d,m)}\|}\right)$

for a positive constant $\kappa$.

Once we incorporate this contrastive loss, the final objective function
becomes:

(14)
$\mathcal{L}^{(d,\cdot)}=-\frac{1}{N^{(d,\cdot)}}\sum_{n=1}^{N^{(d,\cdot)}}ELBO^{(d,\cdot)}_{n}+\lambda_{\beta}\|\boldsymbol{\beta}^{(d,\cdot)}\|_{2}+\frac{1}{N^{(d)}}\sum_{n=1}^{N^{(d)}}\mathcal{L}^{(d,\cdot)}_{cont}(n)$

###  4.5. Model Inference

For any sample $\mathbf{x}^{(d,\cdot)}_{*}$ in domain $d$, we can compute its
integrated feature representation and impute missing modalities.

#### Integration Representation

We first feed the features $\mathbf{x}^{(d,\cdot)}_{*}$ to the encoders and
aggregate the result as in Eq. (7)-(8) to obtain the cell topic embedding
$\boldsymbol{\delta}_{n}^{(d,\cdot)}$. We then compute the topic mixture
probability $\boldsymbol{\theta}_{n}^{(d,\cdot)}$ by applying a softmax
transformation on $\boldsymbol{\delta}_{n}^{(d,\cdot)}$ to obtain the final
integrated feature representation:

(15)
$\mathbf{z}^{(d,\cdot)}_{*}=\boldsymbol{\theta}_{n}^{(d,\cdot)}(\boldsymbol{\alpha}+\boldsymbol{\beta}^{(d,\cdot)})$

#### Missing Modality Imputation

Our framework supports missing modality imputation, even if a modality is
missing across all samples in a domain. We can impute the missing features
$\hat{\mathbf{X}}^{(i,\bar{m})}$ for any $\bar{m}\notin\mathbf{m}^{(d)}$ as
follows:

(16)
$\hat{\mathbf{x}}^{(d,\bar{m})}_{n}=\text{softmax}\left(\boldsymbol{\theta}_{n}^{(d,\cdot)}(\boldsymbol{\alpha}+\boldsymbol{\beta}^{(d,\cdot)})\boldsymbol{\rho}^{(\cdot,\bar{m})}\right)$

##  5\. Experiments

###  5.1. Dataset and Experiment Settings

We use the 2021 NeurIPS single-cell challenge dataset with inherent missing
modality settings (unmeasured modalities) and simulated missing modality
settings (masked modalities) to evaluate the methods (Luecken et al., 2021).
This dataset contains bone marrow mononuclear cells profiled Multiome (gene
expression & chromatin accessibility) (Buenrostro et al., 2015) and CITE-seq
(gene expression & protein abundance) (Stoeckius et al., 2017), respectively.

The NeurIPS dataset is suitable for evaluating integration robustness under
missing modality settings as it inherently contains missing modalities
(chromatin accessibility is missing from CITE-seq samples while protein
abundance is missing from Multiome samples). In order to evaluate the
imputation performance, we must simulate missing modalities where we have
access to ground-truth features.

We consider three different scenarios with the NeurIPS single-cell challenge
dataset, as shown in Figure 4. We select half of the domains in each
experimental scenario as imputation target domains, where we simulate all
imputation possibilities using our ground-truth features. We formulate these
scenarios as follows:

Figure 4. Domains and modality availability under each experiment scenario.

  1. (1)

CITE-seq. We consider the data of CITE-seq from four different sites. We treat
each site as one domain. In this setting, two modalities (gene expression
(GEX) & protein abundance (ADT)) are considered. We assume the availability of
both modalities in domains 1 and 2, protein abundance in domain 3 and gene
expression in domain 4.

  2. (2)

Multiome. We consider the data of Multiome from four different sites. We treat
each site as one domain. In this setting, two modalities (gene expression
(GEX) & chromatin accessibility (ATAC)) are considered. We assume the
availability of both modalities in domains 1 and 2, gene expression in domain
3 and protein abundance in domain 4.

  3. (3)

Combine. We combine the above data of CITE-seq and Multiome. We treat each
site as one domain. In this setting, three modalities (GEX, ADT, ATAC) are
considered. We assume the availability of gene expression & protein abundance
in domains 1 and 2, gene expression & chromatin accessibility in domains 3 and
4. For domains 5-8, only one modality is available.

Table 2. Results of embedding-based clustering. Two metrics ARI and NMI are used to evaluate the performance. Averaged performance over different imputation target domains are reported. | Method | Cite-seq | Multiome | Combine
---|---|---|---|---
ARI | CCA | $0.497$ | $0.432$ | $0.360$
moETM | $0.546$ | $0.561$ | $0.522$
$\rm SC^{5}$ | $\mathbf{0.590}$ | $\mathbf{0.578}$ | $\mathbf{0.572}$
NMI | CCA | $0.694$ | $0.307$ | $0.523$
moETM | $0.731$ | $0.632$ | $0.693$
$\rm SC^{5}$ | $\mathbf{0.746}$ | $\mathbf{0.641}$ | $\mathbf{0.706}$
Table 3. Results of cell type classification. Accuracy is used to evaluate the performance. Averaged performance over different imputation target domains are reported. Method | Cite-seq | Multiome | Combine
---|---|---|---
CCA | $0.657$ | $0.589$ | $0.539$
moETM | $\mathbf{0.796}$ | $0.637$ | $0.748$
$\rm SC^{5}$ | $0.794$ | $\mathbf{0.671}$ | $\mathbf{0.751}$
Figure 5. Distribution of topic assignment scores in the $Combine$ setting.
The top 20% of assigned topics are selected across 500 sampled cells from 2
domains. Color intensity values correspond to the cell-topic feature value
before normalization into a topic mixture probability. Boxed topic features
correspond to cell topic features with strong association to unique cell
types.

###  5.2. Evaluation Metric

We evaluate the performance of all methods from two aspects:

####  5.2.1. Integrated feature quality from dataset combination

We follow previous works (Zhou et al., 2023) and assess the self-supervised
multi-modal integration using the following downstream tasks:

  1. (1)

Self-supervised cell type clustering. We apply k-means clustering on the
integrated multi-modal embeddings to cluster cells into biologically
meaningful clusters that are useful for identifying cell types. We evaluate
clustering quality using adjusted Rand index (ARI) and normalized mutual
information (NMI) between k-means cell clusters and cell type labels.

  2. (2)

Supervised cell type classification. We also use the integrated embeddings to
train a downstream cell type classifier. We train the classifier on the
integrated cell embeddings and evaluate the performance on a hold-out set
using cross-entropy.

####  5.2.2. Imputation quality under missing modality settings

We assess the imputation quality under simulated missing modality settings
using Pearson correlations between the true and imputed features.

###  5.3. Baselines

We compare $SC^{5}$ with the following baseline methods by extending them to
our proposed setting where no full-modality reference is available. As
previous single-cell integration methods are unable to integrate features
across cohorts/domains with conflicting modality availability in this proposed
setting, we use modified baseline methods that are suited for our task.

  1. (1)

”Missing-aware” CCA is an extension of canonical correlation analysis (CCA),
which has been used for single-cell integration (Butler et al., 2018), that
produces multi-modal embeddings using SVD across samples while masking missing
samples.

  2. (2)

”Missing-aware” moETM is an extension of a multi-omic topic embedding
framework (moETM) (Zhou et al., 2023) that learns multi-modal embeddings using
a variational autoencoder across samples while masking missing samples.

Both of these modified baselines can be used to learn integrated
representation after adjusting for missing modalities, while ”missing-aware”
moETM supports imputing missing modalities. Therefore, we compare $SC^{5}$
with both baselines on cell type clustering and classification and compare
$SC^{5}$ with moETM on missing modality imputation.

###  5.4. Cell Type Clustering and Classification

The experimental results on cell type clustering and classification are shown
in Tables 2 and 3. Topic-modelling based methods outperforms CCA by large
margin, which indicates that topic modelling is more suitable for describing
cell activities. $SC^{5}$ outperforms the baselines in all but one of the
settings. As NMI is best suited for evaluating imbalanced clusters, which we
observe in practice (Figure 6), this metric is higher than ARI across all
methods.

We further demonstrate that $SC^{5}$ learns informative topic features
$\boldsymbol{\omega}$ for identifying cell types by visualizing the topic
distribution. In Figure 5, we visualize the distribution of the top 20% of
inferred topics shared among 500 sampled cells across 2 domains. We annotate
cell topic features that have high topic scores for unique cell types.
Notably, $SC^{5}$ infers common topics from separate domains that are
associated with the same cell type. For example, topic 7 is associated with
high topic intensity for monocyte cells in both domain 6 (annotated as
‘monocyte’) and in domain 7 (annotated by the subtypes ‘CD16+ Mono’ and ‘CD14+
Mono’). Additionally, topic 61 is associated with ‘B’ cells in both domains.

In Appendix C, we compare the learned topic embeddings between domains and
find that the integrated features are flexible to domain-dependent
heterogeneity. These results suggest that $SC^{5}$ learns topic features that
are informative while being flexible to domain-specific heterogeneity.

Figure 6. UMAP visualization of cell clustering on imputed modalities under the $Combine$ setting. Each row corresponds the imputation results in one domain. In each subfigure, individual data points represent cells, while different colors are used to denote different cell types. Table 4. Results of missing modality imputation. Pearson correlation between true and imputed features is used to evaluate the performance. Averaged performance over different imputation target domains are reported. Method | Cite-seq | Multiome | Combine
---|---|---|---
moETM | $0.465$ | $0.184$ | $0.326$
$\rm SC^{5}$ | $\mathbf{0.671}$ | $\mathbf{0.278}$ | $\mathbf{0.444}$

###  5.5. Missing Modality Imputation

The experimental results on missing modality imputation are displayed in Table
4. From Table 4, $SC^{5}$ shows a significant performance improvement over the
baselines.

By the results, we demonstrate that $SC^{5}$ can impute modalities that are
completely missing in a domain, even when a full multimodal reference is
unavailable during training. By learning complex relationships between samples
from varying modalities and domains, our method can impute features that match
the modality-specific and domain-specific distributions of the available data.
Unlike previous works, $SC^{5}$ can impute features for a modality that is
completely unseen in a target dataset, as long as the target modality is
available in at least one domain during training.

Table 5. Ablation study on neighborhood contrastive loss. Averaged performance over different imputation target domains are reported. | Method | Cite-seq | Multiome | Combine
---|---|---|---|---
ARI |  $\rm SC^{5}$ wo. NCL | $0.529$ | $0.523$ | $0.554$
$\rm SC^{5}$ | $\mathbf{0.590}$ | $\mathbf{0.578}$ | $\mathbf{0.572}$
NMI |  $\rm SC^{5}$ wo. NCL | $0.738$ | $\mathbf{0.643}$ | $0.699$
$\rm SC^{5}$ | $\mathbf{0.746}$ | $0.641$ | $\mathbf{0.706}$
Classification |  $\rm SC^{5}$ wo. NCL | $\mathbf{0.804}$ | $\mathbf{0.699}$ | $\mathbf{0.757}$
$\rm SC^{5}$ | $0.794$ | $0.671$ | $0.751$
Imputation |  $\rm SC^{5}$ wo. NCL | $\mathbf{0.694}$ | $\mathbf{0.286}$ | $0.435$
$\rm SC^{5}$ | $0.671$ | $0.278$ | $\mathbf{0.444}$

We verify $SC^{5}$’s imputation capabilities by visualizing the imputed
features using Uniform Manifold Approximation and Projection (UMAP) (McInnes
et al., 2018) in Figure 6. Our $SC^{5}$ variants impute features that preserve
finer-grained cell type clusters in contrast to moETM which does not
explicitly model domain-dependent features. In particular, ”CD4+ T cells”
associated with immune response are have greater disentaglement when imputing
protein abundance from gene expression (rows 1, 2). When imputing chromatin
accessibility from gene expression (row 3), $SC^{5}$ better separates
”Erythoroblast”, ”Proerythroblast”, ”Normoblast”, and ”MK/E prog” cells
associated with red blood cell production. We note that on the most difficult
imputation task (imputing dense gene expression features from sparse chromatin
accessibility features), the imputed cell features share similar patterns,
though $SC^{5}$ produces embeddings that are more compact compared to the
baseline. Our results demonstrate that $SC^{5}$ learns robust domain-dependent
and modality factors that preserve cell type clustering patterns when imputing
masked features.

###  5.6. Ablation Study

We now experimentally validate the effectiveness of neighborhood contrastive
loss (NCL) introduced in Section 4.4. We conduct an ablation study by
comparing $SC^{5}$ with its variant obtained by removing the neighborhood
contrastive loss term. As the results show in Table 5, adding the neighborhood
contrastive loss (NCL) term consistently improves the clustering performance.
Notably, the version with neighborhood contrastive loss can achieve the best
or comparable classification and imputation results for all tasks under the
$Combine$ setting.

##  6\. Conclusion

In this paper, we propose a novel framework for joint analysis of single-cell
data across cohorts (domains) with incomplete modalities. The proposed
framework utilizes topic modeling to capture the data generative process,
incorporating shared and domain-specific parameters. Through experiments on
real-world datasets, we demonstrate that our proposed framework can learn
integrated representations and impute missing modalities more efficiently than
the baselines.

## References

  * (1)
  * Ashuach et al. (2023) Tal Ashuach, Mariano I Gabitto, Rohan V Koodli, Giuseppe-Antonio Saldi, Michael I Jordan, and Nir Yosef. 2023.  MultiVI: deep generative model for the integration of multimodal data.  _Nature Methods_ 20, 8 (2023), 1222–1231.
  * Baysoy et al. (2023) Alev Baysoy, Zhiliang Bai, Rahul Satija, and Rong Fan. 2023.  The technological landscape and applications of single-cell multi-omics.  _Nature Reviews Molecular Cell Biology_ (2023), 1–19.
  * Braams (1991) Johannes Braams. 1991.  Babel, a Multilingual Style-Option System for Use with LaTeX’s Standard Document Styles.  _TUGboat_ 12, 2 (June 1991), 291–301.
  * Buenrostro et al. (2015) Jason D Buenrostro, Beijing Wu, Ulrike M Litzenburger, Dave Ruff, Michael L Gonzales, Michael P Snyder, Howard Y Chang, and William J Greenleaf. 2015.  Single-cell chromatin accessibility reveals principles of regulatory variation.  _Nature_ 523, 7561 (2015), 486–490.
  * Butler et al. (2018) Andrew Butler, Paul Hoffman, Peter Smibert, Efthymia Papalexi, and Rahul Satija. 2018.  Integrating single-cell transcriptomic data across different conditions, technologies, and species.  _Nature biotechnology_ 36, 5 (2018), 411–420.
  * Cao et al. (2022) Kai Cao, Yiguang Hong, and Lin Wan. 2022.  Manifold alignment for heterogeneous single-cell multi-omics data integration using Pamona.  _Bioinformatics_ 38, 1 (2022), 211–219.
  * Demetci et al. (2022) Pinar Demetci, Rebecca Santorella, Manav Chakravarthy, Bjorn Sandstede, and Ritambhara Singh. 2022.  Scotv2: Single-cell multiomic alignment with disproportionate cell-type representation.  _Journal of Computational Biology_ 29, 11 (2022), 1213–1228.
  * Ding et al. (2015) Zhengming Ding, Ming Shao, and Yun Fu. 2015.  Missing Modality Transfer Learning via Latent Low-Rank Constraint.  _IEEE Transactions on Image Processing_ 24, 11 (2015), 4322–4334.  https://doi.org/10.1109/TIP.2015.2462023
  * Gayoso et al. (2021) Adam Gayoso, Zoë Steier, Romain Lopez, Jeffrey Regier, Kristopher L Nazor, Aaron Streets, and Nir Yosef. 2021.  Joint probabilistic modeling of single-cell multi-omic data with totalVI.  _Nature methods_ 18, 3 (2021), 272–282.
  * Gong et al. (2021) Boying Gong, Yun Zhou, and Elizabeth Purdom. 2021.  Cobolt: integrative analysis of multimodal single-cell sequencing data.  _Genome biology_ 22, 1 (2021), 1–21.
  * Haghverdi et al. (2018) Laleh Haghverdi, Aaron TL Lun, Michael D Morgan, and John C Marioni. 2018.  Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors.  _Nature biotechnology_ 36, 5 (2018), 421–427.
  * Kingma and Welling (2013) Diederik P Kingma and Max Welling. 2013.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ (2013).
  * Konwer et al. (2023) Aishik Konwer, Xiaoling Hu, Joseph Bae, Xuan Xu, Chao Chen, and Prateek Prasanna. 2023.  Enhancing modality-agnostic representations via meta-learning for brain tumor segmentation. In _Proceedings of the IEEE/CVF International Conference on Computer Vision_. 21415–21425.
  * Kriebel and Welch (2022) April R Kriebel and Joshua D Welch. 2022.  UINMF performs mosaic integration of single-cell multi-omic datasets using nonnegative matrix factorization.  _Nature communications_ 13, 1 (2022), 780.
  * Li et al. (2020) Xiangjie Li, Kui Wang, Yafei Lyu, Huize Pan, Jingxiao Zhang, Dwight Stambolian, Katalin Susztak, Muredach P Reilly, Gang Hu, and Mingyao Li. 2020.  Deep learning enables accurate clustering with batch effect removal in single-cell RNA-seq analysis.  _Nature communications_ 11, 1 (2020), 2338.
  * Lotfollahi et al. (2022) Mohammad Lotfollahi, Anastasia Litinetskaya, and Fabian J. Theis. 2022.  Multigrate: single-cell multi-omic data integration.  _bioRxiv_ (2022).  https://doi.org/10.1101/2022.03.16.484643 arXiv:https://www.biorxiv.org/content/early/2022/03/17/2022.03.16.484643.full.pdf
  * Luecken et al. (2021) Malte D Luecken, Daniel Bernard Burkhardt, Robrecht Cannoodt, Christopher Lance, Aditi Agrawal, Hananeh Aliee, Ann T Chen, Louise Deconinck, Angela M Detweiler, Alejandro A Granados, et al. 2021\.  A sandbox for prediction and integration of DNA, RNA, and proteins in single cells. In _35th Conference on Neural Information Processing Systems (NeurIPS 2021) Track on Datasets and Benchmarks_.
  * Luecken et al. (2022) Malte D Luecken, Maren Büttner, Kridsadakorn Chaichoompu, Anna Danese, Marta Interlandi, Michaela F Müller, Daniel C Strobl, Luke Zappia, Martin Dugas, Maria Colomé-Tatché, et al. 2022\.  Benchmarking atlas-level data integration in single-cell genomics.  _Nature methods_ 19, 1 (2022), 41–50.
  * Ma et al. (2021) Mengmeng Ma, Jian Ren, Long Zhao, Sergey Tulyakov, Cathy Wu, and Xi Peng. 2021.  Smil: Multimodal learning with severely missing modality. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 35. 2302–2310.
  * McInnes et al. (2018) Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger. 2018.  UMAP: Uniform Manifold Approximation and Projection.  _Journal of Open Source Software_ 3, 29 (2018), 861.
  * Rajendran et al. (2024) Suraj Rajendran, Weishen Pan, Mert R Sabuncu, Yong Chen, Jiayu Zhou, and Fei Wang. 2024.  Learning across diverse biomedical data modalities and cohorts: Challenges and opportunities for innovation.  _Patterns_ (2024).
  * Rajendran et al. (2023) Suraj Rajendran, Weishen Pan, Mert R Sabuncu, Jiayu Zhou, and Fei Wang. 2023.  Patchwork Learning: A Paradigm Towards Integrative Analysis across Diverse Biomedical Data Sources.  _arXiv preprint arXiv:2305.06217_ (2023).
  * Stoeckius et al. (2017) Marlon Stoeckius, Christoph Hafemeister, William Stephenson, Brian Houck-Loomis, Pratip K Chattopadhyay, Harold Swerdlow, Rahul Satija, and Peter Smibert. 2017.  Simultaneous epitope and transcriptome measurement in single cells.  _Nature methods_ 14, 9 (2017), 865–868.
  * Swamy et al. (2023) Vinitra Swamy, Malika Satayeva, Jibril Frej, Thierry Bossy, Thijs Vogels, Martin Jaggi, Tanja Käser, and Mary-Anne Hartley. 2023.  MultiModN- Multimodal, Multi-Task, Interpretable Modular Networks.  (2023).
  * Wu and Goodman (2018) Mike Wu and Noah Goodman. 2018.  Multimodal generative models for scalable weakly-supervised learning.  _Advances in neural information processing systems_ 31 (2018).
  * Xu et al. (2022) Yang Xu, Priyojit Das, and Rachel Patton McCord. 2022.  SMILE: mutual information learning for integration of single-cell omics data.  _Bioinformatics_ 38, 2 (2022), 476–486.
  * Zhang et al. (2023) Ziqi Zhang, Haoran Sun, Ragunathan Mariappan, Xi Chen, Xinyu Chen, Mika S Jain, Mirjana Efremova, Sarah A Teichmann, Vaibhav Rajan, and Xiuwei Zhang. 2023.  scMoMaT jointly performs single cell mosaic integration and multi-modal bio-marker detection.  _Nature Communications_ 14, 1 (2023), 384.
  * Zhou et al. (2023) Manqi Zhou, Hao Zhang, Zilong Baii, Dylan Mann-Krzisnik, Fei Wang, and Yue Li. 2023.  Single-cell multi-omic topic embedding reveals cell-type-specific and COVID-19 severity-related immune signatures.  _Cell Reports Methods_.  https://pubmed.ncbi.nlm.nih.gov/36778483/

##  Appendix A Data processing

All datasets were processed into the format of cells-by-features matrices. For
each dataset, read counts for each feature were first normalized per cell by
total counts within the same omic. This normalization was executed using the
scanpy. Subsequently, a log1p transformation was applied to the normalized
data.

Following normalization and transformation, we identified highly variable
features by scanpy.pp.highly_variable_genes. In the case of protein features,
all measured surface proteins were included due to the inherently limited
quantity of proteins measured by the scADT-seq assay in comparison to genes or
chromatin regions, and their substantial elucidation of cellular functions.

##  Appendix B Derivations for Loss Functions

###  B.1. Evidence Lower Bound

(17) $\displaystyle\log p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n})$
$\displaystyle=\log\int_{\boldsymbol{\delta}}p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})p^{(d,\cdot)}(\boldsymbol{\delta})d\boldsymbol{\delta}$
$\displaystyle=\log\int_{\boldsymbol{\delta}}q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n}))\frac{p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})p^{(d,\cdot)}(\boldsymbol{\delta})}{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}d\boldsymbol{\delta}$
$\displaystyle=\log\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\frac{p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})p^{(d,\cdot)}(\boldsymbol{\delta})}{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}$
$\displaystyle\geq\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\log\frac{p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})p^{(d,\cdot)}(\boldsymbol{\delta})}{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}$
$\displaystyle=\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\log
p^{(d,\cdot)}(\mathbf{x}^{(d,\cdot)}_{n}\mid\boldsymbol{\delta})$
$\displaystyle\quad\quad-\textsc{KL}(q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})\parallel
p^{(d,\cdot)}(\boldsymbol{\delta}))$
$\displaystyle=\mathbb{E}_{q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})}\sum_{m\in\mathbf{m}^{(d)}}\log
p^{(d,\cdot)}(\mathbf{x}^{(d,m)}_{n}\mid\boldsymbol{\delta})$
$\displaystyle\quad\quad-\textsc{KL}(q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x}^{(d,\cdot)}_{n})\parallel
p^{(d,\cdot)}(\boldsymbol{\delta}))$

###  B.2. KL-Divergence

The KL-divergence in Equation 6 between the univariate Gaussian posterior
$q^{(d,\cdot)}(\boldsymbol{\delta}\mid\mathbf{x})$ and univariate Gaussian
learnable priors $p_{\theta}^{(d,\cdot)}(\boldsymbol{\delta})$ for a domain
$d$, which we express as $\textsc{KL}\left(q(\boldsymbol{\delta})\parallel
p_{\theta}(\boldsymbol{\delta})\right)$ for simplicity, is derived as follows:

(18) $\displaystyle\textsc{KL}\left(q(\boldsymbol{\delta})\parallel
p_{\theta}(\boldsymbol{\delta})\right)$ $\displaystyle=\int
q(\boldsymbol{\delta})\log\frac{q(\boldsymbol{\delta})}{p_{\theta}(\boldsymbol{\delta})}d\boldsymbol{\delta}$
$\displaystyle=\int q(\boldsymbol{\delta})\left[\log
q(\boldsymbol{\delta})-\log
p_{\theta}(\boldsymbol{\delta})\right]d\boldsymbol{\delta}$
$\displaystyle=\int\frac{1}{\sqrt{2\pi\boldsymbol{\Sigma}_{q}}}\text{exp}\left[\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{q})^{2}}{2\boldsymbol{\Sigma}_{q}}\right]$
$\displaystyle\quad\times\left[\log\frac{\sqrt{\boldsymbol{\Sigma}_{p}}}{\sqrt{\boldsymbol{\Sigma}_{q}}}-\left(\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{p})^{2}}{2\boldsymbol{\Sigma}_{p}}+\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{q})^{2}}{2\boldsymbol{\Sigma}_{q}}\right)\right]$
$\displaystyle=\mathbb{E}_{q(\boldsymbol{\delta})}\left[\log\frac{\sqrt{\boldsymbol{\Sigma}_{p}}}{\sqrt{\boldsymbol{\Sigma}_{q}}}-\left(\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{p})^{2}}{2\boldsymbol{\Sigma}_{p}}+\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{q})^{2}}{2\boldsymbol{\Sigma}_{q}}\right)\right]$
$\displaystyle=\log\frac{\sqrt{\boldsymbol{\Sigma}_{p}}}{\sqrt{\boldsymbol{\Sigma}_{q}}}+\mathbb{E}_{q(\boldsymbol{\delta})}\left[\frac{(\boldsymbol{\delta}-\boldsymbol{\mu}_{p})^{2}}{2\boldsymbol{\Sigma}_{p}}\right]-\frac{1}{2}$
$\displaystyle=\log\frac{\sqrt{\boldsymbol{\Sigma}_{p}}}{\sqrt{\boldsymbol{\Sigma}_{q}}}+\frac{\boldsymbol{\Sigma}_{p}+(\boldsymbol{\mu}_{q}-\boldsymbol{\mu}_{p})^{2}}{2\boldsymbol{\Sigma}_{p}}-\frac{1}{2}$

##  Appendix C Topic embedding analysis

We assess $SC^{5}$’s ability to infer latent topic features that are flexible
to domain-specific heterogeneity. In figure 7, we visualize the distribution
of the top 20% of topics inferred by $\rm SC^{5}$ averaged within each test
domain. Notably, our method is able to learn cell topic feature distributions
that are unique to each domain.

Figure 7. Distribution of scores for top 20% of assigned topics averaged
across samples within select domains. Topic features are learned by under the
Combine setting.