Robust and interpretable prediction of gene markers and cell types from spatial transcriptomics data

Domain: physics
Published: 2026-01-16
Editorial outcome: major_revision

Abstract
Spatial transcriptomics (ST) links tissue morphology with gene expression values, opening new avenues for digital pathology. Deep learning models are used to predict gene expression or classify cell types directly from images, offering significant clinical potential but still requiring improvements in interpretability and robustness. We present STimage as a comprehensive suite of models to predict spatial gene expression and classify cell types directly from standard H&E images. STimage enhances robustness by estimating gene expression distributions and quantifying both data-driven (aleatoric) and model-based (epistemic) uncertainty using an ensemble approach with foundation models. Interpretability is achieved through attribution analysis at single-cell resolution integrated with histopathological annotations, functional genes, and latent representations. We validated STimage across diverse datasets, demonstrating its performance across various platforms. STimage-predicted gene expression can stratify patient survival and predict drug response. By enabling molecular and cellular prediction from routine histology, STimage offers a powerful tool to advance digital pathology.

Introduction
Haematoxylin and Eosin (H&E) staining of tissue samples has been used by pathologists for more than a century. However, microscopic examination of H&E images is often time-consuming and variable1. Moreover, a range of molecularly distinct cell types may not be distinguishable by eye2. Recent FDA approvals3 of whole-slide imaging (WSI) for primary diagnostic use have encouraged more widespread adoption of digital WSI4. However, the development and optimisation of automated cancer classification models requires a large number of pathologist-annotated tissue images5,6. HE2RNA7 was a pioneering model trained without pathological annotation. Instead, HE2RNA uses matched H&E images and RNA-seq data from the TCGA database for 8725 patients across 28 different cancer types, to predict gene expression at the tile level (each WSI H&E image is split into multiple tiles). However, due to the nature of the data input, HE2RNA lacks ground truth at the tile level, relying instead on the aggregated gene expression value from the whole image (i.e., all tiles in the image have the same ground truth for one gene), thus limiting its ability to assess tumour heterogeneity.

Spatial transcriptomics (ST-seq) produces both imaging and sequencing information and has the unique capability of measuring over 20,000 genes without tissue dissociation, preserving tissue anatomy and the microenvironmental context8. However, spatial transcriptomics remains prohibitively expensive for clinical applications. To improve scalability, researchers have attempted to predict spatial gene expression using deep learning methods, including STnet9, Hist2Gene10, Hist2ST11, and DeepSpace12. However, these methods are still limited in interpretability, robustness, and accuracy. Crucially, clinical applications require models to be understandable to clinicians and patients for safety, trust, validation, and therapeutic decision-making13.

Here, we report STimage, a deep learning probabilistic framework with reduced uncertainty and improved robustness and interpretability, applicable even in cases where training datasets are limited. We tested STimage across three different cancer datasets, including breast, skin, and kidney cancer, as well as a liver immune disease dataset, assessing the performance of predicting the most variable genes or functional genes. The shared biological relationships between these genes are learned simultaneously, allowing us to assess the performance of gene panels. We introduced a geometric assessment metric that better represents the model’s ability to predict spatial patterns.

Results
STimage model-distribution prediction, ensemble, interpretability and uncertainty estimates
First, in the STimage pipeline, we implement the image preprocessing step (Fig. 1a). STimage performs stain normalisation to account for technical variation between image inputs used during model training. To remove uninformative tiles, we calculate tissue coverage and filter out tiles with low coverage. In STimage, image tiles are either used to fine-tune the ResNet50 model (as a mini batch of 64 tiles by default) or passed through pretrained foundation models to be transformed into feature embeddings (Fig. 1b). All latent features are then connected to a layer of neurons representing the parameters of the negative binomial (NB) distribution for each gene, which is optimised using a log-likelihood loss function. Alongside the distribution estimates, we also quantify model prediction uncertainty (Fig. 1b). This is a critical feature not available in existing tools, yet essential given the high variability expected in cancer datasets. Multiple genes can be trained simultaneously, with each gene corresponding to a pair of NB parameters. By sampling from the NB distribution, STimage can learn to predict the expression of individual genes or gene groups. Importantly, we incorporate model interpretability via Local Interpretable Model-agnostic Explanations (LIME) (Fig. 1b).

Several existing programs for gene expression prediction, such as STnet9, Hist2ST11, Hist2Gene10, DeepSpace12, and Xfuse14, produce fixed-point estimates without any uncertainty quantification. In contrast, STimage introduces an approach to quantify and mitigate uncertainties in gene expression predictions, distinguishing between uncertainty due to model knowledge gaps (epistemic) and inherent variability in the data (aleatoric) (Fig. 1b).

Prediction performance metrics and STimage accuracy
To assess model performance, we first visually compared the prediction of an abundant cancer marker (COX6C), broadly expressed across the tissue, with a less abundant marker, KRT5, expressed specifically in myoepithelial cells surrounding Ductal Carcinoma In Situ (DCIS) regions. The model trained on nine Visium breast cancer datasets was used to predict expression in a new Visium test dataset, showing highly consistent spatial patterns between predicted and observed values (Fig. 1c). Quantitatively, using the spatial autocorrelation metric Moran’s I, we observed a strong positive correlation (i.e., high-high or low-low clustering) for COX6C expression across most spatial spots (Fig. S2a). We applied the same model, trained on the nine fresh frozen Visium samples, to two external, previously unseen out-of-distribution FFPE samples: one measured using the version-1 FFPE protocol and the other with the FFPE CytAssist protocol (Fig. 1d). Results for one representative cancer-immune marker, ATP1A1, are shown for both FFPE test samples (Fig. 1d). Results for all 14 cancer-immune markers, evaluated using PCC and Moran’s I metrics, demonstrate consistently positive scores, although performance varies across samples (Fig. S2c–e). This spatially-aware quantification offers a distinct advantage over non-spatial models such as HE2RNA7, which do not provide per-tile gene expression predictions.

Prediction performance-more robust across diverse datasets and sample types
To assess the robustness of the approach, we increased heterogeneity in the training and testing data by including samples from different protocols, laboratories, and populations. We trained a model using nine Visium datasets from fresh frozen tissues, with three datasets generated by 10x Genomics and six by Swarbrick’s Laboratory15. Using an ensemble approach, the model trained on fresh frozen samples was applied to two external FFPE samples processed using Visium protocols. The top predicted gene, ATP1A1, showed high PCC values (0.51 and 0.52) and consistent gene expression patterns (Fig. 1d). Assessing model performance across 1522 functional genes (trained simultaneously), we found that predictions were positively correlated with measured values for 1136 of 1522 genes. The top 10 genes achieved PCC scores around 0.4 (Fig. 1e). Notably, model performance was consistently higher when our image normalisation strategies were applied (Fig. S3). We also thoroughly assessed the impact of image tile resolution and found that the default 299×299 tile size in STimage yielded the highest predictive performance (Fig. S4). These findings highlight options for improving robustness.

The trained model incorporated both the public dataset from 10x Genomics and Swarbrick’s dataset, while testing was performed on an independent dataset generated by a European laboratory, consisting of 33 tissue samples. STimage was benchmarked against six other approaches: five software tools, THItoGene, IGI_DL, STnet, His2gene, and Hist2ST, and an alternative model, ViT_NB. STimage consistently outperformed all other methods across the majority of test samples, except for six (out of 33), where its performance was comparable to the top-performing models (Fig. 2e). To further enhance performance, we developed an ensemble model, FSTimage, which integrates multiple foundation-STimage models. Each model combines a foundation model as the upstream image extractor with the STimage negative binomial module. The foundation models used include GigaPath16, H-Optimus17, Virchow218, Phikon19, CTransPath20, CONCH21, UNI22, and Virchow23. Final FSTimage predictions were obtained by averaging results across these models. Our findings indicate that STimage, Virchow2_STimage, and FSTimage consistently outperformed all other models compared.

To further evaluate robustness and increase sample diversity, we incorporated additional datasets, including kidney cancer (fresh frozen, n = 615,632 spots), non-cancer liver disease (FFPE, n = 419,967 spots), and skin cancer (FFPE, n = 1014,405 spots) (Fig. S7, S6, S5). We used a leave-one-out model training and testing strategy to evaluate model performance on individual patient samples. Among the top 100 most predictable genes, Gene Ontology GO enrichment analysis revealed terms corresponding to relevant tissue biology across the three GO databases (Fig. S7a–c, S6a–c, S5a–c). The PCC between predicted and ground truth gene expression values among the top 100 genes showed a positive correlation, ranging from 0.2 to 0.8 (Fig. S7d, S6d, S5d). Performance also varied across patients, suggesting tissue heterogeneity. Representative spatial expression plots are shown in Fig. S7e–g, S6e–g, S5e–g.

Prediction uncertainty can be quantified for each gene
In addition to improving robustness, ensemble models allow us to quantify uncertainty arising from both the data (i.e., aleatoric uncertainty) and the model itself (i.e., epistemic uncertainty). Existing methods, such as STnet9 and HisToGene10 use a CNN or ViT as an image feature extractor. Both models rely on multi-layer neural networks to predict fixed gene expression values, without estimating variability in their predictions. Hist2ST11 employs a zero-inflated negative binomial distribution to predict gene expression distributions, but it does not account for epistemic uncertainty. To evaluate these capabilities in STimage, we used the model trained on the nine Visium fresh frozen breast cancer datasets to test its performance on two independent, unseen FFPE tissues from separate cohorts (Fig. 2a–d and S1). Two representative genes were selected for comparison: ATP1A1, a breast cancer marker, and CD63, an immune marker. For each gene, aleatoric uncertainty, epistemic uncertainty, and predicted gene expression were visualised in spatial tissue plots (Fig. 2a, b). Aleatoric uncertainty (i.e., the average of predicted NB variance values across single models) was found to be higher than epistemic uncertainty (i.e., the variance of NB mean estimates) in both FFPE samples. As expected, both uncertainty types showed a positive correlation with gene expression levels. Interestingly, lower epistemic uncertainty often coincided with higher aleatoric uncertainty, as observed for ATP1A1 (Fig. 2b), thereby increasing confidence in the predicted NB parameters (i.e., mean and variance).

In addition to lower epistemic uncertainty, ATP1A1 also showed higher predictive performance (PCC: 0.52; Fig. 2b) compared to CD63 (PCC: 0.34; Fig. 2a). This finding aligns with interpretations of epistemic uncertainty, which reflects the confidence in model predictions24. In this case, the higher epistemic uncertainty observed for CD63 may be partially attributed to the fact that the test samples were generated using different technological platforms compared to those in the training data. As a result, the input data were considered out-of-distribution.

An ensemble approach improves uncertainty performance
To increase model robustness across a broad range of diverse, unseen datasets, we trained an ensemble of STimage regression models. Ensemble models have been shown to improve generalisation in the transcriptomics context25. We observed variation in model predictions, even when using the same training dataset and model architecture, but with different random seeds. To assess and mitigate this variability, we analysed the distribution of performance across individual models within 12 different groups and for two distinct tissue slides (a total of 30 models). The results are shown as boxplots in Fig. 2c. Each box represents the average PCC from independent predictions made by 10 single models for each of 12 genes (out of a total of 14 trained genes; 2 were not presented in the additional CytAssist FFPE dataset). The lines represent ensemble models obtained by averaging predictions from all 10 single models within each group, with error bars indicating the standard deviation across different genes (12 genes). Although there is considerable variation among individual models, their performance remains consistent across ensemble groups. Furthermore, the ensemble models outperform the single models and exhibit reduced variability across all groups (Fig. 2c). We further demonstrate that this ensemble strategy also applies to the FSTimage_ens model, where embeddings from different foundation models are combined to further enhance performance (Fig. 2e).

Out-of-distribution performance for models trained on a small dataset
To assess model performance on out-of-distribution (OOD) data, we designed an experiment where the STimage model and three other models, STnet, HisToGene, and Hist2ST, were trained on fresh frozen samples (Visium poly-A captured data) and tested on FFPE data (Visium probe-hybridisation data). The probe-hybridisation data detected more genes and provided more information per gene. As the two datasets were derived from different tissue types and technologies, they were considered OOD (Fig. S1).

Previous studies, such as Hist2Gene10 and Hist2ST11, showed limited performance when trained on small datasets. In contrast, we found that the STimage model with a fine-tuning option improved performance on small training sets compared with models using ViT for feature extraction. For a robust comparison, we predicted expression levels for the top 1000 highly variable genes (HVGs) and computed both PCC and Moran’s I scores for each gene. Across all 1000 genes, STimage consistently outperformed the other three models.

We further evaluated robustness by testing the prediction performance for each of the 1000 HVGs, comparing STimage to the three other methods (Fig. S1). Notably, methods based on ViT (e.g., HisToGene and Hist2ST), which generate attention maps across tissue sections, require substantially more memory. These models become challenging to run on datasets with more spatial measurements, a scenario increasingly common with recent advances in spatial transcriptomics. STimage uses a CNN-based image extractor that treats each spatial measurement as an independent data point, enabling scalability to datasets with high spatial resolution. Additionally, STimage leverages pre-trained models for feature extraction, an approach not used by ViT-based methods.

We further evaluated robustness and increased sample diversity by incorporating additional datasets, including kidney cancer (fresh frozen, n = 6, 15,632 spots), non-cancer liver disease (FFPE, n = 4, 19,967 spots), and skin cancer (FFPE, n = 10, 14,405 spots) (Figs. S7, S6, S5). We used a leave-one-out training and testing strategy on tissue samples to evaluate model performance for each patient. Among the top 100 most predictable genes, GO enrichment analysis revealed terms consistent with the underlying tissue biology across the three GO databases (Figs. S7a–c, S6a–c, S5a–c).

The PCC between predicted and ground truth gene expression values among the top 100 genes ranged from 0.2 to 0.8 (Figs. S7d, S6d, S5d), demonstrating a positive correlation. Performance also varied among patients, highlighting tissue heterogeneity. Example spatial gene expression plots are shown in Figs. S7e–g, S6e–g, S5e–g.

Out-of-distribution performance on independent spatial transcriptomics datasets
In the following sections, we further evaluated the robustness of our model by testing it under more challenging scenarios, involving independent datasets from different cancer types (skin and breast cancer), diverse data types, and various experimental protocols, including low-resolution legacy ST, Visium, single-cell resolution Xenium, non-spatial TCGA images, and PhenoCycler Fusion data.

The model, initially trained on nine Visium breast cancer datasets using a functional gene list of 1522 genes, was applied to the HER2ST26 legacy ST dataset. In this experiment, H&E images from the HER2ST dataset were tiled at the spatial location of each spot, using the Visium spot size of 55 μm, which is smaller than the original legacy spot size of 100 μm. Gene expression predictions were generated for these smaller tiles, visualised, and evaluated against the ground truth measurements. Figure 3a displays one of the top predictable genes, ESR1, while additional examples, ERBB2 and B2M, are presented in Fig. S8. Figure S8a and b show samples with low and high prediction performance, respectively, while Fig. S8c provides pathological annotations for these samples. Even in low-performance predictions, we observed consistent spatial expression patterns when compared with both ground truth and pathological annotation. A quantitative assessment of a large number of genes in the HER2ST dataset is presented in Fig. S8d. This includes PCC values for 1146 predicted genes, derived by intersecting the 1522 predictable genes with the gene list available in the legacy ST data. We found that a model trained on the high-resolution Visium dataset was able to effectively predict data from the legacy Visium protocol. Prediction performance was highest for samples B1-B6 (median PCC values ranging from 0.4 to 0.6), moderate for samples C1-D6, G1-G3, and H1-H3 (median PCC values from 0.2 to 0.4), and lowest for samples A1-A6 and E1-F3. Notably, samples A1-A6 and E1-F3 consistently showed low prediction performance, even in in-distribution datasets, as shown in Fig. 2e.

In contrast to the previous OOD analysis using a reduced-resolution test dataset (Visium legacy), here we applied the model to higher-resolution data from Xenium. The same model, trained on Visium H&E images, was applied to images from the single-cell resolution Xenium dataset (Figs. 3c and S9). Subcellular gene expression measured in the Xenium breast cancer dataset was grouped to match the Visium spot size, serving as ground truth for comparison with model predictions. We again observed highly consistent spatial expression patterns for gene markers such as ESR1 (Fig. 3c), although the quantitative correlation was lower than for datasets of the same resolution (Fig. S9).

Out-of-distribution performance on another cancer type assessed by spatial proteomics datasets
The same model, trained on nine Visium breast cancer datasets and a functional gene list of 1522 genes, was applied to H&E images from Visium skin cancer data (Figs. 3b and S10), i.e., for a different cancer type. We demonstrate that a model trained on breast cancer can predict gene marker expression in skin cancer. Predictions for key melanoma markers such as SOX10 and CD34 were highly consistent with protein data measured by PhenoCycler Fusion (Fig. S10). The subcellular protein staining signal from the PhenoCycler Fusion skin cancer dataset was grouped to match Visium spot size, serving as ground truth for comparison with predictions from H&E images on the same slide (Figs. 3d, S10). The ability to predict expression patterns that match protein markers highlight the model’s promise for diagnostic applications.

STimage performs well when trained and tested on other cancer types
Using the same model architecture, we trained the STimage regression model on a skin cancer dataset (Fig. S5), a kidney cancer dataset (Fig. S7), and a non-cancer liver immune disease dataset (Fig. S6). A Leave-One-Out Cross Validation (LOOCV) strategy was used to assess model performance. To further evaluate the biological applicability of the gene expression prediction model on H&E images, we performed clustering on the predicted gene expression data to analyse tissue heterogeneity. Overall, the clustering results identify the structural organisation of skin layers (Fig. S18a) with expected gene markers. Notably, Cluster 1 corresponds to the melanoma region. We further utilised a set of marker genes associated with epidermal genes, dermal fibroblast genes, and sweat gland genes to characterise clusters. We found that epidermal gene markers were enriched in Cluster 0, while sweat gland genes were enriched in Cluster 5 (Fig. S18b). The spatial patterns of the clusters generally align with the H&E morphology.

The classification model provides complementary information on tissue regions with high vs low gene expression
In addition to predicting continuous gene expression values, STimage also implements a classification model that enables the prediction of high vs low gene expression directly from H&E images (Fig. S11c). Refer to “STimage classification model” for details on the classification model. In this model, continuous gene expression values were discretised into two categories: high vs low (Fig. S11a, b).

We evaluated the classification model’s performance on the top 100 predictable genes in two test samples, FFPE and 1160920F, using Moran’s I, AUC, and IoU metrics (Fig. S11c). As the aim of this model is to classify tissue regions into high vs low expression zones, IoU is the most informative metric. The top predictable genes showed consistent classification performance across both test samples (Fig. S11d, e). Since the classification model predicts at the tile level, we found that it performed well on tiles with high expression, but was noisier for tiles with intermediate expression levels. Examples of top-predictable genes include CD24, GNAS, and VEGFA (Fig. S11d, e). These genes are highly relevant to cancer; for instance, CD24 is known to be overexpressed in HER2+ tumours and is associated with poor prognosis27. To assess the classification model’s robustness, we compared two genes with contrasting expression distributions, CD52 and CD24 (Fig. S11). CD52 is an immune-related gene and a favourable biomarker in breast cancer28. Class probabilities from the model were consistent with the ground truth continuous gene expression values (Fig. S11g). All three performance metrics supported the model’s predictive accuracy for both genes.

Interpretability analysis
We investigated the interpretability of the model by assessing feature importance using LIME (Fig. S12; see Methods, “Model interpretation of STimage”). As examples, we selected VEGFA and CD52, as they were ranked among the top predictable genes and are known cancer and immune markers, respectively. We calculated LIME scores for the regression models (FFPE and 1160920F), focusing on the top five tiles with high VEGFA expression from cancer regions and the top five tiles with high CD52 expression from immune regions (Fig. S12c, d). LIME scores overlaid on H&E images, used for both qualitative and quantitative assessment, are shown in Fig. S12a, b. Two representative examples are displayed for two different test datasets: FFPE (Fig. S12a) and 1160920F (Fig. S12b).

The interpretability analysis shows that different nuclei contribute differently to gene expression prediction (Fig. S12a, b). We observed that nuclei with positive LIME scores for VEGFA have a larger area than those with positive LIME scores for CD52 in cancer regions, consistent with the biological expectation that cancer cells are typically larger than immune cells. When integrated with pathological annotations, we also found more nuclei with high LIME scores for VEGFA in cancer tiles (Fig. S12a, b) than in immune tiles. This observation is consistent with the known biology that VEGFA is overexpressed in cancer cells in HER2+ samples29. To support user engagement and exploration of interpretability, we developed an interactive web application that allows visualisation of LIME scores from the STimage model on arbitrary tiles.

Moreover, analysis of the latent space revealed that the learned features from ResNet50 could effectively distinguish cancer tiles from non-cancer tiles. We overlaid the true expression and predicted probability of CD24, one of the top predictable genes across both datasets, on the first two principal components of the latent space (Fig. S11f). CD24 is highly expressed in cancer regions and aligns well with pathologist annotations in Fig. S12a, b. In addition, saliency mapping confirmed that ResNet50 focuses on the most informative features (nuclei) in the images (Fig. S12).

Xenium cell type classification
The Xenium platform from 10x Genomics enables RNA transcript detection with subcellular resolution. Compared to Visium, Xenium-generated data is particularly well-suited for deep learning-based cell classification due to its single-cell resolution profiling. We obtained Xenium data for five primary breast cancer tissue sections and trained a deep learning model to classify nuclei using whole-slide H&E images. Figure 4 shows the classification metrics (Fig. 4a, b), indicating that cancer epithelial cells were the most accurately predicted, while other cell types also showed reasonable accuracy. Qualitatively, predictions closely matched the ground truth labels at both the macroscopic scale (Fig. 4d), where distinct tissue structures were well reconstituted, and the microscopic scale (Fig. 4e), showing individual nuclei segmentation and classification.

We also visualised the feature space learned by the model to interpret inter-cell type similarity and dissimilarity. The 256×256×64 feature map from the final layer of the classification branch was associated with the predicted cell types. Feature vectors were sampled for each cell type and the background class, and visualised using 2D t-SNE embeddings (Fig. 4c). The background was clearly separated from cells, while stromal, immune, and endothelial cells showed substantial overlap, suggesting potential classification challenges among these types.

Furthermore, we assessed model uncertainty by computing the entropy of predicted class probabilities across the image. We overlaid this uncertainty map with misclassified instances. Figure 4f shows that regions of high entropy were enriched for misclassified cells, while low-entropy regions largely corresponded to correct predictions.

We performed leave-one-out cross-validation (LOOCV) across the five Xenium datasets. Results presented in Figure S14 show consistent performance across test samples, demonstrating the robustness of our classification model. Among the cell types, cancer epithelial cells were most accurately predicted, likely due to their abundance and distinctive morphology. Immune cells (including B cells, T cells, and myeloid cells) formed the second most predictable group. This suggests potential clinical applications in profiling immune cell distribution in the tumour microenvironment, which is critical for cancer diagnosis. With more training data, particularly for rarer cell types, predictive performance could be further improved.

We benchmarked our model against the recently published CellViT model30, a Vision Transformer-based architecture that supports both single-cell segmentation and classification (Fig. S15). At the tile level, both CellViT and the STimage classification model accurately predicted cell segments. However, the STimage model outperformed CellViT in cell type prediction (Fig. S15a). We also compared mean IoU scores across LOOCV test samples (Fig. S15b), demonstrating that our model achieved superior performance.

Notably, the STimage classification model was trained using Xenium data in which training labels were generated automatically via image registration. This process aligns the DAPI-derived cell segmentation mask with cell type annotations obtained from marker gene expression measurements in the Xenium dataset, which are considered more accurate than those inferred from H&E alone. This automated approach reduces the need for manual labelling and minimises human error.

Evaluating clinical and translational potential of STimage using the TCGA dataset and a drug response dataset
To evaluate whether STimage could be applied to tissue images from a non-spatial dataset, we first tested two randomly selected breast cancer samples from The Cancer Genome Atlas (TCGA) (Fig. 1f–h). As ground truth for comparison, we obtained pathological annotations for these two H&E images. The predictions for cancer markers were consistent across the two TCGA replicates (Fig. 1f–h). We observed strong agreement between the predicted expression and pathological annotation, with high expression of cancer markers in tumour regions and low expression outside tumour areas.

We further demonstrate that the model-predicted gene expression can cluster cell types from H&E images, even in the absence of spatial measurements. This was evaluated using five tissue sections from three randomly selected TCGA patients. Louvain clustering based on the expression of 1522 predicted genes per tile is shown in Fig. S13a. The resulting clusters correlated with pathological annotations for each sample (Fig. S13a, b). Spatial plots of two known breast cancer markers, ESR1 and GATA331, also aligned with annotated tumour and normal regions (Fig. S13c, d).

For a broader analysis, we predicted expression levels of 1522 highly variable genes for 1034 H&E images from 670 TCGA patients. The average Pearson correlation coefficient between predicted and bulk-measured gene expression was 0.48 (Fig. 5a). STimage performed consistently well for both lowly expressed genes (e.g., MPPED1) and highly expressed genes (e.g., VEGFA), as illustrated in Fig. 5b. Notably, the predicted gene expression was also informative for patient survival (see Methods, “Validation using non-spatial data and assessing predictive values for survival and drug response”). Using the top five genes identified from univariate Cox regression models, we predicted gene expression for 1034 images with matched survival data. Survival curves stratifying patients into low- and high-risk groups were significant and consistent for both true and predicted expression values (Fig. 5c, d).

Top genes with the highest PCCs from gene-wise correlation analysis for each subtype were: FGFR3 (PCC = 0.46) and PAK1 (PCC = 0.47) for HER2; IL5RA (PCC = 0.19) and CCL13 (PCC = 0.20) for Luminal; and KLF4 (PCC = 0.47) and BDNF (PCC = 0.48) for TNBC. These genes are known to be associated with treatment resistance or breast cancer progression. For example, PAK1 copy number is increased in HER2-positive cases32, and its overexpression is linked to poor prognosis and metastasis33. KLF4 is a known prognostic marker in TNBC34.

For the survival benchmarking analysis, among the four comparisons, including: (1) all samples, (2) HER2-positive, (3) Luminal, and (4) TNBC, the STimage outperforms all the other four models in 1, 2, and 3, while it ranked the second best in the fourth comparison. (Fig. S16b)

To further evaluate the clinical utility of STimage, we independently assessed its ability to predict responses to drugs. Here, we analysed another cohort of breast cancer patients with treatment response information35. Average gene-wise Pearson’s correlation coefficient was found to be 0.47 for 1139 common genes (Fig. S17b). We then evaluated the performance of true and predicted gene expression of 1139 common genes for classifying pathological complete response (pCR) vs residual disease (RD) using an radial basis function kernel support vector machine classifier. Given the different distribution of the data, the predicted gene expression data was log-transformed, while the true raw gene expression data were CPM-normalised and standardized.

We used a two-level stratified cross-validation approach to test how well our model could predict new data. First, we divided the full dataset into five equal folds. Each time, we used four folds to train the model and the remaining fold as a held-out set. Within the four folds, we did another round of stratified 5-fold splitting to get the best-trained model to be used for predicting the 5th held-out set. We repeated this process five times so that every patient was part of the held-out set. This was done to ensure the results are robust and reliable. STimage could stratify patient responses with AUC = 0.70, compared with using the bulk RNA-seq with AUC = 0.78 (Fig. S17c).

Methods
Datasets
The training step requires spatial transcriptomics data as input, with large histology images cropped into 299 × 299 pixel tiles. Spatial coordinates encoded by spatial barcodes for each spot, and gene expression values for each matched spot, are used as labels to train the model. The number of microns corresponding to 299 pixels depends on the image resolution and the size of the Visium spots. The tiling process ensures that each tile contains one spot of 55 μm for 10x Visium and 100 μm for the Visium legacy. Therefore, the 299 × 299 window represents a square tile of 55 μm in dimension for the 10x Visium, and a separate model is trained for this resolution compared to a model trained for a resolution of 100 μm. With this setup, a well-trained model can be applied to histological images (non-ST data) to predict gene expression at the spot level.

First, we focused on breast cancer datasets, including an in-house Xenium single-cell spatial transcriptomics dataset and multiple publicly available breast cancer datasets that differ in resolution (33 samples from the legacy ST platform and nine samples from the 10x Visium platform) and sample preservation types, with two fresh frozen and one Formalin-Fixed Paraffin-Embedded (FFPE) sample among the nine Visium samples. To benchmark STimage, we compared it with three existing methods: STnet9, Hist2Gene10, and Hist2ST11, using the 33 samples from the legacy ST platform9. Model training, evaluation, and testing were conducted using a leave-one-out cross-validation (LOOCV) strategy to benchmark the prediction accuracy of the regression model against the existing models.

To test the generalisability of the model to datasets generated from different laboratories and independent cohorts, we used the nine-sample Visium dataset, including six Visium samples for HER2+ breast cancer15 and three from the 10x Genomics public dataset. The dataset was divided into 70% training data and 30% validation data. The test data included one FFPE sample from 10x Genomics and one fresh frozen sample from the HER2+ dataset15. For independent validation, we randomly selected three whole-slide histological images from the TCGA HER2+ cohort. In this way, the trained models from the STimage software on the Visium dataset were applied to assess their applicability to non-spatial data.

We also produced subcellular-resolution data using Xenium for five samples, and we demonstrated the capability of using this data to validate gene expression predictions as an out-of-distribution (OOD) test. We then further extended the model to predict cell types for hundreds of thousands of individual cells. In all model training scenarios, the data were split at the sample level (WSI), not at the tile level, to avoid overfitting. We also performed OOD evaluations in several different scenarios. The regression model trained on the Visium breast cancer dataset was used to make predictions on H&E images from the legacy ST HER2+ dataset9.

Beyond testing for breast cancer data, we also assessed STimage on two other cancer types (skin and kidney cancer datasets) and a non-cancer liver immune disease dataset (primary sclerosing cholangitis). The in-house skin cancer dataset includes 13 melanoma patients samples, processed using the Visium FFPE protocol. The kidney cancer dataset comprises six Visium samples37, and the liver dataset includes four Visium samples38. LOOCV was used to evaluate the STimage model performance for each dataset separately. To evaluate cross-platform prediction, and on H&E images from in-house PhenoCycler Fusion spatial proteomics data of skin cancer, to compare the predicted gene expression with the protein modality. The detailed dataset information can be found in Supplementary Table S4

Xenium in situ gene expression analysis
The in-house FFPE breast cancer samples were sectioned at 5 μm thickness and mounted on Xenium slides (10x Genomics). In situ gene expression profiling was performed using the Xenium In Situ Platform (10x Genomics) according to the manufacturer’s protocol with the Pre-designed Xenium gene panel targeting 280 genes. In brief, tissue sections underwent deparaffinisation and decrosslinking, followed by probe hybridization overnight. The correct-binding probes were circularised by ligation and amplified by rolling cycle amplification. Autofluorescent quench and nuclei staining were performed prior to Xenium Instrument, software version 1.4.2.0. Raw fluorescence images were processed using Xenium onboard Analysis version 1.4.0.6 with default parameters for cell segmentation, transcript decoding, and assignment.

Skin cancer Visium dataset
The in-house FFPE tissue blocks of 13 melanoma patients samples from 10 patients were sectioned at 5 μm by rotary microtome, and the sections were processed for spatial sequencing library preparation following the Visium Spatial Gene Expression for FFPE User Guide (CG000407, CG000408, CG000409 - 10x Genomics). SpaceRanger V1.0 was used to map FASTQ reads to the human reference genome (version GRCh38-3.0.0) to generate the gene expression count matrix.

Skin cancer Phenocycler Fusion dataset
Single-cell spatial phenotyping of melanoma Formalin-Fixed Paraffin-Embedded (FFPE) slides was performed using the Phenocycler Fusion platform in collaboration with Akoya Biosciences (Marlborough, MA, USA). Tissue sections were stained in a single step with a 33-plex antibody panel. The imaging protocol was established with the Fusion Experiment Designer software according to the manufacturer’s instructions. Following the acquisition of all cycles, a final composite QPTIFF image file was exported from Phenocycler Fusion software and aligned to the H&E image from the same tissue section.

Preprocessing of the input image
Technical variation between H&E images can be caused by various factors, such as different staining procedures, imaging hardware and settings, or operators. Such technical artefacts can affect model training and testing, negatively impacting generalisability. In STimage, we perform stain normalisation for each image, such that the mean R, G, and B channel intensities of the normalised images are similar to those of the template images, while preserving the original colour distribution patterns. STimage uses StainTool V2.1.3 to perform Vahadane39 normalisation as the default option. We show that this processing step improves model performance.

In addition, the nature of tissue sectioning will inevitably result in some tiles containing low tissue coverage. These tiles should be removed. STimage uses OpenCV2 for tissue masking and removes tiles with tissue coverage lower than 70% (default).

STimage regression model
We utilised both sequencing data and H&E image pixel intensity data to train STimage regression models for predicting spot-level gene expression based solely on imaging data.

We adopted the NB-distributed likelihood NB(r, p) for the STimage regression model, where r2NG denotes the “number of successes”, p∈(0,1)G the “probability of each trial’s success”, and G is the number of (assumed) independent genes. The parameters of the NB distribution, r and p, are estimated from the spatial spot images using an image feature extractor (ResNet50 by default), fθ, with learnable pre-trained parameters θ from the ImageNet dataset. The functions hω and hη represent parameter spaces containing the weights and biases of the neural network layers that connect the latent space bz to the two output neurons r and p.

bz= fθ (s)
(1)
br=Softplus(hω (bz))

pb=Softmax(hη (bz)) (2)

The original STimage regression entails the ResNet50 CNN fθ mapping a H&E spot tile, s, to a latent representation bz2Rn features, which is then fed into separate output layers (1) and (2), estimating br and pb, respectively. We denote estimates for the g-th gene’s likelihood parameters by br and pb . Therefore, the STimage regression model defined g-th (random) gene expression X ∣s, as a function of H&E spots, by
g g
X js(cid:2) NB(sjbr ,pb ) (3)
g g g
STimage estimates a NB distribution to model each gene’s (random) expression X ∈ {0, 1, …} conditional on the H&E spot s, i.e.,
μ g (s)=E[X g js]=br g (1(cid:4)pb g )=pb g (4)
g g g
p(X js)=NB(sjbr,pb).

σ g(s)= V[X js]= br (1(cid:4)pb )=pb , (5)
g g g g
where μ (s) in (4) is the mean prediction for the g-th gene’s expression with expectation operator E[:], and σ (s) in (5) is the standard deviation function with variance operator V[:] g .

Overall, the ResNet50 CNN layers are optimized to take H&E images as input to estimate each gene expression’s predictive distribution with NB’s parameters r and p (for all g∈{1,…,G}, Eq. (3)). The ResNet50 base, fθ, is used to perform convolution operations and spatial filtering to extract a three-dimensional feature map (i.e., tensor) from input image tiles. The pretrained fθ of the ResNet50 includes the GlobalAveragePooling reduction operation40 to convert the three-dimensional tensor to the feature vector z, whereby each z corresponds to a single tile/spot within a H&E image.

The fully-connected output layers hω and hη were fine-tuned by maximising the NB’s negative log-likelihood,

XG XB Γ(br )Γ(x +1)
log
Γ
g
(br
,i
+x
g,i
)
(cid:4)br
g,i
log(pb
g,i
)(cid:4)x
g,i
log(1(cid:4)pb
g,i
), (6)
g=1i=1 g,i g,i
with batch size B, Gamma function Γ(. ), and x representing the observed RNA count of the g-th gene (i.e., gene expression).

To assess the model’s predictive capability, we first trained our model on a panel of 14 marker genes, which included three immune marker genes (CD74, CD63, CD81) and eleven cancer-associated genes (COX6C, TP53, PABPC1, GNAS, B2M, SPARC, HSP90AB1, TFF3, ATP1A1, FASN, VEGFA). Two 10x fresh frozen breast cancer samples were used (replicate 1 for training and replicate 2 for testing) to evaluate the model’s performance in biomarker prediction. We then extended the model to nine breast cancer datasets, initially using the top 1000 most variable genes and subsequently applying it to a comprehensive set of the most relevant functional genes (1522 genes). The top 100 genes were selected as the most predictable based on their Pearson’s correlation coefficient (PCC). During training, the model also learned the relationships between these genes.

We further improved STimage by integrating the latest foundation models for image feature extraction. The H&E image accompanying the ST data is preprocessed in the same way as described earlier. Instead of using ResNet50 as the backbone of STimage, we utilise a foundation model to convert H&E tile images into feature embeddings. These embeddings are then connected to STimage’s neural network layers to estimate gene expression distributions. To demonstrate the improvement, we incorporated Virchow2, one of the most advanced foundation models, which has been trained on 3.2 million whole-slide H&E images and contains 632 million parameters. Due to its extensive training dataset and large-scale parameterisation, Virchow2 serves as a more suitable image feature extractor compared to traditional models such as ResNet50. Further, we created an ensemble model named FSTimage, where we integrated multiple foundation models into the STimage base model. Each foundation model contains one foundation model as the upstream image extractor and an STimage negative binomial module. The foundation models used are: GigaPath, H-Optimus, Phikon, cTranspath, CONCH, UNI, and Virchow2. We then aggregated by averaging the results from these models to achieve a final FSTimage ensemble (FSTimage_ens) prediction. We found that STimage, Virchow2_STimage, and FSTimage_ens performed better than any other models that we compared to, including THItoGene and IGI_DL.

Uncertainty assessment
STimage estimates a NB distribution to model each gene’s (random) expression X ∈ {0, 1, …} conditional on the H&E spot s, i.e., X js(cid:2) NB(br g ,pb ). An ensemble of STimage models produces a distribution p(X g ∣s,Θ). Therefore, following the law of total variance, we can obtain two terms accounting for uncertainty:

V[X js]=V[E[X js,Θ]]+E[V[X js,Θ]],
|fflfflfflfflfflfflfflg{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflffl{gzfflfflfflfflfflfflffl} (7)
μ g (s) σ g (s)

where the variance of the mean estimates, V[μ g (s)], defines epistemic uncertainty, and the average of variances, E[σ g (s)], defines aleatoric uncertainty. Aleatoric uncertainty reflects the noise inherent in the data. Epistemic uncertainty reflects the model’s limitations (e.g. lack of information due to small sample sizes or suboptimal model parameters) and can be reduced with more data and/or improved models.

Applying the formula (7), we quantify total uncertainty for the prediction of each gene in each spot. Briefly, we trained 30 STimage models for approximately 1522 functional genes separately, using the same training dataset, each with a unique random seed. Corresponding values of X g were predicted for individual spots s using each trained model with model parameters Θ. Based on the 30 individual STimage models, we computed five ensemble models. Each STimage ensemble model was derived by averaging NB estimates μ g (s) from five randomly pooled STimage single models. For each of the single models and ensemble models, PCC was then calculated as the prediction accuracy metric.

Model interpretation of STimage
STimage implements model explainability in order to highlight important segments/features (nuclei) from the image that contribute to the high expression of a gene. We used a local model-agnostic approach (LIME)41, in which different segments of the images are repeatedly perturbed to measure how the predictions made by the trained model change, in order to score each individual segment.

Morphological diversity among cells in the tumour is an important hallmark that provides clues for histopathological assessment. In order to obtain a score for each nucleus, we used Cellpose-3 segmentation42 to segment the nuclei from H&E images. We performed an image transformation, converting RGB images to Haematoxylin-Eosin-DAB (HED) format before applying the Cellpose segmentation. Inputs required by LIME are a trained regressor/classifier, an image tile, a gene name, and a segmentation function. Using LIME, we identified nuclei segments corresponding to high or low feature importance in the regression model. Given below are the formulas used to compute the LIME values.

ξ(x)=argminL(f,g,π x)+Ω(g),
g2G (8)
where
L(f,g,π x)= π x(z)(f(z)(cid:4)g(z0))2 (9)
z,z02Z

LIME finds a simple model g that can approximate our complex neural network model f. The total loss in LIME is given in Eq. (8). The term L(f,g,π) in (8) is a weighted sum of squared errors (SSE), where the weights are given by π(z), depending on the proximity of the local neighbourhood data point z to x. This ensures that points close to x are given greater weight.

Eq. (9) creates a local approximation of our trained model in the neighbourhood of a given input image x using a simple, interpretable regression model g, where G is a set of sparse linear regression models and π x denotes the kernel function that defines the local neighbourhood around data point x. Here, f(z) is the prediction made by the neural network for features z, and g is a simple model that approximates the behaviour of f in the locality of x, as defined by π. The vector z0 is an explainable subset of z. The second term in (8) is a regularisation term applied to g, which assigns zero weight to non-important input features.

STimage Web Application. To provide an interactive, visual interpretability option, we have developed the STimage WebApp, where users can select image tiles and genes of interest to visualise the nuclei important for the prediction of gene expression. The WebApp allows users to either use a pre-trained model or upload their own trained model for interpretability visualisation. We provide trained STimage regression models for breast cancer, kidney cancer, liver disease, and skin cancer. The Web App can be accessed at: https://gml-stimage-web-app.streamlit.app/.

STimage classification model
We reasoned that categorical prediction of gene expression, mapping whether genes are present or absent, or highly or lowly expressed in specific tissue regions, could already have significant clinical potential.

z= x(cid:4)μ
: (10)
σ

We transformed the raw gene expression values (continuous) to z-scores using (10), and categorised them as “Low” for z-scores below zero and “High” for z-scores above zero. For the classification model, we used a pre-trained ResNet50 (optionally with a fine-tuning layer) as a feature extractor for all tiles, and then applied a regularised logistic regression classifier, optimising the log-loss (binary cross-entropy) function.

yb= 1 , i2f1,...,Ng, (11)
i 1+exp(cid:4)(β0X+b)

We used the “elastic net”43 to balance model performance and complexity by applying both L1 and L2 penalties. In Eq. (11), the coefficient set β=(β,…,β ) and X=(x,…,x ) denote the latent space features (N=2048 features).
1 n 1 n

L logistic (β)= (cid:4) 1 XN XN y i log(p yi )+ (1(cid:4)y i )log(1(cid:4)p yi ), (12)
N i=1 i=1

L enet (β)=L logistic (β)+ω (1(cid:4)α) 2X048 j j 2+α 2X048 (cid:4)β(cid:4) (13)
j=1 j=1

The elastic loss L in Eq. (13) is the sum of the binary cross-entropy loss (log loss) from Eq. (12) and the regularisation terms. Here, α in Eq. (12) is a mixing parameter defining the relative contribution of RIDGE and LASSO to the penalty term. For α = 1, the Elastic Net function reduces to LASSO regularisation (~logistic loss + sum of absolute values of weights), while for α = 0, the function reduces to RIDGE regularisation (~logistic loss + sum of squares of weights). ω is the regularisation parameter, and we set the regularisation strength to ω=0.1 (a smaller value implies weaker regularisation). The solver “SAGA”44 was used for optimisation, as it is faster for large datasets.

The classification model was trained on seven breast cancer samples for the top 100 highly predictable breast cancer genes (from the STimage regression model) and tested on 10x FFPE and fresh frozen HER2+ (1160920F) samples.

Benchmarking and performance assessment
We applied a LOOCV strategy on 33 samples (from a total of 36, with 3 excluded due to low quality) from the HER2+ ST dataset to evaluate STimage regression model performance and benchmark it against existing methods such as STnet9, HisToGene10, and Hist2ST11. Briefly, all models were evaluated using the same data-splitting LOOCV strategy, where models were trained on 32 samples, and the spatial gene expression was predicted on the held-out sample.

For each gene, Pearson’s correlation coefficient (PCC) was calculated between the predicted and experimentally measured spatial gene expression values. PCC was used as the primary performance metric to assess prediction accuracy. We also computed the autocorrelation Moran’s I index as an additional metric, a weighted correlation measure that incorporates the local spatial context.

To assess generalisation, we evaluated model performance on the top 1000 highly variable genes (HVGs). After removing genes expressed in fewer than 10% of total spots, 737 HVGs were retained for further analysis. Each gene received a score from applying the model to a patient, resulting in 737 scores per patient for comparison across models.

To ensure fair comparisons, both the HisToGene10 and Hist2ST11 models were run using the authors’ recommended hyperparameters. However, STnet failed to execute due to technical issues in the original source code, so we re-implemented the STnet model according to the original paper9. Once all predictions from all models were collected, we obtained 97,284 values (33 models × 4 methods × 737 genes) for comparative analysis.

We also benchmarked STimage against the other methods using an out-of-distribution (OOD) test dataset. In this scenario, all models were trained on two fresh frozen breast cancer tissue samples from 10x Genomics and tested on an FFPE tissue sample, which we considered as an OOD sample due to the substantial differences in data generation protocols between these two sample types, poly-A capture for fresh frozen versus probe-capture for FFPE protocols. Using the same gene filtering approach as described above, we obtained 982 of the top 1000 HVGs in both the training and prediction datasets.

Hist2ST and HisToGene include a vision transformer (ViT) component as part of their neural networks, which requires loading all spot tile images and expression values from a single tissue into the model in one batch. Consequently, running these models requires greater computational resources (particularly memory) when applied to Visium data, which contain approximately 10 times more data than legacy ST datasets. Training of Hist2ST and HisToGene failed on a computing cluster equipped with an NVIDIA Tesla V100 SXM2 32 GB GPU, necessitating a workaround: we split the Visium tissue into multiple smaller subsets using a sliding window approach to reduce memory load.

Each subset contained a square tissue region of 25 × 25 spots with no overlap. Each subset was treated as an individual dataset, which reduced the number of spots loaded per batch while preserving the spatial neighbourhood structure. All subsets from the two fresh frozen samples were used to train Hist2ST and HisToGene models using their default hyperparameters. The trained models were then applied to predict gene expression in the FFPE dataset. Model performance was assessed using both Pearson’s correlation coefficient (PCC) and Moran’s I index.

A geometric metric to better assess model predictions of spatial patterns
The commonly used metric in gene expression prediction tasks, PCC, is known to be influenced by outlier noise and is unable to capture spatial patterns. To address this, we developed a metric for evaluating model performance in terms of how well predictions match ground truth spatially. For a given gene, we standardise the predicted and ground truth values to z-scores and then categorise them as binary values, based on whether the values are greater or less than zero. These binary values are then organised into a spatial matrix, creating a binary image where gene expression is plotted according to spot coordinates. To minimise noise and retain reliable prediction patterns, we apply Gaussian filtering with a sigma value of 1 to both rows (spatial X) and columns (spatial Y). After obtaining the refined binarised prediction and ground truth images, we compute the Intersection over Union (IoU) scores to assess the overlap between the spatial patterns of the ground truth and the prediction. This score effectively reflects the accuracy of predicting spatial expression patterns.

Cell type classification model for Xenium single-cell ST data
We modified the Hover-Net model45, a CNN-based instance classification and segmentation model, to enable cell segmentation and cell type prediction for Xenium spatial transcriptomics data at single-cell resolution. First, we generated cell type labels by annotating the cells in the Xenium data with individual cell types using Seurat reference-based integration with the human breast cancer atlas from ref. 15. We merged similar cell type labels into four categories: immune cells, stromal cells, cancer epithelial cells, and endothelial cells. The Xenium labels were then mapped to the H&E image by applying Scale-Invariant Feature Transform (SIFT) registration46 between the Xenium DAPI channel and the H&E image. This enabled us to obtain, for each tissue sample, segmentation labels with unique cell types for nearly every cell in the H&E images (i.e., 1,294,600 cells).

To generate the training and test data, the H&E images were divided into 256 × 256-pixel tiles, and tiles containing fewer than two nuclei were filtered out. The model was trained using four of the five breast cancer tissue samples and evaluated on one held-out test sample. We also calculated the mean IoU score for each LOOCV test sample. We compared our model with the recent CellViT model30, a vision transformer-based model that was trained to predict single-cell segmentation and cell type classification.

Validation using non-spatial data and assessing predictive values for survival and drug responses
To perform external validation, assess robustness, and demonstrate the clinical applicability of STimage, we applied the model to the TCGA dataset. For this analysis, H&E images were downloaded from TCIA47, and bulk gene expression and clinical data for corresponding patients were downloaded from https://www.cancer.gov/ccg/research/genome-sequencing/tcga. We obtained 3066 images from TCIA. A total of 1532 images without matched bulk and clinical data were excluded from the analysis, and the remaining 1534 images were processed for quality control using HistoQC. We filtered out low-quality images using the following parameters: blurry_removed_percent≥0.9, small_tissue_removed_num_regions>97.5th percentile, and pixels_to_use<5th percentile. These thresholds were chosen based on the empirical distribution of quality parameters, removing extreme values. Using the mentioned filtering criteria, 500 poor-quality images were excluded. We finally used 1034 breast cancer images from 670 patients and predicted the gene expression of 1522 predictable genes. Predicted gene expression for each image was averaged across all spots to compare with the ground truth TCGA bulk gene expression. We measured the performance of STimage by calculating Pearson Correlation (PCC) between the average predicted expression across all the spots with the true bulk gene expression for the matched samples. We also computed the performance of STimage for predicting top-300 predictable genes for three subtypes (HER2+, Luminal and TNBC). Patient subtypes were defined according to the criteria described in ref. 48.

STimage predictions were further used to perform survival analysis and stratify patients into low- or high-risk groups. We first built univariate Cox regression models for the top 300 predictable genes using the function ‘coxph’ from the R package survival49. Cox regression models were built separately for true and predicted gene expression. Days to death and days to last follow-up (for patients still alive) were used as survival times, and survival status (alive/dead) was used as the event indicator. The univariate models were evaluated using the concordance index (c-index) for both predicted and true survival outcomes.

The top five genes with the highest c-index were used to build multivariate Cox regression models. We used a 3-fold cross-validation strategy with 100 repetitions to evaluate the multivariate models. The median survival risk score from the multivariate model was used to stratify patients into low- or high-risk groups. Univariate and multivariate models were constructed separately for each subtype. For both bulk and predicted data, the average c-index for the multivariate models was found to be over 0.6 across all subtypes. Model performance was assessed using the ‘concordance’ function from the survival package. Univariate Cox regression models and Kaplan-Meier survival curves were computed using the survival package. Multivariate models were constructed using the R package ClassifyR50.

Moreover, we also compared STimage’s survival performance with the existing gene expression prediction models on TCGA dataset. To conduct the survival benchmarking analysis, we utilised the publicly available code from https://zenodo.org/records/1460248951 and integrated the results of STimage into it. To ensure fair benchmarking, for each method four-fold cross-validation strategy was used to select the best model trained on the Her2ST spatial transcriptomics dataset for the top 1000 HVGs (737 after QC for (gene less than 1000 spots)).

To further assess the clinical utilities, we added to our previous survival analysis an independent assessment of the ability to predict responses to drugs. Here, we analysed a cohort of 168 breast cancer patients treated with neoadjuvant and chemotherapy with (n=65) or without HER2-targeted therapy35. 161 patients were assessed at surgery using the residual cancer burden (RCB) classification. 42 (26%) had a complete response (pCR), 25 (16%) had a good response, 65 (40%) had a moderate response and 29 (18%) had extensive residual disease. We evaluated the performance of classification models using true or predicted gene expression for classifying pCR vs RD using an RBF kernel support vector machine. We used 1139 genes in the functional gene list that were present in this data. The continuous predicted gene expression data were log-transformed, while the raw count gene expression data were CPM-normalised and standardized. We applied a nested 5-fold stratified cross-validation: each outer fold used four folds for training and one for testing, with an inner 5-fold split for model selection. This process ensured every patient was tested once. AUC was calculated for each fold, and the average AUC was reported as the final performance metric.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Tan et al. used very few training and test datasets (training: n=7; test: n = 2) to build a weak model for BRCA, aiming to predict the spatial expression distribution of a series of genes in WSI at the spot level. The logic of the whole manuscript is confusing, and there are the following problems:\n1. The selection of training and test dataset:\na) All samples were derived from only one cancer type: BRCA, and were from the 10X Visium platform. Because of the differences in histological characteristics among cancer types, does this model only predict the gene expression distribution of BRCA WSI samples? The sequencing depth of different sequencing technologies is different, and the lack of ST platform data in the training dataset will cause the problem of domain bias.\nb) The number of train/test datasets is only 9, so it is necessary to increase the number of train datasets for establishing prediction models.\n2. In the Materials and Methods section, the Datasets section mentioned that ”The model applied to this dataset predicts 14 cancer and immune markers and corresponds to all results\". In the STimage regression model, the author indicated that \"We trained our model on a panel of 14 marker genes which include four immune marker genes (CD74, CD24, CD63, CD81) and ten cancer genes (COX6C, TP53, PABPC1, GNAS, B2M, SPARC, HSP90AB1, TFF3, ATP1A1, FASN)\". However, in the legend of Fig. 2e, the author pointed: \"Pearson correlation and Moran's I score of predicted 1000 HVGs against ground truth for each model were compared \", given the above, there are the following questions:\na) Does the CNN-NB model predict the expression distribution of the 14 marker genes listed in the text or the expression of 1000 HVGs?\nb) Why select these 14 marker genes? They are neither classical immunological marker genes nor BRCA-characteristic marker genes, was the selection of these 14 marker genes supported by the papers related to their function in BRCA progression?\nc) Why were PCC of only 12 target genes detected in Figure 1e instead of the 14 listed above? Among them, TTLL12 is not present in the 14 marker genes listed above. The authors need .\n3. In the Preprocessing of the image input section:\na) The author's description is unclear, only indicating the tools of stain normalization and filtering, without clarifying the tiling method.\nb) The authors repeatedly emphasized in the article that the GEP distribution prediction is spot level, how many microns correspond to 299 pixels? Is it the spot (55um) level of 10x Visium? Moreover, this paper involves both 10X Visium and ST data with different spatial resolution data. Can 299 pixels also represent the spot level for ST data?\nc) In Fig S3, the author selected several tile sizes of 224, 299, 450, 600, and 900 to try to explain the reasons for selecting 299 pixels, but firstly, the selection of these sizes is not continuous, and secondly, the spacing among tile sizes is different, with a gap of 75 between 224 and 299, 151 between 299 and 450, 200 between 600 and 900, it is recommended to start with a smaller value and do tile size scale tests every 75 intervals.\n4. Figure 1c-d uses the RNA expression of 10X Visium data as the ground truth, which cannot truly reflect the spatial gene expression. It is recommended to use IF/IHC staining to correlate fluorescence intensity with GEP. For example, the 10X Visium BRCA-IDC dataset, or the IHC staining picture in the TCGA database.\n5. Genes compared with \"ground truth\" in the manuscript are COX6C, KRT5, ATP1A1, S100A14, and C3... S100A14 and C3 are not in the 14 marker genes mentioned above. The comparison results between the predicted values of all 14 marker genes mentioned above and the ground truth are not presented in this manuscript, so it is difficult to avoid the suspicion that the suitable genes are selected for display, which means that the model cannot well predict the spatial distribution of all 14 marker genes. It is recommended that a comparison of ground truth fluorescence intensity at the protein level for all 14 marker genes be added to the supplementary figure.\n6. In the STimage regression model section of the method, the authors claim that \"Our assessment found that the ZINB-based model did not perform as well as the NB-based model\". And \"the Poisson likelihood was also deemed problematic for spatial sequencing data, which is typically over-scattered and sparse.\" In both cases, the authors do not provide data or cite literature to prove their accuracy.\nMinor\n1. In the Materials and Methods, Datasets section, paragraph 3, line 9, the PPFE sample should be FFPE sample.\n2. Materials and methods part, the first paragraph in the Uncertainty assessment section, line 1, Negative Binomial (NB), lack of space. The full name of the Negative Binomial has already appeared once before, there should be an abbreviation.\n3. Fig1 f-h: Fig f-g lacks the annotation of gene name, h lacks the pathologist annotation of tumor/normal region, and all feature plots lack the scale bar.", "reviewer_2": "Summary:\nSpatial transcriptomics is a technique that can predict gene expressions from H&E slides while preserving tissue anatomy. However, this approach is costly, highlighting the necessity of developing an AI-based alternative. Current deep learning methods for this task have limitations in terms of interpretability, robustness, and performance, which are essential for clinical applications. To address this issue, the authors developed a deep-learning probabilistic framework that overcomes these limitations, providing improved robustness and interpretability, even in cases where the training datasets are limited. The developed method's performance was measured against existing methods and found to be superior.\nConcerns:\nThe authors use LIME and SHAP to analyze which regions were important to the predictions. They state that nuclei regions contributed more to the predictions, but this visual analysis was only done for two images. To make such a statement, this visual analysis needs to be done for many more images. Also, this analysis isn’t quantitative, so it’s subject to interpretation and confirmation bias. I recommend removing that sentence altogether since two images aren’t enough to claim this.\nThe caption of Fig 1g states that the predictions are highly correlated with tumor morphology, which again, requires many more samples to make such a statement. Recommend removing.\nSince the interpretability aspect is a core contribution of the paper, it requires more thorough evaluation. Evaluations of explanation faithfulness need to be performed and reported. This may be added as a supplementary table or figure.\nSince no statistical tests were performed, it is difficult to assess the reliability of the results. Where appropriate, statistical tests need to be conducted and their p-values should be reported.\nFig 2f is hard to read. It should be split into two rows for better readability.\nThe text could be worded better and contains many spelling and grammar errors, making it difficult to read. Should be put through a grammar checker tool.\nRecommendation:\nThe method appears to improve upon the currently existing methods in terms of predictive performance. It was also tested on different protocols, labs, and populations. Therefore, I recommend accepting after the concerns are resolved.", "reviewer_3": "In this study the authors present a new algorithm STimage which uses machine learning to predict gene expression on H&E images after learning the spatial gene expression patterns from spatial transcriptomics data. STimage and its performance are demonstrated on breast cancer datasets.\nThe concept is not novel and despite seemingly good performance of STimage, there are from my point of view many parts that are immature analytically which are accompanied with several typos and inconsistency through which make the manuscript difficult to assess.\nMAIN CONCERNS\n- The algorithm is tested only on a very limited breast cancer datasets of breast cancer. The authors should clearly state which data are used for training and testing throughout the manuscript to enable easier comparisons between different sections.\n- There are different ST-seq dataset used across the manuscript Visium 1OX, Visium-FFPE and Visium legacy. While the use of the larger dataset which : Visium legacy is mentioned in the abstract and material and method, it is not seen in the results. The authors should provide more information on how the different ST-seq datasets were used in the study and consider building a training and testing dataset that includes a mix of the three methodologies to better evaluate the algorithm's performance on different types of data.\n- Building on the previous comment in the material and methods the authors claim to use 1 FFPE tissue sample and when coming to the results they report testing 2 tissue FFPE samples. Such inconsistencies across the different sections in terms of how the datasets are used and reported highlight the immaturity of the manuscript.\n- The authors should provide a clearer description of the output of the algorithm, both in the abstract and throughout the manuscript, to ensure that readers understand what the algorithm does and what the results mean.\n- The definition of the 14 marker genes used in the study could be clarified. It is unclear how these genes were selected and whether they represent the top genes predicted by STimage. Additionally, the choice of these genes is not well-justified, and it is not clear why no stromal or endothelial markers were included. Providing more context on the selection criteria and the rationale behind these choices would help readers better understand the approach. Moreover, it could be useful to consider incorporating stromal and endothelial markers as these compartments are important in breast tumors and can also be annotated by pathologists. This would improve the algorithm's ability to predict gene expression in all compartments of breast tumors.\n- The Swarbrick’s lab dataset has been annotated in the original publication, could this be used for validation, interpretation here?\n- The non spatial tile labelling and testing is not convincing, first it is not clearly explained how the training was performed. Second testing is performed only on two samples which are closely related since they are successive sections from the same tumor. One needs to really extend the testing to prove applicability.\n- The preprocessing, processing use of the ST-seq data is not described.\n- It is unclear what is the impact of OpenCV2 on tiles removal , especially the tiles situated in the edges that may contain less than 70% tissue coverage, it seems to me that some tiles in the 'middle' of the tissue with bad coverage may be removed, but also tiles on the edges or tiles containing adipose or less dense tissue? How many tiles are removed per image and where? Are important tiles removed? Is it necessary to remove them? Are these tiles creating so much noise that they need to be removed?\n- Stain normalization seems to be performed to look like the RGB intensity of a template image. What is the influence of the choice of the template image. How was the template image chosen?\n- The uncertainties are assessed spatially. Have the authors assessed if the uncertainties seem to be mainly driven by spatial location or the genes assessed? In other words, what is the highest source of uncertainty, the quality of the image? Or the capacity of the algorithm to predict a specific gene expression?\n- Breast cancer is a very heterogeneous disease, when assessing only 14 genes across the manuscript it may be that the performances are different from subtype to subtype. No track of the clinicopathological features of the sample used is given or hint towards the model being associated with subtypes is given. Some have shown that subtype can be predicted from H&E images, it is therefore important to assess how this can affect the model.\n- The interpretability analysis in the end of the ‘prediction performance’ paragraph is disappointing, the only feature which the author mention in their interpretability section is nuclei which has been reported before, the author attempt to conclude on other feature in a gene-specific manner… perhaps… but I am not managing to understand the last sentence which involve COX6C and KRT5 genes. It would be interesting to understand further which feature drive the prediction of different genes\n- Can uncertainties by training on different ST-seq platform/ methodologies be assessed?\n- What is the scalability of STimage, can more genes be assessed? It is important to provide the performance of STimage for all genes tested at least the 14 mention all along ideally all genes.\n- The interpretability, last chapter, is not clear, suddenly in this chapter a new gene for immune marker appear C3 is it a typo? do the author mean CD63? or the C3 gene is used, C3 is not a typically marker gene for immune cells.\nMINOR comments\n- No p values are shown for the PCC (pearson correlation)\n- Typo: for indepenedent validation (p2)\n- Typo: PPFE sample (p2)\n- Typo: A classification model can be added top produce complementary information (p8)"}, "author_response": "We appreciate the reviewer's comment. In our original submission, we utilized the largest and most diverse collection of data available, which included HER2+ low-resolution samples (n=23 patients, 36 tissues) and three datasets with higher-resolution Visium data (n=6 patients, n=2 for Fresh Frozen samples, and n=1 for an FFPE sample). Our results demonstrated that our model not only outperformed three existing models in accuracy but also introduced significant improvements in robustness and interpretability, which were the primary focus of our work.\nWe observed that with the rich content provided by spatial transcriptomics data, a neural network can be effectively trained on a relatively small dataset. For instance, the STnet model (He et al., 2020, Nature Biomedical Engineering) was developed using the same 36 samples, comprising a total of 13,620 spots (legacy ST). In spatial transcriptomics, each sample includes thousands of individual spots, each corresponding to an image tile with thousands of gene expression labels. In our original submission, the nine Visium samples spanned three datasets: 16,456 spots (Wu et al.), 7,785 spots (10x Genomics, Fresh Frozen), and 2,518 spots (10x Genomics FFPE).\nIn this revised version, to further evaluate robustness and increase sample sizes, we incorporated additional datasets, including kidney cancer (Fresh Frozen, n=6, 15,632 spots), non-cancer liver disease (FFPE, n=4, 19,967 spots), and skin cancer (FFPE, n=10, 14,405 spots). Additionally, we applied our model to the latest subcellular resolution Xenium spatial transcriptomics data (n=5). This subcellular data enabled us to predict cell types at single-cell resolution and assess the generalization of the model trained on low-resolution data when applied to single-cell spatial data. Furthermore, we validated our models using a non-spatial TCGA dataset (n=450). To ensure robustness and validity, we compared our predictions with spatial proteomics data (CODEX). We believe our work represents the most comprehensive assessment of different spatial datasets for this type of model.\nWe apologize for any lack of clarity in the logical flow of the manuscript. In response to the reviewer's concerns, we have extensively revised the paper as detailed point-by-point below.\nThis work addresses a significant challenge in digital pathology by leveraging the rich content of spatial transcriptomics data while tackling the high cost barrier of this technology. We developed machine learning models that use H&E image inputs to predict gene marker expression and cell types, both crucial tasks in digital pathology. Unlike existing methods, our approach focuses on two key aspects often overlooked: robustness and interpretability. We demonstrate that our model performs well across various sample types, disease types, resolutions, and diverse cohorts, including non-spatial data. Robustness was achieved by predicting data distribution, estimating uncertainty, and employing an ensemble approach, while interpretability was ensured through feature attribution methods (LIME), latent space explanation, and integration with pathological annotations and the biological functions of the genes. Additionally, we applied our model to breast cancer survival prediction, showcasing its potential clinical benefits.\n1. The selection of training and test dataset:\na) All samples were derived from only one cancer type: BRCA, and were from the 10X Visium platform. Because of the differences in histological characteristics among cancer types, does this model only predict the gene expression distribution of BRCA WSI samples? The sequencing depth of different sequencing technologies is different, and the lack of ST platform data in the training dataset will cause the problem of domain bias.\nWe appreciate the reviewer's question. Addressing potential domain bias is a key focus of this work, and we have thoroughly evaluated the model's performance on both in-distribution and out-of-distribution samples. Our model was trained using various protocols, including the 10x Visium poly-A capture protocol (offering comprehensive gene coverage), the probe-capture protocol (more sensitive), the legacy ST protocol (lower resolution), and the latest Xenium platform (single-cell resolution). The 10x spatial platform is unique in generating both H&E images and gene expression data by default, which is crucial for training imaging-to-expression models. Additionally, we tested the model trained on the spatial BRCA dataset by predicting gene expression in non-spatial H&E images from the TCGA dataset (n=5).\nIn the revision, we extended our model's evaluation to multiple cancer types, including kidney cancer (Fresh Frozen, n=6, 15,632 spots), non-cancer liver disease (FFPE, n=4, 19,967 spots), and skin cancer (FFPE, n=10, 14,405 spots). The extensive analyses have been added to the Supplementary Figures S6, S7, S8.\nWe also conducted a comprehensive assessment of the model's robustness using multiple training and testing strategies. Overall, we found no significant domain bias when the ST platform was excluded from the training dataset. In addition to the previously applied leave-one-out strategies, we evaluated the model on both in-distribution and out-of-distribution datasets, considering variations in resolution and different cancer types. We further analyzed which genes were most predictable by the model.\nb) The number of train/test datasets is only 9, so it is necessary to increase the number of train datasets for establishing prediction models.\nAs mentioned above, we observed that with the rich content provided by spatial transcriptomics data, a neural network can be effectively trained on a relatively small dataset. For instance, the STnet model (He et al., 2020, Nature Biomedical Engineering) was developed using the same 36 samples, comprising a total of 13,620 spots (legacy ST). In spatial transcriptomics, each sample includes thousands of individual spots, each corresponding to an image tile with thousands of gene expression labels. In our original submission, the nine Visium samples spanned three datasets: 16,456 spots (Wu et al.), 7,785 spots (10x Genomics, Fresh Frozen), and 2,518 spots (10x Genomics FFPE).\nIn this revised version, to further evaluate robustness and increase sample sizes, we incorporated additional datasets, including kidney cancer (Fresh Frozen, n=6, 15,632 spots), non-cancer liver disease (FFPE, n=4, 19,967 spots), and skin cancer (FFPE, n=10, 14,405 spots). Additionally, we applied our model to the latest subcellular resolution Xenium spatial transcriptomics data (n=5). This subcellular data enabled us to predict cell types at single-cell resolution and assess the generalization of the model trained on low-resolution data when applied to single-cell spatial data. Furthermore, we validated our models using a non-spatial TCGA dataset (n=450). To ensure robustness and validity, we compared our predictions with spatial proteomics data (CODEX). We believe our work represents the most comprehensive assessment of different spatial datasets for this type of model.\n2. In the Materials and Methods section, the Datasets section mentioned that ”The model applied to this dataset predicts 14 cancer and immune markers and corresponds to all results\". In the STimage regression model, the author indicated that \"We trained our model on a panel of 14 marker genes which include four immune marker genes (CD74, CD24, CD63, CD81) and ten cancer genes (COX6C, TP53, PABPC1, GNAS, B2M, SPARC, HSP90AB1, TFF3, ATP1A1, FASN)\". However, in the legend of Fig. 2e, the author pointed: \"Pearson correlation and Moran's I score of predicted 1000 HVGs against ground truth for each model were compared \", given the above, there are the following questions:\na) Does the CNN-NB model predict the expression distribution of the 14 marker genes listed in the text or the expression of 1000 HVGs?\nWe apologize for the lack of clarity. We trained two models separately: one with 1,000 HVGs and another with the 14 genes. This has been clarified in the current version. Additionally, we included an analysis of the most predictable genes among approximately 1,522 biologically relevant genes, representing cell type markers, immuno-oncology genes, and others. Our findings suggest that this approach can identify highly predictable gene sets that are biologically related to the cancer types and the organs of origin.\nb) Why select these 14 marker genes? They are neither classical immunological marker genes nor BRCA-characteristic marker genes, was the selection of these 14 marker genes supported by the papers related to their function in BRCA progression?\nWe selected the 14 markers to demonstrate the model's performance on a panel of genes, providing a relevant analysis for exploring potential clinical applications of predicting gene panels rather than individual gene markers. These genes include a set of cancer and immune marker genes identified as spatially variable in the Visium training dataset.\nSpecifically, the genes are COX6C (cancer metabolic processes, such as the Warburg effect in cancer glycolysis), SPARC (extracellular matrix regulation, associated with cancer metastasis and progression), HSP90AB1 (association with angiogenesis and proliferation, a potential therapeutic target), TFF3 (associated with breast cancer proliferation, hormone regulation and survival), ATP1A1 (cell proliferation, survival, adhesion, migration, and interaction with the tumor microenvironment in breast cancer), CD74 (enhance invasiveness and interact with macrophage migration inhibitory factor) and CD81 (involved in cell adhesion, migration, and signal transduction, cancer proliferation, migration and metastasis). The remaining genes are CD24 (regulate cell adhesion and migration; high CD24 is associated with increased metastasis, poorer prognosis), CD63 (may act as both tumour promoter or suppressor depending on the context), TP53 (a well established cancer marker), PABPC1 (globally manage mRNA, may contribute to uncontrolled cell growth), GNAS (contribute to cell signaling and metastasis), B2M (essential for MHC-I, can act as either tumour suppressor or promoter), FASN (upregulated in breast cancer, associated with cancer metabolism, implicated in resistance to cancer chemotherapy).\nIt is important to note that the purpose of testing this gene set is to illustrate the performance of the model with a gene panel, while STimage can be trained on any gene set. In the revised version, we conducted a thorough analysis of a large set of 1,522 known gene markers, including those related to cell types, immune pathways, oncology pathways, and other essential pathological processes. This approach contrasts with the common analysis of the 1,000 most variable genes, many of which may not be relevant to cancer. Based on the prediction results, we selected the genes with the highest predictability, as indicated by their Pearson correlation scores. These top predictive genes appear to be specific to cancer types and are enriched for the organs from which the cancer originates, such as kidney cancer.\nc) Why were PCC of only 12 target genes detected in Figure 1e instead of the 14 listed above? Among them, TTLL12 is not present in the 14 marker genes listed above. The authors need .\nWe apologize for the inconsistency as among the 14 markers, two genes were not present in the training datasets. In our revised version, we have extended the analyses to thousands of genes, beyond the initial 14 markers. The 14 genes include seven known markers (COX6C, SPARC, HSP90AB1, TFF3, ATP1A1, CD74, and CD81) and seven genes used in STnet (CD24, CD63, TP53, PABPC1, GNAS, B2M, and FASN) (He et al., 2020). However, since two of these genes could not be removed from the test dataset, we replaced them with two new genes, including the TTLL12. In the revised version, we have shifted our focus from the top 14 marker genes to the hundreds of most predictable and functionally relevant genes.\n3. In the Preprocessing of the image input section:\na) The author's description is unclear, only indicating the tools of stain normalization and filtering, without clarifying the tiling method.\nWe apologize for not describing in detail the tiling method. Briefly, the spatial x and y coordinates for each spot were used to determine the center of each tile. Depending on the WSI resolution, the cropping size uses pixels/µm to convert the pixel values for 100 µm (legacy ST) and 55 µm (Visium). The cropped tiles were then resized to 299x299 to ensure that tiles from different samples with different resolutions have comparable area coverage. We have added this description to the method session.\nb) The authors repeatedly emphasized in the article that the GEP distribution prediction is spot level, how many microns correspond to 299 pixels? Is it the spot (55um) level of 10x Visium? Moreover, this paper involves both 10X Visium and ST data with different spatial resolution data. Can 299 pixels also represent the spot level for ST data?\nAs described in the session above, the number of microns corresponding to 299 pixels are dependent on the image resolution and the size of the Visium spots. The tiling process ensures that each tile contains one spot of 55um for 10x Visium and 100um for the Visium legacy. Therefore, the 299x299 window represents a spot of 55 micron, a squared tile of 55 micron in dimension; a model was trained for this resolution separately from another model that was trained for a resolution of 100 micron.\nAs described above, the number of microns corresponding to 299 pixels depends on the image resolution and the size of the Visium spots. The tiling process ensures that each tile contains one spot: 55 µm for 10x Visium and 100 µm for Visium legacy. Thus, a 299x299 pixel window represents a 55 µm or 100 µm square tile. A model was trained specifically for a 55 µm resolution, separate from another model trained for the 100 µm resolution.\nc) In Fig S3, the author selected several tile sizes of 224, 299, 450, 600, and 900 to try to explain the reasons for selecting 299 pixels, but firstly, the selection of these sizes is not continuous, and secondly, the spacing among tile sizes is different, with a gap of 75 between 224 and 299, 151 between 299 and 450, 200 between 600 and 900, it is recommended to start with a smaller value and do tile size scale tests every 75 intervals.\nWe thank the reviewer for the suggestion. It is noted that 224 and 299 are two commonly used image input sizes for neural networks. As suggested by the reviewer, we have added a 750 tile size to ensure 150 intervals from 450 to 900. There is no change in our conclusions and it remains that 299 x 299 is still the best. We have added the results to the Supplementary Figure 4 (Supplementary Figure S4).\n4. Figure 1c-d uses the RNA expression of 10X Visium data as the ground truth, which cannot truly reflect the spatial gene expression. It is recommended to use IF/IHC staining to correlate fluorescence intensity with GEP. For example, the 10X Visium BRCA-IDC dataset, or the IHC staining picture in the TCGA database.\nOur model is designed to predict gene expression, not protein expression. Therefore, the ground truth for our predictions is the gene expression values obtained from experimental spatial transcriptomics technology. Other ground truths considered in this work include bulk gene expression data (for TCGA data) and pathological annotations for comparing cell type clustering results based on predicted gene expression values with those based on pathological annotations. We evaluated the clustering results using predictions versus the original spatial transcriptomics data, relative to the annotation ground truth (Supplementary Figure S5).\nFurthermore, we utilized the extensive non-spatial H&E data from the TCGA database, where bulk RNA-seq data and survival data serve as the ground truth (Figure 1d, Figure 5, Supplementary Figure S9).\nFollowing the reviewer’s suggestions, we also tested a case where we used protein expression from CODEX data as a “ground truth” to compare to predicted gene expression values based on H&E of the same section (Supplementary Figure S11). We wish to emphasize that our models are to predict gene expression, not protein expression. It is established from the literature that, most of the time, the correlation between these two modalities is mostly not high. Therefore, we expected that the performance will be lower for the proteins that are not well correlated with gene expression. Our results, nevertheless, show high consistency between two modalities. We added the following figure to the manuscript (Also shown below).\nFollowing the reviewer’s suggestions, we also tested a case where we compared predicted gene expression values based on H&E images with protein expression from CODEX data as a “ground truth” (Supplementary Figure S11). We emphasize that our models are designed to predict gene expression, not protein expression. Literature indicates that the correlation between these two modalities is often low (Pearson correlation ranges from 0.53 to 0.71, PMID:27951527). Therefore, we anticipated that performance might be lower for proteins not well correlated with gene expression. However, our results demonstrate high consistency between the two modalities. We have included this additional figure in the manuscript (also shown below).\nSupplementary Figure S11.\n5. Genes compared with \"ground truth\" in the manuscript are COX6C, KRT5, ATP1A1, S100A14, and C3... S100A14 and C3 are not in the 14 marker genes mentioned above. The comparison results between the predicted values of all 14 marker genes mentioned above and the ground truth are not presented in this manuscript, so it is difficult to avoid the suspicion that the suitable genes are selected for display, which means that the model cannot well predict the spatial distribution of all 14 marker genes. It is recommended that a comparison of ground truth fluorescence intensity at the protein level for all 14 marker genes be added to the supplementary figure.\nAs mentioned in response to comment #4, the use of 14 gene markers was intended to demonstrate how the model performs with a gene panel, similar to those used in current sequencing-based diagnostics. The selection of these 14 genes was based on literature rather than our experiments. In the revised version, we moved beyond the 14 genes and expanded our analysis to 1,522 well-curated functional genes. We present the performance of all 1,000 highly variable genes (Figure 1, Figure S1) and the top 100 most predictable genes (Figure S1).\nFor the task of predicting gene RNA expression, the ground truth is the measured RNA levels, not protein expression. Previous research indicates that protein levels are generally not highly correlated with RNA levels, with Pearson correlations ranging from 0.53 to 0.71 (PMID:27951527). Nonetheless, we assessed the predicted RNA expression against measured protein expression levels using spatial proteomics CODEX data. These new results are included in Figure 3d and Figure S11. Additionally, we have added multiple independent evaluation analyses to assess model performance. For example, we compared predictions with out-of-distribution expression data for all 106 shared genes in single-cell resolution Xenium data (Figure S10) and for all 1,146 shared genes in lower-resolution legacy spatial transcriptomics data (Figure S9). We also compared predictions for all 1,522 genes with bulk RNA sequencing data from 450 patients in the TCGA dataset and with independent pathological annotations (Figure 5).\n6. In the STimage regression model section of the method, the authors claim that \"Our assessment found that the ZINB-based model did not perform as well as the NB-based model\". And \"the Poisson likelihood was also deemed problematic for spatial sequencing data, which is typically over-scattered and sparse.\" In both cases, the authors do not provide data or cite literature to prove their accuracy.\nWe provided analysis results for the ZINB-based model performance here. We have also added results for Poisson likelihood for comparison.\nMinor\nWe thank the reviewer again for the comments and suggestions.\n1. In the Materials and Methods, Datasets section, paragraph 3, line 9, the PPFE sample should be FFPE sample.\nAmended\n2. Materials and methods part, the first paragraph in the Uncertainty assessment section, line 1, Negative Binomial (NB), lack of space. The full name of the Negative Binomial has already appeared once before, there should be an abbreviation.\nAmended\n3. Fig1 f-h: Fig f-g lacks the annotation of gene name, h lacks the pathologist annotation of tumor/normal region, and all feature plots lack the scale bar.\nAmended"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Although the author has done more analysis on the original version, the performance of the model needs to be further improved. Some similar works in the past two years seem to have better performance than STimage, such as THItoGene and IGI_DL. The following are the additional comments for the author's modification:\n1. The process of selecting the variable genes and 1,522 functional genes is unclear. Please provide a clear explanation of how these genes were chosen, including the specific methods and thresholds utilized. The rationale for selecting the highly variable genes (HVGs) and the genes used for prognosis prediction should be clearly articulated.\n2. The rationale for selecting only 450 samples for the prognosis prediction analysis, rather than using the full TCGA cohort, is unclear. If any sample selection was performed, the criteria and justification must be provided. Please present the results for all TCGA samples.\n3. The correlation analysis presented was conducted on a per-sample basis. Please provide the results of the correlation between the predicted gene values and the true bulk gene expression values for each gene. Are the genes with a high correlation to bulk values the ones that your gene prediction model predicts with greater accuracy?\n4. The presence of outliers in the data may affect the accuracy of the Pearson correlation coefficient. The Spearman correlation coefficient results require further explanation to interpret result significance in the context of your study.\n5. Moran's I is used to analyzing the spatial distribution pattern of a single variable. However, the author uses it to assess the accuracy of gene predictions without comparing it to the Moran index of actual gene values, which is confusing.\n6. The details of the prognosis prediction model are missing, particularly the method used to derive the survival score from the predicted gene expression values. More thorough explanation is needed.\n7. The poorer performance of the STImage model in some cases of the Her2ST LOOCV benchmark is concerning. Whether STImage is biased to the spatial pattern of a specific slice?\n8. While multiple models were trained using different datasets, the benchmark analysis was limited to the Her2ST dataset and 785 genes instead of 1522 genes that used in prediction. This raises questions about the ability to accurately reflect the models' predictive performance and generalization capabilities.\n9. The authors only benchmarked HisTogene, STnet, and Hist2ST released in 2021 or 2022. How do the more recently released THItoGene and IGI_DL perform?\n10. The median PCC of BRCA FFPE ST data predicted by Stimage is less than 0.1 in Fig. 2E, this suggested a very weak correlation.\n11. In Figure S8, the PCCs around different slides show quite different results, with some slides showing median PCCs below 0.\n12. IoU values in Fig.S11g are abnormally high. From the figure, it seems there aren't many overlapping areas between prediction and ground truth.\n13. The section on Xenium cell type classification does not present the results of the leave-out test comprehensively. It lacks comparative results with other classification models and evaluation metrics, such as mIoU.\n14. Is there a cell type marker in the genes predicted by the model? This can be compared with the classification results.\n15. The authors didn’t assess prognosis performance for STimage, such as C-index.", "reviewer_3": "Thank you to the authors for addressing the questions and comments previously raised. Although I appreciate the technical efforts to improve robustness, reproducibility, and interpretability, the manuscript remains highly technical without demonstrating the biological application of the developed tool suite. The tool suite's applicability appears complex and not easily accessible to researchers.\nThe main concerns are the lack of easy access and the absence of user-friendly software, as well as insufficient biological interpretations, which limit the study's interest. While the authors highlighted key novel aspects and technical improvements in their response, these improvements should be benchmarked against other algorithms. More critically, there are no new biological findings that enhance our understanding of diseases or physiological processes, which raises serious questions about the usefulness of the presented tool suite and algorithm.\nIn the response, the authors mentioned inter and intra-tumor heterogeneity of cancer. However, this aspect is critically missing in the manuscript. There are no measurements of intra-tumor heterogeneity or its association with gene expression and/or clinically relevant parameters. The algorithm's ability to detect phenotypically different tumor clones is neither assessed nor validated. Intra-tumor heterogeneity could be compared to matched WGS for cancer cell heterogeneity or scRNA-seq data for immune cell heterogeneity.\nAdditionally, with the few samples analyzed for each disease, it is probably not relevant to discuss inter-tumor heterogeneity. The authors also mention \"highly relevant predictable genes for each cancer type,\" but it is unclear in the manuscript what these genes predict or how they reveal relevant biological insights.\nIn short, while I appreciate the technical improvements and the additional data used in the revised version, the manuscript still lacks critical elements such as cell type specificity, clinically relevant analysis, spatial ecosystem generalization, and specific spatial gene expression. These missing aspects result in a lack of novel biological insights, which limits the overall interest and impact of the study.", "reviewer_4": "All comments of the reviewer 2 has been addressed."}, "author_response": "We have compared STimage with THItoGene and IGI_DL. In both cases, STimage shows superior performance (Figure R1). Notably, in this revised version, we have markedly improved the STimage performance by integrating state-of-the-art foundation models into STimage. Foundation Models, such as Virchow2, were trained on 3.2 million whole slide H&E images, with 632 million parameters, and therefore are more suitable image extractors than traditional methods, such as ResNet50. Our advanced Virchow2_STimage model outperformed all other models. Further, we created an ensemble model named FSTimage_ens, where we integrated multiple foundation models into the STimage base model. Each foundation STimage model contains one foundation model as the upstream image extractor and an STimage negative binomial module. The foundation models used are: GigaPath, H-Optimus, Phikon, cTranspath, CONCH, UNI, and Virchow2. We then aggregated by averaging the results from these models to achieve a final FSTimage ensemble (FSTimage_ens) prediction. We found STimage, Virchow2_STimage, and FSTimage_ens performed better than any other models that we compared to, including THItoGene and IGI_DL. We introduce the FSTimage in our revised Figure 1b.\nOur new results are shown in the figure below (and have been added to the main Figure 2e in the revised version).\nFigure R1. Benchmarking STimage with IGI-DL, THItoGene and other models. Introducing two new models Virchow2_STimage and FSTimage_ens.\n1. The process of selecting the variable genes and 1,522 functional genes is unclear. Please provide a clear explanation of how these genes were chosen, including the specific methods and thresholds utilized. The rationale for selecting the highly variable genes (HVGs) and the genes used for prognosis prediction should be clearly articulated.\nWe consider a set of all 1630 functional genes (1,522 remained after removing lowly expressed genes) combined from panel gene lists curated by main spatial transcriptomics platforms available for Xenium and CosMX. These genes are the union of 1253 genes from the human pan-cancer panel and 1056 genes from the human immunology panel (10xGenomics). The pan-cancer panel contains: 1)21 cancer-related pathways such as Hedgehog/Ras/p53/MAPK/Notch/mTOR/TNF/TGF-beta/TLR/Wnt signalling pathways;2)markers for 28 cancer cell/tissue types, such as colon cancer, breast cancer; 3) 17 cancer-related cellular processes such as immune response, chemotaxis; and 4) 37 functional annotation(e.g., nucleus, extracellular matrix) and processes (e.g., angiogenesis, mitosis). The immunology panel contains: 1) 20 immune pathways, 2) markers for 29 cell/tissue types (e.g., liver, lung, NK cells, macrophage, DC), 3) 41 immune-related cellular processes (e.g., inflammatory responses, innate immune responses), 4) 25 immune-related functional processes. Although these lists are not exhaustive, the genes represent an unbiased selection of functionally relevant genes, not just focusing on 1000 highly variable genes as in previous reports.\n2. The rationale for selecting only 450 samples for the prognosis prediction analysis, rather than using the full TCGA cohort, is unclear. If any sample selection was performed, the criteria and justification must be provided. Please present the results for all TCGA samples.\nWe initially performed random sampling to analyse a subset of the TCGA dataset. To address the reviewer’s request, we downloaded all 3,066 images for TCGA breast cancer and applied the pre-trained STimage model to all samples. The conclusion from this analysis remains the same compared to the initial random subset analysis (also refer to responses to comment 3, and Figure 4b).\nH&E images that didn’t have matched bulk RNAseq data or clinical data were excluded. We used HistoQC to remove low quality images from the analysis, removing 215 images using filtering criteria below:\n● blurry_removed_percent ≥ 0.9\n● small_tissue_removed_num_regions > 97.5th percentile\n● pixels_to_use < 5th percentile\nAll images that passed the quality filtering criteria were used for further analysis. The results show that our predicted data can significantly stratify overall survival across the three subtypes of breast cancer (Figure R2).\nThe sample-wise correlation was performed for the samples that had matched bulk RNA-seq gene expression and H&E tissue. The gene-wise correlation was performed for the samples with matched bulk RNA-seq expression data, H&E tissue and breast cancer subtype information (three breast cancer subtypes, Her2, Lum and TNBC), (Figure R2). Survival analysis was performed for the samples with matched H&E, bulk RNA-seq gene expression, and survival data.\nFigure R2. Kaplan Meier survival curves for all patients in each subtype. The first row is the results using the measured bulk RNAseq data and the second row is for the predicted data from just H&E (at spot levels and then averaged across all spots to get one value per gene per image). The genes used for plotting these curves are top 5 genes with the highest c-index from the univariate cox-regression analysis.\n3. The correlation analysis presented was conducted on a per-sample basis. Please provide the results of the correlation between the predicted gene values and the true bulk gene expression values for each gene. Are the genes with a high correlation to bulk values the ones that your gene prediction model predicts with greater accuracy?\nWe performed analysis on a per-sample basis for the whole transcriptome rather than per gene across samples because, in the TCGA dataset, the bulk sequencing data and the image data for a patient often come from different tissue blocks of the same patient. The two tissue blocks from a breast cancer patient can be markedly different and so the information encoded by the bulk sequencing and the image data could be relatively large. Given a gene, changes in the expression value from one patient to another patient was measured by RNAseq data, but not by imaging data. The deviation of the correlation value between RNAseq data and image-based prediction is therefore the composite of three sources of variation: between patients, between imaging sections and tissue blocks used for sequencing, and true difference between predicted values of the genes from the imaging sections and the sequencing results of the same imaging section (which does not exist). Therefore, the correlation for one gene across samples is affected by multiple factors not just how accurately the model predicts the gene expression in an image, but also affected by the variation between patients in terms of image quality, the technical sequencing variation across samples, and biological variation between tissues used for imaging and different tissue blocks used for sequencing. Therefore, the sample-wise correlation becomes a metric that is not possible to be decomposed to assess the model performance. Even the gene-wise correlation for each sample is less ideal, but is a better metric to use.\nWe also re-run the gene-wise correlation analysis on a per-sample basis for all TCGA samples and the results are shown in the FigR4 (added to the revised version as Figure 5). The average correlation for the new results is similar to previous results i.e., 0.48.\nFigure R4. Sample-wise and Gene-wise Pearson correlation analysis using the TCGA breast cancer dataset. This figure shows the gene-wise correlation for each sample and across 1034 samples. The number on the top right indicates the number of samples (each dot in the plot represents a sample).\n4. The presence of outliers in the data may affect the accuracy of the Pearson correlation coefficient. The Spearman correlation coefficient results require further explanation to interpret result significance in the context of your study.\nWe agree with the reviewer and therefore we added more assessment metrics, including IoU as discussed in more details below.\n5. Moran's I is used to analyzing the spatial distribution pattern of a single variable. However, the author uses it to assess the accuracy of gene predictions without comparing it to the Moran index of actual gene values, which is confusing.\nHere, Moran’s I was used to assess how the spatial autocorrelation in the predicted values correlates with that of the original spatial transcriptomics data. A gene with a higher Moran’s I value suggests that the spatial pattern is better preserved than a gene with a lower Moran’s I. This is an additional metric to PCC, a commonly used metric that does not account for the spatial context.\nThe interpretation of spatial autocorrelation is as follows: if the score is high, it means that the prediction matches the ground truth and that the prediction of the gene across the spots also preserve the spatial pattern of that gene between neighbor spots. If the score is low, it could mean either that the prediction is inaccurate or that the gene does not have a spatial pattern.\n6. The details of the prognosis prediction model are missing, particularly the method used to derive the survival score from the predicted gene expression values. More thorough explanation is needed.\nWe have performed more thorough analysis as described in our updated text of the Methods section as below:\nTo perform external validation, assess the robustness and demonstrate clinical application of STimage, we applied the model on the TCGA dataset. For this analysis, H&E images were downloaded from TCIA and bulk gene expression and clinical data for corresponding patients were downloaded from “https://www.cancer.gov/ccg/research/genome-sequencing/tcga”. We downloaded 3,066 images. 1,532 images that didn't have matched bulk and clinical data were excluded from the analysis and remaining 1,534 images were used for quality check using HistoQC. For 285 images HistoQC results were NA or not generated as the images were corrupted and couldn’t be read. For the remaining 1,249 images with HistoQC results, we filtered out low quality images using the following parameters:\nblurry_removed_percent ≥ 0.9, small_tissue_removed_num_regions > 97.5th percentile and pixels_to_use < 5th percentile. These thresholds were chosen based on the empirical distribution of quality parameters, removing extreme values. Using these filters, 215 images were excluded. Only images that passed all the filtering criteria were then used for further analysis. We finally used 1,034 breast cancer images from 670 patients and predicted the gene expression of 1,522 highly variable genes.\nPredicted gene expression for each image was averaged across all spots to compare with the ground truth TCGA bulk gene expression. We measured the performance of Stimage by calculating Pearson Correlation (PCC) between the average predicted expression across all the spots with the true bulk gene expression for the matched samples. We also computed the performance of STimage for predicting top-300 predictable genes for three subtypes (HER2+, Lumnial and TNBC). Patient subtypes were defined according to the criteria described in Xu et. al. (2022).\nSTimage prediction was then used to perform survival analysis and stratify the patients into low- or high-risk groups. We first built a univariate cox-regression model for the top 300 predictable genes using the function ‘coxph’ from R package “survival”. Using the top-300 predictable genes, we built the univariate Cox regression model for each gene individually. Cox regression models were built separately for true and predicted gene expression. Days to death and days to last follow-up (for patients still alive) were used as survival times, and survival status (alive/dead) was used as the event indicator. The univariate models were evaluated using the concordance index (c-index) for both predicted and true survival outcomes. The top-5 genes that had the highest c-index were then used to build the multivariate cox-regression model. We used a 3-fold cross validation strategy with 100 repetitions to evaluate the multivariate models. The median survival risk score from the multivariate model was used to stratify patients into low- or high-risk groups. Univariate and multivariate models were built separately for each subtype. For the bulk and the predicted data, the average c-index was found to be 0.75 for all subtypes. We used the function ‘concordance’ from the R package “survival” to assess the concordance performance. Univariate Cox regression models and Kaplan-Meier curves were computed using R package “survival”. The multivariate regression models were built using the R package “ClassifyR”.\n7. The poorer performance of the STImage model in some cases of the Her2ST LOOCV benchmark is concerning. Whether STImage is biased to the spatial pattern of a specific slice?\nSTimage performed better than other methods in the legacy dataset (Her2ST data) in most cases. We found that the performance was only equal or slightly lower in the samples that the quality of the gene expression matrices was spare (i.e., too many 0s values across the whole matrix). Notably, for the sparse data, the spatial feature appears to be less reliable, even with the ground truth data, suggesting that the comparisons for those samples with low quality should be interpreted with caution, or should even be excluded. We found that the STimage performance is consistently better for less sparse samples.\nImportantly, our new FSTimage_ens performed better than all other models in all cases, including the ones with over-sparse data, although we suggest that for sparse samples, they should not be used for assessing models.\nFigure R5. The PCC for all models on the Her2ST dataset in LOOCV settings. The PCC for each test sample is grouped and ordered by the sparsity of the real measurement of ST data. The sparsity is calculated by the number of zeros in the ST gene expression data, ordered from high sparsity (left) to low sparsity (right). Refer to the Figure R1 above for a clearer presentation of the Figure (not ordered by sparsity).\n8. While multiple models were trained using different datasets, the benchmark analysis was limited to the Her2ST dataset and 785 genes instead of 1522 genes that used in prediction. This raises questions about the ability to accurately reflect the models' predictive performance and generalization capabilities.\nWe appreciate the reviewer’s comment regarding the benchmarking analysis. To address this, we have added additional benchmark experiments across 10 samples from 5 Visium slides covering other cancer types, including skin melanoma datasets, to further evaluate the predictive performance and generalization capabilities of our model. The 785 genes mentioned by the reviewer were from our test for the top 1000 highly variable genes (HVG), not for the 1522 functional genes. The reason for using 758 genes, but not 1000 HCG, was that we removed the genes that expressed in less than 10% of the total spots, an unbiased way to remove noise. In addition, we introduced the analysis of 1522 functional genes.\nThe results demonstrate that STimage consistently outperformed the models suggested by the reviewer, including THItoGene and LGL_DL (Figure R6). Moreover, our new model FSTimage_ens more markedly outperformed all other models (Refer to Figure R1).\nFor the generalization, we assessed the models for different cancer types, including skin cancer, kidney cancer, and a non-cancer liver dataset and across different experimental protocols, patient cohorts and data modalities (Table 1 below).\nFigure R6. PCC for benchmarking STimage, THItoGene, and IGL_DL on 10 melanoma samples across 10 Visium slides. Each group of boxes shows the PCC for the top 300 predictable genes for each model across test Visium slides.\nTable 1. Datasets used for evaluation of STimage performance.\n9. The authors only benchmarked HisTogene, STnet, and Hist2ST released in 2021 or 2022. How do the more recently released THItoGene and IGI_DL perform?\nAs discussed above in the Comments #7 and #8, and the general comment with Figure R1, we have added new benchmarking analyses on HER2ST and the Lagrange ST dataset for breast cancer, as well as Visium data for skin cancer. We consistently found higher performance for STimage baseline across datasets. Moreover, our FSTimage_ens markedly improves the performance of the STimage baseline model.\n10. The median PCC of BRCA FFPE ST data predicted by Stimage is less than 0.1 in Fig. 2E, this suggested a very weak correlation.\nThe low performance for this particular setting was expected, as we trained the model on a small Visium fresh-frozen dataset (2 samples) and tested it on the FFPE dataset (1 sample). This experiment evaluated the model's performance on an out-of-distribution (OOD) test, as the image and gene expression data quality differed significantly between the fresh-frozen and FFPE datasets. Therefore, given this stringent setting to test both OOD and the small-size training datasets, a weak correlation is expected. In most applications, the training dataset is often more than 10 slides, a sample size that consistently resulted in predictable genes with PCC ranges from 0.25 to 0.8 (Figure 2c, e). Moreover, PCC is not the accurate metric to reflect model performance regarding how well the spatial pattern can be predicted. In this case, intersection of union (IoU) is a more relevant metric. We show that we can achieve IoU consistently higher than 0.7 (Figure 2d and Figure S4b).\n11. In Figure S8, the PCCs around different slides show quite different results, with some slides showing median PCCs below 0.\nPlease refer to our responses to comment #7. The low performance for these samples was due to poor data quality in some samples (sparse data with too many 0 values). In Figure S8a and S8b, we show that even with a very low PCC value, the prediction of cancer markers is still consistent with the pathological annotation and, in some cases, even better than the ground truth spatial transcriptomics measurements, possibly due to the model also learning from imaging data not just from the sparse sequencing data.\n12. IoU values in Fig.S11g are abnormally high. From the figure, it seems there aren't many overlapping areas between prediction and ground truth.\nWe thank the reviewer for identifying this issue and we apologise for this unintentional mistake. The high IoU scores for classification model were due to applying the same IoU function from the regression model to the binary classification model. The parameter ‘average’ in the function “jaccard_score” from sklearn module was set to “binary” for the classification model rather than “None” for the regression model, where the predicted values are continuous. We have now fixed that issue. The correction has been made for the figure S11 and the IoU remains relatively high (0.3 to 0.6).\nFigure R7. Recomputed IoU scores for the classification model. The Jaccard scores were obtained for both “High” and “Low” class for a gene (high and low refer to the binary classes of the categorised gene expression). The average of scores of both classes were computed as a final IoU score of a gene. An example of the gene CD52 for the two test-set samples is shown. We have corrected the IoU scores in the Figure S11 panels C and G.\n13. The section on Xenium cell type classification does not present the results of the leave-out test comprehensively. It lacks comparative results with other classification models and evaluation metrics, such as mIoU.\nWe performed a leave-out test using five Xenium datasets, and the results are presented in Figure R8 (we added this Figure as the Figure S14 in the revised version). The performance remained consistent across all test samples, demonstrating the robustness of our classification model. Among the different cell types, cancer epithelial cells were the most accurately predicted (82.0%-92.6% accuracy). The second most predictable group was immune cells, which included B cells, T cells, and myeloid cells (58.6%-83.1%). Ability to predict cancer cells and immune cells from H&E images enables the prediction of cancer-immune interactions within the tumour microenvironments, suggesting diagnosis and prognosis values. Given that much more data will be available in the future, our model has the potential to achieve even better predictive performance, particularly for rarer cell types.\nWe also calculated the mean IoU score for each LOOCV test sample (Figure R9; added as Figure S15 in the revised version). In benchmarking against the newly published CellViT (Hörst et al. 2024) model, which is a Vision Transformer-based model that can also predict single-cell segmentation and cell type classification, our performance in terms of IoU score is 2 to 4 times more accurate than CellViT (Figure R10; added as Figure S15 in the revised version). Notably, our classification model utilises Xenium data from a small set of samples, different to CellViT that utilised hundreds of thousands of images. With our approach, the training data is generated automatically using an image registration method to align the DAPI cell segmentation mask (which is more accurate than HE) with cell types derived from marker gene expression (experimental measurement) from Xenium data. This approach reduces the need for manual labelling and minimises human error.\nFigure R8. Leave-one-out cross-validation for single-cell segmentation and cell type classification in the breast cancer Xenium dataset. Each row represents one test sample in each cross-validation experiment. The first column shows the predicted cell type, the next column contains the ground truth cell type labels from the Xenium data, and the next row presents the HE image, which is the only data required in the prediction step. The last column shows the confusion matrix for all five cell types. The numbers indicate the number of cells, and the percentages represent the proportion of each prediction relative to all ground truth values.\nFigure R9. Mean IoU score for each sample in the LOOCV experiment. The blue bar represents CellViT, while the orange bar represents STimage Classification.\nFigure R10. Randomly selected predictions at the tile level to compare the performance of the STimage classification model, CellViT. Each tile has three columns: from left to right, the STimage classification model, CellViT, and the ground truth.\n14. Is there a cell type marker in the genes predicted by the model? This can be compared with the classification results.\nWe thank the reviewer for their comment. We performed cell type prediction on one of the held-out Xenium test samples to demonstrate that the model can accurately identify cell types (Figure R11a, b). Using the same H&E image, we applied the STimage gene expression prediction model that was trained on the low-resolution Her2st dataset to predict marker expression in the Xenium sample (Figure R11d). CD79A, an immune cell marker, and ESR1, a breast cancer marker, were selected for visualisation. The predicted expression patterns of CD79A and ESR1 closely align with the corresponding immune and cancer cell types predicted by the classification model. This consistency across models suggests that the predictions are biologically meaningful. Importantly, the classification and regression models were trained on different datasets, yet both produced coherent and spatially consistent results in for same HE image. This highlights the robustness of the models and indicates their ability to generalise and reflect real biological signals.\nFigure R11. Comparison of cell type prediction from the classification model at the single-cell level and marker gene prediction from the regression model at the spot level. (a) Predicted cell types for one leave-out test Xenium sample. (b) Ground truth labels annotated from Xenium gene expression data. (c) H&E image for the same slide used as input for prediction. (d) Predicted markers for CD79A and ESR1 at the spot level using the model trained on the Her2st dataset. (e) Pathological annotation for the same slide.\n15. The authors didn’t assess prognosis performance for STimage, such as C-index.\nWe thank the reviewer for the suggestion to assess the prognostic performance of STimage. We have calculated the c-index for true and predicted data for both univariate and multivariate models. We performed the survival analysis for each subtype separately. The top-5 genes for each subtype resulted in a better performance following the multivariate analysis; giving consistent results with a c-index of 0.75 (Figure R12; added to the revised version as Figure 5d).\nUnivariate Analysis: Using the top-300 predictable genes, we run the univariate cox-regression model for each gene individually.\nMultivariate Analysis: Using the top-5 genes with highest c-index obtained from univariate analysis, we built a multivariate cox-regression model using the function CrossValidate from ClassifyR to predict the risk/survival score.\nFigure R12. We assessed the prognostic performance of the STimage predictions for TCGA data as well as the bulk data for the three subtypes (Her2, Lum and TNBC). Each dot in the univariate analysis box-plot represents the c-index of a gene. For the multivariate analysis, we used a 3-fold cross-validation strategy with 100 repetitions to evaluate the performance. Each dot in the multivariate analysis is a c-index of the test set. The red dotted line indicates the mark of 0.5 c-index."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have improved the study, the reviewer still has following comments based on their revised manuscript.\n1. The authors suggested the relatively low Pearson correlation coefficients (PCC) observed in several benchmark datasets due to insufficient training sample size. However, they proceed to perform downstream analyses using these potentially inaccurate predictions. As the task of inferring gene expression levels from H&E images is inherently a regression problem, the use of Intersection over Union (IoU) — which binarizes continuous gene expression values — may not appropriately reflect the accuracy of regression outputs. It is strongly recommended that the authors report additional regression metrics such as Mean Absolute Error (MAE) or Mean Squared Error (MSE) to provide a more comprehensive and quantitative assessment of prediction accuracy.\n2. This study predominantly highlights genes with strong prediction performance, which appear to coincide with spatial patterns of H&E staining intensity. This raises the concern of potential model bias, specifically, whether STimage disproportionately favors genes with highly localized expression patterns that align with histological contrast. It would be important for the authors to clarify whether STimage struggles with genes exhibiting broadly distributed expression patterns. Furthermore, is there a correspondence between high Pearson correlation coefficients and high spatial overlap of predicted versus true gene expression domains? For genes with low PCC, is inaccuracy primarily due to incorrect prediction of spatial distribution or due to inaccurate expression value estimation?\n3. Multiple previous models incorporate spatial relationships among tissue tiles either through architectural design or as part of the input features, yet STimage appears to omit such explicit spatial encoding. Despite this, the model achieves better predictive performance across multiple benchmarks. It would be valuable for the authors to discuss the implications of this design choice. Does the model implicitly capture spatial dependencies through other means, or is such information not critical for the tasks addressed?\n4. This study lacks a direct comparison between STimage's survival prediction performance and existing H&E-based models, such as IGIDL survival model. The authors primarily compare their own univariate and multivariate models, which limits the reader’s ability to contextualize the improvements achieved. Furthermore, the rationale for performing subtype-specific survival modeling in breast cancer cohorts is not clearly articulated. Would a pan-subtype survival model (without stratifying by subtype) result in inferior performance, and if so, could this be quantified and discussed?", "reviewer_3": "I appreciate the thorough responses provided by the authors addressing my concerns regarding the biological applications and user accessibility of the STimage tool suite.\nOverall, I am satisfied with the responses and the revisions made by the authors."}, "author_response": "We appreciate the reviewer’s suggestion for adding more regression metrics for each dataset. Previous works such as STNet and Hist2ST use MAE or MSE as loss functions to optimise model parameters. Thus, these cannot be used as performance metrics for benchmarking. In line with the reviewer’s recommendation, we will report MAE and MSE for all datasets as additional regression metrics (see Figure R1 and the accompanying supplementary CSV file). These metrics evaluate point-wise accuracy within each sample but do not capture gene-wise trends across spots within a sample. The PCC provides a reasonable measure for this aspect and is widely used in previous works. However, PCC, MAE and MSE do not account for spatial context, which is central to our use case. We therefore further assess the IoU, which quantifies the overlap of high-expression domains. One of the reasons we use IoU is that this metric is intuitively more similar to tissue assessment by pathologists, who classify tissue regions in a manner akin to IoU rather than at a single tissue location. MWe also report Moran’s I (see Figure S2), which measures concordance in spatial autocorrelation between predicted and observed gene expression.\nFigure R1. Comparison of regression performance metrics (MAE, MSE, PCC) across multiple spatial transcriptomics datasets. Boxplots show mean absolute error (MAE), mean squared error (MSE), and Pearson correlation coefficient (PCC) for predicted versus observed gene expression values across samples within each dataset. Panels correspond to (top row) HER2ST breast cancer samples (Legacy ST), (second row) skin melanoma samples (Visium), (third row) kidney samples (Visium), and (bottom row) liver samples (Visium).\n2. This study predominantly highlights genes with strong prediction performance, which appear to coincide with spatial patterns of H&E staining intensity. This raises the concern of potential model bias, specifically, whether STimage disproportionately favors genes with highly localized expression patterns that align with histological contrast. It would be important for the authors to clarify whether STimage struggles with genes exhibiting broadly distributed expression patterns. Furthermore, is there a correspondence between high Pearson correlation coefficients and high spatial overlap of predicted versus true gene expression domains? For genes with low PCC, is inaccuracy primarily due to incorrect prediction of spatial distribution or due to inaccurate expression value estimation?\nWe appreciate the reviewer’s concern about potential bias towards genes with high histological contrast. We evaluated STimage on three gene sets chosen to be independent of H&E intensity patterns: (i) functional genes drawn from immuno-oncology and cell-type marker panels that were selected for biological relevance rather than morphology; (ii) the top 1,000 highly variable genes defined by expression variance in the matrix, not by spatial variability; and (iii) a clinically used 14-gene panel.\nTo address the concern regarding potential bias towards spatial patterns of H&E staining intensity, we first examined the genes with high prediction accuracy and low correlation with H&E image intensity (Figure R2). The results show that the model is not learning H&E intensity-related features for those highly predictable genes.\nFigure R2. examining relationship between gene expression predictions and H&E image intensity, along with spatial expression maps for genes in the HER2ST dataset. For each gene with high PCC and low correlation with H&E image intensity (CCL19, HLA-DPA1, JCHAIN, LYZ), the left panels show scatter plots of predicted expression versus median RGB intensity per tile, annotated with evaluation metrics (IoU, MAE, MSE and PCC). The middle and right panels show spatial expression plots for predicted and ground truth expression, respectively. Warmer colours indicate higher expression levels.\nWe further performed a more systematic and comprehensive analysis of the relationship between prediction accuracy and the correlation of both true and predicted gene expression with median H&E intensity across four datasets (her2st, skin, kidney, liver), evaluating four metrics (MAE, MSE, PCC, IoU), as shown in Figure R3. We observed dataset-specific trends. In the HER2ST and liver datasets, prediction accuracy tends to be higher for genes with moderate to high correlation with H&E intensity. This could indicate that H&E intensity dominates the histomorphological features. Conversely, for the skin and kidney datasets, prediction performance appears largely independent of H&E intensity correlation, indicating that STimage is not systematically biased towards HE intensity-linked genes. For genes with low PCC, this can be due to both inaccurate spatial distribution and inaccurate expression-value estimation.\nWe also observe the key marker gene KRT5 for myoepithelial cells in Figure 1C. The predicted pattern matches the ground truth, is expressed in the basal cell layer of DCIS, and shows no correlation with H&E image intensity. It is well recognised that not all genes show distinct spatial structures or a strong relationship to histological morphology, which inherently constrains their predictability.\nFigure R3. Relationship between prediction accuracy, gene expression, and histology image intensity across datasets for top 100 most predictable genes. Each scatter plot shows the correlation of predicted gene expression with the median H&E intensity (y-axis) versus the correlation of true gene expression with the median H&E intensity (x-axis) for each top predictable genes in four datasets (HER2ST, skin, kidney, and liver). Point colours indicate different regression metrics(MAE, MSE, PCC and IoU).\n3. Multiple previous models incorporate spatial relationships among tissue tiles either through architectural design or as part of the input features, yet STimage appears to omit such explicit spatial encoding. Despite this, the model achieves better predictive performance across multiple benchmarks. It would be valuable for the authors to discuss the implications of this design choice. Does the model implicitly capture spatial dependencies through other means, or is such information not critical for the tasks addressed?\nWe thank the reviewer for this thoughtful comment. We agree that spatial relationships are important. In the gene expression prediction task, gene expression values are spot-level aggregates at relatively coarse scales (for example, ~55 µm in Visium and ~100 µm in legacy ST). Much of the morphologic context and local cell organisation relevant to these gene expression values are already contained within each spot-level image tile. So explicit cross-spots graphs can mix signals over large distances that may not correspond to the target spot and can introduce noise unless very large training sets are available. STimage therefore prioritises robust estimation by modelling gene counts with a negative binomial likelihood rather than using a fixed point estimate, and by employing deep ensembles, which we demonstrate improves performance and stability across benchmarks.\nAlthough we do not impose an explicit spot-level graph to capture the nearby spatial information, in the updated FSTimage variant, we replace image features with foundation-model embeddings (ViT based) trained self-supervised on millions of whole slide images; these embeddings are known to encode multi-scale contextual structure and further improve robustness.\n4. This study lacks a direct comparison between STimage's survival prediction performance and existing H&E-based models, such as IGIDL survival model. The authors primarily compare their own univariate and multivariate models, which limits the reader’s ability to contextualize the improvements achieved. Furthermore, the rationale for performing subtype-specific survival modeling in breast cancer cohorts is not clearly articulated. Would a pan-subtype survival model (without stratifying by subtype) result in inferior performance, and if so, could this be quantified and discussed?\nWe thank the reviewer for this constructive comment. The IGI-DL survival model implements a graph-based deep learning approach that takes as input both IGI-DL–predicted gene expression and patient-level clinical variables (e.g., age, race, gender, AJCC pathological stage). The scope of this model is therefore broader than ours, which intentionally evaluates survival prediction performance based solely on predicted gene expression to isolate and quantify the prognostic value of gene expression features derived from histopathology.\nMoreover, we attempted to reproduce the IGI-DL survival model following the authors’ GitHub repository. While we were able to perform the initial steps (whole slide image pre-processing, patch tiling, stain colour normalisation, and nuclei segmentation), we encountered reproducibility issues at the graph construction stage. Specifically, the pipeline is built for Visium spatial transcriptomics data and relies on Visium spot coordinates, log transform gene expression matrix and hard-coded parameters (e.g., the neighbour filtering distance threshold is 20 µm, corresponding to 0.5 µm per pixel for graph construction. The Visium pipeline uses 55 µm tiling (matching the Visium spot size), whereas TCGA uses 100 µm tiling). These constraints prevented us from benchmarking with the IGI-DL model on the TCGA dataset.\nInstead, we benchmarked four models, representing the three most widely used architectures: CNN-based (ST-Net), ViT-based (HisToGene, Hist2ST), and capsule-net/graph-based (THIToGene). To ensure fair benchmarking, we used four-fold cross-validation for each method to select the best model trained on the Her2ST spatial transcriptomics dataset for the top 1,000 HVGs (737 remaining after quality control, i.e., genes expressed in fewer than 1,000 spots were removed). The model was then used to make predictions on 1,253 TCGA samples (all of TCGA-BRCA with survival information and bulk RNAseq data). Among the four comparisons, including: 1) all samples, 2) HER2-positive, 3) Lumnimal, and 4) TNBC, the STimage outperforms the all the other four models in 1, 2, and 3, while it ranked the second in the 4 comparison (Figure R4). Notably, in all cases, STimage can still stratify patients into high- and low-risk groups across (Figure R4).\nRegarding the reviewer’s question on breast cancer subtype-specific modelling, we stratified the survival models by molecular subtypes, which reflect the well-established prognostic categories to display information that is more clinically relevant for each subtype (ER-positive, HER2-positive and triple-negative). A pan-subtype model would aggregate biologically distinct cohorts, potentially diluting subtype-specific prognostic signals. To address the reviewer’s request, we have added a pan-subtype model for comparison and we still observed significant stratification between high- and low-risk groups for the combined category (Figure R4). We have updated the main Figure X accordingly.\nFig R4. Benchmarking of survival analysis. a) Kaplan-Meier survival curves for true and predicted gene expression from different models for all samples and three subtypes of breast cancer using the predicted gene expression from STimage. b) C-index comparison between different models for all samples and three subtypes of breast cancer.\nWe have also updated the figure 5 in the manuscript, to have a pan-subtype survival model."}, {"round": 4, "reviewer_comments": {"reviewer_1": "The authors have significantly improved the manuscript. The authors did not perform a comparison between the STimage and IGI-DL, but they provided explanations for this decision. I am satisfied with the responses and the revisions made by the authors. Congrats."}, "author_response": ""}]

Decision letter