Model Architecture and Analysis Pipeline
SpaCEy is an explainable method that detects spatial patterns in tissue samples which are linked to clinical outcomes. The model input consists of protein-abundance measurements at single-cell level represented as spatial graphs; clinical outcomes are used as supervised labels, while clinical annotations are only retained for descriptive and post hoc interpretation (Methods, Data preprocessing and Training settings) (Fig. 1a). Cell-type and compartment annotations are not used as SpaCEy node features, but are retained as metadata for descriptive and post hoc biological interpretation. The spatially resolved samples are transformed into a spatial graph by Delaunay triangulation method24. This approach ensures that each cell is connected to its spatially nearest neighbours, preserving the tissue architecture in the graph structure. In these spatial graphs, each node represents a cell, and the features of each node represent the abundance of each measured protein marker in the corresponding spatial location. A GNN model is then trained on this spatial graph representation, leveraging the connectivity and spatial relationships between cells to learn meaningful latent representations of the input samples (Methods, Model architecture) (Fig. 1a). To interpret the learned model, an explainer model generates feature and edge masks that highlight the most relevant components contributing to the GNN’s predictions for the clinical outcome of interest (Methods, Explainer model) (Fig. 1b)23. The input to the explainer model consists of the trained model and constructed spatial graph, and its output is the learned edge masks (i.e., edge importance values) for the input graph. A key novelty of SpaCEy lies in its ability to generate interpretable spatially contiguous regions (i.e., connected subgraphs) that uncover spatially localised molecular drivers of clinical outcomes. Node importance values are computed by aggregating edge-mask values from k-hop neighbours of each node, enabling the identification of key spatially contiguous, compact subgraphs that mark important regions in the tissue (Methods, Node importance and important regions). Here, node and edge importance scores produced by the explainer quantify unsigned relevance magnitudes. Thus, they identify the graph components that are most influential for the model prediction, but they do not by themselves encode the direction of effect, i.e., whether a given node, edge, or spatial motif shifts the prediction toward or away from a clinical outcome of interest. Using these importance values together with the latent representations from the trained model, we perform unsupervised analysis of the latent space representations of samples, and perform differential analysis to extract biologically relevant insights (Fig. 1c). Based on the identified molecular and spatial features, SpaCEy provides a data-driven framework for understanding the heterogeneity of structure and function in the tumour microenvironment and enables an improved patient stratification (Fig. 1c).
Fig. 1: SpaCEy model architecture and analysis workflow.
a Samples from multiple cohorts were collected as the source of high-dimensional spatial data for this study. These samples were converted into spatial graphs using the Delaunay triangulation method, and a graph neural network (GNN) was trained on these graphs, integrating spatial information and clinical data to predict relevant clinical outcomes. b The trained GNN model is then analysed using the explainer model, which outputs learned edge masks. Node importance values are computed by aggregating edge-mask values within the k-hop neighbourhood of each node, enabling the identification of compact connected subgraphs that define important spatial regions in the tissue. c Latent representations of samples and identified important regions were then used for visualisation, differential analysis, and biological interpretation based on node importance values, ultimately enabling patient stratification and downstream analysis. The female icon in panel a was obtained from Servier Medical Art (https://smart.servier.com/) and is licensed under CC BY 4.0. The breast icon in panel a was obtained from Bioicons (https://bioicons.com/?query=breast) and is licensed under CC BY 3.0 Unported. Source data are provided as a Source Data file.
SpaCEy identifies spatially organised protein and cellular patterns predictive of tumour progression
We first applied SpaCEy to a recently published lung adenocarcinoma (LUAD) imaging mass cytometry data by Sorin et al.25 comprising samples from 416 patients. In the LUAD cohort, we modelled progression as a supervised binary classification task using the progression annotations provided with the dataset. Summary statistics and cohort characteristics of LUAD dataset is provided in Supplementary Table 1. Fig. 2a shows Kaplan–Meier estimates of overall survival stratified by progression status, with censoring marks indicated on each curve, revealing no significant difference in overall survival between the two progression groups (log-rank p = 0.0767). To investigate the relationship between the clinical outcomes and the latent representations generated by our model, we visualised the embeddings across all nodes for each sample in UMAP space (Fig. 2b), with each point representing a sample and coloured by progression status. Together, these visualisations demonstrate that the learned embedding is associated with clinical progression.
Fig. 2: Progression-associated spatial and cellular patterns in lung adenocarcinoma.
a Kaplan–Meier estimates of overall survival for patients stratified by progression status, with censoring marks shown on each curve. Survival distributions were compared using a non-directional omnibus log-rank test. No adjustment for multiple comparisons was applied because a single global comparison was performed. b A two-dimensional UMAP projection of sample-level embeddings generated from spatial proteomic features, showing the organisation of samples in a low-dimensional manifold based on the clinical outcome of patients c. Differential abundance analysis identifying distinct protein markers between progressor and non-progressor patients. d Cell type proportions of important nodes in the LUAD dataset (top panel) and significant proportional differences between progression and no-progression groups (bottom panel). Positive values indicate enrichment in the progression group, while negative values indicate enrichment in the no-progression group. e Representative progressor LUAD sample showing important cancer and endothelial cell nodes. The zoomed region highlights direct edges in the input spatial graph between important cancer and endothelial cells f. Comparison of SpaCEy with the state-of-the-art methods. Individual points represent results from each cross-validation fold (n = 5 folds per method). Bars represent the arithmetic mean. Error bars are centred on the arithmetic mean and extend ± 1 sample standard deviation across folds. The same fixed, progression-stratified folds were used for all methods. Source data are provided as a Source Data file.
Analysis of marker abundance in the important nodes revealed distinct immune microenvironments and tumour-associated profiles between patients with and without progression Fig. 2c. The T-cell labels follow the source annotations from Sorin et al.25: Th cells refers to CD4-enriched helper T cells, Tc refers to CD8a-enriched cytotoxic T cells, Treg refers to FOXP3-enriched regulatory T cells, and T other denotes T cells not assigned to these specific subsets. Here, important nodes refer to a group of cells selected by the explanation module as part of model-informative spatial regions. In contrast, non-progressing tumours exhibited areas with elevated numbers of Th cells, Tc cells, and B cells. Specifically, FOXP3, CD163, Histone H3, and DNA1 were enriched in the important nodes of progression group. FOXP3 is a canonical regulatory T-cell marker, and CD163+ macrophages form part of an immunosuppressive axis previously associated with aggressive LUAD architecture and poor outcomes25. In contrast, non-progressing patients exhibited areas with higher expression of HLA-DR and CD4, and a slight upregulation of epithelial cancer cell marker pan-cytokeratin. HLA-DR marks primarily antigen-presenting cells that are part of an active anti-tumour immune response26,27. In the original Sorin et al. study, Th and Tc populations showed stronger cancer-cell interactions in lower-grade tumours than in high-grade solid LUAD25. Therefore, SpaCEy specifically highlights these regions of cancer-immune cell interactions that are characteristic of tumours with better prognosis. Analysis of cell-type proportions revealed differences between progression and non-progression groups Fig. 2d which shows cell-type enrichment among important nodes. Tumours from patients who progressed showed areas with significantly increased frequencies of macrophage populations, including CD163+ macrophages, as well as higher cancer and endothelial cell content. This is consistent with the original publication demonstrating that CD163+ tumour-associated macrophages, typically linked to immunosupression and tumour progression, are enriched in high-grade tumours and strongly co-occur with FOXP3+ regulatory T cells25. The increased co-localisation and interaction between cancer and endothelial cells observed in higher-grade histological patterns reflects their aggressive nature and high metastatic potential25. To provide an image-level example, we selected a representative progressor LUAD sample and visualised direct graph contacts between explainer-important cancer and endothelial cell nodes, highlighted with bold edges (Fig. 2e). This example illustrates local cancer–endothelial adjacency within an explainer-highlighted tissue region. Interestingly, SpaCEy also identifies a previously not described enrichment of CD163− macrophages in regions associated with progression. In contrast, non-progressing tumours exhibited areas with elevated numbers of Th cells, Tc cells and B cells. To assess whether this immune-cell signal was located in direct proximity to tumour cells, we defined tumour adjacency as a direct graph edge between an important immune cell and at least one cancer cell. Using this definition, 49.9% of pooled important Tc cells (3,019 of 6,053 cells) and 40.9% of pooled important Th cells (4,416 of 10,803 cells) in the no-progression group were cancer-adjacent, indicating that this no-progression-associated T-cell signal is partly tumour-localised. For comparison with the entire Tc and Th cell populations in the same no-progression group, 46.8% of all Tc cells (38,074 of 81,330 cells) and 38.6% of all Th cells (48,467 of 125,557 cells) in the no-progression group were cancer-adjacent. Thus, important Tc and Th cells were only modestly more tumour-adjacent than the corresponding global Tc/Th populations. CD4, CD8a, and FOXP3 are measured marker channels, whereas Th, Tc, Treg, and T other are source cell-type annotation labels shown in Fig. 2d. Tc and Th cells play key roles in anti-tumour immunity, whereas B cells have previously been associated with improved overall survival in LUAD25. Overall, these patterns indicate that SpaCEy identifies clinically relevant regions predictive of tumour progression and also shows potential to reveal emergent spatial trends.
We evaluated SpaCEy against two state-of-the-art baselines, Ali et al.14 and SPACE-GM12, using five-fold cross-validation across all performance metrics (accuracy, F1-score, and AUC) Fig. 2f. We performed a random split stratified by disease status and fixed the splits to ensure a fair comparison and an equal distribution of disease status across folds. SpaCEy achieved the highest overall performance and the most stable behaviour across folds. In terms of accuracy, SpaCEy reached a mean of 0.68, outperforming both Ali et al. (0.61) and SPACE-GM (0.55). Similar trends were observed for F1-score, where SpaCEy got 0.68 on average, whereas the method from Ali et al. and SPACE-GM achieved 0.59 and 0.52, respectively. SpaCEy also demonstrated higher discriminative capability, yielding the highest mean AUC (0.62), compared with 0.57 for Ali et al. and 0.52 for SPACE-GM. Notably, SpaCEy consistently performed strongly across folds, including peak values of 0.77 accuracy, 0.75 F1-score, and 0.78 AUC.
SpaCEy stratifies breast cancer patients by subtype-agnostic localised spatial patterns
The second dataset used to train SpaCEy was an imaging mass cytometry dataset comprising 720 samples from two patient cohorts at the University Hospital of Basel and the University Hospital of Zurich, which includes breast cancer samples across a range of clinical subtypes and cancer stages28. We refer to this dataset as the JacksonFischer dataset throughout the remainder of the text. Here, our objective was to model overall survival (hereafter referred to as survival) in a censoring-aware time-to-event setting. Summary statistics and cohort characteristics are provided in Supplementary Table 2. For the JacksonFischer dataset, Kaplan–Meier curves (Fig. 3a) stratified by clinical subtype (HR+HER2+, HR+HER2-, HR-HER2+, and TripleNeg) showed significant differences among clinical subtypes in overall survival (log-rank p = 0.012), where patients with triple-negative breast cancer (TNBC) exhibit the lowest survival rates due to the aggressive nature of the subtype, lack of targeted therapies, and higher likelihood of early metastasis29,30. The corresponding UMAP of the embeddings learned by the SpaCEy is shown in (Fig. 3b). The embeddings revealed a structured distribution of samples, with similarity grouping associated primarily with survival.
Fig. 3: Survival-associated spatial patterns in the JacksonFischer breast cancer cohort.
a Kaplan-Meier estimates of overall survival across clinical subtypes in the JacksonFischer cohort, with censoring marks shown on each curve and a global log-rank test b UMAP visualisation of sample embeddings derived from the trained model, showing overall survivability distributions across clusters and highlighting potential cluster-specific survivability trends. c Differential abundance analysis identifying distinct protein markers between clusters, with Cluster 0 representing high-survival patients and Cluster 2 representing low-survival patients. A representative lower-survival patient is shown alongside the analysis, illustrating overlap between regions identified as important by the model and local abundance of Slug and c-Myc proteins. d Cell type proportions of important nodes in the JacksonFischer dataset and proportional differences between Cluster 2 and Cluster 0. Positive values indicate enrichment in Cluster 2, while negative values indicate enrichment in Cluster 0, with enrichment scores shown for each cell type. e Distribution of important regions identified by the model. f Cell type annotations one representative sample: first column shows overall annotation, second highlights tumour core, peritumoural, and extratumoural regions, and third indicates important regions. g UMAP visualisation of latent representations for HR+HER2- patients, showing clustering structure and overall survivability distributions per cluster. h Differential expression analysis identifying distinct protein markers between clusters, with Cluster 0 corresponding to high-survival and Cluster 2 to low-survival patients. Source data are provided as a Source Data file.
We then performed post hoc clustering after model training by applying the Leiden algorithm to the learned sample-level embeddings31. Based on the resulting clusters, patients were stratified into three subtype-agnostic groups exhibiting distinct survival distributions (Fig. 3b). Patient groups obtained by clustering of sample embeddings showed significant global survival separation in JacksonFischer dataset (global log-rank p = 0.00253; Supplementary Fig. 5). Clinical subtype composition was summarised only post hoc to contextualise these embedding-derived groups: Cluster 0 was enriched for high-survival patients with high variance, showing the highest proportion of HR+HER2+ samples (45.7%) and the lowest proportion of TNBC samples (22.7%) (Fig. 3b, Supplementary Fig. 1). In contrast, Cluster 2 was enriched for low-survival patients, showing the highest proportion of TNBC samples (41%) and the lowest proportion of HR+HER2+ samples (12.4%). Next, we isolated the important regions and examined differentially abundant proteins between Cluster 2, enriched for low-survival samples, and Cluster 0, enriched for high-survival samples. These regions are not restricted to a single biological compartment, but may include tumour, stromal, immune, and vascular cells. Next, we identified spatially contiguous protein expression patterns that distinguished high- and low-survival clusters in the overall dataset (Fig. 3c). In particular, proteins associated with high-survival patients (Cluster 0) differed from those in low-survival clusters (Clusters 2). The identification of localised expression of Slug, Ki-67, CD3, and c-Myc in lower survival patients from the overall dataset suggests their association with breast cancer progression. For example, Slug is a transcription factor involved in epithelial-to-mesenchymal transition (EMT), which enhances tumour invasiveness and metastasis32. Ki-67, a widely used proliferation marker, indicates aggressive tumour growth, and elevated Ki-67 levels are generally associated with worse survival in breast cancer33,34. Finally, c-Myc, an oncogene involved in cell cycle regulation and proliferation, is frequently overexpressed in aggressive breast tumours and is associated with poor prognosis35. The combination of these markers indicates a tumour biology characterised by high proliferation (Ki-67, c-Myc), EMT-driven invasion (Slug). To illustrate these marker level findings at the tissue level, we visualised one representative lower-survival sample from Cluster 2. SpaCEy highlighted spatially localised important regions within the tissue section, and zoomed views of these regions showed locally elevated Slug and c-Myc abundance within the model-selected subgraphs (Fig. 3c, bottom) which indicates that the associated marker signals are not only detected at the aggregate level, but are also spatially organised within important regions of individual tissue samples. Notably, detecting CD3 in lower survival patients in the overall dataset suggests heterogeneity in immune cell infiltration, with interpretation depending on whether CD3-positive cells are located in tumour or stromal important regions. To determine whether survival-associated marker signals were localised to specific tissue compartments rather than reflecting global tissue-wide changes, important nodes were first separated into broad cellular compartments, including tumour, stroma, immune, and vessel. Marker abundance was then summarised as the sample-level mean within each compartment (Supplementary Table 6). For example, a higher CD3 signal in stromal-context important nodes means that the ROI-level mean CD3 abundance among important nodes assigned to the stromal context is higher in lower-survival ROIs. This suggests the presence of stromal regions with a high abundance of CD3-enriched cells in patients with poorer survival, whereas such regions are largely absent in the better-survival group. Given the immune-excluded phenotype, therapeutic strategies aimed at stromal remodelling or reprogramming, such as combination therapy with anti-TGF-β agents, or at normalisation of the tumour vasculature may promote T-cell infiltration into the tumour and potentially improve survival in patients with poor prognosis36. In the JacksonFischer cohort, Ki-67 had its strongest lower survival shift in tumour-context important nodes, whereas CD3, Slug, and c-Myc had their strongest lower survival shifts in stromal-context important nodes. For each marker, the reported compartment corresponds to the context with the largest positive standardised mean difference between lower- and higher-survival samples. We therefore interpret these markers as compartment-resolved signals within important regions rather than as global tissue-wide marker changes or direct evidence of cell-intrinsic mechanisms.
To further dissect cellular composition within these key clusters, we performed cell type proportion analysis and identified the highest proportional differences based on the nodes identified as important in the JacksonFischer dataset (Fig. 3d. Specifically, Cluster 0 exhibited important nodes with a higher abundance of hormone receptor-positive cells known to respond to endocrine therapy and, which may contribute to favourable survival37. In contrast, Cluster 2 was enriched with p53+ EGFR-positive, hypoxic, basal, and highly proliferative cancer cell phenotypes. Breast tumours with alteration in TP53 have a worse prognosis38, and EGFR overactivation drives proliferation and survival signalling39, and has been associated with poor clinical outcomes40. Both hypoxia and high proliferation indicate rapid tumour growth and aggressive clinical behaviour41,42. Basal cells, also known as myoepithelial cells, express markers of EMT and contribute to invasive behaviour43. Collectively, areas with a higher proportion of cells exhibiting upregulation of oncogenic proteins, signatures and processes linked to rapid tumour growth, metastasis and therapy resistance is expected in the poor survival group. Similar associations between phenotypic composition and patient outcomes have also been reported in previous studies28,44. We further characterised the localisation of important regions within the tissue (Fig. 3e). We found that 55.2% of these regions were in tumour areas suggesting that the abundance values in the tumour areas are the most informative for the survival prediction. The second most important regions are the peritumoural areas suggesting that tumour microenvironment plays a critical role in disease progression and patient outcomes. For example, hypoxia in the tumour microenvironment can lead to genetic instability and more aggressive tumour phenotypes. In breast cancer specifically, hypoxic regions are linked to metastasis, treatment resistance, and poor prognosis and these factors may contribute to the observed spatial distribution of important regions and their association with patient survival45,46. An example case is given in Fig. 3f where these important regions predominantly align with the tumour core and peritumoural areas. Finally, we assessed the composition of important nodes across the JacksonFischer dataset and found that they were distributed across multiple tissue compartments, with mean proportions of 55.2% tumour, 29.3% stroma, 13.5% immune, and 2.0% vascular compartments. This indicates that identified regions by the explainer are not confined to a single compartment, but instead capture outcome-associated local tissue structures spanning mainly tumour and stromal contexts, with additional immune and vascular contributions. Given this mixed composition, we focused on the stromal compartment by quantifying the cellular composition of important nodes at the sample level. To directly assess whether stromal-associated marker signals were driven by a higher stromal fraction among important nodes, we compared proportions of important nodes per compartment between lower- and higher-survival groups. In JacksonFischer, lower-survival samples showed a lower stromal fraction of important nodes than higher-survival samples (27.7% vs. 30.9%; Mann–Whitney p = 0.018), indicating that stromal marker enrichment should not be interpreted as evidence that a larger stromal area alone explains the lower-survival phenotype. Rather, these markers represent marker-level signals within heterogeneous important regions, which we further interpret using context-stratified marker summaries where annotations are available.
SpaCEy explains the heterogeneity of the survival of patients within a clinical subtype
Hormone receptor-positive, HER2-negative (HR+HER2-) breast cancer represents the most prevalent molecular subtype, accounting for the majority of breast cancer diagnoses. Although patients with HR+HER2- tumours generally exhibit more favourable prognoses compared to more aggressive subtypes such as TNBC, substantial heterogeneity remains in disease progression, therapeutic response, and overall clinical outcomes47. This variability presents a significant challenge for precision oncology, as current classification systems may not adequately capture the biologically relevant distinctions within this subgroup. To address this gap, we investigated the latent molecular architecture of HR+HER2- tumours using spatially resolved proteomic data. As shown in the UMAP embedding, patient samples with similar survival tend to cluster together, specifically, reduced survival in cluster 2 and higher survival values in cluster 0 in Fig. 3g. Notably, samples within HR+HER2- subtype show varying survival, underscoring the biological heterogeneity that exists within this molecular subtype and can be group into three clusters. Kaplan-Meier curves for the clusters in HR+HER2- clinical subtype (global log-rank p = 0.03132) are shown in Supplementary Fig. 7.
Beyond patient re-stratification, SpaCEy allows for improved insights into previously well-established clinical subtypes and heterogeneity of outcome. When comparing explainer features between the entire dataset and the HR+HER2-subgroup (Fig. 3h), we observed a significant overlap in the survival-associated proteins. To test whether the overlap between HR+HER2- and the full cohort was driven by the prevalence of HR+HER2- samples, we performed a targeted subtype analysis restricted to the identified markers associated with lower survival patients: Slug, c-Myc, CD3, and Fibronectin. For each clinical subtype and for the combined non-HR+HER2- group, patients were split by the group-specific median overall survival. Marker abundance was then summarised as the patient-level mean abundance among hard-mask SpaCEy-important nodes (Supplementary Table 7). All four reported markers showed positive lower-survival-associated shifts in both the full cohort and the HR+HER2- subgroup. Importantly, the same lower-survival-associated trend was also observed in the combined non-HR+HER2- group (Supplementary Data 3). In this group, Slug, c-Myc, and CD3 were FDR-supported, while Fibronectin showed the same direction of change but was not statistically significant. These results indicate that the marker overlap is not solely explained by HR+HER2- dominance. Rather, SpaCEy captures a shared lower-survival-associated component across receptor-defined subtypes while retaining subtype-specific variation. The high expression of Ki-67, Slug, and c-Myc in the lower-survival samples of both the whole data set and HR+HER2- subset highlights the central role of these proteins in tumour proliferation and progression in breast cancer in general. The elevated fibronectin signal in important regions from lower-survival HR+HER2- patients is consistent with prior METABRIC-based evidence linking high fibronectin expression in primary breast tumours to decreased patient survival, while also pointing to context-dependent variation in epithelial-mesenchymal state and metastatic potential48. More broadly, we observed that high expression of apoptotic markers, such as cleaved caspase and PARP, was enriched in the high survival group. The spatial expression patterns of these proteins enhance our comprehension of the biological diversity within HR+HER2- breast cancer, and can be potentially used as a guide to personalised treatment strategies.
SpaCEy generates generalisable insights across cohorts
We used the METABRIC dataset as a second independent breast cancer cohort to train a separate model49. This allowed us to further investigate the relationship between patient survival and spatially resolved proteomic features. The METABRIC dataset comprises 460 imaging mass cytometry samples from 405 patients, capturing diverse clinical features. The preprocessing steps are described in the Data and Preprocessing subsection, and summary statistics of the METABRIC dataset are presented in Supplementary Table 3. For the METABRIC imaging cohort, Kaplan–Meier curves stratified by clinical subtype (HR+HER2+, HR+HER2-, HR-HER2+, and TripleNeg) also showed significant differences among clinical subtypes in overall survival (log-rank p = 0.0175). Similar to our previous observations, UMAP embeddings of the METABRIC dataset (Fig. 4b) demonstrated relevant patterns associated with overall survivability, with Cluster 0 being predominantly associated with high survival and Cluster 2 with lower survival. Patient groups obtained by clustering of the sample embeddings showed significant global survival separation in METABRIC dataset (global log-rank p = 0.00371; Supplementary Fig. 5). Differential protein abundance analysis identified distinct molecular markers distinguishing these survival-associated clusters (Fig. 4c). In Cluster 0 (higher survival), we observed enrichment of HER2 and CK8/18. Conversely, Cluster 2 (lower survival) was characterised by increased expression of fibronectin, vimentin, beta-catenin, and SMA. Because METABRIC contains a higher fraction of HR+HER2- patients than JacksonFischer, we performed a descriptive subtype-composition analysis. HR+HER2- patients accounted for 71% of METABRIC patients versus 63% of JacksonFischer patients. Within METABRIC, fibronectin abundance in important nodes was not enriched in HR+HER2- patients relative to other subtypes (mean 5.62 vs. 6.34; Mann–Whitney U test, p = 0.84). We therefore interpret fibronectin as an important feature associated with the lower-survival group in METABRIC dataset, rather than as a standalone prognostic biomarker uniformly conserved across cohorts. This observation is consistent with the clinical characteristics of patients within this cluster50. Notably, PR and PanCK were highly expressed across the two clusters, with PR showing relatively higher abundance in the higher survival group, suggesting a shared molecular signature among patient subgroups. The occurrence of these patterns across datasets indicates that the SpaCEy captures underlying biological signals.
Fig. 4: Survival-associated spatial patterns in the METABRIC imaging cohort.
a Kaplan–Meier estimates of overall survival across clinical subtypes in the METABRIC imaging cohort, with censoring marks shown on each curve and a global log-rank test. b UMAP visualisation of patient samples coloured by survival values and clustering results along with the distribution of survival values across different clusters. c Marker expression profiles across lower- and higher survival clusters. d Cell type proportions of important nodes in the METABRIC dataset and cell type proportion differences between Cluster 2 and Cluster 0. e UMAP visualisation of HR+HER2- samples coloured by survival values and clusters and distribution of survival values across different clusters in HR+HER- samples f Differential expression analysis identifies distinct protein markers between clusters, with Cluster 0 representing high-survival patients and Cluster 2 representing low-survival patients in the dataset. g Representative METABRIC sample showing important regions and the corresponding spatial distribution of Fibronectin. The upper panel shows node importance across the sample, with explainer-highlighted nodes and their spatial graph connections overlaid. The boxed region is enlarged in the lower panel, showing Fibronectin abundance. Source data are provided as a Source Data file.
To further dissect the cellular composition underlying these clusters, we compared the proportions of cell types across clusters and calculated significant cell type proportional differences between clusters (Fig. 4d). As in JacksonFischer dataset, we observed areas with a higher proportion of hormone receptor-positive cells in Cluster 0, while hypoxic and myoepithelial/basal were more abundant in Cluster 2, which is associated with worse survival outcomes. The METABRIC dataset also allowed analysis of tumour microenvironment and revealed areas with higher proportion of myofibroblasts and fibroblasts in patients with lower survival, demonstrating the role of fibrotic stroma and cancer-associated fibroblasts in promoting malignant growth and tumour invasion51. Across METABRIC samples, important nodes were also compositionally mixed, with mean important-node proportions of 58.0% tumour, 30.4% stroma, 9.0% immune, 1.6% vessel, and 1.0% unclassified regions. Similarly, in METABRIC, stromal important-node fractions were comparable between lower- and higher-survival ROIs (30.0% vs 30.8%; p = 0.263). Focusing on the HR+HER2- clinical subtype, the visualisation of the learned latent representations in UMAP space (Fig. 4e) reveals clustering of samples according to patients’ overall survival profiles. Kaplan-Meier curves for the clusters in HR+HER2- clinical subtype (global log-rank p = 0.01498) are shown in Supplementary Fig. 7. The distribution of overall survival across these clusters is consistent with our previous observations, with Cluster 2 enriched for patients with shorter survival and Cluster 0 comprising those with longer survival times. Next, we explored differential protein abundance across HR+HER2- patient clusters (Fig. 4f). In Cluster 0 (higher survival), we observed higher localised expression of PR while Cluster 2 (lower survival) showed an enrichment of fibronectin, β-catenin and CK19. Higher abundance of fibronectin in the lower survival group is consistent with the overall analysis above. Increased β-catenin expression in breast cancer is significantly associated with poor prognosis, as demonstrated by a study where patients exhibiting nuclear and/or cytoplasmic β-catenin expression had reduced disease-specific survival rates52. The association between PR (progesterone receptor) expression and favourable prognosis in hormone receptor-positive breast cancer is widely acknowledged. A meta-analysis demonstrated that patients with high PR expression had significantly better survival compared to those with low PR expression53. These results highlight partially overlapping spatial and molecular patterns in the METABRIC and JacksonFischer datasets, reinforcing the association between protein abundance, tumour cellular composition and patient survival outcomes, and demonstrating SpaCEy’s ability to consistently identify prognostic patterns.
SpaCEy improves prediction performance using spatially resolved data
We also benchmarked SpaCEy against other machine-learning methods that do not incorporate spatial information, using the JacksonFischer and METABRIC datasets. This comparison enabled a direct assessment of whether spatial context improves predictive performance. To this end, we compared our results with four related methods: Fast Survival SVM9, Random Survival Forest10, Gradient Boosting Survival Analysis11, and a Cox proportional hazards (CoxPH) model54. To perform this comparison, we generated pseudobulk profiles of the samples based on protein abundance values and cell type annotations. For protein abundance, we applied four aggregation functions (minimum, maximum, mean, and sum) resulting in four distinct training datasets. For cell type-based aggregation, we created two additional training datasets using mean and sum aggregators (Fig. 5a).
Fig. 5: Comparison with non-spatial baseline methods.
a Baseline comparison pipeline. b Performance comparison on JacksonFischer Dataset. Comparison of test-set concordance indices for SpaCEy and non-spatial survival baselines using pseudobulk marker summaries and cell-type-composition features. Each box summarises n = 10 fixed patient-separated train/test splits. The centre line indicates the median; box limits indicate the 25th and 75th percentiles; whiskers extend to the smallest and largest observations within 1.5 times the interquartile range; and values outside the whiskers are shown as individual points. c Similar performance comparison on METABRIC Dataset. Source data are provided as a Source Data file.
We performed a comprehensive hyperparameter optimisation for each method-aggregation combination (Supplementary Data 1), optimising a total of 24 models. To ensure methodological rigor and a fair comparison across models, we employed consistent cross-validation splits in which folds were stratified at the patient level to prevent data leakage. The results for the JacksonFischer and METABRIC datasets are presented in Fig. 5b, c, respectively. Model performance was evaluated using the concordance index (C-index) as the primary metric. The choice of aggregation functions was motivated by their ability to capture different aspects of the protein abundance values, while mean aggregation smooths out individual variations, maximum and minimum values preserve extreme signals that might be biologically relevant. Sum aggregation, on the other hand, retains total protein abundance levels, which could be important for capturing global trends in tumour heterogeneity.
JacksonFischer dataset, SpaCEy achieved a mean performance of 0.720 ± 0.061, outperforming all 24 non-spatial baseline configurations evaluated through the same fixed 10-fold protocol. The best-performing baseline method was Cell Type Composition with Gradient Boosting SurvivalAnalysis using sum aggregation, which achieved 0.617 ± 0.092. For the Cox proportional hazards (CoxPH) baseline, the best-performing variant used pseudobulk marker summaries with sum aggregation and achieved a mean test C-index of 0.600 ± 0.095. SpaCEy demonstrated a substantial improvement of 0.102, representing a 16.5% relative improvement over this best baseline method.
On the METABRIC dataset, SpaCEy achieved a mean performance of 0.688 ± 0.038, again significantly outperforming all baseline methods. The best-performing baseline method was Cell Type Composition with Gradient Boosting Survival Analysis using sum aggregation, which achieved 0.570 ± 0.071. For the CoxPH baseline, the best-performing variant again used pseudobulk marker summaries with sum aggregation and achieved a mean test C-index of 0.562 ± 0.048. SpaCEy demonstrated an even larger improvement of 0.118, representing a 20.7% relative improvement over this best baseline method. The lower standard deviation of SpaCEy’s performance compared to baseline methods indicates more stable and reliable predictive capability across both datasets.
Importantly, our proposed model demonstrated the highest and most robust performance, with a clear improvement in median C-index compared to all alternative methods in both datasets. This gain was particularly evident in the METABRIC cohort, where traditional models showed more variable performance. These results indicate that our model generalises well across cohorts and outperforms existing approaches by leveraging aggregated cellular features more effectively.
To compare SpaCEy with a supervised non-spatial representation in terms of embedding-level stratification, we trained a non-spatial multilayer perceptron (MLP) baseline using pseudobulk marker features from the same samples and the same censoring-aware Cox loss used for SpaCEy. We extracted sample-level MLP embeddings, clustered them into three groups, and evaluated patient-level survival using Kaplan–Meier curves and a global log-rank test across all clusters. The non-spatial MLP sample-embedding baseline did not show significant global separation in JacksonFischer (global log-rank p = 0.23210) or METABRIC (global log-rank p = 0.12239; Supplementary Fig. 6). We therefore use these cluster-level survival analyses as descriptive summaries of supervised representations, not as independent validation of predictive performance.
SpaCEy learns predictive spatial representations from CODEX data and retains signal across batches
To evaluate SpaCEy representations under batch variation, we performed a batch-focused analysis using a colorectal cancer cohort, measured with CODEX, comprising 140 tissue regions from 35 patients55. This dataset provides a same-panel setting in which cross-batch transfer can be assessed directly, avoiding the marker-panel incompatibility that limits comparisons across some spatial proteomics datasets. We used this dataset for two related purposes: first, to test whether SpaCEy can learn predictive spatial representations from CODEX data, and second, to evaluate its sensitivity to batch effects under a controlled same-panel transfer setting. We selected the CODEX colorectal cohort because matched batch samples share the same marker panel, enabling a controlled same-panel cross-batch test. A direct JacksonFischer-to-METABRIC breast cancer transfer was not used for this analysis because marker-panel mismatch (Supplementary Data 4).
The prediction task was formulated as a supervised binary classification problem using immune phenotype annotations, distinguishing Crohn’s-like reaction (CLR) from diffuse inflammatory infiltration (DII). We used a patient-stratified 80/10/10 train/validation/test split, ensuring that no patient appeared in more than one split. After hyperparameter optimisation, we repeated the selected model configuration across five random seeds while keeping the patient-level split fixed. Across these five runs, SpaCEy achieved a held-out test AUPRC of 0.846 ±0.026. These results indicate that SpaCEy can be trained on an independent CODEX multiplex imaging cohort while preserving patient-level separation between training, validation, and test sets.
We next assessed batch sensitivity using the matched tissue microarray batches A and B from the same CODEX cohort. Because each patient was represented in both batches, a naive batch holdout would leak patient information across splits. We therefore separated train/validation and test partitions at the patient level. Specifically, the model was trained and validated on batch A and evaluated on held-out samples from different patients in batch B. This A-to-B transfer setting is more stringent than the randomly patient-stratified analysis above, because the model is exposed during training only to the source batch and must generalise to samples acquired in a different batch.
SpaCEy retained predictive performance on held-out batch B in this cross-batch setting, achieving an AUPRC of 0.801 ± 0.016 across random seeds. As expected, performance on held-out batch B was slightly lower than in the randomly patient-stratified analysis, consistent with the additional challenge imposed by batch transfer. To provide direct visual evidence for this A-to-B transfer setting, we visualised first fully connected layer sample-level embeddings from the cross-batch model using UMAP, with the same samples coloured by technical batch and by CLR/DII immune phenotype (Supplementary Fig. 4). The embeddings did not show strong separation by technical batch alone, with a technical-batch silhouette score near zero for A-to-B (0.0009), supporting the conclusion that the representation was not dominated solely by batch identity.
Scalability analysis results
To ensure our model remains computationally practical for large-scale spatial datasets, we assessed the computational scalability of the model across 30 configurations varying in both marker dimensionality and sample size (Fig. 6). The benchmarking workflow comprised four main stages: synthetic spatial data generation, model training and simulation, computation of performance metrics, and downstream analysis (Fig. 6a). This systematic analysis enabled quantitative evaluation of runtime and GPU memory usage under diverse experimental conditions. The complete results of scalibility analysis are given in Supplementary Data 2.
Fig. 6: Computational performance and scalability analysis.
a Overview of the computational pipeline comprising synthetic spatial data generation, model training and simulation, computation of summary metrics, and downstream analysis. b Average epoch time (in seconds) as a function of the number of samples and markers, showing the scaling behaviour of training time with increasing dataset size. c Heatmap of average GPU memory utilisation (%) across varying numbers of samples and markers, highlighting resource demands for large-scale simulations. d Scaling of average epoch time as a function of sample number for different graph sizes, illustrating the expected increase in runtime with both the number of samples and the number of cells per sample. Source data are provided as a Source Data file.
Model optimisation was performed over multiple epochs, where one epoch corresponds to a complete pass of the model through the entire training dataset. Each epoch therefore represents a single iteration of parameter updates across all samples. Average epoch time ranged from 0.63 to 22.86 s (mean = 6.89 ± 6.12 s), highlighting favourable performance across all tested conditions (Fig. 6b). For small datasets (50 samples), all marker configurations exhibited comparable performance (0.95 ± 0.43 s per epoch), enabling rapid model prototyping. Epoch time increased approximately linearly with dataset size, reaching 7.59 ± 1.63 s at 500 samples, 11.60 ± 3.56 s at 1000 samples, and 16.50 ± 3.77 s at 1500 samples. Interestingly, a non-monotonic trend was observed between marker count and training time, with intermediate dimensionality (100 markers) resulting in slower performance at large sample sizes (20.32 s vs. 12.48 s for other configurations with ≥1000 samples), potentially due to suboptimal GPU memory access patterns. Overall, the model exhibited sub-linear scaling behaviour, with per-sample processing time decreasing at larger batch sizes, suggesting improved computational efficiency for large-scale, population-level analyses.
GPU memory utilisation remained modest across all configurations, ranging from 4.36 to 41.37% (mean = 14.14 ± 9.90%), with a maximum absolute memory usage of only 235.18 MB for the largest configuration (1000 markers, 500 samples) (Fig. 6c). Memory consumption increased primarily with marker dimensionality (30 markers: 8.84%; 1000 markers: 25.67%) rather than sample size, reflecting the feature-centric architecture of the model. Notably, memory efficiency improved by 9.7-fold at larger batch sizes (0.072 MB/sample at 1500 samples vs. 0.701 MB/sample at 50 samples), confirming that the model scales favourably in both compute time and memory usage.
To address whole-slide-scale settings directly, we updated the scalability benchmark to vary both the number of cells per sample and the number of samples. One synthetic sample was represented as a spatial point cloud (cells) with 1500, 10,000, 50,000, 100,000, 250,000, 500,000, or 1,000,000 cells and 50 marker channels (Fig. 6d, e). Samples exceeding a specified cell-count threshold were partitioned into cell tiles and processed; the learned features from these tiles were then aggregated to generate the final sample-level representation. The calculated per-sample processing time increased from 0.014 s for 1500-cell samples to 7.42 s for 1,000,000-cell samples. Scaling these measured per-sample times to epochs containing 50, 100, 200, 500, 1000, or 1500 samples gave estimated epoch times from 0.69 s for 50 samples with 1500 cells each to approximately 3 h for 1500 samples with 1,000,000 cells each. The largest setting, therefore, corresponds to 1.5 billion cells per epoch. Peak allocated GPU memory remained bounded by the tile size, reaching approximately 1.06 GB, while peak resident CPU memory was approximately 1.45 GB.
Ablation and sensitivity analyses
To assess the contribution of individual SpaCEy components and evaluate the sensitivity and robustness of the model, we performed a series of ablation studies focusing on key architectural and input-design choices. To assess the sensitivity of SpaCEy to partial field of view, we performed an inference-only crop-sensitivity analysis on the JacksonFischer dataset using the trained survival model. Each full region of interests (ROIs) was compared against five fixed 50% crops defined from the ROI bounding box (centre, upper-left, upper-right, lower-left, and lower-right). Across full ROIs and 1578 valid crop comparisons, crop-based predictions remained highly concordant with full-ROI predictions (mean Spearman ρ = 0.943), while risk-tertile agreement remained 0.823 on average. Here, risk-tertile agreement means that the crop and its parent full ROI were assigned to the same low/mid/high risk tertile using cutpoints defined from the full-ROI risk-score distribution. Crop embeddings also remained highly similar to the corresponding full-ROI embeddings (mean cosine similarity 0.920). We further assessed whether crop-derived embeddings preserved the latent-space cluster assignment of the corresponding full ROI. In this analysis, full-ROI embeddings were clustered globally with Leiden, and each crop embedding was assigned to the nearest full-ROI cluster centroid by cosine similarity with the aim of measuring whether the crop remains in the same full-ROI embedding cluster as its parent ROI. On a stratified explanation subset of full ROIs and valid crop explanations, the cell-type composition of important nodes was stable (mean Pearson correlation 0.83–0.97). The cell-type composition correlation was calculated by restricting the full-ROI explanation to cells present in the crop, converting the important nodes from the restricted full ROI and crop explanation into normalised cell-type proportion vectors and computing the Pearson correlation between the two vectors. These results indicate that SpaCEy is moderately robust to partial observation, but ROI restriction still remains an important limitation for both prediction and interpretation. The crop-level prediction, embedding, and explanation-stability summaries are shown in Supplementary Table 5.
In a second ablation study, we asked whether explicit removal of identified important regions alters SpaCEy predictions and survival-ranking performance more strongly than removal of matched non-important regions. In the JacksonFischer cohort, we progressively removed 1%, 5%, 10%, 20%, 50%, and 80% of important nodes, together with their incident edges, and compared the performance of the trained model. Across all ablation levels, removal of important nodes produced consistently larger decreases in model performance than removal of non-important nodes and led to a stronger degradation in concordance index, indicating that explainer-identified regions had a greater effect on both individual risk estimates and survival-ranking performance. Specifically, mean decreases in model output were 0.0094 versus 0.0071 at 1%, 0.0509 versus 0.0179 at 5%, 0.0818 versus 0.0281 at 10%, 0.1034 versus 0.0455 at 20%, 0.1553 versus 0.0913 at 50%, and 0.2175 versus 0.1289 at 80% for important-node versus matched non-important-node removal, respectively. This consistent monotonic ablation pattern indicates that explainer-identified regions are more important for SpaCEy predictions than matched non-important regions, supporting the functional relevance of the highlighted nodes.
The aim of the third ablation study was to validate the explainer model in a setting with a known ground-truth spatial signal. For this, we generated a synthetic, standalone toy tumour-stroma dataset in which positive samples contained high intratumoural CD8 abundance, whereas negative samples had matched CD8 abundance shifted to the stromal compartment. The model received only synthetic marker-channel features and spatial graphs as input; the synthetic marker channels represented tumour, CD8, stromal, immune-activation, proliferation, fibroblast/SMA-like, and noise-control signals. Simulated cell-type labels (tumour cells, CD8 T cells, fibroblasts, and macrophage-like cells) and compartment labels were used only for ground-truth evaluation and post hoc interpretation, not for model prediction. The final synthetic dataset comprised 300 graphs with varying numbers of cells per graph, and a model was trained to predict the corresponding graph-level labels. When explainer model was applied to held-out test graphs, intra-tumoral CD8 nodes received higher mean importance than stromal CD8 nodes in explained graphs, with mean group importance 0.0467 for intra-tumoral CD8 versus 0.0162 for stromal CD8 (paired one-sided Wilcoxon signed-rank p = 8.29 × 10−14), supporting that the explainer selectively highlights the intended intra-tumoral CD8 spatial motif rather than extra-tumoral CD8 cells. Representative synthetic samples further showed that high-importance regions localised to intra-tumoral CD8 cells in positive graphs, while stromal CD8 cells outside the tumour compartment received lower importance scores (Supplementary Fig. 3).
