Accurate imputation of pathway-specific gene expression in spatial transcriptomics with PASTA

Domain: neuroscience
Published: 2025-12-16
Editorial outcome: major_revision

Abstract
Mapping the entire transcriptome at single-cell resolution under its natural spatial context is essential for investigating the oncogenesis and progression of diseases. The recently emerged targeted in-situ technologies retain the spatial organization of cells at high resolution, although they remain limited in the number of genes that can be simultaneously measured. To overcome this obstacle, numerous computational methods have been developed to predict unmeasured gene expression in spatial transcriptomics data by leveraging scRNA-seq data. Most of these methods focus on the expression of individual genes and usually generate highly variable predictions. In this study, we introduce PASTA (PAthway-oriented Spatial gene impuTAtion), a spatial pathway expression imputation method that leverages cell type and spatial proximity to enhance prediction accuracy. PASTA assumes that nearby cells and cells of the same type exhibit similar expression patterns, along with pathway information integrated into the imputation process, which improves prediction robustness and enhances biological relevance in spatial transcriptomics data. We demonstrate PASTA’s superior performance across both simulated and real-world datasets, highlighting its ability to impute pathway gene expression with improved stability and biological significance.

Introduction
Named “Method of the Year 2020” by Nature Methods, spatial transcriptomics (ST) opened the way to quantify gene expression with spatial context1. Unlike traditional RNA sequencing methods, which average gene expression across dissociated cells, ST preserves the spatial relationships between cells, enabling the study of gene expression patterns in situ. This advancement is particularly significant for understanding tissue architecture, cellular heterogeneity, and the molecular micro-environments of complex tissues2–4.

ST is becoming more accessible due to the development of commercialized technologies. Different ST technologies cover spatial content using different methods. For example, Visium and Stereo-seq are sequencing-based, which rely on spatially barcoded DNA; CosMx and MERFISH are imaging-based and rely on in situ hybridization5. However, most of the technologies still suffer from low sensitivity and resolution by limiting their detection to a set of pre-selected genes, which leads to the emergence of computational methods for gene imputation, especially the methods integrating the whole transcriptomic profiling from single-cell RNA-sequencing (scRNA-seq) data. For example, k-nearest neighbors (k-NN) has been used in several existing methods to predict uncovered genes with direct alignment between ST and scRNA-seq data, such as Spatial Gene Enhancement (SpaGE)6 and stPlus7. Another method, Tangram, adopted a different approach through a deep learning framework to map scRNA-seq data directly to ST spots8. This strategy minimizes the distance at both feature and sample levels and improves the imputation accuracy. Along this thread, other existing works also used diffusion model or other deep learning approaches to borrow information from scRNA-seq and predict the gene expression coverage in ST data, including stDiff9, SpotDiff10, and SPRITE11.

The current ST imputation methods share several limitations, such as overestimation of spatial smoothness, uncertainty of the imputed biological accuracy, and reliance on the additional data. None of the methods aims on or is driven by biological pathway information. Pathway analysis has gradually become a standard way to discover and interpret biological functions based on omics datasets. For example, pathway enrichment analysis has been used in understanding the genesis of diseases12,13, detecting druggable targets14,15, and improving the therapeutic effects in patients with insomnia16, myocardial infarction17, and cerebral ischemia18. Many existing tools are available to help researchers move beyond individual gene-level changes to understand the broader biological processes and molecular functions affected by differential gene expression, providing insight into how these genes collectively contribute to disease mechanisms, cellular responses, or physiological changes19. Thus, obtaining pathway expression from ST data can facilitate the translation of spatial findings into clinical and therapeutic applications. Moreover, existing methods generally suffer from high noise level and low accuracy in their imputation results. Combining multiple functionally similar genes from a pathway can improve imputation robustness. Additionally, current imputation methods rarely integrate cell type information, although it has been widely acknowledged that cells belonging to the same cell type have more similar expression levels. Missing this critical knowledge in the imputation process introduces extra noise in the prediction results.

In this paper, we present PAthway-oriented Spatial gene impuTAtion (PASTA), a pathway expression imputation tool for ST data that simultaneously leverages cell type and spatial coordinate information while borrowing information from scRNA-seq data (Fig. 1A). PASTA constructs a mapping matrix between the scRNA-seq reference and ST data by assuming nearby cells and cells from the same cell types share similar expressions. During the mapping, pathway information is provided to guide the training process. In general, gene expression is highly sparse in both scRNA-seq and ST data. In contrast, incorporation of pathway aggregates multiple genes that have more similar behaviors and provides more stable and precise prediction for further analysis. In addition, the pathway-level expression across the ST landscape can offer more biologically related interpretations and conclusions, which reduces the efforts for downstream pathway analysis. PASTA shows its consistent performance in simulated datasets with different pathways and reference data settings. Applied to several real ST datasets, PASTA demonstrates its ability on pathway expression imputation compared with existing methods that serve similar purposes.

Results
Simulation
We carefully performed comprehensive simulations to compare PASTA with several other existing popular ST imputation methods20,21. The simulation datasets were generated based on Stereo-seq datasets using the Axolotl Regenerative Telencephalon22. Compared with other ST technologies, Stereo-seq can measure the whole genome transcriptomic information at single-cell resolution. To begin, we focused on a relatively straightforward scenario: Each simulated ST dataset contains 3000 cells and 1000 genes, which are randomly selected from the Stereo-seq dataset; ScRNA-seq data was obtained from the unselected cells in the same Stereo-seq dataset. This setting ensured that the scRNA-seq and ST datasets maintained similar gene expression distributions. The pathways were constructed by randomly selecting genes from the simulated ST datasets. To summarize pathway expression, we calculated the summation of the expression of the genes in each pathway. In our initial simulations, we opted for a simple but relatively uncommon case where all genes within a pathway were observed in the ST dataset. We then compared the performance of our proposed method against several other established methods, using metrics Mean Squared Error (MSE), Pearson correlation, and Spearman correlation to assess the prediction accuracy. Our results showed that PASTA consistently outperformed the alternatives, achieving the highest correlations and the lowest MSE (Fig. 1B). However, we did observe a slight decline in performance as the number of genes within each pathway increased, suggesting that pathway size can influence imputation accuracy. However, the decline did not continue if we increased the pathway size (Fig. S1). This observation aligns with expectations, as more complex pathways naturally introduce additional variability, making imputation more challenging. Next, we extended our analysis to a more realistic scenario, where only a subset of randomly selected pathway genes were observed in the ST dataset, with the remaining genes coming from the unmeasured portion of the scRNA-seq data. In this partially observed pathway setting, performance across all methods declined, as expected. Nevertheless, our proposed method continued to demonstrate superior performance compared to the other approaches, maintaining the best correlations and lowest MSE (Figs. 1B and S2), no matter whether the genes was observed or not originally (Fig. S3). The ability to handle missing information and still achieve strong results highlights the robustness of PASTA.

We then evaluated the proposed method under some more challenging conditions. In one scenario, we simulated a situation where the reference scRNA-seq data and the ST data were collected from the same subject but different regions (Fig. 1C). The number of genes in a pathway was fixed at 35, and all of them were observed in the simulated ST dataset. This was designed to introduce more variability in cells between the reference and query datasets, as cells from different regions are likely to exhibit distinct expression patterns. Despite this increased complexity, PASTA demonstrated notable stability, consistently providing accurate predictions even in the presence of this variability. In practical applications, collecting data from the same tissue is common, but obtaining both ST and scRNA-seq data from the same subject can be logistically challenging. To address this, we further constructed another scenario where the scRNA-seq data originated from a different subject. After generating the simulated ST data, we down-sampled scRNA-seq data from different individuals’ Stereo-seq datasets. Our proposed method continued to exhibit superior performance, maintaining accuracy and stability despite the increased complexity (Fig. 1D). We also conducted another situation where the scRNA-seq data and the ST data contain cells belong to different cell types. Under this condition, several existing methods, such as Seurat, struggled to produce stable predictions, reflecting their sensitivity to variability between cell types. However, although compared to the other scenario PASTA saw a bit of decline on the performance, it still shows satisfactory prediction precision (Fig. 1E). These results highlight the robustness of PASTA in dealing with the variability inherent in real-world data collection.

We also observed that the performance of our proposed model is not substantially affected by the size of the scRNA-seq reference dataset (Fig. S4), except when the number of reference cells is relatively small (e.g., <2000), in which case increasing the reference size improves imputation accuracy. However, this improvement stops once a sufficient number of reference cells is reached. With respect to the ST data, major cell types tend to exhibit slightly better imputation performance, although the differences are modest (Fig. S5). Since our method borrows information from neighboring cells, a higher number of cells can provide more context and improve the imputation process. Notably, we found that over-clustering of cells introduces unnecessary heterogeneity, which can obscure shared biological signals and reduce the stability of pathway-level imputation (Fig. S6).

Analyzing xenium dataset from human lung cancer tissue
We first collected an ST dataset from human lung cancer tissue with 480 genes processed by 10× Genomics Xenium23. We annotated the dataset manually referring to the PanglaoDB24 for the seven clusters that Xenium provides. The corresponding scRNA-seq dataset was collected from the UCSC cell Browser website25, which is from human lung tissue and contains 26485 genes.

To demonstrate the effectiveness of PASTA in predicting pathway expression, in addition to the methods we used in simulation studies, we expanded our comparison to a few more recently developed tools, including SpaGE, gimVI26, stDiff, stPlus, and stAI27. We examined 44 pathways from the Hallmark gene set collection, which is known to be associated with cancer progression28. We summarized pathway expression by aggregating the expression levels of genes within each pathway and used the overlapping genes between the pathways and the ST data as the ground truth. The predicted pathway expression was visualized by projecting the results into the spatial profiles (Figs. 2 and S7). We did not include stAI in Fig. 2 due to its inability to predict the entire tissue slice; specifically, it requires over 100 GB of GPU memory, which exceeds our available computational resources.

Overall, PASTA significantly achieved higher pathway correlations (Figs. 3B and S8): Compared to SpaGE, predictions of our model showed higher accuracy in 79.5% of the pathways. Compared to Tangram, our model generated better predictions in all pathways. Notably, for certain pathways, such as the Hedgehog signaling pathway involved in regulating cancer stem cells, the correlation reached as high as 0.9. Despite the overall pathway distribution among all cells, within each cell type PASTA shows homogeneous high accuracy in terms of the pathway expression distribution.

We further checked how the predicted pathways connected to different cell types. Figure 3C illustrates the enrichment of hallmark pathways across various cell types, where each dot represents a significant pathway enrichment in a particular cell type. The size of the dot indicates the significance level, with larger dots corresponding to more significant enrichment. PASTA predicted that cancer cell development-related pathways, MYC TARGETS V1 and V2, are enriched in tumor cells. Additionally, PASTA identified immune-related pathways, including INFLAMMATORY RESPONSE, IL2-STAT5 SIGNALING, and ALLOGRAFT REJECTION, as highly expressed in T cells and macrophages, while using only the measured genes failed to capture correct enrichment patterns across cell types. In addition to pathway analysis, we applied PASTA to predict meta-programs (MPs) – gene sets defined by Gavish et al.29 that are upregulated in specific subpopulations of tumor cells. These MPs capture diverse cellular processes underlying intratumor heterogeneity. PASTA predicted that nearly all MPs are highly expressed in basal cells within the tumor area (Fig. S9).

Pathways provide more comprehensive and biologically meaningful information than single genes. We illustrated this by clustering the dataset using pathway expression. We first applied Integrative and Reference-Informed Tissue Segmentation (IRIS)30 to detect the tumor region within the dataset. IRIS integrates scRNA-seq data and biological knowledge to identify domains in the ST dataset, and we set the number of domains as two during the application. Next, we performed K-means clustering on a pathway-by-cell matrix, where each row represents pathway-level expression. For the K-means cluster, we also set the number of clusters to two aiming to distinguish tumor cells from the rest of the dataset. To evaluate clustering performance, we selected three representative regions and compared cell type annotation, pathway-based clustering results, and IRIS-defined clusters (Fig. 3D, from top to bottom). While IRIS grouped all cells within the tumor region into a single cluster, pathway-based clustering achieved a finer separation, effectively distinguishing tumor cells within a tumor-infiltrating lymphocyte environment. In addition, we also observed that pathway expression shows more spatial structure compared to the single gene expression (Fig. S10): It can accurately show the tissue structure, while single gene may fail on capturing it.

Among the two domains detected by IRIS, one mainly contained fibroblasts and immune cells, while the other primarily comprised basal cells. We defined the latter as the “tumor domain." The HALLMARK pathway expression aligned well with the tumor cells distribution (Fig. 3E). We extracted the cells from the tumor domain, and constructed trajectory analysis using the R package Monocle331, which can reflect the sequence of biological expression changes among cells. Instead of using the original gene-by-cell expression matrix, we still used the pathway-by-cell matrix as we used for K-means clustering. The pathway-by-cell matrix reduces the computational complexity. It also removes potential gene noises, which can make the trajectory clearer and more robust, providing a more coherent picture of biological processes over time that can make it easier to interpret shifts in cell states along the trajectory. Following the standard process provided by Monocle3, we generated a principal graph and ordered the cells by their pseudo-development time (pseudotime) using the first principal component as the root to get the trajectory. We then classified the cells into “early stage" and “late stage" based on their pseudotime. The predictions from our proposed method and Tangram showed distinct patterns (Fig. S11): Our method revealed that more developed tumor cells were concentrated in the center, whereas Tangram’s predictions placed them at the right edge. To further validate these results, we examined the expression of the MET gene, which is associated with abnormal cell proliferation. The distribution of the gene expression aligned closely with our method’s predictions, demonstrating that our approach effectively captures the underlying biological development among the tumor cells. K-means cluster results of the other existing methods are presented in Fig. S13.

Analyzing merfish mouse brain ST data
We collected an adult mouse brain dataset processed by Merfish technology (Zhuang-ABCA-3, brain section: Zhuang-ABCA-3.010)32. The dataset contains 82075 cells with a panel of 1122 genes. The cell types were identified by the authors manually, including mostly neurons. The corresponding mouse primary visual cortex scRNA-seq data is collected from the Allen Institute.

We first examined how PASTA and the other spatial feature imputation methods perform on pathway imputation. The ground truth was defined as the summation of the genes’ expression common to both the pathway and the ST dataset. We chose the GOBP pathways related to mouse brain development. Overall, we observed that PASTA achieved the highest correlations between predicted pathway expression and the ground truth (Figs. 4B and S12). Particularly, for the BRAINSTEM DEVELOPMENT pathway, the correlation achieved 0.958 (Fig. 4A).

The trajectory analysis revealed that pathway expressions can be referred for cell development information. We extracted the inhibitory neurons from the data, and then only kept the layer cells for the trajectory analysis. Four groups of cells are shown in Fig. 4C, which are layer 2 cells (L2), layer 2 or 3 cells (L2/3), layer 4 or 5 cells (L4/5), and layer 6 cells (L6). These layers are formed during a specific developmental process known as inside-out neurogenesis, where neurons in deeper layers (L5 and L6) are generated first, followed by neurons in the more superficial layers (L4, L3, and L2)33. We used the predicted GOBP pathways to construct a pathway-by-cell matrix to infer the cells’ development. Using the pathways predicted by PASTA, the estimated pseudotime follows the theoretical inside-out development pattern (Fig. 4D). However, trajectory analysis with pathways predicted by Tangram failed on predicting the time shift. (Fig. 4D).

We then clustered the cells using the pathway expression matrix predicted by different methods. The cluster and cell distributions from PASTA are shown in Fig. 4E, F, along with the cell distributions generated from Tangram and Seurat. We found that pathways predicted by PASTA can generate clear cell distribution in the clusters: Cluster 1 consists of astrocyte and ependymal-like cells, oligodendrocyte progenitor cells, and vascular cells; Cluster 2 consists of the other cells. Astrocytes and oligodendrocytes play a vital role in providing metabolic support to neurons. Both astrocytes and vascular cells are crucial components of the blood-brain barrier. However, clusters generated from pathways predicted by Tangram and Seurat have cell mixture.

Analyzing mouse primary visual cortex data
Another dataset that we examined is collected from mouse primary visual cortex34,35. There are three tissue slices, and the largest contains 119 genes profiled from 82785 cells. We clustered the dataset using Seurat and annotated each cluster manually referring to the PanglaoDB. The corresponding primary visual cortex scRNA-seq data is collected from the Allen Institute, including 15413 cells with the expressions of thousands of genes.

We evaluated the performance of PASTA against other spatial feature imputation methods, including Tangram, SpaGE, gimVI, stDiff, stPlus, and Seurat, to assess its ability to predict pathway expression through checking the correlation between predictions and the ground truth. During the training process, we did not incorporate cell type information for PASTA, as the majority of the cells were neurons that are distributed across the spatial domain. We summarized pathway expression by aggregating the expression levels of pathway genes. We focused on GOBP pathways, particularly those involved in neuron development. The ground truth was defined as the overlapping genes between each pathway and the ST data. Overall, our model achieved the highest correlations in predicting pathway expression, especially for pathways related to brain activity (Fig. 5A, B). We did not include stAI in Fig. 5B due to its inability to predict the entire tissue slice. We observed that some methods, such as SpaGE and Seurat, struggled to predict the pathway expression. In contrast, our model and Tangram showed better performance among all the methods. Our model showed better performance in around 70% of the checked pathways compared to Tangram. We found that our method generated non-satisfactory predictions in some pathways, such as actin-mediated cell contraction and adult behavior pathways, which were expected as the pathways elves did not have a strong relevance with the dataset source.

To check whether the predicted spatial gene expression is biologically meaningful, we first examined how the predicted pathways are enriched in the cells. We found that a variety of neuron-related pathways are highly expressed among the neuron cells, which is accurately detected by the imputation methods. However, some methods, such as Seurat, generated false positive enrichment among the Epithelial cells (Fig. 5C). Additionally, we constructed developmental trajectories using the R package Monocle3 to reveal the progression of cell states. The proposed method and Tangram produced different developmental patterns: our method highlighted early and late stages across the spatial domain, while Tangram predictions were impacted by spatial coordinates, resulting in a uniformly distributed pseudotime in UMAP (Fig. 5D, E). It is rare to have uniformly distributed pseudotime, which might be led by overfit or over-relies on the spatial coordinates during the imputation. If the model prioritizes spatial smoothness, it may impose coordinate-based trends, especially when information from actual gene expression data is limited. In contrast, our method balanced spatial smoothing with biologically informed patterns, capturing the biological characteristics of the dataset in a more meaningful way.

We also applied the imputation methods on the other two additional ST slices. The proposed method still shows advantages in predicting the pathway information (Fig. S14).

Analyzing human frontal cortex ST data
We collected a large dataset that contains 18,686 cell profiles from human frontal cortex tissue. The corresponding scRNA-seq data are collected from a study related to the developmental atlas of the human prefrontal cortex36 with 24153 genes leveled from 2394 cells. The cells in the dataset were identified to 13 clusters, including oligodendrocytes and layer cells. We evaluated the performance of our proposed model against other spatial feature imputation methods: Tangram, SpaGE, gimVI, stDiff, stPlus, and Seurat. The ground truth was defined as the summation of the expression levels of genes common to both the pathways and the ST dataset. With more than 6000 genes present in the spatial dataset, almost all pathway genes were represented. We would like to note that imputation may not be necessary for high-resolution ST datasets, although sparse ST datasets can still benefit from imputation (Fig. S15). We focused on GOBP pathways related to brain development. Our proposed method demonstrated superior performance, in around 70% of the examined pathways, achieving the highest correlations between predicted pathway expression and the ground truth, followed by Tangram and GimVi (Fig. 6A, B). We did not include stAI in Fig. 6B due to its inability to predict the entire tissue slice. The brain renin angiotensin system pathway from GOBP database aligns with the astro cells correctly than highly expressed genes among astro cells (Fig. S16). To further illustrate how the predicted pathways capture biological processes, we analyzed cortical layer cells, specifically L2 cells (outer granular layer), L3 cells (outer pyramidal layer), L4 cells (inner granular layer), and L6 cells (polymorphic layer). The development of the layer cells in the brain follows an inside-out pattern, from L6 to L2. We integrated the predicted expression levels from over 200 GO pathways into a matrix and applied the R package Monocle3 for trajectory construction. By selecting L6 cells as the development root, we ordered the cells according to their pseudotime. The predicted cell development order aligned with biological expectations: L4 cells developed after L6 cells, while L2 and L3 cells developed at the very end (Fig. 6E).

We evaluated the predictions using the recently developed calibration score presented by TISSUE35. The calibration score measures imputation uncertainty and is defined as the ratio between the absolute prediction error and the variability within cell neighborhoods, which integrates both spatial and gene expression similarity. Following the standard procedure outlined by TISSUE, we examined the overall distribution of calibration scores as well as the distribution within each of the 13 cell subsets. Across all cells, our proposed method achieved a calibration score comparable to Tangram and StDiff (Fig. 6C, D). Within individual cell subsets, our method consistently showed lower calibration scores. Overall, the lower median calibration score from PASTA indicates more stable spatial gene expression estimation.

Methods
Method overview
In this section, we will introduce the algorithm of PASTA. We would like to impute the missing pathway expression, i.e., the summation of the gene expression within the pathway, through borrowing gene information in ST data. The gene expression in ST data is generally sparse with high variance. Pathway expression, in contrast, aggregates multiple genes that are more functionally similar. Thus, incorporating the pathway information in the imputation process would stabilize the prediction and provide more accurate estimations. Additionally, pathway information adds known biological knowledge to the results, allowing for easier interpretation and more straightforward downstream analysis.

Let T ∈ RnT×GT denote the ST gene expression matrix measured at nT cells for GT genes (denoted as GT), and S ∈ RnS×GS be the scRNA-seq gene expression matrix measure for the (usually larger) GS genes (denoted as GS) in nS single cells. The two assays share the subset of genes O = GT ∩ GS. We focus on a curated collection of P pathways, for each pathway l its full gene list is Gl ⊆ GS. Within Gl we distinguish the part observed in ST as Gobs l from the remaining part missing in ST but present in scRNA-seq as Gmiss l.

Our goal is to predict at every spatial location, the sum of expression of the entire pathway geneset Gl even when many of its genes are unmeasured in ST. To do so, we first learn a trainable alignment matrix A ∈ RnS×nT, after row-wise soft-max transformation, its entries are converted into Aik = exp(Aij) / ∑n j0=1 exp(Aij0), which form a probability distribution that maps each scRNA-seq cell i onto spatial location j in ST. and S* = ATS is therefore the single cell informed gene expression estimation in ST.

For a given pathway l, we define its imputed expression as Tl = [S*]l 1, where Tl is an nT by 1 pathway activity vector on the tissue slide. Here we are not simply imputing genes one-by-one, we aim to aggregate real ST knowledge and scRNA-seq information into a robust pathway activity measurement.

At each time we focus on one pathway l, the alignment matrix is learned by minimizing a weighted sum of four loss functions that each serves a distinct role. The first reconstruction loss is:
Lrecon = ∑nT j=1 MSE(ATS[:,Gobs l]N(j), T[:,Gobs l]N(j)), (1)
which aligns the pathway expression within each spatial location’s neighboring cells from the same cell type. This is because spatial expression changes smoothly within the same cell type but can jump sharply across different cell types; therefore, for each focal spatial location j, we consider its k-nearest neighboring cells from the same cell type as N(j), k=10 by default.

The global similarity loss is:
Lglobal = ∑g2O Cos sim(S[:,g], T[:,g]), (2)
which enforces transcriptome-wide agreement on all shared genes to enable mapping between scRNA-seq with ST.

The pathway similarity loss is:
Lpathway = ∑nT j=1 Cos sim(S*j,Gobs l, Tj,Gobs l), (3)
which narrows down to the current pathway’s own genes, preventing them from being washed out by highly expressed genes in the global similarity loss, and ensuring the alignment is still optimized for pathways that are under-represented in ST.

The neighborhood similarity loss is:
Lneighbor = ∑nT j=1 Cos sim(S*N(j),Gobs l, TN(j),Gobs l), (4)
which encourages cell-type-coherent expression patterns in each cell’s immediate cellular neighborhoods.

Taking together, the total loss is
L = λ1 Lrecon + λ2 Lglobal + λ3 Lpathway + λ4 Lneighbor. (5)

We minimized the loss function through gradient-based optimization using the PyTorch library. After model training, we obtain the predicted pathway Tl.

As demonstrated by the ablation analysis, all components of the loss function contribute meaningfully to the overall performance of PASTA (Fig. S17). Sensitivity analysis further revealed that the algorithm is particularly responsive to the pathway loss weight: assigning a higher value to this parameter generally leads to improved imputation accuracy. Additionally, performance benefits from a reasonably large learning rate and moderate values for the other loss regularization terms (Fig. S18). By default, the algorithm sets λ1 and λ3 to 1, λ2 to 2, and λ4 to 10.

Our method is inspired by Tangram, a well-established tool for spatial transcriptomics by forcing gene expression in individual cells and in ST spots. However, we improved it by incorporating neighbor cells and cell type information, aiming to produce more accurate and interpretable imputed pathway expression.

Real datasets
Xenium human lung cancer dataset. We collected 161,000 ST cells from human lung cancer tissue processed by 10× Genomics’ high-plex Xenium In Situ technology. The dataset includes 480 genes. Cell annotations were performed manually based on features identified by K-means clustering (K=7) referring to the PanglaoDB (Human, including tumor/cancer samples and cell lines). Annotation details can be found in the next section. The dataset includes 31249 basal cells, 28866 fibroblasts, 27453 T cells, 21975 macrophages, 18946 plasma cells, 18046 epithelial cells, and 13913 endothelial cells. The corresponding reference scRNA-seq data was obtained from the UCSC Cell Browser. This reference dataset includes all cell types observed in the ST dataset, along with 30 additional cell types.

Merfish mouse brain data. The dataset is downloaded through the Allen Brain Cell Atlas. We extracted the cells with brain section number brain section labeled “Zhuang-ABCA-3.010". Cell type annotations were provided by the original authors, and we filtered out cell types with fewer than 40 cells. The final dataset contains 82,075 cells from 27 cell types. The corresponding mouse primary visual cortex scRNA-seq data was collected from the Allen Institute, comprising 15,413 cells with the expression of ten thousands of genes.

ISS mouse primary visula cortex dataset. We analyzed an in situ sequencing (ISS) dataset from the mouse primary visual cortex34,35. This dataset includes three slices: the largest contains 119 genes profiled from 82,785 cells, while the other two contain 136 genes profiled from 14,066 and 15,823 cells, respectively. We clustered the dataset using Seurat and manually annotated each cluster based on the PanglaoDB. Details can be found in the next section. The corresponding scRNA-seq data was collected from the Allen Institute, comprising 15,413 cells with the expression of ten thousands of genes.

CosMx Human frontal cortex dataset. We analyzed a large dataset containing 188,686 cells from human frontal cortex tissue. Using the CosMx Spatial Molecular Imager technique, over 6000 genes were characterized from human frontal cortex FFPE tissue. The dataset mainly contains 7 cell types, which are astrocytes, inhibitory neuron, endothelial, cortical layer cells, microglia, oligodendrocyte, oligodendrocyte progenitor cells (OPCs). The corresponding scRNA-seq data was obtained from a study on the developmental atlas of the human prefrontal cortex36, which includes 24,153 genes from 2394 cells, including six major cell types: neural progenitor cells, excitatory neurons, interneurons, astrocytes, OPCs, and microglia, which are intersected with the cell types in the spatial dataset.

Celltype annotation for the real datasets
Xenium human lung cancer dataset. The clusters and their features were provided by Xenium (https://www.10xgenomics.com/datasets/ffpe-human-lung-cancer-data-with-human-immuno-oncology-profiling-panel-and-custom-add-on-1-standard). We used the K-means cluster with K=7. For each cluster, we ranked features by p-value and log2 fold change, selecting the top five to annotate each cluster using PanglaoDB. The annotation was performed with species set to human, including tumor/cancer samples, cell lines, as well as non-adult and non-primary samples.

ISS mouse primary visual cortex dataset. We annotated this dataset using the Seurat R package, following its standard workflow for normalization and feature selection, with the number of variable features set to 2000. Clustering was performed at a resolution of 0.2, and marker genes for each cluster were identified using the Wilcoxon rank-sum test. Markers were ranked by p-value and log2 fold change, and the top five features were used for cluster annotation, referring to PanglaoDB with species set to mouse, including non-adult and non-primary samples.

Pathway selection
Selecting appropriate pathways is crucial for obtaining accurate and interpretable results. In general, pathway selection should consider dataset characteristics such as source tissue and disease context. We recommend using well-curated pathway databases, such as KEGG and Reactome, and selecting pathways with an appropriate level of detail. Broad pathways may be less informative for datasets focused on a specific biological process. To improve robustness, pathway relevance can be validated using independent datasets.

For the human lung cancer dataset, we used 49 hallmark pathways which has overlapped genes with the ST dataset. For the mouse brain and primary cortex dataset, we used GOBP pathways, particularly, we included the pathways related to brain development. For the human cortex dataset, we checked 40 brain development-related pathways collected from GOBP sets. An example table of the pathways we examined was included in the supplementary file.

Benchmarking
Compared methods. We compared our proposed model with several state-of-art methods, including Tangram, StDiff, gimVI, Seurat, SpaGE, stPlus, and stAI. These methods are all imputation tools for gene expression. After we get the gene expression, we will compute the pathway expression by aggregating the pathway genes for each cell.

Tangram aligns scRNA-seq profiles with ST data to impute the missing genes in it. For different levels of resolution, it voxels the spatial volume differently. For example, a voxel for Visum technology is an individual spots, but for MERFISH is an individual cell. It minimizes an imputation loss function by incorporating gene-by-voxel and voxel-by-gene expressions, as well as the cell density in each voxel. In our implementation, we used the default parameter settings, except for adjusting the voxel-gene regularization term to 2.

StDiff introduces a diffusion model for ST gene imputation referring to scRNA-seq data. Instead of seeking similarity between the scRNA-seq and ST data, the model leverages the rules behind the scRNA-seq data to impute the missing genes in ST data, which can reduce the potential technical noise introduced by the two different experimental platforms. In summary, the model uses two interconnected Markov chains: the forward chain adds noise to the scRNA-seq data. All genes in scRNA-seq data undergoes noise perturbation, then the unique genes in scRNA-seq data will be added a time-dependent noise. Then the shared genes and unique gene part are concatenated to train the model to predict the noise. After the model is trained, unique genes in scRNA-seq data will be set as zero values in ST data and combine with the shared genes’ values in the ST data to be fed into the model to do the prediction. The usage of diffusion model is under-explored in bioinformatics analysis. We expect that the model can capture local dependencies in ST data where gene expression patterns are often highly dependent on location. For real-world datasets, we prepared the inputs of spatial and scRNA-seq data in H5AD format and ran the program with default settings.

GimVI uses a variational inference framework to optimize the joint likelihood of scRNA-seq data and ST data. For scRNA-seq data, it uses zero-inflated negative binomial or negative binomial distribution to sample gene expressions, assumed conditionally independent on a latent variables. For ST data, it uses Poisson or negative binomial distribution to sample the gene expressions, assumed conditionally independent on a latent variable. Then, variational inference approximates the posterior distributions, from which missing genes will be imputed. We trained and ran Gimvi with default parameters.

Seurat, a widely used R package initially designed for scRNA-seq analysis and visualization, has been adapted to impute missing genes in ST datasets. It first jointly reduces the dimension of both scRNA-seq and ST datasets using canonical correlation analysis (CCA) with a L2-normalization on the canonical correlation vectors. Mutual nearest neighbors were used in the low-demensional space to get the “anchors", which refer as the cell pairs between scRNA-seq data and ST data. Then, a weighted matrix is constructed between each query cell (ST dataset) and each anchor to impute the missing genes. We followed the standard Seurat workflow, normalizing the data and using CCA by setting the dimension argument as 30 for anchor detection. Then principal component analysis (PCA) was set to impute the missing genes, as projecting a query dataset onto an existing PCA structure is more efficient in cases working with a large reference dataset and can eliminate the potential batch differences.

SpaGE follows a two-step process: it first applies domain adaptation to align ST and scRNA-seq data in a shared space. Then, for each cell in the ST data, neighbor cells from the scRNA-seq data are defined. Then missing genes’ expressions of that cell are predicted as a weighted average of the expressions of the neighbor cells, where the weight is associated with the cosince distance between cells. We trained the SpaGE model using the ST datasets with their corresponding scRNA-seq reference and imputed the missing genes by setting the number of principal vectors in SpaGE equal to 30.

StPlus jointly embeds ST and scRNA-seq data into a common gene space using an autoencoder. It then predicts the missing gene expression of a spatial cell by computing the weighted average of its neighborhood scRNA-seq cells. We performed pathway gene imputation using its default parameters, setting the random seed to 10.

stAI is a deep learning-based model which extracted sharing features of ST and scRNA-seq data through an encoder-decoder-based backbone framework. Missing genes in ST data were imputed through a weighted linear combination of the genes in the near cells from scRNA-seq data. The method can also perform cell annotation for ST data. We performed pathway gene imputation using its example parameters.

Evaluations. After the gene imputation, we performed the downstream analysis to explore the rationality of the results.

Correlation. To evaluate the performance of pathway imputation, we computed the Pearson correlation between the observed and imputed pathway expression in the ST dataset. For a given pathway, the ground truth was defined as the expression summation of its genes measured in the ST dataset, while the predicted values corresponded to the imputed expression summation of these same genes. By quantifying the correlation between observed and imputed pathway expression, we assessed the extent to which the imputation methods preserved pathway expression patterns within the spatial context.

TISSUE calibration score. The goal of the TISSUE framework35 is to calibrate the predictions of gene expression in ST, incorporating a measure of uncertainty to reflect the confidence in these predictions. The formula of the score is
sij = |Xij - X^ij| / Uij
Uij = 1 + √∑k2NiP Wik (X^kj - X^ij)2 / ∑k2NiP Wik
Wik = exp ||X^i: - X^k:|| / ||X^i: - X^k:|| (6)
The X denotes the gene expression for cell i and gene j, and X^ij is the corresponding predicted expression. Uij is a cell-centric variability that measures the deviation between the predicted expression of a cell and its neighborhood cells, with a weight of the exponential cosine similarity between that cell and its neighbors.

During our implementation, the input of the algorithm was pathway-by-cell matrix. The calibration score could reflect how well the model’s predicted probabilities (and uncertainties) align with the true spatial pathway expression data, which can be a fair measurement for the imputation methods.

Trajectory analysis. We aim to use pathway information to infer the temporal dynamics of biological processes through trajectory analysis. Since pathway expression reflects underlying biological activity, we hypothesize that trajectory analysis based on pathway profiles can effectively capture the progression of biological states in ST data. We perform this analysis using the R package Monocle3, which is designed to model cell state transitions in scRNA-seq data, particularly in contexts such as differentiation, development, and disease progression. In our analysis, the input to Monocle3 is a cell-by-pathway expression matrix, where each row corresponds to a pathway’s expression. We originally included 49 pathways for the Xenium lung cancer dataset and 292 pathways for the CosMx frontal cortex dataset. The number of pathways used can vary depending on tissue complexity and heterogeneity; incorporating a broader and biologically relevant set of pathways generally improves the resolution and robustness of trajectory inference. In practice, including several dozen to around a hundred well-chosen pathways often yield more meaningful and stable trajectories. Following the standard Monocle3 workflow, we preprocess the data and construct developmental trajectories. The root cells of the trajectory are assigned for the theoretically early-developed cells.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript introduces a novel pathway-oriented imputation approach for spatial transcriptomics (ST) data, leveraging\nsingle-cell RNA sequencing (scRNA-seq) as a reference. The motivation is well-articulated: limited gene detection remains\na common technical limitation of in situ spatial transcriptomic platforms. By imputing at the pathway level, the proposed\nmethod provides both biological interpretability and methodological innovation.\nThat said, I have several questions and comments that I hope the authors can address:\nMajor:\n1. Does the number of scRNA-seq cells/genes impact the results of imputation? A brief discussion or experiment showing\nhow varying the size of the reference impacts imputation would be helpful.\n2. The pathways in the simulation are not real pathways. The authors extracted a part of the genes in the ST data as a\n“pathway”. Do the pathways represent known biological modules? Have the authors evaluated performance using real\npathway sets (e.g., from KEGG, Reactome, or MSigDB)?\n3. Related with the above question, the rationale for why incorporating pathway data improves expression imputation\nremains unclear. It is uncertain whether the benefit arises from treating multiple genes collectively or from leveraging the\nfunctional similarity among genes within known pathways. To clarify this, a simulation study comparing the imputation\nperformance using true biological pathway data versus randomly constructed “pathway” groups could provide valuable\ninsights.\n4. Does the cell type impact the results of imputation? For example, does the method perform better on immune cells than\nthe other types of cells?\n5. In the algorithm, there are different regularization parameters corresponding to different loss terms. The authors need to\nprovide some explanations about the parameter tuning. For example, what are the effects of assigning a larger weight to the\npathway loss term? A sensitivity analysis figure showing how performance changes across a range of parameter values\nwould enhance transparency and usability.\n6. Has the author tried how the proposed method performs on very sparse ST data? If the ST data has many dropouts, can\nthe proposed method still impute the pathway expressions accurately?\n7. The ST data of the real studies in the paper are all single cell leveled. Can the proposed methods be applied on spot-level\nST data? If so, the authors should briefly discuss any adaptations needed for this use case.\n8. For the trajectory analysis, the authors used a pathway expression matrix. How does the number of pathways included\naffect the results? For example, is it feasible to use only a few (e.g., 2–5) pathways, or does the method require a larger set\n(e.g., hundreds) to produce meaningful results?\n9. The authors interpreted the biological interpretability from pathway expression using clustering method: The cluster\nresults from PASTA do not have cell type mixtures. However, the authors should explain if the cluster results are meaningful\nbiologically. Were the cells with the same functions grouped together? Or were the cells belong to the same pathways\ngrouped together?\n10. Will the cluster number in ST dataset impact the imputation results?\nMinor:\n1. The letter “B” is missing in Figure 3.", "reviewer_2": "The authors have developed PASTA, a novel spatial pathway expression imputation method that leverages cell type and\nspatial proximity to enhance prediction accuracy. This method is the first to incorporate pathway information into the\nimputation of unmeasured genes. However, there are multiple concerns that need to be addressed:\n1. In practical scenarios, for the task of predicting unmeasured gene expression, we do not have access to the true\nexpression of the predicted genes. Therefore, the genes being predicted should not be used in the training process. The\nauthors need to clarify this point in both the simulated and real-world experiments. It is recommended that the authors\nperform cross-validation experiments at the gene level, similar to the approach used by SpaGE, to address this issue.\n2. Considering that this article incorporates pathway information, the authors could also add cross-validation experiments at\nthe pathway level. For each individual pathway, the following three real-world scenarios should be considered: (1) all genes\nwithin the pathway are measured, (2) some genes within the pathway are measured, and (3) none of the genes within the\npathway are measured. The authors should discuss the pathway prediction performance under these different scenarios.\n3. In the simulation experiments, the authors only mention the variation in the number of genes in each pathway. However,\nthe total number of pathways used is not specified. Does the method differentiate between pathways, or does it integrate all\nthe genes across pathways into a single list for input? These points should be clarified in the manuscript to aid\nunderstanding. Moreover, does the number of pathways used affect the performance of the method? Pathways typically\ncontain dozens to hundreds of genes, yet the article only compares scenarios with 35 and 70 genes. This aspect needs\nfurther elaboration.\n4. In the quantitative evaluation of this study, the authors only assessed performance at the pathway level, i.e., by evaluating\nthe average expression of all genes within a pathway. However, we believe that for the task of imputing unmeasured genes,\nevaluation at the individual gene level is indispensable. The authors should add evaluations of the predictions at the\nindividual gene level, not just for pathway genes, but for all unmeasured genes that are predicted.\n5. In the simulation experiments, the authors should specify the exact proportion of unmeasured genes in the “partial\nobserverd” scenario. How does varying this proportion impact the results?\n6. In the real-world experiments, did the authors predict the genes in pathways that were unmeasured in the ST data? If so,\nthis should be explicitly mentioned in the manuscript, along with a discussion on whether predicting unmeasured pathway\ngenes provided more biological insights compared to the genes already measured in the ST data. If not, we recommend\nadding this analysis.\n7. In situations where nearly all pathway genes are already measured in the ST data (e.g., Section 2.5), such as when\npathway expression is fully captured, does the use of predicted gene expression reveal any biological phenomena that the\nspatially measured gene expression does not? If not, does predicting gene expression in such datasets, where most genes\nare already measured, provide any added value?\n8. The manuscript’s formulas need to be clarified. For example, in the first term of the objective function, the MSE calculation\nis not specified. Based on the definition of 'pathway expression' in Section 4.1, the reader might be misled to think this term\ncalculates the MSE of the average expression of genes within the pathway. However, the implementation in the code differs.\nAdditionally, the use of cell neighborhood information is not represented in the formula, but is only described in the text. The\nformulas should be consistent with the code. For example, the parameters lambda_g1 and lambda_g2 used in the code are\nnot reflected in the formula. Also, the third term in the objective function corresponds to gv_term2 in the code, but this term is\nunrelated to cell type—does the Σ index for k in the formula apply correctly?\n9. Are the MSE term and the cosine similarity term on the same scale? If not, how does the variation of each parameter affect\nthe results? Additionally, the authors should conduct ablation studies to assess the impact of each individual term in the\nobjective function.", "reviewer_3": "In this manuscript, Li et al. proposed a method PASTA for imputation of gene expression in spatial transcriptomics. They\nused gene pathway information to enhance imputation. Here are some comments.\nreferences are missed. Some examples are shown below.\nTISSUE: uncertainty-calibrated prediction of single-cell spatial transcriptomics improves downstream analyses. Nature\nMethods\nstDiff: a diffusion model for imputing spatial transcriptomics through single-cell transcriptomics. Briefings in Bioinformatics\nstAI: a deep learning-based model for missing gene imputation and cell-type annotation of spatial tr anscript omics. Nucleic\nAcids Research\nSpotDiff: Spatial Gene Expression Imputation Diffusion with Single-Cell RNA Sequencing Data Integration. AAAI 2025\ninformation provided in the manuscript, it is hard to evaluate PASTA.\nabove.\nvalues on results should be shown\nBased on the comments above, I do not think this manuscript is suitable for publication in Nature Communications."}, "author_response": "We sincerely thank the associate editor and reviewers for their thorough reviews, thoughtful\ncomments, and constructive suggestions for our manuscript “Accurate Imputation of Pathway-\nspecific Gene Expression in Spatial Transcriptomics with PASTA”. Based on the reviewers’\ncomments, we have carefully revised the manuscript (modifications in the manuscript are\nhighlighted in red font).\nIn this revision, we made the following major improvements:\n• Performed additional simulation experiments to explore news scenarios, including varying\nreference cell and gene numbers, varying pathway sizes, dropout sparsity, and different\nproportions of observed genes in the spatial transcriptomics (ST) data.\n• We included an ablation study to examine the sensitivity of PASTA to key hyperparameters.\n• We provided more explanations/discussions in our real data applications, such as the use\nof pathway data in trajectory inference. We also revised our Method section (page13-14)\nfor clarity, and revised captions for Figs 2-6, to provide better explanations and\nclarifications.\n• We improved the public python package, we now provide a step-by-step vignette that walks\nusers through installation, data import, execution of all core functions, and interpretation\nof outputs.\nA detailed point-to-point response to the reviewers’ comments is attached below. The reviewers’\ncomments are in bold black, and our replies are in blue.\nResponse to Reviewer 1’s comments:\nThis manuscript introduces a novel pathway-oriented imputation approach for spatial\ntranscriptomics (ST) data, leveraging single-cell RNA sequencing (scRNA-seq) as a reference.\nThe motivation is well-articulated: limited gene detection remains a common technical\nlimitation of in situ spatial transcriptomic platforms. By imputing at the pathway level, the\nproposed method provides both biological interpretability and methodological innovation.\nThat said, I have several questions and comments that I hope the authors can address:\nThank you for your positive feedback and constructive comments. Our detailed responses are listed\nbelow.\nMajor:\n1. Does the number of scRNA-seq cells/genes impact the results of imputation? A brief\ndiscussion or experiment showing how varying the size of the reference impacts imputation\nwould be helpful.\nThanks for the comments. Following your comment, we have now added new experiments to\nexplore how scRNA-seq cells/genes impact the results (figures below) through simulation using\nsimulated data generated from Stereo-seq Axolotl Regenerative Telencephalon datasets.\nSpecifically, we evaluated the Pearson’s correlation, Spearman’s correlation and MSE between the\npredicted and the true pathway expressions. We found that increasing the number of cells in\nscRNA-seq reference data improves imputation performance (Fig R1.A), especially when the\nnumber is relatively small (e.g., <2000 cells). We found the median Pearson’s correlation increased\nfrom 0.87 to 0.90 when the cell number increased from 1000 to 3000. However, the improvement\nstops once the reference cell number reaches a sufficient size. In contrast, varying the number of\ngenes in the reference data did not have a noticeable impact on imputation results, the Pearson’s\ncorrelation ranges from 0.71 to 0.94 (sd=0.04) (Fig R1.B). We added a brief discussion of these\nfindings in the manuscript to clarify how the size of the reference influences performance. The\nnew results are shown in Figure S3 and are described in the updated Simulation section (line 110\non pages 8) in the updated manuscript.\nFigure R1. Evaluation of PASTA’s performance with different numbers of genes and cells in\nthe single cell reference dataset. A. The Pearson’s correlation and Spearman’s correlation\nbetween the true and predicted pathway expressions and mean squared error (MSE) of them at\ndifferent number of cells in the reference data. B. Pearson’s correlation, Spearman’s correlation,\nand MSE at different number of genes in the reference data.\n2. The pathways in the simulation are not real pathways. The authors extracted a part of the\ngenes in the ST data as a “pathway”. Do the pathways represent known biological modules?\nHave the authors evaluated performance using real pathway sets (e.g., from KEGG,\nReactome, or MSigDB)? 3. Related with the above question, the rationale for why\nincorporating pathway data improves expression imputation remains unclear. It is uncertain\nwhether the benefit arises from treating multiple genes collectively or from leveraging the\nfunctional similarity among genes within known pathways. To clarify this, a simulation study\ncomparing the imputation performance using true biological pathway data versus randomly\nconstructed “pathway” groups could provide valuable insights.\nWe appreciate the reviewer’s thoughtful suggestion regarding the underlying rationale for\nincorporating pathway data into expression imputation. Our method leverages gene sets (i.e.,\npathways) to aggregate expression signals. Based on our benchmarking, assembling genes into\nsets that reflect the tissue biology, whether or not they coincide with canonical pathways,\nconsistently enhances imputation accuracy.\nTo further investigate this point, we conducted a simulation study using 30 random pathways from\nGOBP (Gene Ontology Biological Process) database, based on the COSMX frontal cortex dataset.\nPASTA achieves a higher Pearson’s correlation (median value = 0.87) among real pathways\ncompared to the simulated pathways generated by randomly selecting genes from the dataset.\nThese results suggest that the pathway imputation improvement not only from grouping genes\ntogether, but more importantly from the biological coherence of those gene sets with respect to the\nbiological context.\nFigure R2. Performance comparison of PASTA of simulated pathways experiments and real\nbiological pathways from GOBP database. A simulation study using real/simulated pathways\nfrom human GOBP (Gene Ontology Biological Process) database. Simulated ST data and their\ncorresponding scRNA-seq data was generated based on the COSMX frontal cortex dataset.\nSimulated pathways are generated by randomly selected genes from the ST data.\n4. Does the cell type impact the results of imputation? For example, does the method perform\nbetter on immune cells than the other types of cells?\nThis is an interesting question. Following your comment, we investigated whether imputation\nperformance varies by cell type. We extracted nine distinct cell types from the Merfish brain\ndataset, including: Astro-Epen (astrocyte and ependymal-like), CB (cerebellum) GABA (Gamma-\naminobutyric acid;), CNU-LGE (cerebral nuclei - lateral ganglionic eminence) GABA, Immune,\nOB-IMN (olfactory bulb - immature neurons) GABA, OEC (lfactory ensheathing cell), OPC-\nOligo (oligodendrocyte progenitor cell), and Vascular. We then conducted additional simulations\nby generating paired scRNA-seq and ST data for each cell type (simulation details can be found in\nSimulation section in the manuscript, line 70). Among them, OEC and Vascular cells showed\nsignificantly lower and higher imputation accuracy, respectively (Wilcoxon test, P < 0.01), while\nthe remaining cell types showed no substantial differences. We believe this variation is related to\nthe number of cells available for each type: OEC had the fewest cells, while Vascular and OPC-\nOligo had the most. Since our method borrows information from neighboring cells, a higher\nnumber of cells can provide more context and improve the imputation process. We have added\nthese observations and the corresponding figure to the revised manuscript. (page 8, line 114; Figure\nS4).\nFigure R3. Evaluation of imputation performance across five cell types in 30 simulation\ndatasets. Nine distinct cell types from the Merfish brain dataset, including: Astro-Epen (astrocyte\nand ependymal-like), CB (cerebellum) GABA (Gamma-aminobutyric acid;), CNU-LGE (cerebral\nnuclei - lateral ganglionic eminence) GABA, Immune, OB-IMN (olfactory bulb - immature\nneurons) GABA, OEC (olfactory ensheathing cell), OPC-Oligo (oligodendrocyte progenitor cell),\nand Vascular. Among them, OEC and Vascular cells showed significantly lower and higher\nimputation accuracy, respectively (Wilcoxon test, P < 0.01), while the remaining cell types showed\nno substantial differences based on Wilcoxon test.\n5. In the algorithm, there are different regularization parameters corresponding to different\nloss terms. The authors need to provide some explanations about the parameter tuning.\nFor example, what are the effects of assigning a larger weight to the pathway loss term?\nA sensitivity analysis figure showing how performance changes across a range of parameter\nvalues would enhance transparency and usability.\nThank you for your suggestion. Following your comments, we have conducted a comprehensive\nsensitivity study to evaluate the sensitivity of our method to key hyperparameters. In one set of\nexperiments, we fixed the learning rate at 0.05 and varied the parameters for pathway loss from\n0.1 to 5, pathway regulation from 0.1 to 5, and neighborhood loss at 1, 2, and 5 (Figure R4).\nAcross all panels (A–E), we observed that pathway prediction accuracy, measured by Pearson’s\ncorrelation, was highly sensitive to the choice of pathway loss. Specifically, pathway loss of 2 and\n5 consistently resulted in higher accuracy (typically > 0.8), while a low pathway loss of 0.01 led\nto poor performance across all neighborhood and reconstruction regulation values (Pearson\ncorrelation around 0.6). When reconstruction loss parameter was fixed at moderate values (e.g., 2),\nPASTA achieved robust performance across the range of neighborhood loss parameter values. For\nexample, in panel C (pathway loss = 1), the median Pearson correlation was consistently around\n0.82. Neighborhood loss parameter between 1 and 2 yielded the best results in most cases,\nsuggesting that moderate weighting of the pathway loss term helps balance the influence of\nindividual gene expression and biological structure. A more detailed description of this sensitivity\nanalysis is provided in the Supplementary Materials (Figure S18).\nFigure R4. Sensitivity analysis of model performance across different regulation levels. Each\npanel corresponds to a fixed pathway loss parameter (from left to right: 0.1, 0.5, 1, 2, and 5), while\nthe x-axis represents increasing reconstruction parameter. The y-axis shows the Pearson\ncorrelation between the predicted and true gene expression values. Lines represent different\nneighborhood loss parameter: red = 1, blue = 2, and green = 5.\n6. Has the author tried how the proposed method performs on very sparse ST data? If the\nST data has many dropouts, can the proposed method still impute the pathway expressions\naccurately?\nThanks for the question. Following your comments, we evaluated our method on 30 simulated\nsparse ST datasets by randomly masking gene expression values to mimic dropout events\ncommonly observed in real ST data. Specifically, each simulated ST data has 3000 cells and 1000\ngenes, and contains 95% zero values after masking the gene expression values. Based on Figure 1,\nwe found that Liger does not performs well, so we did not include it in this simulation. Despite\nthe increased sparsity, our method maintained good imputation performance with Pearson’s\ncorrelation ranging from 0.79 to 0.92 (sd=0.03) (Figure R5). This robustness is likely because our\nmethod aggregates gene expression, which helps buffer against noise and missing values in\nindividual genes. By leveraging shared information across related genes and incorporating the\nspatial structure our method can still recover meaningful pathway activity even under sparse input\nconditions, highlighting the benefits of incorporating pathway knowledge into the imputation\nprocess. The new results are shown in Figure S15 and are described in the updated Result section\n(lines 255 on page 12).\nFigure R5. Evaluation of imputation performance under sparse ST dataset. Boxplots display\nthe distribution of correlation (left) and mean squared error (MSE) (right) for each method. Our\nmethod shows better performance than competing methods (Seurat and Tangram) under the sparse\nsettings. Each dot summarizes the result from a simulation run.\n7. The ST data of the real studies in the paper are all single cell leveled. Can the proposed\nmethods be applied on spot-level ST data? If so, the authors should briefly discuss any\nadaptations needed for this use case.\nThank you for the thoughtful comment. Yes, the proposed methods can potentially be applied to\nspot-level spatial transcriptomics (ST) data. However, spot-level data often capture transcriptomes\nfrom multiple neighboring cells, resulting in a mixture of cell types within each spot. This cellular\nheterogeneity can obscure true cell-type-specific expression patterns and introduce challenges for\ndownstream analyses such as imputation or trajectory inference. Thus, focusing on high resolution\ndata at single cell level is essential for reliable results. To apply our method effectively in this\ncontext, additional preprocessing steps, such as cell type deconvolution or spot refinement could\nbe necessary to estimate the underlying cellular compositions and improve resolution.\nTherefore, following your comment, we have also included a couple of sentences in the Discussion\nsection to emphasize that future method development for robust de-noise tools that specifically\ntargeted for larger/mixed spatial transcriptomics spots. Pairing PASTA with these future methods\nwill likely have added benefits and improved accuracy for enhancing the biological interpretability\nof the imputed pathway signals (page 12). Moreover, we would like to note that in general spot\nlevel data has majority of the genes measured. In such condition, imputation might not be a\nnecessary step as the benefits would be naturally limited.\n8. For the trajectory analysis, the authors used a pathway expression matrix. How does the\nnumber of pathways included affect the results? For example, is it feasible to use only a few\n(e.g., 2–5) pathways, or does the method require a larger set (e.g., hundreds) to produce\nmeaningful results?\nThank you for the insightful question. In our study, we originally used 49 pathways for the Xenium\nlung cancer dataset and 292 pathways for the CosMx frontal cortex dataset. The number of\npathways included can vary depending on the biological complexity and heterogeneity of the tissue.\nIn general, incorporating a larger and biologically relevant set of pathways helps to better capture\nthe diversity of cellular states and transitions, leading to a more comprehensive and stable\ntrajectory. While it is technically feasible to perform trajectory inference using only a few\npathways (e.g., 2–5), doing so may overlook important biological processes. In practice, using\nseveral dozen to a few hundred well-chosen pathways is typically sufficient and more effective for\nmeaningful trajectory inference. We suggest users start from pathways that are related to the organ\ntissue. For example, we choose the brain related pathways for the Merfish brain dataset and cancer-\nrelated hallmark pathways for the Xenium lung cancer dataset. Users can search the keywords,\nsuch as brain, in the pathway database to get them. We have added a brief discussion of this point\nto the manuscript to clarify the influence of pathway number on trajectory analysis (page 19, line\n484).\n9. The authors interpreted the biological interpretability from pathway expression using\nclustering method: The cluster results from PASTA do not have cell type mixtures. However,\nthe authors should explain if the cluster results are meaningful biologically. Were the cells\nwith the same functions grouped together? Or were the cells belong to the same pathways\ngrouped together?\nThank you for the thoughtful comment. In our study, we used pathway-level expression profiles\nto cluster cells, aiming to capture functional similarities beyond transcript-level noise. The\nresulting clusters from PASTA are indeed biologically meaningful. For the Merfish dataset,\nastrocyte and ependymal-like cells (Astro-Epen), oligodendrocyte progenitor cells (OPC-Oligo),\nand vascular cells were grouped together. This grouping is supported by biological knowledge:\nastrocytes and oligodendrocytes play a vital role in providing metabolic support to neurons [1].\nBoth astrocytes and vascular cells are crucial components of the blood-brain barrier [2]. In addition,\nAstro-Epen and OPC-Oligo cells both enriched for gliogenesis and myelination pathways [3].\nIn summary, our cluster revealed biologically interpretable groupings of Astro-Epen, OPC-Oligo,\nand vascular cells based on shared cellular functions. The new results are described in the updated\nResult section (line 210 on page 11).\n10. Will the cluster number in ST dataset impact the imputation results?\nThank you for the question. To assess how the number of clusters in the ST dataset affects\nimputation performance, we conducted a targeted simulation study. Specifically, we extracted a\nsingle cell type and manually grouped the cells into 1, 2, and 3 clusters.\nWe then applied our imputation method to each of these scenarios and evaluated the results.\nWe observed that using 1 cluster (i.e., treating all cells of the same type as a homogeneous group)\nresulted in the best imputation performance; Using 2 clusters led to slightly worse performance;\nAnd using 3 clusters resulted in the lowest accuracy (Figure R6 below and Figure S5 in manuscript).\nThese results suggest that over-clustering introduces unnecessary heterogeneity, which may\nobscure shared biological signals and reduce the stability of pathway-based imputation. In addition,\nwe also checked how under-clustering will affect the results. We treated two cell types as one group\nand checked the single cell type’s performance, which we found the performance dropped a bit as\nno correlation higher than 0.9.\nTherefore, choosing the correct or biologically meaningful number of clusters is important and can\npositively impact imputation quality.\nFigure R6. Effect of clustering granularity on imputation performance in a spatial\ntranscriptomics dataset. A single cell type was manually grouped into 1, 2, and 3 clusters. Under-\nclustering was generated by treating two cell types as one. Over-clustering introduces unnecessary\nheterogeneity.\nMinor:\n1. The letter “B” is missing in Figure 3.\nWe appreciate your insights. We updated the figures accordingly.\nReviewer #1 (Remarks on code availability):\nI have reviewed the code, and it is runnable and reproducible.\nResponse to Reviewer 2’s comments:\nThe authors have developed PASTA, a novel spatial pathway expression imputation method\nthat leverages cell type and spatial proximity to enhance prediction accuracy. This method\nis the first to incorporate pathway information into the imputation of unmeasured genes.\nHowever, there are multiple concerns that need to be addressed:\nThank you for your constructive comments. Our detailed responses are listed below.\n1. In practical scenarios, for the task of predicting unmeasured gene expression, we do not\nhave access to the true expression of the predicted genes. Therefore, the genes being predicted\nshould not be used in the training process. The authors need to clarify this point in both the\nsimulated and real-world experiments. It is recommended that the authors perform cross-\nvalidation experiments at the gene level, similar to the approach used by SpaGE, to address\nthis issue.\nWe appreciate the reviewer’s suggestion. Although in all of our original analyses, the model was\nnot trained on the expression values of the genes being predicted, we recognize that this was not\nstated clearly. We have now added explicit sentences in the updated Methods section (page 14-15)\ndescribing how the measured gene panel is partitioned. In response to the reviewer’s cross-\nvalidation recommendation, we performed a new cross-validation experiment to further validate\nour approach. We first performed the leave-one-out simulation. Across 20 spatial simulation\ndatasets, the median Pearson’s correlation between predicted and observed expression for the\nheld‑out genes remained comparable with Tangram under identical masking (Fig. R7.A).\nWe hope to clarify that our method is not designed for gene-level prediction like SpaGE and\nTangram. Instead, the primary objective of our approach is to impute pathway expression, which\nis computed as an aggregated measure (e.g., the summation) of the expression values of all genes\nwithin a given pathway. The focus is thus on recovering biologically meaningful pathway activity,\nrather than on predicting individual gene expression values.\nThus, we adopted another 2‑fold gene‑mask cross-validation strategy within each pathway: in\nevery fold we randomly hide 50% of the genes belonging to the pathway (treated as unmeasured\nin ST), train the model on the remaining pathway genes and shared genes between scRNA-seq and\nST, and pathway score is then computed leveraging predicted values for the masked genes. Masked\ngenes are excluded from all input features and loss terms thus preventing information leakage.\nAcross the 20 simulation runs, the median Pearson’s correlation between predicted and observed\nheld-out pathway expression of PASTA (0.26) remained higher than Tangram (0.18).\nFigure R7. Cross-validation. A. Leave-one-out cross-validation at gene level. Pearson’s\ncorrelation between predicted and observed expression for the held‑out genes were summarized in\nthe boxplot from 20 simulations. B. Within‑pathway cross‑validation. For every pathway, genes\nwere randomly split 50/50 in each fold. The first half was retained as observed in the simulated\nspatial‑transcriptomics (ST) data and used to train the model; the second half was masked and later\nimputed. Each dot represents one simulation run; horizontal bars indicate the median performance\nacross runs.\n2. Considering that this article incorporates pathway information, the authors could also add\ncross-validation experiments at the pathway level. For each individual pathway, the\nfollowing three real-world scenarios should be considered: (1) all genes within the pathway\nare measured, (2) some genes within the pathway are measured, and (3) none of the genes\nwithin the pathway are measured. The authors should discuss the pathway prediction\nperformance under these different scenarios.\nThank you for the helpful suggestion. In our original simulations, we included two realistic\nscenarios: (1) all genes within a pathway are measured and (2) only a subset of genes within a\npathway are measured. To further address the reviewer's point, we have now extended our\nsimulation to include a third scenario in which (3) none of the genes within a pathway are measured.\nWe observed that when all or some of the genes in a pathway are measured, as shown in Figure 1\nin the manuscript, the imputation performance remains high with a median Pearson’s correlation\nof around 0.8. However, when none of the pathway genes are observed, all methods experience a\nsignificant drop of median Pearson’s correlation to around 0.3 in performance (Figure R8 in this\nresponse letter, Figure S2 in manuscript). Despite this, our method still achieves better accuracy\ncompared to the alternatives with Wilcoxon’s test P-value less than 0.01.\nWe note that the scenario where none of the genes in a pathway are measured is likely rare in\npractical settings. This is because spatial transcriptomics platforms typically prioritize the\ninclusion of biologically meaningful genes, and key pathways relevant to the study design are\nusually at least partially captured in the gene panel. For example, when we examine the Xenium\nlung cancer dataset with only 480 genes, among all hallmark pathways, only 5 of them have no\noverlapping genes with the observed data, while the other pathways have all or partial overlap with\nthe observed gene list, suggesting the rarity of pathways with none overlapping genes. We have\nadded this analysis and discussion to the revised manuscript to provide a more complete picture of\nour method’s robustness under various real-world conditions (page 6, line 89).\nFigure R8. Impact of pathway gene coverage on pathway prediction accuracy. Performance\nis shown for three evaluation metrics, Pearson correlation (left), Spearman correlation (middle),\nand mean‑squared error (MSE, right), between imputed and true expression of the masked genes.\nThe x‑axis denotes the proportion of each pathway’s genes that are observed in the ST panel (0,\n1/2, 1/3, and 1/4). Box‑plots represent the distribution across 30 independent simulation runs.\n3. In the simulation experiments, the authors only mention the variation in the number of\ngenes in each pathway. However, the total number of pathways used is not specified. Does\nthe method differentiate between pathways, or does it integrate all the genes across pathways\ninto a single list for input? These points should be clarified in the manuscript to aid\nunderstanding. Moreover, does the number of pathways used affect the performance of the\nmethod?\nThank you for the insightful comment. Our algorithm is pathway‑specific by design, during\ntraining we focus on one predefined pathway at a time, supplied as a simple gene‑list, and optimize\nthe alignment matrix exclusively for that pathway. Genes that belong solely to other pathways\nparticipate only through the transcriptome‑wide similarity term and do not influence the\npathway‑specific losses. For advanced use cases, we also updated our software package to support\ninput of multiple pathways. A user can pass a list of pathways at a time, and the pipeline will then\niterate internally and return an imputed activity vector for each pathway. This multi‑pathway mode\noffers convenience and enables parallel computing, but each pathway is still processed in an\nindependent loop, so inference accuracy is unaffected. We have clarified these points in the revised\nDiscussion section (page 12).\n4. Pathways typically contain dozens to hundreds of genes, yet the article only compares\nscenarios with 35 and 70 genes. This aspect needs further elaboration.\nThank you for the observation. We have extended our simulation experiments to include pathways\ncontaining a larger number of genes, reflecting more realistic biological scenarios in which\npathways often consist of dozens to hundreds of genes. Figure R9.A presents the gene numbers in\nthe pathways we investigated, ranging from 35 to 300 genes per pathway. In Figure R9.B (also\navailable as Figure S1 in manuscript), we observed that as the number of genes per pathway\nincreases, imputation accuracy shows a slight decline. This is likely due to greater heterogeneity\nand increased noise associated with larger gene sets. We clarified this point in the manuscript and\nnote that our framework remains effective for pathways containing dozens to hundreds of genes,\nas commonly observed in real biological datasets (page 6, line 83).\nFigure R9. Evaluation of imputation performance with different size of pathways.\nA. Histogram of the pathway sizes checked in our simulation study. B. Imputation accuracy for\nthe two largest pathway sizes (100‑gene pathways in red and 200‑gene pathways in blue). For each\nmetric, Pearson correlation, Spearman correlation, and mean‑squared error (MSE), box‑plots show\nthe distribution across simulation replicates for the three competing methods (PASTA, Seurat,\nTangram).\n4. In the quantitative evaluation of this study, the authors only assessed performance at the\npathway level, i.e., by evaluating the average expression of all genes within a pathway.\nHowever, we believe that for the task of imputing unmeasured genes, evaluation at the\nindividual gene level is indispensable. The authors should add evaluations of the predictions\nat the individual gene level, not just for pathway genes, but for all unmeasured genes that\nare predicted.\nWe appreciate the reviewer’s insightful suggestion.\nUsing the MERFISH mouse brain dataset, we randomly selected 70 genes from the GOBP (Gene\nOntology Biological Process) pathways and then imputed them with PASTA, Tangram, and\nSeurat. PASTA achieved a median Pearson correlation of 0.93 for single gene prediction, compared\nwith 0.71 for Tangram and 0.3 for Seurat (Spearman and MSE yielded the same ranking; Wilcoxon\nsigned‑rank tests P < 0.001 for all pairwise comparisons, Figure R10).\nNevertheless, we would like to emphasize again that our method does not produce gene-level\npredictions as an output. Its strength lies in providing robust and biologically interpretable\nestimates of pathway activity across spatial locations. The primary output of our model is pathway\nexpression, computed as the aggregate (e.g., summation) of gene expressions within a predefined\npathway. Our focus is on accurately imputing these pathway-level signals, rather than predicting\nindividual gene expression values. We believe that pathways aggregating multiple genes will\nreduce noise compared to the single genes. In addition, pathways expression also provides more\nprobabilities on downstream analysis which can be functionally interpretable. As a result, we\ncurrent still maintain all the pathway-level evaluation but leave out the gene-level evaluation to\navoid further confusion.\nFigure R10. Single gene imputation results in the Merfish mouse brain dataset. A total of 70\ngenes from GOBP pathways were randomly selected and imputed their spatial expressions.\nBox‑plots show the Pearson correlations between imputed and ground‑truth expression across the\n70 genes; points are individual genes.\n5. In the simulation experiments, the authors should specify the exact proportion of\nunmeasured genes in the “partial observed” scenario. How does varying this proportion\nimpact the results?\nThank you for the helpful comment. In our original simulation, the “partially observed” scenario\nincluded 50% unmeasured genes within each pathway (i.e., 35 out of 70 genes were unobserved).\nTo further investigate the impact of this proportion, we expanded our simulations to include\nadditional conditions where 1/3 and 1/4 of the pathway genes were unmeasured.\nWe observed that when all or some of the genes in a pathway are measured, as shown in Figure 1\nin the manuscript, the imputation performance remains high with a median Pearson’s correlation\nof around 0.8. The median of Pearson’s correlation between the imputed and true pathway\nexpressions increased from 0.7 to 0.8 for our method as the partial portion dropped from 1/2 to\n1/4. Our method achieves better accuracy compared to the alternatives with Wilcoxon’s test P-\nvalue less than 0.01. This trend is intuitive, as having more observed genes provides stronger\ncontext for estimating pathway-level activity. We have added these observations and the\ncorresponding results to the manuscript (page 6, line 89).\nFigure R11. Impact of pathway gene coverage on pathway prediction accuracy. Performance\nis shown for three evaluation metrics, Pearson correlation (left), Spearman correlation (middle),\nand mean‑squared error (MSE, right), between imputed and true expression of the pathway genes.\nThe x‑axis denotes the proportion of each pathway’s genes that are observed in the ST panel (1/2,\n1/3, and 1/4). Box‑plots represent the distribution across 30 independent simulation runs.\n6. In the real-world experiments, did the authors predict the genes in pathways that were\nunmeasured in the ST data? If so, this should be explicitly mentioned in the manuscript,\nalong with a discussion on whether predicting unmeasured pathway genes provided more\nbiological insights compared to the genes already measured in the ST data. If not, we\nrecommend adding this analysis.\nThanks for pointing this. Yes, in our real-world experiments, we predicted genes in pathways that\nwere unmeasured in the ST data. All downstream analyses including clustering and differential\nexpression were performed using the complete set of pathways, regardless of whether they were\noriginally measured in the ST data. Although the ST data does not observe the genes originally,\nthe imputation process learns the genes expression from the scRNA-seq reference, which we think\nwill provide us more information regarding the true pathway expression.\nFor example, the Xenium lung cancer ST data measures 289 genes from all the hallmark pathways,\nwhereas our curated Hallmark pathways span almost 4000 genes. Consequently, about 89 % of\npathway genes were unmeasured in ST but were imputed by our model. Specifically, individual\ngene maps often fail to identify tumor from non-tumor structures, while our imputed hallmark\npathway expressions clearly show the separation between tumor and non-tumor structures, which\nis also shown by differential expression analysis results (Figure 2E and Figure S10). Beyond visual\nstratification, pathway expression provides more possibilities on downstream analysis. For\nexample, after imputing the hallmark pathways in the Xenium lung cancer ST data, we performed\npathway enrichment testing across different cell types, generating interpretable functional\nsignatures (Figure 1D), and also reconstructed a plausible tumor evolution trajectory (Figure S11).\nWe revised the manuscript to explicitly clarify this point and included a discussion on how\nincorporating unmeasured pathway genes can enhance biological interpretation by providing a\nmore comprehensive view of pathway activity.\n7. In situations where nearly all pathway genes are already measured in the ST data (e.g.,\nSection 2.5), such as when pathway expression is fully captured, does the use of predicted\ngene expression reveal any biological phenomena that the spatially measured gene expression\ndoes not? If not, does predicting gene expression in such datasets, where most genes\nare already measured, provide any added value?\nThanks for the comments. In the Cosmx frontal cortex dataset (Section 2.5), although most\npathway genes are already measured, using pathway-level imputed values still offers advantages.\nFrom pathway level, our reconstruction loss, ∑𝑛𝑇 𝑀𝑆𝐸((𝐴𝑇𝑆 ) ,(𝑇 ) ) (Method\n𝑗=1 :,𝒢𝑜𝑏𝑠 :,𝒢𝑜𝑏𝑠\nl 𝒩(𝑗) l 𝒩(𝑗)\nsection), aligns the pathway expression within each spatial location’s neighboring cells from the\nsame cell type, reducing possible noise of gene expression measured in ST data across the cell\ntypes and cell locations. For example, predicted pathway expression enables accurate pseudo time\nestimation while reducing the dimensionality of the input expression matrix. Specifically, we\nestimated the pseudo time correctly using around 300 pathways instead of using the large gene\nexpression matrix with more than 6,000 genes while reducing the computation time. In addition,\npathway expression can highlight some brain structures while single gene might have difficulties\nregarding of that. For example, the brain renin angiotensin system pathway from GOBP database\naligns with the astro cells correctly than other highly expressed genes among astro cells (Figure\nS16).\nIn reality, if almost all genes were observed, the added value of gene imputation is naturally limited.\nIn such cases, we would not expect users to apply imputation tools. However, many ST platforms,\nsuch as Xenium, Merfish, Cosmx, are sparse due to technical limitations and targeted panels, where\nimputation tools can still be applied on the datasets to extract additional information. We performed\na new simulation study to mimic the impact of data sparsity on imputation performance using ST\ndata. Specifically, we randomly masking gene expression values to mimic dropout events in the\nsimulated ST data. Despite the increased sparsity, our method maintained good imputation\nperformance, suggesting the possible application under such cases (Figure R12 in this response\nletter, Figure S15 in manuscript).\nFigure R12. Evaluation of imputation performance under sparse ST dataset. Boxplots show\nthe distribution of performance across 30 simulation runs for the three evaluation metrics: Pearson\ncorrelation (left), Spearman correlation (center), and mean‑squared error (MSE) (right) between\nimputed and ground‑truth expression. Sparsity was simulated by retaining only 5 % of the original\nST gene panel. PASTA (red) consistently achieves the highest correlations and the lowest MSE,\nindicating superior robustness under extreme gene‑dropout conditions.\n8. The manuscript’s formulas need to be clarified. For example, in the first term of the\nobjective function, the MSE calculation is not specified. Based on the definition of 'pathway\nexpression' in Section 4.1, the reader might be misled to think this term calculates the MSE\nof the average expression of genes within the pathway. However, the implementation in the\ncode differs. Additionally, the use of cell neighborhood information is not represented in the\nformula, but is only described in the text. The formulas should be consistent with the code.\nFor example, the parameters lambda_g1 and lambda_g2 used in the code are\nnot reflected in the formula. Also, the third term in the objective function corresponds to gv_term2 in the code, but this term is\nunrelated to cell type—does the Σ index for k in the formula apply correctly?\nThanks for the insights regarding our formula and code. We have rewritten the Method section\n(page 13-14).\nWe redefine the ST and scRNA-seq matrix: Let 𝑇 ∈ ℝ𝑛𝑇×𝐺𝑇 be the ST gene expression matrix\nmeasured at 𝑛 spatial locations for 𝐺 genes (denoted as 𝒢 ), and 𝑆 ∈ ℝ𝑛𝑠×𝐺𝑆 be a scRNA-seq\n𝑇 𝑇 𝑇\nmatrix measured for the (usually larger) 𝐺 genes (denoted as 𝒢 ) in 𝑛 single cells.\n𝑆 𝑆 𝑆\nWe assume the two assays share the subset of genes 𝒪 = 𝒢 ∩𝒢 . We focus on a curated collection\n𝑇 𝑆\nof 𝑃 pathways, for each pathway 𝑙 its full gene list is 𝒢 ⊆ 𝒢 ). Within 𝒢 we distinguish the part\nl 𝑆 l\nobserved in ST as 𝒢𝑜𝑏𝑠 from the remaining part missing in ST but present in scRNA-seq as 𝒢𝑚𝑖𝑠𝑠.\nl l\nWe first learn a trainable alignment matrix 𝐴̃ ∈ ℝ𝑛𝑆×𝑛𝑇, after row-wise soft-max transformation,\nexp(𝐴̃ )\nits entries are converted into 𝐴 = 𝑖𝑗 , which form a probability distribution that maps each\n𝑖𝑗 ∑ 𝑛𝑇 𝐴̃\n𝑗′=1 𝑖𝑗′\nscRNA‑seq cell 𝑖 onto spatial location 𝑗 in ST, and 𝑆∗ = 𝐴𝑇𝑆 is therefore the single cell informed\ngene expression estimation in ST.\nFor a given pathway 𝑙, we define its imputed expression as 𝑇 = [𝑆∗ ]𝟏, where 𝑇 is a 𝑛 by 1\nl :,𝒢 l 𝑇\nl\npathway activity vector on the tissue slide. At each time we focus on one pathway 𝑙, the alignment\nmatrix is learned by minimizing a weighted sum of four loss functions that each serves a distinct\nrole. The first reconstruction loss is:\n𝑛𝑇\n𝐿 = ∑ 𝑀𝑆𝐸((𝐴𝑇𝑆 ) ,(𝑇 ) )\n𝑟𝑒𝑐𝑜𝑛 :,𝒢𝑜𝑏𝑠 :,𝒢𝑜𝑏𝑠\n𝑗=1 l 𝒩(𝑗) l 𝒩(𝑗)\n, which aligns the pathway expression within each spatial location’s neighboring cells from the\nsame cell type. This is because spatial expression changes smoothly within the same cell type but\ncan jump sharply across different cell types, therefor for each focal spatial location 𝑗, we consider\nits k-nearest neighboring cells from the same cell type as 𝒩(𝑗), k= 10 by default.\nThe global similarity loss is:\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑔𝑙𝑜𝑏𝑎𝑙 𝑠𝑖𝑚 :,𝑔 :,𝑔\n𝑔∈𝒪\n, which enforces transcriptome-wide agreement on all shared genes to enable mapping between\nscRNA-seq with ST.\nThe pathway similarity loss is:\n𝑛𝑇\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑝𝑎𝑡ℎ𝑤𝑎𝑦 𝑗=1 𝑠𝑖𝑚 𝑗,𝒢 l 𝑜𝑏𝑠 𝑗,𝒢 l 𝑜𝑏𝑠\n, which narrows down to the current pathway’s own genes, preventing them from being washed\nout by highly expressed genes in the global similarity loss, and ensuring the alignment is still\noptimized for pathways that are under‑represented in ST.\nThe neighborhood similarity loss is:\n𝑛𝑇\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑟\n𝑗=1\n𝑠𝑖𝑚 𝒩(𝑗),𝒢\nl\n𝑜𝑏𝑠 𝒩(𝑗),𝒢\nl\n𝑜𝑏𝑠\n, which encourages cell‑type‑coherent expression patterns in each cell’s immediate cellular\nneighborhoods.\nTaking together, the total loss is\n𝐿 = 𝜆 𝐿 +𝜆 𝐿 +𝜆 𝐿 +𝜆 𝐿\n1 𝑟𝑒𝑐𝑜𝑛 2 𝑔𝑙𝑜𝑏𝑎𝑙 3 𝑝𝑎𝑡ℎ𝑤𝑎𝑦 4 𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑟\nWe have unified the symbols between the code and the formulas in the manuscript.\n9. Are the MSE term and the cosine similarity term on the same scale? If not, how does the\nvariation of each parameter affect the results? Additionally, the authors should conduct\nablation studies to assess the impact of each individual term in the objective function.\nThank you for your suggestion. Following your comments, we have conducted a comprehensive\nsensitivity study to evaluate the sensitivity of our method to key hyperparameters. In one set of\nexperiments, we fixed the learning rate at 0.05 and varied the parameters for pathway loss from\n0.1 to 5, pathway regulation from 0.1 to 5, and neighborhood loss at 1, 2, and 5 (Figure R13).\nAcross all panels (A–E), we observed that pathway prediction accuracy, measured by Pearson’s\ncorrelation, was highly sensitive to the choice of pathway loss. Specifically, pathway loss of 2 and\n5 consistently resulted in higher accuracy (typically > 0.8), while a low pathway loss of 0.01 led\nto poor performance across all neighborhood and reconstruction regulation values (Pearson\ncorrelation around 0.6). When reconstruction loss parameter was fixed at moderate values (e.g., 2),\nPASTA achieved robust performance across the range of neighborhood loss parameter values. For\nexample, in panel C (pathway loss = 1), the median Pearson correlation was consistently around\n0.82. Neighborhood loss parameter between 1 and 2 yielded the best results in most cases,\nsuggesting that moderate weighting of the pathway loss term helps balance the influence of\nindividual gene expression and biological structure. A more detailed description of this sensitivity\nanalysis is provided in the Supplementary Materials (Figure S18).\nFigure R13. Sensitivity analysis of model performance across different regulation levels. Each\npanel corresponds to a fixed pathway loss parameter (from left to right: 0.1, 0.5, 1, 2, and 5), while\nthe x-axis represents increasing reconstruction parameter. The y-axis shows the Pearson\ncorrelation between the predicted and true gene expression values. Lines represent different\nneighborhood loss parameter: red = 1, blue = 2, and green = 5.\nReviewer #2 (Remarks on code availability):\nThe provided code appears to cover the main components of the proposed method; however,\nit would significantly enhance the reproducibility and utility of the resource if the authors\ncould include detailed code for all experiments reported in the paper. In particular, I\nrecommend adding:\n(1) Code for simulation experiments, including data generation and evaluation metrics.\n(2) Scripts used to reproduce the real-data experiments presented in the main paper and\nsupplementary materials.\n(3) A clear structure separating different experimental settings, with corresponding\ninstructions for running each.\n(4) Additionally, while a README file is provided, it would be helpful to include more\ndetailed instructions on installing dependencies, setting up the environment, and running\neach experiment end-to-end. These improvements would greatly benefit readers who wish to\nreplicate the results or build upon this work.\nWe appreciate the valuable suggestions. As suggested, we built a new website and have now added\ncode for the simulation experiments, scripts used to reproduce the real data experiments, and a\nclear structure separating different experimental setting. We have also included more detailed\ninstructions to our README file regarding the installing dependencies, setting up the environment,\nand running each experiment end to end. We agree with this reviewer that these additions greatly\nimprove the user experience who would like to apply our method.\nFigure R14. Screen-shot of our newly updated website. We added code for the simulation\nexperiments and scripts used for the real data experiments.\nResponse to Reviewer 3’s comments:\nIn this manuscript, Li et al. proposed a method PASTA for imputation of gene expression in\nspatial transcriptomics. They used gene pathway information to enhance imputation. Here\nare some comments.\nThank you for your constructive comments. Our detailed responses are listed below.\ntranscriptomics were cited. Too many references are missed. Some examples are shown below.\nTISSUE: uncertainty-calibrated prediction of single-cell spatial transcriptomics improves\ndownstream analyses. Nature Methods\nstDiff: a diffusion model for imputing spatial transcriptomics through single-cell\ntranscriptomics. Briefings in Bioinformatics\nstAI: a deep learning-based model for missing gene imputation and cell-type annotation of\nspatial tr anscript omics. Nucleic Acids Research\nSpotDiff: Spatial Gene Expression Imputation Diffusion with Single-Cell RNA Sequencing\nData Integration. AAAI 2025\nThanks for your insights. We have checked the citation in the manuscript accordingly and added\nall the missing ones.\nin PASTA? Using the method information provided in the manuscript, it is hard to evaluate\nPASTA.\nThanks for the comments. We rewrite the method section to provide more systematic information\nfor our method. In summary, the core idea of PASTA is to align ST cells to scRNA-seq cells via\na soft alignment matrix, then use this mapping to impute spatial pathway expression. This is\nachieved by learning a probabilistic mapping between spatial cells and scRNA-seq reference cells,\nusing a multi-term loss function that incorporates pathway consistency, local neighborhood\nagreement, and gene-level similarity. The details are as below.\nLet 𝑇 ∈ ℝ𝑛𝑇×𝐺𝑇 be the ST gene expression matrix measured at 𝑛 spatial locations for 𝐺 genes\n𝑇 𝑇\n(denoted as 𝒢 ), and 𝑆 ∈ ℝ𝑛𝑠×𝐺𝑆 be a scRNA-seq matrix measured for the (usually larger) 𝐺\n𝑇 𝑆\ngenes (denoted as 𝒢 ) in 𝑛 single cells. The two assays share the subset of genes 𝒪 = 𝒢 ∩𝒢 .\n𝑆 𝑆 𝑇 𝑆\nWe focus on a curated collection of 𝑃 pathways, for each pathway 𝑙 its full gene list is 𝒢 ⊆ 𝒢 ).\nl 𝑆\nWithin 𝒢 we distinguish the part observed in ST as 𝒢𝑜𝑏𝑠 from the remaining part missing in ST\nl l\nbut present in scRNA-seq as 𝒢𝑚𝑖𝑠𝑠.\nl\nOur goal is to predict at every spatial location, the sum of expression of the entire pathway gene\nset 𝒢 even when many of its genes are unmeasured in ST. To do so, we first learn a trainable\nl\nalignment matrix 𝐴̃ ∈ ℝ𝑛𝑆×𝑛𝑇, after row-wise soft-max transformation, its entries are converted\nexp(𝐴̃ )\ninto 𝐴 = 𝑖𝑗 , which form a probability distribution that maps each scRNA‑seq cell 𝑖 onto\n𝑖𝑗 ∑ 𝑛𝑇 𝐴̃\n𝑗′=1 𝑖𝑗′\nspatial location 𝑗 in ST, and 𝑆∗ = 𝐴𝑇𝑆 is therefore the single cell informed gene expression\nestimation in ST.\nFor a given pathway 𝑙, we define its imputed expression as 𝑇 = [𝑆∗ ]𝟏, where 𝑇 is a 𝑛 by 1\nl :,𝒢 l 𝑇\nl\npathway activity vector on the tissue slide. Here we are not simply imputing genes one‑by‑one, we\naim to aggregate real ST knowledge and scRNA-seq information into a robust pathway activity\nmeasurement.\nAt each time we focus on one pathway 𝑙, the alignment matrix is learned by minimizing a weighted\nsum of four loss functions that each serves a distinct role. The first reconstruction loss is:\n𝑛𝑇\n𝐿 = ∑ 𝑀𝑆𝐸((𝐴𝑇𝑆 ) ,(𝑇 ) )\n𝑟𝑒𝑐𝑜𝑛 :,𝒢𝑜𝑏𝑠 :,𝒢𝑜𝑏𝑠\n𝑗=1 l 𝒩(𝑗) l 𝒩(𝑗)\n, which aligns the pathway expression within each spatial location’s neighboring cells from the\nsame cell type. This is because spatial expression changes smoothly within the same cell type but\ncan jump sharply across different cell types, therefor for each focal spatial location 𝑗, we consider\nits k-nearest neighboring cells from the same cell type as 𝒩(𝑗), k= 10 by default.\nThe global similarity loss is:\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑔𝑙𝑜𝑏𝑎𝑙 𝑠𝑖𝑚 :,𝑔 :,𝑔\n𝑔∈𝒪\n, which enforces transcriptome-wide agreement on all shared genes to enable mapping between\nscRNA-seq with ST.\nThe pathway similarity loss is:\n𝑛𝑇\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑝𝑎𝑡ℎ𝑤𝑎𝑦 𝑗=1 𝑠𝑖𝑚 𝑗,𝒢 l 𝑜𝑏𝑠 𝑗,𝒢 l 𝑜𝑏𝑠\n, which narrows down to the current pathway’s own genes, preventing them from being washed\nout by highly expressed genes in the global similarity loss, and ensuring the alignment is still\noptimized for pathways that are under‑represented in ST.\nThe neighborhood similarity loss is:\n𝑛𝑇\n𝐿 = −∑ 𝐶𝑜𝑠 (𝑆∗ ,𝑇 )\n𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑟\n𝑗=1\n𝑠𝑖𝑚 𝒩(𝑗),𝒢\nl\n𝑜𝑏𝑠 𝒩(𝑗),𝒢\nl\n𝑜𝑏𝑠\n, which encourages cell‑type‑coherent expression patterns in each cell’s immediate cellular\nneighborhoods.\nTaking together, the total loss is\n𝐿 = 𝜆 𝐿 +𝜆 𝐿 +𝜆 𝐿 +𝜆 𝐿\n1 𝑟𝑒𝑐𝑜𝑛 2 𝑔𝑙𝑜𝑏𝑎𝑙 3 𝑝𝑎𝑡ℎ𝑤𝑎𝑦 4 𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑟\nWe minimized the loss function through gradient-based optimization using the PyTorch library.\nWe hope this can help readers to better understand our method.\nlist of methods can be found above.\nThank you for your comments. We had already included stDiff in our original method evaluation,\npreviously, following your suggestion, we have now added a comparison with stAI as well. All\nmethods are evaluated across four different platforms that covers probe-based, in-situ sequencing,\nand imaging-based technologies (Xenium, Merfish, ISS, and Cosmx). The implementation details\nand hyper-parameter choices are provided in the updated Methods section (pages xxx, lines xxx-\nxxx). Across datasets, PASTA consistently achieved the highest accuracy for pathway expression\nprediction. For example, in the Xenium lung cancer dataset containing 480 genes, we calculated\nthe Pearson’s correlation among the predicted and the true pathway expression. PASTA achieved\nthe highest correlation value (median = 0.7), followed by SpaGE (median = 0.48) and stAI (median\n= 0.43). In the ISS frontal cortex data, PASTA also achieves the highest correlation score (median\n= 0.80) compared with the other methods (stAI median = 0.49). We attempted to include SpotDiff,\nbut no public implementation (GitHub, PyPI, or Docker image) was available at the time of\nrevision, we have nevertheless cited the paper in the Introduction section. These results suggest\nthat the pathway expression imputed by PASTA more accurately reflects the underlying biological\nsignal in the tissue, at least relative to the other evaluated methods. The new results are shown in\nthe updated Figures 3-6 and are described in detail in the updated Results section accordingly.\nof the choice of parameter values on results should be shown\nThank you for your suggestion. Following your comments, we have conducted a comprehensive\nsensitivity study to evaluate the sensitivity of our method to key hyperparameters. In one set of\nexperiments, we fixed the learning rate at 0.05 and varied the parameters for pathway loss from\n0.1 to 5, pathway regulation from 0.1 to 5, and neighborhood loss at 1, 2, and 5 (Figure R15).\nAcross all panels (A–E), we observed that pathway prediction accuracy, measured by Pearson’s\ncorrelation, was highly sensitive to the choice of pathway loss. Specifically, pathway loss of 2 and\n5 consistently resulted in higher accuracy (typically > 0.8), while a low pathway loss of 0.01 led\nto poor performance across all neighborhood and reconstruction regulation values (Pearson\ncorrelation around 0.6). When reconstruction loss parameter was fixed at moderate values (e.g., 2),\nPASTA achieved robust performance across the range of neighborhood loss parameter values. For\nexample, in panel C (pathway loss = 1), the median Pearson correlation was consistently around\n0.82. Neighborhood loss parameter between 1 and 2 yielded the best results in most cases,\nsuggesting that moderate weighting of the pathway loss term helps balance the influence of\nindividual gene expression and biological structure. A more detailed description of this sensitivity\nanalysis is provided in the Supplementary Materials (Figure S18).\nFigure R15. Sensitivity analysis of model performance across different regulation levels. Each\npanel corresponds to a fixed pathway loss parameter (from left to right: 0.1, 0.5, 1, 2, and 5), while\nthe x-axis represents increasing reconstruction parameter. The y-axis shows the Pearson\ncorrelation between the predicted and true gene expression values. Lines represent different\nneighborhood loss parameter: red = 1, blue = 2, and green = 5.\nThanks for your suggestion. To assess the contribution of each component in PASTA's loss\nfunction, we conducted an ablation study by systematically removing individual terms and\nevaluating the resulting imputation performance. Specifically, we compared the full model\n(Baseline) against three ablated versions: (1) No neighborhood, which removes the neighborhood\nsimilarity loss; (2) No pathway, which removes the pathway loss; and (3) No reconstruction, which\nremoves the reconstruction loss from the mode. Our model requires global similarity, which cannot\nbe eliminated. Performance was evaluated using three metrics: Pearson correlation, Spearman\ncorrelation, and mean squared error (MSE) between the imputed and ground-truth pathway\nexpression (Figure R16 and Figure S17). Across 30 simulation runs, removing any of the loss terms\nled to a decline in performance. Notably, omitting the pathway loss caused the largest drop: median\nof Pearson’s correlation from 0.91 to 0.5. The reconstruction loss and the neighborhood loss have\na more moderate effect, while they still led to performance drop demonstrating their roles in\nrefining the imputation stability.\nThese results confirm that all components of the loss function contribute meaningfully to the\noverall performance of PASTA, and that their combined effect yields robust and accurate pathway\nimputation.\nFigure R16. Ablation analysis. Boxplots show the distribution of performance across 30\nsimulation runs for the three evaluation metrics (Pearson correlation (left), Spearman correlation\n(center), and mean‑squared error (MSE) (right)) between imputed and ground‑truth expression\nunder different loss function settings. Baseline includes all the loss functions. When pathway loss\nmissed, the imputation accuracy saw a big drop: median of Pearson’s correlation from 0.91 to 0.5.\nNeighborhood loss and the reconstruction loss also helps on the imputation stability.\nloss function is modified.\nBased on the comments above, I do not think this manuscript is suitable for publication in\nNature Communications.\nWe appreciate the reviewer’s feedback and the opportunity to clarify the novelty of our work. To\nbetter highlight the contributions of PASTA, we would like to outline the key innovations and\nadvantages of our method.\n1. Methodological innovation\nPASTA introduces a pathway-oriented design that fundamentally differs from traditional\ngene-level imputation methods. Unlike traditional methods that impute individual gene\nexpressions, PASTA is designed to directly optimize pathway-level expression. To our\nknowledge, no existing methods optimize for pathway-level expression directly or tailor\nthe learning process toward functional, pathway-based biological interpretation.\nWe chose to model pathway expression rather than individual gene expression because\npathway scores aggregate multiple genes with related biological functions, reducing noise\nand enhancing interpretability. In doing so, PASTA shifts the focus from individual gene\nfluctuations to broader, more stable biological signals.\nFor each biological pathway, PASTA estimates pathway activity by aggregating both the\nobserved expression of spatially measured genes and the imputed expression of\nunmeasured genes derived from scRNA-seq, guided by the learned alignment. The\nalignment matrix is trained by minimizing a multi-term loss function, each component of\nwhich encodes a specific biological prior. Uniquely, PASTA focuses the optimization on\nthe pathway of interest by aligning the expression of pathway genes between the two\nmodalities and enforces cell type-aware local coherence by comparing the imputed and\nobserved pathway gene expressions in the neighborhood of each spatial cell. By integrating\nboth global and local constraints, PASTA’s architecture aligns well with biological\nprinciples such as cell type homogeneity and spatial continuity principles often overlooked\nin existing methods (as summarized in Table R1).\n2. Improved imputation accuracy in terms of pathway expression.\nWe evaluated PASTA across four spatial transcriptomics datasets generated by diverse\ntechnologies (Xenium, MERFISH, ISS, and CosMx). By aligning the imputation process\npathway expressions, our method consistently achieves higher accuracy in terms of\npathway expression. For example, in the CosMx frontal cortex dataset, PASTA achieved a\nmedian Pearson correlation of 0.80 between predicted and ground-truth pathway\nexpression—substantially higher than existing methods such as SpaGE and stPlus. These\nimprovements are consistent across other datasets as well (see Figures 2–6 in the\nmanuscript)\n3. Enhanced Downstream Biological Insight\nBeyond improved accuracy, the pathway-based design of PASTA enables deeper biological\ninsights in downstream analyses such as clustering, pseudo time inference, and differential\nexpression.\nFor example, in the Merfish brain dataset, we used the imputed pathway expression to\ncluster the cells, which can reduce the cell type mixture in the clusters (Figure 4 in the\nmanuscript). In the Xenium lung cancer dataset, the hallmark pathways can highlight\ntumor/non-tumor area based on their differential expressions (Figure 2E in the manuscript).\nThese outcomes represent a conceptual and practical advance over prior approaches.\nIn conclusion, the novelty of our method lies in its biologically grounded, pathway-focused design,\nwhich offers accurate pathway expression imputation and more possibilities on the downstream\nanalysis that uncover cell states and trajectories. We hope this clarification helps highlight the\nnovelty and broader relevance of our contribution. We have added a more detailed discussion of\nthe methodological and biological innovations in the revised manuscript.\nTable R1. Summary of several existing imputation tools for spatial transcriptomics data.\nTools Model ST cell type ST cell Pathway\ninformation neighborhood information\nPASTA Trainable\nalignment matrix\nTangram Deep learning\nalignment\nSeurat CCA+MNN\nStDiff Diffusion model\nStplus Autoencoder\nGimVI Variational\ninference\nSpaGE Alignment +\nKNN\nstAI Deep learning-\nbased model\nReferences:\n1. Philippot, Camille, et al. \"Astrocytes and oligodendrocytes in the thalamus jointly maintain\nsynaptic activity by supplying metabolites.\" Cell Reports 34.3 (2021).\n2. Manu, Doina Ramona, et al. \"Astrocyte involvement in blood–brain barrier function: a\ncritical update highlighting novel, complex, neurovascular interactions.\" International\njournal of molecular sciences 24.24 (2023): 17146.\n3. Jin, Kelly, et al. \"Brain-wide cell-type-specific transcriptomic signatures of healthy ageing\nin mice.\" Nature 638.8049 (2025): 182-196."}, {"round": 2, "reviewer_comments": {"reviewer_1": "1.You noted OEC and vascular cells differed due to cell number availability. Could intrinsic biological variability of those cell\ntypes also play a role? In addition, in a small dataset setting, when the number of cells are all small how does it perform?\n2.You recommend using dozens to hundreds of pathways. Could you show a quantitative stability analysis — e.g., clustering\nconsistency or pseudotime correlation — when using fewer vs. more pathways?", "reviewer_2": "I acknowledge the authors’ revisions and the improvements made in clarifying methodological details and providing\nadditional analyses. However, several critical issues remain insufficiently addressed, particularly regarding the validity of\npathway-level aggregation, the meaningful contexts for pathway prediction, and the rigor of the evaluation strategy. Without\nfurther clarification, these concerns limit the strength of the methodological claims and the scope of biological interpretations.\nMy detailed comments are as follows.\nMajor:\n1. The manuscript directly predicts pathway-level expression by aggregating the expression of genes within each pathway.\nHowever, genes within the same pathway are often heterogeneous, with some acting as positive regulators and others as\nnegative regulators. Simply summing gene expression to represent pathway expression may obscure such gene-specific\ncontributions and potentially lead to misleading interpretations. A more rigorous justification for this aggregation strategy\nshould be provided, together with a discussion of its potential limitations.\n2. The manuscript should further clarify the circumstances under which pathway-level prediction is truly meaningful. When a\npathway contains measured genes, its expression can already be estimated from these genes, as the evaluation itself\ndefines the sum of measured gene expression as the ground truth. In such cases, the added value of pathway prediction\nremains unclear. A more convincing design would be: demonstrate utility in pathways without any measured genes, or\nevaluate predictions at the single-gene level for unmeasured genes in partially observed pathways, or show that predicted\npathway expression reveals additional biological insights compared to pathway expression calculated directly from\nmeasured genes, rather than compared to single-gene measurements as currently presented.\n3. The current evaluation strategy is problematic. Pathway expression is defined as the sum of ST-observed genes, which\nare also used for model training. This conflates training and evaluation, essentially fitting on the training set. To ensure\nindependence and rigor, all experiments should adopt gene-mask cross-validation. Gene-mask better reflects the\nunmeasured-gene scenario and provides a fairer test of generalization. If pathway integrity is a concern, three\ncomplementary evaluations could be reported: held-out genes (generalization), observed genes (reconstruction fidelity), and\ncombined genes (overall consistency and utility).\n4. Given that PASTA is a pathway-oriented framework, the authors should provide more details regarding pathway selection\nand coverage. Specifically: for each dataset, how many pathways are present in the database, and how many were included for\nprediction? Were unselected pathways mainly irrelevant to the tissue or disease context? How does PASTA perform on\nthose excluded pathways? Provide statistics on the overlap between pathway genes and spatially measured genes, and\nanalyze the relationship between overlap degree and pathway prediction accuracy. Such analyses would help readers\nunderstand the scope of PASTA’s applicability across different biological contexts, and provide practical guidance for\npathway selection in future studies.\nMinor:\n1. The authors should clearly specify the default parameter values used for each component of the loss function.\n2. On line 341, there is an extra comma that should be removed.", "reviewer_3": "The revision has addressed my concerns."}, "author_response": "Introduction\nWe thank the editor and reviewers for another round of thoughtful review for our manuscript\n“Accurate Imputation of Pathway-specific Gene Expression in Spatial Transcriptomics with\nPASTA”. In this revision, we provided additional clarification for the rationale behind the score\nthresholds in the recommendation chart. The reviewers’ comments are in bold black, and our\nreplies are in blue.\nResponse to Reviewer 1’s comments:\n1. You noted OEC and vascular cells differed due to cell number availability. Could\nintrinsic biological variability of those cell types also play a role? In addition, in a small\ndataset setting, when the number of cells are all small how does it perform?\nThanks for the insightful comments. To our knowledge, when pathways are closely linked to\ncell activity, biological variability is unlikely to be the main driver of prediction performance\ndifferences. In the small dataset setting, we observed minimal variation in prediction\nperformance across most cell types. The only exception was the CB GABA cells, which\nshowed a lightly lower correlation. However, this appears to be related to the sparse spatial\ndistribution of these cells across the spatial domain rather than biological variability.\nFigure R1. Evaluation of imputation performance across different cell types in 30\nsimulation datasets. Nine distinct cell types from the Merfish brain dataset, including:\nAstro-Epen (astrocyte and ependymal-like), CB (cerebellum) GABA (Gamma-aminobutyric\nacid;), CNU-LGE (cerebral nuclei - lateral ganglionic eminence) GABA, Immune, OB-IMN\n(olfactory bulb - immature neurons) GABA, OEC (olfactory ensheathing cell), OPC-Oligo\n(oligodendrocyte progenitor cell), and Vascular.\n2. You recommend using dozens to hundreds of pathways. Could you show a quantitative\nstability analysis — e.g., clustering consistency or pseudotime correlation — when\nusing fewer vs. more pathways?\nThanks for your comments. We compared clustering and pseudotime performance under\nconditions with fewer versus more pathways. For clustering, using 10 pathways resulted in\nmixed cell populations (Figure R2.B), whereas increasing the number to 36 pathways enabled\nclear cell separation (Figure R2.A). In the CosMx frontal cortex dataset, using 40 pathways\nfailed to capture the correct developmental trajectory (Figure R2.D), while including 200\npathways successfully reconstructed the expected temporal progression (Figure R2.C). The\noptimal number of pathways is related to the cell number and data complexity. A brief\ndiscussion has been found in the Methods section (page 19): “In our analysis, the input to\nMonocle3 is a cell-by-pathway expression matrix, where each row corresponds to a pathway’s\nexpression. we originally included 49 pathways for the Xenium lung cancer dataset and 292\npathways for the CosMx frontal cortex dataset. The number of pathways used can vary\ndepending on tissue complexity and heterogeneity; incorporating a broader and biologically\nrelevant set of pathways generally improves the resolution and robustness of trajectory\ninference. In practice, including several dozen to around a hundred well-chosen pathways\noften yield more meaningful and stable trajectories.”.\nFigure R2. Downstream analysis with fewer and less pathways using Merfish brain\nand Cosmx frontal cortex dataset. A. 10 pathways were used to generate the clusters in\nthe Merfish brain dataset. B. 35 pathways were used to generate the clusters in the Merfish\nbrain dataset. C. 200 pathways were used to calculate the pseudotime in the Cosmx frontal\ncortex dataset. D. 40 pathways were used to calculate the pseudotime in the Cosmx frontal\ncortex datase. More pathways tend to provide more reliable information.\nResponse to Reviewer 2’s comments:\nMajor\n1. The manuscript directly predicts pathway-level expression by aggregating the\nexpression of genes within each pathway. However, genes within the same pathway\nare often heterogeneous, with some acting as positive regulators and others as\nnegative regulators. Simply summing gene expression to represent pathway\nexpression may obscure such gene-specific contributions and potentially lead to\nmisleading interpretations. A more rigorous justification for this aggregation strategy\nshould be provided, together with a discussion of its potential limitations.\nThanks for the comment. We understand the concerns about the aggregating the expression\nof genes within each pathway. We would like first to clarify that aggregating genes from a\npathway to represent the pathway-level expression is a commonly used method in terms of\nthe pathway analysis. For example, the “AddModuleScore” function in the Seurat package\ndefines a module as a set of genes (i.e., a pathway) and, for each cell, aggregates (averages)\nthe expression levels of the genes within the module. This function has been broadly applied\nin bioinformatics studies, including tumor progression prediction [1], cancer prognosis [2],\nand disease subtype identification [3].\nWe also acknowledge that aggregation may obscure heterogeneity among genes. Therefore,\nwe recognize the importance of studying single-cell behavior using existing imputation\nmethods such as Tangram, Seurat, StDiff, and STplus. However, as demonstrated in our\nexperiments, when ST data contain only a limited number of genes, single-gene imputation\nmethods often produce noisy predictions. In such cases, our proposed pathway-level\nimputation approach offers a more robust alternative by leveraging pathway-level information\nto improve the accuracy of inferred expression patterns. We emphasize in the discussion\nsection that our method is not intended to replace existing single-gene imputation tools. Rather,\nit is designed to complement them—particularly in scenarios where pathway-level analysis is\nthe focus and the ST data have limited genomic coverage.\n2. The manuscript should further clarify the circumstances under which pathway-level\nprediction is truly meaningful. When a pathway contains measured genes, its\nexpression can already be estimated from these genes, as the evaluation itself\ndefines the sum of measured gene expression as the ground truth. In such cases, the\nadded value of pathway prediction remains unclear. A more convincing design would be:\ndemonstrate utility in pathways without any measured genes, or evaluate predictions\nat the single-gene level for unmeasured genes in partially observed pathways, or show\nthat predicted pathway expression reveals additional biological insights compared to\npathway expression calculated directly from measured genes, rather than compared to\nsingle-gene measurements as currently presented.\nThanks for your comments. We agree that evaluating the prediction accuracy for unobserved\ngenes is essential for demonstrating the validity of the proposed method. To this trend, we\nconducted a cross-validation experiment using a leave-one-gene-out strategy. Specifically, we\npredicted the expression of each gene after excluding it from the training data. Across 20\nspatial simulation datasets, the median Pearson’s correlation between predicted and\nobserved expression for the held‑out genes was comparable to Tangram under the same gene\nmasking conditions (Figure R3.A).\nTo demonstrate the advantage of PASTA, we implemented an additional 2‑fold gene‑mask\ncross-validation strategy within each pathway. In every fold, we randomly masked 50% of the\ngenes belonging to the pathway (treating them as unmeasured in the ST data), trained the\nmodel using the remaining pathway genes along with shared genes between the scRNA-seq\nand ST datasets, and then computed pathway scores using the imputed values for the masked\ngenes. These masked genes were excluded from both the input features and the loss\nfunctions, ensuring that no information leakage occurred. Across the 20 simulation runs,\nPASTA consistently outperform Tangram, achieving a higher median Pearson’s correlation\nbetween predicted and observed held-out genes (Figure R3.B).\nThe motivation for developing pathway-based imputation originated from a collaborative\nproject focused on evaluating TCF pathway expression in T cells. In this project, only a limited\nnumber of genes were profiled, resulting in a limited number of genes were profiled, the\nexpression data were extremely sparse in the real data (Figure R4.A), making further analysis\ndifficult. After applying our imputation method, we observed that the TCF pathway expression\nbecame smoother across cells (Figure R4.B). Since this dataset has not yet been published,\nwe did not include its details in the manuscript. At this stage, we did not conduct additional\nwet-lab experiments to validate our findings for this dataset. In order to assess the validity of\nour imputation results, we utilized published spatial transcriptomics datasets to perform\ndownstream analyses, including pathway enrichment, clustering, and trajectory analyses.\nIn the real data applications presented in the manuscript, we evaluated how predicted pathway\nexpression provides additional biological insights compared with pathway expression\ncalculated directly from the measured genes in the Xenium lung cancer dataset. Using only\nthe measured genes failed to capture correct enrichment patterns across cell types. By\ncontrast, PASTA successfully predicted enrichment of cancer cell development related\npathways (MYC TARGETS V1 and V2) in tumor cells. Furthermore, PASTA identified immune-\nrelated pathways, including INFLAMMATORY RESPONSE, IL2–STAT5 SIGNALING, and\nALLOGRAFT REJECTION as highly expressed in T cells and macrophages (Figure R5).\nRelated descriptions and discussions have been provided in the manuscript Result section,\npage 9, paragraph 3.\nFigure R3. Cross-validation. A. Leave-one-out cross-validation at gene level. Pearson’s\ncorrelation between predicted and observed expression for the held‑out genes were summarized\nin the boxplot from 20 simulations. B. Within‑pathway cross‑validation. For every pathway, genes\nwere randomly split 50/50 in each fold. The first half was retained as observed in the simulated\nST data and used to train the model; the second half was masked and later imputed. Each dot\nrepresents one simulation run; horizontal bars indicate the median performance across runs.\nFigure R4. Pathway expression examination across T cell. A. Sparse expression data was\nobserved in the original spatial dataset. B. Predicted pathway expression among the T cells.\nPathway expression became smoother across cells after the imputation process.\nFigure R5. Enrichment analysis across cell types. Using only the measured genes failed to\ncapture correct enrichment patterns across cell types. The predicted pathways showed\nenrichment of cancer cell development related pathways in tumor cells and enrichment of immune\nrelated pathways in T cells.\n3. The current evaluation strategy is problematic. Pathway expression is defined as the\nsum of ST-observed genes, which are also used for model training. This conflates\ntraining and evaluation, essentially fitting on the training set. To ensure independence\nand rigor, all experiments should adopt gene-mask cross-validation. Gene-mask better\nreflects the unmeasured-gene scenario and provides a fairer test of generalization. If\npathway integrity is a concern, three complementary evaluations could be reported:\nheld-out genes (generalization), observed genes (reconstruction fidelity), and\ncombined genes (overall consistency and utility).\nThanks for your comments. As described in the previous response, we conducted cross-\nvalidation at both the gene-level and the pathway level (Figure R3). We first performed the\nleave-one-out simulation to examine the prediction accuracy at the gene level. Across 20\nspatial simulation datasets, the median Pearson’s correlation between predicted and\nobserved expression for the held‑out genes achieved comparable performance with Tangram\nunder identical masking. We then adopted another 2‑fold gene‑mask cross-validation strategy\nwithin each pathway. In every fold, we randomly masked 50% of the genes belonging to the\npathway (treating them as unmeasured in the ST data), trained the model using the remaining\npathway genes along with shared genes between the scRNA-seq and ST datasets, and then\ncomputed pathway scores using the imputed values for the masked genes. These masked\ngenes were excluded from both the input features and the loss functions, ensuring that no\ninformation leakage occurred. Across the 20 simulation runs, PASTA consistently outperform\nTangram, achieving a higher median Pearson’s correlation between predicted and observed\nheld-out genes.\nWe further conducted a simulation to compare the held-out genes, observed genes, and\ncombined genes. Under 30 simulation datasets, PASTA consistently showed higher\ncorrelations between prediction and the true pathway expression among held-out genes\n(Pearson correlation median 0.5 by PASTA VS. 0.18 by Tangram), overlapped genes (Pearson\ncorrelation median 0.89 VS. 0.60), and combined genes (Pearson correlation median 0.85 VS.\n0.57) (Figure R6). For the real dataset, since we do not have ground truth for the unobserved\ngenes (also combined genes) in ST dataset, we present the correlation of the observed genes\nand used the downstream analysis to validate the whole pathway (combined genes). A\ndiscussion of these results has been added to the Simulation section (page 6 paragraph 1).\nFigure R6. Evaluation of imputation performance for different part of pathways. Boxplots\ndisplay the distribution of correlation for each method. Our method shows better performance\nthan Tangram among heldout genes, overlap genes, and the whole pathways.\n4. Given that PASTA is a pathway-oriented framework, the authors should provide more\ndetails regarding pathway selection and coverage. Specifically: for each dataset, how\nmany pathways are present in the database, and how many were included for\nprediction? Were unselected pathways mainly irrelevant to the tissue or disease\ncontext? How does PASTA perform on those excluded pathways? Provide statistics on\nthe overlap between pathway genes and spatially measured genes, and\nanalyze the relationship between overlap degree and pathway prediction accuracy. Such analyses\nwould help readers understand the scope of PASTA’s applicability across different\nbiological contexts, and provide practical guidance for pathway selection in future\nstudies.\nThank you for your comments. Given that pathway databases contain thousands of pathways,\nwe focused our evaluation on subsets that are biologically relevant to each dataset’s tissue\ncontext. Specifically,\n• For the Xenium lung cancer data, we examined 49 hallmark pathways (Table S1.1).\n• For the MERFISH brain data (36 pathways, Table S1.2) and ISS visual cortex data (16\npathways, Table S1.3), we selected Gene Ontology Biological Process (GOBP)\npathways related to mouse brain development.\n• For the CosMx dataset, we evaluated 40 GOBP pathways associated with human\nbrain development (Table S1.4).\nThe details of the pathways were summarized in the Supplementary material Table S1.\nPathways outside these curated sets were largely irrelevant to the corresponding tissue\ncontext and exhibited uniformly low expression across cells, as expected. For example, when\ntesting a GOBP pathway unrelated to cancer or lung tissue development (Coronary artery\nmorphogenesis, Figure R7.B), the predicted pathway expression values remained\nconsistently below 0.01 across all cells.\nWe have also included the summary statistics on the overlap between pathway genes and\nspatially measured genes as a supplement table (Table S1). In order to analyze the\nrelationship between overlap degree and pathway prediction accuracy, we conducted a\nsimulation study based on the overlap proportion between the pathway and spatial dataset.\nWe observed that when all or some of the genes in a pathway are measured, as shown in\nFigure 1 in the manuscript, the imputation performance remains high with a median Pearson’s\ncorrelation of around 0.8. In contrast, when none of the pathway genes are observed, all\nmethods experience a significant drop of median Pearson’s correlation to around 0.3 in\nperformance (Figure R8 in this response letter, Figure S2 in manuscript). Despite this, our\nmethod consistently outperforms alternative approaches, with Wilcoxon’s test P-value less\nthan 0.01, demonstrating its robustness even under limited gene coverage.\nFigure R7. Pathway predictions. A. HEDGEHOG SINGALING from hallmark pathway database\nthat are related to cancer development. B. CORONARY ARTERY MORPHOGENESIS from\nGOBP pathway database that are related heart vessel organization. The irrelevant pathways in\ngeneral show low expression across all the cells.\nFigure R8. Impact of pathway gene coverage on pathway prediction accuracy. Performance\nis shown for three evaluation metrics, Pearson correlation (left), Spearman correlation (middle),\nand mean‑squared error (MSE, right), between imputed and true expression of the masked genes.\nThe x‑axis denotes the proportion of each pathway’s genes that are observed in the ST panel (0,\n1/2, 1/3, and 1/4). Box‑plots represent the distribution across 30 independent simulation runs.\nMinor\n1. The authors should clearly specify the default parameter values used for each\ncomponent of the loss function.\nThanks for your suggestion. We have clarified this in the method section (page 15, paragraph\n5).\n2. On line 341, there is an extra comma that should be removed.\nThanks for your insights. We removed the extra comma.\nReference\n1. Luo, Linjie, et al. \"Single-cell RNA sequencing identifies molecular biomarkers predicting late\nprogression to CDK4/6 inhibition in patients with HR+/HER2-metastatic breast cancer.\"\nMolecular cancer 24.1 (2025): 48.\n2. Zhang, Zhanshuo, et al. \"Five-gene prognostic model based on autophagy-dependent cell\ndeath for predicting prognosis in lung adenocarcinoma.\" Scientific Reports 14.1 (2024): 26449.\n3. Zethoven, Magnus, et al. \"Single-nuclei and bulk-tissue gene-expression analysis of\npheochromocytoma and paraganglioma links disease subtypes with tumor\nmicroenvironment.\" Nature Communications 13.1 (2022): 6262.\n4. Barbie, David A., et al. \"Systematic RNA interference reveals that oncogenic KRAS-driven\ncancers require TBK1.\" Nature 462.7269 (2009): 108-112.\n5. Aibar, Sara, et al. \"SCENIC: single-cell regulatory network inference and clustering.\" Nature\nmethods 14.11 (2017): 1083-1086."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have addressed all my concerns and significantly improved the paper quality. Thanks for their diligent efforts.", "reviewer_2": "The authors have satisfactorily addressed all my previous comments. I have no additional comments at this stage."}, "author_response": ""}]

Decision letter