SOTIP is a versatile method for microenvironment modeling with spatial omics data

The rapidly developing spatial omics generated datasets with diverse scales and modalities. However, most existing methods focus on modeling dynamics of single cells while ignore microenvironments (MEs). Here we present SOTIP (Spatial Omics mulTIPle-task analysis), a versatile method incorporating MEs and their interrelationships into a unified graph. Based on this graph, spatial heterogeneity quantification, spatial domain identification, differential microenvironment analysis, and other downstream tasks can be performed. We validate each module’s accuracy, robustness, scalability and interpretability on various spatial omics datasets. In two independent mouse cerebral cortex spatial transcriptomics datasets, we reveal a gradient spatial heterogeneity pattern strongly correlated with the cortical depth. In human triple-negative breast cancer spatial proteomics datasets, we identify molecular polarizations and MEs associated with different patient survivals. Overall, by modeling biologically explainable MEs, SOTIP outperforms state-of-art methods and provides some perspectives for spatial omics data exploration and interpretation.

With simulation 4 (see Methods), we attempted to use different embedding methods to compute the distance in the embedded space as the ground distance of earth mover's distance, which is further used to construct the MECN graph. These embedding methods include most popular ones, such as PCA (a), UMAP (b), diffusion map (c), ForceAtlas (d), and PHATE (e). All these embedding methods bended the simulated manifold of cells, thus leading to simply utilizing Euclidean distance problematic. The PAGA methods only capture the connectivity between closely related clusters, and simply treat far clusters disconnect (f). We also extended the diffusion pseudo-time (DPT) method by iterating DPT by assigning the root cell to each cluster center. This practice is time consuming, and fails to obtain comparable distance values obtained by each run of DPT (g). Our proposed CGMGD method best approximated the mutual manifold distance between cell clusters (h).  c). d-f, SDM performed by SOTIP shown in 3D (d) and in silico sliced sections along zaxis (e-f). g, SOTIP successfully stratifies different tissue regions (i.e. tumor and stroma region), consistent with the marker enrichment.

Supplementary Figure 11. Simulation 5.
This simulation is similar with Fig. 2e-f, in which the aim is to detect sample specific boundaries. The difference is that, in simulation 3 ( Fig. 2e-f), the three cell clusters share the same distance in gene expression manifold, while in simulation 5 (this figure), the C1-C2 are closer, while C3 is further from C1 and C2 (a).     Randomly positioned these three cluster of cells as three sequential bands shaped tissue bulk on a two dimensional plane.
In order to verify the SHN methods could not only identify the existence of SHN peaks, but also distinguish the differences in SHN values.

Supplementary Notes
Supplementary Note 1. Parameter discussion.
The parameter setting is consistent across this manuscript. There are three parameters for user to choose according to their applications.
The first one is the resolution of Leiden clustering. This parameter controls the initial clustering resolution before MECN definition. Since SOTIP could incorporate the distances between cell clusters into MECN representation, when the resolution is set too high (overclustering), the algorithm can still remedy the error with its characteristics, but when the resolution is set too low (under-clustering), the reliability of constructed MEG could be compromised. As a result, users should rather set this parameter too high than set it too low. Across our manuscript, for datasets with true number of cell clusters, we provided a function for searching the modest resolution, and for datasets without true number of cell clusters ,we set the Leiden resolution to 2. We also tested SOTIP's SHN quantification performance with different clustering resolutions on osmFISH ( Supplementary Fig. 2a,b), seqFISH+ ( Supplementary Fig. 2e,f), and Visium datasets (Fig. 2i, Supplementary Fig. 2i,j), the results showed that SOTIP performed better than other methods, even with inappropriate setting of parameters (e.g. over clustering).
The second one is how to define MECNs. We don't use the Voronoi tessellation graph because it is reported that cells resided in microenvironment could perform interplay within broader scales than nearest cells. Since different tissues have different ranges of cell-cell interactions (for example the neuron cells have wider spatial range of interplays because of the presence of axons), using the KNN spatial graph would allow users to freely adjust the scale of MECN according to applications or tissue prior. Across the manuscript, we used kNN search for each cell to define k nearest cells/spots as a microenvironment. For all datasets in our manuscript, k is set to 10 consistently, except for the SHN quantification in cerebral cortex FISH dataset for evaluating SOTIP's robustness to MECN size. SOTIP's SHN quantification performance is also validated to be better than other methods, with various choices of MECN sizes ( Supplementary Fig. 2a,b,e,f).
The third one is the number of clusters in SDM identification. Following the widely used methodology 5,6 , we set the number of clusters to the true number of clusters. In the real applications, where the true number of clusters is not known a prior, SOTIP can provides users with a multi-level tissue partitioning hierarchy, from which tissue regions with different resolutions can be investigated.

Supplementary Note 2. Further explanations on SHN.
SOTIP's SHN (Spatial heterogeneity quantification) module is used to quantify the gene expression variance of cells within a spatial neighborhood. Since each cell has a spatial neighborhood, containing a number of cells, the output of SHN is a scalar value per cell. If a cell's associated spatial neighborhood contains cells with very different gene expression profiles, the SHN value of this cell is high, meaning spatially more heterogeneous. The only control method of SHN is NUCC (IGD mentioned in the original manuscript is not a previous method, but a variant of SOTIP), which simply regards the number of unique cell clusters within each neighborhood as spatial heterogeneity. SOTIP-SHN considers the gene expression variation within neighborhood, and NUCC considers the cell cluster variation within neighborhood. That being said, SOTIP puts the cell relationships into a gene expression space so that distance between cells could be computed by distance metrics between gene expression profiles, but NUCC puts the cell relationships into a onehot cell type representation space so that the relative distances between cells on the gene expression manifold are missed. So, SOTIP-SHN has higher resolution than NUCC. Because SOTIP-SHN considers gene expression relationships between cells, it is not restricted by the accuracy of cell clustering, but the reliability of NUCC could be easily influenced by the accuracy of cell clustering. SOTIP-SHN is especially useful when continuity pronounced in gene expression space (as demonstrated in Figure 2a,b), or when it is difficult to choose a cell clustering resolution (as demonstrated in Figure 3i and Supplementary Figure 2a-d).

Supplementary Note 6. Comparison of SOTIP-SDM with other related methods.
For clustering procedure, SOTIP followed a hierarchical merging procedure, and the intermediate clustering result could be maintained and evaluated in just one-time running. While other three methods do not enjoy this benefit due to their different principles. Specifically, SpaGCN needs to initialize the cluster labels with predefined number of clusters, then applying the iterative clustering procedure to stochastically optimize parameters and cluster labels. As a probabilistic graphical model, BayesSpace also needs predefined number of clusters to formulize the model. It regarded the labels as hidden variable, then building a fully Bayesian model with Gaussian as conditional distribution and Markov random field as smooth prior. As to stLearn, it adopts either k-means or Louvain to perform clustering after generating the normalized graph adjacency matrix. Because the true number of clusters is not accessible, which is the common situation in real applications, all these three methods need to constantly re-run entire program to test different number clusters or resolution parameters.
For the extra information usage, as with BayesSpace, SOTIP does not require incorporate histological information (e.g., H&E) to obtain better performance than SpaGCN and stLearn. In the common situation, producing spatial omics data together with paired histological image needs additional efforts on experimental design and finer operations, so that many mainstream spatial omics researches would not actually provide paired or vertically adjacent histological images 13,18,21,22,[36][37][38] . What's more, the discontinuity of adjacent sections and the uncontrollable noise further complicate the utility of histological information.

Supplementary Note 7. SOTIP delineates microenvironment trajectory of human cerebral cortex.
Inspired by single cell trajectory inference, SOTIP provided a microenvironment trajectory (MET) inference module. The MECN graph built by SOTIP effectively represented the relationships among microenvironments. We applied PAGA 39 to delineate the topological structure of MECNs, with SpatialLIBD 20 human cortex 10X Visium datasets. We used this dataset to demonstrate MET because this dataset provides the cortical layer annotation and the cerebral cortex naturally displays an ordinal spatial pattern as ground truth. We adopted SOTIP to perform MET inference on all 12 samples, the results showed high consistency between inferred trajectory and the cortex layer order (Supplementary Fig. 3). Figure 7.

Supplementary Note 8. Literature validations for MECN 24 found in
Above analysis implied that MECN cluster 24 (CD8 T, macrophage, and Keratin+ tumor) may drive the difference between compartmentalized and mixed subtypes of TNBC. Tumor-infiltrating lymphocytes (TILs) have aroused wide interest because of its prognostic significance 40 . As an important part of TILs, CD8 T cells have been proved to be associated with lots of immune checkpoint molecules as well as the survival rate of patients 41,42 . However, a pooled analysis 43 containing 3,771 patients showed that increased TILs were not always a positive prognostic factor for breast cancer patients' survival. How TILs related to the prognosis of patients differed in breast cancer subtypes. For example, high-level TILs contradicted with the survival benefit in patients of luminal-HER2-negative breast cancer. On the other hand, a recent study revealed that macrophages played a role in excluding CD8 T cells from reaching the tumor cells, thus causing a poor clinical outcome 44 . Combining these studies, an assumption came up that TILs should help to fight against tumor but it may be hindered by other cells in the microenvironment. Once the hindering happened, it would indicate a condition where immunocompetent hosts failed to defend against the antigenic tumor. For example, the mixture of CD8 T cells and macrophages could be regarded as a malignant signal indicating bad survival chance. This assumption was supported by our finding that MECN cluster 24 specifically occurred more in the mixed group (reported to be more malignancy and output worse survival 22 ) than in the compartmentalized group across the patient cohort. Altogether, by fully integrating the spatial information and quantitatively measuring the microenvironment difference between TNBC samples, our analysis demonstrated a solid case to show the power of our MEbased computational framework. Our findings suggested a potential theory to explain why compartmentalized-pattern TNBC patients would response better than the mix-pattern ones and exhibited the potential of MECN as a powerful clinical indicator in a biologically meaningful way.

Supplementary Note 9. SOTIP is scalable with 3D spatial omics datasets.
Since there are currently no method for SDM identification in 3D spatial datasets, we next applied SOTIP on two recently published 3D spatial omics datasets, i.e., EASI-FISH 24 and 3D IMC 25 . EASI-FISH 24 dataset contains spatial transcriptomics measurement of 36,422 single cells in thick brain slice of lateral hypothalamus region. The results showed an agreement between the anatomical annotation and the low-dimensional UMAP embedding computed with MEG ( Supplementary Fig. 9a). The within-cluster-between-cluster ratio (WB ratio) also showed significant agreement between the MEG and anatomical annotation ( Supplementary Fig. 9b).
On 3D IMC 25 dataset, we applied SOTIP to perform spatial heterogeneity quantification and spatial domain identification on a Her2 breast carcinoma model, processed by 3D IMC.
The results are shown in Supplementary Figure 10. Specifically, supplementary Figure 10a shows the spatial heterogeneity results in 3D. For visualization, supplementary Figure  10b,c shows in silico sliced sections along z-axis. These results shows consistently higher spatial heterogeneity around the tumor-stroma boundary. Supplementary Figure 10d shows the spatial domain identification results in 3D. For visualization, supplementary Figure 10e,f shows in silico sliced sections along z-axis. These results shows SOTIP successfully stratifies different tissue regions (i.e. tumor and stroma region), consistent with the marker enrichment ( Supplementary Fig. 10g).