soScope frameworkA generative model for resolution enhancement of spatial omicssoScope utilizes three modalities from the same tissue: spot-level omics profiles, spatial neighboring relations of spots, and subspot-level morphological features. We aim to divide the observed expression profile \({{{{\bf{x}}}}}^{(n)}\in {{\mathbb{R}}}^{G}\) at the \(n\)th spot \({s}^{(n)}\) into \(K\) subspots \({s}_{1:K}^{(n)}\):$${{{{\bf{x}}}}}^{{(}n{)}}={\sum}_{k=1}^{K}{\widehat{{{{\bf{x}}}}}}_{k}^{{(}n{)}}+{{{{\mathbf{\epsilon }}}}}^{{(}n{)}},\, {n}=1\ldots N$$
(1)
where \({\hat{{{{\bf{x}}}}}}_{1:K}^{(n)}\in {{\mathbb{R}}}^{G}\) are the latent expression to be inferred at \(K\) subspots, \({{{{\mathbf{\epsilon }}}}}^{(n)}\) is drawn from a Gaussian noise with a mean of \({{{\mathbf{0}}}}\) and variance of \({\sigma }^{2}{{{\boldsymbol{I}}}}\), and \(N\) is the total number of spots in a spatial dataset. We formulate the task of spatial omics enhancement into a probabilistic generative model:$${{{{\bf{z}}}}}^{(n)} \sim {{{\rm{Gaussian}}}}\left({{{\mathbf{0}}}},\, {\sigma }^{2}{{{\boldsymbol{I}}}}\right)\\ {{{{\mathbf{\omega }}}}}_{k}^{(n)}=f\left({{{{\bf{y}}}}}_{k}^{(n)},\, {{{{\bf{z}}}}}^{(n)}\right),\, k=1\ldots K\\ {\widehat{{{{\bf{x}}}}}}_{k}^{(n)}{{{\rm{|}}}}{{{{\bf{z}}}}}^{(n)},\, {{{{\bf{y}}}}}_{k}^{(n)} \sim P\left({{{{\mathbf{\omega }}}}}_{k}^{(n)}\right),\, k=1\ldots K\\ {{{{\bf{x}}}}}^{(n)}{{{\rm{|}}}}{\widehat{{{{\bf{x}}}}}}_{1:K}^{(n)} \sim {{{\rm{Gaussian}}}}\left({\sum }_{k}{\widehat{{{{\bf{x}}}}}}_{k}^{(n)},\, {\sigma }^{2}{{{\boldsymbol{I}}}}\right)$$
(2)
where \({{{{\bf{z}}}}}^{(n)}\in {{\mathbb{R}}}^{D}\) is the latent representation of the underlying spatial omics state at spot \({s}^{(n)}\), and \({{{{\bf{y}}}}}_{k}^{(n)}{{{\boldsymbol{\in }}}}{{\mathbb{R}}}^{L}\) is the feature extracted from the corresponding image region at subspot \({s}_{k}^{(n)}\). A neural network \(f(\cdot )\) combines subspot-level image feature \({{{{\bf{y}}}}}_{k}^{(n)}\) and the spot-level latent state \({{{{\bf{z}}}}}^{(n)}\) to learn the parameter \({{{{\mathbf{\omega }}}}}_{k}^{(n)}\) for the probability \(P(\cdot )\) that models expression profiles at the subspot level. The generative probability \(P(\cdot )\) is selected with respect to spatial omics data types to best describe its property. Finally, the profile \({{{{\bf{x}}}}}^{(n)}\) of a spot can be obtained by aggregating profiles at its corresponding subspots jittered by the zero-mean Gaussian noise with \({\sigma }^{2}\) as the variance.For simplicity, we use matrix notations \({{{\boldsymbol{X}}}}{{{\boldsymbol{=}}}}{\left[{{{{\bf{x}}}}}^{(1)}\,{{\ldots }}\,{{{{\bf{x}}}}}^{(N)}\right]}^{T}{{\in }}{{\mathbb{R}}}^{N\times G}\), \({{{\boldsymbol{Y}}}}{{{\boldsymbol{=}}}}{\left[{{{{\bf{y}}}}}_{1}^{(1)}\,{{\ldots }}\,{{{{\bf{y}}}}}_{K}^{(N)}\right]}^{T}{{\in }}{{\mathbb{R}}}^{{NK}\times L}\), \(\widehat{{{{\boldsymbol{X}}}}}{{=}}{\left[{\widehat{{{{\bf{x}}}}}}_{1}^{(1)}\,{{\ldots }}\,{\widehat{{{{\bf{x}}}}}}_{K}^{(N)}\right]}^{T}{{\in }}{{\mathbb{R}}}^{{NK}\times G}\), \({{{\boldsymbol{Z}}}}{{{\boldsymbol{=}}}}{\left[{{{{\bf{z}}}}}^{(1)}\,{{\ldots }}\,{{{{\bf{z}}}}}^{(N)}\right]}^{T}{{\in }}{{\mathbb{R}}}^{N\times D}\) and obtain the log-likelihood form for the resolution enhancement process in Eq. 2:$$\log p\left({{{\boldsymbol{X}}}} | {{{\boldsymbol{Y}}}}\right)=\log \left\{{\sum}_{\widehat{{{{\boldsymbol{X}}}}}}{\sum}_{{{{\boldsymbol{Z}}}}}p\left({{{\boldsymbol{X}}}} | \widehat{{{{\boldsymbol{X}}}}}\right)p\left(\widehat{{{{\boldsymbol{X}}}}} | {{{\boldsymbol{Y}}}},\, {{{\boldsymbol{Z}}}}\right)p\left({{{\boldsymbol{Z}}}}\right)\right\}$$
(3)
Variational inference of soScopeThe summation over latent variables makes calculating \(\log p\left({{{\boldsymbol{X}}}} | {{{\boldsymbol{Y}}}}\right)\) intractable. Therefore, we solve the log-likelihood in Eq. 3 via the variational Bayesian inference40. We hypothesize that latent state \({{{\boldsymbol{Z}}}}\) is related to spatial omics profile \({{{\boldsymbol{X}}}}\) and spatial neighboring relations \({{{\boldsymbol{A}}}}\). Then, we estimate its variational distribution in the form of \(q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)\), where \({{{\boldsymbol{A}}}}\in {[{{\mathrm{0,1}}}]}^{N\times N}\) is the spatial neighboring relation graph41 and is calculated using profiles between neighboring spots in \({{{\boldsymbol{X}}}}\) (Supplementary Note 1). The evidence lower bound (ELBO) is formulated as (Supplementary Note 2):$${{{\rm{ELBO}}}}={{\mathbb{E}}}_{q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)}\left\{{{\mathbb{E}}}_{p({\widehat{{{{\boldsymbol{X}}}}}} | {{{\boldsymbol{Y}}}},\,{{{\boldsymbol{Z}}}})}\left[\log p\left({{{\boldsymbol{X}}}} | \widehat{{{{\boldsymbol{X}}}}}\right)\right]\right\}-{D}_{{KL}}\left[q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}} \,{{{\boldsymbol{A}}}}\right)\parallel p\left({{{\boldsymbol{Z}}}}\right)\right]$$
(4)
where \({D}_{{KL}}\left[q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}},\,{{{\boldsymbol{A}}}}\right)\parallel p({{{\boldsymbol{Z}}}})\right]\) is the Kullback-Leibler (KL) divergence between \(q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)\) and \(p({{{\boldsymbol{Z}}}})\).Image regularizationTo ensure that the enhanced profiles accurately reflect morphological similarities observed in the images, we introduce an additional image regularization term to enforce their consistency. In detail, we calculate the similarity matrices among subspots based on enhanced omics profiles \({{{\boldsymbol{\Lambda }}}}\in {[{{\mathrm{0,1}}}]}^{{NK}\times {NK}}\) and image features \({{{\boldsymbol{W}}}}\in {[{{\mathrm{0,1}}}]}^{{NK}\times {NK}}\), respectively (Supplementary Note 1). Then, we use a soft-label cross entropy42 to evaluate their consistency:$${L}_{{{{\rm{Image}}}}}= {{{\rm{CE}}}}\left({{{\boldsymbol{\Lambda }}}},{{{\boldsymbol{W}}}}\right)\\= -\frac{1}{\left({NK}\right)^{2}}{\sum}_{i,\, j=1}^{{NK}}\left[{w}^{(i,\, j)}\log {\lambda }^{(i,\, j)}+\left(1-{w}^{(i,\, j)}\right)\log \left(1-{\lambda }^{(i,\, j)}\right)\right]$$
(5)
where \({w}^{(i,j)}\) and \({\lambda }^{(i,j)}\) are the entry at the \(i\)th row and the \(j\)th column in \({{{\boldsymbol{W}}}}\) and \({{{\boldsymbol{\Lambda }}}}\), respectively. By minimizing the \({L}_{{{{\rm{Image}}}}}\), enhanced omics profiles are encouraged to reflect the morphological patterns from the images.soScope objective optimization and inferenceWe formulate the overall optimization function by combining ELBO and image regularization:$${LL}={{{\rm{ELBO}}}}-\beta {L}_{{{{\rm{Image}}}}}$$
(6)
Here, \(\beta\) is a hyperparameter balancing the level of constraint (we set \(\beta=1\) for all experiments). In the implementation of soScope (Fig. 1), we choose \(q\left({{{\boldsymbol{Z}}}} | {{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)\) as a Gaussian distribution, with its mean and covariance matrix given by a graph encoder \(h\left({{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)\). The parameters modeling \(p\left(\widehat{{{{\boldsymbol{X}}}}} | {{{\boldsymbol{Y}}}},\, {{{\boldsymbol{Z}}}}\right)\) are given by an enhancement decoder \(f({{{\boldsymbol{Y}}}},\, {{{\boldsymbol{Z}}}})\) (Supplementary Fig. 1).The model is optimized with the Adam optimizer43 using a two-step strategy:(1) Initialize the network \(h\left({{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)\) using a simple variational graph auto-encoder framework44 without resolution enhancement (Supplementary Note 3).(2) Maximize the overall target function (Eq. 6) to infer subspot profiles.After optimization, we use the expectation of subspot-level profiles \({\widehat{{{{\boldsymbol{X}}}}}}^{*}\) as the enhanced profiles:$$\begin{array}{c}{{{{\boldsymbol{Z}}}}}^{*}={{\mathbb{E}}}_{q\left({{{\boldsymbol{Z}}}}{{{\rm{|}}}}{{{\boldsymbol{X}}}},\, {{{\boldsymbol{A}}}}\right)}\left[{{{\boldsymbol{Z}}}}\right]\\ {\widehat{{{{\boldsymbol{X}}}}}}^{*}={{\mathbb{E}}}_{p\left(\widehat{{{{\boldsymbol{X}}}}}{{{\rm{|}}}}{{{{\boldsymbol{Z}}}}}^{*},\, {{{\boldsymbol{Y}}}}\right)}\left[\widehat{{{{\boldsymbol{X}}}}}\right]\end{array}$$
(7)
Extension to different omics typesTo properly model the generative process, \(P{(\cdot )}\) is determined by the spatial omics type. For transcript data, we use the negative binomial distribution (NB)45; for protein and histone modification, we use the Poisson distribution37,46; for other spatial omics profiles without conclusive probabilistic distribution knowledge (such as principal components (PCs) of DNA, and normalized ATAC peaks), we use the Gaussian distribution (Supplementary Table 1) as a general choice.Extension to spatial multiomics resolution enhancementTo enable soScope to simultaneously enhance the resolution of multiomics profiles, we modify the original model framework with three major changes: (1) constructing the spatial neighboring relation graph \({{{\boldsymbol{A}}}}\) based on multiomics similarities; (2) mapping latent representations \({{{\boldsymbol{Z}}}}\) using a graph encoder with multiomics inputs; and (3) generating distribution parameters for multiomics subspot profiles simultaneously from the decoder network. The elaborated description of the multiomics soScope framework is provided in Supplementary Note 4.DatasetsHuman intestine dataset from Visium platformThe human intestine tissue section was collected from the colon of a 66-year-old male individual (labeled as A1 in the original data). The tissue was profiled by Visium platform and contained 2,649 sequencing spots. Tissue region annotations and their regional markers were obtained from the previous publication34. For the enhancement analysis, 9 marker genes from 3 regions (epithelium, muscularis, and immune) were included as input. For the image, we took the coordinate information for each subspot and segmented the corresponding H&E image region, with the radius of the region sett to 5 times the spot radius. Then, we resized image patches to 299×299 pixel size and fed them into a pretrained Google Inception-v333 on ImageNet to obtain image features (in 2,048 dimensions).Mouse head dataset from Xenium platformThe mouse head section was obtained from a one-day-old mouse pup and profiled using the Xenium platform. This region contains 457,781 cells, with each expressing 379 genes. For biological variability analysis, expression counts were normalized against total counts, scaled by median gene expressions, and log1p transformed. These data were then processed using the modelGeneVar function in the scran R package (v1.20.1). Only genes with positive biological variability (n = 194) were selected for gene correlation analysis and categorized into high (top 60 genes), medium (61–120), and low (121–194) variability groups. For the imaging modality, deep features were extracted following the pipeline used in the intestine data.Mouse kidney dataset from Xenium platformThe mouse kidney dataset is measured from the same mouse pup section using Xenium platform. This tissue region covers 1538 cells and 379 genes. For the enhancement analysis, we selected the top 5 abundant genes and the top 5 marker genes. For the image modality, we segmented the image for each cell from the tissue H&E and extracted deep features using the same method as in our intestinal data analysis. As the cell locations were not arranged in a regular array, we defined the neighborhood of each cell as the five spatially nearest cells using the NearestNeighbors function from the scikit-learn.neighbors Python package (v1.2.0) for the purpose of graph construction.Human ductal carcinoma dataset from Xenium platformThe dataset is measured from the breast section of a human female using Xenium platform. This tissue region covers 37,894 cells and 313 genes. We identified three distinct gene groups for targeted analyses. For highly correlated genes, we selected the top 10 tumor marker genes as reported in the prior study5. To identify lowly correlated genes, we employed a linear regression model using the LinearRegression function from the scikit-learn.linear_model Python package (v1.2.0), comparing the top 250 abundant genes against the top 50 principal components of image features. We then screened the bottom 10 genes with the lowest \({R}^{2}\), determined via the r2_score function from the scikit-learn.metrics Python package (v1.2.0). For immune cell analysis, we screened 11 marker genes across three cell types (four for T cells, two for B cells, and five for macrophages) based on the previous publication5, and quantified cell density as the summation of log1p-normalized expression of these marker genes following the setting of BayesSpace30. Image modality features were processed using the same pipeline as utilized in the intestinal data analysis, while for soScope without image modality, standard Gaussian noise was employed as a substitute.Mouse embryo dataset from spatial-CUT&Tag platformThe mouse embryo dataset was collected from a mouse embryo at embryonic day 11. The section was profiled by the spatial-CUT&Tag platform and contained 1,974 sequencing spots. In the original publication18, 11 subpopulations were reported. We selected four organ regions with clear spatial structures (liver, heart, forebrain, and spinal cord) and identified the top 15 variable peaks from each region (getMarkers function in archR v1.0.1 R package) for the resolution enhancement analysis. For the image modality, we followed the established pipeline and learned image features from the H&E image of the embryo. To estimate gene activities from peak counts, we mapped peak regions to genes using bedtools R package (v2.28.0) and performed a negative log transformation on peak data. A higher activity score represents a lower suppression for a gene.Human embryonic heart dataset from spatial transcriptomics platformThe human heart tissue section was collected from a human embryo at 6.5 post-conception weeks (PCW). The tissue was profiled by the ST platform and contained 186 sequencing spots. Gene expressions were scaled using total transcript counts and log1p-transformed. Regional markers FHL2 and LDHA were identified and reported in the original study25.Mouse liver dataset from slide-DNA-seq and slide-RNA-seq platformsThe dataset was obtained from a mouse liver metastasis section, encompassing two distinct tumor clone regions labeled as clone A and clone B. The dataset included 24,679 spots for slide-DNA-seq data, and 21,902 spots for slide-RNA-seq data. We used the processed data provided in the original work19. For the resolution enhancement purpose, we divided the tissue slide into 60×60 squared regions and averaged original data in each region into spot-level profiles. Then, we included the top 2 DNA PCs and three marker genes (Hmga2, Tm4sf1, and Aqp5) highlighted in the original work, respectively, in the analysis. For the H&E image data, we extracted deep features from the H&E image for subspots as previously described.Human skin dataset from spatial-CITE-seq platformThe human skin tissue section was collected from an adult donor who received an early immune activation due to the administration of a Coronavirus Disease 2019 (COVID-19) mRNA vaccine. The tissue contained 1,618 sequencing spots, with each spot including 15,486 genes and 283 proteins. Both omics were normalized using the SCTransform function from R package sctransform (v0.3.5). After that, all proteins and 114 top variable genes (the union of top 20 overexpressed genes from each subpopulation identified using analytical pipeline in the original study15) were included for the enhancement analysis. For the image modality, we extracted deep features from the bright-field image following the same pipeline described above.Mouse embryo dataset from spatial ATAC-RNA-seq platformThe mouse embryo tissue section was from an E13 (embryonic day 13) mouse embryo. The section included 1,874 sequencing spots, each containing data for 15,748 genes and 24,071 ATAC peaks. We selected the 25 most variable genes (modelGeneVar function in scran v1.20.1 R package) and their corresponding peaks for the enhancement analysis. For the ATAC modality, we performed normalization using the SCTransform function from the R package sctransform (v0.3.5) and subsequently applied min-max normalization. For the image modality, deep features were extracted from the bright-field images using the previously described pipeline.Experiment setupIn silico “low-resolution” profile simulationIn performance benchmarks, we implemented a simulation strategy wherein we combined neighboring spots to create “low-resolution” spots. Then, we assessed the ability of resolution enhancement techniques to accurately restore the original spot profiles. In the case of Visium datasets, where spots were arranged in a hexagonal pattern, we aggregated the profiles of seven spots within the same hexagon to form a single “low-resolution” spot, with the centroid of the hexagon as its position. For the Xenium datasets, to simulate ST data, we divided the whole tissue regions into the small squared regions (in total 240×240 for mouse head, 50×50 for mouse kidney, and 39×39 for human ductal carcinoma) and integrated the gene expressions of all cells within these square regions to generate regional aggregated expression profiles. For datasets from other spatial platforms, we merged neighboring spots within a square region to produce these “low-resolution” spots. For spatial-CUT&Tag, we merged spots in every 2×2 region. In the multi-scale enhancement experiments with slide-DNA-seq and slide-RNA-seq datasets, we varied the region size from 2×2 to 6×6.Subpopulation identificationFor human skin dataset, we followed the analysis in the original study15 and employed the Louvain algorithm47 (resolution=0.75) to identify subpopulations from the top 10 PCs of normalized proteins or genes.Comparing methodsAll compared methods were benchmarked with the same input features on the same Linux server (Ubuntu 20.04.3 LTS operation system) with Xeon(R) 6226 R CPU and NVIDIA GeForce RTX 3090 GPU (Driver Version: 470.63.01, CUDA Version: 11.4).Spatial-based enhancement approaches
BayesSpace
(v1.2.1R package)30 is a statistical method for resolution enhancement of ST data. It uses the neighboring expression information to estimate the subspot-level gene expressions. We used the spatialEnhance function from BayesSpace R package for resolution enhancement.
Linear interpolation
(Spatial linear) calculates subspot-level profiles by averaging profiles from the nearest neighboring spots. In the case of Visium data in a hexagonal array, we considered three neighbors within the same triangle. For platforms with a square array configuration, we selected four neighbors from the top, bottom, left, and right directions. This approach was implemented using the interpolate function from the SciPy Python package (v1.4.1).
Spatial gaussian process
(Spatial GP) employs a radial basis function (RBF) to model the covariance matrix of spot coordinates48. This method was trained on spot-level coordinates and profiles and then applied to estimate subspot-level profiles based on their spatial coordinates. We utilized the GaussianProcessRegressor from the SciPy Python package (v1.4.1) for implementation. The choice of the RBF kernel weight was determined by testing across a wide range of values for the best performance. The final parameter settings were provided in Supplementary Table 2a.
Image-based enhancement approachesiStar32 is a deep regression model designed to enhance the resolution of ST by minimizing the MSEs between the aggregated inferred gene expressions within a spot and the observed gene expression. For our study, the Python implementation of iStar was used, and the model was run using its default configuration settings.XFuse31 is a deep generative model designed to infer expression maps at enhanced resolution while concurrently reconstructing histological images and spatial gene expressions. We used the Python implementation of XFuse (v0.2.1) and fine-tuned the model under default parameters. Both XFuse and iStar infer profiles at the pixel level. We collected and combined all pixel profiles within the region of a subspot for comparison.Linear Regression Model from Image (Image linear) employs a linear mapping to predict molecular profiles based on image features. For its implementation, we utilized the LinearRegression module from the scikit-learn Python package (v1.2.0). We trained the model at spot-level data and applied it to the subspot image features to predict subspot profiles.Image Gaussian Process (Image GP) utilizes image features as the input for the Gaussian Process model to predict molecular profiles. We implemented Image GP using GaussianProcessRegressor from the SciPy Python package (v1.4.1). The weight of the RBF kernel and white noise used for each spatial omics profile was summarized in Supplementary Table 2b.Image Multilayer Perceptron (Image MLP) models the mapping between image features and molecular profiles through a multilayer perceptron. We implemented the model using MLPRegressor function from scikit-learn Python package (v1.2.0). For each spatial omics dataset, we defined the Image MLP model with the same set of parameters (hidden layer sizes = (128, 32), activation=relu, solver=adam, learning rate init=1E-3, max iter=1E4).Joint enhancement approachJoint Multilayer Perceptron (Joint MLP) model has the same framework as Image MLP but takes in the concatenated feature of the image and spatial position for each spot/subspot for omics profile prediction. Implementations and hyperparameters were following the same setting as in Image MLP.Spatial omics scope (soScope) is implemented in Python with PyTroch (v1.8.0) and PyG (v1.7.2) packages. We provided detailed instructions and a demonstration to run the model on GitHub (details were provided at https://github.com/deng-ai-lab/soScope).EvaluationsRecovery accuracy evaluation with mean square errorMSE measures the difference between two vectors under \({l}_{2}\) norm. An MSE close to 0 indicates an accurate expression recovery; in the implementation, we first min-max normalized the input vectors and then used the mean_squared_error function from sklearn.metrics Python package (v1.2.0).Recovery accuracy evaluation with Pearson correlation coefficientsPearson correlation coefficient is a metric to quantify the consistency between two variables. A coefficient close to 1 indicates a strong positive linear correlation of spatial patterns between reconstructed results and ground truth. In the implementation, we used the pearson function from scipy.stats Python package (v1.10.0).Differential expression measurement with Kolmogorov–Smirnov distanceThe Kolmogorov-Smirnov distance (KS distance) measures the distance between two distributions by calculating the maximum distance between two cumulative distribution functions. A KS distance close to 1 indicates a better separation between two distributions. In the implementation, we used the kstest from scipy.stats Python package (v1.10.0).Gene variability measurementWe used the modelGeneVar function from the scran R package (v1.20.1) to estimate the biological variability of genes. The biological variability is defined as the difference between the log-normalized expression variance and the technical component. A positive biological variability indicates that the observed variability of a gene exceeds the non-informative variation predicted by the model. Conversely, a negative biological component suggests that the observed variability is less than the non-informative variation expected by the model.Gene abundance measurementThe abundance of a gene is the total number of expression counts in a tissue. In the implementation, we used the sum function in numpy from the Python package (v1.19.2).Cluster compactness measurement with average spatial distanceWe first min-max normalized the coordinates of each spot. To assess the compactness of clusters, we computed the average distance between every pair of spots within a cluster. A spatial distance of 0 indicates that all spots within the cluster are closely co-localized in space. For implementation, we used the cdist function from the scipy.spatial.distance Python package (v1.10.0).Software used in this studyThe software used for generating images can be accessed via the following link:Adobe Illustraor: https://www.adobe.com/products/illustrator.htmlSoftware packages used in the study can be accessed via following links:iStar: https://github.com/daviddaiweizhang/istarXFuse: https://github.com/ludvb/xfuseBayesSpace: https://github.com/edward130603/BayesSpaceSeruat: https://satijalab.org/seurat/ArchR: https://www.archrproject.com/sctransform: https://github.com/satijalab/sctransformLinear interpolation: https://docs.scipy.org/doc/scipy/tutorial/interpolate.htmlGaussian Process: https://scikit-learn.org/stable/modules/gaussian_process.htmlLinear Regression Model:https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.htmlMultilayer perceptron: https://scikit.learn.org/stable/modules/neural_networks_supervised.htmlComputational efficiencyThe detailed runtime and memory consumption are reported in Supplementary Table 3. The experiments were conducted on an NVIDlA GeForce RTX 3090 GPU.Statistics and reproducibilityThe experiments were not randomized. No statistical method was used to predetermine the sample size. For the human intestine dataset, we excluded spots that could not be accommodated within the set of hexagonally arranged “low-resolution” spots. For the mouse liver dataset, we excluded spots located outside the inscribed circular area of the slide-DNA-seq and slide-RNA-seq, as our experimental design involved merging the sampling spots to a lower resolution, and the square region facilitated our experimental procedures. For the mouse head, mouse kidney, human ductal carcinoma, mouse embryo (E11), human heart, human skin, and mouse embryo (E13) datasets, no samples were excluded from the analysis. The investigators were not blinded to allocation during experiments and outcome assessment.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.