Multi-scale spatial modeling of immune cell distributions enables survival prediction in primary central nervous system lymphoma

Summary To understand the clinical significance of the tumor microenvironment (TME), it is essential to study the interactions between malignant and non-malignant cells in clinical specimens. Here, we established a computational framework for a multiplex imaging system to comprehensively characterize spatial contexts of the TME at multiple scales, including close and long-distance spatial interactions between cell type pairs. We applied this framework to a total of 1,393 multiplex imaging data newly generated from 88 primary central nervous system lymphomas with complete follow-up data and identified significant prognostic subgroups mainly shaped by the spatial context. A supervised analysis confirmed a significant contribution of spatial context in predicting patient survival. In particular, we found an opposite prognostic value of macrophage infiltration depending on its proximity to specific cell types. Altogether, we provide a comprehensive framework to analyze spatial cellular interaction that can be broadly applied to other technologies and tumor contexts.


INTRODUCTION
The interaction between tumor cells and their surrounding non-malignant cell populations, the so-called tumor microenvironment (TME), is considered one of the major drivers in oncogenesis and tumor evolution of almost all malignancies. 1 The TME consists of a wide variety of non-malignant cells, including various classes of regulatory and effector lymphoid cells of T-and B-lineage, and various other non-hematological cell populations such as accessory cells of fibroblastic, histiocytic, and dendritic nature. To adequately determine single cells in the TME within their architectural context, various cutting-edge multidimensional technologies have been developed, such as multiplex immunofluorescence (mIF) techniques (Vectra Polaris) 2,3 and more recently higher-throughput techniques including DNA-barcoded multiplex IHC (CODEX) 4 and imaging mass cytometry (Hyperion). 5,6 Now that these technologies allow for accurate quantification of TME markers and record their spatial distribution, the challenge is to handle the vast amount of information that is obtained from each experiment. Therefore, advanced computational frameworks for spatial analysis are needed to translate complex cellular interplays into clinically relevant biomarkers. Primary central nervous system lymphoma (PCNSL) is a very rare type of aggressive B cell lymphoma that most frequently involves the brain parenchyma and more rarely the spinal cord with or without involvement of the vitreoretinal space of the eyes. 7,8 Even though intensified treatment protocols have improved survival, 5-year survival has stagnated at only 30% which is obtained at the cost of severe medical and neuro-psychological early and late side effects. 9,10 PCNSL has a unique underlying biology characterized by a defined genetic-and immunological context, in which immune evasion is a key mechanism. 11 On the one hand, immune evasion is mediated via genetic inactivation of antigen presentation molecules (Beta2microglobulin, MHC class I and II) by the tumor, with subsequent loss of expression of these proteins, escaping recognition by cytotoxic CD8 + T cells. 12,13 On the other hand, immune evasion is mediated by immune checkpoint deregulation, most frequently by translocations and chromosomal copy number alterations (CNAs) involving the PD-L1/PD-L2 locus at 9p24.1. 13 Deregulation of the PD-1 pathway, by increased expression of PD-L1/PD-L2 leads to reversible inhibition of T cell proliferation and activation. Studies

Characteristics of the PCNSL patients from the HOVON105 study
Of 134/202 registered patients in the HOVON 105/ALLG NHL24 clinical trial 20 biopsy material was available for pathology review and a diagnosis of PCNSL was confirmed. For the present study, for 88/134 patients biopsy material of sufficient amount and quality was selected, and processed for PD-L1/PD-L2 fluorescence in situ hybridization (FISH) analysis and mIF analysis. This patient cohort can be considered representative of the entire cohort of patients included in the clinical trial with minor, non-significant enrichment of male sex, worse performance status (WHO 3), and non-significant difference in survival (Table 1 and Figure S1).
We found frequent genetic alterations of the 9p24.1 region that contains PD-L1/PD-L2 in tumor cells of PCNSL, determined by FISH analysis. The types of alteration include polysomy (n = 3), copy number gain (n = 33), amplification (n = 3), and rearrangement (n = 2) of the 9p24.1 region in a total of 41 out of 88 samples (47%; Figure S2). However, genetic deregulation of PD-L1 and PD-L2 as analyzed by FISH was not significantly associated with outcome in PCNSL patients (log rank p value of 0.74; Figure S2E). iScience Article Immune-TME profiling by mIF of PCNSL cohort Formalin-fixed, paraffin-embedded (FFPE) biopsy samples of 88 PCNSL patients were used for mIF analysis using an optimized combination of primary antibody and fluorescent dyes as indicated in Table 2. From each slide, on average 15 representative regions both within tumor areas (range 4-49, total tumor regions: 1,013) and at the interface between the infiltrating tumor border and the surrounding pre-existent cerebral tissues (range 0-17, total border regions: 380; Figure 1A) were taken for multispectral imaging (MSI) analysis. For 21/88 patients, no obvious tumor border was present in the available biopsy material, and hence only tumor regions were included. Detailed analysis to determine single cells was done using InForm software ( Figure 1). In brief, raw mIF images ( Figure 1B) were segmented to identify single cells ( Figure 1C), and individual cells were phenotyped into the following 10 phenotypes ( Figures 1D and 1E): Tumor cells (PAX5 + cells) with (1) and without (2) PD-L1 expression, T cells (CD3 + cells) and their subtypes classified by expression of CD8 and PD-1 (3-6), macrophages (CD163 + cells) with (7) and without (8) PD-L1 expression, and other cell types with (9) and without (10) PD-L1 expression ( Figure 1B). Cells that did not stain for any of the markers in the panel were not included in the downstream analysis.
Extraction of multi-scale TME features from mIF data In addition to the standard non-spatial metrics such as cell counts and densities, we extracted diverse types of quantitative spatial features which collectively capture spatial characteristics of cells in PCNSL. See feature extraction in STAR Methods for a detailed list of all spatial features. Among the spatial features, the statistics that capture spatial associations between two cell types are illustrated in Figure 2A. We first used standard local and global spatial features that summarize spatial characteristics focusing on adjacent and all cells, including median distance between two cell types measured using only the most adjacent cells (i.e., median of minimum distances; local spatial feature) or all possible pairs of cells (global spatial feature).
To further dissect spatial characteristics at multiple scales, we introduced radius-based features that describe spatial patterns in the neighborhood within a specific radius around cells. These radius-based spatial features are derived from marked Poisson point processes that enable us to calculate several statistics for (1) gaps between the cells, (2) clustering behavior of individual cell types, and (3) attraction or repulsion between two cell types. An example is Ripley's L function which measures the average amount of one cell type within a specific radius from another cell type ( Figure 2B). We calculated these features based on small radii (5, 10, or 25 mm) and large radii (50, 75, or 100 mm) to cover both short-and longdistance spatial interactions between cells. Representative examples of mIF images with no spatial association, attraction, and repulsion between macrophages and T cells as captured by Ripley's L function are presented in Figures 2C and 2D.
By extracting the aforementioned local, global, and radius-based spatial features together with the standard non-spatial features, we obtained a total of n = 4,494 TME features that comprehensively characterize the TME of PCNSL patients. The high number of TME features is due to the inclusion of both coarse and fine-grind classification of total of n = 13 cell types and all possible pairs among them which added up to a total n = 58 cell type pairs ( Figure 3A). Pairing a fine-grind rare cell type (e.g., PD-1+CD8 + T cells) with an abundant coarse cell type (e.g., macrophage) makes it less prone to missing observations than pairing it with another rare cell type (e.g., PD-L1+ macrophage). To overcome persistent missing observations, we imputed missing values by the median value per feature if the feature is observed in > 50% cases (> 44 samples), otherwise the observation was excluded ( Figure S3). After filtering, we retained 2,980 TME features for the downstream analysis. We managed to obtain more features for rare cell types, partly due to more counterpart cell types available to be paired with ( Figure 3B). The spatial features are more abundant iScience Article than non-spatial features, particularly radius-based spatial features since they are obtained for all cell types (n = 13), cell type pairs (n = 58 pairs), and all six radii ( Figure 3C). Cumulatively, we comprehensively capture each cell type in the TME of PCNSL patients and the interplay among these cell types at multiple scales.

Unsupervized analysis of TME features identified subgroups of patients with distinct survival outcome
To evaluate the clinical utility of the obtained TME features, we first performed clustering analysis based on cell counts only and explored the added value of spatial TME features. Unsupervized clustering analysis of immune cell counts revealed two clusters mainly separated by the presence of tumor cells, macrophages and T cells with PD-L1 (tumor, macrophages) and PD-1 (T cells) expression, which can further be subdivided into five clusters by the abundance of macrophages and T cells ( Figure 4A). The clusters revealed a trend for survival difference (p-value of 0.064 for two clusters and 0.1 for five clusters); with a better survival in the group with higher expression of PD-1 and PD-L1 ( Figure 4B). However, a better prognostication could be achieved by abundance of either macrophages (p = 0.0041) or T cells (p = 0.008) without accounting for the expression of PD-L1 and PD-1 ( Figure S4). The clustering analysis of border regions yielded similar results with two clusters by PD-1 and PD-L1 status, but the five clusters from the border regions were not correlated with the five clusters from the tumor regions ( Figure S5). There was no apparent correlation of these clusters with the outcome of the patients ( Figure S5). Therefore, clustering analysis based on cell count alone, regardless of border and tumor regions, did not benefit from assessing PD-1/PD-L1 status.
Next, we explored the impact of spatial information in addition to cell counts and detailed immunophenotypic information using non-negative matrix factorization (NMF) clustering using the entire feature set. A  iScience Article consensus clustering analysis with a diverse number of clusters (defined by parameter K) could identify up to four stable clusters. Less stable results were obtained as of K = 5 ( Figures 4C and S6). We noted that with K = 2, both clusters identified non-spatial features as the top contributing feature. However, the majority of the clusters derived with a higher K were exclusively defined by spatial features, particularly radius-based spatial features (e.g., clusters 3-4 with K = 4; Figures 4D and S7 for K = 2 and K = 3). Significant survival differences between the clusters were observed with clustering outcomes, particularly with higher K (p = 0.069, p = 0.019, p = 0.00086 for K = 2, K = 3, and K = 4, respectively; Figure 4E), suggesting that combined spatial and to a lesser extent non-spatial information most adequately describes the biological impact of the TME. On the contrary, NMF clustering analysis of the border regions yielded a less stable clustering outcome with a less apparent correlation with the survival outcome. This may attribute to the higher heterogeneity of border regions as compared to the tumor regions ( Figure S8). All in all, spatial TME features play a prominent role in identification of prognostic subgroups in PCNSL.

Supervised analysis revealed the most predictive TME features on patient outcome
To identify which spatial TME features had the most significant impact on survival, we used supervised analysis. To this end, a Random Forest (RF) classifier was applied using the TME features from the tumor regions iScience Article to predict 12-month event free survival (EFS) as the clinically most meaningful dichotomizing survival cut point that is used in clinical practice. Thereby, the RF classifier achieved an out-of-bag AUC of 0.72 in models including all local, global and radius-based spatial TME features ( Figure 5A). The performance declined by training the same model with non-spatial features only (cell counts and densities; AUC of 0.62), further substantiating the findings on the importance of spatial features. Limiting inclusion of classes of spatial features did not improve the model (local and global spatial features only, AUC of 0.66).
For the interpretation of the RF classifier, we next analyzed the group feature importance per feature category classified by type of statistics (Figures 5B and 5C) and cell types ( Figure 5D). Among the types of statistics, radius-based spatial features with a small radius up to 25 micro-meters were identified as the most important, followed by local spatial features ( Figure 5B). Among those features, median minimum distance was the most important, followed by L statistics with radius of 10 micro-meters ( Figure 5C). In general, most of the small radius-based features were important, except for the F function that measures empty space between cells. These findings indicate that close interactions between cell types are the most informative to predict patient outcome. Overall, interactions between tumor cells and non-malignant TME cell types were associated with a good outcome, except for the interaction between PD-L1+ tumors with (PD-L1+) macrophages ( Figure 5D). Furthermore, it was apparent that macrophage content was associated with both good and poor outcomes, depending on its interaction partner. Patients had a poor outcome when the macrophages in the TME interact with T cells, and a good outcome when (PD-L1-) macrophages interact with tumor cells).
From the RF classifier, we identified the top features with the highest importance for a subset of cell type pairs. First, top features for interactions between PD-L1+ tumors and PD-L1+ macrophage were mostly enriched with radius-based features with radius of 5 and are high in poor outcome patients ( Figures 5D and  5E). In the NMF clustering results with K = 2 and K = 4, these features were also identified as the top features iScience Article in the clusters with poor survival (Cluster 2, Cluster 1, and Cluster 3 from NMF clustering with K = 2, K = 3, and K = 4, respectively; Figure S9). The top image according to the G function with radius of 5 shows a tight interaction between PD-L1+ tumors and PD-L1+ macrophages ( Figure 6A). The top features for interaction between T cells and tumors, specifically PD-L1-tumors and CD8 À T cells, reflected close interactions (radius of 5) and were high in good outcome patients ( Figures 5D, 5F, and 6B). In particular, these features may explain the survival difference of Cluster 3 and Cluster 4 from NMF clustering with K = 4, which are mostly derived from the poor surviving group from NMF clustering with K = 2 ( Figures 4E and S9). However, these features individually were not predictive of patient outcome ( Figure S10).
Of interest, we identified important roles of the spatial context of T cells and macrophages to predict patient outcome. First, spatial features reflecting close interactions between CD8 + T cells and macrophages are high in poor outcome patients, among which median minimum distance being most important iScience Article according to RF (Figures 5G and 6C). Furthermore, interaction between CD8 + and CD8 À T cells as determined by the L function with a radius of 10 was associated with poor outcome (Figures 5H and 6D). These features individually were not predictive of survival, unlike multivariate classification by RF ( Figure S10). Finally, the distribution of macrophages is an important predictor for patient survival (Figures 5H and  6E). A high interaction between macrophages determined by L statistics with radius of 25, suggesting clustering of this cell type, was the most important according to RF and was significantly associated with a poor outcome (p value of 0.00069 and FDR of 0.048; Figures 5H, 6E, and S10). In contrast, a similarly high abundance of macrophages is not associated with a poor survival, when their interaction is low and thus presents a diffuse distribution (example in Figure S11). Taken together, a multi-scale computational framework combining frequency, architectural and spatial interaction parameters is optimal to determine the prognostic impact of TME characteristics in PCNSL.  iScience Article subgroups that were specifically dictated by spatial TME characteristics rather than mere cell abundance. The high concordance between the top features from independent analysis by NMF and RF underpinned the robust prognostic value of TME characteristics. Particularly, the most important features for predicting patients' outcomes were close spatial interactions within a radius of 25 mm ( Figure 5B), but we identified varied importance between different scales. For instance, G and K functions were the most important at a radius of 5 and 25, whereas the L function was the most important at a radius of 10 ( Figure 5C). This observation confirms the complementary roles of the features at diverse scales, which overall improved the performance of the RF model ( Figure 5A).
Previous studies on the TME of PCNSL have been mostly based on non-spatial techniques, such as transcriptome analysis 16 or counts of single-or multiple markers by immunohistochemistry of single or few markers. 17,18,[21][22][23] These studies suggested that a high number of PD-L1 positive macrophages may be associated with a better overall survival, 16,21,22 while the presence of tumor-associated macrophages overall rather contribute to a worse survival. 15 Furthermore, presence of perivascular and central tumor T cell infiltrates may be predictive for a favorable outcome, although the implication of their interaction with other cell types are largely unresolved. 17,23 In the present study, we confirmed significantly worse survival in patients with a high number of macrophages and found a trend for worse prognosis in patients with a low number of T cells. However, we showed that spatial context matters more than the abundance of these cell types. For instance, the abundance of macrophages overall was deemed a poor prognostic characteristic, but this impact did not apply in a spatial context with close spatial interaction between macrophages and T cells. In fact, our RF model showed that the top features to predict clinical outcome were close spatial interaction of tumor cells and PD-L1+ macrophages for a good clinical outcome and close spatial interaction of macrophages and T cells for a poor patient outcome. Importantly, we showed that close spatial interactions between different cell types were most predictive for clinical outcome in this dataset, while spatial interactions at large distances (e.g., global spatial feature and radius-based spatial features with a large radius) between cells and non-spatial features (counts and densities) were less predictive and thus not selected in the RF model. These results underscore the importance for taking spatial characteristics into account when studying the TME, not readily feasible by standard histopathological evaluation.
Although our method was developed in a very specific type of malignant lymphoma, spatial omics has a much broader interest and is currently a major focus of attention in the field of solid oncology. For instance, in melanoma, the close proximity of specific interacting myeloid cell populations 24 and specific CD8 + T cell subset of tumor-infiltrating lymphocytes (TIL) were highly predictive of response to immune checkpoint inhibition. 25 Furthermore, studies in gastric and colorectal carcinoma have shown that spatial architecture at the tumor/stromal interface as well as specific cellular interactions at the interface is associated with metastatic behavior. [26][27][28] And indeed, such parameters may be predictive of outcome. 29,30 Further machinelearning-based integration of spatial context with other clinical information has shown to be able to predict prognosis, for instance, in colorectal cancer, breast cancer, 31 and muscle-invasive bladder cancer. 32 These examples underpin the significance and broad applicability of developing dedicated algorithms for spatial analysis of TME to support designing cutting-edge personalized therapeutic strategies.
The sequential steps of extracting a large number of TME features at multiple scales followed by collapsing information by unsupervised/supervised analysis to reveal the key TME features with the highest biological or clinical impact may provide a basis for focused functional studies. The computational framework presented can also be applied to more immune cell populations, such as CODEX and mass-tag based approaches that can profile more immune cell populations. 33 These imaging-based technologies provide an integral method to deepen our understanding on complex cellular interplays within the TME of malignancies of various nature, including malignant lymphoma.

Limitations of the study
Our study has several technical limitations, including manual steps involved in image analysis, a limited number of cell types analyzed, and the analysis of only small regions of the tissue section instead of the entire tissue.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:   d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Patients sample collection
Formalin-fixed paraffin-embedded (FFPE) tumor samples and clinical data were obtained from PCNSL patients enrolled in the HOVON105/ALLG NHL 24 intergroup, multicenter, open-label, randomized phase 3 study (NTR2437 and ACTRN12610000908033) 20 through the HOVON Pathology Facility and Biobank. The pathology on all cases was reviewed according to HOVON guidelines (see Table 1). The HOVON105/ALLG NHL 24 study was conducted in accordance with the Declaration of Helsinki and Good Clinical Practice, and approved by all relevant institutional review boards. Written informed consent was obtained from all patients, including use of biopsy material for research purposes.

Multiplex immunofluorescence (mIF)
Multiplex immunofluorescence was performed on 4-mm-thick formalin-fixed, paraffin-embedded whole tissue sections using the Opal 7-color fluorescence immunohistochemistry (IHC) kit (Akoya biosciences, USA), as previously described. 34 Specifically, slides were deparaffinized and rehydrated, followed by a blocking step for endogenous peroxidase using 0.3% H 2 O 2 /methanol and fixation with 10% neutral buffered formalin (Leica Biosystems, Germany). Slides were washed in Milli-Q water and 0.05% Tween20 in 1x Tris-Buffered Saline (TBS-T). Antigen retrieval was done by placing the slides in 0.05% ProClin300/Tris-EDTA buffer pH 9.0 in a microwave at 100% power until boiling, followed by 15 min at 30% power. Slides were cooled in Milli-Q water, washed in 1x TBS-T and blocked with Antibody Diluent (Agilent, USA). The slides were then incubated with primary antibody diluted in Normal Antibody Diluent, followed by incubation with the broad spectrum HRP from the SuperPicture Polymer Detection Kit (Life Technologies, USA). Next, the slides were incubated with Opal TSA fluorochromes diluted in an amplification buffer (Akoya biosciences, USA). The primary and secondary antibody complex was stripped by microwave treatment with 0.05% ProClin300/Tris-EDTA buffer at pH 9.0. The combination of primary antibody and fluorescent dyes is indicated in Table 2, in order of staining. Finally, DAPI working solution (Akoya biosciences, USA) was applied and the slides were mounted with Prolong Diamond Anti-fade mounting medium (#P36965; Life Technologies).

Image acquisition and quantification
Stained slides were scanned using the Vectra Polaris Automated Quantitative Pathology Imaging System (Akoya biosciences, USA). From each slide, representative tumor regions and regions at the junction of tumor and surrounding cerebral tissue were selected and multispectral imaging (MSI) images were acquired at 40x resolution ( Figure 1A). After image capture, the images were spectrally unmixed ( Figure 1B) and analyzed, using supervised machine learning algorithms within Inform 4.2.2. (Akoya biosciences). In brief, the Tissue Finder software integrated within inForm utilizes automated segmentation to identify and locate individual cells. This is followed by quantification of labels and stains ( Figure 1C) and eventually allows for phenotyping on a per-cell basis ( Figure 1D). Cells were assigned into ten different phenotype categories:  (Figures 1D and  1E). In the training phase, cells representative for each category were manually selected after which the InForm software predicts the phenotype for all remaining cells. The InForm software utilizes a multionomial logistic regression classifier with softmax output 35 for training purposes. To increase the accuracy, the phenotypes were verified manually, and adjustments were made to the training process by feeding additional cells with given phenotypes into the training process. This algorithm was trained using one image per sample and subsequently applied to all images within that same sample. Data tables for each patient were exported for subsequent analyses.

Feature extraction
For each cell type (pair) several features were calculated which represent different aspects of the nonspatial or spatial distribution of the cells. Marked Poisson point processes in two-dimensional space were used for these calculations, where points represent the cells and marks the phenotypes. Median spatial score. For each Tumor cell, calculate the spatial score, i.e. the distance to the nearest Tcell divided by the distance from that Tcell to the Macrophage that is nearest to that Tcell. Let s 1 ; .; s N denote these spatial scores, where N is the number of Tumor cells. Define the median spatial score by MSS = median s 1 ; .; s N Median absolute deviation of spatial scores. Define the median absolute deviation of spatial scores by MADSS = median absolute deviation s 1 ; .; s N

Radius based spatial features
Radius based spatial features describe spatial patterns in the data at a specific spatial scale, i.e. a specific radius around cells. In each of the following subsections we define a function which is evaluated at 6 different radii (r = 5, 10, 25, 50, 75, 100 micro-meters) yielding 6 features per function per cell type (pair). We refer to radii r = 5, 10, 25 as small radii, and to radii r = 50, 75, 100 as large radii. Four distinct functions were evaluated at each of six radii, yielding 24 features per cell type (pair). Nearest neighbor function (G function). The nearest neighbor function is closely related to the empty space function. The difference is that the nearest neighbor function looks for disk-shaped empty space centered at a cell of a certain type, while the empty space function looks for disk-shaped empty space in general. For each pair of cell types i; j (possibly i = j) and a point x in the MSI, let d j ðxÞ be the distance from x to the nearest cell of type j. For each pair of cell types i; j and radius r,

Empty space function (F function
where the mean is over all cells x of type i. If i = j then d j ðxÞ is the distance from x to the nearest cell of type i not equal to x, so that d j ðxÞ is always strictly larger than 0. If the cells of type i and the cells of type j are distributed as two independent Poisson point processes with densities l i and l j , respectively, then the theoretical expected value of G ij ðrÞ equals G theo ij ðrÞ = 1 À exp À À l j pr 2 Á For each pair of cell types i; j and radius r, define the nearest neighbor function bỹ where sðG theo ij ðrÞÞ is the theoretical standard deviation as calculated in Baddeley et al. 19 In addition, we calculateG iÃ ðrÞ for each cell type i, radius r, and cell type j replaced by all cells of type unequal to i. Ripley's K function. Ripley's K function measures the average amount of cells of type j in the neighborhood of (i.e. less than a certain radius from) cells of type i. A normalization is applied to make the value independent of the density of the cells. High values of Ripley's K function indicate attraction between cells, while low values indicate repulsion. For a pair of cell types i; j (possibly i = j) and radius r, K ij ðrÞ = l À 1 j mean i N j ðBðx; rÞÞ where l j is the density of cells of type j, the mean is over all cells x of type i, and N j ðBðx; rÞÞ is the number of cells of type j that are at distance less than r from the cell x. If i = j then the cell x itself is excluded from the calculation of N j ðBðx;rÞÞ. If the cells of type i and the cells of type j are distributed as two independent Poisson point processes with densities l i and l j , respectively, then the theoretical expected value of K ij ðrÞ equals Finally, average feature values per patient for tumor and border regions were obtained for the extracted features while ignoring missing values when taking the average. Features with more than 50% cases of persisting missing values were excluded, while the rest of the missing values were imputed by the median value per feature. Median imputation is a simple and effective method that has little impact on downstream clustering and classification analyses. The feature extraction script was written in R and is available on GitHub (https://github.com/tgac-vumc/spatstat_vectra). Part of the script uses the package spatstat. 19

Unsupervized analysis
Hierarchical clustering analysis was performed using the R packages pheatmap (v1.0.12) with Pearson correlation distance measure and Ward linkage. Consensus NMF clustering was performed using the R package NMF (v0.21.0) using Brunet update rules, with 20 iterations to determine a consensus adherence of samples per cluster and the number of clusters between 2 and 5. Top contributing features per cluster are selected by those with relative basis contribution above 0.8 compared to the maximum contributing feature (also implemented in the NMF R package).

Supervised analysis
Supervised analysis using random forest (RF) was performed in R (v3.6.1) using the package randomFor-estSRC (v2.9.3). The number of trees was set to 100,000. Standard settings were used for the number of variables randomly selected as candidates for splitting a node (square root of number of features), and for the average number of unique data points in a terminal node (1). Predictive performance was measured by the area under the receiver operating characteristic curve (AUC) based on out-of-bag predictions. The importance of individual features and groups of features was determined by calculating (joint) permutation variable importance. 36 ll OPEN ACCESS