## Distance-Preserving Spatial Representations in Genomic Data

Wenbin Zhou and Jin-Hong Du

Abstract -The spatial context of single-cell gene expression data is crucial for many downstream analyses, yet often remains inaccessible due to practical and technical limitations, restricting the utility of such datasets. In this paper, we propose a generic representation learning and transfer learning framework dp-VAE , capable of reconstructing the spatial coordinates associated with the provided gene expression data. Central to our approach is a distance-preserving regularizer integrated into the loss function during training, ensuring the model effectively captures and utilizes spatial context signals from reference datasets. During the inference stage, the produced latent representation of the model can be used to reconstruct or impute the spatial context of the provided gene expression by solving a constrained optimization problem. We also explore the theoretical connections between distance-preserving loss, distortion, and the bi-Lipschitz condition within generative models. Finally, we demonstrate the effectiveness of dp-VAE in different tasks involving training robustness, out-of-sample evaluation, and transfer learning inference applications by testing it over 27 publicly available datasets. This underscores its applicability to a wide range of genomics studies that were previously hindered by the absence of spatial data.

Index Terms -Representation learning, spatial transcriptomics, single-cell RNA sequence, variational inference, and dimensionality reduction.

✦

## 1 INTRODUCTION

U NDERSTANDING the tissue spatial context from which the gene expressions are extracted is crucial to many applications, deepening the understanding of cellular organizations and their interactions within tissues, as well as illuminating biological insights in neuroscience, developmental biology, precise diagnostics, personalized therapies, and a range of disease mechanisms such as cancer [57]. As a result, measurement technologies such as spatial transcriptomics have emerged as a powerful tool in modern biology, bridging the gap between single-cell genomics and the larger physiological processes at play in living tissues.

However, various practical and technical constraints limit the range of spatial context and the granularity that these measurement tools can directly capture. For example, spatial transcriptomics (ST) technologies such as Visium by 10x Genomics and NanoString GeoMx [10] typically do not have single-cell resolution [7], making the extracted spatial coordinates within the tissue space too coarse for many precision-focused downstream analyses. Enhancing resolution to this degree requires specialized equipment, intricate protocols, and substantial computational infrastructure, thereby limiting its availability to highly specialized institutions and facilities. On the other hand, the usual single-cell RNA sequencing (scRNA-seq) is capable of extracting gene expression at a localized single-cell level but lacks spatial contexts. Given observations of the complementary nature of the pro-and-cons of these two existing mainstream gene

- · W. Zhou is with the Heinz College of Information Systems and Public Policy and Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA, 15213.
- · J. Du is with the Department of Statistics and Data Science and Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA, 15213.

expression extraction technologies, the following question can be asked:

Does there exist an alternative solution to reconstruct or impute spatial context information from gene expression data?

An approach towards solving this question can potentially avoid the need for the deployment of spatial transcriptomic, making the tasks we intend to carry out with spatial context to be more accessible.

As an initial step toward tackling this question, we propose a generic representation learning and transfer learning framework enabling the inference of spatial context from purely gene expression data ( e.g. , collected from scRNAseq) by drawing reference from datasets where gene expression is paired with spatial information ( e.g. , collected from spatial transcriptomics). This procedure requires little or even no reliance on ST technologies depending on the type of reference data that the user intends to use, which can even come from existing publicly available databases. Specifically, our representation learning framework is based on the variational autoencoder (VAE) models [23], consisting of an encoder and a decoder network and trained on the reference datasets. Instead of using the standard ELBO loss, we propose to regularize it with a distance-preserving loss function, which utilizes the spatial context as auxiliary information to enforce the learned representation to be geometrically similar to the reference dataset. This regularization helps the model capture useful marker genes and signals from gene expressions that are representative of their spatial contexts. During inference, since spatial context is not used as the input to the model, the model can be applied to gene expression datasets without spatial context information but may share similar expression and geometric patterns with the reference dataset, thereby effectively extrapolating to

Fig. 1: The model architecture of dp-VAE with vanilla VAE as its backbone. During the training stage, the distancepreserving loss ( L DP) enforces distance preservation between the spatial domain ( S ) and the latent space ( Z ), enforcing the encoder network ( e.g. , parameterized as a fully connected neural network) to capture meaningful signals of the tissue spatial characteristics of the provided gene expression ( X ). During the inference/test stage, the extracted latent representations are fed to an optimizer, such as multidimensional scaling (MDS) if it is a spatial reconstruction task instead of a spatial context imputation task.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-24d6e73f-b2ce-4a1e-8197-897297305861.jpg)

the task of spatial reconstruction and imputation of the missing spatial context in these inference datasets. Next, the spatial reconstruction task can be formulated as an optimization that reversely solves for the optimal coordinates sets aligning with the distance geometry of extracted latent representations from the input gene expression data. This optimization can be bypassed when we set the latent representation space to match the dimensionality of the tissue's spatial domain. An illustration of the pipeline is provided in Figure 1.

Contributions. The contributions of this paper are summarized as follows:

- 1) We propose dp-VAE , a novel variational-inferencebased representation learning framework, capable of extracting useful signals of the spatial context from gene expression data. Its uniqueness lies in leveraging spatial context information exclusively during the training phase, thereby enabling spatial reconstruction and imputation of gene expression datasets, which is a task underexplored by existing studies.
- 2) We conduct an extensive empirical evaluation of the proposed method's robustness over 27 publically available datasets on our method's training stability, out-ofsample performance and evaluation, and transfer learning (out-of-distribution) inference capacity, providing a comprehensive analysis of the effectiveness of our method.
- 3) Additionally, we establish the theoretical connection between the distance-preserving loss and the distortion or bi-Lipschitz condition, contributing to the scarce literature on characterizing the mathematical notion of smoothness of distance-preserving mappings in the

context of generative models.

Our approach enables researchers to bridge the gap between standard gene expression profiles and spatial analyses, paving the way for more scalable and cost-effective research that integrates spatial insights into a variety of downstream applications.

## 2 RELATED WORKS

Spatial transcriptomics analysis has consistently garnered interest from the statistical machine learning community. It includes different tasks such as cellular deconvolution, clustering, visualization, cell classification, batch correction, gene imputation, and differential expression [35, 31, 33, 13]. Classical statistical models such as marked point process [14], Gaussian process [46], generalized linear spatial models [45, 61], and autoregression [35] have been deployed for these modeling tasks, where spatial context information are typically integrated into these methods in the form of covariance matrices. In the meantime, machine learning models empowered with deep neural networks can instead directly leverage the spatial context information. For example, convolution neural network-based methods such as CoSTA [54], and/or graph-based methods using graph neural networks such as GLISS [62], SpaGCN [20], GraphST [31]. See [58] for a survey of these methods.

Our work aligns more closely with literature that explores how various generative models ( e.g. autoencoders) can be applied to learning efficient representation from ST and scRNA-seq, such as STAGATE [9], scVSC [50], scVI [33] and gimVI [32]. Some studies have also proposed spatial reconstruction methodologies of tissue spatial structure,

such as Tangram [5], and NovoSpaRc [38]. However, to the best of our knowledge, our work is the first to propose to bridge representation learning with tissue spatial reconstruction drawing reference from out-of-sample and outof-distribution datasets. Additionally, the previously mentioned statistical and machine learning models dedicated to spatial transcriptomics are also inapplicable to such tasks, as the spatial contexts are required as a conditional input to the models during the inference stage.

The distance-preserving loss used in our study has (concurrently) emerged in parallel domains within the broader machine-learning literature. For instance, distancepreserving mappings have been introduced for GANs for unsupervised domain mapping [2], image-to-audio transformation [40, 22]. It has also appeared as regularizers for VAEs [6] and adversarial learning [18, 37]. A stronger notation compared to distance-preserving, named isometry [56, 26, 4] has also been extensively studied in machine learning models such as VAEs. Building on top of these works, our derivation of distance-preserving loss' connection to distortion or the bi-Lipschitz condition contributes to formalizing this notion.

## 3 PROBLEM SETUP

At the high level, the goal of representation learning is to learn a mapping that extracts latent embeddings z ∈ R d from given gene expression x ∈ R D . In our setting, we additionally assume that we have access to some spatial context information s ∈ R 2 associated with each gene expression indicating where the cell is located in the tissue slice. It is only available in the training stage, but not when during inference. We expect that the learned representation learning model can also undergo some procedure to reconstruct the spatial domain given only gene expression during the inference stage. Formally, given the gene expression a set of gene expressions denoted as x i ∈ R D , i = 1 , . . . , N , the goal of spatial reconstruction is to obtain a mapping to infer the spatial coordinates s i ∈ R 2 for each cellular location i .

In some cases, this task can be framed as an imputation task, where the spatial information of some cell has already been observed, indexed by i ∈ I ⊆ { 1 , . . . , N } . Therefore, imputation aims to infer the rest of the spatial coordinates locations s i , where i ∈ { 1 , . . . , N } / I . We remark that, in practice, the spatial reconstruction tasks here are useful for many downstream analyses, such as spatial domain identification, region-region interaction analysis, and spatial trajectory analysis. The inference task is similar in spirit to the generic task of transfer learning [39].

## 4 METHOD

To tackle the aforementioned challenges, we propose a gene expression representation learning framework. The model, named distance-preserving variation autoencoder ( dp-VAE ), incorporates the spatial information into the learning objective as a penalty for the error of the pairwise similarity between the model's learned latent variable with the spatial information of the training data. In this case, the spatial information can also be viewed as an auxiliary variable during the training stage rather than a direct input of the model, therefore it is not required during the inference stage . During inference, spatial reconstruction, and spatial 1 imputation, this penalty can be used reversely to decipher the spatial coordinates from the extracted latent representations of the data.

In this section, we give a thorough description of our method by first introducing the background knowledge of variational autoencoders in Section 4.1, then we describe the proposed dp-VAE in Section 4.2. Next, we formulate the spatial reconstruction and imputation tasks as a tractable optimization problem. Finally, we comment on the theoretical connections between this distance-preserving loss and the distortion and bi-Lipschitz condition in Section 4.4.

## 4.1 Background: Variational Autoencoder

Variational autoencoders (VAEs) are representation-learning models that have been proven to be capable of effectively extracting representations from gene expression data. The architecture of a traditional VAE includes an encoder-decoder model designated by the posterior probability p θ ( z x | ) and the likelihood q ϕ ( x z | ) , where x is the gene expression, and z is the latent variable that lives in low-dimensional space R m ′ ( m ′ ≤ m ). The objective of VAE is to estimate the parameterized marginal distribution p ϕ ( x ) of the gene expression

$$\Phi _ { \phi } ( y ) = \int p _ { \phi } ( y | z ) p ( z ) d z,$$

where p z ( ) serves as a predefined prior distribution for the latent variable. Given the computational inefficiencies linked to the normalizing constant in likelihood calculations, VAEs prioritize maximizing the evidence lower bound (ELBO) obtained via variational inference as opposed to directly optimizing (1):

$$\underbrace { \mathbb { E } _ { z \sim q _ { \theta } ( z | x ) } \log p _ { \phi } ( x | z ) } _ { - \mathcal { L } _ { r e c o n } } - \underbrace { D _ { K L } \left ( q _ { \theta } ( z | x ) | | p ( z ) \right ) } _ { \mathcal { L } _ { K L } }. \quad ( 2 )$$

This ELBO format is proven to be a lower bound to the marginal likelihood log p ϕ ( y ) ; the derivation is included in the appendix. We refer to the negative of the first term as the reconstruction loss 2 L recon , and the negative of the second term as the KL divergence L KL . Optimizing (2) involves tuning parameters θ and ϕ , essentially conducting representation learning. This process enables the encoder network p θ ( x ) to effectively translate signals from an ambient space R D into a fundamental subspace R d of lower dimensions.

However, the existing encoding paradigm presented in (2) reveals implicit biases, undermining its efficacy. Two primary concerns highlight this limitation: ( ) The likelihoodi focused optimization objective, while capturing some meaningful signals, underscores probabilistic aspects that could be incorrectly specified due to cross-dependencies or nonidentical distributed genomic samples. This could lead to entangled and biased representations of spatial context to be learned from the training data; ( ii ) In the meantime, directly

- 1. Whereas conditional models such as CVAE [43, 52] do require the spatial context as inputs, therefore cannot be applied in this setting.
- 2. To see why it is called the reconstruction term, if setting the likelihood and the posterior distribution to Gaussian distributions, the term is simplified to ∥ ̂ x i -x i ∥ 2 , which is the standard autoencoder reconstruction loss.

Fig. 2: Pearson correlation between the pairwise distance of gene expression ( X ) and the pairwise distance ( S ) of spatial location, computed on 27 spatial transcriptomics datasets. The red line indicates zero correlation, and the error bars are the 95% confidence intervals of the correlation coefficient calculated using Fisher's Z transformation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-8f744b91-87a8-4421-934d-8d97a9cd2f0f.jpg)

inserting the spatial information as part of the input variable violates our assumption about the setting where spatial context is not available as input during the inference stage, but rather should be the output of the model. Therefore, there is a need to develop a new method that incorporates likelihood-based generative modeling strategies with spatial information.

## 4.2 Distance-Preserving Variational Autoencoder

We propose to leverage the spatial information by incorporating it in the objective function as a distance-preserving regularization term of the VAE model. First, we define the distance-preserving loss

$$\mathcal { L } _ { \text{DP} } = \frac { 1 } { n ^ { 2 } } \sum _ { i, j } | \| z _ { i } - z _ { j } \| - \lambda \| s _ { i } - s _ { j } \| | \,, \quad \ ( 3 ) \quad \ s _ { \text{f} } ^ { \text{p} ^ { 1 } }$$

in attributes like gene expression and biological functions. Representation-based generative models, which aim to capture such spatial relationships and tissuelevel patterns, should utilize such geometric structures in learning disentangled representations that are meaningful. Prior works [20, 45] have illustrated that proximity in physical distance between cells has a strong association with their gene expressivity. This is also the initial of many spatial transcriptomics machine learning models to incorporate spatial context as inputs to enhance the modeling of cellular behavior and tissue organization, enabling more robust, context-aware insights into biological systems. If viewed reversely, it is also reasonable to claim that the learned gene expression representations can also reveal the associated spatial context from where they are located in the tissue.

where ∥·∥ is the Euclidean norm. The intuition behind the design of this loss is to ensure that the latent space preserves the relative pairwise distances between data points in the original space, while the scaling factor λ can be adjusted appropriately to account for differences in the variance of the underlying prior distributions of the latent and original spaces. This encourages the model to maintain geometric consistency while learning meaningful representations of gene expression aligned with the objective of VAE. Next, we incorporate this loss as a weighted regularize in the objective function of the training model

$$\alpha _ { 1 } \left ( \mathcal { L } _ { \text{recon} } + \beta \mathcal { L } _ { \text{KL} } \right ) + \alpha _ { 2 } \mathcal { L } _ { \text{DP} }, \quad \ \ ( 4 ) \quad \text{rever}$$

where the reconstruction term L recon and KL divergence term L KL are as defined in (2). The outer-level regularization coefficients α , α 1 2 ≤ 0 control for the numerical balance of the distance-preserving loss and the ELBO loss, and β is the temperature coefficient for the KL-divergence term [19]. The model is trained by learning the parameters of the encoder and decoder networks θ and ϕ as well as the scaling parameter λ to minimize (4).

The purpose of introducing the distance-preserving penalty term can be viewed from two perspectives:

- ( i ) From the biological perspective, neighboring cells, often from similar tissue types, tend to exhibit homogeneity
- ( ii ) From the modeling perspective, the penalization term directs the model to focus on capturing useful signals from the gene expression, which might serve as marker genes for its spatial information.

It is noted that the first point is also supported by a simple experiment we conduct, illustrated in Figure 2. It can be seen that out of the 27 datasets that we are using (see the experiment section for a detailed description of them), 23 datasets show a statistically significant ( 95% ) positive correlation between the spatial distance of the sample locations and the gene expression similarity. This highlights the potential of seeing the spatial context as a useful source of information to aid the learning of gene expression representations, or reversely, gene expression representations are capable of revealing information about their associated spatial context.

In practice, we recommend choosing a balanced penalization coefficient α 1 and α 2 to strike a good balance between the statistical accuracy (controlled by the reconstruction and KL term) and the spatial signal preservation (distance preserving term), achieving strong interpolation and generalization across different datasets. Also, we recommend a masked version of the distance-preserving:

$$\widetilde { \mathcal { L } } _ { \text{DP} } = \frac { 1 } { n ^ { 2 } } \sum _ { i, j } \text{G} _ { i, j } \left | \left \| z _ { i } - z _ { j } \right \| - \lambda \right \| s _ { i } - s _ { j } \right \| \right |,$$

where G is a { 0 1 , } n × n mask matrix that denotes the spatial connectivity of the cells, prespecified by the user. This mask

Fig. 3: An illustrative example of the masking procedure is shown. Suppose we have three spatial locations, so the mask matrix G is 3 × 3 . One user-specified connection in S is masked during training by setting the corresponding entries in G . Consequently, the distance-preserving regularization in L DP is not enforced between those two locations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-e7195963-ee3e-4f78-98f0-d07f60a88c27.jpg)

allows the incorporation of biological knowledge to selectively emphasize spatial relationships that are most relevant for model training. For example, if one has prior knowledge about the association between cells indexed by i and j , then G i,j can be set to one. The matrix G can also be specified by algorithms such as gene expression clustering, allowing it to generalize to different tasks where spatial information is present for model training. An illustration of the masking procedure is provided in Figure 3.

## 4.3 Spatial Reconstruction and Imputation

In the next step, we aim to carry out the spatial reconstruction and imputation task, where we only have inference data consisting of gene expression but not spatial context information. Assuming that we have access to some gene expression data that are paired with spatial context information and are statistically similar to the inference data, then dp-VAE can be trained on reference dataset, but transfered to the inference data for application. Specifically, since dp-VAE encodes the spatial context patterns in the latent representation via a distance-preserving regulated mapping, we can reversely find the optimal spatial coordinates in the spatial domain that align with the distance patterns in the extracted latent space. This can be formulated as an optimization problem regarding the estimated spatial coordinates. Additionally, if our task is imputation, meaning that some spatial coordinates may be exactly available, the solution to the estimation is, therefore, an optimization problem with constraints that these exact coordinates must match.

Formally, the procedure is as follows: First, we obtain the latent embeddings z i for each gene expression in the inference dataset indexed i = 1 , . . . , N by passing them through the encoder network p θ of the fitted dp-VAE . Then, solve the following optimization problem to obtain the estimated spatial information ̂ s i ,

$$\min _ { \widehat { s } _ { i } } \ \sum _ { i, j } | \| z _ { i } - z _ { j } \| - \| \widehat { s } _ { i } - \widehat { s } _ { j } \| |$$

$$s. t. \ \widehat { s } _ { i } = s _ { i }, \forall i \in \mathcal { I },$$

$$\min _ { \widehat { s } _ { i } } & \sum _ { i, j } | \| z _ { i } - z _ { j } \| - \| \widehat { s } _ { i } - \widehat { s } _ { j } \| | \quad \quad ( 6 ) \quad \text{ of t} \\ s. t. & \widehat { s } _ { i } = s _ { i }, \forall i \in \mathcal { I }, \quad \quad ( 7 ) \quad \text{ uses} \\ & \text{char}$$

where I denotes the index set of the entries whose exact spatial coordinate is known in the imputation task. Intuitively, the objective (6) enforces the pairwise distance of the estimated spatial information to be consistent with the latent embeddings, and the constraint (7)

requires that the estimated spatial information be the same for all known spatial coordinates. The problem can be relaxed to an unconstrained optimization problem by modifying the objective function to a penalized version: ∑ i,j |∥ z i -z j ∥ - ∥ ̂ s i -̂ s j ∥| + γ ∑ i ∈I | ̂ s i -s i | for easier solving, where γ denotes the regularization coefficient ( i.e. Lagrange multiplier).

The objective function (6) is highly similar to the distance-preserving loss defined in (3). Observe (3) encodes distance-preserving property into the learned latent embedding, while (6) extracts the embedded distancepreserving property to reconstruct spatial information. This can be viewed as an implicit second set of encode-decodestructure, which we specifically utilize, aside from the first set of encoder-decoder-structure in VAE. However, a subtle difference between the two losses is that the parameters to be optimized in (6) and (3) are fundamentally different, where (6) is parameterized by the spatial information s i , and (3) is parametrized by the parameters in the encoder network θ .

Fig. 4: Illustration of forward and inverse mappings. The mapping that recovers the coordinates S from pairwise distances D S is not unique or non-deterministic when |I| &lt; 3 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-f0c29f58-0a7c-4c05-b99e-e6c6c092e7bd.jpg)

The cardinality of the observed spatial information index set I plays a critical role in the solution of the optimization problem. The optimal solution to (6) is only unique when |I| ≥ 3 but otherwise not, as outlined in Figure 4. When |I| = 0 , the problem boils down to the classical multidimensional scaling (MDS) [41] objective, which is a well-studied dimensionality reduction problem and can be solved via, e.g. , the SMACOF algorithm [27].

Additionally, when |I| = 0 , the optimization procedure can also be skipped: By directly setting the dimension of the latent space to d = 2 , the learned latent embedding can then be naturally regarded as the estimated spatial information, i.e. z i = ̂ s i . However, when the extracted latent representation is intended for other uses, directly setting d = 2 is not recommended. This is because the dimension of the latent embedding is physically interpreted as the intrinsic dimension of the manifold where the gene expression features lie, which is typically greater than two. Deliberately setting it to two may result in over-compression of the gene expression feature, resulting in information loss and making the learned embeddings ineffective for other uses. Tuning for the latent dimension could also require a changed set of training hyperparameters, which should be cautious.

The whole algorithm is summarized in Algorithm 1. In the ideal case where the algorithm is capable of fully capturing the distribution, the recover the real spatial information matrix invariant under translation, scaling, and rotation .

## Algorithm 1 Spatial Reconstruction with dp-VAE

Require: Training datasets X ,S tr tr , test dataset X te .

- 1: Train dp-VAE model via maximizing the distancepreserving ELBO (4).
- 2: Compute latent encoding ̂ Z te from the trained encoder using X te , and obtain its pairwise distance matrix ̂ D te .
- 3: Solve optimization problem in (6)-(7) to obtain coordinate matrix ̂ S te .
- 4: return Coordinate matrix ̂ S te .

This invariance is sufficient in practice because the relative spatial relationships between cells are preserved, which is the key information needed for most spatial transcriptomics analyses. The absolute coordinates are often arbitrary and less important than the relative positions and distances between cells in the tissue context. We will describe how such similarity can be evaluated in the experiment section.

## 4.4 Theoretical Property

We briefly comment on our finding on the theoretical connection between the proposed distance-preserving loss and the distortion [8] or bi-Lipshitz condition[36, 47]. Specifically, the enforcement result can be viewed as a probably approximately correct (PAC) learning version of the two, summarized in the following theorem.

Theorem 1. For any ϵ &gt; 0 , there exist some constants λ &gt; 0 and

$$L \lesssim \mathcal { O } \left ( \frac { \mathbb { E } \left [ \mathcal { L } _ { \text{DP} } \right ] } { \lambda \epsilon } \right ). \quad \quad \quad \ ( 8 ) \text{ \ \ will } \\ \exp \epsilon$$

such that the following condition holds with probability greater than 1 -ϵ :

$$\lambda \| s - s ^ { \prime } \| \leq \| z - z ^ { \prime } \| \leq L \cdot \lambda \| s - s ^ { \prime } \|, \quad \ ( 9 ) \, \overset { \text{ NVIDIA} } { \text{Fav.} }$$

where the probability arises from the following generation process

$$( x, s ), ( x ^ { \prime }, s ^ { \prime } ) \stackrel { i i d } { \sim } p, \, z \sim f _ { \theta } ( y ), \, z ^ { \prime } \sim f _ { \theta } ( y ^ { \prime } ). \quad ( 1 0 ) \quad \text{and} \quad \text{for two}$$

In Theorem 1, expression (15) resembles the distortion or bi-Lipschitz condition defined in prior works, except that it only holds probabilistically. This probability is incurred by the underlying data generation process and the learned VAE generative model. The distortion constant or bi-Lipschitz constant L is porportional to the expected value of the distance-preserving loss L DP , which in the ideal case ( i.e. perfect generalization) scales with the value of its training error. This means that the encoder network becomes 'smoother' as training progresses, and the distancepreserving loss is what regulates such behavior of this generative model. We note that this result summarizes and can be potentially generalized to characterize the theoretical properties of the generative networks regulated with distance-preserving loss, as it is not unique to the dp-VAE framework we are using in this paper. We provide the exact bound to Theorem 1 in the appendix.

## 5 EXPERIMENTS

In this section, we conduct experiments to validate the effectiveness of our method and shed light on some of our method's practical implications, including training stability and in-sample behavior of dp-VAE (Section 5.2), its comparison against vanilla VAE in out-of-sample evaluation (Section 5.3) and out-of-distribution inference applications, i.e., transfer learning (Section 5.4). We use a total of 27 processed Visium Spatial Gene Expression data from the 10x Genomics database provided by Scanpy [51]. We only select the top 100 most variable genes for all datasets to ensure model stability and learning feasibility 3 The names of the datasets are the initials of the datasets, whose full names and correspondence can be found in the appendix. Before we proceed, we first describe our experiment setups in the next subsection.

## 5.1 Experiment Setups

Model description. Since the VAE model is a special case of dp-VAE when the distance-preserving regularization coefficient is set to zero. Therefore, we will only describe our set-up for the dp-VAE model. The dp-VAE model consists of an encoder and a decoder network, each comprising a one-layer fully connected neural network with 124 hidden layers and a latent dimension of 2 ; therefore, the procedure of solving for an optimization problem can be bypassed by directly using the latent representations as the estimated spatial coordinates. The prior, posterior, and likelihood distributions are all specified as Gaussian distributions. We apply no mask between any pair of entries. The regularization coefficients α , β 1 for the ELBO objective are set to one throughout all experiments, therefore for simplicity, we will only refer to α 2 as the regularization coefficient in our experiment section. The learning rate is set to 5 × 10 -4 or 10 -4 depending on the experiment, and the total number of training epochs is set to 2000 . All experiments are run on NVIDIA GeForce RTX 4070Ti GPU.

Evaluation protocals. We evaluate our model based on the Procrustes distance [15, 16, 25] between the reconstructed and the ground-truth spatial information matrix. Formally, for two datasets denoted by matrices D 1 , D 2 ∈ R n × m , the procedure of calculating the Proscrutes distance between D 1 and D 2 can be defined as the procedure of finding the optimal translation, rotation, and scaling transformation from R m to R m , such that the Frobenius norm between the two matrices is minimized. Formally, it can be described as

$$\min _ { \mathbf A, b, c } \ \| c \cdot \mathbf A \mathbf D _ { 1 } + b - \mathbf D _ { 2 } \| _ { F }.$$

Here A ∈ R m,m is an orthogonal matrix, b ∈ R m , and c is a scalar. We will refer to it as the spatial reconstruction error in the following sections when clear from context.

## 5.2 In-Sample Training Behaviors

First, we evaluate the training stability of dp-VAE under various choices of the regularization coefficient. This analysis is crucial because the proposed modification of the original ELBO objective may introduce uncertainty about whether the new objective ensures convergence or risks numerical instability.

- 3. Most datasets only contain roughly around three thousand entries, making the problem high-dimensional and numerically unstable.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-cae06f3b-0df6-43ef-bd5b-4b0102d7bbbb.jpg)

Fig. 5: Latent space visualization over different training epochs. Each dot represents the 2-dimensional embedding of a cell. The blue and orange stars represent two randomly selected triplets of spatial locations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-39fe7569-f00f-49ab-8f93-98d400891e8f.jpg)

Fig. 6: Out of sample (leave-three-out) spatial location imputation. Stars of different colors represent a randomly selected fixed triplet. The top row is the ground truth spatial context, and the bottom row is the reconstructed spatial context.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-86b5bfc5-6a99-44fd-9548-65c2cd2ed646.jpg)

Fig. 7: Training stability for different choices of regularization coefficients for the distance-preserving loss over 1 × 10 3 epochs of training, averaged over 27 datasets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-ac68aa2d-031b-4a24-83a9-2b691f642bbb.jpg)

Figure 7 shows the normalized training loss curve of the dp-VAE model with different choices of the regularization coefficient α 2 averaged over 27 datasets used as training data. Our results show that increasing α 2 may introduce small abrupt oscillations in the loss curve, globally the loss manages to converge at the same point for all choices of α 2 . This means that different choices of α 2 do not significantly impact the convergence speed and stability of the training, making its training behavior almost identical to the standard VAE model. In practice, dp-VAE can be tuned and trained very similarly to a standard VAE. However, it is noted that using large learning rates may amplify the oscillation in the

Fig. 8: Heatmap of spatial reconstruction error of four Visium mouse brain datasets obtained from dp-VAE , where the training data is selected to be only VMBSP.

training, making it numerically unstable.

Figure 5 reports the reconstructed spatial coordinates at different stages of training (every 400 epochs) of dp-VAE compared with VAE on the Mouse Brain Sagittal Anterior Section 2 dataset. We randomly pick two sets of three anchor points to help visualize how the distance-preserving property is gradually enforced. It can be seen that the dp-VAE

Fig. 9: Comparions of out-of-sample spatial reconstruction error between VAE and dp-VAE , evaluated on 27 datasets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-360e94e8-4b3e-4259-a764-bbac408926ba.jpg)

Fig. 10: Sensitivity analysis of spatial reconstruction error with different regularization coefficient choices, evaluated over 27 datasets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-5c9e5658-9f00-4d1b-95c5-609d6587b6f0.jpg)

model quickly converges within 800 -th to accurately capture the global shape of the spatial context in the training data. At more refined resolutions, the relative locations of the anchor points have been accurately recovered by dp-VAE , where the reconstructed triangular shapes are almost identical to the ground-truth data. On the other hand, the reconstructed spatial context by the vanilla VAE model does not converge correctly to the training data. This implies that the distancepreserving loss in the dp-VAE model enforces disentangled

Fig. 11: The matrix of out-of-distribution spatial reconstruction error incurred by dp-VAE . The model uses each column data as the training data to infer four test datasets labeled on the rows.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-46f3a4b0-5237-422f-b666-662e5c16e23e.jpg)

latent representations from gene expressions, and captures meaningful biological signals helpful for the task of spatial reconstruction.

## 5.3 Out-of-Sample Evaluation

To evaluate the out-of-sample generalization ability and verify that the dp-VAE does not simply memorize the training data, we carry out comparative studies comparing VAE and dp-VAE , as well as different choices of the regularization coefficient α 2 on out-of-sample spatial reconstruction evaluation tasks. We obtain our training and test data by dividing each one of 27 datasets by a ratio of 80% -20% , yielding 27 training datasets paired with 27 test datasets.

Figure 9 reports the normalized spatial reconstruction error comparing the VAE and dp-VAE models on the 27 test datasets. It can be seen that the dp-VAE significantly outperforms vanilla VAE on spatial reconstruction tasks across almost all datasets, demonstrating the generalization capability of the distance-preserving regularizer. This result aligns with our previous theoretical analysis. Furthermore, Figure 2 indicates a statistical dependence between gene expression and its spatial context, suggesting that spatial context can be inferred from signals derived from gene expressions, even under mild extrapolation.

Figure 6 follows Figure 9 and reports the detailed reconstructed spatial context (bottom row) with the real spatial context (top row). The three colored spots are random anchor points chosen within the test dataset. It can be seen that the three imputed anchor points are placed in the correct relative position within the reconstructed spatial context. In the meantime, the relative locations and distances of the three points are also accurately represented. This confirms the nuanced rationale underlying the strong performance on the aggregated metric of spatial reconstruction error.

Figure 10 presents the spatial reconstruction error averaged over 27 test datasets, evaluated under varying regularization coefficient values. As the regularization coefficient increases, the test (out-of-sample) spatial reconstruction error decreases logarithmically. This finding further validates that applying stronger distance-preserving regularization in dp-VAE enhances its ability for spatial reconstruction extrapolation.

However, we note that using a larger regularization coefficient is not always optimal in practice. With larger networks, excessively strong regularization can increase the risk of overfitting the training data, ultimately impairing extrapolation performance. Furthermore, as previously mentioned, a high regularization coefficient can lead to instability during training, potentially causing the learning process to fail. Therefore, selecting an appropriate regularization coefficient requires careful consideration of various practical factors.

## 5.4 Out-of-Distribution Evaluation and Inference

In previous sections, we have been evaluating the performance of our methods, where the training and test data are split from the same dataset. We further apply our method to a more challenging scenario, where the train and test data come from separate datasets, and biological structure may

be potentially different in the two populations. We experiment on four mouse brain tissue spatial transcriptomics datasets, creating 4 × 4 = 16 pairs of training and test datasets. We aim to apply the proposed dp-VAE to extract insights about the transfer learning capability between these four datasets.

Figure 11 visualize the normalized spatial reconstruction error matrix, where the columns are the test data where the error is calculated, and the rows are the training data where the dp-VAE model is trained. It can be seen that the train-test pairs VMBSPS2-VMBSP, VMBSPS2-VMBSP, and VMBSA-VMBSAS2 incur relatively small error scales compared to other train-test pairs. This may indicate that the joint distribution of gene expression and spatial context between these sets of tissues are similar, and this shared information may be used to infer one another.

Figure 8 shows the detailed view of the spatial reconstruction error at the coordinate level of four test datasets using VMBSP as the training data, corresponding to the bottom row of Figure 11. The results reveal distinct patterns of spatial reconstruction error between different tissue regions, suggesting that heterogeneous gene expression carries varying levels of biological signals closely linked to its spatial context. Further research is anticipated to illuminate the underlying causes of this heterogeneity.

## 6 CONCLUSION AND DISCUSSION

In this paper, we proposed a novel representation learning framework based on a variational autoencoder structure to enable tissue spatial reconstruction and imputation from gene expression data. Note that our approach relies minimally on spatial transcriptomics in the sense that they are only used during the training stage as an auxiliary variable, which allows for diverse downstream applications that rely on having spatial context to be carried out when such information is unavailable during the inference stage, opening up the possibility of a broad spectrum of genomics-related studies that might be once unavailable on many datasets. The technical novelty lies in the distancepreserving loss that we introduce to regulate the training process, which helps the model to effectively capture gene expression signals oriented towards describing its spatial context. Spatial reconstruction can be achieved by solving an optimization problem or automatically achieved when the latent space dimension is set to be the dimension of the spatial context. We also point out the theoretical gap between the distance-preserving loss or bi-Lipschitz condition in the context of generative modeling. Finally, we study the training robustness, out-of-sample inference, and transfer learning performance of our proposed method on 27 spatial transcriptomics datasets to evaluate and demonstrate its empirical effectiveness.

A potential future research direction is to more extensively evaluate the method in more challenging transfer learning settings, such as examining the performance of the method in inferring the spatial context of gene expressions that are at a larger scale or with higher spatial resolution, studying the effect of different masking procedures [6], performance in different downstream analysis when integrated with other algorithms, or examine the performance of the method when using reference datasets with higher dissimilarity to the inference datasets. This comparative study might need to incorporate more biological domain expert insights and is more resource-demanding. Nevertheless, we believe that this general framework has strong potential to be modified to various tailored settings.

## REFERENCES

- [1] Deep generative models for spatial networks , 2021.
- [2] Sagie Benaim and Lior Wolf. One-sided unsupervised domain mapping. Advances in neural information processing systems , 30, 2017.
- [3] Tristan Bepler, Ellen Zhong, Kotaro Kelley, Edward Brignole, and Bonnie Berger. Explicitly disentangling image content from translation and rotation with spatial-vae. Advances in Neural Information Processing Systems , 32, 2019.
- [4] Kosio Beshkov, Jonas Verhellen, and Mikkel Elle Lepperød. Isometric representations in neural networks improve robustness. arXiv preprint arXiv:2211.01236 , 2022.
- [5] Tommaso Biancalani, Gabriele Scalia, Lorenzo Buffoni, Raghav Avasthi, Ziqing Lu, Aman Sanger, Neriman Tokcan, Charles R Vanderburg, Asa Segerstolpe, Meng ˚ Zhang, et al. Deep learning and alignment of spatially resolved single-cell transcriptomes with tangram. Nature methods , 18(11):1352-1362, 2021.
- [6] Nutan Chen, Patrick van der Smagt, and Botond Cseke. Local distance preserving auto-encoders using continuous knn graphs. In Topological, Algebraic and Geometric Learning Workshops 2022 , pages 55-66. PMLR, 2022.
- [7] Tsai-Ying Chen, Li You, Jose Angelito U Hardillo, and Miao-Ping Chien. Spatial transcriptomic technologies. Cells , 12(16):2042, 2023.
- [8] Leena Chennuru Vankadara and Ulrike von Luxburg. Measures of distortion for machine learning. Advances in Neural Information Processing Systems , 31, 2018.
- [9] Kangning Dong and Shihua Zhang. Deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder. Nature communications , 13(1):1739, 2022.
- [10] Ruben Dries, Jiaji Chen, Natalie Del Rossi, Mohammed Muzamil Khan, Adriana Sistig, and GuoCheng Yuan. Advances in spatial transcriptomic data analysis. Genome research , 31(10):1706-1718, 2021.
- [11] Ruben Dries, Qian Zhu, Rui Dong, Chee-Huat Linus Eng, Huipeng Li, Kan Liu, Yuntian Fu, Tianxiao Zhao, Arpan Sarkar, Feng Bao, et al. Giotto: a toolbox for integrative analysis and visualization of spatial expression data. Genome biology , 22:1-31, 2021.
- [12] Jin-Hong Du, Zhanrui Cai, and Kathryn Roeder. Robust probabilistic modeling for single-cell multimodal mosaic integration and imputation via scvaeit. Proceedings of the National Academy of Sciences , 119(49):e2214414119, 2022.
- [13] Jin-Hong Du, Tianyu Chen, Ming Gao, and Jingshu Wang. Joint trajectory inference for singlecell genomics using deep learning with a mixture prior. Proceedings of the National Academy of Sciences , 121(37):e2316256121, 2024.

- [14] Daniel Edsg¨ ard, Per Johnsson, and Rickard Sandberg. Identification of spatial expression trends in single-cell gene expression data. Nature methods , 15(5):339-342, 2018.
- [15] Colin Goodall. Procrustes methods in the statistical analysis of shape. Journal of the Royal Statistical Society: Series B (Methodological) , 53(2):285-321, 1991.
- [16] John C Gower. Generalized procrustes analysis. Psychometrika , 40:33-51, 1975.
- [17] Tiantian Guo, Zhiyuan Yuan, Yan Pan, Jiakang Wang, Fengling Chen, Michael Q Zhang, and Xiangyu Li. Spiral: integrating and aligning spatially resolved transcriptomics data across different experiments, conditions, and technologies. Genome Biology , 24(1):241, 2023.
- [18] Raia Hadsell, Sumit Chopra, and Yann LeCun. Dimensionality reduction by learning an invariant mapping. In 2006 IEEE computer society conference on computer vision and pattern recognition (CVPR'06) , volume 2, pages 1735-1742. IEEE, 2006.
- [19] Irina Higgins, Loic Matthey, Arka Pal, Christopher P Burgess, Xavier Glorot, Matthew M Botvinick, Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with a constrained variational framework. ICLR (Poster) , 3, 2017.
- [20] Jian Hu, Xiangjie Li, Kyle Coleman, Amelia Schroeder, Nan Ma, David J Irwin, Edward B Lee, Russell T Shinohara, and Mingyao Li. Spagcn: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network. Nature methods , 18(11):13421351, 2021.
- [21] Yaofeng Hu, Kai Xiao, Hengyu Yang, Xiaoping Liu, Chuanchao Zhang, and Qianqian Shi. Spatially contrastive variational autoencoder for deciphering tissue heterogeneity from spatially resolved transcriptomics. Briefings in Bioinformatics , 25(2):bbae016, 2024.
- [22] Chelhwon Kim, Andrew Port, and Mitesh Patel. Faceto-music translation using a distance-preserving generative adversarial network with an auxiliary discriminator. arXiv preprint arXiv:2006.13469 , 2020.
- [23] Diederik P Kingma and Max Welling. Auto-encoding variational bayes, 2022.
- [24] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. Advances in neural information processing systems , 29, 2016.
- [25] Wojtek Krzanowski. Principles of multivariate analysis , volume 23. OUP Oxford, 2000.
- [26] Yonghyeon Lee, Hyeokjun Kwon, and Frank Park. Neighborhood reconstructing autoencoders. Advances in Neural Information Processing Systems , 34:536-546, 2021.
- [27] DE LEEUW. Application of convex analysis to multidimensional scaling. Recent developments in statistics , pages 133-145, 1977.
- [28] Hongfei Li, Catherine A Calder, and Noel Cressie. Beyond moran's i: testing for spatial dependence based on the spatial autoregressive model. Geographical analysis , 39(4):357-375, 2007.
- [29] Shang Li, Kuo Gai, Kangning Dong, Yiyang Zhang,
- and Shihua Zhang. High-density generation of spatial transcriptomics with stage. Nucleic Acids Research , page gkae294, 2024.
- [30] Teng Liu, Zhao-Yu Fang, Zongbo Zhang, Yongxiang Yu, Min Li, and Ming-Zhu Yin. A comprehensive overview of graph neural network-based approaches to clustering for spatial transcriptomics t. liu et al. overview of spatial transcriptomics' spatial clutering. Computational and Structural Biotechnology Journal , 2023.
- [31] Yahui Long, Kok Siong Ang, Mengwei Li, Kian Long Kelvin Chong, Raman Sethi, Chengwei Zhong, Hang Xu, Zhiwei Ong, Karishma Sachaphibulkij, Ao Chen, et al. Spatially informed clustering, integration, and deconvolution of spatial transcriptomics with graphst. Nature Communications , 14(1):1155, 2023.
- [32] Romain Lopez, Achille Nazaret, Maxime Langevin, Jules Samaran, Jeffrey Regier, Michael I Jordan, and Nir Yosef. A joint model of unpaired data from scrnaseq and spatial transcriptomics for imputing missing gene expression measurements. arXiv preprint arXiv:1905.02269 , 2019.
- [33] Romain Lopez, Jeffrey Regier, Michael B Cole, Michael I Jordan, and Nir Yosef. Deep generative modeling for single-cell transcriptomics. Nature methods , 15(12):10531058, 2018.
- [34] Calvin Luo. Understanding diffusion models: A unified perspective, 2022.
- [35] Ying Ma and Xiang Zhou. Spatially informed celltype deconvolution for spatial transcriptomics. Nature biotechnology , 40(9):1349-1359, 2022.
- [36] Sepideh Mahabadi, Konstantin Makarychev, Yury Makarychev, and Ilya Razenshteyn. Nonlinear dimension reduction via outer bi-lipschitz extensions. In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing , pages 1088-1101, 2018.
- [37] Chengzhi Mao, Ziyuan Zhong, Junfeng Yang, Carl Vondrick, and Baishakhi Ray. Metric learning for adversarial robustness. Advances in neural information processing systems , 32, 2019.
- [38] Noa Moriel, Enes Senel, Nir Friedman, Nikolaus Rajewsky, Nikos Karaiskos, and Mor Nitzan. Novosparc: flexible spatial reconstruction of single-cell gene expression with optimal transport. Nature protocols , 16(9):4177-4200, 2021.
- [39] Sinno Jialin Pan and Qiang Yang. A survey on transfer learning. IEEE Transactions on knowledge and data engineering , 22(10):1345-1359, 2009.
- [40] Andrew Port, Chelhwon Kim, and Mitesh Patel. Earballs: Neural transmodal translation. arXiv preprint arXiv:2005.13291 , 1, 2020.
- [41] Nasir Saeed, Haewoon Nam, Mian Imtiaz Ul Haq, and Dost Bhatti Muhammad Saqib. A survey on multidimensional scaling. ACM Computing Surveys (CSUR) , 51(3):1-25, 2018.
- [42] Elizaveta Semenova, Yidan Xu, Adam Howes, Theo Rashid, Samir Bhatt, Swapnil Mishra, and Seth Flaxman. Priorvae: encoding spatial priors with variational autoencoders for small-area estimation. Journal of the Royal Society Interface , 19(191):20220094, 2022.
- [43] Kihyuk Sohn, Honglak Lee, and Xinchen Yan. Learning structured output representation using deep condi-

tional generative models. Advances in neural information processing systems , 28, 2015.

- [44] Casper Kaae Sønderby, Tapani Raiko, Lars Maaløe, Søren Kaae Sønderby, and Ole Winther. Ladder variational autoencoders. Advances in neural information processing systems , 29, 2016.
- [45] Shiquan Sun, Jiaqiang Zhu, and Xiang Zhou. Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies. Nature methods , 17(2):193-200, 2020.
- [46] Valentine Svensson, Sarah A Teichmann, and Oliver Stegle. Spatialde: identification of spatially variable genes. Nature methods , 15(5):343-346, 2018.
- [47] Alexandre Verine, Benjamin Negrevergne, Yann Chevaleyre, and Fabrice Rossi. On the expressivity of bilipschitz normalizing flows. In Asian Conference on Machine Learning , pages 1054-1069. PMLR, 2023.
- [48] Xiaomeng Wan, Jiashun Xiao, Sindy Sing Ting Tam, Mingxuan Cai, Ryohichi Sugimura, Yang Wang, Xiang Wan, Zhixiang Lin, Angela Ruohao Wu, and Can Yang. Integrating spatial and single-cell transcriptomics data using deep generative models with spatialscope. Nature Communications , 14(1):7848, 2023.
- [49] Zhengyang Wang, Hao Yuan, and Shuiwang Ji. Spatial variational auto-encoding via matrix-variate normal distributions. In Proceedings of the 2019 SIAM International Conference on Data Mining , pages 648-656. SIAM, 2019.
- [50] Zile Wang, Haiyun Wang, Jianping Zhao, Junfeng Xia, and Chunhou Zheng. scvsc: Deep variational subspace clustering for single-cell transcriptome data. IEEE/ACM Transactions on Computational Biology and Bioinformatics , 2024.
- [51] F Alexander Wolf, Philipp Angerer, and Fabian J Theis. Scanpy: large-scale single-cell gene expression data analysis. Genome biology , 19:1-5, 2018.
- [52] Shenghao Wu, Wenbin Zhou, Minshuo Chen, and Shixiang Zhu. Counterfactual generative models for timevarying treatments. arXiv preprint arXiv:2305.15742 , 2023.
- [53] Hang Xu, Huazhu Fu, Yahui Long, Kok Siong Ang, Raman Sethi, Kelvin Chong, Mengwei Li, Rom Uddamvathanak, Hong Kai Lee, Jingjing Ling, et al. Unsupervised spatially embedded deep representation of spatial transcriptomics. Genome Medicine , 16(1):12, 2024.
- [54] Yang Xu and Rachel Patton McCord. Costa: unsupervised convolutional neural network learning for spatial transcriptomics analysis. BMC bioinformatics , 22(1):397, 2021.
- [55] Yi Yang, Xingjie Shi, Wei Liu, Qiuzhong Zhou, Mai Chan Lau, Jeffrey Chun Tatt Lim, Lei Sun, Cedric Chuan Young Ng, Joe Yeong, and Jin Liu. Scmeb: spatial clustering with hidden markov random field using empirical bayes. Briefings in bioinformatics , 23(1):bbab466, 2022.
- [56] LEE Yonghyeon, Sangwoong Yoon, Minjun Son, and Frank C Park. Regularized autoencoders for isometric representation learning. In International Conference on Learning Representations , 2021.
- [57] Qichao Yu, Miaomiao Jiang, and Liang Wu. Spatial transcriptomics technology in cancer research. Frontiers
- in Oncology , 12:1019111, 2022.
- [58] Zexian Zeng, Yawei Li, Yiming Li, and Yuan Luo. Statistical and machine learning methods for spatially resolved transcriptomics data analysis. Genome biology , 23(1):83, 2022.
- [59] Xinyi Zhang, Xiao Wang, GV Shivashankar, and Caroline Uhler. Graph-based autoencoder integrates spatial transcriptomics with chromatin images and identifies joint biomarkers for alzheimer's disease. Nature Communications , 13(1):7480, 2022.
- [60] Xiang Zhou, Kangning Dong, and Shihua Zhang. Integrating spatial transcriptomics data across different conditions, technologies and developmental stages. Nature Computational Science , 3(10):894-906, 2023.
- [61] Jiaqiang Zhu, Shiquan Sun, and Xiang Zhou. Spark-x: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies. Genome biology , 22(1):184, 2021.
- [62] Junjie Zhu and Chiara Sabatti. Integrative spatial single-cell analysis with graph-based feature learning. Biorxiv , pages 2020-08, 2020.

## APPENDIX A OMITTED PROOFS IN SECTION 4.4

In this section. we provide the detailed proof for Theorem 1 in Section 4.4. We begin by describing a technical lemma that will be used in the main proof.

Lemma2. For any given distribution defined on some space X ⊆ R n , given arbitrary 0 &lt; δ &lt; 1 , there exist some 0 &lt; M 1 &lt; M &lt; 2 + ∞ , such that:

$$\mathbb { P } \left [ M _ { 1 } \leq \| x - x ^ { \prime } \| \leq M _ { 2 } \right ] \geq 1 - \delta, \quad \ \ ( 1 1 ) \quad \ \ \ \ O t _ { \ }.$$

where x, x ′ are drawn i.i.d. from this distribution.

Proof. If X is a finite space or bounded space, then simply take:

$$M _ { 1 } = \inf _ { x, x ^ { \prime } \in \mathcal { X } } \| x - x ^ { \prime } \|, \, M _ { 2 } = \sup _ { x, x ^ { \prime } \in \mathcal { X } } \| x - x ^ { \prime } \|.$$

If X is an infinite space or an unbounded space, such as X = R n . In this case, since we know that the Euclidean space R n can be exhausted by a nested sequence of compact balls centered at the origin:

$$\mathcal { X } \subseteq \mathbb { R } ^ { n } = \bigcup _ { r = 1 } ^ { + \infty } B ( 0, r ),$$

where B (0 , r ) := { x ∈ R n : ∥ x ∥ ≤ r } . Therefore, by the total probability property, we know that:

$$\mathbb { P } ( \mathcal { X } ) \leq \mathbb { P } \left [ \bigcup _ { r = 1 } ^ { + \infty } B ( 0, r ) \right ] = 1. \quad \quad ( 1 2 ) \quad \quad \text{The }!$$

By the dominated convergence theorem:

$$\mathbb { P } \left [ \lim _ { R \to + \infty } \bigcup _ { r = 1 } ^ { R } B ( 0, r ) \right ] = \lim _ { R \to + \infty } \mathbb { P } \left [ \bigcup _ { r = 1 } ^ { R } B ( 0, r ) \right ]. \quad ( 1 3 )$$

By definition, given any 0 &lt; δ &lt; 1 , there exists some R &gt; ′ 0 such that:

$$\mathbb { P } \left [ \bigcup _ { r = 1 } ^ { R ^ { \prime } } B ( 0, r ) \right ] \geq 1 - \delta,$$

which implies:

$$\mathbb { P } \left [ - 2 R ^ { \prime } \leq \| x - x ^ { \prime } \| \leq 2 R ^ { \prime } \right ] \geq 1 - \delta.$$

Therefore, we can set M 1 = -2 R ′ and M 2 = 2 R ′ , and this concludes the proof.

Theorem 3. Denote δ, M , M 1 2 be given by Lemma 2 with X set to S , let L DP be defined in (3) . Then, for any ϵ &gt; 0 , there exist some constants λ &gt; 0 and

$$L \leq \frac { M _ { 2 } } { M _ { 1 } } + \frac { \mathbb { E } \left [ \mathcal { L } _ { \text{DP} } \right ] } { \lambda M _ { 1 } ( \epsilon - \delta ) }. \quad \quad \ ( 1 4 ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

such that the following condition holds with probability greater than 1 -ϵ :

$$\lambda \| s - s ^ { \prime } \| \leq \| z - z ^ { \prime } \| \leq L \cdot \lambda \| s - s ^ { \prime } \|, \quad ( 1 5 ) \, \stackrel { \text{work} } { \text{variat} }$$

where the probability arises from the following generation process

$$( x, s ), ( x ^ { \prime }, s ^ { \prime } ) \stackrel { i i d } { \sim } p, \, z \sim f _ { \theta } ( y ), \, z ^ { \prime } \sim f _ { \theta } ( y ^ { \prime } ). \quad ( 1 6 ) \quad \text{graph}$$

Proof. Denote the following events in the probability space

$$\begin{array} { c c c c c c } \cdot & \cdot & \cdot & \cdot & \cdot \\ \mathcal { A } = \{ M _ { 1 } \leq \| s - s ^ { \prime } \| \leq M _ { 2 } \} \,, \\ \mathcal { B } = \{ \lambda \| s - s ^ { \prime } \| - C \leq \| z - z ^ { \prime } \| \leq \lambda \| s - s ^ { \prime } \| + C \} \\ \mathcal { C } = \{ \lambda \| s - s ^ { \prime } \| \leq \| z - z ^ { \prime } \| \leq L \cdot \lambda \| s - s ^ { \prime } \| \} \,, \end{array}$$

we will set L = M /M 2 1 + C/ λM ( 1 ) for now, which then by the definition of the three events, there is

$$\mathbb { P } ( \mathcal { B } | \mathcal { A } ) \leq \mathbb { P } ( \mathcal { C } | \mathcal { A } ).$$

On the other hand, using Markov's inequality and by the definition of the distance preserving loss (3), for any C &gt; 0 , we have:

$$\frac { \mathbb { E } \left [ \mathcal { L } _ { \text{DP} } \right ] } { C } \geq \mathbb { P } ( | \| z - z ^ { \prime } \| - \lambda \| s - s ^ { \prime } \| | \geq C ).$$

Rearranging terms on both sides, we get

$$\det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det\det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \dotdet \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det
\det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \cot \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \dot
\det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det \det
\end{cases}$$

On the other hand, using the law of total probability, we can bound the probability of event C by:

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \delta$$

The first line uses the law of total probability, the second line uses (17) and Lemma 2, and the third line uses the law of total probability and Lemma 2:

$$\mathbb { P } ( \mathcal { B } ) & = \mathbb { P } ( \mathcal { B } | \mathcal { A } ) \mathbb { P } ( \mathcal { A } ) + \mathbb { P } ( \mathcal { B } | \mathcal { A } ^ { c } ) \mathbb { P } ( \mathcal { A } ^ { c } ), \\ & \leq \mathbb { P } ( \mathcal { B } | \mathcal { A } ) ( 1 - \delta ) + \delta$$

and the fourth line uses (18). Set C = E [ L DP ] / ϵ ( -δ ) and plug it back to the formula of L , then we can organize the previous expression and into concluding when:

$$L \leq \frac { M _ { 2 } } { M _ { 1 } } + \frac { \mathbb { E } \left [ \mathcal { L } _ { \text{DP} } \right ] } { \lambda M _ { 1 } ( \epsilon - \delta ) },$$

there is

$$\mathbb { P } \left ( \lambda \| s - s ^ { \prime } \| \leq \| z - z ^ { \prime } \| \leq L \cdot \lambda \| s - s ^ { \prime } \| \right ) \geq 1 - \epsilon.$$

This proves what we desire.

## APPENDIX B

## ADDITIONAL RELATED WORKS

Spatial Data Modeling Autoencoders have been widely utilized in the analysis of spatial data due to their powerful representation learning capabilities. Traditional autoencoders (AE) have been explored in several studies [29, 48]. For instance, SEDR [53] employs a deep auto-encoder network to learn gene representations while incorporating a variational graph auto-encoder to embed spatial information simultaneously, highlighting the integration of spatial and genetic data. The combination of autoencoders with graph neural networks (GNNs) has also proven effective,

as shown in works by [31, 60, 17, 30, 59], where such models enhance spatial data representation by leveraging graph structures. Additionally, advanced probabilistic models like the hidden Markov random field have been applied to spatial data, as demonstrated by Giotto [11], which identifies spatial domains using an HMRF model with a spatial neighbor prior. This approach has been extended in studies like [55], further showcasing the potential of probabilistic methods in spatial data analysis.

Spatial Variational Autoencoders Compared to the aforementioned approaches, variational autoencoders (VAEs) are recently attracting more attention in representation learning of single-cell data analysis [12, 13]. They are preferred over traditional autoencoders because they provide a structured, continuous, and interpretable latent space that facilitates robust generative capabilities and handles uncertainty by explicitly modeling distributions. This results in more meaningful data generation and easier manipulation of latent features compared to traditional autoencoders, which can then be used for data integration [12], trajectory inference [13], etc. Furthermore, the study of VAEs lays the groundwork for extending methodologies to other variational inferencebased models, such as diffusion models and variational graph autoencoders.

A series of works have been conducted to study the interplay of VAEs with spatial information in the data, and different variants of spatial VAEs have been proposed. [3] enables learning VAE models of images that separate latent variables encoding content from rotation and translation. This can be done by directly constraining them as part of the latent space to be learned. [49] proposed to instead directly modify the latent space from the standard multivariate Gaussian distributions to low-ranked matrix-variate normal distributions to allow for the better capturing ability of the spatial information present in image data. [1] considered using VAE to learn disentangled representations of spatial networks via careful designs of the model structure and latent variable factorization. They proposed a tractable optimization algorithm to carry out variational learning. More recently, [42] considered the setting of learning spatial Gaussian priors via variational autoencoder. [21] combined graph variational autoencoder with contrastive learning.

However, a critical distinct difference to be noted is that while most existing works studying spatial VAE aim to capture the representation of the spatial information, our work aims to refine the representation learned by incorporating the spatial information given as context. This is a different but equally challenging task with broad applications in fields where spatial information is easily accessible, such as spatial transcriptomics data analysis.

## APPENDIX C ELBO DERIVATION

Below, we present the derivation of lower bound for the log probability density function (PDF) of VAE, namely the evidence lower bound (ELBO). To begin with, the marginal probability of data x can be written as

$$\log f _ { \theta } ( x ) = \log \int p _ { \theta } ( x, z ) d z,$$

where z is a latent random variable. This integral has no closed form and can be usually estimated by Monte Carlo integration with importance sampling, ie,

$$\int p _ { \theta } ( x, z ) d z = \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \frac { p _ { \theta } ( x, z ) } { q ( z | x ) } \right ].$$

Here q z x ( | ) is the proposed variational distribution, where we can draw sample z from this distribution given x and a . Therefore, by Jensen's inequality, we can find evidence lower bound (ELBO) of the conditional PDF:

$$\cdots \cdots \cdots \cdots \cdots \\ \log f _ { \theta } ( x ) & = \log \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \frac { p _ { \theta } ( x, z ) } { q ( z | x ) } \right ] \\ & \geq \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log \frac { p _ { \theta } ( x, z ) } { q ( z | x ) } \right ]. \\ \mathbb { m } \quad \cdots \cdots \cdots \cdots \cdots \cdots$$

Using Bayes rule, the ELBO can be equivalently expressed as:

$$\cdot \quad \text{as} \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \\ \cdot \quad \cdot \quad \begin{array} { c c c c } & \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log \frac { p _ { \theta } ( x, z ) } { q ( z | x ) } \right ] \\ \cdot \quad \quad = \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log \frac { p _ { \theta } ( x | z ) p _ { \theta } ( z ) } { q ( z | x ) } \right ] \\ \cdot \quad \quad = \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log \frac { p _ { \theta } ( z ) } { q ( z | x ) } \right ] + \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log p _ { \theta } ( x | z ) \right ] \\ \cdot \quad \quad = - \, D _ { \text{KL} } ( q ( z | x ) | | p _ { \theta } ( z ) ) + \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log p _ { \theta } ( x | z ) \right ]. \\ \, \text{This concludes the derivation}. \end{array}$$

This concludes the derivation.

## APPENDIX D KL DIVERGENCE

Lemma 4 (KL divergence of multivariate normals) . Suppose both p and q are the pdfs of multivariate normal distributions N d ( µ 1 , Σ 1 ) and N d ( µ 2 , Σ 2 ) , respectively. The Kullback-Leibler distance from p to q is given by

$$\Phi _ { K L } ( p | | q ) & = \frac { \mu } { 2 } \left [ \log \frac { \mu } { | \Sigma _ { 2 } | } - d + \text{tr} \{ \Sigma _ { 2 } ^ { - 1 } \Sigma _ { 1 } \} \\ & \quad + ( \mu _ { 2 } - \mu _ { 1 } ) ^ { T } \Sigma _ { 2 } ^ { - 1 } ( \mu _ { 2 } - \mu _ { 1 } ) \right ].$$

When Σ 1 = diag( σ , . . . , σ 2 1 2 d ) and Σ 2 = I d , it reduces to

$$\underset { \text{to cap} } { \overset { \dots \dots \dots } { \text{ing} } } \text{ } \ D _ { K L } ( p | | q ) = - \frac { 1 } { 2 } \left [ d + \sum _ { j = 1 } ^ { d } \log \sigma _ { j } ^ { 2 } - \sum _ { j = 1 } ^ { d } \sigma _ { j } ^ { 2 } - \| \mu _ { 2 } - \mu _ { 1 } \| _ { 2 } ^ { 2 } \right ].$$

Proof of Theorem 4. By definition, we have

$$\i \, \text{to cap} ^ { - } _ { \ } \, \text{for theorem $4$} \, \text{By definition, we have} \\ \text{ordering} \\ \text{different} \\ \text{in} \\ \text{such as} \, \text{in} \\ \text{such as} \, \text{in} \\ \text{in} \\ \text{for the} \\ \text{nely the} \\ \text{marginal} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \text{in} \\ \end{array} \, \text{for the} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \\ \text{in} \\ \text{nely the} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \, \text{in} \\ \text{neq} \,\text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{in} \\ \text{neq} \, \text{in} \, \text{neq} \, \text{in} \\ \text{neq} \, \text{in} \, \text{neq} \, \text{in} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq } \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{neq} \, \text{neq} \\ \text{nege} \, \text{neq} \\ \text{nege} \, \text{neq} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{\text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \, \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{peq} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \text{nege} \\ \end{cases}$$

The last expression on the right-hand side translates to:

$$\frac { 1 } { 2 } \left [ \log \frac { | \Sigma _ { 2 } | } { | \Sigma _ { 1 } | } - d + \text{tr} \{ \Sigma _ { 2 } ^ { - 1 } \Sigma _ { 1 } \} + ( \mu _ { 2 } - \mu _ { 1 } ) ^ { T } \Sigma _ { 2 } ^ { - 1 } ( \mu _ { 2 } - \mu _ { 1 } ) \right ],$$

$$& \frac { 1 } { 2 } \left [ \log \frac { | \Sigma _ { 2 } | } { | \Sigma _ { 1 } | } - d + \text{tr} \{ \Sigma _ { 2 } ^ { - 1 } \Sigma _ { 1 } \} + ( \mu _ { 2 } - \mu _ { 1 } ) ^ { T } \Sigma _ { 2 } ^ { - 1 } ( \mu _ { 2 } - \mu _ { 1 } ) \right ], \\ & \text{which finishes the proof.} \quad \Box \quad \text{where $v,\,$s$}$$

## APPENDIX E ADDITIONAL EXPERIMENT RESULTS

In this section, we show an additional earlier version of experiment result evaluating the spatial autocorrelation and isometry preserving property of our proposed dp-VAE model. Below we introduce the two main evaluation metrics used.

- · Moran's I [28] is a correlation coefficient that measures the overall spatial autocorrelation of a dataset. Intuitively, for a given gene, it measures how one spot is similar to others surrounding it. If the spots are attracted (or repelled) by each other, it implies they are not independent. Thus, the presence of autocorrelation indicates the spatial pattern of gene expression. Moran's I value ranges from -1 to 1, where a value close to 1 indicates a clear spatial pattern, a value close to 0 indicates random spatial expression, and a value close to -1 indicates a chessboard-like pattern. To evaluate the spatial variability of a given gene, we calculate the Moran's I using the following formula,

$$I = \frac { N } { W } \frac { \sum _ { i } \sum _ { j } \left [ w _ { i j } \left ( x _ { i } - \bar { x } \right ) \left ( x _ { j } - \bar { x } \right ) \right ] } { \sum _ { i } \left ( x _ { i } - \bar { x } \right ) ^ { 2 } },$$

where x i and x j are the gene expression of spots i and j , ¯ x is the mean expression of the feature, N is the total number of spots, w ij is the spatial weight between spots i and j calculated using the 2D spatial coordinates of the spots, and W is the sum of w ij . We select the k nearest neighbors for each spot using spatial coordinates. Moran's I statistic is robust to the choice of k and is set at 5 in our analysis. We assign w ij = 1 if spot j is in the nearest neighbors of spot i , and 0 otherwise.

- · Geary's C is another commonly used statistic for measuring spatial autocorrelation. It is calculated as:

$$C = \frac { N - 1 } { 2 W } \frac { \sum _ { i } \sum _ { j } \left [ w _ { i j } \left ( x _ { i } - x _ { j } \right ) ^ { 2 } \right ] } { \sum _ { i } \left ( x _ { i } - \bar { x } \right ) ^ { 2 } },$$

The value of Geary's C ranges from zero to two, where zero indicates perfect positive autocorrelation.

Evaluation Protocals. Moran's I [28] is a correlation coefficient that measures the overall spatial autocorrelation of a dataset. Intuitively, for a given gene, it measures how one spot is similar to others surrounding it. If the spots are attracted (or repelled) by each other, it implies they are not independent. Thus, the presence of autocorrelation indicates the spatial pattern of gene expression. Moran's I value ranges from -1 to 1, where a value close to 1 indicates a clear spatial pattern, a value close to 0 indicates random spatial expression, and a value close to -1 indicates a chessboard-like pattern. To evaluate the spatial variability of a given gene, we calculate the Moran's I using the following formula,

$$I = \frac { N } { W } \frac { \sum _ { i } \sum _ { j } \left [ w _ { i j } \left ( x _ { i } - \bar { x } \right ) \left ( x _ { j } - \bar { x } \right ) \right ] } { \sum _ { i } \left ( x _ { i } - \bar { x } \right ) ^ { 2 } },$$

where x i and x j are the gene expression of spots i and j , ¯ x is the mean expression of the feature, N is the total number of spots, w ij is the spatial weight between spots i and j calculated using the 2D spatial coordinates of the spots, and W is the sum of w ij . We select the k nearest neighbors for each spot using spatial coordinates. Moran's I statistic is robust to the choice of k and is set at 5 in our analysis. We assign w ij = 1 if spot j is in the nearest neighbors of spot i , and 0 otherwise.

Geary's C is another commonly used statistic for measuring spatial autocorrelation. It is calculated as:

$$C = \frac { N - 1 } { 2 W } \frac { \sum _ { i } \sum _ { j } \left [ w _ { i j } \left ( x _ { i } - x _ { j } \right ) ^ { 2 } \right ] } { \sum _ { i } \left ( x _ { i } - \bar { x } \right ) ^ { 2 } }, \\. \quad. \quad. \quad. \quad.$$

The value of Geary's C ranges from zero to two, where zero indicates perfect positive autocorrelation.

The experiment results are shown in Figure 12, Figure 13, Table 1. Figure 12 shows the gene expression levels of six marker genes extracted by VAE. For each gene in each panel, the left subplot shows the gene expression levels in the spatial coordinates, and the right subplot shows the 2DUMAP of the latent representation of the corresponding model. The cells are dispersed in the latent space. The olfactory bulb region, with genes Calb2, Gng4, S100a5, and Febp7 enriched, is well separated from other regions. Table 1 shows that enforcing distance-preserving property induces stronger spatial autocorrelations on out-of-sample datasets. Figure 13 shows that the proposed method can improve isometric property of the learned latent space.

## APPENDIX F ADDITIONAL BACKGROUND OF VAE

CVAE. A conditional variational encoder (CVAE) [43, 52] is an extension of VAE, where an additional conditional input y ∈ Y can be specified in the encoder and decoder network

$$q _ { \theta } & \colon \mathcal { X } \times \mathcal { Y } \to \mathcal { Z }, \\ p _ { \phi } & \colon \mathcal { Z } \times \mathcal { Y } \to \mathcal { X }.$$

Similar to VAE, a typical parameterization of CVAE specifies the encoder and decoder functions as conditional Gaussian distributions, where the mean and variance are parameterized by neural networks.

β -VAE. β -VAE [19] incorporates a temporature term of the KL loss term in the VAE model:

$$\Gamma _ { s } ^ { r } - \beta \cdot D _ { \text{KL} } ( q ( z | x ) | | p _ { \theta } ( z ) ) + \mathbb { E } _ { z \sim q ( \cdot | x ) } \left [ \log p _ { \theta } ( x | z ) \right ].$$

Intuitively, a larger KL weight forces the VAE to confine to the prior distribution more, and a smaller KL weight allows the VAE to focus more on minimizing the reconstruction loss.

scVI (single-cell Variational Inference). scVI [33] is a leading deep generative model for analyzing single-cell RNA sequencing (scRNA-seq) data. It employs a CVAE architecture

Spatial feature of Calb2

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-5adbb557-12f8-4c19-afa5-71a1b170251c.jpg)

Spatial feature of Gng4

Spatial feature of Nrgn

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-29d01557-3fd9-4a72-b334-75878758e4a9.jpg)

Latent representation of Calb2

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-e67ff7df-ef6c-4074-967f-caa3c2779579.jpg)

Latent representation of Gng4

Latent representation of Nrgn

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-0fccb63d-8dae-4984-be0d-26a20e5ae579.jpg)

Spatial feature of Ttr

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-3915aa8c-7f2c-4972-8e3d-d0b1fd11ac8a.jpg)

Spatial feature of S100a5

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-d434c69f-bcba-4425-8e29-1b9e124a8a47.jpg)

Spatial feature of Fabp7

9

8

7

6

5

4

3

2

1

0

4.0

3.5

3.0

2.5

2.0

1.5

1.0

0.5

0.0

Latent representation of Ttr

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-db03c914-1299-4451-b99b-52487f10cd77.jpg)

Latent representation of S100a5

Latent representation of Fabp7

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-30a802bf-d3f5-4c90-8be6-2c02a3652b0c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-69919e32-af0c-4a9e-8b8b-84324379774d.jpg)

Fig. 12: Latent representation extracted from VAE trained on the Mouse Brain Serial Datasets anterior 1. The olfactory bulb region, with genes Calb2, Gng4, S100a5, and Febp7 enriched, is well separated from other regions

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-7b91bd4c-63e4-4418-b746-6ba708ce3f23.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-a149439d-bb1e-451f-9414-3c374e8e0bd7.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-47339555-74a9-41c9-bd9c-575ecf76d9a3.jpg)

TABLE 1: Moran's I and Geary's C of the latent representation extracted by scVI and VAE using 4 mouse brain cell test datasets, with and without distance-preserving penalty, averaged over 5 repeated trials to account for the randomness of the training process.

|            | Method   | Moran's I   | Moran's I   | Moran's I   | Moran's I   | Geary's C   | Geary's C   | Geary's C   | Geary's C   |
|------------|----------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|
|            | Method   | A2          | A1          | P2          | P1          | A2          | A1          | P2          | P1          |
| Vanilla    | VAE      | 0.62(0.07)  | 0.55(0.05)  | 0.52(0.05)  | 0.52(0.03)  | 0.36(0.06)  | 0.41(0.03)  | 0.49(0.05)  | 0.43(0.03)  |
|            | scVI     | 0.43(0.03)  | 0.52(0.04)  | 0.37(0.02)  | 0.45(0.04)  | 0.57(0.03)  | 0.48(0.04)  | 0.62(0.02)  | 0.55(0.04)  |
| Distance   | VAE      | 0.64(0.02)  | 0.60(0.03)  | 0.56(0.06)  | 0.49(0.06)  | 0.35(0.02)  | 0.37(0.02)  | 0.45(0.07)  | 0.46(0.05)  |
| Preserving | scVI     | 0.45(0.04)  | 0.52(0.04)  | 0.43(0.02)  | 0.47(0.03)  | 0.55(0.05)  | 0.48(0.04)  | 0.57(0.02)  | 0.53(0.03)  |

to model gene expression patterns across individual cells. In scVI, the input consists of gene expression counts for each cell, while the model aims to learn a latent representation that captures biological variability. The model accounts for technical factors such as library size and batch effects, which are treated as observed variables. Specifically, scVI uses the normalized gene expression counts as the target variable, with the library size included as a nuisance variable to be inferred. Batch information, when available, is incorporated as conditional information to help disentangle technical from biological variability. This approach allows scVI to perform various downstream tasks, including normalization, batch correction, imputation, and dimensionality reduction, making it a versatile tool for single-cell transcriptomics analysis. Markovian Hierarchical VAE. A Markovian Hierarchical Variational Autoencoder (MHVAE) [24, 44] is a generalization of a VAE that extends to multiple hierarchies over latent variables. Under this formulation, latent variables are interpreted as generated from other higher-level, more abstract patterns. Let T &gt; 1 denote the total number of hierarchical levels; the objective of an MHVAE is similar to VAE, except that there are now T latent variables whose joint distribution can be decomposed into Markov transitions:

$$p _ { \theta } ( x ) = \int p _ { \theta } ( x | z _ { 1 \colon T } ) d z _ { 1 \colon T }$$

$$= \overset { \kappa } { \int } p _ { \theta } ( x | z _ { 1 } ) \prod _ { t = 2 } ^ { T } p _ { \theta } ( z _ { t } | z _ { t + 1 } ) d z _ { 1 \colon T }. \quad ( 2 0 )$$

Deriving for the ELBO of the MHVAE can follow a similar

Fig. 13: Visualization of latent representation space obtained from scVI (left) and scVI regularized with distancepreserving loss (right). More isotropic and homogeneous ellipses indicate more distance-preserving [56].

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-6bd3e60e-5a3c-4fef-81b9-275f6b6922a5.jpg)

## procedure as of VAE:

$$\begin{array} { l } \cdot \\ p _ { \theta } ( x ) \geq \mathbb { E } _ { q _ { \phi } ( z _ { 1 } | x ) } \left [ \log p _ { \theta } ( y | z _ { 1 } ) \right ] \\ & - \mathbb { E } _ { q _ { \phi } ( z _ { T - 1 } | x ) } \left [ D _ { \text{KL} } \left ( q _ { \phi } \left ( z _ { T } | z _ { T - 1 } \right ) \| p ( z _ { T } ) \right ) \right ] \\ & - \sum _ { t = 1 } ^ { T - 1 } \mathbb { E } _ { q _ { \phi } ( z _ { t - 1 }, z _ { t + 1 } | x ) } \left [ D _ { \text{KL} } ( q _ { \phi } ( z _ { t } | z _ { t - 1 } ) \| p _ { \theta } ( z _ { t } | z _ { t + 1 } ) \right ) \right ]. \\ \text{with } z _ { \lambda } \ = \ 0 \text{ for notational communication.                                                                                                                                                                                                        
                                                                                                                                                                                                        }                                                                                                                                                                                                        \end{array}$$

$$\cdot$$

with z 0 ≡ 0 for notational convenience. The first term (reconstruction term) and the second term (prior matching term) can be interpreted similarly. The third term (consistency term) is new. It endeavors to make the distribution at z t consistent from both forward and backward processes. All terms are taken with respect to distributions that we can sample from, and Monte Carlo estimation can be obtained.

However, a key caveat is that the expectation over two random variables { z t -1 , z t +1 } for every timestep, the variance of its Monte Carlo estimate could potentially be higher than a term that is estimated using only one random variable per timestep. As it is computed by summing up T -1 consistency terms, the final estimated value of the ELBO may have a high variance for large T values. Instead, we can derive an equivalent version that does not pose this issue:

$$p _ { \theta } ( y ) & \geq \mathbb { E } _ { q _ { \phi } ( z _ { 1 } | x ) } \left [ \log p _ { \theta } ( x | z _ { 1 } ) \right ] - D _ { \text{KL} } \left ( q _ { \phi } \left ( z _ { T } | x \right ) \| p ( z _ { T } ) \right ) \\ & \quad - \sum _ { t = 2 } ^ { T - 1 } \mathbb { E } _ { q _ { \phi } ( z _ { t } | x ) } \left [ D _ { \text{KL} } ( q _ { \phi } ( z _ { t } - 1 | z _ { t }, x ) \| p _ { \theta } ( z _ { t - 1 } | z _ { t } ) \right ) \right ]. \\ \longrightarrow \quad \cdot \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

The changed last term is called denoising matching term following the variational diffusion model literature [34]. We learn desired denoising transition step p θ ( z t -1 | z t ) as an approximation to tractable, variational denoising transition step q ϕ ( z t -1 | z , x t ) .

Fig. 14: 10-lengthed-genes aggregated expression heatmap of the mouse brain cell. Gene expressions show strong spatial trends and local correlations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00911v2-figure-2c11ab47-e5bd-4f9e-9496-4fab24822b2c.jpg)
Wenbin Zhou Jin-Hong Du
Distance-Preserving Spatial Representations in Genomic Data
The spatial context of single-cell gene expression data is crucial for many
downstream analyses, yet often remains inaccessible due to practical and
technical limitations, restricting the utility of such datasets. In this paper,
we propose a generic representation learning and transfer learning framework
dp-VAE, capable of reconstructing the spatial coordinates associated with the
provided gene expression data. Central to our approach is a distance-preserving
regularizer integrated into the loss function during training, ensuring the
model effectively captures and utilizes spatial context signals from reference
datasets. During the inference stage, the produced latent representation of the
model can be used to reconstruct or impute the spatial context of the provided
gene expression by solving a constrained optimization problem. We also explore
the theoretical connections between distance-preserving loss, distortion, and
the bi-Lipschitz condition within generative models. Finally, we demonstrate
the effectiveness of dp-VAE in different tasks involving training robustness,
out-of-sample evaluation, and transfer learning inference applications by
testing it over 27 publicly available datasets. This underscores its
applicability to a wide range of genomics studies that were previously hindered
by the absence of spatial data.