Single-cell spatial immune landscapes of primary and metastatic brain tumours – Nature

Clinical samples for IMC

A cohort of 185 patients underwent surgical resection for primary brain tumours or BrM between 2006–2019. A breakdown of the patients whose tumour samples were used in this study can be found in Extended Data Fig. 1a. The clinical data for all patients was obtained from surgical and pathological reports. All tumour samples obtained from patients with glioblastoma were previously untreated and classified by a certified neuropathologist (M.C.G.) following primary surgical debulking. A subset of tumours was removed during a second follow-up surgery (residual tumour) or following tumour recurrence (recurrent tumour). In accordance with World Health Organization 2021 guidelines²³, tumours formerly classified as IDH mutant glioblastoma are now considered grade IV IDH mutant astrocytoma; therefore, in this study, only grade IV IDH wild-type tumours were designated as glioblastoma. We further distinguished glioblastomas that resulted from progression from grade II/III tumours (Extended Data Fig. 1a). LTSs were defined as patients with an overall survival greater than three years (much longer than expected survival time), and STSs were defined as patients with an overall survival less than one year (shorter than expected survival time). Brain metastases samples⁵⁴ were surgically removed from patients bearing lung (n = 51 images), breast (n = 29 images), and melanoma (n = 19 images) primary tumours as well as a small number of bladder, colorectal, gastric, gastrointestinal and ovarian tumours, collectively called ‘other’ (n = 20 images) in this study. Pre-operative MRI images were used to determine the extent of peritumoural oedema (scored 0–3) by a neuroradiologist (S.L.). Leptomeningeal disease⁵⁵ was determined by contrast-enhancing lesions in the subarachnoid or ventricles as determined on MRI by a neuroradiologist (S.L.). All patients underwent standard of care (SOC) following surgery, unless otherwise specified. Cores (1–1.5 mm in diameter) were removed from formalin-fixed, paraffin embedded (FFPE) tissue blocks and assembled into tissue microarrays (TMAs). Within the glioma cohort, we included tumour-adjacent ‘normal’ tissues as well as primary brain tumour samples extracted from the tumour bulk, which were confirmed using corresponding haematoxylin and eosin (H&E) staining by a neuropathologist (M.C.G.). BrM samples were extracted from the tumour bulk and/or tumour–brain interface (termed BrM-cores and BrM- margins, respectively). A total of 242 tissue regions were sampled across the 185 patients including 139 high-grade glioma, 18 glioma-adjacent normal, 41 BrM-cores and 44 BrM-margins. Of these 242 regions, 142 were sampled in duplicate and 2 were sampled in triplicate for a total of 389 cores. Additionally, 39 patients with BrM had matched core-margin pairs. All surgical specimens and clinical information were obtained following written informed patient consent. Clinical information was de-identified and used in accordance with the institutional review boards of McGill University and Montreal Neurological Institute-Hospital (REB: NEU-10-066, 2018-4150).

Antibody optimization

Antibodies were optimized on control tissues including spleen, tonsil, lymph node, liver, kidney, normal lung, normal brain, lung cancer, glioblastoma and/or BrM. In Extended Data Fig. 3a and Supplementary Figs. 1 and 2, we show representative optimization images of both immunohistofluorescence (IHF) and IMC staining for all markers in our panel, with some exceptions: IHF was not performed for antibodies that were commercially available with conjugated metal isotopes (except CD20 and CD45; unconjugated forms were used for IHF), and IHF was not performed for Ki67 as the B56 clone is routinely used. A list of all antibodies can be found in Supplementary Table 1. For IHF staining, FFPE sections underwent deparaffinization and heat-mediated antigen retrieval using the Ventana Discovery Ultra auto-stainer platform (Roche Diagnostics) according to manufacturer instructions. FFPE slides were incubated at 70 °C in pre-formulated EZ Prep solution (Roche Diagnostics), followed by incubation at 95 °C in pre-formulated Cell Conditioning 1 solution (Roche Diagnostics) for a total run time of ~2.5 h. Slides were rinsed in 1× PBS and incubated for 1 h at room temperature in Dako Serum-free Protein Block solution (Agilent). An antibody cocktail was prepared in Dako Antibody Diluent and slides were incubated with primary antibodies overnight at 4 °C. Slides were rinsed with 1× PBS and incubated with secondary antibody cocktail prepared in Dako Antibody Diluent for 1 h at room temperature. Slides were counterstained with DAPI for 5 min at room temperature and mounted using Dako Mounting Medium. An AxioScan Z1 scanner was used to capture tissue images.

Immunostaining and IMC

FFPE TMA slides underwent deparaffinization and heat-mediated antigen retrieval using the Ventana Discovery Ultra auto-stainer platform (Roche Diagnostics) according to the manufacturer’s instructions. FFPE slides were incubated at 70 °C in pre-formulated EZ Prep solution (Roche Diagnostics), followed by incubation at 95 °C in pre-formulated Cell Conditioning 1 solution (Roche Diagnostics) for a total run time of ~2.5 h. Slides were rinsed in 1× PBS and incubated for 45 min at room temperature in Dako Serum-free Protein Block solution (Agilent). An antibody cocktail containing metal-conjugated antibodies was prepared in Dako Antibody Diluent at optimized dilutions. Slides were stained with primary antibodies at 4 °C overnight and subsequently washed with 0.2% Triton X-100 and 1× PBS. A secondary antibody cocktail containing metal-conjugated anti-biotin was prepared in Dako Antibody Diluent at the optimized dilution. Slides were incubated with anti-biotin for 1 h at room temperature and subsequently washed with 0.2% Triton X-100 and 1× PBS. Slides were counterstained with Cell-ID Intercalator-Ir (Fluidigm) diluted at 1:400 in 1× PBS for 30 min at room temperature, rinsed for 5 min with distilled water, and air-dried prior to IMC acquisition. IMC acquisition was performed using the Hyperion Imaging System (Fluidigm).

Data transformation and normalization

All IMC data presented were not transformed and analyses were based on raw measurements. Single-cell marker expressions are summarized by mean pixel values for each channel. For heat map visualization, expression data were normalized to the 95th percentile and z-scored cluster means were plotted.

Cell segmentation and lineage assignment

All lineage and functional markers underwent a staining quality check prior to cell segmentation. A subset of functional markers passed initial quality control, but did not stain consistently with IMC, and were subsequently removed from analysis (GM-CSF-R, M-CSF-R, PD-1, PD-L1 and CTLA-4; see Supplementary Fig. 1). Cell segmentation was done using a combination of classical and modern machine learning-based computer vision algorithms. This pipeline enables high-throughput segmentation and accurately resolves individual cells across diverse tissues and structures. Importantly, this algorithm fully automates the detection of cells, thus eliminating subjective bias. The DNA channel is pre-processed for nuclei segmentation to obtain foreground regions of interest using mixtures of generalized Gaussian distributions (MoGG). The channels are also tiled for segmentation so we can pass them as inputs for inference to the MaskRCNN model. A detailed description of our segmentation and image analysis pipeline is available⁵⁶. To assign cell phenotypes, we established a supervised approach based on canonical lineage markers, expected population abundance, staining quality, and maturity of cell lineage (Extended Data Fig. 1b). We first used k-means clustering⁵⁷ and a mixture of generalized Gaussian models⁵⁸ to create multi-level image stacks based on the staining intensity of each marker. Masks were curated for each lineage marker in the panel based on consideration of 6 levels using the following procedure. (1) Greyscale image channel is convolved with a median filter with a 3 × 3 window size. (2) Each pixel is clustered into 6 groups of intensity levels using the k-means algorithm. (3) For each channel we select all groups up to a particular level as foreground (1) and the rest are designated as background (0). (4) We apply morphological blob removal to obtain smoother binary masks, where binary blobs of a particular area are removed from masks to avoid noisy regions. (5) To further refine the accuracy of select markers, additional channel-specific morphological operations are applied by computing an additional binary mask obtained using the adaptive binarization method with a sensitivity of 0.4. This mask is then amalgamated with the mask obtained in step 4. (6) To enhance the image intensity of select channels, we apply a simple contrast enhancement filter by saturating the bottom and top intensity levels of pixels in particular channels. This process enables us to capture more accurate masks from channels when phenotyping cells within our cores.

The method of lineage assignment is represented in the following formula: for each cell c_i we consider the curated mask for each lineage marker M_k, where k = 1,…,n and n is the number of lineage markers. Let us assume ({p}_{{c}_{i}}^{j}) be the j^th pixel that lies in the surrounding of c_i and each pixel has the following presence vector based on the lineage markers:

where ({p}_{{M}_{i}}={0,text{or},1}), which determines whether the pixel ({p}_{{c}_{i}}^{j}) is positive for a particular marker. Next, to determine whether each pixel within a cell is positive or negative for a given marker, we determine the majority vector by summing over the presence of all vectors as:

where ({N}_{{c}_{i}}) is the number of pixels in the cell c_i. The maximum value in vector ({M}_{{c}_{i}}) determines the cell type assignment. Cell lineages are assigned in rank priority order (Extended Data Fig. 1b).

All code used to perform these analyses is available at https://github.com/walsh-quail-labs/IMC-Brain.

Cell–cell pairwise interaction analysis

To identify significant pairwise interaction and avoidance behaviours between cell types, we performed permutation tests of single-cell interactions as previously described^11,59. Cells within a 6 pixel radius (6 µm) were considered interacting. Significant interaction or avoidance behaviours were defined as having a P-value of less than 0.01.

Cellular neighbourhood identification

To identify spatial cellular neighbourhoods, we first computed neighbour windows, which we defined as being the number (N) of nearest cells to each cell (as indicated), as previously described⁶⁰. Each window is a frequency vector consisting of the types of N closest cells to a given cell. Neighbour windows were clustered. Cellular neighbourhood discovery on glioblastoma and BrM-cores combined (performed in 2021) was performed using Scikit-learn, a software machine learning library for Python. Clustering was performed using MiniBatchKMeans clustering algorithm version 0.24.2 with default batch size = 100 and random_state = 0. BrM-margins were excluded from cellular neighbourhood discovery due to their variable mix of tumour versus stromal content. Cellular neighbourhood analysis on glioblastoma cores alone (performed in 2022) used MiniBatchKMeans clustering algorithm version 1.1.2 with default batch size = 1,024 and random_state = 0. Every cell was then assigned to a particular cellular neighbourhood based on their neighbour window. Cellular neighbourhood prevalence in each core was normalized so the sum of cellular neighbourhood prevalence for that core was 100%. Values were then z-scored and cores with z-score above or equal to 0 and below 0 were compared for survival outcome.

IMC survival analysis

Glioblastoma survival analysis was conducted using a clinically controlled cohort of patients that received gross total resection of the tumour prior to treatment, as confirmed by post-surgical MRI, and were treated with SOC (Extended Data Fig. 9f and Supplementary Table 2). Overall survival was calculated from the date of surgery to date of death. For patients with glioblastoma whose date of death was not specified, overall survival was estimated using the date of their last known follow-up. For BrM survival analyses, local recurrence-free survival was assessed in previously untreated lesions with complete macroscopic gross total resections as confirmed by post-surgical MRI. For all Kaplan–Meier analyses, images were averaged when multiple cores were collected from the same patient’s tumour (that is, each patient had only one survival value represented in the analysis).

t-SNE

Using default parameters, t-SNE dimensionality reduction plots were generated in MATLAB (version 2019b). Clustering was performed using a customized high-dimensional spectral-based clustering algorithm, due to the curse dimensionality of our cells (order of million number of cells). In our customized algorithm, we first use the DBSCAN to isolate clusters that have a similar density (with fixed parameters of a maximum distance of 3 pixels minimum number of 30 points per cluster). This approach produces some small and some big clusters with densities that are similar to each other. The big cluster group is then re-clustered using a spectral clustering algorithm. To be able to achieve a spectral clustering result on our massive dimensional data, we do a subsampling of the data (with a subsampling rate of 10), which gives us the overall shape of the data. Next, we assign each cluster with its cluster labels obtained from spectral clustering. Finally, we fit a k-nearest neighbour classifier (with k = 5) to our labelled subsampled data, to identify the cluster labels of all samples. Markers used for t-SNE analysis include CD14, CD16, CD68, CD163, P2Y12, CC3, Ki67, CD40, CD206, HIF1α, MMP9, MPO, OX40L and pSTAT3. For visualization, expression data were normalized to the 95th percentile.

Immunohistofluorescence co-staining

FFPE tissue sections were deparaffinized and underwent heat-mediated antigen retrieval in citrate buffer pH 6.0 or EDTA buffer pH 9.0. Slides were blocked with Power Block for 5 min at room temperature and incubated with the primary antibody for 30 min at room temperature. Slides were rinsed with TBS-T and subsequently incubated with secondary antibody–horseradish peroxidase (HRP) for 30 min at room temperature. Slides were rinsed with TBS-T and stained with Opal fluorophore working solution for 10 min (AKOYA Biosciences; Opal 520: FP1487001KT, lot 202212718; Opal 570: FP1488001KT, lot 20212821). This was followed by heat-mediated antibody stripping to remove primary and secondary antibodies. These steps were repeated for each primary antibody for a total of two rounds of labelling: MPO, Abcam, EPR20257, ab208670, lot GR3390666-13, 1:500; IBA1, Fujifilm Wako Pure Chemicals, polyclonal, 019-19741, lot 41375175, 1:400; and Horse Anti-Rabbit IgG HRP Polymer Kit, Vector Laboratories, MP-7801, lot ZH0611, 1:1.

Antibody specificities and dilutions were optimized individually before multiplexing was performed. Tissue images were captured using the AxioScan Z1 scanner and processed using HALO software (version 3.5).

Clinical samples for IHF

Tissue microarrays containing glioblastoma primary tumour samples from n = 135 patients were consolidated from McGill University (Quebec, Canada), University of Calgary (Alberta, Canada) and the Netherlands Cancer Institute (NKI, Amsterdam, The Netherlands). All patients were previously untreated and classified as IDH wild-type glioblastoma by a certified neuropathologist following primary surgical debulking, and later treated with SOC. All patient information and tissues were obtained after written informed consent and used in accordance with the following ethics oversight.

McGill University cohort

n = 70 patients underwent surgical resection between 2006–2019; McGill University Health Centre and the Montreal Neurological Institute and Hospital institutional review boards (NEU-10-066, 2018-4150); a neuropathologist reviewed all cases and provided the TMA (M.C.G.). These samples represent a subset of our original IMC cohort based on tissue availability on the TMA (an independent section was used). Survival <1 year, n = 25 patients; survival 1–3 years, n = 18 patients; survival >3 years, n = 27 patients.

University of Calgary cohort

n = 58 patients underwent surgical resection between 2002 and 2020; Health Research Ethics Board of Alberta, Cancer Committee (HREBA.CC-16-0762), Clark Smith Tumour Biobank; a neuropathologist reviewed all cases and provided the TMA (J.A.C.). Survival <1 year, n = 31 patients; survival 1–3 years, n = 23 patients; survival >3 years, n = 4 patients.

NKI cohort

n = 7 patients; Institutional Review Board of the NKI-AvL and NKI-biobank (CFMPB541); TMA provided by D.B. Survival <1 year, n = 5 patients; survival 1–3 years, n = 1 patient; >3 years, n = 1 patient.

Publicly available RNA-sequencing data

Single-cell RNA-sequencing datasets were downloaded from the following.

1. (1)

   GSE154795 (ref. ⁵¹) (GSE154795_GBM.AllCell.Integrated.Scaled.ClusterRes.0.1.rds.gz). Using the Seurat object file GSE154795_GBM.AllCell.Integrated.Scaled.ClusterRes.0.1.rds, a new Seurat object was created (Seurat 4.1.1), with the RNA assay counts from the subset of the 14 new patients with glioblastoma and was normalized with the default parameters of the Seurat function NormalizeData.

2. (2)

   GSE162631 (ref. ⁵³) (GSE162631_raw_counts_matrix.zip.gz). A Seurat object was created using Seurat 4.1.1 from the expression matrix count files with the parameters min.cells = 0 and min.features = 200. The counts were normalized with the default parameters of the Seurat function NormalizeData.

3. (3)

   https://doi.org/10.17605/OSF.IO/4Q32E⁵². Using the Seurat object file seurat.obj_MNN_ref.RDS, a new Seurat object was created using Seurat 4.1.1 with the RNA assay counts of the source Seurat object and were normalized with the default parameters of the Seurat function NormalizeData.

MDMs from each dataset were characterized as CD68^high (normalized expression >0) and P2RY12^low (normalized expression <0.1) and were isolated for further downstream analysis. MDMs were subdivided by MPO^high (normalized expression >0.05) or MPO^low (normalized expression <0.05). For each individual patient, an average expression matrix was generated from the MPO^high and MPO^low MDMs. The FindMarkers function in Seurat was used to generate a list of differentially expressed genes between the MPO^high and MPO^low MDMs. Pathway enrichment was assessed using Ingenuity Pathway Analysis software v.01–13 (Qiagen). Differentially expressed genes (adjusted P < 0.05) were selected for each dataset and ‘Core Analysis’ was used with all default parameters.

The following data were also used. Transcriptomic data from human immune cells in blood⁶¹, accessed via Human Protein Atlas (https://www.proteinatlas.org/ENSG00000005381-MPO/immune+cell); glioblastoma RNA-sequencing data from The Cancer Genome Atlas (TCGA Firehose Legacy), accessed via the cBioPortal for cancer genomics (https://www.cbioportal.org); and bulk RNA-sequencing data from sorted immune cells from brain tumours¹⁶, accessed via https://joycelab.shinyapps.io/braintime/.

Statistics and reproducibility

All image analysis steps were performed in MATLAB (version 2019b) and Python (version 3.7.12). Statistical analyses were performed using GraphPad Prism 9 statistical software. P-values of <0.05 were considered significant and data were expressed as mean ± s.e.m. unless indicated otherwise in the figure legends. Normal distribution was examined via the Shapiro–Wilk test. Parametric data were analysed by Student’s t-test, one-way or two-way ANOVA. Non-parametric data were analysed by Mann–Whitney test; for large sample size comparisons, Student’s t-test was used⁶². Survival data were analysed by log-rank (Mantel–Cox) tests, as indicated in the relevant figure legends. Contingency tables were analysed by Fisher’s exact test. Tukey’s test was used for multiple comparisons. For all analyses related to survival, including Kaplan–Meier analysis and the LTS and STS cohort, images were averaged when multiple cores were collected from the same patient’s tumour to prevent biasing results toward individuals with more images. All other analyses unrelated to survival (for example, population dynamics) were performed using individual images to appropriately capture heterogeneity within the TME. Area analysis of IMC images was performed using ImageJ (version 1.53k). All antibody optimization was repeated at least two times by IHF and an additional two times by IMC, using a broad variety of tissues as shown in Extended Data Fig. 3 and Supplementary Fig. 1. Additional representative IMC images (including Fig. 1b, Extended Data Figs. 6a and 10b and Supplementary Fig. 3) were selected from 389 total images and depict the statistical changes and/or staining quality as described; similar results in staining quality were obtained for all samples included in analysis. All other representative immunostaining (Fig. 4b, Extended Data Fig. 10g and Supplementary Fig. 2) was performed on at least five full tissue samples with similar results.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.