Prediction of outcome from spatial Protein profiling of triple-negative breast cancers.

Introduction
The tumor microenvironment is a complex amalgamation of diverse interacting components, such as individual cells, structural components, and signaling proteins. These components impact pro- and anti-tumor immune activity, progression, metastasis, and survival risk. While recent studies have shown the importance of tumor-immune interactions in response to therapy and risk stratification, common descriptors of immune activity, such as the abundances of tumor-infiltrating lymphocytes (TILs) and tumor-associated macrophages, explain only a limited subset of the behaviors that pervade the tumor microenvironment1,2.
Multiplex antibody-based imaging techniques, such as imaging mass cytometry (IMC), now enable simultaneous measurement of dozens of epitopes at subcellular resolution, but computational approaches for interpretation of these data types are lacking. Most current approaches (e.g., refs. 3–5) are predicated on cell segmentation and classification, but these steps are challenging and have fundamental limitations in 2D. Furthermore, categorizing cells into a limited number of classes obfuscates the heterogeneity of differentiated and regulated states within each cell class. To address these challenges, we have developed SparTile, a computational approach for analyzing multiplex protein images. Inspired by deep learning approaches in digital pathology, SparTile breaks protein images into smaller tiles and uses sparse non-negative matrix factorization to generate low-dimensional tile representations. These tile representations are interpretable and can be used to robustly identify microenvironments (MEs) within and across tumors.
To demonstrate the power of SparTile, we have studied tissue samples from patients with triple-negative breast cancer (TNBC), which account for 15–20% of all breast cancers and carry the worst prognosis across all breast cancers. TNBCs are heterogeneous, and an improved understanding of the tumor-immune ecosystem is necessary to develop prognostic markers of risk and response to therapy1. Using SparTile, we identify and confirm MEs previously shown to correlate with survival, identify and externally validate MEs associated with survival, and distinguish fine spatial arrangements of cells that can be quantitatively associated with patient outcome. Our work demonstrates how segmentation-free analysis of multiplex antibody-based imaging of tumors can be used effectively for clinically relevant outcome prediction.

Methods

Data collection
Tissue Microarrays (TMAs) were obtained from Yale School of Medicine and were generated under Yale IRB number 0003011706, allowing construction and annotation of TMAs with waiver of patient consent. The IRB at The Jackson Laboratory approved the use of the TMA for the current study under IRB number 2019–088, and the study was determined as non-human subject research. The following process was used for sample selection. A search for TNBC was done from 2000–2012. All cases that were available in the archives were collected. The grade and stage information were not used as inclusion criteria. Expert review identified the block with the most tumor present and annotated a representative region on a hematoxylin and eosin (H&E)-stained slide from that block. The final core was extracted from the pathologist-annotated region. TMA spots were 600 to 800 microns in diameter. Two adjacent 5-micron cuts were obtained for each spot. IMC was performed on TMA spots. Images corresponding to cell lines, relapsing or metastatic patients, or patients who received neoadjuvant therapy were removed from downstream analysis. 5 patients had two microarray spots, while all other patients had 1 spot. The images corresponding to 88 primary TNBC patients were used for downstream analysis. These patients were monitored, received adjuvant chemotherapy, adjuvant radiotherapy, or both adjuvant chemotherapy and radiotherapy.
Additionally, two publicly available breast cancer IMC datasets were used for analysis: METABRIC3 and Basel4, available freely at6 and7, respectively. In the METABRIC study3, “all samples were obtained with written, informed patient consent, and the study protocol was approved by the NRES Cambridgeshire 2 Research Ethics Committee (REC ref. 07/H0308/161)” (ref. 3). In the Basel study4, the reporting summary states that “samples from the University Hospital Basel were approved by the Ethikkommission Nordwest- und Zentralschweiz (EKNZ-2014-397) and for samples to use from the Institute of Pathology and Molecular Pathology, University Hospital Zurich were approved by the Ethikkommission Kanton Zürich (KEK-2012-553).” (ref. 4). It is clarified in ref. 8 that “[a]ccording to the Swiss Federal Law for research, a positive vote of an ethical committee in a retrospective study is sufficient for using tissue for research purposes without further need of individual informed consents for the samples.” (ref. 8). These datasets are composed of IMC images of TMAs of several breast cancer subtypes. Preprocessed IMC images of these datasets, with no further preprocessing, were used for downstream analysis.

IMC data acquisition
TMA slides were incubated for 15 minutes at 58 °C in a dry oven, deparaffinized in fresh histoclear, and rehydrated through a series of graded alcohols. Antigen retrieval was performed in a decloaking chamber (BioSB TintoRetriever) for 15 minutes at 95 °C in citrate buffer, pH 6.0. After blocking in SuperBlock buffer (ThermoFisher Scientific), slides were incubated overnight at 4 °C with a cocktail of metal-conjugated IMC-validated primary antibodies and described in Supplementary Data 1. The following day, slides were washed twice in DPBS and counterstained with iridium intercalator (0.25 μmol/L) for 5 minutes at room temperature to visualize the DNA. After a final wash in Milli-Q water, the slides were air-dried for 20 minutes. The slides were then loaded on the Fluidigm Hyperion imaging mass cytometer. Each TMA spot was selected as a separate region of interest using Fluidigm CyTOF Software (7.0) and ablated by the Hyperion. The resulting images were exported as 16-bit “.tiff” files using the Fluidigm MCDViewer software.
For IMC panel design, where available, pre-conjugated antibodies were obtained from Standard BioTools. In Panel 1, we utilized all available isotopes except for 172Yb, assigning each to a distinct marker. To incorporate additional markers—particularly those focused on myeloid cell populations—we designed Panel 2, which includes 11 unique markers not present in Panel 1. While Panel 2 shares 20 markers with Panel 1 to ensure continuity and comparability, the 11 additional markers could not be accommodated in Panel 1 due to the unavailability of 11 distinct isotopes. Therefore, the two-panel approach was a necessary compromise to achieve broad cellular and phenotypic coverage while working within the technical constraints of the IMC platform. For custom-conjugated antibodies, isotope selection followed a set of strategic guidelines to optimize signal quality and minimize spillover: (a) Markers with high expression were assigned to isotopes known to produce minimal spillover. (b) Markers with low expression were paired with isotopes that generate a strong signal intensity. (c) Markers with low expression were also prioritized for isotopes that are less susceptible to receiving spillover. (d) When necessary, low-expression markers were assigned to isotopes that may produce higher spillover, balancing overall panel performance.

Data preprocessing
The IMC images were preprocessed and denoised. Each IMC image was log-normalized at the pixel level using the log(1+x) transform. For each channel, pixels with values less than a threshold (T1) and connected components (8-connectivity) with size less than T2 were removed. The thresholds used for each channel are provided in Supplementary Data 1. PNG images of all IMC markers across all spots are provided at9.

Visualization of IMC images
MCD Viewer was used to visualize raw IMC images from mcd files. For each IMC image and selected markers, the visualization parameters of MCD Viewer were selected to provide a reliable representation of marker distribution. These parameters are part of the MCD viewer software and were not used in the computational analysis via SparTile.

Virtual H&E generation from IMC markers
DNA1, DNA3, CD3, CD20, and Ki76 were added to construct the “nucleic” channel, and aSMA, HLA-ABC, CD31, Collagen, PD1, and PDL-1 channels were added to construct the “cytoplasmic” channels. The nucleic and cytoplasmic channels were multiplied by a gain coefficient, and were fed to the falsecolor package10 to generate virtual H&E images. A 2D grid search was used to optimize gains of nucleic and cytoplasmic channels by minimizing L1 distance between median RGB values of virtual H&E images and a template image (Supplementary Fig. 1, obtained from: www.pathologyoutlines.com/topic/liveradenovirus.html). While the input channels are called nucleic or cytoplasmic in FalseColor10, the goal is to mimic hematoxylin and eosin staining. We observed that adding protein markers found in both the cytoplasm and extracellular matrix improves the virtual H&E images. Furthermore, we observed non-zero expression of membrane proteins CD3 and CD20 on or immediately adjacent to pixels in the 2D space that stained positive for DNA. This signal arises from a membrane signal directly above or below the nucleus in the z-axis and is notable especially for small cells such as lymphocytes. Due to this effect, we found that treating CD3 and CD20 as nuclear markers generates darker lymphocyte nuclei in the virtual H&Es and makes them appear more similar to real H&E images (Supplementary Fig. 2). PNG images of all virtual H&E’s are provided at9.

Registering two IMC spots
We used the similarity transform of the Airlab package (v0.2.111) to register grayscale virtual H&E images to each other. The IMC image of the second panel was transformed to match the coordinates of the first image. Max projection was used for IMC markers common between the two panels.

Region type separation
IMC images were tiled using 20-micron overlapping tiles (75% overlap). Pixels with total IMC marker intensities of less than 1 were labeled as non-tissue, and tiles with less than 50% tissue pixels were removed. The average expression of IMC markers on each tile was used to compute the probability of each tile belonging to a tissue region type (tumor, stroma, immune). The following process, similar to ref. 12, was used to compute the probability of each tile belonging to each class. First, each region type was associated with its indicative markers. For example, the tumor was associated with panCK and Epcam-Ecadherin markers. The full list of markers indicative of each region type is provided in Supplementary Data 1. A quadratic link was used to compute the probability of each region type, i.e., for tile , we havewhere is the probability tile belongs to class , is the set of markers associated with class C, is the weight for marker , and is the average expression of marker m in tile . Weights are optimized by an iterative process. At each iteration, tiles were classified, and the gradient of the weight vector was used to update the weights. After training, given tile-level probabilities, the class probability vector of each pixel was calculated by averaging the class probabilities of all tiles containing the pixel.

Tile representation learning
Given a set of tiles of a given size, the pixel-level tissue region type probability maps were used to obtain tile-level probabilities by averaging pixel-level probabilities over each tile. Then, each tile was associated with the class with the largest probability. For tiles belonging to each region type, the average expression of IMC markers and IMC marker pair colocalizations was calculated. To compute colocalization of two IMC markers, first, each marker was dilated by a disk (K = 5). Then, dilated images of the two IMC markers were multiplied at the pixel level, and the average value of the product channel over the tile was saved. Given the average expression of IMC markers and marker pairs, sparse NMF (sNMF) was used to reduce the dimensionality of encoded tiles to D = 10,15, 20, 30, 40, and 50 dimensions.
IMC images were tiled using overlapping tiles with the following sizes: S = 20, 30, and 40 microns. The sNMF analysis described above was performed for each tile size. For each tissue region type, 18 sNMF models corresponding to the 3 tile sizes (S = 20, 30, and 40) and 6 dimensionality reduction parameters (D = 10, 15, 20, 30, 40, and 50), were computed. The distribution of the number of cells per tile for each S is provided in Supplementary Fig. 3.

Clustering tiles in each region type
First, tissue pixels, i.e., those with total IMC intensity larger than 1, were classified as a tumor, stroma, or immune. Given tile-level representations, pixel-level representations were obtained by averaging over the tiles containing the pixel. For each pixel, the sNMF model corresponding to the pixel region type was used for all tiles containing it. Finally, each pixel was associated with the sNMF component with the largest coefficient.

Importance scores of each cluster
An importance score was computed for IMC markers of each sNMF cluster:where is the importance of marker in sNMF cluster , is the weight of marker in sNMF component , is the weight of marker pair and in sNMF component , and is the total number of IMC markers used for analysis.

Interpretation of SparTile clusters as MEs
Top IMC markers and marker pairs of an sNMF cluster, importance scores of IMC markers, cross-referencing clusters across representations of different sNMF models, and confirmation of marker presence in raw IMC data of a cluster were used to interpret sNMF clusters and identify the MEs they encode. These interpretations were additionally used to select the hyperparameters of SparTile. We observed that using a large number of dimensions D (e.g., D = 50) reduces sNMF stability and results in noise determining many sNMF clusters (Supplementary Data 2). Therefore, we used S = 30 and D = 20 for tumor and stroma regions, and used S = 20 and D = 10 for immune regions.

Association of SparTile clusters with survival
The fraction of tissue pixels corresponding to each SparTile cluster was computed for each patient. Unless otherwise stated, a ratio threshold T = 0.05 was used to separate patients into high and low cluster ratios. Lower thresholds were used for less prevalent SparTile clusters. The log-rank test of the lifelines Python package 13 with default parameters (two-sided) was used to compare Kaplan-Meier (KM) plots of patients in high and low ratio groups. Additionally, patients with overall survival less/more than 3 years were labeled as short/long survival. The Mann-Whitney U test was used to compare the ratio of sNMF clusters between the two groups as an additional test. The statsmodel python package14 was used for multiple test correction using the “fdr_bh” method.

Spatial distance of tumors and myeloid cell markers
PanCK was used as the tumor marker, and the sum of CD68, CD11b, CD14, and CD33 markers was used as the myeloid marker. First, patients with low average expression of tumor or myeloid markers were removed (min 0.05). The remaining 80 patients were used in downstream analysis. The tumor and myeloid markers were blurred using a Gaussian Kernel with a bandwidth of 50 pixels. Afterwards, the blurred images were normalized to have a unit volume, representing the 2D distributions of tumor and myeloid markers. The total variation norm between these two distributions was used as the distance:where and are the distributions of tumor and myeloid markers, and denotes the pixels in an IMC image. A threshold of 0.25 was used to separate images with high/low distances.
For the METABRIC IMC images, we used panCK as the tumor marker, the sum of CD68, CD163, and CD11c as myeloid markers, used a threshold of 0.1 for selecting images staining positive for tumor and myeloid cell markers, and a threshold of 0.15 to separate images with high/low distance between tumor and myeloid markers. For the Basel IMC dataset we used sum of panCK, CK7, CK8/18, CK19, CK5, and KRT14 as the tumor marker, CD68 as myeloid marker, used a threshold of 0.2 for selecting images staining positive for tumor and myeloid cell markers, and a threshold of 0.15 to separate images with high/low distance between tumor and myeloid cell markers. The log-rank test of the lifelines python package11 with default parameters (two-sided) was used to compare KM plots of patients in high/low distance groups in all datasets.
The relative spatial distance can similarly be computed for any two arbitrary sets of proteins, but the number of potential distances to compute is computationally infeasible. The ME identification process can be used to narrow down the list of marker sets to use.

Multivariate Cox proportional hazard model analysis
The lifelines Python package (v0.27.4)13 was used to train multivariate Cox models for a set of independent variables. The p-values and confidence intervals of the trained Cox model are reported. For Multivariate Cox models using ratios of SparTile MEs, we used elastic-net regularization (penalty = 0.01, l1_ratio = 0.5) for stable parameter estimation. No regularization was used for other Cox models.

Correlation between spatial distance and the ratio of MEs
The difference between t0 and i2 was used for analysis. Patients with a low fraction of t0 or i2 were removed. Patients were labeled as distance high/low based on the spatial distance, and difference high/low based on the difference between t0 and i2. Fisher's exact test was used to test if there is a correlation between spatial distance and the difference between t0 and i2.

Segmentation-based analysis
Mesmer15 was used to perform cell segmentation. The area-integrated preprocessed IMC marker intensities, normalized by the square root of cell area (in square microns), were used to quantify each cell. We observed that normalizing by the square root of cell area improves the separation of cell types, deriving from membrane markers usually having stronger expression on the cell boundary. We considered cell classification and microenvironment inference under several hyperparameters. Hierarchical clustering was used to identify the n = 20, 25,...,60 cell types. We allowed a larger number of cell types to avoid missing phenotypes affected by segmentation errors or spillage of ME markers over segmented cells. The number of cells within a distance d = 10, 20,…,50 microns of each center cell was used as the neighborhood embedding for each center cell, and hierarchical clustering was used to identify c = 20, 25,…,50 MEs with different cellular compositions. Comparing different hyperparameters, we selected n = 50, d = 30, and c = 25 for subsequent analysis. Other hyperparameters worsened cell phenotyping, did not robustly encode cellular neighborhoods, or resulted in ME ratios that had weaker correlations with survival. For the distance analysis between tumor and myeloid cells, only patients with at least 10 tumor cells and 10 myeloid cells were considered. For each tumor cell, we computed the average distance to the K = 5 nearest myeloid cells (T → M), and for each myeloid cell we computed the average distance to the K = 5 nearest tumor cells (M → T). These distances were averaged across all tumor and myeloid cells, respectively. The mean of these two distances, denoted by T ↔ M, was also computed.

Statistics and reproducibility
All statistical analyses were performed in Python. The following packages were used: imctools (v2.1.8)16, numpy (v1.26.0)17, scipy (v1.8.1)18, falsecolor10, Airlab (v0.2.1)11, scikit learn (v1.1.2)19, matplotlib (v3.5.2)20, seaborn(v0.12.0)21, lifelines (v0.27.4)13, statsmodels (v0.12.2)14, and pandas (v1.4)22. Data from 88 TNBC patients were used for downstream analysis. All analyses were done at the patient level. For patients with more than one spot, the ratio of each ME across the total area of IMC spots was used for downstream analysis.

Results

Multiplex imaging of TNBC tissue microarrays
To investigate tumor-immune interactions in primary TNBC tumors, we generated IMC data from tissue microarray (TMA) spots of 88 tumors. We developed panels enabling imaging of DNA and 45 epitopes. Two (5 um) adjacent slices were registered for each spot to enable profiling of a larger marker set (Fig. 1A). This allowed identification of diverse cellular and microenvironmental markers and thorough characterization of the immune landscape (Fig. 1A, B). The panel for the first slice provided a general characterization of cell composition, while the second panel focused on immune markers. To provide a common feature space for registration of adjacent slices, we generated virtualized hematoxylin and eosin (H&E) images from the IMC data of each slice (methods, Fig. 1C). Registration of adjacent H&E images enabled spatial comparison across all assayed features, though tissue differences between slices constrained cell segmentation-based comparisons (Supplementary Fig. 4). For example, we frequently observed cells within only one of the two slices. We also observed differences in nucleus or cell membrane shapes across slices, not unexpected considering a typical cell size is 10 um and reflecting the fact that slices do not capture whole cells but instead reflect cell cross-sections (Supplementary Fig. 4).

SparTile robustly identifies microenvironments
Given the limitations of segmentation-based analysis, we developed an analysis approach to identify MEs directly from antibody image data. This approach, SparTile (SPARse non-negative matrix factorization for TILE analysis), was designed to robustly identify MEs within and across spatial antibody images. Such MEs can then be studied more deeply for interactions and spatial relationships, as well as for identification of MEs that are predictive of survival (Fig. 1C, D). SparTile first clusters tiles into three general region types: tumor, stroma, and immune (methods, Supplementary Fig. 4). Region type identification occurs at the tile scale (methods), though there may be finer structure within each tile, e.g. low densities of immune cells within predominantly tumor regions (Supplementary Fig. 4). The second step is to distinguish MEs with more specific differences. SparTile uses sparse non-negative matrix factorization (sNMF) to learn low-dimensional representations of tiles of size SxS microns within each region type (methods, Fig. 2A). Using these representations, we then learned D = 20 clusters in the tumor (S = 30), D = 20 clusters in the stroma (S = 30), and D = 10 clusters in the immune (S = 20) region type (methods, Figs. 1D, 2B, and 2C). Different markers characterize each cluster, indicating that SparTile effectively distinguishes different MEs (methods, Figs. 1D and 2B).
We found that SparTile clusters are interpretable and correspond to distinct MEs. The set of top IMC markers and marker pairs (Supplementary Data 2), as well as marker importance scores (Fig. 2B) were used to interpret the ME corresponding to each cluster (Fig. 2D and 3A–C). For example, tumor cluster 0 (t0) is an MX1+ tumor-myeloid cell ME (Fig. 2D). PanCK, MX1, and myeloid markers, such as CD14 and CD68, are frequent among the top markers and marker pairs of t0, and are observed in raw IMC images in regions labeled as t0 (Fig. 2D). Furthermore, PanCK and MX1 have high importance in t0 (Fig. 2A).

Individual tumor microenvironments correlate with survival
We observed that several MEs and expression patterns correlate with patient survival. t0 is an MX1+ tumor-myeloid cell ME. Remarkably, we found that the fraction of the spot classified as t0 correlates with worse survival (log-rank test p-value = 0.04), even though the TMA spots are only 0.6 to 0.8mm-diameter tissue samples. t0 fraction was also predictive of survival risk in a Cox model accounting for clinical variables (methods, Fig. 3A, HR = 5.25, CI = [1.32–9.18], p-value = 0.01). This ME was robustly identified and separated survival curves under other hyperparameters (e.g., D = 15, S = 40, T = 0.01, log-rank p-value = 0.015, D = 20 and S = 20, log-rank p-value = 0.05, Supplementary Fig. 5, methods). This finding is consistent with previous studies showing that the presence of myeloid cells and higher expression of MX1 in tumor ME is indicative of higher survival risk1,23. While MX1 in the tumor ME was associated with poor survival, the fraction of spots classified as i2, which is an MX1+ myeloid immune ME (Fig. 2B, Supplementary Data 2), was not predictive of survival (log-rank p-value = 0.66). Interestingly, the difference in t0 and i2 fraction (Supplementary Fig. 4, log-rank p-value = 0.004) was a stronger predictor of survival than t0.
SparTile also revealed survival associations for t4, a vimentin+ tumor ME. Vimentin and panCK have high importance in t4 (Fig. 2B) and are prevalent in the list of top markers and top marker pairs (Supplementary Data 2). We also frequently observed vimentin in manual inspection of t4 regions (Fig. 3B). Prior studies have shown that the epithelial to mesenchymal transition is characterized by high expression of vimentin on tumor cells, suggesting t4 can potentially identify tumors in this state24–27. Higher ratios of t4 correlate with poor survival (Fig. 3B, log-rank p-value = 1.1e−5). The ratio of t4 was predictive of survival risk in a Cox model accounting for clinical variables (methods, Fig. 3B, HR = 14.84, CI = [6.86−22.83], p-value < 0.005). This ME was robustly identified and separated survival curves under other hyperparameters (e.g., D = 20 and S = 2,0 log-rank p-value = 4e−4, D = 15 and S = 40, log-rank p-value = 2e−5, Supplementary Fig. 5).
Ki67 is frequently used as a marker of proliferation in tumor cells and is an indicator of poor outcome28. Importance scores (Fig. 2B), top IMC markers (Supplementary Data 1), raw IMC images (Fig. 3C), and comparison of SparTile representations with different hyperparameters (Supplementary Fig. 6) show that SparTile clusters t6, t17, and s13 are contributors to the Ki67+ MEs. Higher fractions of Ki67+ MEs correlate with poor survival (Fig. 3C, log-rank p-value = 0.052). This ME was robustly identified and separated survival curves under other hyperparameters (e.g. D = 15 and S = 20 for both tumor and stroma region types, T = 0.01 log-rank p-value = 0.06, D = 20 and S = 40 for both tumor and stroma region types, T = 0.01 log rank p-value = 0.02, Supplementary Fig. 5). Although Ki67+ status is an established marker of poor survival, we observed that Ki67+ ME fraction was not statistically significant in a multivariate Cox model accounting for clinical variables (methods, Fig. 3C, p-value = 0.4).
Recent studies have emphasized the role of AXL in resistance to chemotherapy through promotion of proliferation and migration of cancer cells29. t7 and s15 are the AXL+ MEs, and higher AXL + ME fraction correlates with poor survival (log-rank p-value = 0.03, Supplementary Fig. 6). The AXL + ME fraction was not statistically significant in a multivariate Cox model accounting for clinical variables (methods, Supplementary Fig. 6, p-value = 1). This ME was robustly identified and separated survival curves under other hyperparameters (e.g., D = 20 and S = 40 log-rank p-value = 0.06, D = 20 and S = 20 log rank p-value = 0.04, Supplementary Fig. 5).
While we observed that several tumor MEs correlate with survival, from a data-driven perspective, only t4 was statistically significant after correcting for multiple testing (methods). We believe this is due to the heterogeneity of TNBCs that results in a large number of distinct MEs with complex interactions.

Individual immune microenvironments correlate with survival
A higher prevalence of lymphocytes has generally been considered to be positive for survival1. Consistently, we observed that several lymphocyte-associated MEs within the stromal and immune region types correlate with longer survival. s9 is composed of CD8+ and CD45RO + T-cells in collagen+ and aSMA+ stromal regions, and it is associated with longer survival (methods, T = 0.03, log-rank p-value = 0.07, Mann-Whitney U-test p-value = 0.003, Supplementary Fig. 7). This ME was robustly identified and separated survival curves under other hyperparameters (e.g. D = 20, S = 20, and T = 0.01 log-rank p-value = 0.06, D = 15, S = 20, and T = 0.01 log rank p-value = 0.02, Supplementary Fig. 5). s9 fraction was a weak predictor of longer survival in a multivariate Cox model accounting for clinical variables (methods, Supplementary Fig. 6, HR = −21.98, CI = [−53.98, 9.94], p-value = 0.18). i4 is a T-cell immune ME composed of both CD4+ and CD8 + T-cells (Supplementary Data 2, Fig. 2B). Higher i4 fraction correlates with longer survival (Supplementary Fig. 7, T = 0.03, log-rank p-value = 0.01, Mann-Whitney U-test p-value = 2.1e-5). i4 fraction was statistically significant in a multivariate Cox model accounting for clinical variables (methods, Supplementary Fig. 7, HR = −31.62, CI = [−61.05, −2.19], p-value = 0.04). This ME was robustly identified and separated survival curves under other hyperparameters (e.g., D = 15, S = 20, T = 0.0,1 log-rank p-value = 0.003, D = 20, S = 30, T = 0.01 log rank p-value = 0.02, Supplementary Fig. 5).
In contrast to the T-cell microenvironments, we observed a smaller effect size for B-cells (i1, Fig. 2B). i1 fraction weakly correlates with longer survival (methods, Mann-Whitney U-test p-value = 0.03, log-rank p-value = 0.26, Supplementary Fig. 7). We observed that B-cell MEs usually comprise a small fraction of the sample area (Supplementary Fig. 7). This exacerbates difficulties in representative quantification of B-cell ME prevalence from the small TMA spots, which may contribute to the weaker associations. The B-cell ME was robustly identified across hyperparameters, i.e., tile sizes 20-40 microns and NMF dimensions D < 40, but survival associations were not statistically significant.

Spatial distribution of tumors and myeloid cell markers is prognostic of survival
Because we observed distinct associations of myeloid protein prevalence with survival in tumor (t0 ME) and immune (i2 ME) regions, we hypothesized that relative spatial distributions between myeloid and tumor cells (Fig. 4A) might contain information to further improve survival predictions. To test this, we computed the total variation norm of the distribution of tumor and myeloid markers in each TMA spot, i.e., quantifying how separated tumor and myeloid cells are (methods). We observed a negative correlation (methods, Pearson correlation coefficient = −0.25) between large distance and large difference between ratios of t0 and i2 (methods, Fisher's exact test p-value = 0.036). Lower distance between tumors and myeloid markers is associated with poor survival (methods, Fig. 4B, p-value = 0.0025). Interestingly, we observed a higher expression of MX1 in images with low tumor-myeloid cell distance (Supplementary Fig. 8, Mann-Whitney U-test p-value = 0.02). While previous studies have shown that intratumoral presence of myeloid cells, such as macrophages and neutrophils, correlates with poor outcome30,31, our results further quantify how myeloid infiltration affects survival. Spatial distance was statistically significant (p-value = 0.02) in a multivariate Cox model accounting for several clinical variables (Fig. 4B). While stage has negative risk in this Cox model, it incurs no risk in a multivariate model without lymph node status (Supplementary Fig. 8, HR = −0.07, CI = [−0.73, 0.48], p-value = 0.84), and has a weak positive risk in a univariate Cox model (HR = 0.33, CI = [−0.28, 0.75], p-value = 0.29). These results suggest that multicollinearity between variables may impact their hazard ratio in Cox models and emphasize the robustness of our spatial distance in predicting risk.
We used two external IMC datasets, the Basel IMC dataset7 and the METABRIC IMC dataset6, to further validate if the spatial distance between myeloid cells and tumors is indicative of survival in breast cancers (methods). Smaller distance correlates with poor survival in both datasets (Fig. 4c, d, log rank test, Basel IMC p-value = 0.004, METABRIC IMC p-value = 0.003). As both datasets are composed of multiple subtypes, we used a multivariate Cox model to account for subtype differences and other clinical variables. SparTile’s spatial distance was a weak independent prognostic marker of survival (Basel IMC p-value = 0.11, METABRIC IMC p-value = 0.17) in both datasets. These results are interesting as both Cox models account for a wide range of clinical variables.

Comparison with segmentation-based analysis
To further assess the robustness and interpretability of SparTile we compared it to cell segmentation-based analysis. For segmentation and classification, we analyzed only the tissue slice corresponding to panel 1. This is because individual cells are sometimes not large enough to span the adjacent slices (Supplementary Fig. 4), and panel 2 lacks markers for diverse cell phenotyping on its own (methods, Supplementary Data 3). We used Mesmer15 for cell segmentation and hierarchical clustering to identify cell phenotypes (methods, Supplementary Data 3, Supplementary Fig. 6). We then defined cell segmentation-based microenvironments based on the cell types around each individual cell. For each center cell, the prevalence of cell types within radius R = 30 was calculated. Hierarchical clustering of these prevalence vectors was then applied, yielding N = 25 MEs (methods, Supplementary Fig. 9).
For a matched comparison to the cell segmentation analysis, we re-applied SparTile to markers of only panel 1 (MEs as a function of hyperparameters are provided in Supplementary Data 2). Most MEs correlating with survival in the two-panel analysis were re-identified. For instance, the MX1+ tumor myeloid ME, Ki67+ ME, vimentin+ tumor ME, T-cell MEs, and B-cell MEs were re-identified using panel 1 alone. This is not surprising as these MEs are described by markers present in panel 1. Furthermore, the SparTile distance between tumor and myeloid markers (Fig. 4) robustly separated survival curves using the panel 1 data alone (Supplementary Fig. 10a, log-rank p-value < 0.01).
Segmentation-based analysis consistently identified several MEs also found by SparTile, such as the tumor ME with epithelial to mesenchymal transition (cluster 7), Ki67+ tumor MEs (clusters 2, 6,15, and 19) and the B-cell dominated ME (clusters 9 and 11) (Supplementary Fig. 9). The segmentation-based tumor ME with epithelial to mesenchymal transition (cluster 7, T = 0.01, log rank p-value = 0.0009) and CD8 + T-cell ME (cluster 14, T = 0.01, log rank p-value = 0.05) each showed strong correlation with survival risk, but contrary to expectations the Ki67+ tumor MEs showed weak survival associations (clusters 2, 6, 15, and 19, T = 0.01, log rank p-value = 0.25) (Supplementary Fig. 10). In contrast, with SparTile the p-value is 0.052. Moreover, the MX1+ tumor-myeloid signal found by SparTile was split among several MEs when analyzed by segmentation, regardless of hyperparameters, hindering identification of the survival risk for this ME (methods, Supplementary Fig. 9, Supplementary Data 3).
We next tested if a segmentation-based spatial distance between tumor and myeloid cells could separate survival curves (methods). We measured the average distance from each tumor cell to the 5 nearest myeloid cells (denoted by T → M, methods) and vice versa (denoted by M → T). We additionally computed the average of both distances (denoted by M ↔ T). We only observed a significant association between distance and improved survival in the M → T case (Supplementary Fig. 11). None of these metrics were significantly associated with survival in a multivariate Cox model (Supplementary Fig. 11), although M → T correlated with the SparTile distance (Pearson Correlation coefficient=0.54, p-value = 7e−8). We externally tested our segmentation-based distances on the Basel7 and METABRIC IMC6 datasets (Supplementary Fig. 11). We observed a lack of consistency across the 3 datasets. For example, for the M → T metric, in our JAX dataset and the Basel datasets, smaller M → T correlated with higher survival risk, while in METABRIC it correlated with lower survival risk, though weakly (Supplementary Fig. 11). For the T → M metric, the direction of effect was also inconsistent among the three datasets. This was the case for the T ↔ M metric as well. Moreover, even when the effects were in the expected direction (higher tumor-myeloid distance associated with better survival), p-values were often weak. In only the T → M JAX and M ↔ T METABRIC cases were the p-values < 0.01. No cases showed statistically significant effects in multivariate Cox models (Supplementary Fig. 11).

Discussion
The tumor microenvironment is a complex and dynamic network of interactions involving cells, structural components, and signaling proteins. These spatial interactions are critical for understanding tumor evolution, therapeutic responses, and survival risks. SparTile offers a segmentation-free framework for studying these interactions by leveraging the spatial distribution of membrane, structural, and signaling proteins to distinguish diverse MEs. For example, t0 is an ME characterized by tumor and myeloid membrane markers and MX1 signaling protein.
While prior studies have demonstrated that marker prevalence can predict survival independent of spatial relations32, SparTile advances this field by incorporating spatial structure at the tile level. This approach enables detailed analyses of multi-celltype-driven survival associations. For example, SparTile summarized the spatial distribution of myeloid and tumor cell markers into a single metric, which effectively separated survival curves and demonstrated prognostic value in multivariate Cox models across three breast cancer datasets.
Segmentation-based methods are valuable but have challenges, including error-prone cell segmentation, signal spillage from adjacent cells and MEs, and the computational burden of identifying and representing cellular neighborhood16,33. Recent cell phenotyping tools aim to correct for signal spillage33. This improves cell phenotyping but may filter out information important to ME inference. Graph neural networks have emerged as potential tools for understanding microenvironments from segmented cells5 but their interpretability is limited, and incorporating structural and signaling proteins into their node representations remains complex. Additionally, the influence of network architecture and hyperparameter choices on the resulting representations is not fully understood.
SparTile addresses these limitations by employing a tile-based approach that encodes spatial relationships and facilitates the direct identification of distinct microenvironments. Compared to single-pixel analyses (e.g., ref. 34), tiles provide richer contextual information, reducing the need for extensive preprocessing and hyperparameter tuning. SparTile has conceptual advantages over cell-cell contact-based metrics (e.g., ref. 12), which fail to account for dynamic cellular states, such as the persistence of T-cell activation post-interaction with cancer cells. Because SparTile integrates such dynamics, it provides a more complete characterization of spatial information.
We also note that SparTile identifies MEs in an unbiased fashion and encodes the relative spatial distribution of arbitrary sets of protein markers. As the number of combinations to test for spatial distribution analysis is computationally intractable, SparTile’s ME identification step can be used to narrow down the search. Another strength of SparTile is that MEs are more directly computed from pixel data, while segmentation-based models have multiple steps involving complex biological assumptions, including but not limited to cell segmentation, cell phenotyping, and ME identification. This likely explains the lower robustness of segmentation approaches in identifying survival associations. Nonetheless, cell-type information remains essential for many downstream applications, such as designing experiments using organoid models or patient-derived xenografts for targeted therapies. Therefore, segmentation-based methods and SparTile are not mutually exclusive, but complementary. Although we have shown that the cellular and spatial composition of SparTile’s MEs can be estimated post hoc, segmentation-based methods provide an essential complementary approach for clarifying cell-level interpretability.
In summary, SparTile’s ability to incorporate membrane, structural, and signaling proteins into a cohesive analytical framework enhances the biological interpretability of multiplex imaging data. This capability enables more precise and comprehensive analyses, advancing our understanding of tumor biology and informing prognostic assessments.

Supplementary information