Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions

Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions View ORCID Profile Saad M. Khan , View ORCID Profile Anthony Z. Wang , View ORCID Profile Rupen R. Desai , View ORCID Profile Colin R. McCornack , Rui Sun , View ORCID Profile Sonika M. Dahiya , View ORCID Profile Jennifer A. Foltz , Ngima D. Sherpa , Lydia Leavitt , View ORCID Profile Timothy West , Alexander F. Wang , Aleksandar Krbanjevic , Bryan D. Choi , Eric C. Leuthardt , View ORCID Profile Bhuvic Patel , Al Charest , View ORCID Profile Albert H. Kim , Gavin P. Dunn , View ORCID Profile Allegra A. Petti doi: https://doi.org/10.1101/2025.04.04.647096 Saad M. Khan 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Saad M. Khan Anthony Z. Wang 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 3 Department of Pathology and Immunology, Washington University in St. Louis School of Medicine , St. Louis, Missouri Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Anthony Z. Wang Rupen R. Desai 4 Department of Neurosurgery, Medical College of Wisconsin , Milwaukee, Wisconsin Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Rupen R. Desai Colin R. McCornack 5 Department of Neurological Surgery, Washington University School of Medicine , St. Louis, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Colin R. McCornack Rui Sun 6 Department of Immunology, School of Basic Medical Sciences, Anhui Medical University , Hefei, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sonika M. Dahiya 7 Division of Neuropathology, Department of Pathology and Immunology, Washington University School of Medicine , St. Louis, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sonika M. Dahiya Jennifer A. Foltz 8 Department of Medicine, Division of Oncology, Washington University School of Medicine , St. Louis, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jennifer A. Foltz Ngima D. Sherpa 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 9 Biological and Biomedical Sciences Graduate Program, Harvard University , Cambridge, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site Lydia Leavitt 10 Department of Neurological Surgery, University of Louisville , Louisville, KY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Timothy West 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Timothy West Alexander F. Wang 11 The University Program, Case Western Reserve University School of Medicine , Cleveland, OH, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Aleksandar Krbanjevic 12 Department of Pathology, Microbiology, and Immunology, Vanderbilt University , Nashville, TN, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Bryan D. Choi 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site Eric C. Leuthardt 5 Department of Neurological Surgery, Washington University School of Medicine , St. Louis, MO, USA 13 The Brain Tumor Center, Siteman Cancer Center, Washington University School of Medicine , St. Louis, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Bhuvic Patel 5 Department of Neurological Surgery, Washington University School of Medicine , St. Louis, MO, USA 13 The Brain Tumor Center, Siteman Cancer Center, Washington University School of Medicine , St. Louis, MO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Bhuvic Patel Al Charest 14 Department of Medicine, Beth Israel Deaconess Medical Center, Harvard Medical School , Boston, MA, USA 15 Cancer Research Institute, Beth Israel Deaconess Medical Center , Boston, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Albert H. Kim 5 Department of Neurological Surgery, Washington University School of Medicine , St. Louis, MO, USA 12 Department of Pathology, Microbiology, and Immunology, Vanderbilt University , Nashville, TN, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Albert H. Kim Gavin P. Dunn 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: apetti{at}mgh.harvard.edu Allegra A. Petti 1 Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts 2 Brain Tumor Immunology and Immunotherapy Program, Department of Neurosurgery, Massachusetts General Hospital, Harvard Medical School , Boston, Massachusetts Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Allegra A. Petti For correspondence: apetti{at}mgh.harvard.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Treatment resistance in glioblastoma (GBM) is largely driven by the extensive multi-level heterogeneity that typifies this disease. Despite significant progress toward elucidating GBM’s genomic and transcriptional heterogeneity, a critical knowledge gap remains in defining this heterogeneity at the spatial level. To address this, we employed spatial transcriptomics to map the architecture of the GBM ecosystem. This revealed tumor cell states that are jointly defined by gene expression and spatial localization, and multicellular niches whose composition varies along the tumor core-edge axis. Ligand-receptor interaction analysis uncovered a complex network of intercellular communication, including niche- and region-specific interactions. Finally, we found that CD8 positive GZMK positive T cells colocalize with LYVE1 positive CD163 positive myeloid cells in vascular regions, suggesting a potential mechanism for immune evasion. These findings provide novel insights into the GBM tumor microenvironment, highlighting previously unrecognized patterns of spatial organization and intercellular interactions, and novel therapeutic avenues to disrupt tumor-promoting interactions and overcome immune resistance. INTRODUCTION Glioblastoma (GBM), the most common primary malignant brain tumor in adults, continues to be a devastating disease, with a median survival of 15-20 months 1 – 3 . Of the many factor impeding the development of new therapies, a major challenge is the multi-level heterogeneity of the GBM ecosystem, which has been extensively studied from the perspective of genomics 4 – 6 , transcriptomics 7 , 8 , and single-cell RNA-sequencing (scRNA-seq) 9 – 15 . In particular, exome and genome sequencing studies have catalogued extensive intratumoral genetic heterogeneity in GBM. More recently, scRNA-seq studies have elucidated complex intratumoral transcriptional heterogeneity, with the emerging view that GBMs are mixtures of cells in different transcriptionally-defined cell states, although there is no unified classification of cell states for this disease. These studies have also provided insight into the GBM immune microenvironment, which is unique in two key ways that hinder the development of effective immunotherapies. First, it is “immune-cold” in that it contains few infiltrating cytotoxic T cells, and those T cells are dysfunctional, with only a small subset exhibiting a canonical exhaustion profile 16 – 20 . Second, it is enriched for immunosuppressive myeloid cells 21 that influence other immune compartments, for example, by driving T cell exhaustion 22 . Spatial heterogeneity in GBM is particularly multifaceted, ranging from microscopic histopathological hallmarks (e.g. regions of microvascular proliferation and necrotic foci surrounded by palisading cells) to larger spatial domains visible using MRI (necrotic core, viable tumor core, infiltrative margin). Recently, spatial transcriptomic (ST) technologies have transformed our ability to understand the spatial organization and in situ biological programs of tumor ecosystems. In GBM, ST has been used to show that regions of hypoxia in the tumor core drive higher-level tissue organization 23 , genomic instability 24 , and macrophage cell states 25 that promote vascular hyperpermeability 26 and influence prognosis 27 . It has also been shown that the tumor edge and core differ in their distribution of predefined tumor cell states and gene expression profiles 23 , 28 , with infiltrating tumor cells exhibiting Notch signaling and synaptic gene expression, for example. Such differences may have profound clinical implications, given that most recurrences arise from the tumor margin 29 . Despite recent progress in describing the spatial heterogeneity of diverse cell types within GBM, the underlying intercellular interactions and their spatial organization remain largely unexplored. Here, we leveraged ST, scRNA-seq, and multi-region tissue sampling to comprehensively define the spatial architecture of GBM. We used the Visium V1 platform to survey the spatial landscape of 22 samples along the edge-core axis of nine patients with primary isocitrate dehydrogenase (IDH) wildtype (WT) GBM, and the Xenium platform to analyze nine samples from six patients with primary IDH WT GBM. We integrated these data with scRNA-seq data from 13 GBM samples and a nonmalignant brain atlas 30 to derive a novel and unique set of tumor cell states defined jointly by transcriptional profile and spatial localization. We then identified spatially segregated multicellular communities, or niches, which vary along the edge-core axis and display distinct tumor cell state compositions. To identify potential therapeutic targets, we identified ligand-receptor pairs that sustain intercellular communication within and between these niches as well as larger spatial domains. Using the clinically relevant interaction between hypoxia and vasculature as a case study, we identified niche-specific ligand-receptor interactions promoting angiogenesis, highlighting the complex, and highly redundant, network of interactions that drive this critical histopathological pathway in GBM. Finally, we harnessed the single cell resolution of the Xenium platform and found that the majority of CD8+GZMK+ T cells, the most predominant CD8+ T cell population within GBM TIL 20 , tend to colocalize with LYVE1+CD163+ myeloid cells in vascular regions, raising the possibility that these heterotypic interactions may attenuate T cell presence and function in GBM. Together, these results provide a spatially resolved portrait of GBM across multiple resolutions, from single cells to broad tissue domains, and identify potentially targetable ligand-receptor interactions underlying this organization. RESULTS Identification of tumor core, edge, and core-to-edge transition zones GBM differs from many solid tumors in that it is highly infiltrative and lacks clear margins, complicating both surgical resection and data analysis. To reliably identify and collect tumor specimens from locations along the core-edge axis, we employed MRI and/or intraoperatively visualized 5-aminolevulinic acid fluorescence gradients from core to edge 31 (Table S1). Based on such intraoperative localization, twenty-two samples were obtained from nine patients undergoing standard-of-care surgery for primary GBM, and designated (1) tumor core ( n =7), (2) tumor edge ( n =6), or (3) core-to-edge transition zone ( n =9). Each sample was preserved in OCT or FFPE, stained with H&E, and analyzed using the Visium V1 spatial transcriptomic platform ( Fig. 1A , Fig. S1, Table S1). Download figure Open in new tab Figure 1: Data Overview. (A) Summary of samples, patients, data types, and analyses. (B) Clockwise from top left: H&E image for core tumor sample from B189, neuropathologist’s tumor content estimates, neuropathologist’s identification of microvascular proliferation (MVP), and unsupervised graph-based clusters depicting MVP and additional variation. (C) Louvain gene expression clusters in representative core samples, with cluster IDs and Gene Ontology enrichments as indicated. (D) Clusters in representative edge samples. (E) Clusters in representative transitional 5ALA+ samples. (F) Hypoxia-response cluster, C3, in core samples. (G) MVP cluster, C6, in core samples and transitional 5ALA+ samples. (H) Glial (C0, C1) and neuronal (C2) clusters in edge samples. Several methods were first employed to distinguish spatially distinct domains within our GBM ST data. Histopathologic tumor core and edge designations by a board-certified neuropathologist (S.D.) were consistent with those based on 5-ALA fluorescence, MRI, and surgical evaluation, but indicated that most samples, even those from the tumor edge, contained some regions of high tumor cell density. Within each sample, manual histologic annotation indicated regions of high tumor cell density (consistent with tumor “core”), low tumor cell density (infiltrative tumor “edge”), microvascular proliferation (MVP), palisading necrosis, predominantly non-malignant white matter, and predominantly non-malignant gray matter ( Fig. 1B , Fig. S2). Spot-level analysis of Visium data indicates that the tumor edge and core are transcriptionally distinct Spot-level analysis of Visium ST data yields limited insight on its own, but here served as a preliminary step in more nuanced downstream analyses. Graph-based clustering of the spot-level data (Methods) revealed features that were apparent in the H&E images, such as microvascular proliferation, and also captured extensive transcriptional heterogeneity that was not evident in those images. Low-resolution clustering (e.g. 0.1) captured high-level features of GBM shared across many samples (Figs. S3 A-C), but missed more subtle signatures present in subsets of samples. At higher clustering resolutions (e.g. 0.9), more extensive transcriptional variation was apparent, including shared and sample-specific clusters, as well as specific features enriched at either the tumor core ( Figs. 1C,E ) or edge ( Figs. 1D,E ) (Fig. S4). Tissue slices resected from the transition zone contained a mixture of clusters but typically resembled the core more closely than the edge at a transcriptional level ( Fig. 1E ). Each cluster was classified as edge-associated, core-associated, interface-associated (between core and edge), or spatially non-specific (“N.S.”) based on predominant histologic location as defined by the neuropathology tumor content assessments. The “Interface” region was not specifically noted by the pathologist, but interface-associated clusters typically appeared as bands between the core-associated and edge-associated clusters, and were thus assigned their own category. Core-associated clusters (e.g. C3, C6, etc.) were enriched for transcripts associated with several pathways such as hypoxia (C3, Fig. 1F ), microvascular proliferation (C6, Fig. 1G ), cilia function (C8), inflammation and antigen presentation (C9), myeloid leukocyte migration and activation (C12), cellular respiration (C17), response to virus (C19), TGF response (C20), and RNA metabolism (C22). In contrast, edge-associated clusters were enriched for genes associated with oligodendrocyte-related functions such as myelination (C0, C1), and neuron-related functions such as synaptic signaling (C2, C13) ( Fig. 1H ). These findings highlight large-scale differences in gene expression between the tumor core and edge, presumably reflecting spatial heterogeneity in cellular composition. Single-cell mapping reveals distinct cell type composition at the tumor core and edge To analyze spot-level transcriptomic data at higher resolution, the cellular composition of each spot must be computationally inferred. We tested many deconvolution tools (Methods) and chose a random forest-based approach called CellTrek 32 that maps cells from a single-cell reference to x and y coordinates within the Visium data. We created a single-cell reference by generating scRNA-seq data from eleven GBM samples, including two previously published samples 33 , and combining it with published data 30 from non-malignant brain cell types that are often under-represented in GBM scRNA-seq data, including neurons, glia, endothelial cells, and immune cells (Methods) ( Fig. 2A ). The resulting deconvolved cell mappings revealed distinct regions enriched for tumor, immune, endothelial, neuronal, and glial cells. These were associated with the location of the tissue sample—core samples were enriched for tumor cells and infiltrating immune cells, particularly macrophages and often distinct pockets ( Fig. 2B ), while edge samples were enriched for neurons and glia, as expected ( Fig. 2C ), and transitional samples displayed a mixture of cell types but were generally more core-like ( Fig. 2D ). Thus, the deconvolved ST data was consistent with both pathology annotations and functional enrichment of spot-based clusters. Download figure Open in new tab Figure 2: Cell mapping and derivation of spatial tumor cell states. (A) UMAP projection of scRNA-seq data used as a reference for cell mapping. (B) Cell mapping results for representative core tumor sections showing dense tumor cellularity, vascularization, and pockets of immune infiltration. (C) Cell mapping results for representative edge sections showing minimal tumor content and abundant neurons or glia. (D) Cell mapping results for representative transitional 5-ALA+ sections showing dense tumor cellularity, vascularization, and pockets of immune infiltration. (E) Consensus clustering of marker genes for spatial tumor cell states derived from mapped tumor cells. (F) Spatial tumor cell states in core tumor slices showing abundant Hypoxic (S4) and Mesenchymal (S3) tumor cells. (G) Spatial tumor states in edge samples, showing enrichment for Proliferating (S2) and Stem-like (S5) states in regions of high neuronal/glial content. (H) Spatial tumor cell states in transitional 5-7866ALA+ tumor slices showing a gradient from core Mesenchymal/Hypoxic to edge (Proliferative/Stem-like) states. GBM tumor cell states are defined by both gene expression and spatial localization Several studies have employed elegant computational methods to classify GBM cells into distinct transcriptional cell states based on scRNA-seq data 9 , 12 – 14 , 34 . Although common themes have emerged, each study has yielded a different set of cell states. ST offers potential benefits over scRNA-seq for cell state inference because it does not introduce dissociation-associated artifacts, captures a more representative set of cells, and is more likely to reflect the true state of cells that are programmed by their spatial context 9 . We therefore asked whether the spatial organization of the mapped tumor cells could refine our understanding of tumor cell states. To this end, we clustered the mapped tumor cells in each sample using a spatially-weighted consensus clustering approach (Methods, Fig. S5). Each cluster was associated with unique marker genes, and we found that similar clusters recurred across samples with only subtle differences in marker genes. We then used consensus clustering and hierarchical clustering to derive “consensus” clusters defined by the shared features in these sample-specific clusters. This yielded five gene sets, or “spatial tumor cell states,” that defined similar clusters across samples ( Fig. 2E , Methods). Functional enrichment of the five consensus clusters indicated association with: (1) histone methylation and Sox2 interactions (e.g. KMT2A , KMT2C , NR2F1 , SPAG9 ); (2) cell cycle and proliferation (e.g. STMN1 , CDKN2C , TUBA1B ); (3) apoptosis, wound-healing, antigen presentation, and mesenchymal cell state (e.g. VIM , CD99 , CHI3L1 , CD63 , SPP1 , ANXA2 , HLA* ); (4) response to hypoxia (e.g. EGLN1/3 , PGK1 , SOD2 , VEGFA , MIF ); and (5) regulation of gliogenesis and neuronal stem cell differentiation (e.g. SOX4/6/8/11 , PTPRZ1 , NRXN1 , VCAN , BCAN , PDGFRA , OLIG3 , DLL3 ) (Table S3). We refer to these states as Histone Methylation (S1), Proliferating (S2), Mesenchymal (S3), Hypoxic (S4), and Stem-like (S5). Although this classification is unique by incorporating both transcriptional and spatial elements, it combines aspects of previously published classifications 9 , 23 , which have included mesenchymal-like states, stem-like states, neuronal states, proliferating states, and hypoxic states, as well as many of the marker genes listed above. Next, each tumor cell mapped by CellTrek was assigned to one of the five states S1-S5 based on marker gene expression, revealing mesenchymal regions, hypoxic regions, and mixed regions containing stem-like and proliferating tumor cells ( Fig. 2F-H , Methods). Tissue samples from the tumor core, edge, and transition zone revealed unique distributions of tumor cell states. The same scoring procedure was also performed on the original single-cell reference data, demonstrating that these states were also present in scRNA-seq data. Through this analysis, we provide a new classification of GBM cells that is based on gene expression and spatial organization, and show that these states reside preferentially in different regions of the tumor. This novel approach to defining tumor cell states highlights the importance of ST technologies in understanding cell state and intratumoral heterogeneity. The tumor core and edge are composed of distinct multicellular niches with characteristic tumor cell compositions We next sought to understand how these tumor cell states co-localized with other cell types by defining cellular “niches,” or co-localized multicellular communities. To do so, we first used CellTrek to infer the cellular composition of each oligonucleotide spot (and thus each spot-based cluster) (Methods). This enabled us to interpret each spot cluster as a “niche” composed of cells of different types ( Fig. 3A ). Based on the most statistically enriched cell type in each niche, we identified four major niche categories: Tumor, Vascular, Immune, and Neuronal/Glial. For each major category, we observed clear differences in the spatial distribution of these niches. The tumor core contained Tumor niches (e.g. C3, C20), Immune niches (e.g. C16, C9), one Vascular niche (C6), and one Glial niche (C22) ( Fig. 3A ). Core Tumor niches were more highly enriched for Mesenchymal (S3) and Hypoxic (S4) tumor cell states than other regions, and Hypoxic tumor cells were primarily found in regions of highest tumor cell density ( Figs. 3B-C ). The two vascular niches displayed markedly different patterns of spatial organization, suggesting two types of vasculature in GBM: the core-associated vascular niche (C6) matched with histologically-identifiable regions of MVP and also enriched for Mesenchymal (S3) tumor cells; meanwhile, the other vascular niche (C15) corresponded to isolated spots that were dispersed throughout the tissue, not identifiable in the H&E image, and not enriched for any particular tumor cell state. The most immune-rich cluster in the tumor core (C12) was markedly enriched for diverse myeloid and lymphoid cells and strikingly Hypoxic (S4) cells as well, demonstrating co-localization of immune cells and hypoxic tumor cells. Clusters at the edge/core interface (C4, C7, C19, and C21) showed a broad distribution of cell types. As expected, edge-associated niches were enriched for glia (C0, C1, C23) or neurons (C2, C13). Unlike the tumor core, the tumor edge was enriched for Stem-like (S5), Proliferating (S2), and Histone Methylation (S1) tumor states. The enrichment of stem-like and proliferating tumor cells within the edge was supposed by scRNA-seq performed on matched core and edge samples from two GBM patient samples ( Figs. 3D-F ). Download figure Open in new tab Figure 3: Tumor niches defined using cell mapping and clustering. (A) Composition of each graph-based cluster (columns, ordered by tumor location and tumor cell percentage) with respect to mapped cell types (rows). Edge/Core/Transition classifications were made by the surgeons and refined by the pathologist based on tumor cell density. Edge/Core/Interface/N.S. classifications were based on cluster distribution within the tissue. Circle size represents cell-type enrichment, and black borders represent statistically significant enrichment. Enrichment was calculated with respect to all cell types. The niche designation is indicated above each column based on the most-enriched cell type. (B) Heatmap showing statistical enrichment of each spatial tumor cell state (S1-S5, columns) in each niche (rows). Here, enrichment was calculated with respect to tumor cells only. (C) Fisher exact test for enrichment of Tumor states between domains. Dot color represents the bonferroni corrected p-value of enrichments and dot-size is the fold change of proportion of tumor states (w.r.t to tumor cells) between domains. (D) Heatmap showing spatial cell states and regional enrichment in scRNA-seq data from the core and edge of two tumors. (E) Left, UMAP layout of edge/core scRNA-seq data with cells colored by sample. Right, UMAP layout of edge/core scRNA-seq data with cells colored by region. (F) Distribution of scores for each cell state in the edge (blue) or core (red). (G) Spatial proximity map of niches, where edge width is proportional to Weighted Mean Adjacency and edge color is proportional to the fraction of samples in which each niche pair is adjacent. We next quantified the spatial proximity of the niches in order to build a niche “map” ( Fig. 3G ) (Methods). This analysis showed that Hypoxic (S4) and Mesenchymal (S3) tumor niches, TAM niches, and MVP niches were mutual neighbors within the tumor core, while neuronal and glial niches were mutual neighbors at the tumor edge. Niches at the core-edge interface were near core and edge niches, but not necessarily near each other. In particular, the mixed Stem/Mesenchymal (S5/S3) tumor niche C7 connected the edge and core, reflecting the gradient of tumor state enrichment from core to edge. Unexpectedly, the spatially non-specific niches (C5, C15, and C10), which were abundant in immune cells and non-MVP blood vessels, tended to be closer to edge niches. These results highlighted the enrichment of pathologic blood vessel development within the core, but relatively normal vascular development at the tumor edge. Overall, we found that the tumor core and edge are comprised of distinct multicellular niches. Moreover, Mesenchymal (S3) and Hypoxic (S4) tumor cells were enriched in core-associated niches, while Histone Methylation (S1), Proliferating (S2), and Stem-like (S5) tumor cells were enriched in edge-associated niches. Due to the infiltrative nature of GBM, it will likely be vital to develop therapies that target such tumor cells in the edge in order to prevent recurrence. Ligand-receptor (LR) interactions mediating communication within and between niches Similar spatial relationships in GBM, including highly organized hypoxic regions and stem-like/proneural tumor states at the tumor edge, have been observed in several separate studies, but with key differences in either relative co-localization or composition 23 , 28 , 35 , 36 . However, the mechanisms driving such organization have yet to be fully elucidated. Ligand-receptor (LR) interactions are a powerful determinant of spatial organization, and there is significant clinical interest in leveraging them therapeutically. Therefore, we used the stLearn algorithm to identify LR interactions underlying the spatial organization at spot-level resolution. We first assessed contact-dependent LR interactions occurring within the same niche (intra-niche), which revealed highly niche-specific communication patterns ( Fig. 4A ). The vascular niches (C6, C15) were strikingly enriched for the pro-angiogenic NOTCH signaling pathways, and each used distinct combinations of DLL , JAG , and NOTCH genes. Meanwhile, edge-associated neuron- and glia-rich niches (e.g. C2, C13) were enriched for diverse NRXN -mediated signaling pathways, which were recently found to mediate pro-tumorigenic electrophysiological communication between neurons and GBM tumor cells 37 . More surprisingly, the core niche enriched for Proliferative (S2) tumor cells (C8) was significantly involved in NRXN2 -mediated communication, suggesting that pro-tumorigenic tumor/neuron interactions occur even in tumor-dense regions, and may preferentially involve S2 tumor cells. The TAM-enriched core niche (C9) was highly enriched for MHC Class II interactions, suggesting potential CD4 driven responses within the tumor core. Download figure Open in new tab Figure 4: Ligand-Receptor analysis of Visium Samples. (A) Enrichment dotplot of significant contact-dependent (within-spot) ligand receptor interactions in visium data. Size of dot indicates enrichment score, bold circle around spot indicates significance based on adjusted p.values. Row Barplots show celltype diversity for the said LR whereas column barplots indicate celltype diversity in the cluster. The clusters are ordered from edge rich clusters to core-rich clusters. (B) LR connectivity (C i , C j ) for diffusion based LR’s is shown as a cytoscape network between clusters, indicating a higher connectivity between core-rich clusters than edge-rich clusters. Circos plots on the network nodes indicate sender and receiver celltype frequency. (C) Enrichment dot plot based on hypergeometric test for diffusion-based LR interactions between groups of clusters, where dot color indicates values for -log10(p.adj+1) for LR’s that are significant in at least 1 group. The dot size indicates enrichment score and the black boundary around dots indicate significant LRs. The attached heatmap on the right indicates celltype frequency (from celltrek) for the significant spots in each group combination (arranged from edge to core similar to the dotplot). (D) Similar to B but only inter-niche interactions are shown here in 2 directions. The first direction being from edge towards core and the second being from core towards edge. (E) Summary of the inter-niche LR signalling results. Circles are colored by groups. Next, to investigate LR interactions involving diffused/secreted ligands occurring between different niches (inter-niche), we generated a high-level map which largely reflected the spatial relationships in Fig. 3G and indicated that Mesenchymal (S3) tumor cells were a hub of intercellular communication ( Fig. 4B , Methods). Predicted inter-cellular communication was abundant within the core, with extensive crosstalk within and between Mesenchymal (S3) niches, TAM niches, Hypoxic (S4) niches, and regions of MVP. The majority of the LR interactions involved in this crosstalk included those known to drive angiogenesis in GBM, specifically ANGPTL4- , SDC- , TNC- , THBS1- , ADM- , FGF- , VEGFA- , and VEGFC- mediated interactions, CXCL12 : CXCL4 , and WNT4 : FZD * (Supplementary Data). A given biological process (e.g. angiogenesis) can be driven by different ligand-receptor interactions depending on spatial context We were then interested in how these diffusion-based LR interactions may change as a function of localization within the broad domains defined above – tumor edge, core, edge/core interface, and spatially non-specific regions ( Figs. 4C-D , Supplementary Data). Beginning with diffusion-based LR interactions that occurred within the same domain, this analysis revealed extensive pro-angiogenic signaling in the tumor core that utilized diverse LR pairs, including some ligands canonically associated with angiogenesis (e.g. ANGPTL4 , PDGFB, and ADM ). Meanwhile, the interface was enriched for interactions including RELN:LRP8, RELN:VLDLR, and SEMA3B:PLXNA3 , all of which mediate important neuronal signaling interactions and whose presence suggest tumor-neuron interaction in these specific regions 38 , 39 . Major pro-angiogenic signaling interactions observed in the core were not statistically enriched in the interface, notably, supporting our previous observations that pathologic microvascular proliferation was not observed in these regions. Surprisingly, we found a unique pro-angiogenic signaling pathway in the edge, APLN:APLNR , which was not enriched in any other major domain. Taking pro-angiogenic pathways as an example, these results highlight the complexity of ligand-receptor interactions within the GBM tumor microenvironment and show that different ligands and receptors may drive a common biological process in different regions of the tumor. In addition, LR interactions that were present in multiple tumor regions sometimes involved different cell types in each region. Ligand-receptor (LR) interactions mediating communication among broader spatial domains Finally, we sought to understand intercellular communication between different spatial domains ( Figs. 4D-E ). This analysis showed that the majority of interactions were bidirectional, e.g. the core signaled to the interface via PTN:PTPRZ1 and vice versa. Several of these interactions included ADM:RAMP3 , a pro-angiogenic signaling interaction, and PDGFB:PDGFBR and FGF1:FGFR3 , two important growth signaling interactions. Unidirectional interactions were also observed, such as C3:C3AR1 , in which the interface signaled unidirectionally to the edge. Meanwhile, in a reciprocal fashion, the edge signaled unidirectionally to the interface via SEMA3B:PLXNA3 . Overall, these results highlight the abundant interactions among all domains ( Fig. 4E ) and the complexity of interactions mediating important biological processes within the tumor microenvironment. As exemplified by pro-angiogenic intercellular interactions, different ligands and receptors are utilized depending on spatial localization, highlighting the need for multi-target drug regimens, rather than single target therapies. Xenium in situ transcriptomics recapitulates spatial tumor states derived from Visium data and highlights the presence of key immune populations To supplement our ST results, we employed Xenium, a panel-based technology that quantifies transcript abundance at high resolution (<30 nm XY precision and <100 nm Z localization precision). We extended the pre-designed Brain/GBM Xenium panel (10X Genomics) with custom oligos, ultimately targeting 362 cell markers (Methods and Table S4). We used this platform to generate high-quality data for nine sections of primary GBM tumors predominantly from the tumor core (Table S1& S5). We generated Xenium data for two patient cohorts that were analyzed using different cell segmentation technologies (Methods): the primary cohort contained five samples from four different patients; the secondary cohort contained four samples from one patient. Cell typing was performed on the primary cohort at single-cell resolution (Materials and Methods), and cell types were superimposed on the secondary cohort using transfer learning ( Fig. 5A-B , Fig. S7 A-B). Our custom probe panel enabled us to annotate numerous key cell populations, including tumor cells, immune cells, vascular-related cells, neuronal cells, and glia (Fig. S7 C, Table S4-S6). We identified several myeloid populations that were consistent with macrophages, but due to the limited number of probes in our Xenium panel, we conservatively referred to these as “myeloid” cells rather than “macrophages.” Two of these myeloid populations (C7, C13) expressed low levels of MHC class II and were consistent with E-MDSC and M-MDSCs, respectively, observed in a recently published study 25 . Both CD8 + and CD4 + GZMK-expressing T cells were also observed and consistent with our previously published data, expressing low levels of IL7R and GNLY 20 . Several neuronal clusters were also detected, with gene expression profiles consistent with several different types of excitatory and inhibitory neurons. Notably, sections from different regions of the same tumor showed heterogeneous cell composition in both cohorts (Fig. S8 D). Download figure Open in new tab Figure 5: Xenium-based Analysis of GBM Tumor Samples. (A) UMAP visualization of all cells from primary cohort samples. (B) Spatial plots of GBM026I and GBM030 with cell types colored as per Fig. 5A . (C) Hierarchical clustering of tumor clusters from both cohorts (primary and secondary) based upon expression of Visium-based tumor signatures. (D) Hierarchical clustering of BANKSY clusters from both cohorts (primary and secondary) based upon normalized values of cell type frequency. (E) Spatial plots of GBM026I and GBM030 with niches colored. (F) Dot plot of enrichment scores for each cell type per respective niche. The enrichment score was calculated by dividing the cell type frequency per niche by the cell type frequency within the entire data set and signified by the size of each circle. Significance was calculated using the hypergeometric test and denoted by a black border around each respective circle. (G) Node-edge plot visualizing the average distance between niches across all samples (in both cohorts). Distances between niches were calculated on a per sample basis and then averaged across samples if applicable. The number of samples with a given niche-niche distance is signified by the circled number. In (B) and (E) , specific regions of interest (ROIs) are highlighted by dashed white boxes with the respective ROI matched by the marked number in the top left or bottom right corner. Next, we verified the presence of the Visium spatial tumor cell states in the orthogonal Xenium data set. Briefly, Louvain-based clusters of tumor cells were scored for expression of Visium-based signatures, then hierarchically clustered. All Visium-based tumor cell states were recapitulated in the Xenium data except for the Histone Methylation state (S1), which was the rarest cell state in both the Visium data and scRNA-seq data ( Fig. 5C ). We also detected a Xenium-specific “Astrocytic” tumor cell state (S6) marked by high RYR3 expression, and two “meta states,” M1 and M2, which were defined by combinations of Visium-based signatures. M1 represented a combination of the Proliferative (S2) and Stem-like (S5) tumor cell states, and M2 a combination of the Mesenchymal (S3) and Hypoxic (S4) tumor cell states, which we interpreted as a “Hypoxic Mesenchymal” cell state and has been previously reported 9 . Thus, tumor cell states derived from Visium data were largely recapitulated by the single-cell resolution Xenium data. Xenium facilitates single-cell analysis of niche composition and suggests additional complexity within hypoxic and vascular niches Tumor sections analyzed using Xenium exhibited complex, visually striking spatial organization ( Fig. 5B ). To characterize these patterns quantitatively, we used BANKSY 40 to define niches. We first defined niches on a per-sample basis, then unified them across samples by clustering the sample-specific niches for both cohorts ( Fig. 5C-D , Table S8, Materials and Methods). Overall, we characterized eight niches marked by variable enrichment for tumor, immune, vascular, and glial/neuronal cells. Each niche exhibited characteristic tumor cell state enrichment: Hypoxic.1 (enriched for the M2 state), Hypoxic.2 (M2- and S4), Stem-like (M1/S5), Mixed.1 (M1/S4), Mixed.2 (S3/S5), Immune-Infiltrated (S3), and Neuronal (S6). The Vascular niche was defined by high levels of endothelial, mural, and immune cells, and contained very few tumor cells. While some niches were found in most samples, such as the Hypoxic and Vascular niches, others were more sample-specific, such as the Immune-Infiltrated and Mixed niches (Fig. S7 E). We next calculated the cell-type enrichment of each niche, which captured their unique cellular “fingerprints” ( Fig. 5F ). Notably, the Immune-Infiltrated niche was enriched for several immune cell types in addition to Mesenchymal (S3) tumor cells, consistent with the Visium data. Although colocalization of mesenchymal tumor cells with macrophages has been previously described 23 , 28 , this data suggests that mesenchymal tumor cells also colocalize with T cells. We also observed T cell enrichment in the Hypoxic.1 niche, which is notable because GBM tumors are typically considered “immune-cold” 41 , with low T cell infiltration and high T cell dysfunction. These data suggest that T cells are present but sequestered within specific niches, which may significantly impact their function. Furthermore, we found Monocyte-like (C13) cells, similar to M-MDSCs, were enriched in the Immune-Infiltrated niche while Hypoxic.Myeloid.1 (C7) cells, similar to E-MDSCs, were enriched in both hypoxic niches. These observations were consistent with previous findings 25 and highlight an additional level of immunosuppression within the tumor microenvironment, particularly in tumor-rich niches enriched in T cells. Both hypoxic niches were enriched for Hypoxic (S4) tumor cells (by definition) and multiple myeloid populations. The Hypoxic.1 niche was also enriched for KIT-expressing endothelial cells, which have been shown to play an important role in driving angiogenesis in GBM 42 , 43 . These results suggest that KIT-endothelial cells may function as the progenitor cell from which microvascular proliferation may arise either near, or even within, hypoxic regions. Stem-like and Mixed niches were primarily enriched for tumor cells, and secondarily for myeloid populations (e.g., CCL4 + myeloid cells in Mixed.1). Meanwhile, the Neuronal niche was enriched for glutamatergic and inhibitory neurons, glial cells including oligodendrocytes and astrocytes, M1 (i.e. Proliferating (S2) and Stem-like (S5)) tumor cells, and the unique Xenium-specific population of astrocytic tumor cells (S6). The only niche not significantly enriched for any particular tumor cell state was the Vascular niche, populated by diverse immune cells, endothelial cells, and mural cells. Overall, these observations corroborated three important findings from our Visium-based analysis: 1) Different tumor cell states reside in different multi-cellular niches within the tumor, 2) Mesenchymal (S3) tumor cells preferentially co-localize with immune cells, and 3) Stem-like (M1/S5) tumor cells preferentially co-localize with neuronal cells Xenium analysis identifies niches along a spatial gradient defined by proximity to vasculature We next sought to quantify the spatial relationships among the Xenium-derived niches by quantifying the average minimum distances between niches across samples ( Fig. 5G , Table S9, Methods). This analysis revealed striking differences in the proximity of the tumor-enriched niches to the Vascular niche. The organization of tumor niches around the vascular niche were as follows from closest to farthest: the Mixed.2 niche (which was enriched for Mesenchymal (S3) and Stemlike (S5) tumor cells), the Stem-like niche (enriched for Stem/Proliferating (M1) and Stemlike (S5) tumor cells), the Immune Infiltrated niche (enriched for Mesenchymal (M3) tumor cells), the Mixed.1 niche (enriched for Stem/Proliferating (M1) and Hypoxic (S4) tumor cells), and finally the hypoxic niches (Hypoxic.1 which was enriched for just Hypoxic Mesenchymal (M2) or Hypoxic.2 which was enriched for Hypoxic Mesenchymal (M2) as well as Hypoxic (S4) tumor cells)—i.e., niches enriched in stem-like or mesenchymal tumor cells were closest to the vasculature while hypoxic tumor cells were further. Although both hypoxic niches were far from the vascular niche, they were also localized away from each other, on the magnitude of several hundred micrometers. These results further highlight the different hypoxic niches that exist within GBM – not only do Hypoxic.1 and Hypoxic.2 niches have different cell type enrichment profiles, but they were also spatially localized in different regions within the tumor microenvironment. Single-cell LR analysis reveals a central role for TAMs, vasculature, and Mesenchymal tumor cells in intratumoral communication Next, we performed spatially-informed LR analysis using CellChat to determine which LR pairs mediate interactions among (a) cell types and (b) niches. This revealed the overall structure of the GBM ecosystem from the perspective of intercellular interactions, cell-state-specific LR interactions, and niche-specific cell states not previously observed. Significant LRs between every cell type were aggregated across samples and both unique and total number of LR interactions was quantified (Table S10, Materials and Methods). This highlighted the striking diversity of LR interactions utilized by non-tumor cells such as vascular, myeloid, glial, and neuronal cells ( Fig. 6A,B ). While vascular and myeloid cells utilized these diverse LRs primarily in the context of heterotypic interactions, glial and neuronal cells utilized them primarily in homotypic interactions. Within the tumor compartment, Mesenchymal (S3) tumor cells utilized the greatest number of unique LRs, as both senders and receivers, and in both heterotypic and homotypic interactions. Download figure Open in new tab Figure 6: Ligand-Receptor Analysis of Xenium Samples. (A) Interaction network illustrating unique ligand receptor interactions between general cell types. Only unique ligand receptors detected in more than 3 samples were visualized.The color, opacity, and thickness of arrows are proportional to the number of unique LR interactions connecting any two cell types. (B and C) Heatmaps of the (B) unique and (C) total number of ligand receptor interactions between general cell types. (D) Ligand-receptor interactions of tumor cells (Hypoxic: S2 or S4, Stem-like: S1 or S5) as senders and myeloid, neuronal/glial, and vascular cells as receivers. Ligand-receptor pairs unique to a given sender-receiver pair were italicized and labeled with an asterisk. (E) Interaction network of cell type sender-receiver pairs of the ADM:CALCRL ligand-receptor pair. The color and density of arrows were proportional to the number of samples the ADM:CALCRL interaction was detected as significant for a given sender-receiver pair. (F) Heatmap of ADM and CALCRL expression per cell type per niche. (G) Expression of ADM and CALCRL based upon localization in tumor region. Data derived from IVY-GAP. We then considered the total number of LR interactions between pairs of cell types ( Fig. 6C , Fig. S8 A), which is influenced by cell-type frequency. For each pair of cell types and each LR, we defined an interaction as above, then added those interactions across instances within the same sample, across LR pairs, and across samples. As above, this highlights the primacy of the vascular cells, myeloid cells, and Mesenchymal (S3) tumor cells on the communication landscape. In particular, Mesenchymal tumor cells sent more signals to other cells (tumor and non-tumor) than any other tumor cell state, particularly myeloid and vascular cells (Fig. S8 B). LR analysis reveals hypoxic niche-specific cell states and interactions among tumor, immune, vascular, and neuronal/glial cells Recent studies have highlighted the specialized roles of different tumor transcriptional states in shaping the tumor microenvironment, especially with respect to myeloid cells, vasculature, and neuronal cells 23 , 28 , 44 – 46 . This prompted us to ask whether the different tumor cell states in our data utilized different LRs to communicate with vascular and myeloid cells ( Fig. 6D ). We found that Mesenchymal (S3) tumor cells were unique in utilizing CNTN1:NOTCH1, GJA1:GJA1, and LAMA2:CD44 interactions to communicate with vascular cells; the LAMA2:CD44, CD99:PILRA, and CD99:CD99 interactions to communicate with immune cells; the CNTN1:NRCAM interactions to communicate with neuronal/glial cells. Meanwhile, the hypoxic tumor states (M2 and S4) were unique in utilizing ADM:CALCRL to communicate with both vascular and neuronal/glial cells. Finally, stem-like tumor cells (M1 and S5) only uniquely interact with neuronal/glial cells through MPZL1:MPZL1 and NRXN1:NLGN1 . Several LRs were utilized by all tumor cell states to communicate with vascular or myeloid cells, including ANXA1:FPR1 and CSF1:CSF1R interactions between tumor and myeloid cells, and VEGFA:VEGFR1 and PTN:NCL interactions between tumor and vascular cells. Meanwhile, several previously reported interactions between tumor and neuronal/glial cells were shared, such as NRXN1:NLGN3 and SEMA6A:PLXNA4, though only by Mesenchymal (S3) and stem-like tumor cells. Thus, some LR pairs were broadly utilized by all tumor cell states to interact with the surrounding microenvironment, while others were utilized by specific tumor states, which provide potentially unique targets. Given that our Visium results highlighted abundant pro-angiogenic interactions within the tumor core that included ADM:CALCRL, we decided to analyze this pathway in greater detail within the Xenium data set because recent work 26 , 47 , 48 suggested that (1) ADM expression by myeloid cells can lead to neovascularization and vascular dysfunction in GBM, and (2) therapeutic targeting of ADM improves drug delivery and efficacy. Through the single cell resolution of the Xenium technology, our analyses have added significant nuance to this story. First, the data suggested that hypoxic myeloid cells and hypoxic tumor cells use ADM:CALCRL to communicate in a unidirectional manner with vascular cells, especially KIT+ endothelial cells ( Fig. 6 D-E). The expression patterns of ADM and CALCRL were also recapitulated in our scRNA-seq data set (Fig. S8 C-D). This finding was consistent with the known role of hypoxia in inducing ADM expression via the transcription factor HIF1A. To understand how these interactions depended on the broader spatial context, we then investigated the ADM:CALCRL interaction at the niche level, and found that only one hypoxic niche (“hypoxic.1”) interacted with the vascular niche via ADM:CALCRL in multiple samples (Fig. S8 E). Importantly niche-specific expression of ADM and CALCRL seemed to underlie the unidirectionality of this inter-niche interaction ( Fig. 6F ). For example, Hypoxic Mesenchymal (M2) tumor cells expressed ADM more highly when they were located in the hypoxic.1 niche than in the vascular niche; similarly, endothelial.2 cells expressed CALCRL more highly when they were in the vascular niche than in the hypoxic.1 niche. We observed similar region-specific expression of ADM and CALCRL in two orthogonal data sets: our Visium data, and bulk RNA-sequencing data from the Ivy Glioblastoma Atlas Project (Ivy-GAP) 49 . In the Ivy-GAP data, ADM was higher in the perinecrotic zone and pseudopalisading cells, regions associated with hypoxia, while CALCRL was higher in hyperplastic blood vessels and microvascular proliferation ( Fig. 6G ). With these analyses, we were able to show that domain-specific ligand-receptor interactions observed within the Visium data were recapitulated at a single cell level and likely driven due to spatially-restricted gene expression profiles. GZMK + CD8 + T cells co-localize with LYVE1 + CD163 + myeloid cells in vascular niches In recent work, we and others found that GZMK + CD8 + T cells were the most abundant T cell population in several tumor types 20 , 50 , especially GBM, and were selectively clonally expanded in the tumor compared to matched peripheral blood. GZMK+ CD8 + T cells exhibited a highly non-uniform spatial distribution in our ST data, so we sought to further characterize the cellular environment surrounding GZMK + CD8 + T cells, or GZMK+ CD8+ T cell niches. To this end, we defined an enrichment score to identify cell types enriched within a given radius surrounding each GZMK + CD8 + T cell ( Fig. 7A ) (Methods). This showed that the most highly enriched cell types within a 0.5-20 µm radius of GZMK + CD8 + T cells included CD8 + T cells, CD4 + T cells, NK cells, and LYVE1 + CD163 + myeloid cells. The proximity of LYVE1 + CD163 + myeloid cells was striking, because these cells are thought to regulate migration of T cells into tumors, although the directionality of their influence is ambiguous 51 – 53 . In the homeostatic brain, LYVE1 + CD163 + myeloid cells resided within the perivascular space 54 – 56 and given T cells accumulate within the perivascular space in gliomas 57 , 58 , we hypothesized that the colocalization of LYVE1 + CD163 + myeloid cells and GZMK + CD8 + T cells may be specific to the vascular niche in GBM. To test this, we used several approaches to quantify the distance of each GZMK + CD8 + cell to the nearest LYVE1 + CD163 + myeloid cell, and found that these cells were significantly closer to each other in the vascular niche compared to non-vascular niches (22.45 µm vs. 156.7 µm, p.value < 0.0001, Mann-Whitney U-test) ( Fig. 7B-E , Methods). This relationship was not seen with CD4 + memory or CD4 + GZMK + T cells across samples, and was not due to differences in the distribution of lymphoid cells and LYVE1 + CD163 + myeloid cells across the various niches (Figs. S9 A-C). We finally used an orthogonal data set to support the spatial relationship between GZMK + CD8 + T cells and LYVE1 + CD163 + myeloid cells. In bulk RNA-seq data from IVY-GAP, GZMK and LYVE1 expression were significantly correlated (R 2 =0.3418, p<0.001) in vascular regions only ( Fig. 7 F ). Notably, GZMB was not correlated with LYVE1 in either vascular or non-vascular regions (Fig. S9 D). Download figure Open in new tab Figure 7: CD8 + GZMK + T cells preferentially localize with LYVE1 + CD163 + myeloid cells in vascular niches. (A) Top 6 cell types enriched within 20 micron radius of CD8+ GZMK+ T cells. Enrichment score was calculated by dividing the frequency of a target cell type within a given radius of a root state (CD8+GZMK+ T cells) by the frequency of the target cell type within the sample. Each dot represents a patient sample. (B) Distance to the nearest target myeloid cell of each CD8 +. GZMK + T cell separated by localization within a vascular or non-vascular niche. Each dot represents a single CD8 +. GZMK + T cell, with all CD8 +. GZMK + T cells across all samples visualized. (C) Distances from Fig. 6B averaged on a pers sample basis. Each dot represents an individual sample. Significance was calculated using the Wilcoxon matched-pairs signed rank test. (D) Frequency of CD8 +. GZMK + T cells in the vascular niche given proximity to the target myeloid cell (either within 30 μm or further than 30 μm). (E) Frequency of CD8+GZMK+ T cells within 30 μm of target myeloid cell given presence of T cell in vascular or non-vascular niche. (F) Correlation of GZMK and LYVE1 expression based upon localization in vascular or non-vascular niche. Data derived from IVY-Gap. (G) Confocal imaging of immunofluorescence staining of human GBM tumor sample. Surprisingly, LYVE1 + CD163 + myeloid cells were not detected in our scRNA-seq data set nor many other published scRNA data sets, though a recent publication highlighted LYVE1 as a contributing gene of the “scavenger suppressive” profile in glioma myeloid cells 59 (Fig. S9 E). As a result, to validate the presence of these myeloid cells, we performed immunofluorescence imaging on human GBM tumor sections and found that LYVE1 and CD163 were present in human GBM tumors at the protein level ( Fig. 7G ) and co-expressed by the same cell ( Fig. 7H ). Given the potentially immunosuppressive nature of LYVE1 + CD163 + myeloid cells in other tumor types, this interaction may present one possible mechanism by which T cells are sequestered within the perivascular space in the GBM TME. DISCUSSION It is increasingly appreciated that solid tumors function as multicellular ecosystems in which the spatial organization and interactions among diverse malignant and non-malignant cells influence the development and behavior of the tumor, including its resistance to therapy. Using a combination of complementary single-cell and novel applications of spatial transcriptomic technologies, we characterized the spatial organization of GBM and the ligand-receptor interactions that govern intercellular communication across the tumor landscape. This revealed two key features of the GBM ecosystem that have implications for the development of new therapies. First, we found marked differences in cellular composition along a radial axis from the tumor core to the tumor edge. In initial analyses, we derived six base tumor cell states jointly defined by gene expression and spatial proximity, as well as two “meta-states” that are combinations of these base states. Specifically, we found that the distribution of these cell states varied along the core-edge axis, with the tumor edge being comparatively enriched for states associated with Proliferation (S2), Stemness (S5), and Histone Methylation (S1). Meanwhile, we found that the tumor core was comparatively enriched for Mesenchymal (S3) and Hypoxic (S4) tumor cell states. These results support previous studies that have found similar gradients of tumor cell states along the GBM core-edge axis 23 , 28 , 35 , 60 and emphasize the need to understand GBM tumor cell states not only from a transcriptional, but also a spatial perspective. In particular, GBM stem cells have been implicated in treatment resistance and recurrence 61 – 63 . Their enrichment at the tumor edge, which is impossible to fully identify and excise during surgery, motivated us to analyze the intercellular interactions that mediate their communication with other cell types and maintain a distinct ecosystem at the tumor edge. With this in mind, we next identified spatially coherent communities of cells, or “niches.” Most niches were enriched for one or more specific tumor cell states, and co-enrichment of immune cells with Mesenchymal (S3) and Hypoxic (S4) tumor states was particularly evident. Critically, we observed a gradient of niche types (and tumor states) across the core-edge axis: the tumor core exhibited high tumor cell density with pockets of vascular, hypoxic, and immune-infiltrated niches that were typically enriched for Mesenchymal (S3) or Hypoxic (S4) tumor states. In contrast, the tumor edge was enriched for glial and neuronal niches that were often enriched for Stem-like (S5), Proliferating (S2), and Methylation (S1) tumor cell states. Interestingly, we identified two types of vasculature with distinct spatial profiles: abnormal microvascular proliferation (MVP) which was concentrated in the tumor core, and a distinct vascular niche, not associated with MVP, at the tumor edge. Characterizing and understanding the cellular composition of such niches is pivotal to understanding how cell states may drive, result from, and interact with the surrounding environment, providing key insights into potential therapeutic targets. In our study, we highlight several examples: First, the enrichment of T cells within hypoxic niches presents several levels of immunosuppression that may benefit from combination therapies as hypoxia itself can lead to T cell dysfunction 64 – 66 and hypoxic regions are associated with both immunosuppressive macrophages 22 , 27 , 67 and E-MDSCs 25 . As a result, immune checkpoint blockade (ICB) therapies, such as anti-PD-1, may need to be supplemented with therapies that inhibit such immunosuppressive myeloid cells. Hypoxic tumor signatures themselves have also been associated with worse prognosis in patients with GBM 35 . Second, our data provides further evidence that Proliferating (S2) and Stem-like (S5) tumor cells interact with glial and neuronal cells, building upon previous work that has shown that tumor-neuron signaling can drive tumorigenesis 46 . Such niches would therefore require targeting of not only these stem-like tumor cells, but potentially the surrounding neurons or glia as well. Finally, from a strictly immunological perspective, we found significant interaction between CD8+GZMK+ T cells and LYVE1+CD163+ myeloid cells within the vascular niche. Coupled with evidence that LYVE1+CD163+ myeloid cells can inhibit T cell migration into the tumor core 51 , 53 , and such genes define a “scavenger” suppressive profile linked to worse survival and therapy resistance in patients with glioma 59 , this suggests that delivery of targeted therapies to vascular regions of the tumor may improve response to ICB therapies. Taken together, these observations suggested that the tumor edge represents a unique ecosystem that may respond to different therapies from the tumor core; understanding why the edge responds differently to treatment and how it may contribute to recurrence will ultimately allow us to develop multiple arms of therapy that can ultimately eradicate all tumor cells. Second, our data highlighted the extraordinary complexity of ligand-receptor (LR) interactions underlying the spatial organization of the tumor. On one hand, we found that LR interactions vary spatially, and the same cell types may utilize different LR interactions depending on tumor localization. Most prominently, we found that there were several spatially variable LR interactions that all have been shown previously to drive angiogenesis either in GBM or other tumor models. Within the tumor core, VEGFA:VEGFR1 , ADM:CALCRL , PDGFB:PDGFRB, and ANGPTL4 -mediated interactions were significantly enriched. Our data suggested even further nuance within the core, such that ADM:CALCRL was limited to hypoxic regions. Meanwhile, within the tumor edge, APLN:APLNR was specifically and significantly enriched. Such results provide a potential explanation for the ineffectiveness of bevacizumab in patients with primary GBM 68 and limited response in patients with recurrent GBM, 69 given that bevacizumab only targets VEGF-A-mediated interactions. An additional layer of complexity suggested by the Visium data was that the same LR may be utilized by different cell types at different positions along the core-edge axis. For instance, the Visium data suggested that APLN:APLNR mediates interactions between glial and Stemlike (S5) tumor cells at the tumor edge, but among S3/S4/Vascular cells at the tumor core. We observed similar interaction patterns for THBS2:CD47 . In contrast, Visium data suggested that PTN:PTPRZ1 interactions were mediated by Stemlike (S5) and glial cells at the edge; Hypoxic (S4) and Stemlike (S5) cells at the interface; and Mesenchymal (S3), Stemlike (S5), and myeloid cells at the core. In the context of cancer, disrupting ligand-receptor (LR) interactions holds significant clinical potential. Therapies targeting intercellular interactions, such as bevacizumab, are already in use, with numerous additional antibodies and compounds currently in development for this purpose. However, GBM remains a highly lethal disease, and the extensive diversity and spatial specificity of LR interactions observed in our data may explain the limitations of existing therapies. Our findings suggest that targeting a single LR interaction is unlikely to completely disrupt communication between specific cell types or niches. According to the Drug/Gene Interaction database 70 , many of the LR interactions we identified could potentially be targeted by at least one existing drug or compound (Supplementary Data). In addition to spatially restricted LRs, we identified several ubiquitous interactions, including those between SPP1 and the receptors CD44, ITGAV , and ITGA4/5/8/9 , as well as several LR interactions significantly enriched across multiple tissue domains, such as THBS2:CD47 and CD99:CD99 . These interactions may represent particularly promising drug targets due to their broad relevance, and several approved or investigational compounds are believed to inhibit their activity (Supplementary Data). Overall, our study highlights the need for combinatorial therapies that incorporate spatial variation in tumor cell states and heterotypic intercellular interactions which may ultimately offer a more effective approach to improving treatment outcomes for patients with GBM. METHODS Patient recruitment and sample collection Adult patients undergoing neurosurgical intervention were screened for the following criteria: (1) age > 18 years and (2) likely presence of primary glioblastoma with clinical indications for surgical resection. Final diagnosis of primary IDH WT GBM was later confirmed by histopathological analysis and molecular testing. Informed consent from patients was obtained prior to surgery following the IRB Protocol #201409046 (Washington University in St. Louis) or the Dana Farber/Harvard Cancer Center IRB Protocol #10-417 and Massachusetts General Hospital Secondary IRB Protocol #2022P001982. Patient characteristics were summarized in Table S1. Following surgical resection, samples were placed in normal saline and maintained on ice until further processing. All procedures and experiments were performed in accordance with the Helsinki Declaration. Visium data generation and processing Using image-guided surgical biopsies and 5-aminolevulinic acid (5-ALA) fluorescence, we collected 22 samples from tumor regions ranging from the core to the margin in nine patients with IDH WT GBM undergoing surgical resection of radiographic, primary high-grade glioma. Each section (i.e. sample) was preserved, stained with hematoxylin and eosin (H&E), imaged, and processed to generate a spatial transcriptomic library. The first 13 samples were fresh-frozen and OCT-embedded, and cDNA libraries were constructed according to the Visium Spatial 3’ Gene Expression User Guide Rev. D. The next nine samples were preserved in FFPE, and cDNA libraries were generated according to the Visium Spatial Gene Expression for FFPE User Guide Rev. E. All libraries were sequenced on an Illumina Novaseq S2 flow cell to a target depth of 50,000 reads/spot and yielded high quality transcriptome data (median values: 2300 spots under tissue, 55917 mean reads under tissue spot, 2275 panel genes per spot, 4544 median panel UMIs/spot, and 22508 (OCT) or 16334 (FFPE) total genes detected) (Table S2). Demultiplexing, image and transcript alignment, and transcript quantification was performed using the Spaceranger pipeline (10x Genomics, v1.2, GRCh38 reference). For comparability across OCT-preserved and FFPE-preserved samples, all samples were analyzed using the genes in the Visium FFPE gene panel. Spots were filtered to exclude those with fewer than 500 unique transcripts, and initial downstream analysis— including normalization for per-spot sequencing depth using the SCTransform algorithm, principal component analysis, UMAP layout, graph-based clustering, differential gene expression analysis, and spatially variable gene identification — was performed using the R package Seurat (V4) 71 . Previously defined gene expression modules associated with tumor and immune cell states were incorporated into this analysis using the AddModuleScore Seurat function 9 , 72 . Spatially variable biological processes and cell types were also inferred using clusterProfileR 73 , 74 . Annotation and analysis of histological regions H&E images of all tissue sections were manually annotated with guidance from a board-certified neuropathologist (S.D.). In FFPE-preserved samples, regions of microvascular proliferation, necrosis, palisading necrosis, hemorrhage, and normal brain parenchyma were identified (e.g., Fig. 1 , Fig S2). In OCT-preserved samples, H&E quality permitted only the identification of gross features such as regions of high and low tumor content. Tumor cell percentages were visually estimated, with extensive variation noted within and between sections. Importantly, Visium data does not have single-cell resolution; by overlaying the Visium oligonucleotide “spots” on the H&E images, we estimated that each spot contained 5-8 cells. Single-cell RNA-sequencing (scRNA-seq) data generation and analysis We used the 10x Genomics 3’ v2 and 5’ v2 Chromium Gene Expression Solution to generate high-quality scRNA-seq data for 11 GBM samples from nine patients, three of whom also contributed Visium data. Because the size of tumor samples was often limiting, we could not obtain scRNA-seq samples from every patient analyzed using Visium. Sequencing reads were demultiplexed, aligned, and counted using the Cellranger 5.0 pipeline 75 . Downstream processing steps, including filtering, doublet removal, normalization, scaling, identification of variable genes, dimensionality reduction, clustering, and cell typing were performed using a WDL pipeline we developed for tailored analysis of GBM scRNA-seq data ( https://dockstore.org/organizations/PettilabMGH ). Briefly, doublets were identified using multiple tools, including “doubletCells” 76 ,“cxds”, “bcds”, “hybrid” 77 , “scDblFinder” 78 , and “DoubletFinder” 79 . Cells predicted to be doublets by a majority of these were removed from the data. Next, cells were retained for further analysis if they met the following filtering criteria: mitochondrial gene percentage < 5%, nCount (UMIs/cell) 700. Data was then log-normalized, scaled, and corrected for cell cycle gene expression using functions from Seurat V4 80 . The optimal number of Principal components (PCs) for graph-based clustering 71 and UMAP layout was determined based on the Jackstraw test and statistical significance after Bonferoni correction. Cell types were then determined using a combination of reference-based (SingleR) 81 and marker-based (scSorter) 82 tools. To this end, cells were first partitioned into CD45+ (immune) and CD45-(non-immune) clusters based on PTPRC expression. Immune cells were further classified by training these tools on the single-cell Brain immune atlas 83 and the immune component of the STAB 84 atlas (adult samples only). Non-immune cells were further classified by training these tools on a combination of published 9 GBM data and the non-immune component of STAB 84 . Putative tumor cells were further classified based on the cell states in Neftel et al 9 using a marker-based prediction method called scSorter 82 . The malignant status of putative tumor cells was further confirmed using copy number variant (CNV) predictions obtained using CONICSmat 85 . Cells with ambiguous or noisy predictions (e.g. tumor cells without CNVs and nonmalignant cells with CNVs) were excluded from the data set. Generation of a brain/GBM single-cell reference for Visium data deconvolution We combined this data with normal brain cells from the cortical dissections reported in a recent state-of-the-art whole-brain cell atlas 30 and T cells from our recent single-cell study of GBM T cells 20 . Because downstream analysis steps are computationally intensive, we downsized this data by randomly sampling from within each data set. The resulting single-cell reference contained 40,000 cells and 48 cell types. Integration of spatial transcriptomic and scRNA-seq data There are two classes of tools for inferring Visium spot composition, “deconvolution” tools, which estimate cell type proportions for each spot, and “cell charting tools,” which estimate tissue coordinates for individual cells based on a reference single cell dataset. We tested and compared a variety of tools in each class including SPOTlight 86 , Seurat, RCTD 87 , CellTrek 32 , and CytoSPACE 88 ). We considered (1) the extent to which the tools’ predictions are consistent with marker gene expression, graph-based clustering, H&E annotations by pathologist, and other tools; (2) robustness to composition of the single-cell reference; and (3) ability to identify tumor cells in GBM mouse models in which the location of tumor cells is unambiguous (data not shown). Within each category, CellTrek and SPOTlight (V1) worked best. CellTrek enabled more nuanced downstream analysis of the mapped single cells, and was selected for subsequent analyses. Celltrek also allowed us to interpolate cells that may exist between spots and thus was used in two modes “with interpolation” and “without”. Superimposing Visium spots on mapped cells In addition to using Celltrek with interpolation, Celltrek source code was modified to obtain cell types based on the single cell reference for each spot (without interpolation). We limited the maximum number of spots to which one cell can be charted to 10 as well as the maximum number of cells that can be mapped to a spot were limited to 8. The single cell to spot mappings were then used to calculate cell type enrichment ( Fig. 3A,G , Fig. 4 ). Derivation of consensus spatial tumor cell states Spatially variable tumor cell clusters were identified using Celltrek output. Each cluster had a unique set of differentially expressed genes (DEGs). For each sample, the DEGs for each tumor cluster were combined in a binary matrix of genes across samples. This matrix was clustered using a consensus clustering algorithm called cola 89 (version 2.7.1). The ATC method was used for consensus partitioning and spherical k-means was used for sub-group classification. This approach yielded 5 major groups of DEGs that were further used to define the consensus spatial tumor cell states. To further standardize the number of genes in resulting tumor cell states, genes not present in at least 3 patients overall and at least 1 patient in both techniques (OCT & FFPE) were removed from the final set. The resulting gene set was ranked using cumulative mutual information, and the top 100 genes for each program were used to further classify tumor cells in the single cell reference dataset. Ligand receptor interaction analysis in scRNA-seq data CellChat 90 was used to define ligand receptor interactions to better understand how tumor cells (and cell states) communicate. CellChat defines the directionality of the interaction by calculating a communication probability based on a mass-distance algorithm. We identified interactions between various cell-type pairs, namely myeloid-tumor and tumor-myeloid, based on the direction of signal between the two cell types. Batch correction and integration of Visium data using Reciprocal Principal Component analysis Visium data from all samples was combined and batch-corrected (treating each patient as a “batch”) using Seurat’s 80 reciprocal principal component analysis (RPCA) method. In summary, unique variable features were identified across samples using the ‘SCT’ assay, and the merged dataset was normalized with ‘SCTransform.’ For each patient, ‘SelectIntegrationFeatures’ was applied with ‘nfeatures = 2000’. RPCA-based batch correction was then performed using FindIntegrationAnchors, leveraging the previously selected variable features. The data was integrated using ‘IntegrateData’, with k-nearest neighbors (‘k.weight’) set to 50. Statistically significant Principal Components (nPC) were selected using the Jackstraw test. The “integrated” assay of the merged data was then scaled and ‘RunPCA’ was used with nPC=63. Clustering was then performed on the resulting object using multiple resolutions. Cluster resolution 0.9 was used for further downstream analysis as it was able to capture spatial diversity as indicated in pathologist annotations (Fig. S2). Cluster Adjacency We defined a metric, Cluster Adjacency (C A ), that reflects the probability that spots from two clusters, i and j, were “adjacent” given the frequency of each cluster and the total number of adjacent spot pairs. Here, two spots were considered adjacent if the Euclidian distance between them was ∼100μm (which is the minimum distance between 2 spots on a Visium slide). Below, P(i)P(j) is the probability of drawing one spot from cluster i and one spot from cluster j given the number of spots in each cluster; P(i, j) is the observed frequency at which spots from i and j are adjacent. Then, Where and (similar for P(j) ) To summarize the CA score across samples, we defined the following: weighted mean , where w s = adjacent spot pairs for i, j in sample s , and Σ s w s = adjacent spot pairs for i, j across samples The weighted mean adjacency score was calculated for each ( i, j ) combination across all RPCA cluster pairs. The cluster adjacency network ( Fig. 3G ) was then drawn using Cytoscape, where circos nodes represent the cell type distribution (from Celltrek) for each cluster. The edge thickness indicates the Weighted mean adjacency( i, j ) between a cluster pair and the edge color/transparency indicate proportion of samples in which the cluster pair was present with minimum sample proportion set to 20%. Spatially adjacent spots for selected cluster pairs are shown in Fig. S6. Ligand receptor interaction analysis in Visium data To analyze the spatial variation in ligand-receptor (LR) interactions of Visium data at per-spot resolution, stLearn 91 was utilized with ligand-receptor pairs obtained from the Cellchat V2 92 database. For contact-dependent interactions, calculations were performed within spot mode , while for "Secreted Signaling" and "Extracellular Matrix Remodeling" interactions (diffusion-dependent interactions), statistically significant spot pairs were identified over a given distance. It has been shown that different LR pairs interact across varying distance ranges 93 . To identify the optimal set of significant spot pairs, multiple distance parameters were explored, and a distance of 2000 μm was selected. Statistically significant contact-dependent LR pairs for each spot were identified using stLearn, and their enrichment within each cluster was assessed using the hypergeometric distribution. For contact dependent interactions if they were significant for at least 10 spots in a cluster and 100 spots over all were used in further analysis.Hypergeometric test was done for calculating significance of enrichment for each LR pair in each cluster, where: P-values were corrected using the Bonferroni correction. The top five contact-dependent LR pairs for each cluster are shown in Fig. 4A , while all significant LR pairs meeting above threshold are provided in the Supplementary Data. For LR analysis over a distance, stLearn gives us significant spots for each LR pair. In order to distinguish a sender spot with Ligand L k from a receiver spot with Receptor R k (where k indicates known interaction based on cellchat V2 database) for significant spots over a distance, Spots with L k expression ≥ 2 were classified as sender spot/s similarly, spot with expression R k ≥ 2 were designated as receiver spot/s. For all directionally unique spot pairs within/between any pair of batch corrected clusters we calculated LR connectivity score. Where U = Unique sender and receiver spot pairs significant for at least one LR pair between Ci and Cj T = Total possible spot pairs involving Ci (sender) and Cj (receiver) LR connectivity is shown in Fig. 4B with thickness indicating the connectivity score and color indicating no of samples where the cluster pairs interact. To further understand enrichment of certain LR pairs in inter-niche (inter-Cluster) and intra-niche (intra-Cluster) LR signalling between a cluster-pair or between 2 groups (with groups being edge/interface/non-specific/core) we conducted hypergeometric test of enrichment for each LR pair (L k R k ) such that Resulting p -values for each analysis were corrected using the Bonferronni correction. Significant LRs and their enrichment scores meeting the threshold criteria and present in at least 3 samples and/or 2 patients are shown ( Fig. 4 , Supplementary Data). Hexagon analysis Hexagons were superimposed on Celltrek output and the pseudo-bulk expression profile of each hexagon was calculated. The resulting “hexagon-by-gene” expression matrix was analyzed using a standard Seurat pipeline, with ‘vst’ normalization and Louvain clustering. The resulting clusters were classified into niches based on the predominant broad celltype (i.e. Tumor/Neuron/Glia/Vascular/Myeloid/Lymphoid) recapitulating the niches identified through integrated spot-based clustering and Celltrek mapping. Xenium panel design We designed a panel of oligonucleotides targeting transcripts from 362 genes, comprising 266 genes from the predesigned 10X Genomics human brain-and-glioblastoma panel and 96 custom genes indicated by our scRNA-seq analyses (Table S4). Together, this customized panel captures diverse neurons, glia, endothelial cells, microglia, immune cells, GBM cells (including both spatial GBM cell states derived from our Visium data and published scRNA-seq-derived GBM cell states), and ligand-receptor pairs that putatively mediate tumor/neuron, tumor/T cell, and tumor/macrophage interactions. Xenium In Situ Workflow Sample Preparation - Primary Cohort Formalin-fixed, paraffin-embedded (FFPE) tumor samples were sectioned into 5-um sections using a microtome onto Xenium slides according to Protocol CG000578 published by 10X Genomics. The Xenium slides were then deparaffinized and decrosslinked according to the Protocol CG000580. Both a pre-designed gene expression panel (Human Brain Probe set) and a custom probe set were used during probe hybridization, ligation, and amplification according to Protocol CG000749, which also contained cell segmentation staining. The Xenium slide was then loaded into the Xenium Analyzer using Protocol CG000584 and regions of interest (ROIs) were initially selected using the on-instrument software. Sample Preparation - Secondary Cohort Samples were processed similarly as previously published 20 . In brief, formalin-fixed, paraffin-embedded (FFPE) tumor samples were sectioned into 5-um sections using a microtome onto Xenium slides according to Protocol CG000578 published by 10X Genomics. The Xenium slides were then deparaffinized and decrosslinked according to the Protocol CG000580. Both a pre-designed gene expression panel (Human Brain Probe set) and a custom probe set were used during probe hybridization, ligation, and amplification according to Protocol CG000582. The Xenium slide was then loaded into the Xenium Analyzer using Protocol CG000584 and regions of interest (ROIs) were initially selected using the on-instrument software. Sample Filtering and Cell Type Identification - Primary Cohort Raw FASTQ files were analyzed in Xenium Ranger (version 2.0.0.10) using the Xenium Multi-Tissue Stain cell segmentation method. Data were then processed using the Seurat R package (5.0.1). Initial filtering, and subsequent analysis, was performed on a per-sample basis. Counts were then log-normalized and PCA was performed with the optimal number of PCs determined based upon elbow plots, jackstraw resampling, and PC expression heatmaps. Dimensionality reduction and visualization were performed with the UMAP algorithm (Seurat implementation) followed by unsupervised graph-based clustering. Clusters with overall low nCount and nFeature were removed, and PCA and UMAP visualization was performed again on the filtered object. Cell types in each sample were individually characterized using marker genes and differentially expressed genes (Table S6-S7). Samples were then merged using SCTransform (Seurat implementation) and PCA and UMAP visualization was performed again on the merged object. Cell types were then standardized by performing hierarchical clustering on different groups of cell types, such as immune cells, tumor cells, etc. Cell types with consistent gene expression profiles were labeled similarly and differentially expressed genes of each cell type were analyzed to confirm correct identification. Visium-based tumor signature scores were generated using AddModuleScore (Seurat implementation) and applied to patient-specific clusters within the primary cohort to identify tumor cell identities. All patient-specific clusters were subsequently hierarchically clustered according to Visium-based tumor signature scores and appropriately classified. Sample Filtering and Cell Type Identification - Secondary Cohort Raw FASTQ files were reprocessed in Xenium Ranger (version 1.7.0.2) with an expansion distance of 0um. Data were then processed using the Seurat R package (5.0.1). Initial filtering excluded cells that had any of the following: an nFeature count less than 5, an nFeature count greater than 125, or an nCount value greater than 400. Samples were then merged and down-sampled to 50,000 cells using SketchData (Seurat). Counts were then log-normalized and PCA was performed with the optimal number of PCs determined based upon elbow plots, jackstraw resampling, and PC expression heatmaps. Dimensionality and visualization were performed with the UMAP algorithm (Seurat implementation) followed by unsupervised graph-based clustering. Clusters with overall low nCount and nFeature were removed, and PCA and UMAP visualization was performed again on the filtered object. To label cells, the primary cohort Seurat object was first log-normalized and PCA was applied. Then, FindTransferAnchors and TransferData was performed using the primary cohort Seurat object as a reference and applied to the secondary cohort Seurat object. Similarity between gene expression profiles of primary and secondary cohort cell types was assessed by performing a correlation analysis. Xenium Niche Analysis Niche analysis was performed using the BANKSY R package 40 . On a per sample basis, RunBanksy was performed on log-normalized data with a lambda value of 1.0 for both primary and secondary cohort samples. Relative proportions of each cell type per sample-specific niche were then aggregated and log-normalized. Hierarchical clustering was then performed and niches were appropriately classified based upon cell type composition. Statistical significance for cell type enrichment per niche was calculated using the hypergeometric test. Xenium Enrichment Score Enrichment score was calculated by summing the relative frequency of a target cell state within a specific radius (or radial range) of a source cell for all source cells, scaled to the overall frequency of the particular cell state. Xenium Niche Distance Analysis Average minimum distances between niches were quantified by first calculating the distance between a cell of a given primary niche (Niche A) to its nearest neighbor in a different, secondary niche (Niche B) within a given sample. This distance was labeled as the “minimum distance” and calculated for every cell in a primary niche (Niche A), on a per-sample basis, and averaged to calculate the minimum distance from Niche A to Niche B for a given sample. This was then averaged across samples to give the final average minimum distance. The node-edge map was then generated using layout_with_kk (igraph package 94 ), which generates a layout using the Kamada-Kawai algorithm. Xenium Ligand Receptor Analysis Ligand receptor analysis was performed using the CellChat v2 R package 92 on a per sample basis. Significant ligand-receptor (LR) pairs between every cell type across samples were aggregated. For cell types A and B and LR i , we considered A to interact with B via i if the interaction was observed in more than 3 samples. Then, for each cell type pair, we counted unique i and total i .The total number of ligand-receptor interactions were quantified by summing the number of statistically significant ligand-receptor interactions across samples between two given pairs of cell types. The number of unique ligand-receptor interactions was quantified by summing the number of unique interactions across samples between two given pairs of cell types. Immunofluorescence Analysis Formalin-fixed, paraffin-embedded (FFPE) tumor samples were sectioned into 5-um sections and baked at 55°C for 2 hours. Slides were then deparaffinized in: xylene for 10 minutes (twice), 100% ethanol for 5 minutes (twice), 90% ethanol for 3 minutes (once), 70% ethanol for 3 minutes (once), and PBS for 5 minutes (twice). Slides were then incubated in 1X Tris-EDTA buffer, pH 9 (Abcam) at 95C for 30 minutes and then cooled in antigen retrieval buffer for 25 minutes. Slides were then washed in PBS for 5 minutes and subsequently incubated in blocking buffer (1X PBS with 0.1% Triton X-100, 5% BSA, and 5% donkey serum). Primary antibodies (diluted in blocking buffer) were used as follows: goat anti-LYVE1 antibody 1:30 (AF2089, RND Systems), rabbit anti-CD163 1:500 (93498, Cell Signaling Technology). Tissues were stained with primary antibody solution overnight at 4°C and then washed in PBS for 10 minutes thrice. Secondary antibodies (diluted in blocking buffer) were used as follows: AlexaFluor 488 donkey anti-goat 1:500 (Jackson ImmunoResearch), AlexaFluor 647 donkey anti-rabbit 1:500 (Jackson ImmunoResearch). Slides were incubated in secondary antibody solution at room temperature for 1 hour and then washed in PBS for 10 minutes thrice. Slides were then stained with DAPI at a final concentration of 5 μg/mL for 30 minutes and then mounted with VectaShield Vibrance mounting medium (VectorLabs). Images were acquired by fluorescence microscopy with a Nikon AXR microscope. DATA AVAILABILITY Visium and scRNA-seq FASTQ files will be made available in the NCBI Sequence Read Archive (SRA) (PRJEXXXX). Raw Xenium data will be made available on Gene expression Omnibus (GEO) under accession code GSEXXXX.Processed Xenium and Visium Seurat objects and intermediate files used for all analyses are available on the open-access data sharing platform Zenodo at https://zenodo.org/uploads/15014838 . CODE AVAILABILITY Software used for analysis is public and described in detail in the Methods section. All the code used in the analysis is available at https://github.com/Petti-Lab/GBM_Spatial_paper_figures_code/ . Citable code for this study is available in the Zenodo repository https://zenodo.org/uploads/15014838 . Competing Interests Authors declared no Competing interests. Supplementary Figures and Table Legends Fig. S1: H&E visium slides of 22 samples. H&E Visium slides from the 22 samples are presented for both OCT and FFPE specimens. In some cases, distinct core and edge regions were available, as illustrated; in others, transition regions—annotated by the surgeons or pathologists—were obtained, as shown. Fig. S2: Pathologist annotations of the Visium slides. Annotations by the resident pathologist for the H&E slides of ST073021.04 (Core), ST073021.02 (Edge), WU1220.01 (Edge), WU1220.02 (Edge), and WU1227 (Core), with Visium spots superimposed, are presented. Key histopathological features characteristic of GBM such as microvascular proliferation (MVP), palisading necrosis, myelin-rich regions, and tumor-rich regions are illustrated. Fig. S3. Batch-corrected RPCA clustering across patients at resolution 0.1. (A) Visium spots colored by cluster assignment across patients and samples. (B) Bar plots depicting the distribution of patients and samples within each cluster. (C) Gene Ontology (GO) biological processes enriched among the top 200 differentially expressed genes for each cluster (adjusted p ≤ 0.05, ranked by average log fold change). The analysis indicates that low-resolution batch integration effectively identifies biologically distinct clusters, including regions associated with microvascular proliferation (MVP), as well as glia-rich and neuron-rich areas. Fig. S4. Batch-corrected RPCA clustering across patients at resolution 0.9. (A) Visium spots colored by cluster assignment at resolution 0.9 across patients. (B) Bar plots showing the distribution of patients and samples within each cluster. (C) Gene Ontology (GO) (biological processes) enriched among the top 200 differentially expressed genes for each cluster (adjusted p ≤ 0.05, ranked by average log fold change). Higher-resolution (0.9) batch-integration clustering captured more extensive transcriptional variation, revealing both shared and sample-specific clusters, as well as features specifically enriched in either the tumor core or edge. Fig. S5. Visium analysis workflow. (A) Single-cell RNA-seq data was used as a reference for each patient, and cells were mapped onto the corresponding Visium slides using CellTrek (with interpolation) (1) . Spatially variable tumor cell clusters were identified for each patient using CellTrek. A binary matrix representing these tumor sub-clusters across patients was constructed (2,3) . Consensus partitioning and ranking methods were then applied to this matrix to define spatially variable tumor cell states (4) (see Methods.). (B) Visium data was integrated across patients using RPCA-based batch clustering at a resolution of 0.9. CellTrek (without interpolation) was then used to assign cell types to each Visium spot (1) . Cell type frequencies and cell type enrichments were calculated for each cluster, as shown in Fig. 3A (2) . Subsequently, spatial adjacency analysis of cellular niches (3) and ligand–receptor interaction analysis (4) were performed on the Visium spots, and cell type frequencies were visualized. Fig. S6. Spatial adjacency between transcriptional clusters. As shown in Fig. 3G , multiple clusters display spatial co-localization. Spatial adjacency between the Hypoxia cluster (C3) and the microvascular proliferation (MVP) cluster (C6) is demonstrated on the Visium spots across multiple samples, consistent with the adjacency map in Fig. 3G . Additional examples of spatially adjacent clusters include MVP (C6) with S3 Mesenchymal (C17), MVP (C6) with S3 Mesenchymal (C9), and MVP (C6) with S5 Stem-like (C7), emphasizing the central role of microvascular proliferation in GBM tumor architecture. In edge-enriched regions, spatial adjacency between Oligo-rich (C0) and Neuron-rich (C13) clusters is also observed. Fig. S7: Analysis of Cohort 1 and Cohort 2 samples. (A) UMAP visualization of cohort 2 samples alone. (B) Correlation analysis of gene expression of respective cell type annotations in cohort 1 and cohort 2. (C) Dot plot of expression of select genes used for cell type annotation. (D) Bar plot of frequency of each cell type per sample in both cohorts 1 and 2. (E) Heatmap of niche occurrence within every sample with the total occurrence of each niche across all samples summarized at the top. Fig. S8: Cell chat analysis of Xenium samples. (A) Bar plot showing the number of ligand-receptor interactions with the specified cell type either as the sender or receiver. (B) Interaction network showing ligand receptor interactions with tumor cells as senders and all other cell types as receivers. The density and size of arrows were proportional to the number of interactions detected as significant for a given sender-receiver pair. (C) UMAP visualization of scRNA-seq data including both GBM samples and normal brain cell types. (D) UMAP visualization of both ADM and CALCRL gene expression within scRNA-seq data set. (E) Interaction network of niche sender-receiver pairs of the ADM:CALCRL ligand-receptor pair. The color and density of arrows were proportional to the number of samples the ADM:CALCRL interaction was detected as significant for a given sender-receiver pair. Fig. S9: Analysis of T cell and myeloid cell co-localization within Xenium samples. (A) Distance to the nearest target myeloid cell of each CD4 +. GZMK + T cell separated by localization within a vascular or non-vascular niche averaged on a per sample basis. Each dot represents an individual sample. (B) Distance to the nearest target myeloid cell of each CD4 +. Memory T cell separated by localization within a vascular or non-vascular niche averaged on a per sample basis. Each dot represents an individual sample. (C) Frequency of lymphoid and myeloid cell types based upon vascular or non-vascular localization calculated on a per sample basis. (D) Correlation of GZMB and LYVE1 expression based upon localization in vascular or non-vascular niche. Data derived from IVY-Gap. (E) UMAP visualization of LYVE1 and CD163 expression within scRNA-seq data set. Significance in (A-C) was calculated using the Wilcoxon matched-pairs signed rank test. Table S1: Clinical details of the Samples sequenced for Visium and Xenium Table S2: QC metrics of Visium data Table S3: Gene sets of tumor states (S1-S5) identified after Consensus Clustering. Table S4: Xenium Probe sets. Table S5: Xenium QC metrics. Table S6: Gene markers for cell type identification. Table S7: DEGs for each cell type (both CD45+ and CD45-cells). Table S8: Sample-specific BANKSY clusters associated with Fig. 5C . Table S9: Individual distances used to calculate node-edge plot in Fig. 5G . Table S10: All individual significant LR interactions per cell type pair per sample. Supplementary Data (Contact_dependent_LR) : Within-spot Hypergeomtric test enrichment results of LR (per cluster) for Contact-dependent ligand-receptor pairs only. Supplementary Data (diffusion_LR_enrichment_byDomain) : Hypergeomtric test results of LR enrichment for diffusion-based ligand-receptor pairs. Supplementary Data (DGIDB_visium) : Ligand receptor candidates from Visium and their associated drug interactions from DGIDB database. Supplementary Data (DGIDB_xenium) : Ligand receptor candidates from Xenium and their associated drug interactions from DGIDB database. ACKNOWLEDGEMENTS This work was supported by a Cancer Research Foundation Young Investigator’s Award (to A.A.P.), National Institutes of Health grants R01 NS094670, R01 NS128470 (to A.H.K.), the Christopher Davidson and Knight Family Fund (to A.H.K.); and the Duesenberg Research Fund (to A.H.K.). Footnotes Some signs and symbols were missing. https://zenodo.org/uploads/15014838 REFERENCES 1. ↵ Weller , M. , Butowski , N. , Tran , D.D. , Recht , L.D. , Lim , M. , Hirte , H. , Ashby , L. , Mechtler , L. , Goldlust , S.A. , Iwamoto , F. , et al. ( 2017 ). Rindopepimut with temozolomide for patients with newly diagnosed, EGFRvIII-expressing glioblastoma (ACT IV): a randomised, double-blind, international phase 3 trial . Lancet Oncol . 18 , 1373 – 1385 . OpenUrl CrossRef PubMed 2. ↵ Liau , L.M. , Ashkan , K. , Brem , S. , Campian , J.L. , Trusheim , J.E. , Iwamoto , F.M. , Tran , D.D. , Ansstas , G. , Cobbs , C.S. , Heth , J.A. , et al. ( 2023 ). Association of autologous tumor lysate-loaded dendritic cell vaccination with extension of survival among patients with newly diagnosed and recurrent glioblastoma: A phase 3 prospective externally controlled cohort trial: A phase 3 prospective externally controlled cohort trial . JAMA Oncol . 9 , 112 – 121 . OpenUrl CrossRef PubMed 3. ↵ Omuro , A. , Brandes , A.A. , Carpentier , A.F. , Idbaih , A. , Reardon , D.A. , Cloughesy , T. , Sumrall , A. , Baehring , J. , van den Bent , M. , Bähr , O. , et al. ( 2023 ). Radiotherapy combined with nivolumab or temozolomide for newly diagnosed glioblastoma with unmethylated MGMT promoter: An international randomized phase III trial . Neuro. Oncol . 25 , 123 – 134 . OpenUrl CrossRef PubMed 4. ↵ Cancer Genome Atlas Research Network ( 2008 ). Comprehensive genomic characterization defines human glioblastoma genes and core pathways . Nature 455 , 1061 – 1068 . OpenUrl CrossRef PubMed Web of Science 5. Brennan , C.W. , Verhaak , R.G.W. , McKenna , A. , Campos , B. , Noushmehr , H. , Salama , S.R. , Zheng , S. , Chakravarty , D. , Sanborn , J.Z. , Berman , S.H. , et al. ( 2013 ). The somatic genomic landscape of glioblastoma . Cell 155 , 462 – 477 . OpenUrl CrossRef PubMed Web of Science 6. ↵ Schaettler , M.O. , Richters , M.M. , Wang , A.Z. , Skidmore , Z.L. , Fisk , B. , Miller , K.E. , Vickery , T.L. , Kim , A.H. , Chicoine , M.R. , Osbun , J.W. , et al. ( 2022 ). Characterization of the Genomic and Immunologic Diversity of Malignant Brain Tumors through Multisector Analysis . Cancer Discov . 12 , 154 – 171 . OpenUrl Abstract / FREE Full Text 7. ↵ Verhaak , R.G.W. , Hoadley , K.A. , Purdom , E. , Wang , V. , Qi , Y. , Wilkerson , M.D. , Miller , C.R. , Ding , L. , Golub , T. , Mesirov , J.P. , et al. ( 2010 ). Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1 . Cancer Cell 17 , 98 – 110 . OpenUrl CrossRef PubMed Web of Science 8. ↵ Wang , Q. , Hu , B. , Hu , X. , Kim , H. , Squatrito , M. , Scarpace , L. , Decarvalho , A.C. , Lyu , S. , Li , P. , Li , Y. , et al. ( 2017 ). Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment . Cancer Cell 32 , 42 – 56.e6 . OpenUrl CrossRef PubMed 9. ↵ Neftel , C. , Laffy , J. , Filbin , M.G. , Hara , T. , Shore , M.E. , Rahme , G.J. , Richman , A.R. , Silverbush , D. , Shaw , M.L. , Hebert , C.M. , et al. ( 2019 ). An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma . Cell 178 , 835 – 849.e21 . OpenUrl CrossRef PubMed 10. Patel , A.P. , Tirosh , I. , Trombetta , J.J. , Shalek , A.K. , Gillespie , S.M. , Wakimoto , H. , Cahill , D.P. , Nahed , B.V. , Curry , W.T. , Martuza , R.L. , et al. ( 2014 ). Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma . Science 344 , 1396 – 1401 . OpenUrl Abstract / FREE Full Text 11. Bhaduri , A. , Di Lullo , E. , Jung , D. , Müller , S. , Crouch , E.E. , Espinosa , C.S. , Ozawa , T. , Alvarado , B. , Spatazza , J. , Cadwell , C.R. , et al. ( 2020 ). Outer Radial Glia-like Cancer Stem Cells Contribute to Heterogeneity of Glioblastoma . Cell Stem Cell 26 , 48 – 63 e6 . OpenUrl CrossRef PubMed 12. ↵ Couturier , C.P. , Ayyadhury , S. , Le , P.U. , Nadaf , J. , Monlong , J. , Riva , G. , Allache , R. , Baig , S. , Yan , X. , Bourgey , M. , et al. ( 2020 ). Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy . Nat. Commun . 11 , 3406 . OpenUrl CrossRef PubMed 13. Wang , L. , Babikir , H. , Müller , S. , Yagnik , G. , Shamardani , K. , Catalan , F. , Kohanbash , G. , Alvarado , B. , Di Lullo , E. , Kriegstein , A. , et al. ( 2019 ). The phenotypes of proliferating glioblastoma cells reside on a single axis of variation . Cancer Discov . 9 , 1708 – 1719 . OpenUrl Abstract / FREE Full Text 14. ↵ Garofano , L. , Migliozzi , S. , Oh , Y.T. , D’Angelo , F. , Najac , R.D. , Ko , A. , Frangaj , B. , Caruso , F.P. , Yu , K. , Yuan , J. , et al. ( 2021 ). Pathway-based classification of glioblastoma uncovers a mitochondrial subtype with therapeutic vulnerabilities. Nat . Cancer 2 , 141 – 156 . OpenUrl 15. ↵ Darmanis , S. , Sloan , S.A. , Croote , D. , Mignardi , M. , Chernikova , S. , Samghababi , P. , Zhang , Y. , Neff , N. , Kowarsky , M. , Caneda , C. , et al. ( 2017 ). Single-Cell RNA-Seq Analysis of Infiltrating Neoplastic Cells at the Migrating Front of Human Glioblastoma . Cell Rep . 21 , 1399 – 1410 . OpenUrl CrossRef PubMed 16. ↵ Brooks , W.H. , Netsky , M.G. , Normansell , D.E. , and Horwitz , D.A . ( 1972 ). Depressed cell-mediated immunity in patients with primary intracranial tumors. Characterization of a humoral immunosuppressive factor . J. Exp. Med . 136 , 1631 – 1647 . OpenUrl Abstract 17. Chongsathidkiet , P. , Jackson , C. , Koyama , S. , Loebel , F. , Cui , X. , Farber , S.H. , Woroniecka , K. , Elsamadicy , A.A. , Dechant , C.A. , Kemeny , H.R. , et al. ( 2018 ). Sequestration of T cells in bone marrow in the setting of glioblastoma and other intracranial tumors . Nat. Med . 24 , 1459 – 1468 . OpenUrl CrossRef PubMed 18. Mathewson , N.D. , Ashenberg , O. , Tirosh , I. , Gritsch , S. , Perez , E.M. , Marx , S. , Jerby-Arnon , L. , Chanoch-Myers , R. , Hara , T. , Richman , A.R. , et al. ( 2021 ). Inhibitory CD161 receptor identified in glioma-infiltrating T cells by single-cell analysis . Cell 184 , 1281 – 1298.e26 . OpenUrl CrossRef PubMed 19. Naulaerts , S. , Datsi , A. , Borras , D.M. , Antoranz Martinez , A. , Messiaen , J. , Vanmeerbeek , I. , Sprooten , J. , Laureano , R.S. , Govaerts , J. , Panovska , D. , et al. ( 2023 ). Multiomics and spatial mapping characterizes human CD8 + T cell states in cancer . Sci. Transl. Med . 15 , eadd1016 . OpenUrl CrossRef PubMed 20. ↵ Wang , A.Z. , Mashimo , B.L. , Schaettler , M.O. , Sherpa , N.D. , Leavitt , L.A. , Livingstone , A.J. , Khan , S.M. , Li , M. , Anzaldua-Campos , M.I. , Bradley , J.D. , et al. ( 2024 ). Glioblastoma-Infiltrating CD8+ T Cells Are Predominantly a Clonally Expanded GZMK+ Effector Population . Cancer Discov . 14 , 1106 – 1131 . OpenUrl CrossRef PubMed 21. ↵ Xuan , W. , Lesniak , M.S. , James , C.D. , Heimberger , A.B. , and Chen , P . ( 2021 ). Context-Dependent Glioblastoma-Macrophage/Microglia Symbiosis and Associated Mechanisms . Trends Immunol . 42 , 280 – 292 . OpenUrl CrossRef PubMed 22. ↵ Ravi , V.M. , Neidert , N. , Will , P. , Joseph , K. , Maier , J.P. , Kückelhaus , J. , Vollmer , L. , Goeldner , J.M. , Behringer , S.P. , Scherer , F. , et al. ( 2022 ). T-cell dysfunction in the glioblastoma microenvironment is mediated by myeloid cells releasing interleukin-10 . Nat. Commun . 13 , 925 . OpenUrl CrossRef PubMed 23. ↵ Greenwald , A.C. , Darnell , N.G. , Hoefflin , R. , Simkin , D. , Mount , C.W. , Gonzalez Castro , L.N. , Harnik , Y. , Dumont , S. , Hirsch , D. , Nomura , M. , et al. ( 2024 ). Integrative spatial analysis reveals a multi-layered organization of glioblastoma . Cell 187 , 2485 – 2501.e26 . OpenUrl CrossRef PubMed 24. ↵ Ravi , V.M. , Will , P. , Kueckelhaus , J. , Sun , N. , Joseph , K. , Salié , H. , Vollmer , L. , Kuliesiute , U. , von Ehr , J. , Benotmane , J.K. , et al. ( 2022 ). Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma . Cancer Cell 40 , 639 – 655.e13 . OpenUrl CrossRef PubMed 25. ↵ Jackson , C. , Cherry , C. , Bom , S. , Dykema , A.G. , Wang , R. , Thompson , E. , Zhang , M. , Li , R. , Ji , Z. , Hou , W. , et al. ( 2025 ). Distinct myeloid-derived suppressor cell populations in human glioblastoma . Science 387 , eabm5214 . OpenUrl CrossRef PubMed 26. ↵ Wang , W. , Li , T. , Cheng , Y. , Li , F. , Qi , S. , Mao , M. , Wu , J. , Liu , Q. , Zhang , X. , Li , X. , et al. ( 2024 ). Identification of hypoxic macrophages in glioblastoma with therapeutic potential for vasculature normalization . Cancer Cell 42 , 815 – 832.e12 . OpenUrl CrossRef PubMed 27. ↵ Haley , M.J. , Bere , L. , Minshull , J. , Georgaka , S. , Garcia-Martin , N. , Howell , G. , Coope , D.J. , Roncaroli , F. , King , A. , Wedge , D.C. , et al. ( 2024 ). Hypoxia coordinates the spatial landscape of myeloid cells within glioblastoma to affect survival . Sci Adv 10 , eadj3301 . OpenUrl CrossRef PubMed 28. ↵ Harwood , D.S.L. , Pedersen , V. , Bager , N.S. , Schmidt , A.Y. , Stannius , T.O. , Areškevičiūtė , A. , Josefsen , K. , Nørøxe , D.S. , Scheie , D. , Rostalski , H. , et al. ( 2024 ). Glioblastoma cells increase expression of notch signaling and synaptic genes within infiltrated brain tissue . Nat. Commun . 15 , 7857 . OpenUrl CrossRef PubMed 29. ↵ De Bonis , P. , Anile , C. , Pompucci , A. , Fiorentino , A. , Balducci , M. , Chiesa , S. , Lauriola , L. , Maira , G. , and Mangiola , A. ( 2013 ). The influence of surgery on recurrence pattern of glioblastoma . Clin. Neurol. Neurosurg . 115 , 37 – 43 . OpenUrl CrossRef PubMed 30. ↵ Siletti , K. , Hodge , R. , Mossi Albiach , A. , Lee , K.W. , Ding , S.-L. , Hu , L. , Lönnerberg , P. , Bakken , T. , Casper , T. , Clark , M. , et al. ( 2023 ). Transcriptomic diversity of cell types across the adult human brain . Science 382 , eadd7046 . OpenUrl CrossRef PubMed 31. ↵ Hadjipanayis , C.G. , Widhalm , G. , and Stummer , W . ( 2015 ). What is the surgical benefit of utilizing 5-aminolevulinic acid for fluorescence-guided surgery of malignant gliomas? Neurosurgery 77 , 663 – 673 . OpenUrl CrossRef PubMed 32. ↵ Wei , R. , He , S. , Bai , S. , Sei , E. , Hu , M. , Thompson , A. , Chen , K. , Krishnamurthy , S. , and Navin , N.E . ( 2022 ). Spatial charting of single-cell transcriptomes in tissues . Nat. Biotechnol . 40 , 1190 – 1199 . OpenUrl CrossRef PubMed 33. ↵ Mahlokozera , T. , Taiwo , R. , Patel , B. , Hafez , D. , Paturu , M. , Salehi , A. , Mao , D.D. , Gujar , A. , Dunn , G.P. , Mosammaparast , N. , et al. ( 2021 ). Competitive binding of E3 ligases TRIM26 and WWP2 controls SOX2 proteostasis in glioblastoma . Under Revision . 34. ↵ Wang , L. , Jung , J. , Babikir , H. , Shamardani , K. , Jain , S. , Feng , X. , Gupta , N. , Rosi , S. , Chang , S. , Raleigh , D. , et al. ( 2022 ). A single-cell atlas of glioblastoma evolution under therapy reveals cell-intrinsic and cell-extrinsic therapeutic targets. Nat . Cancer 3 , 1534 – 1552 . OpenUrl 35. ↵ Zheng , Y. , Carrillo-Perez , F. , Pizurica , M. , Heiland , D.H. , and Gevaert , O . ( 2023 ). Spatial cellular architecture predicts prognosis in glioblastoma . Nat. Commun . 14 , 4122 . OpenUrl CrossRef PubMed 36. ↵ Moffet , J.J.D. , Fatunla , O.E. , Freytag , L. , Kriel , J. , Jones , J.J. , Roberts-Thomson , S.J. , Pavenko , A. , Scoville , D.K. , Zhang , L. , Liang , Y. , et al. ( 2023 ). Spatial architecture of high-grade glioma reveals tumor heterogeneity within distinct domains . Neurooncol. Adv . 5 , vdad142 . OpenUrl 37. ↵ Venkatesh , H.S. , Morishita , W. , Geraghty , A.C. , Silverbush , D. , Gillespie , S.M. , Arzt , M. , Tam , L.T. , Espenel , C. , Ponnuswami , A. , Ni , L. , et al. ( 2019 ). Electrical and synaptic integration of glioma into neural circuits . Nature 573 , 539 – 545 . OpenUrl CrossRef PubMed 38. ↵ Jossin , Y . ( 2020 ). Reelin functions, mechanisms of action and signaling pathways during brain development and maturation . Biomolecules 10 , 964 . OpenUrl CrossRef PubMed 39. ↵ Zhang , X. , Shao , S. , and Li , L . ( 2020 ). Characterization of class-3 semaphorin receptors, neuropilins and plexins, as therapeutic targets in a pan-cancer study . Cancers (Basel ) 12 , 1816 . OpenUrl PubMed 40. ↵ Singhal , V. , Chou , N. , Lee , J. , Yue , Y. , Liu , J. , Chock , W.K. , Lin , L. , Chang , Y.-C. , Teo , E.M.L. , Aow , J. , et al. ( 2024 ). BANKSY unifies cell typing and tissue domain segmentation for scalable spatial omics data analysis . Nat Genet 56 , 431 – 441 . OpenUrl CrossRef PubMed 41. ↵ Tomaszewski , W. , Sanchez-Perez , L. , Gajewski , T.F. , and Sampson , J.H . ( 2019 ). Brain tumor microenvironment and host state: Implications for immunotherapy . Clin. Cancer Res . 25 , 4202 – 4210 . OpenUrl Abstract / FREE Full Text 42. ↵ Sun , L. , Hui , A.-M. , Su , Q. , Vortmeyer , A. , Kotliarov , Y. , Pastorino , S. , Passaniti , A. , Menon , J. , Walling , J. , Bailey , R. , et al. ( 2006 ). Neuronal and glioma-derived stem cell factor induces angiogenesis within the brain . Cancer Cell 9 , 287 – 300 . OpenUrl CrossRef PubMed Web of Science 43. ↵ Sihto , H. , Tynninen , O. , Bützow , R. , Saarialho-Kere , U. , and Joensuu , H . ( 2007 ). Endothelial cell KIT expression in human tumours . J. Pathol . 211 , 481 – 488 . OpenUrl CrossRef PubMed 44. ↵ Hara , T. , Chanoch-Myers , R. , Mathewson , N.D. , Myskiw , C. , Atta , L. , Bussema , L. , Eichhorn , S.W. , Greenwald , A.C. , Kinker , G.S. , Rodman , C. , et al. ( 2021 ). Interactions between cancer cells and immune cells drive transitions to mesenchymal-like states in glioblastoma . Cancer Cell 39 , 779 – 792.e11 . OpenUrl CrossRef PubMed 45. Huang-Hobbs , E. , Cheng , Y.-T. , Ko , Y. , Luna-Figueroa , E. , Lozzi , B. , Taylor , K.R. , McDonald , M. , He , P. , Chen , H.-C. , Yang , Y. , et al. ( 2023 ). Remote neuronal activity drives glioma progression through SEMA4F . Nature 619 , 844 – 850 . OpenUrl CrossRef PubMed 46. ↵ Venkatesh , H.S. , Johung , T.B. , Caretti , V. , Noll , A. , Tang , Y. , Nagaraja , S. , Gibson , E.M. , Mount , C.W. , Polepalli , J. , Mitra , S.S. , et al. ( 2015 ). Neuronal activity promotes glioma growth through neuroligin-3 secretion . Cell 161 , 803 – 816 . OpenUrl CrossRef PubMed 47. ↵ Metellus , P. , Voutsinos-Porche , B. , Nanni-Metellus , I. , Colin , C. , Fina , F. , Berenguer , C. , Dussault , N. , Boudouresque , F. , Loundou , A. , Intagliata , D. , et al. ( 2011 ). Adrenomedullin expression and regulation in human glioblastoma, cultured human glioblastoma cell lines and pilocytic astrocytoma . Eur. J. Cancer 47 , 1727 – 1735 . OpenUrl CrossRef PubMed 48. ↵ Garayoa , M. , Martínez , A. , Lee , S. , Pío , R. , An , W.G. , Neckers , L. , Trepel , J. , Montuenga , L.M. , Ryan , H. , Johnson , R. , et al. ( 2000 ). Hypoxia-inducible factor-1 (HIF-1) up-regulates adrenomedullin expression in human tumor cell lines during oxygen deprivation: a possible promotion mechanism of carcinogenesis . Mol. Endocrinol . 14 , 848 – 862 . OpenUrl CrossRef PubMed Web of Science 49. ↵ Puchalski , R.B. , Shah , N. , Miller , J. , Dalley , R. , Nomura , S.R. , Yoon , J.-G. , Smith , K.A. , Lankerovich , M. , Bertagnolli , D. , Bickley , K. , et al. ( 2018 ). An anatomic transcriptional atlas of human glioblastoma . Science 360 , 660 – 663 . OpenUrl Abstract / FREE Full Text 50. ↵ Schalck , A. , Sakellariou-Thompson , D. , Forget , M.-A. , Sei , E. , Hughes , T.G. , Reuben , A. , Bai , S. , Hu , M. , Kumar , T. , Hurd , M.W. , et al. ( 2022 ). Single-Cell Sequencing Reveals Trajectory of Tumor-Infiltrating Lymphocyte States in Pancreatic Cancer . Cancer Discov . 12 , 2330 – 2349 . OpenUrl CrossRef PubMed 51. ↵ Anstee , J.E. , Feehan , K.T. , Opzoomer , J.W. , Dean , I. , Muller , H.P. , Bahri , M. , Cheung , T.S. , Liakath-Ali , K. , Liu , Z. , Choy , D. , et al. ( 2023 ). LYVE-1+ macrophages form a collaborative CCR5-dependent perivascular niche that influences chemotherapy responses in murine breast cancer . Dev. Cell 58 , 1548 – 1561.e10 . OpenUrl CrossRef PubMed 52. Nalio Ramos , R. , Missolo-Koussou , Y. , Gerber-Ferder , Y. , Bromley , C.P. , Bugatti , M. , Núñez , N.G. , Tosello Boari , J. , Richer , W. , Menger , L. , Denizeau , J. , et al. ( 2022 ). Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer . Cell 185 , 1189 – 1207.e25 . OpenUrl CrossRef PubMed 53. ↵ Elfstrum , A.K. , Rumahorbo , A.H. , Reese , L.E. , Nelson , E.V. , McCluskey , B.M. , and Schwertfeger , K.L . ( 2024 ). LYVE-1-expressing macrophages modulate the hyaluronan-containing extracellular matrix in the mammary stroma and contribute to mammary tumor growth . Cancer Res. Commun . 4 , 1380 – 1397 . OpenUrl CrossRef PubMed 54. ↵ Siret , C. , van Lessen , M. , Bavais , J. , Jeong , H.W. , Reddy Samawar , S.K. , Kapupara , K. , Wang , S. , Simic , M. , de Fabritus , L. , Tchoghandjian , A. , et al. ( 2022 ). Deciphering the heterogeneity of the Lyve1+ perivascular macrophages in the mouse brain . Nat. Commun . 13 , 7366 . OpenUrl CrossRef PubMed 55. Faraco , G. , Park , L. , Anrather , J. , and Iadecola , C . ( 2017 ). Brain perivascular macrophages: characterization and functional roles in health and disease . J. Mol. Med . 95 , 1143 – 1152 . OpenUrl CrossRef PubMed 56. ↵ Wen , W. , Cheng , J. , and Tang , Y . ( 2024 ). Brain perivascular macrophages: current understanding and future prospects . Brain 147 , 39 – 55 . OpenUrl CrossRef PubMed 57. ↵ Alanio , C. , Binder , Z.A. , Chang , R.B. , Nasrallah , M.P. , Delman , D. , Li , J.H. , Tang , O.Y. , Zhang , L.Y. , Zhang , J.V. , Wherry , E.J. , et al. ( 2022 ). Immunologic features in DE Novo and recurrent Glioblastoma are associated with survival outcomes . Cancer Immunol. Res . 10 , 800 – 810 . OpenUrl CrossRef PubMed 58. ↵ Mu , L. , Yang , C. , Gao , Q. , Long , Y. , Ge , H. , DeLeon , G. , Jin , L. , Chang , Y.E. , Sayour , E.J. , Ji , J. , et al. ( 2017 ). CD4+ and perivascular Foxp3+ T cells in glioma correlate with angiogenesis and tumor progression . Front. Immunol . 8 , 1451 . OpenUrl CrossRef PubMed 59. ↵ Miller , T.E. , El Farran , C.A. , Couturier , C.P. , Chen , Z. , D’Antonio , J.P. , Verga , J. , Villanueva , M.A. , Gonzalez Castro , L.N. , Tong , Y.E. , Saadi , T.A. , et al. ( 2025 ). Programs, origins and immunomodulatory functions of myeloid cells in glioma . Nature , 1 – 11 . 60. ↵ Ren , Y. , Huang , Z. , Zhou , L. , Xiao , P. , Song , J. , He , P. , Xie , C. , Zhou , R. , Li , M. , Dong , X. , et al. ( 2023 ). Spatial transcriptomics reveals niche-specific enrichment and vulnerabilities of radial glial stem-like cells in malignant gliomas . Nat. Commun . 14 , 1028 . OpenUrl CrossRef PubMed 61. ↵ White , J. , White , M.P.J. , Wickremesekera , A. , Peng , L. , and Gray , C . ( 2024 ). The tumour microenvironment, treatment resistance and recurrence in glioblastoma . J. Transl. Med . 22 , 540 . OpenUrl CrossRef PubMed 62. Lathia , J.D. , Mack , S.C. , Mulkearns-Hubert , E.E. , Valentim , C.L.L. , and Rich , J.N . ( 2015 ). Cancer stem cells in glioblastoma . Genes Dev . 29 , 1203 – 1217 . OpenUrl Abstract / FREE Full Text 63. ↵ Piper , K. , DePledge , L. , Karsy , M. , and Cobbs , C . ( 2021 ). Glioma stem cells as immunotherapeutic targets: Advancements and challenges . Front. Oncol . 11 , 615704 . OpenUrl CrossRef PubMed 64. ↵ Park , J.H. , and Lee , H.K . ( 2022 ). Current understanding of hypoxia in glioblastoma multiforme and its response to immunotherapy . Cancers (Basel ) 14 , 1176 . OpenUrl PubMed 65. Scharping , N.E. , Menk , A.V. , Whetstone , R.D. , Zeng , X. , and Delgoffe , G.M . ( 2017 ). Efficacy of PD-1 blockade is potentiated by metformin-induced reduction of tumor hypoxia . Cancer Immunol. Res . 5 , 9 – 16 . OpenUrl Abstract / FREE Full Text 66. ↵ Kim , A.-R. , Choi , S.J. , Park , J. , Kwon , M. , Chowdhury , T. , Yu , H.J. , Kim , S. , Kang , H. , Kim , K.-M. , Park , S.-H. , et al. ( 2022 ). Spatial immune heterogeneity of hypoxia-induced exhausted features in high-grade glioma . Oncoimmunology 11 , 2026019 . OpenUrl CrossRef PubMed 67. ↵ Sattiraju , A. , Kang , S. , Giotti , B. , Chen , Z. , Marallano , V.J. , Brusco , C. , Ramakrishnan , A. , Shen , L. , Tsankov , A.M. , Hambardzumyan , D. , et al. ( 2023 ). Hypoxic niches attract and sequester tumor-associated macrophages and cytotoxic T cells and reprogram them for immunosuppression . Immunity 56 , 1825 – 1843.e6 . OpenUrl CrossRef PubMed 68. ↵ Gilbert , M.R. , Dignam , J.J. , Armstrong , T.S. , Wefel , J.S. , Blumenthal , D.T. , Vogelbaum , M.A. , Colman , H. , Chakravarti , A. , Pugh , S. , Won , M. , et al. ( 2014 ). A randomized trial of bevacizumab for newly diagnosed glioblastoma . N. Engl. J. Med . 370 , 699 – 708 . OpenUrl CrossRef PubMed Web of Science 69. ↵ Fu , M. , Zhou , Z. , Huang , X. , Chen , Z. , Zhang , L. , Zhang , J. , Hua , W. , and Mao , Y . ( 2023 ). Use of Bevacizumab in recurrent glioblastoma: a scoping review and evidence map . BMC Cancer 23 , 544 . OpenUrl CrossRef PubMed 70. ↵ Cannon , M. , Stevenson , J. , Stahl , K. , Basu , R. , Coffman , A. , Kiwala , S. , McMichael , J.F. , Kuzma , K. , Morrissey , D. , Cotto , K. , et al. ( 2024 ). DGIdb 5.0: rebuilding the drug-gene interaction database for precision medicine and drug discovery platforms . Nucleic Acids Res . 52 , D1227 – D1235 . OpenUrl CrossRef PubMed 71. ↵ Butler , A. , Hoffman , P. , Smibert , P. , Papalexi , E. , and Satija , R . ( 2018 ). Integrating single-cell transcriptomic data across different conditions, technologies, and species . Nat. Biotechnol . 36 , 411 – 420 . OpenUrl CrossRef PubMed 72. ↵ Tirosh , I. , Izar , B. , Prakadan , S.M. , Wadsworth , M.H. , 2nd . , Treacy , D. , Trombetta , J.J. , Rotem , A. , Rodman , C. , Lian , C. , Murphy , G. , et al. ( 2016 ). Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq . Science 352 , 189 – 196 . OpenUrl Abstract / FREE Full Text 73. ↵ Yu , G. , Wang , L.-G. , Han , Y. , and He , Q.-Y . ( 2012 ). clusterProfiler: an R package for comparing biological themes among gene clusters . OMICS 16 , 284 – 287 . OpenUrl CrossRef PubMed Web of Science 74. ↵ Xu , S. , Hu , E. , Cai , Y. , Xie , Z. , Luo , X. , Zhan , L. , Tang , W. , Wang , Q. , Liu , B. , Wang , R. , et al. ( 2024 ). Using clusterProfiler to characterize multiomics data . Nat. Protoc . 19 , 3292 – 3320 . OpenUrl CrossRef PubMed 75. ↵ Zheng , G.X.Y. , Terry , J.M. , Belgrader , P. , Ryvkin , P. , Bent , Z.W. , Wilson , R. , Ziraldo , S.B. , Wheeler , T.D. , McDermott , G.P. , Zhu , J. , et al. ( 2017 ). Massively parallel digital transcriptional profiling of single cells . Nat Commun 8 , 14049 . OpenUrl CrossRef PubMed 76. ↵ Lun , A.T.L. , McCarthy , D.J. , and Marioni , J.C . ( 2016 ). A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor . F1000Res . 5 , 2122 . OpenUrl 77. ↵ Bais , A.S. , and Kostka , D . ( 2019 ). scds: Computational Annotation of Doublets in Single Cell RNA Sequencing Data . bioRxiv , 564021 . doi: 10.1101/564021 . OpenUrl Abstract / FREE Full Text 78. ↵ Germain , P.-L. , Lun , A. , Garcia Meixide , C. , Macnair , W. , and Robinson , M.D . ( 2021 ). Doublet identification in single-cell sequencing data using scDblFinder . F1000Res . 10 , 979 . OpenUrl 79. ↵ McGinnis , C.S. , Murrow , L.M. , and Gartner , Z.J . ( 2019 ). DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors . Cell Syst 8 , 329 – 337.e4 . OpenUrl CrossRef PubMed 80. ↵ Hao , Y. , Hao , S. , Andersen-Nissen , E. , Mauck , W.M. , 3rd . , Zheng , S. , Butler , A. , Lee , M.J. , Wilk , A.J. , Darby , C. , Zager , M. , et al. ( 2021 ). Integrated analysis of multimodal single-cell data . Cell 184 , 3573 – 3587.e29 . OpenUrl CrossRef PubMed 81. ↵ Aran , D. , Looney , A.P. , Liu , L. , Wu , E. , Fong , V. , Hsu , A. , Chak , S. , Naikawadi , R.P. , Wolters , P.J. , Abate , A.R. , et al. ( 2019 ). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage . Nat. Immunol . 20 , 163 – 172 . OpenUrl CrossRef PubMed 82. ↵ Guo , H. , and Li , J . ( 2021 ). scSorter: assigning cells to known cell types according to marker genes . Genome Biol . 22 , 69 . OpenUrl CrossRef PubMed 83. ↵ Pombo Antunes , A.R. , Scheyltjens , I. , Lodi , F. , Messiaen , J. , Antoranz , A. , Duerinck , J. , Kancheva , D. , Martens , L. , De Vlaminck , K. , Van Hove , H. , et al. ( 2021 ). Single-cell profiling of myeloid cells in glioblastoma across species and disease stage reveals macrophage competition and specialization . Nat. Neurosci . 24 , 595 – 610 . OpenUrl CrossRef PubMed 84. ↵ Song , L. , Pan , S. , Zhang , Z. , Jia , L. , Chen , W.-H. , and Zhao , X.-M . ( 2021 ). STAB: a spatio-temporal cell atlas of the human brain . Nucleic Acids Res . 49 , D1029 – D1037 . OpenUrl CrossRef PubMed 85. ↵ Müller , S. , Cho , A. , Liu , S.J. , Lim , D.A. , and Diaz , A . ( 2018 ). CONICS integrates scRNA-seq with DNA sequencing to map gene expression to tumor sub-clones . Bioinformatics 34 , 3217 – 3219 . OpenUrl CrossRef PubMed 86. ↵ Elosua-Bayes , M. , Nieto , P. , Mereu , E. , Gut , I. , and Heyn , H . ( 2021 ). SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes . Nucleic Acids Res . 49 , e50 . OpenUrl CrossRef PubMed 87. ↵ Cable , D.M. , Murray , E. , Zou , L.S. , Goeva , A. , Macosko , E.Z. , Chen , F. , and Irizarry , R.A . ( 2022 ). Robust decomposition of cell type mixtures in spatial transcriptomics . Nat. Biotechnol . 40 , 517 – 526 . OpenUrl CrossRef PubMed 88. ↵ Vahid , M.R. , Brown , E.L. , Steen , C.B. , Zhang , W. , Jeon , H.S. , Kang , M. , Gentles , A.J. , and Newman , A.M . ( 2023 ). High-resolution alignment of single-cell and spatial transcriptomes with CytoSPACE . Nat. Biotechnol . 41 , 1543 – 1548 . OpenUrl CrossRef PubMed 89. ↵ Gu , Z. , Schlesner , M. , and Hübschmann , D . ( 2021 ). cola: an R/Bioconductor package for consensus partitioning through a general framework . Nucleic Acids Res . 49 , e15 . OpenUrl CrossRef PubMed 90. ↵ Jin , S. , Guerrero-Juarez , C.F. , Zhang , L. , Chang , I. , Ramos , R. , Kuan , C.-H. , Myung , P. , Plikus , M.V. , and Nie , Q . ( 2021 ). Inference and analysis of cell-cell communication using CellChat . Nat. Commun . 12 , 1088 . OpenUrl CrossRef PubMed 91. ↵ Pham , D. , Tan , X. , Balderson , B. , Xu , J. , Grice , L.F. , Yoon , S. , Willis , E.F. , Tran , M. , Lam , P.Y. , Raghubar , A. , et al. ( 2023 ). Robust mapping of spatiotemporal trajectories and cell-cell interactions in healthy and diseased tissues . Nat. Commun . 14 , 7739 . OpenUrl CrossRef PubMed 92. ↵ Jin , S. , Plikus , M.V. , and Nie , Q . ( 2025 ). CellChat for systematic analysis of cell-cell communication from single-cell transcriptomics . Nat. Protoc . 20 , 180 – 219 . OpenUrl CrossRef PubMed 93. ↵ Cang , Z. , Zhao , Y. , Almet , A.A. , Stabell , A. , Ramos , R. , Plikus , M.V. , Atwood , S.X. , and Nie , Q . ( 2023 ). Screening cell-cell communication in spatial transcriptomics via collective optimal transport . Nat. Methods 20 , 218 – 228 . OpenUrl CrossRef PubMed 94. ↵ Csárdi , G. , and Nepusz , T . ( 2006 ). The igraph software package for complex network research . 1695 . View the discussion thread. Back to top Previous Next Posted April 22, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions Saad M. Khan , Anthony Z. Wang , Rupen R. Desai , Colin R. McCornack , Rui Sun , Sonika M. Dahiya , Jennifer A. Foltz , Ngima D. Sherpa , Lydia Leavitt , Timothy West , Alexander F. Wang , Aleksandar Krbanjevic , Bryan D. Choi , Eric C. Leuthardt , Bhuvic Patel , Al Charest , Albert H. Kim , Gavin P. Dunn , Allegra A. Petti bioRxiv 2025.04.04.647096; doi: https://doi.org/10.1101/2025.04.04.647096 Share This Article: Copy Citation Tools Mapping the spatial architecture of glioblastoma from core to edge delineates niche-specific tumor cell states and intercellular interactions Saad M. Khan , Anthony Z. Wang , Rupen R. Desai , Colin R. McCornack , Rui Sun , Sonika M. Dahiya , Jennifer A. Foltz , Ngima D. Sherpa , Lydia Leavitt , Timothy West , Alexander F. Wang , Aleksandar Krbanjevic , Bryan D. Choi , Eric C. Leuthardt , Bhuvic Patel , Al Charest , Albert H. Kim , Gavin P. Dunn , Allegra A. Petti bioRxiv 2025.04.04.647096; doi: https://doi.org/10.1101/2025.04.04.647096 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Cancer Biology Subject Areas All Articles Animal Behavior and Cognition (7618) Biochemistry (17635) Bioengineering (13859) Bioinformatics (41846) Biophysics (21401) Cancer Biology (18534) Cell Biology (25422) Clinical Trials (138) Developmental Biology (13352) Ecology (19860) Epidemiology (2067) Evolutionary Biology (24285) Genetics (15582) Genomics (22463) Immunology (17700) Microbiology (40298) Molecular Biology (17141) Neuroscience (88424) Paleontology (666) Pathology (2825) Pharmacology and Toxicology (4813) Physiology (7633) Plant Biology (15107) Scientific Communication and Education (2042) Synthetic Biology (4284) Systems Biology (9808) Zoology (2267)