Abstract

   Glioblastoma (GBM) is an aggressive primary brain cancer with few
   effective therapies. Stereotactic needle biopsies are routinely used
   for diagnosis; however, the feasibility and utility of investigative
   biopsies to monitor treatment response remains ill-defined. Here, we
   demonstrate the depth of data generation possible from routine
   stereotactic needle core biopsies and perform highly resolved
   multi-omics analyses, including single-cell RNA sequencing, spatial
   transcriptomics, metabolomics, proteomics, phosphoproteomics, T-cell
   clonotype analysis, and MHC Class I immunopeptidomics on standard
   biopsy tissue obtained intra-operatively. We also examine biopsies
   taken from different locations and provide a framework for measuring
   spatial and genomic heterogeneity. Finally, we investigate the utility
   of stereotactic biopsies as a method for generating patient-derived
   xenograft (PDX) models. Multimodal dataset integration highlights
   spatially mapped immune cell-associated metabolic pathways and
   validates inferred cell-cell ligand-receptor interactions. In
   conclusion, investigative biopsies provide data-rich insight into
   disease processes and may be useful in evaluating treatment responses.

   Subject terms: CNS cancer, Translational research, Surgical oncology,
   Cancer microenvironment
     __________________________________________________________________

   Stereotactic needle biopsies are currently performed for diagnosis only
   in glioblastoma; this study highlights the depth of information
   available through multi-omics analysis.

Introduction

   Glioblastoma (GBM) is an aggressive primary brain cancer with limited
   treatment options^[121]1. Despite numerous clinical trials to date, the
   standard of care remains maximal surgical resection followed by
   combination chemotherapy and radiotherapy^[122]2–[123]4. Apart from
   tumor-treating fields^[124]4, few novel therapies have been approved.
   Unfortunately, despite aggressive multimodal treatment, the prognosis
   remains poor. While immunotherapies in other solid cancers have shown
   significant improvements in overall survival, trials in GBM involving
   targeted and immunotherapies have been met with limited
   success^[125]5,[126]6.

   Standard clinical endpoints such as overall survival, progression-free
   survival or radiological evidence of progression are coarse markers of
   treatment response, and there is a strong need to improve the current
   understanding of how the tumor, host immune system and tumor
   microenvironment respond to therapies. Currently, there are few
   opportunities to obtain additional tissue samples to understand tumor
   responses during treatment other than at initial surgery in the form of
   either biopsy or resection. This lack of visibility into tumor and
   immune co-evolution under treatment is a fundamental limitation in our
   ability to develop and advance new GBM therapies.

   Stereotactic biopsies are a common neurosurgical procedure performed to
   obtain tissue diagnosis through basic histological and molecular
   evaluation, with a well-defined role and safety profile in the
   diagnosis of brain tumors and brain lesions^[127]7–[128]12. The trend,
   however, toward obtaining additional molecular data besides traditional
   hematoxylin and eosin (H&E), and immunohistochemical (IHC) staining has
   begun in earnest, and many academic medical centers have incorporated
   gene panel sequencing such as MSK-IMPACT^[129]13, DFCI Oncopanel or
   other solid tumor panels^[130]14,[131]15 from needle core biopsies as
   part of the routine clinical diagnostic workflow. The amount of
   clinical material obtained from these core biopsies is small (<50 mg),
   however the uni-directional cutting face of the biopsy needle permits
   multiple cores to be taken from the same site from different directions
   by rotating along its axis (Fig. [132]1a), enabling sampling from a
   microregion of the tumor itself.

Fig. 1. Intra-operative collection and coordination of sample core
distribution and analyses.

   [133]Fig. 1
   [134]Open in a new tab

   a Schematic diagram demonstrating multiple core sampling and proposed
   analysis methods. Note that the biopsy needle can rotate along its axis
   and perform multiple biopsies from the same fixed position. b MR
   imaging of patient undergoing resection of lesion. Left, T1 weighted
   imaging with contrast. Right, MR perfusion imaging demonstrating
   “hotspot” on the anteromedial portion of tumor, destined for resection.
   c Intra-operative neuro-navigated imaging demonstrating the trajectory
   and position of biopsy needle in 3 orthogonal planes and inline view. d
   Representative pictures of multiple biopsy cores (scale bar = 10 mm). e
   Core sampling locations geospatially mapped back to original MR imaging
   using 3D slicer software. QC quality control, PDX patient-derived
   xenografts. (a) was generated with Biorender.

   The role of investigative biopsies—i.e., biopsies where the primary
   goal is to collect systematic, detailed multimodal tumor-related
   information for the purposes of evaluating treatment response, has been
   under-explored^[135]16. Immune cell population flux and neoplastic cell
   responses to therapy vary over the disease course, and insight into
   this process may have important implications for decisions regarding
   treatment, disease progression and prognosis. The ability to track
   tumor and stromal populations and resolve cell-cell signaling pathways
   via spatial and single-cell technologies could potentially offer a more
   comprehensive overview of the tumor fitness landscape under different
   therapeutic modalities^[136]17.

   The application of multi-omics technologies to the GBM tumor
   microenvironment (TME) represents a powerful and unbiased approach to
   understand tumor-stromal interactions and tumor organizational
   structure^[137]18; single-cell RNA sequencing combined with spatial
   metabolomics analysis has shown that metabolites in the TME can shape
   the tumor-immune environment through purinergic signaling pathways,
   involving ATP and adenosine in tumor-microglial interaction via CD73
   and CD39 respectively^[138]19 or direct immunosuppressive effects via
   oncometabolites on T-cells^[139]20, these findings are supported by GBM
   single-cell tumor cell transcriptome analysis which suggests pathway
   convergence along neurodevelopmental and metabolic axes^[140]21.
   Phosphoproteomics and immunopeptidomics techniques have shown how the
   signaling pathways and the antigenic peptide repertoire can alter under
   treatment^[141]22–[142]25 thus shaping the systemic immune response to
   tumor. We therefore selected 6 assays as part of our multimodal
   analysis panel: (1) scRNA sequencing, (2) spatial transcriptomics, (3)
   spatial proteomics, (4) spatial metabolomics, (5) phosphoproteomics and
   (6) HLA Class I immunopeptidomics. The feasibility of adapting and
   applying multi-omics technologies to accommodate low input volumes in
   the “real world” setting is unknown, therefore the primary goal was to
   demonstrate safety and utility of applying multi-omics analysis to
   biopsy samples from a single patient. A secondary objective was to
   establish methods to measure intra- and inter-regional variance from
   biopsies within the tumor.

   Here we demonstrate the feasibility of this multi-omics approach on
   needle core biopsies from an individual patient. We developed an
   optimized workflow to obtain information regarding cellular
   composition, tissue architecture, cellular transcriptomes, the MHC
   Class I immunopeptidome and protein activation pathways from needle
   core biopsies collected during surgery. Secondary integration of
   multimodal measurements then allowed us to cross-validate findings and
   make inferences about cell-cell interactions and immune cell metabolic
   activity within the TME. The overall clinical and scientific impetus to
   develop multimodal analysis methods is driven by the need to understand
   disease progression and tumor-immune evolution over the course of
   treatment. The establishment of protocols that can facilitate this has
   important implications for clinical trials and disease monitoring going
   forward.

Results

Intra-operative spatially localized needle core biopsy sampling

   Multiple biopsies were taken from a patient undergoing craniotomy for
   resection of a recurrent GBM after obtaining informed consent
   (Fig. [143]1b). The patient had presented with tumor recurrence after
   completing standard-of-care first-line treatment (Supplementary
   Table [144]1). After re-opening of the craniotomy and exposure of the
   target resection area, we obtained sequential biopsies using a biopsy
   needle registered to a frameless stereotaxy system (BrainLab AG,
   Germany) in a region of tumor destined for resection along a single
   trajectory. We recorded the 3D coordinates using the intra-operative
   neuronavigation software at each depth and biopsy location (iPlan 3.0
   software, BrainLab AG, Germany). Samples were obtained from 4 quadrants
   where feasible before moving to a separate depth and obtaining
   additional samples. Five core needle biopsy samples (out of a total of
   seven cores obtained) were deemed of sufficient quality for downstream
   analysis (Fig. [145]1c–e). Sample volumes obtained from biopsy needle
   cores were in the range of 41–44 mg. Tissue cores were immediately
   placed on labeled saline-soaked non-adherent patties (Telfa®, American
   Surgical Company, USA). All cores were inspected visually to ensure the
   specimen was solid tissue and not blood clot. The two cores that had
   the most viable appearance were placed in tissue storage solution (MACS
   storage buffer, Miltenyi Biotec, Germany) and placed on ice for fresh
   processing (Core 2). The other cores were assigned according to the
   flowchart (Fig. [146]1a) and were either fixed in formalin (Core 1), or
   snap frozen using liquid nitrogen vapor (Cores 3–5). Cores were then
   transported on an expedited basis to partner institutions for
   processing and analysis. In addition, excess tissue from Core 1
   (formalin-fixed) was preserved for bulk T-cell receptor clonal analysis
   (Supplementary Fig. [147]1), excess frozen tissue from Core 4 (frozen)
   was used to perform bulk RNA sequencing.

Single-cell RNA sequencing of GBM biopsy cores

   Single-cell RNA sequencing (scRNA) technology permits high-resolution
   gene expression analysis of tissue samples, and has been published
   extensively in GBM showing glioma cell states^[148]26–[149]28, immune
   cell composition and oncogenic developmental programs^[150]29–[151]31.
   Nearly all studies utilized surgical resection samples, and it is
   unknown whether current tissue processing protocols for resection
   samples are suitable for needle biopsy cores^[152]32. We therefore
   adapted our existing protocols to low input volumes. Similar challenges
   were faced in a study on diffuse intrinsic pontine glioma (DIPG), where
   tissue samples are necessarily limited due to the extreme sensitivity
   of the brain structures involved. Although scRNAseq was shown to be
   feasible with low tissue input volumes^[153]33,[154]34, the number of
   cells sequenced was low (2458 cells from 6 patients)^[155]35. More
   recently, multi-location biopsy sampling was combined with
   scRNAseq^[156]36, and in this study, the investigators sequenced 6148
   cells from 73 biopsy sites in 13 patients. Overall, the combined
   sequenced cell yield per biopsy in these studies was in the region of
   80–400 cells per biopsy. We estimated, however, that we would need to
   aim for 10,000 cells to adequately detect tumor, immune and stromal
   populations, which would require a minimum starting quantity of 50,000
   viable cells. Single-cell dissociation of brain tumor tissue typically
   involves mechanical dissociation followed by enzymatic digestion,
   optionally followed by density gradient separation or myelin removal
   and red cell lysis^[157]30. This process can be standardized through
   automated tissue trituration and dissociation methods (GentleMACS
   system, Miltenyi Biotec, Germany). Automated tissue processing has the
   advantage of providing consistent conditions, however the instrument
   settings are not optimized for low volume samples. We therefore first
   performed a direct comparison between manual dissociation versus
   automated dissociation and showed that manual dissociation yielded
   approximately ten-fold as many live cells as the automated method, with
   the average yield per core of approximately 610,000 live cells after
   manual dissociation compared to 74,000 with automated methods
   (Supplementary Table [158]2).

   A recent study by Marsh et al. suggested that enzymatic dissociation is
   associated with an ex vivo activation signature (ExAM) which was absent
   in mechanical dissociation and is abrogated using protein translational
   inhibitors^[159]37. The original study was performed in normal mouse
   brains, therefore we sought to determine if this held true for human
   tumor tissue. To address this, we first performed single-cell
   dissociation and scRNAseq with and without translational inhibitors. We
   then repeated the experiment with and without enzymatic digestion, with
   and without inhibitors (Fig. [160]2a, b). Corroborating the report by
   Marsh et al., translational inhibitors clearly decreased the ExAM
   signature in the enzymatically dissociated samples, particularly in the
   myeloid population. Comparison of enzyme versus non-enzymatic
   mechanical conditions in simulated biopsy tissue, however, showed the
   ExAM signature was higher in the mechanically dissociated sample
   compared to the enzymatic condition (Fig. [161]2c). This observation
   suggests that non-enzymatic effects also contributed to the ExAM
   signature in the context of GBM tissue. Pure mechanical dissociation
   alone without enzymes also yielded a much lower number of viable cells,
   severely hindering the quality of the downstream analysis.

Fig. 2. Optimization of single-cell workflow for fresh GBM biopsy cores.

   [162]Fig. 2
   [163]Open in a new tab

   a Comparison of manual versus automated tissue dissociation and
   single-cell RNA sequencing with and without translation inhibitors (n =
   viable cell yield). b Ex vivo activation (ExAM) score in myeloid
   population with and without inhibitors under automated and manual
   dissociation conditions. Boxes are bound at the 25th and 75th
   percentile with the center line being at the median. The dots represent
   minima and maxima of the data. c UMAP showing cellular populations
   using mechanical only versus enzymatic dissociation and subsequent ExAM
   score in low volume simulated biopsy tissue. d UMAP showing cellular
   populations and proportions from actual patient biopsy sample using
   manual dissociation technique. Lower boxes represent higher-level
   annotations subpopulation analysis within the lymphocytic and myeloid
   populations. e Direct comparison of glioma-associated macrophage gene
   expression between frozen tissue derived single-cell nuc-seq versus
   fresh cell single-cell sequencing from a common originating sample
   showing significant reduction in key myeloid-related genes in frozen
   nuc-seq data. f Decision tree for single-cell processing workflow for
   stereotactic biopsy samples (Data from 4 different patient samples
   shown, a, b resection tissue, c simulated biopsy tissue, d actual
   patient biopsy tissue used for multimodal analysis, e resection tissue
   used for parallel scRNA and snRNA analysis).

   Another technical and logistical consideration was to determine whether
   scRNAseq should be performed using freshly dissociated single cells or
   frozen sample-derived single nuclei, with the rationale being that
   frozen samples are easier to transport, have a better chance of
   preserving nuclei in sufficient numbers for downstream single nuclei
   RNA sequencing (snRNA), and capture dissociation-susceptible cell types
   such as neurons. Previous studies comparing these two methods have
   shown loss of specific gene expression, particularly in microglial
   cells^[164]38. Given that one of our key study objectives was to study
   the immune myeloid cell composition, we performed pilot experiments in
   a separate patient with IDH-wildtype GBM and compared the yield of CD45
   positive population in both conditions, which corroborated the findings
   by Thrupp et al. We found important differences in gene expression
   between scRNA and snRNA, with key genes associated with
   glioma-associated myeloid cells (GAM) better represented in the scRNA
   data compared to snRNA (Fig. [165]2e). These genes included CD163,
   CX3CR1, AIF1, CD68, P2RY12 and SORL1. In view of our interest in tumor
   microenvironmental changes to treatment, we concluded that enzymatic
   fresh tissue dissociation without the use of translational inhibitors
   for scRNA would be preferable (Fig. [166]2f), and this final protocol
   was applied to the index patient in this study (Fig. [167]2d), which
   demonstrated neoplastic, myeloid and lymphocytic populations from
   surgical patient biopsy samples.

Spatial proteomic profiling of biopsy cores

   GBM is a histomorphologically heterogeneous tumor with regions of
   cellularity, infiltration, microvascular proliferation and
   necrosis^[168]39. Spatial tissue profiling technologies have emerged as
   a method of identifying intratumoral cellular neighborhoods to gain a
   phenotypic understanding of niche-specific molecular heterogeneity of
   GBM^[169]40. We performed highly multiplexed microscopy using
   Co-detection by indexing (CODEX, Akoya Biosciences, USA). A panel of 51
   antibodies was validated for tumor microenvironment-specific protein
   markers on FFPE tissue to identify immune cell types, tumor, vascular
   and stromal markers along with multiple functional markers
   (Fig. [170]3a, Table [171]1). The distance between specific immune
   cells from GFAP+ cells was evaluated by near neighbor analysis for the
   CD3+ T cells and CD68+ macrophages. Macrophages were more abundant and
   closer to the GFAP+ cells suggesting closer interaction between two
   cell types when compared to T-cells (Fig. [172]3b, table [173]1).
   Clustering of different cellular subpopulations showed an abundance of
   GFAP+ cells (33%) closely followed by myeloid cells (18%). Previous
   studies have shown that lymphocytic infiltration tends to be limited to
   the perivascular spaces; however, in our sample, there was a paucity of
   T cells and comprised only ~1% of all cell type population
   (Fig. [174]3c) with no clear relationship to the vasculature. Further,
   a total of 15 cellular subtypes were identified (Supplementary
   Fig. [175]2). These cell clusters were further interrogated to
   characterize cellular neighborhoods (CN) that were then re-mapped to
   the original image to reveal the spatial organization of CN
   (Fig. [176]3d, e). A total of six different CN were identified, with
   tumor cells (CN1) being most frequent and surrounded by proliferative
   tumor and myeloid cell-rich neighborhood (CN2). Additional CNs
   including vasculature associated immune cells (CN3), macrophage and
   microglial cells (CN4), myeloid and B cells (CN5), and GFAP+ cells in
   association with reactive astrocytes (CN6) were also observed. Defining
   spatial organization can provide insight into interactions of the
   different cell subtypes within the tumor ecosystem and the potential
   impact of drug treatment on tumor evolution and cell-cell
   cross-communication.

Fig. 3. CODEX Analysis of core biopsy sample.

   [177]Fig. 3
   [178]Open in a new tab

   a CODEX panel and representative image showing all biomarker staining.
   b Nearest Neighbor analysis showing distance of CD3+ T cells and CD68+
   macrophages to the GFAP+ cells. c Pie chart showing different
   subpopulations present within the pilot tissue sample. d Spatial
   location of Cellular Neighborhood (CN) as seen on the biopsy tissue.
   Stack plot represents the frequency of different CN. e Neighborhood
   analysis shows different cellular clusters contributing to form
   different CN.

Cell signaling with proteomics and phosphoproteomics

   Protein tyrosine phosphorylation (pTyr) can regulate protein activity,
   interactions, localization, and stability. Signaling networks mediated
   by pTyr control a variety of cellular processes in response to
   intracellular and extracellular (e.g., environmental) perturbations.
   Aberrant pTyr signaling is known to be oncogenic, with many tyrosine
   kinases identified as potent oncogenes, and conversely, protein
   tyrosine phosphatases are known to be tumor suppressors. Moreover,
   tyrosine phosphorylation regulates adaptive and innate immune cell
   activity through cell surface receptors, JAK/STAT signaling, and immune
   checkpoint proteins such as SIRPα and the Siglecs. Receptor Tyrosine
   Kinase (RTK) pathways such as Epidermal Growth Factor (EGF) and
   Platelet-Derived Growth Factor (PDGF) alpha and beta family of
   associated receptors are frequently aberrantly activated in GBM and
   result in an oncogenic phenotype^[179]41,[180]42. STAT3 signaling has
   also been implicated in GAM polarization in GBM and recruitment of
   regulatory T-cells^[181]43,[182]44. We utilized mass spectrometry-based
   phosphoproteomics to identify tyrosine phosphorylated proteins in the
   patient tumor and capture the tumor cell-intrinsic and immune signaling
   components. Several 10-micron thick sections from Core 5 were lysed and
   processed to yield 127 µg of tryptic peptides, which were subsequently
   enriched for pTyr-containing peptides by immunoprecipitation and
   phosphopeptide enrichment. Enriched pTyr peptides were then analyzed by
   liquid chromatography-tandem mass spectrometry (LC-MS/MS) on an
   Exploris 480 Orbitrap mass spectrometer. At 1% FDR, the analysis
   yielded 288 unique phosphopeptides containing 270 pTyr, 40 pThr, and 60
   pSer phosphorylation sites. With a more stringent ions score threshold
   (≥20), the analysis identified 258 unique phosphopeptides, comprising
   248 pTyr, 24 pThr, and 36 pSer phosphorylations on proteins ranging
   from receptor tyrosine kinases such as EGFR and their downstream
   signaling networks, including the ERK mitogen-activated protein kinases
   (MAPKs) to immune cell-specific proteins including CD226 (Y322), CD84
   (Y279, Y296, Y316), and Siglec-9 (Y456) (Fig. [183]4a, Supplementary
   Table [184]3). To evaluate the respective abundance of the various
   phosphorylation signals, the area under the curve (AUC) was quantified
   for the chromatographic elution profile for each phosphopeptide. This
   analysis highlighted the most abundant signals, including activation
   loop phosphorylation sites on well-known oncogenes such as the ERK
   MAPKs and Src-family kinases (SFKs) (Lyn/Hck and Yes/Src/Fyn/Lck) as
   well as activating sites on p38 MAPK and the STAT3 transcription
   factor, which has been implicated in driving tumor progression in GBM.
   Although EGFR phosphorylation was detected, only a single site (Y1173)
   was identified at relatively low abundance compared to SFKs or ERK
   MAPKs, suggesting that this tumor was more likely to be driven by SFKs
   than EGFR signaling. In addition to the oncogenic signals, a large
   number of pTyr sites were indicative of innate immune infiltration into
   the tumor, including sites on the immunoglobulin epsilon receptor
   (FcER) and immunoglobulin gamma receptor (FcGR), as well as sites on
   the Tyro receptor binding protein (TYROBP), Siglec-9 and PECAM1. For
   global phosphoproteomic (3 pTyr, 281 pThr, 1614 pSer) and proteomic
   (2746 unique proteins) analysis, supernatant from
   pTyr-immunoprecipitation was fractionated using a high pH gradient
   (Fig. [185]4a, Supplementary Table [186]3). The quantified proteins and
   phosphoproteins were involved in a plethora of biological pathways
   (e.g., apoptosis, translation, neuron-projection development,
   metabolism, etc.) that cover various hallmarks of cancer
   (Fig. [187]4b). The phosphorylation status of critical immune
   checkpoint proteins such as PD-L1 (S283) and VISTA (S235) showed that
   these receptors were present and engaged in the tissue. Additionally,
   68 proteins and 38 phosphorylation sites implicated in several
   metabolic pathways (e.g., pyruvate, ATP, NADP, and fatty acid metabolic
   pathways) were quantified in the biopsied sample. Taken together, we
   demonstrated that comprehensive phosphoproteomic and proteomic analysis
   could be done on a small amount of biopsy to shed light on the complex
   signaling networks formed by tumor and the surrounding microenvironment
   (Fig. [188]4b).

Fig. 4. The depth of proteomic, phosphoproteomic, and immunopeptidomic
analysis on GBM core biopsy.

   [189]Fig. 4
   [190]Open in a new tab

   a Unique proteins and phosphoproteins quantified are ranked by their
   abundance. The names of select proteins and phosphoproteins associated
   with tumor progression and immune infiltration are shown. b Gene
   ontology terms associated with proteins and phosphoproteins that were
   quantified. Groups that had intersection sizes above 3 are shown. c The
   number of 8 to 11-mer peptides identified by immunopeptidomic analysis.
   d The number of proteins that gave rise to MHC peptides. e Motif
   analysis of the quantified immunopeptidome using Gibbs clustering. f
   MHC peptides quantified from the core biopsy were compared to the HLA
   Ligand Atlas of the benign human brain. Previously reported
   cancer-testis antigens (CTNNA2, SPAG17) and predicted CTAs (CTCFL) are
   labeled along the immunopeptides that were ranked based on their
   abundance.

Immunopeptidomics on GBM biopsy cores

   Immunological control of tumor growth is based upon the recognition of
   MHC Class I and Class II peptides by CD8+ and CD4+ T cells,
   respectively. Although immune checkpoint inhibitors have been
   unsuccessful in GBM clinical trials, results from targeted
   immunotherapies with MHC Class I neoantigens or tumor-associated
   antigens have suggested potential therapeutic benefit in selected
   patients, especially if the treatment can be matched to the
   immunopeptide expression profile for a given patient^[191]45. To
   determine the MHC Class I profile for the current patient, a frozen
   core (Core 5) was subjected to mass spectrometry-based immunopeptidomic
   analysis. Following homogenization and cell lysis, which yielded
   1.95 mg of lysate, MHC class I peptide complexes were isolated using
   W6/32 antibody. MHC I peptides were then released by low pH and
   enriched by molecular weight cutoff filter prior to the analysis by
   LC-MS/MS. This analysis led to the identification of 2915 MHC class I
   peptides, of which 1847 were 9-mers and 605 were 10-mers
   (Fig. [192]4c). Additionally, VIM, PTPRZ1, PRPF8, and SPTBN1 were
   highly represented in the immunopeptidome, contributing 16, 9, 8, and 8
   epitopes, respectively, suggesting their high expression in the biopsy
   core (Fig. [193]4d). These proteins have been previously implicated in
   GBM invasion and progression^[194]46–[195]48. Gibbs clustering of the
   data highlighted four MHC peptide motifs (Fig. [196]4e). Of the 2915
   peptides detected, 199 were found within the HLA Ligand Atlas, a
   comprehensive collection of HLA ligands presented on benign tissues,
   suggesting that these peptides may represent those that are enriched in
   malignant tissues^[197]49. In line with this, previously characterized
   and predicted cancer-testis antigens (CTA)^[198]50,[199]51 were
   identified in the patient-specific immunopeptidome, again highlighting
   the ability of this technique to detect antigens presented in GBM
   (Fig. [200]4f).

Spatial metabolite imaging

   The integration of H&E, t-CyCIF and MALDI-MSI data of consecutive
   biopsy tissue sections reveals insights into the spatial context of
   immune cell metabolism in GBM (Fig. [201]5). Adenosine triphosphate
   (ATP), an essential compound in metabolism and pathway
   signaling^[202]52,[203]53 is generated by aerobic glycolysis and the
   tricarboxylic acid (TCA) cycle in GBM^[204]54. Our analysis
   demonstrates ATP-positive pixels acquired with MALDI-MSI in CD45-
   regions shown in t-CyCIF analysis (Fig. [205]5a). The extracellular
   breakdown of ATP into adenosine, called purinergic signaling, plays a
   role in pathogenesis and enhanced chemoresistance in GBM^[206]55. High
   CD73 expression, the ectonuclease responsible for the conversion of AMP
   to adenosine, is found on CD45+ immune cells that suppress immune
   activity^[207]56 and results in decreased levels of ATP. Conversely,
   the catabolism of ATP is exceedingly slow in glioma cells, resulting in
   the accumulation of ATP in proximity to tumor cells^[208]57. On the
   other hand, malate and linoleic acid metabolism are enriched in regions
   of CD45+. Although the exact mechanisms in GBM are unclear, high levels
   of either malate or linoleic acid lead to metabolic reprogramming in
   the mitochondria^[209]58,[210]59. Linoleic acid, a highly consumed
   polyunsaturated fatty acid, has been shown to improve the antitumor
   function of CD8+ T cells through metabolic reprogramming and
   facilitates superior cytotoxic function^[211]60. Integration of spatial
   proteomics and spatial metabolomics show enriched metabolites in CD45+
   and CD45- regions (Fig. [212]5b) that correspond to metabolic pathways
   throughout the biopsy sample (Fig. [213]5c, d). CD45- regions are
   enriched for ubiquinone and other terpenoid-quinone biosynthesis,
   caffeine metabolism, arginine and proline metabolism, aminoacyl-tRNA
   biosynthesis, tryptophan and tyrosine metabolism. Although GBMs
   primarily rely on glucose for metabolic fuel, they also utilize amino
   acids such as glutamine and glutamate^[214]61. Reprogrammed glutamine
   is essential for the proliferation and survival of cancer cells, which
   also yields an increase in amino acids, TCA cycle intermediates and
   fatty acid biosynthesis^[215]62. Our results consistently demonstrate
   enriched amino acid pathways in the CD45+ regions, such as taurine and
   hypotaurine metabolism, valine, leucine and isoleucine biosynthesis,
   arginine biosynthesis, linoleic metabolism and D-glutamine and
   D-glutamate metabolism.

Fig. 5. Integrated spatial proteomics and metabolomics on GBM biopsy core.

   [216]Fig. 5
   [217]Open in a new tab

   a Representative CyCIF zoom-ins displaying CD45, OLIG2, CD14 and GFAP
   (scale bar = 50 µm). b Hematoxylin and Eosin staining, and
   corresponding spatial distribution of immune cells (CD45+), Adenosine
   triphosphate (ATP), Malate, and Linoleic acid metabolism enrichment
   with immune cells (CD45+) overlayed. The distribution of ATP is
   anticorrelated, while Malate and Linoleic acid metabolism show strong
   spatial colocalization to immune cells. c Diagram summarizing the
   integration of CyCIF and MALDI-MSI of consecutive tissue sections. This
   includes image registration, Pearson’s correlation of the number of
   CD45+ cells per pixel to all metabolites, and pathway enrichment
   analysis. d Comparison of metabolism between immune and non-immune
   regions of the tissue. e Pathways enriched in non-immune regions of the
   tissue. f Pathways enriched in immune regions of the tissue.

Integrative analysis of single patient-derived multi-omics data

   A key advantage of multimodal analysis is the provision of orthogonal
   data from the sampled tumor environment. We therefore performed
   cross-integration of multi-omics data to gain insight into the salient
   features of the TME and validate existing approaches to cell-cell
   interaction inference. Preliminary integrative data analyses were
   performed on the datasets: Mapping of 3D spatial sampling coordinates,
   single-cell RNA sequencing and metabolic profiling of tumor and immune
   populations, single-cell and bulk RNA sequencing and immunopeptidome,
   spatial proteomics and transcriptomics.

Core sampling locations

   MR coordinates were obtained via neuro-navigated biopsy needles, with
   3D coordinates recorded intra-operatively then re-mapped onto patient’s
   pre-operative scan to confirm correlate sampling locations to their
   respective assay (Fig. [218]6a, [219]b).

Fig. 6. Integrative analysis of single patient-derived multi-omics data.

   [220]Fig. 6
   [221]Open in a new tab

   a Core sampling locations and proposed analysis methods. b Comparison
   of pathway enrichment analysis of immune vs non-immune cells using
   spatial metabolomics and scRNAseq. r = Pearson correlation coefficient.
   c Correlation between the abundance of MHC peptides and their matching
   transcripts and proteins. R = Pearson correlation coefficient. d
   Correlation of abundance of MHC peptides with corresponding cellular
   subsets as identified by scRNAseq. e Comparison of CODEX protein data
   to GeoMx-WTA data using corresponding regions of interest. Multiplex
   immunofluorescence was performed on the biopsied section showing
   selected ROIs (left panel). Protein expression for phenotypic markers
   by CODEX (middle) and immune cell deconvolution from GeoMx-WTA (right)
   showing the distribution of cell type scores of immune cell subtypes in
   the selected ROIs. f Rank plot displaying the communication probability
   of all predicted cell-cell interactions and ligand-receptor
   interactions as determined by CellChat using scRNAseq data. Points
   highlighted represent the SPP1-CD44 ligand-receptor pair. The bottom
   stripe graph indicates ligand-receptor pairs with activated receptors
   according to phosphoproteomics data. The vertical line indicates the
   transition between highly probable (green) and unlikely (red)
   interactions. g Spatial neighboring frequency analysis using CyCIF.
   Highly specific markers identify cell types: CD8 (T-cells), TMEM119
   (Microglia), CD14 (Monocytes), SOX2 (MES-like), SOX9 (NPC-like), and
   OLIG2 (OPC-like). Cells positive for both OPN+ and CD44+ are considered
   targets, and their neighboring frequency to presumed sources is
   quantified in a network connecting all cells with centroids within 50
   μm of each other. The statistical significance of these frequencies is
   assessed using a Monte Carlo approach by permuting the cell labels 1000
   times. h Representative CyCIF zoom-ins showing cell-cell interactions
   involving OPN and CD44: Microglia-CD8 T-cells (top), Microglia-MES-like
   (middle), and Microglia-Monocyte (bottom). i Correlation plot between
   Communication Probability (scRNAseq + CellChat) and the statistical
   significance of the neighboring frequency analysis in CyCIF. Spearman’s
   ρ = 0.02, p = 0.02. j Main interactions identified by CyCIF neighboring
   frequency analysis within the context of the whole tissue section. DSP
   digital spatial profiling.

Cross-comparison of spatial proteomics, metabolomics and scRNAseq

   Comparison of spatial metabolomics with single-cell RNA sequencing data
   was performed by metabolic pathway enrichment analysis on both
   datasets. We computed the enrichment in immune and non-immune cell
   populations and found a moderate correlation (Pearson’s r = 0.65,
   p = 0.03) (Fig. [222]6b) between the enrichment ratio of selected
   metabolic pathways. This result indicates agreement between the
   transcriptomic and metabolomic analyses performed. The top co-nominated
   pathways for immune cells included D-glutamate metabolism,
   leucine/isoleucine biosynthesis and linoleic acid metabolism
   (Fig. [223]6b).

Correlation between MHC class I immunopeptidome and bulk and single-cell RNA
expression

   The rank-ordered abundance of MHC peptides assessed by
   immunopeptidomics did not positively correlate with that of
   corresponding transcripts (Fig. [224]6c, left) or proteins
   (Fig. [225]6c, right), highlighting the need for immunopeptidomic
   analysis to directly quantify tumor antigen presentation. Single-cell
   RNA sequencing data generated from tissue from the same procedure
   allowed us to delve deeper into the different cell types. The abundance
   of MHC peptides from immunopeptidomics again did not show significant
   correlation with the corresponding standardized average gene expression
   across the four major cell types (Myeloid, Lymphoid, Neoplastic, and
   Vascular cells) detected on scRNAseq (Fig. [226]6d). Although none of
   the cell types showed a strong correlation with the detected
   immunopeptidome, transcript expression of the ~1900 source genes did
   skew toward the neoplastic population, suggesting that the MHC Class I
   immunopeptidomics data are enriched for tumor-associated proteins;
   however, transcript expression remains a poor predictor of peptide
   abundance.

Cross-comparison of spatial proteomics profiling with spatial transcriptomics

   We used the CODEX and GeoMx -WTA platforms and compared immune cell
   distribution obtained from similar regions of interest (ROI)
   (Fig. [227]6e, left panel) from non-adjacent sections of a single
   biopsy for quantitative comparison across spatial proteomics and
   transcriptomics data. Expression of phenotypic markers from CODEX assay
   (Fig. [228]6e, middle panel) was comparable with GeoMx-deconvolution
   data (Fig. [229]6e, right panel) and indicated comparable immune cell
   infiltration in the selected ROIs with the predominance of myeloid
   cells and fewer B and T lymphocytes.

Validation of cell-cell interactions based on multi-omics integrative
analysis

   Tumor-stromal interactions are thought to play a significant role in
   immune suppression and disease progression. Techniques to investigate
   cell-cell interactions have been developed which leverage knowledge of
   known ligand-receptor pairs (L-R) and gene co-expression between
   different cell types to determine the probability of
   interaction^[230]63. Since phosphoproteomic data obtained from the core
   biopsies provides information regarding activated signaling
   pathways^[231]64, we therefore performed integrative analysis by
   filtering inferred putative L-R interactive pairs based on single-cell
   RNA sequencing data based on the abundance of phosphorylated receptors
   and identified osteopontin (SPP1) and CD44 as a top-ranked L-R pair.
   SPP1 is known to be expressed by tumor-associated myeloid cells and has
   been reported to play a role in impairing T-cell responses^[232]65 and
   tumor cell crosstalk^[233]66,[234]67. We next examined the inferred
   cell types where this interaction was predicted to occur based on
   Cellchat analysis, then validated these cell-cell interactions by
   performing CyCIF staining for SPP1, CD44 and predicted cell type
   markers such as TMEM119, CD3 and SOX2 on the same biopsy sample.
   Analysis of cell type and colocalization of SPP1:CD44 showed good
   correspondence (Fig. [235]6h). In conclusion, this layering of
   orthogonal analysis methods affirmed the role that SPP1:CD44 plays in
   the TME and confirmed some of the predicted cell-type interactions
   identified through single-cell analysis.

Estimated intra-core, intra-positional and inter-positional variance between
needle core biopsy samples

   A critical unresolved issue when performing multiple needle core
   biopsies is the underlying sampling variance^[236]68. GBM is recognized
   as a highly heterogenous tumor^[237]69,[238]70. We therefore sought to
   develop a quantitative metric of heterogeneity to compare samples from
   a given region. Specifically, we assessed the variance between (1)
   intra-regional core biopsies—i.e., multiple biopsies taken from the
   same trajectory location within the tumor, and (2) inter-regional core
   biopsies—differences in sample content from distinct trajectories (and
   hence by extension, regions) of the tumor itself. Since needle rotation
   along a 360° plane angle allows capture of up to 6 biopsies, each
   “microregion” of the tumor sampled represents a tissue compartment of
   approximately ~1 cm^3. A simple proposed criterium for determining
   intra-regional sample homogeneity therefore would be that
   intra-regional sample variance should not exceed inter-regional
   variance, and violation of this rule would signify the sample as an
   outlier.

   To apply this principle to a real-world situation, we performed needle
   core biopsies from spatially distinct locations in a separate patient
   with recurrent GBM undergoing surgical resection of recurrent tumor
   (Table [239]S1). A total of 7 cores were obtained from 3 needle passes
   (referred to as regions) of an en-bloc resected recurrent GBM specimen
   (Fig. [240]7a).

Fig. 7. Multi-regional sampling variance analysis.

   [241]Fig. 7
   [242]Open in a new tab

   a A recurrent GBM is resected en-bloc and samples are taken from 3
   separate regions of the tumor (labeled 1–3). b CNV analysis shows low
   tumor purity in samples 1a and 1b, PBMC control sample is also included
   in the last track. c Phylogenetic analysis demonstrating common shared
   events such as chromosome 7 gain and 10 loss followed by branching
   differences separated by sample of origin. d Diagram representing the
   microregions sampled with each needle core biopsy. e Pre-operative T1
   contrast-enhanced and T2 image showing the location of the lesion. f
   Representative H&E sections from each core (top) with corresponding
   segmentation (bottom) performed with a machine learning model trained
   on contours by a neuropathologist (TB). g The segmentation maps are
   randomly subsampled multiple times to derive a distribution of sampled
   imaging features from each sample. This distribution is visualized
   using PCA and statistically compared across microregions using linear
   mixed effects (LME). Microregion 2 is statistically distinct from the
   other two microregions. Sample 1c is an outlier within microregion 1
   due to an increase in necrotic features. (d) was generated with
   Biorender.

   A global view of tumor purity was initially evaluated by determining
   the tumor fraction leading to integer copy number states for gain of
   chromosome 7 and loss of chromosome 10. Both chromosomal alterations
   are prevalent in GBM and were presumed to be clonal. Samples 1a and 1b
   were noted to have <5% tumor purity and were therefore excluded from
   the downstream analysis, whereas the remainder of the samples had tumor
   purity ranging from 26 to 88% purity (Fig. [243]7b). Peripheral blood
   mononuclear cells (PBMCs) from the same patient were included as normal
   control. Phylogenetic analysis demonstrated that most events were
   shared across all samples, such as gains on chromosomes 7, 19, 20
   (including focal amplification of EGFR), as well as loss on 10 and a
   focal loss on 9p around the CDKN2A locus. Samples had one to two
   additional subclonal events, such as a focal gain on chromosome 11p
   (shared between 1c, 2b, 3a and 3b) and sample-specific events such as
   loss of 15 in 2a, gain of 1q in 2b and loss on 6q in 1c. Differences
   within trajectories were pronounced between 2a and 2b, whereas 3a and
   3b appear indistinguishable at the copy number level. In this case, the
   overall similarities between samples overshadowed the differences, as
   samples consistently exhibited a core karyotype comprising alterations
   on 5 chromosomes, accompanied by a maximum of 1–2 subclonal
   alterations.

   To evaluate tissue state variance between cores, each core sample was
   reviewed by an expert neuropathologist (TB), and histological features
   identified were then used to inform a machine learning classifier to
   quantify tissue variance. We then adopted a Monte Carlo subsampling
   method to simulate randomly sampling of a subregion of a given core to
   derive a distribution of sampled imaging features from each biopsy
   core, we could then compare intra-regional an inter-regional
   distributions through principle component analysis (PCA). These results
   demonstrate how closely related intra-regional samples are compared to
   inter-regional samples (Fig. [244]7f, g, Supplementary Fig. [245]4),
   and also indicate potential reasons for outliers, for example, core 1c
   sits further apart from the other cores due to an increase in necrotic
   features (Fig. [246]7f, g, Supplementary Fig. [247]4). Taken together,
   tissue subsampling and image analysis provides a statistical and
   quantitative method of both determining tumor heterogeneity using H&E,
   but also tissue quality control to identify potential outliers which
   can be either excluded or analyzed separately.

Generation of patient-derived xenografts using needle core biopsy samples

   Advances in patient-derived models of GBM have aided our ability to
   test therapeutics and serve as a platform for the emerging field of
   precision medicine. Having models that are matched to multi-omics data
   from the same needle biopsy would allow for functional integration with
   multi-omics results and further mechanistic insights into resistance.
   Prior studies have shown that small tissue samples from several non-CNS
   cancer types are sufficient for generating enough cells to create
   patient models such as 2D cell lines, 3D organoids/neurospheres or in
   vivo PDX models^[248]71–[249]73. However, the specific success rates
   for such models in the challenging context of necrotic/reactive tissue
   available in recurrent glioma stereotactic biopsies are not known.

   To conceptually evaluate the feasibility of patient model generation
   along with multimodal data generation from a single neurosurgical
   needle biopsy procedure, we retrospectively examined data from model
   generation we had performed in parallel to a recently reported clinical
   trial ([250]NCT03152318) (Fig. [251]8a)^[252]74. Thirty patients
   underwent stereotactic biopsy to confirm pathologically the presence of
   high-grade glioma (HGG, WHO grade 3 or 4) on frozen section, and up to
   6 core biopsy samples per patient were obtained. Patients included 28
   IDH-wildtype glioblastoma (WHO grade 4), 2 IDH-mutant astrocytoma (1
   WHO grade 3 and 1 grade 4) and had the expected spectrum of genomic and
   clinical characteristics for these diagnoses.
   Neuropathologist-estimated composition on H&E of each biopsy
   (Fig. [253]8b) revealed a variable degree of tumor and necrosis with
   tumor content ranging from 5 to 80% and all cores having 30% necrosis
   or less. Transcriptional silencing of the
   O6-methylguanine-DNA-methyltransferase (MGMT) gene was inconsistent
   among the group, with 9 patients (30%) demonstrating partial or full
   MGMT promoter methylation. Next-generation targeted exome sequencing of
   these gliomas demonstrated a broad landscape of canonical genomic
   alterations, with EGFR and CDK6 copy number gain as well as copy number
   loss of PTEN and CDKN2A (Fig. [254]8c).

Fig. 8. In vivo and in vitro modeling of glioblastoma from stereotactic
needle core biopsy.

   [255]Fig. 8
   [256]Open in a new tab

   a Stereotactic needle core biopsies were obtained from 30 HGG patients,
   QC metrics determined, and used in 26 PDX and 22 PDCL attempts. b
   Representative tissue sections from two patients demonstrate the small
   size of stereotactic needle biopsy cores and microscopic glioma
   histology, as well as clear evidence of PDX glioma histology via
   positive immunohistochemical staining for clinically reliable glioma
   markers. c 42% (11/26) of PDX attempts and 23% (5/22) of PDCL attempts
   were successful for a cohort of patients that demonstrated a diversity
   of histopathological, clinical, and tumor genomic characteristics. d
   Neuropathologist-estimated tumor percentage and measured biopsy mass,
   and cell viability in biopsy samples with correlations to model
   generation success (n = biological replicates). Boxes are bound at the
   25th and 75th percentile with the center line being at the median. The
   whiskers for plots are the minima and maxima of the data. The
   statistical test used was an unpaired, two-tailed t-test with
   significance determined as p < 0.05. HGG high-grade glioma, PDX
   patient-derived xenograft, PDCL patient-derived cell line. (a) was
   generated with Biorender.

   Tissue from a single core biopsy was used to attempt model generation.
   Cells were either directly implanted into the brain orthotopically
   (PDX) or concurrently cultured in vitro as 3D neurosphere
   patient-derived cell lines (PDCL) as previously described in defined
   stem cell media with supplemented with bFGF and EGF as previously
   described^[257]75. Cores were dissociated to single cells and
   quantified as to the tumor percentage, core biopsy mass, and the viable
   cell count ultimately isolated. The cores varied greatly in the size,
   content, and quality of tumor (Fig. [258]8d). Tumor percentage was a
   mean of 48% in the cores (range 5–80%). Mean biopsy mass was 39 mg
   (range 6.8–191 mg). Mean needle biopsy viable cell count
   post-processing was 1.5 million cells (range 0.13–21 million) with a
   mean viability of 79% (range 35–100). Overall the core biopsies yielded
   sufficient cells therefore to perform 7 PDX or PDCL attempts per core
   using 100,000 cells injection or T25 flask. We generally prioritized
   attempting two PDX (mean of 1.9 mice injected per core biopsy) (1–4
   mice per core) with excess cells going to PDCL attempts. Out of 30
   patients with needle biopsy core procedures, we were able to attempt
   both orthotopic PDX and PDCL (n = 18 patients), PDX only (n = 8) and
   PDCL only (n = 4). In total, we performed 26 PDX attempts and 22 PDCL
   attempts.

   Growth-verified successful PDXs (defined as histologically verified
   tumors within the brain) were generated from 11/26 patients (42.3%),
   which is analogous to the maximal success rates reported^[259]76.
   Successful PDCLs were generated from only 5/22 patients (22.7%) which
   was lower than published rates for in vitro neurosphere generation from
   larger samples^[260]77. Growth-verified models were obtained even from
   samples with lower tumor content or significant necrosis with the
   specific tumor content or necrosis within the range. Successful PDX
   models were generated from tumors ranging between an estimated 10–70%
   tumor content and 0–30% necrosis. Successful PDCL models were generated
   from tumors ranging between 30–70% tumor content and 0–20% necrosis.

   Correlation of the tumor percentage in patient core biopsies (H&E
   estimates of nuclei percent), biopsy mass, viable cell count was
   examined to identify the most important predictors of success rates
   (Fig. [261]8d). For PDX, there was no significant correlation of
   success with these three variables, but viable cell count showed near
   significance (p = 0.07). For PDCL attempts, the biopsy mass and the
   viable cell count were both significantly correlated with success
   (p = 0.01 and 0.03, respectively), while tumor percentage was not
   significant.

   All successful models demonstrated clear features of diffuse high-grade
   glioma histology by hematoxylin & eosin (H&E) staining as well as
   lineage verification by immunohistochemistry for glioma-defining
   markers including SOX2, Nestin, and GFAP and genomic sequencing. No
   donor patient TME was identified in PDX or PDCL by genomics or
   histology review (Fig. [262]8b).

   In summary, both in vitro and in vivo models were readily generated
   from stereotactic core biopsies but with higher success rates for in
   vivo PDX than PDCL. Tumor content did not significantly influence model
   generation success, but the total viable cell number isolated seemed to
   be the most important predictor, though the small size of our cohort
   warrants further investigation into this topic.

Discussion

   GBM remains the most common and lethal primary brain cancer, with
   minimal incremental improvements in patient outcomes despite over a
   century of research and trials. It is now understood that GBM has a
   highly immunosuppressive, immune “cold” tumor microenvironment,
   impervious to current therapeutic efforts. The lack of progress in GBM
   therapeutics can be partially attributed to the lack of understanding
   of disease response and adaptation under therapeutic pressure. This
   lack of visibility into how cellular states alter and respond to
   cytotoxic stress and how tumor-immune and tumor-stromal interactions
   interact to generate an immunosuppressive tumor milieu that favors
   progression and resists treatment has handicapped efforts to design
   more rational therapies.

   Diagnostic needle core biopsies have been the mainstay of neurosurgical
   diagnosis for decades and the indications and safety profile are well
   established^[263]8. However, the use of biopsies to gain additional
   information pertinent to a patient’s disease course beyond diagnosis,
   what we term investigative biopsies, has not been established. The
   central premise of this study is that standard needle biopsy cores
   contain a wealth of information which at present is underutilized.
   Possible reasons for this include concerns about undue risk to
   patients, technical challenges in obtaining high-resolution data, and
   the lack of experience in generating multi-omics data from low-volume
   samples. These reasons result in a lack of rationale and interest in
   investigative sampling, despite a clear need.

   The ability to generate high-resolution, data-rich tissue profiling
   information to gain an understanding of a patient’s disease course
   potentially alters the risk-benefit calculus for investigative
   biopsies. Extending this concept further, longitudinal investigative
   biopsies may permit unprecedented insight into tumor and immune
   responsiveness to therapy. Underpinning this entire concept, however,
   is the need to establish the feasibility of obtaining high quality
   multimodal measurements from a single patient in the same procedure.

   Large-scale proteogenomic approaches integrating genomic data with
   proteomics and post-translational modifications such as protein
   phosphorylation such as those by the National Cancer Institute Clinical
   Proteomic Tumor Analysis Consortium (CPTAC) have revealed new
   biological mechanisms and ushered in new integrative computational
   methods^[264]78–[265]84. The advent of single-cell profiling has the
   potential to refine these associations further with corresponding
   spatial correlations at the intra-operative macroscopic and microscopic
   tissue level scales. Multi-omics analysis of GBM tumors has previously
   identified core transcription factors associated with mesenchymal
   transformation in GBM^[266]85. Proteogenomic analyses have found
   cellular stress and growth factor networks regulating the EGFR
   signaling pathway^[267]86,[268]87. A similar stress response phenomenon
   is also seen when integrating metabolic and proteomic profiles^[269]88.
   Recently, Dekker et al. performed gene expression, proteomics and
   phosphoproteomics analysis on 8 primary and recurrent sample pairs and
   reported that differences found at the transcriptional level were not
   seen at the proteomic level despite more pronounced differences between
   primary and recurrent samples seen in the proteomics data^[270]89. This
   data reinforces the view that proteomics data provides orthogonal
   information that can be integrated into the overall picture of the TME.
   It should be noted, however, that in all studies, the specimen type was
   not specifically sourced from needle core biopsies, which are an
   abundant and accessible specimen type in routine clinical practice. The
   ability to perform integrative multi-omics analysis on needle core
   biopsies alone therefore has significant implications on the
   availability of tissue specimens for study. Single-cell analysis of
   needle core biopsies in brain tumors has been attempted in the past;
   however, the number of reported cell sequences has been
   low^[271]35,[272]36, with the number of cells per biopsy ranging
   between 70 and 400 cells per core. With advances in tissue processing
   and single-cell encapsulation, we have shown that viable cells in the
   region of 250,000–600,000 cells can be generated from 2 cores, which
   represents a difference of several orders of magnitude, and well within
   the range of standard input requirements for single-cell RNA
   sequencing. Remarkably, these cellular yields were obtained from
   previously treated, recurrent tissue, where one would expect lower cell
   yields, affirming the robustness of the protocol, standing in contrast
   to previous reports, which were obtained from newly diagnosed
   cases^[273]35,[274]36. Besides scRNA, we also demonstrate an
   unprecedented breadth of information that can be obtained from standard
   needle biopsy cores through orthogonal assays. Although multi-core
   sampling is part of routine surgical practice, the application of
   systematic sampling for targeted multi-omics analysis combined with
   precise geotagging of sampling coordinates is not routine and provides
   additional macroscopic spatial information (Fig. [275]6a). The analysis
   techniques were selected based on the organizing principle of
   tumor-immune co-evolution. Modalities were specifically selected to
   resolve cell type composition and cell states, cell-cell interactions
   and spatial cellular profiles at both the gene expression and proteomic
   level, as well as intratumoral TCR clonal frequencies sequences
   (Supplementary Fig. [276]1), buttressed by orthogonal phosphoproteomics
   data, spatial metabolomics data and intratumoral antigen diversity via
   MHC Class I immunopeptidomics, all obtained via minor modifications to
   a standard neurosurgical procedure. The data generated is of high
   quality and relevant, and the results of this study stand as a proof of
   principle that needle core biopsies can serve as suitable input
   material for high-resolution downstream analyses.

   Having established the data quality from these different assays. We
   next performed integrative analyses across modalities in addition to
   demonstrating the feasibility of multi-omics technologies. These
   cross-cutting analyses provided insight into the metabolic states of
   specific immune and tumor cell types. It was, for example, noted that
   the linoleic acid pathway was upregulated in CD45+ regions of the TME.
   We were also able to cross-correlate spatial proteomics imaging and
   whole transcriptome regions of interest, highlighting the utility of
   these interrelated modalities. Immunopeptidomics assays also allowed us
   to correlate the presented HLA Class I peptidome at the bulk and
   single-cell RNA expression level, and in concordance with previous
   studies, show that gene expression correlated poorly with the detected
   surface peptides^[277]90,[278]91, again highlighting the need for
   orthogonal multi-omics measurement.

   More advanced integrative approaches were also performed to both infer
   and validate cell-cell interactions between predicted cell types using
   a combination of single-cell RNA sequencing data, phosphoproteomics and
   cyclic immunofluorescence. Cellchat is a popular cell-cell interaction
   algorithm that takes single-cell RNA sequencing data as input and
   assigns probability scores for specific ligand-receptor (L-R) pairs and
   cell-type pairings^[279]63. The algorithm results in a large number of
   potential cell-cell interactions and validation data for this approach
   is currently lacking—we therefore took advantage of our multi-omics
   approach and used phosphoproteomics data as “ground truth” to filter
   out tyrosine kinase receptors which were not phosphorylated, then
   ranked the strength of cell-cell interaction based on the Cellchat
   probability score. This integrative approach not only identified which
   ligand-receptor pair was highly expressed in scRNA and detected on mass
   spectrometry (SPP1:CD44) but also predicted cell types where this
   interaction would occur between. We cross-validated this prediction
   through CyCIF staining of both the L-R pair and predicted cell types.
   These results confirmed SPP1 secretion by glioma macrophages which has
   previously been reported as interacting with T-cells and tumor
   cells^[280]65–[281]67,[282]92,[283]93, and has been borne out by this
   integrative analysis. Future work will be focused on identifying
   further L-R pairs and interacting cell types using these methods.

   GBM tumor heterogeneity is a recognized phenomenon, which in turn
   impacts sampling variance. A perennial concern with any biopsy is how
   representative the sampled tissue is of the underlying tumor. Here we
   present a conceptual and analytic framework for quantifying
   heterogeneity. By first defining the unit of sampling as a
   microregion—an area that is accessible by rotating a stereotactic
   biopsy needle in a specific location, we then compare CNV and
   histological features intra- and inter-regionally. Additionally, we
   propose a simple rule to determine outlier cores—defined as cores that
   have a higher intra-regional variance than inter-regional variance. We
   proceed to show how this variance can be captured and quantified using
   H&E staining through a combination of image feature classification and
   Monte Carlo subsampling methods, and show that in this particular set
   of patient samples, intra-regional variance was lower than
   inter-regional samples.

   Having this wealth of data will also assist in efforts to develop
   peripheral biomarkers. Future efforts will include concurrent analysis
   of sample blood and CSF, with the goal of identifying peripheral
   biomarkers of therapeutic or immune response to treatment which could
   obviate the need for repeated tissue sampling. Future single-cell
   sequencing of samples would also include TCR VDJ sequencing, as doing
   so over multiple timepoints will enable the tracking and identification
   of responding T cells clones based on their frequencies in the
   repertoire. Moreover, from the single-cell expression data, the
   phenotypic state of those clones could also be assessed, validating
   their effector function. We would also be able to correlate the
   adaptive immune activity of T cells with clinical and immunopeptidome
   data, possibly linking T cells and targets. Identification of shared
   and overlapping TCR clones between patients could also be used to find
   broadly reactive clones. This approach would also be supplemented by
   TCR sequencing of blood samples that would serve as a patient-specific
   background, highlighting tumor-specific expanded clones.

   Finally, from a translational standpoint stemming from a need to create
   patient-specific tumor models integrated tightly with multimodal data,
   we show here that needle core biopsies often contain sufficient tissue
   to generate PDX models, with 42% of needle cores going on to generate
   viable tumors in immunocompromised mice. Somewhat counter-intuitively,
   the success rate in generating in vivo tumor avatars was higher than
   the success rate in generating patient-derived neurosphere 3D cell
   lines from the same tissue in parallel, which may be related to the
   conduciveness of the initial seeding environment. Another conclusion of
   our work is that the total number of viable cells able to be isolated
   from a core biopsy (integrative of tumor content and biopsy mass)
   seemed to be potentially the most important factor related to success.
   A caveat of our conclusions for this is that our sample size was low
   and so further clarity may come from an expanded cohort size over time.
   However, a clear message was that successful models were generated from
   extreme low end of each parameter and so no minimum limit of tumor
   content or viable cell number required for success is likely able to be
   set. Furthermore, based on the literature and multiple cell state
   diversity now documented in gliomas, success rate improvement in PDCLs
   from low input cell numbers may require personalized matching of media
   and growth factors to the unique tumor types.

   In conclusion, we have demonstrated that needle biopsy samples taken
   from a single procedure timepoint can be used to generate
   high-resolution cellular and molecular information regarding immune
   cell content, gene expression profiles, and activated protein pathways
   and antigenic diversity and provide an analytic framework to quantify
   and study core to core variance. Although multi-omics analysis was
   performed on a single patient, each individual analysis component was
   developed separately and optimized in a wider cohort of human specimens
   before integration into the final workflow. Data generated from a
   single patient therefore represents the cumulative end product of
   extensive logistical, protocol and platform development. Although data
   from a small number of patients limits the generalizability of
   conclusions, the ability to simultaneously obtain multimodal
   information from clinically accessible samples remains an important
   proof of principle, and as such, it represents a significant leap in
   capability in the interrogation of samples that were previously thought
   to be of limited value outside of standard diagnosis due to low tissue
   volume. This platform permits the introduction of investigative
   biopsies which have a clear translational motivation: To monitor and
   track treatment response in the setting of n-of-1 studies, and
   high-resolution sampling naturally lends itself to extending this
   technique to multiple timepoints. The lack of cellular and tissue level
   data in tumors undergoing treatment is considered a major roadblock to
   therapy development, and serial sampling is a critical step toward
   identifying tissue response phenotypes and tumor-stromal interactions
   that may occur during the course of disease and treatment.
   Additionally, investigative biopsies could also play a role in
   non-cancer diseases such as neuro-inflammatory conditions to monitor
   response to treatment. Future directions include further refinements in
   technology and dissociation techniques which could permit analysis of
   archival FFPE tissue, further widening the scope for retrospective
   analysis in historically treated samples, and integration of additional
   modalities, including data types and analysis tools; Finally, the
   establishment of a robust multi-omics analysis platform for needle core
   biopsy specimens forms the basis of a clinical trial involving
   longitudinal sampling in 12 patients which is currently underway
   ([284]NCT03152318).

Methods

   This study and its associated tissue collection and processing were
   governed under institutional ethical regulations. Institutional Review
   Boards: (MSKCC #06-107 & #09-156, DF/HCC #10-417 & #16-557). All
   procedures performed in studies involving human participants were in
   accordance with the ethical standards of the institutional and/or
   national research committee and with the 1964 Helsinki Declaration and
   its later amendments or comparable ethical standards. Written patient
   consent was obtained prior to surgery in all cases. Patients were not
   compensated for this study.

Experimental design

   This is a single-subject experimental design where multi-modality
   measurements are taken from tissue samples from the same subject. For
   the PDX generation section, the experimental design is that of a
   retrospective correlative study examining tissue from a cohort of
   patients who were enrolled in a prior clinical trial arm. Age, sex
   and/or gender analysis was not carried out due to small patient
   numbers.

Statistics and reproducibility

   This is primarily a feasibility study involving multiple n-of-1
   analyses, which inherently limits statistical analysis involving
   biological replicates. To ensure and enhance reproducibility, the
   methodology and protocols for each individual assay have been detailed,
   and technical replicates were performed where possible.

Sample collection, geotagging and transportation

   Biopsy cores were collected immediately during operation and
   transported in MACS storage solution (Miltenyi, Germany). Samples were
   then sent to the primary processing institution (MDACC) on ice via
   courier.

Sample distribution

   Formalin-fixed samples were processed locally at MSKCC (Core 1). Fresh
   tissue placed in tissue storage solution and was immediately shipped
   for central processing at MDA (Cores 2 and 3). Snap-frozen specimens
   (Cores 4–5) were shipped to BWH/DFCI and MIT, respectively.

Low-volume glioma tissue processing

   At the receiving institution, tissue cores were processed within 24 h
   of the operation. Fresh tissue samples were directly placed in MACS
   tissue storage buffer (Miltenyi Biotec, Germany) on ice before
   transportation to the processing laboratory. Given the potential need
   to transport tissue from different institutions, we also tested the
   performance of tissue storage buffer (results in Supplementary
   Table [285]4). Once received in the laboratory, tissue cores were
   morselized into small pieces under sterile conditions using sterile
   scalpels, triturated using pipettes, then transferred to a 50 ml falcon
   tube (Corning, USA), then incubated with papain (37 °C, 25 min)
   (Worthington, USA) according to manufacturer instructions. After
   incubation, HBSS without magnesium or calcium was then added to make up
   to 50 ml, and the sample was centrifuged (300×g, 5 min, 4 °C),
   supernatant aspirated before being passed through a 70-μm cell strainer
   (BD, USA). Erythrocytes were removed by resuspending and incubating the
   obtained pellet in 5 ml of ACK-lysis buffer (eBioscience, USA) for
   10 min, followed by centrifugation (300×g, 5 min, 4 °C). Residual
   myelin and extracellular debris were eliminated using the Debris
   Removal Kit (Miltenyi Biotech, Germany). Cell counts were quantified
   using an automated cell counter (Nexcelom, USA) after resuspending the
   pellet in PBS. Cells were stored in Bambanker (Wako, Japan). Cell
   suspensions were immediately placed in a storage container (Mr.
   FrostyTM, Thermo Fisher Scientific, USA) and stored in LN2.

Liquid chromatography-tandem mass spectrometry (LC-MS/MS)

   For pTyr proteome, the enrichment of pTyr peptides was performed by
   incubating digested peptides with 60 μl of protein G agarose beads
   (Sigma, USA) pre-conjugated to 24 μg of 4G10 (BioXcell, USA) and 6 µl
   of PT66 antibodies overnight at 4 °C in immunoprecipitation (IP) buffer
   (100 mM Tris-HCl, 1% NP-40, pH 7.4). Subsequently, immunoprecipitated
   peptides were eluted twice with 0.2% trifluoroacetic acid (TFA) for
   10 min each and enriched for phosphopeptides using ferric
   nitrilotriacetate (Fe-NTA) columns (Thermo Scientific, USA). Eluates
   from the Fe-NTA columns were dried using vacuum centrifuge,
   reconstituted in 4 µl of 3% acetonitrile/0.1% formic acid, and directly
   bomb-loaded onto a hand-packed 10 cm analytical column containing 3 μm
   C18 beads. Agilent 1100 Series HPLC (Agilent Technologies USA)
   connected to Orbitrap Exploris 480 Mass Spectrometer (Thermo
   Scientific, USA) was operated at 0.2 ml/min flow rates with a precolumn
   split to attain nanoliter flow rates through the analytical column and
   nano-electrospray ionization tip. Peptides were eluted with the
   increasing concentrations of buffer B using the following 140-min
   gradient settings (Buffer A: 0.1% acetic acid; Buffer B: 70%
   acetonitrile/0.1% acetic acid): 0 min: 0% B; 10 min: 11% B; 110 min:
   32% B; 125 min: 60% B; 130 min: 100% B; 128 min: 100% B; 130 min: 0% B.
   The MS parameters were: ESI spray voltage, 2.5 kV; no sheath or
   auxiliary gas flow; heated capillary temperature, 275 °C,
   data-dependent acquisition mode with the scan range of 380–2000 m/z.
   MS1 scans were acquired at 60,000 resolution, maximum injection time of
   100 ms, normalized AGC target of 300%, and only included the precursor
   charge states of ≥2 and ≤6. For every full scan, MS/MS spectra were
   collected during a 3-s cycle time. For MS/MS, Ions were isolated
   (0.7 m/z isolation width) for the maximum IT of 250 ms, normalized AGC
   target of 100%, HCD collision energy of 30% at a resolution of 60,000,
   and the dynamic exclusion time of 35 s. Half of the supernatants from
   pTyr-IP were subjected to high-pH reverse-phase fractionation on a
   Kromasil® C18 HPLC Column (5 μm particle size, pore size 100 Å, L ×
   I.D. 250 mm × 4.6 mm) (Nouryon, Sweden) using buffer A (10 mM
   triethylammonium bicarbonate (TEAB), pH 8) and buffer B (10 mM TEAB, pH
   8, 99% acetonitrile) over an 85-min gradient (0 min: 1% B, 1 min: 1% B,
   5 min: 5% B, 65 min: 40% B, 75 min: 70% B, 84 min: 70% B, 85 min: 1%)
   into 10 fractions (concatenated) using a Gilson FC 204 Fraction
   Collector (Marshall Scientific, USA). Fractionation was performed at a
   flow rate of 1 ml/min and 1 min per fraction between 10 and 85 min
   portion of the gradient. For each fraction, 1/10 was allocated for
   protein expression analysis, and 9/10 was allocated for global
   phosphoproteomic (pSer/Thr) analysis. Fractions for pSer/Thr analysis
   were resuspended in 50 µl of 0.2% TFA and subjected to Fe-NTA
   column-based phosphopeptides enrichment, after which they were eluted,
   dried, and resuspended in 3% acetonitrile/0.1% formic acid for
   injection by Dionex UltiMate 3000 Autosampler (Thermo Scientific, USA).

   For immunopeptidome, MHC-peptide complexes were isolated by IP
   (anti-HLA-A/B/C; clone W6/32 (BioXcell, USA); conjugated to 20 μl
   FastFlow protein A sepharose bead (GE Healthcare, USA)), eluted with
   10% acetic acid, filtered by size exclusion filter (Nanosep 10 K,
   PALL), cleaned up with zip tip (Thermo Scientific, USA), and dried with
   speed-vac. Samples were resuspended in 3% acetonitrile/0.1% formic acid
   and directly loaded onto an analytical capillary chromatography column
   (about 1 μm orifice) prepared and packed in-house (50 μm ID × 15 cm and
   1.9 μm C18 beads, ReproSil-Pur). The sample was analyzed using an
   Orbitrap Exploris 480 Mass Spectrometer (Thermo Scientific, USA),
   coupled with an UltiMate 3000 RSLC Nano LC system (Dionex, USA),
   Nanospray Flex ion source (Thermo Scientific, USA) and column oven
   heater (Sonation; operated at 45 °C). Peptides were eluted using a
   125-min gradient (Buffer A: 0.1% formic acid, Buffer B: 80%
   acetonitrile/0.1% formic acid) with the following minute:B% profile:
   0:3, 30:3, 31:8, 85:30, 97:50, 100:97, 103:97, 103.1:3. The flow rate
   of 0.4 μl/min was used between 0 and 29.5 min and 0.075 μl/min between
   30 and 125 min. The MS parameters were: spray voltage, 2.5 kV; no
   sheath or auxiliary gas flow; heated capillary temperature, 275 °C,
   data-dependent acquisition mode with the scan range of 350–1200 m/z.
   MS1 scans were acquired at 60,000 resolution, maximum injection time of
   50 ms, standard AGC target, and only included the precursor charge
   states of ≥2 and ≤4. For every full scan, MS/MS spectra were collected
   during a 3-s cycle time. For MS/MS, Ions were isolated (0.7 m/z
   isolation width) for the maximum IT of 250 ms, normalized AGC target of
   100%, HCD collision energy of 30% at a resolution of 60,000, and the
   dynamic exclusion time of 30 s.

Mass spectrometry data analysis

   Mass spectra were analyzed using Proteome Discoverer (v3.0, Thermo
   Fisher Scientific, USA) and searched using Mascot (v2.8) against the
   human Swiss-Prot database (v2022_11). For peptide-bound HLA, peptides
   were searched with no enzyme and variable methionine oxidation. Peptide
   spectrum matches were filtered by an ion score ≥15, length 8–11, search
   engine rank of 1, and aggregated across unique peptides. GibbsCluster
   2.0 was used for motif analysis^[286]94, and the logos were generated
   using ggseqlogo R package. For phosphoproteomic data, raw files were
   searched with two or fewer missed cleavages, precursor and fragment ion
   matched with a 10 ppm and 20mmu mass tolerances, and fixed
   modifications of carbamidomethyl (Cys) and dynamic modifications of
   phosphorylation (Ser, Thr, Tyr) and oxidation (Met). After
   quantification by minora feature detection method, peptide spectrum
   matches were filtered by an ion score ≥20 and search engine rank of 1.
   The ptmRS node^[287]95 in Proteome Discoverer was used to confidently
   assign the phosphorylated and oxidized sites (filtered for peptides
   with localization probability ≥90% after search). All data were
   processed in R studio (v4.1.0).

Single-cell RNA sequencing

   scRNAseq was performed using the 10x Genomics Chromium Single Cell
   Controller. Single-cell suspensions were prepared from tissue biopsy as
   described above. Freshly isolated cells were droplet-separated using
   the Chromium Single Cell 3′ v.3 Reagent Kit (10x Genomics, USA) with
   the 10x Genomics microfluidic system creating cDNA library with
   individual barcodes for individual cells. Barcoded cDNA transcripts
   from patients were pooled and sequenced using the NovaSeq 6000
   Sequencing System (Illumina, USA).

Single-cell RNA sequencing analysis

   Raw single-cell RNA seq reads generated by Illumina sequencer were
   demultiplexed into FASTQ and aligned to GRCh38 reference genome to
   generate count matrices using Cell Ranger v6.0.0 analysis pipelines
   (10x Genomics). The preliminarily processed data was used for
   downstream analysis by Seurat v4.3.0^[288]96. Low-quality and dead
   cells that had <200 or >2500 unique transcripts or >5% of mitochondrial
   transcripts were filtered out. Multiple datasets (with and without
   inhibitors) were combined, normalized, identified variable features,
   and applied principal component analysis (PCA) using Seurat. Batch
   effect correction of combined data was performed using Harmony
   algorithm v0.1.1^[289]97. Following batch correction, FindClusters was
   applied with a resolution of 0.5 to identify cell clusters and uniform
   manifold approximation and projection (UMAP). For visualization, plots
   were generated using a combination of Seurat, ggplot2, and scCustomize
   v2.1.2^[290]98.

Spatial proteogenomics

   Digital Spatial profiling (DSP, Nanostring) experiments were performed
   according to the manufacturer’s protocol. Briefly, slides were stained
   and imaged to visualize morphology markers, GFAP-AF488, 53-9892-82,
   Clone GA5, diluted 1–100, Syto83, S11364, diluted 1:25, CD45-AF594,
   NBP1-44763AF594, clone EM05, diluted 1:100, Ki67-AF647,
   NBP2-22112AF647, clone 8D5, diluted 1:100. Images at ×20 magnification
   were assembled to yield high-resolution regions of interest. As
   previously described, ROIs were illuminated, and released tags were
   collected into 96-well plates. Sequential sections were used for GeoMx
   Protein and RNA profiling. For protein detection GeoMx (v1.0) Human NGS
   Protein Core, (v1.0) Human NGS Cell Death Protein, (v1.0) Human NGS
   Glial Cell Subtyping Protein, (v1.0) Human NGS Immune Activation Status
   Protein, (v1.0) Human NGS Immune Cell Typing Protein, (v1.0) Human NGS
   IO Drug Target Protein, (v1.1) Human NGS MAPK Signaling Protein, (v1.0)
   Human NGS Myeloid, (v1.0) Human NGS Neural Cell Typing Protein, (v1.0)
   Human NGS Pan-Tumor Protein, (v1.0) Human NGS PI3K/AKT Signaling
   modules were used to identify cell types. Detection of mRNA was
   performed for morphology marker–defined compartments according to the
   NanoString GeoMX RNA assay protocol described previously^[291]99 using
   the Whole Transcriptome Atlas probe reagent. Library preparation was
   also performed according to the manufacturer’s protocol. DSP collection
   sample plates were dried, resuspended in nuclease-free water, and
   amplified using PCR according to the manufacturer’s protocol. Purified
   libraries were sequenced by Illumina NovaSeq6000. The FASTQ reads from
   the sequenced DSP library were processed by the GeoMx NGS DnD Pipeline
   to convert sequencing reads into counts per ROI (NanoString,
   MAN-10133_03). After processing, counts were uploaded to the GeoMx DSP
   Data Analysis Suite (NanoString Technologies, USA).

CODEX

   A 4 µm FFPE tissue section was placed on a positively charged glass
   slide and stored at 4 °C until use. Codex staining assay was performed
   following the manufacturer’s protocol. Purified antibodies, barcodes
   and reporters were purchased commercially, as listed in Supplementary
   Table [292]4. Antibodies CD1c, Bcl-6, DC-LAMP, CXCL13, CXCR5, and CD206
   were conjugated through the Spatial Tissue Exploration Program (STEP)
   from Akoya Biosciences. For GFAP, Olig2, FoxP3, CCR7, Lag3, CD138, and
   AIF1, purified antibodies (free of BSA and glycerol) were conjugated
   in-house. The remaining antibodies were commercially tagged with
   PhenoCycler® Barcodes at the time of purchase from Akoya Biosciences.
   The in-house antibody conjugation process followed the manufacturer’s
   guidelines. Purified barcoded-conjugated antibodies were used within 6
   months of conjugation.

Codex imaging, processing, segmentation, and analysis

   The stained tissue section was captured using PhenoCycler-Fusion
   (version 1.0.8) using DAPI, ATTO550, CY5 and AF750 filters at a scan
   resolution of 0.50 µm (20×) and saturation protection was applied on
   DAPI channel. Alignment of images across cycles, stitching of tiles and
   subtraction of auto-fluorescence was performed using CODEX® Processor
   application. Cell segmentation was performed by applying StarDist, a
   deep learning-based algorithm, to the DAPI channel to segment nuclei;
   cells were approximated based on expansion from the nuclei. Scripts for
   performing the segmentation in the QuPath software and a pre-trained
   model were provided by Akoya Biosciences. All mIF analyses were run in
   R-4.2.1.

Cell clustering

   The cell expression values were first Z-score normalized across
   markers, followed by the application of scaled and centered biomarker
   expression value for each cell (Seurat citations). Principal component
   analysis was performed to identify 20 principal components
   (“Seurat::RunPCA, npcs = 20”), which were then used in the Harmony
   algorithm to group cells and correct for dataset-specific and
   experimental conditions (“harmony::RunHarmony”)^[293]96,[294]100. The
   Harmony reduction was used for downstream analysis. The dimensionality
   was reduced with UMAP (“Seurat::RunUMAP, reduction = ‘harmony’”).
   Clustering was performed using a shared nearest neighbor modularity
   optimization-based algorithm. The UMAP plot and a heatmap plot showing
   the average normalized marker expression in each cluster were graphed
   using the Seurat package. Clusters were annotated using their average
   expression to identify cell types, and these annotations were validated
   by manual inspection of multiplexed IF stains on images. The cluster
   spatial location plot was graphed using the ggplot2 package.

Cell neighborhood analysis

   To define cellular neighborhoods (CNs), the nearest 20 neighbors for
   each cell were recorded and set as a window. The number of each cell
   cluster in the window for each cell was counted, resulting in a matrix
   of cells by cell clusters, with each row representing a cell, each
   column representing a cell annotation (cell type) from the clustering
   above, and each value representing the count of neighbors of the given
   annotation. The neighbor cell proportion was computed for each row. The
   resulting matrix was clustered using k-means clustering
   (“stats::kmeans”), where the optimal k was determined empirically by
   maximizing the silhouette score metric (“factoextra:: fviz_nbclust”).
   Each cluster was defined as a CN. Thus, each cell was given both a cell
   type annotation, which depends only on the cell’s own marker
   expression, and a cell neighborhood annotation, which depends on the
   cell’s type and the identities of its nearest neighbors.

MALDI-MSI tissue preparation and imaging

   Biopsy Core Four underwent a preservation process, being rapidly frozen
   on dry ice. For MALDI MSI, the core was then cryo-sectioned into 10
   µm-thick slices. These slices were subsequently thaw-mounted onto an
   indium tin oxide (ITO) slide. Serial sections were prepared for
   microscopy staining. For the mass spectrometry analysis, we utilized a
   1,5-diaminonaphthalene hydrochloride MALDI matrix, with 15 N glutamate
   spiked in as an internal standard. The matrix was prepared at a
   concentration of 4.3 mg/ml in a mixture of 4.5 parts HPLC grade water,
   5 parts ethanol, and 0.5 parts 1 M HCl (v/v/v). A TM-sprayer from HTX
   Imaging was used to spray the matrix, with the following parameters: a
   flow rate of 0.09 ml/min, spray nozzle velocity of 1200 mm/min, spray
   nozzle temperature of 75 °C, nitrogen gas pressure of 10 psi, track
   spacing of 2 mm, and a four-pass spray. The mass spectrometry analysis
   was carried out using a 15 Tesla SolariX XR FT-ICR MS (Bruker
   Daltonics, Billerica, MA). The instrument was set to negative ion mode,
   and the mass range scanned was from m/z 46.07 to 3000, with a sampling
   step size of 30 µm. Each sampling point consisted of 200 laser shots at
   a laser power of 21% (arbitrary scale), with a laser repetition rate of
   1000 Hz. We employed the Continuous Accumulation of Selected Ions
   (CASI) mode, setting Q1 to m/z 150 with an isolation window of 200. The
   mass range was calibrated using a tune mix solution from Agilent
   Technologies with the electrospray source. Additionally, the internal
   standard 15 N glutamate was used for online calibration during the
   acquisition. For data analysis, SCiLS Lab software (version 2023c Pro,
   Bruker Daltonics, Billerica, MA) was used to view and process ion
   images and mass spectra. The dataset was normalized to the total ion
   current (TIC), and peaks were annotated using Metaboscape (2021b,
   Bruker Daltonics, Billerica, MA). Metabolites were putatively annotated
   based on an accurate mass with a Δppm <0.5 and MSMS measurements.
   MetaboanalystR^[295]101 was used to perform pathway enrichment analysis
   at every pixel. The resulting enrichment ratios were displayed
   spatially with an in-house R script. Metabolite differential expression
   was tested with a t-test and FDR correction.

Tissue-cyclic immunofluorescence (t-CyCIF)

   Tissue-based cyclic immunofluorescence (t-CyCIF)^[296]1 of fresh frozen
   tissue was adapted from a detailed protocol available at protocols.io
   (dx.doi.org/10.17504/protocols.io.bjiukkew). Fresh frozen tissue
   samples mounted on SuperFrost slides (Fisher Scientific #12-550-15,
   Hampton NH) were fixed with 4% PFA for 1 h at room temperature,
   followed by membrane permeabilization for 5 min in 0.5% Triton X in PBS
   (Thermo Scientific #28316). All primary and conjugated antibodies were
   incubated overnight at 4 °C in the dark. Secondary antibodies targeted
   to unconjugated primary antibodies in the first cycle were incubated
   for 1 h at room temperature. Hoechst 33342 co-staining was performed by
   combining Hoechst with primary antibodies in SuperBlock Blocking buffer
   (Thermo Scientific #37515) buffer at every cycle. Coverslips were
   mounted on slides with 10% Glycerol in PBS prior to imaging. Following
   imaging, coverslips were removed in 1X PBS, then fluorophores were
   inactivated by incubating slides in a solution of 4.5% H[2]O[2], 24 mM
   NaOH in PBS for 1 h under an LED light source. Slides were processed
   through multiple cycles of antibody incubation, imaging, and
   fluorophore inactivation. CyCIF images were acquired on a CyteFinder II
   slide scanning fluorescence microscope (RareCyte Inc. Seattle WA) with
   CyteFinder software (v3.11.024). t-CyCIF image processing and data
   quantification were performed as previously described^[297]2. Stitching
   and registration of tiles and cycles were done in MCMICRO^[298]3 using
   ASHLAR (v1.10.2)^[299]4 module. Code is available on GitHub
   ([300]https://github.com/labsyspharm/cd73_coy_spatialcorrelation) and
   Zenodo ([301]https://zenodo.org/record/6628875#.YqJVfqjMJD8).
   10.5281/zenodo.6628875. Additional details and code can be found at
   [302]www.cycif.org. A list of antibodies used in CyCIF is available in
   Supplementary Table [303]5.

MALDI-MSI and t-CyCIF integration

   t-CyCIF images were registered onto the MSI space using an in-house
   MATLAB script. Several fiducials were specified manually around the
   contour and specific tissue landmarks. An affine transformation was
   obtained from the points using the MATLAB function fitgeotrans. The
   obtained transformation was applied to all CyCIF channels.

MALDI-MSI and scRNAseq integration

   Metabolite pathway enrichment analysis was performed using
   MetaboAnalyst^[304]102 for spatial metabolomics data and
   PANTHER^[305]103 for scRNAseq data. The enrichment based on spatial
   metabolomics data was used to select a list of top 5 pathways enriched
   in each group (immune cells vs non-immune cells). Due to mismatching
   names in the corresponding libraries KEGG^[306]104 and GO^[307]105 the
   crosslink was performed manually.

PDCL and PDX generation

   Patient samples for model generation DF/HCC protocol #10-417 from
   patients who planned to enroll in DFCI clinical trial #16-557
   [308]NCT03152318 of oncolytic HSV agent CAN-3110. Biopsy cores deemed
   to be in excess of clinically required tissue were allocated for
   research use. Fresh tumor needle biopsy tissue specimens were
   transported in DMEM + 1% penicillin-streptomycin. Fresh tumor biopsy
   tissue was dissociated in NeuroCult NS-A media to a single-cell level
   and viable cells were counted using a TC20 cell counter (Bio-Rad). If
   the total cell count was greater than 100,000 cells, an intracranial
   mouse injection for PDX was performed. If the total cell count was
   greater than 3,000,000 cells, then the intracranial mouse injection for
   PDX was accompanied by a parallel 3D spheroid cell line attempt. Mean
   needle biopsy mass of the samples was 37.8 mg, and mean cell count
   post-processing was 0.76 million cells.

   PDX attempts were conducted via intracranial injection per DFCI
   institutional animal protocol and protections. Nude 6-week-old mice of
   both sexes were sedated using a 0.2 L/min isoflurane administration
   along with oxygen and arranged within a stereotactic apparatus. Tumor
   cells were injected into the right striatum as 100,000 viable cells in
   1 µL of NSA media at coordinates from bregma: D-V −2.5 mm, M-L 2 mm,
   A-P 0 mm. Mice were subsequently monitored at a minimum of 3 times per
   week for symptomatology as the endpoint. Histological confirmation of
   the presence of human tumor was performed by a neuropathologist on H&E
   analysis at the first harvest and passage of all models (KLL). Models
   were evaluated for histological evidence of patient donor immune cell
   infiltration or patient fibroblast or other TME cell types, and no
   models showed evidence of lymphocyte preservation or non-tumor cell
   types within the mouse brain. The PDX tumors showed variable growth
   patterns with highly infiltrative growth and interaction with the TME
   as well as other patient models having circumscribed growth of the
   tumor and limited TME interactions at the edge of the model as
   previously described^[309]76. Additional methods were as previously
   published^[310]75,[311]106.

   PDCL attempts were performed by placing a minimum of 100,000 viable
   cells on laminin-coated T25 flasks and cultured in the NeuroCult NS-A
   Proliferation Kit media (StemCell Technologies) supplemented with
   0.0002% heparin (StemCell Technologies), EGF (20 ng/ml), and FGF
   (10 ng/ml; Miltenyi Biotec) in a humidified environment of 5% CO[2] at
   37 °C^[312]75,[313]107. Models were growth-verified once five serial
   passages were completed indicating long-term growth. Growth-verified
   models were genomically verified via whole exome sequencing for glioma
   diagnostic copy and mutation aberrations.

   Growth-verified models were genomically verified via whole exome
   sequencing for glioma diagnostic copy and mutation aberrations. All
   models and relevant data are available from the DFCI Center for Patient
   Derived Models (models@dfci.harvard.edu). Tumor genomic alteration data
   was obtained by targeted next-generation sequencing via DFCI Oncopanel
   or MSK-IMPACT. These validated measures identify genomic alterations,
   including single-nucleotide variants, insertions and deletions, copy
   number alterations, and structural variants in the exons and selected
   introns of 447 (Oncopanel) and 410 (MSK-IMPACT) genes, including genes
   relevant to the diagnosis and subclassification of HGG. Sequencing
   results were reviewed and reported by a board-certified
   neuropathologist. Sequencing data from these assays, as well as
   clinical information regarding age, trial cohort, diagnosis at the time
   of injection, and MGMT promoter methylation status, was used to
   generate an Oncoprint genomic profile using R 4.2.2, RStudio
   2022.12.0+353, and the Oncoprint function of the ComplexHeatmap 2.14.0
   package. Statistical analysis (Student’s t-test) evaluating the
   influence of tumor percentage on model success was conducted using the
   same version of R (4.2.2) and RStudio (2022.12.0 + 353).

Bulk DNA and RNA sequencing and analysis

   Each core was processed for DNA, RNA, and protein extraction using the
   Nucleospin Tri-kit (Machery-Nagel). Samples were homogenized first in a
   1 ml pipette and then 10 times through a 26 G needle before being
   prepared per the manufacturer’s recommendations. Extracted DNA and RNA
   were evaluated by Fragment Analyzer (Agilent), and 1 ng of DNA was used
   for nexteraXT library preparation (Illumina) before low pass CNA whole
   genome sequencing on either an Illumina NextSeq500 (40nt paired end) or
   Singular G4 (50nt paired end). Bulk RNAseq libraries were prepared from
   10 ng of RNA using the Pico-input strand-specific total RNA-seq for
   mammalian samples v2 kit from Takara. DNA samples were sequenced on an
   Illumina NextSeq500 (40 paired end). RNA samples were sequenced on a
   MiSeq to provide validation of the methodology using a 60nt forward
   read and were analyzed using the RNAseq NextFLow pipeline. DNA reads
   were mapped to the hg19 reference genome using BWA-MEM. We calculated
   the number of reads in 500 kb genomic bins throughout the genome. In
   the case of tumor samples, read counts were normalized based on the
   read counts from the PBMC sample, followed by GC correction using
   HMMcopy^[314]108 to get read depth ratio’s R. These values were then
   logged to get LogR values plotted in the figure. Gains and losses were
   identified using QDNAseq^[315]109 with a probability cutoff of 0.5. To
   estimate tumor fraction we calculated the purity that resulted in
   average purity corrected copy number values of 3 for chromosome 7 and 1
   for chromosome 10. Gains on chromosome 7 and losses on chromosome 10
   are ubiquitous features in GBM and were assumed clonal. Phylogenetic
   tree were constructed by generating a binary event matrix where rows
   were genome segments and columns were samples. Genome segments were
   defined as continuous regions of the genome that had gains or losses in
   any sample, we ensured that segment boundaries were identical across
   samples. Matrix values were 1 if a sample had a copy number alteration
   in that segment and 0 otherwise. Hamming distances were calculated
   between all samples and a tree was constructed using the neighbor
   joining algorithm using the distance matrix as input.

Multi-regional sampling variance analysis

   The annotations from an expert neuropathologist (TB) were used to train
   an H&E segmentation model using Ilastik^[316]110. The labels included
   technical artifacts (i.e., background, folds, bubbles) and
   histopathologically relevant features (i.e., viable tumor, infiltrative
   tumor, necrosis, vasculature…). All technical artifacts were removed
   from further analysis. An in-house R script was used to subsample a
   thousand 50 × 50 μm regions of interest (ROI) within each image in a
   Monte Carlo setting. For each ROI, the percentage area occupied by each
   label was quantified, deriving a distribution of sampled imaging
   features from each image. PCA was used to visualize the resulting
   distributions, and the lme function from the nlme R package was used to
   perform group comparisons.

Reporting summary

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

Supplementary information

   [318]Supplementary Information^ (10.5MB, pdf)
   [319]Transparent Peer Review file^ (390.3KB, pdf)
   [320]Reporting Summary^ (4.2MB, pdf)

Acknowledgements