Spatiotemporal modeling reveals high-resolution invasion states in glioblastoma

Manoharan, Varsha Thoppey; Abdelkareem, Aly; Gill, Gurveer; Brown, Samuel; Gillmor, Aaron; Hall, Courtney; Seo, Heewon; Narta, Kiran; Grewal, Sean; Dang, Ngoc Ha; Ahn, Bo Young; Osz, Kata; Lun, Xueqing; Mah, Laura; Zemp, Franz; Mahoney, Douglas; Senger, Donna L.; Chan, Jennifer A.; Morrissy, A. Sorana

doi:10.1186/s13059-024-03407-3

Research
Open access
Published: 10 October 2024

Spatiotemporal modeling reveals high-resolution invasion states in glioblastoma

Varsha Thoppey Manoharan^1,2,
Aly Abdelkareem^1,2,
Gurveer Gill^1,2,
Samuel Brown^1,2,
Aaron Gillmor^1,2,
Courtney Hall^1,2,
Heewon Seo^1,2,
Kiran Narta^1,2,
Sean Grewal^1,2,
Ngoc Ha Dang²,
Bo Young Ahn²,
Kata Osz²,
Xueqing Lun²,
Laura Mah^2,3,
Franz Zemp^1,2,
Douglas Mahoney^2,3,7,
Donna L. Senger^4,5^na1,
Jennifer A. Chan^2,6,7^na1 &
…
A. Sorana Morrissy^1,2,7^na1

Genome Biology volume 25, Article number: 264 (2024) Cite this article

Metrics details

Abstract

Background

Diffuse invasion of glioblastoma cells through normal brain tissue is a key contributor to tumor aggressiveness, resistance to conventional therapies, and dismal prognosis in patients. A deeper understanding of how components of the tumor microenvironment (TME) contribute to overall tumor organization and to programs of invasion may reveal opportunities for improved therapeutic strategies.

Results

Towards this goal, we apply a novel computational workflow to a spatiotemporally profiled GBM xenograft cohort, leveraging the ability to distinguish human tumor from mouse TME to overcome previous limitations in the analysis of diffuse invasion. Our analytic approach, based on unsupervised deconvolution, performs reference-free discovery of cell types and cell activities within the complete GBM ecosystem. We present a comprehensive catalogue of 15 tumor cell programs set within the spatiotemporal context of 90 mouse brain and TME cell types, cell activities, and anatomic structures. Distinct tumor programs related to invasion align with routes of perivascular, white matter, and parenchymal invasion. Furthermore, sub-modules of genes serving as program network hubs are highly prognostic in GBM patients.

Conclusion

The compendium of programs presented here provides a basis for rational targeting of tumor and/or TME components. We anticipate that our approach will facilitate an ecosystem-level understanding of the immediate and long-term consequences of such perturbations, including the identification of compensatory programs that will inform improved combinatorial therapies.

Background

Glioblastoma is the most common malignant brain tumor in adults—incurable despite multi-modal treatment with maximal safe surgical resection, radiation, and chemotherapy [1]. Part of the treatment challenge in GBM stems from its highly invasive phenotype, wherein individual tumor cells move through normal tissues diffusely or follow perivascular routes or white matter tracts to spread far beyond the main tumor mass [2]. Although surgery removes the bulk of the tumor, these infiltrating cells are left behind as they are not only difficult to detect but also impossible to remove without unacceptable neurologic consequences. The residual cells can then continue to evolve and adapt to the selective pressures of conventional and rational therapies—a process that is multifaceted and that involves genetic heterogeneity, phenotypic plasticity, and the ability to engage with and co-opt the tumor microenvironment (TME) into a pro-tumorigenic state [2,3,4,5,6,7]—to result in recurrence. Understanding of the biology of the invasive front and delineating the mechanisms by which these cells engage with their surrounding normal cells and environment to promote the malignant phenotype is thus of high clinical importance.

Multiple factors have hindered our ability to understand the invasive front, however, and its relationships with the rest of the spatial glioblastoma ecosystem. First, surgery removes the most solid tumor tissue, but the outermost regions of infiltration are often inaccessible and undersampled due to clinical limitations. Second, surgery and tissue banking typically yield small and unoriented tissue fragments with an unknown spatial relationship relative to each other. Thus, most GBM literature describes heterogeneity within only small regions of the highly cellular or highly vascularized resectable portion of the tumor and are unable to capture overall tumor organization. Third, studies to date have implicated multiple factors in invasion, including signaling pathways [8, 9], components of the extracellular matrix (ECM) [10, 11], and interactions with nearby or distal cells of the tumor [12, 13] and TME [14], but there have been limited approaches and tools to contextualize the many processes as part of an ecosystem of interdependent processes within the tumor’s overall organization [15]. These challenges have collectively limited our ability to assess global patterns of adaptation to local tissue contexts, for instance, invasion along white matter tracts or perivascular routes within the same tumor.

Fortunately, new technologies for spatial profiling, which can delineate the organization of tumor transcriptional cell states in relation to each other, to genetic diversity, and to metabolic and cellular diversity of the TME, have the potential to greatly increase our understanding of tumor ecosystems [16,17,18,19,20]. Nevertheless, the spatial profiling platforms available to date that offer global transcriptome coverage—and that could therefore support exploratory studies—do not have single-cell resolution (NanoString GeoMx, 10X Genomics Visium). The resulting data are therefore mixtures of cell types and states that require computational deconvolution. Strategies to deconvolute spatial data employ either supervised approaches (requiring matched single-cell data to infer cellular composition within each profiled region) [21,22,23], unsupervised approaches based on matrix factorization or probabilistic modelling [24,25,26,27] (to identify latent gene expression programs representing cell types or states), or semi-supervised approaches [28, 29]. Regardless of strategy, most current applications of deconvolution have not robustly disambiguated low-frequency signals and have primarily focused on regions where signals of interest are > 20% of the total. This excludes areas of diffuse tumor infiltration [16, 17, 20] as well as lower frequency cells and states of the TME. Altogether, these limitations hold back key aspects of tumor biology from comprehensive study, including components of the TME and cellular interactions that drive tumor growth and invasion.

Here, we aim to address these challenges by coupling transcriptome-wide spatial profiling of xenografted GBM cells with temporal sampling throughout tumor progression and an analytic workflow that enables highly specific and sensitive detection of cellular programs (Fig. 1a). A key strength of the xenograft strategy is the genomic distinction between human tumor cells and mouse-derived TME, which we leverage to boost signal detection in regions with low tumor content. This approach captures transcriptional phenotypes across the whole tumor including areas of invasion. The GBM lines selected represent multiple common genetic drivers, allowing the identification of both universal and genotype-enriched transcriptional programs that can be linked back to external patient cohorts. To identify molecular programs, we developed a computational workflow centered on unsupervised deconvolution that does not require matched single-cell data and enables de novo discovery of both known and novel expression programs corresponding to cell types, cell activities, or combinations thereof. This led to the deconvolution of the mouse brain and TME with far greater granularity than previously achieved for spatial data, distinguishing 71 cell types and anatomic structures, and 19 TME-related programs in this cohort. Many of the TME programs (including astrocytic, myeloid, and vascular cell types and states) showed spatiotemporal variation in abundance, reflecting dynamic changes during tumor growth and invasion. We further identified multiple highly resolved human tumor cell programs and explored within-tumor and tumor-TME crosstalk across regions of variable tumor density, identifying unique interactions along distinct routes of invasion. Altogether, our work provides insight into previously elusive aspects of the GBM ecosystem and enables systematic exploration of the invasive front.

Results

Spatial transcriptome profiling and unsupervised program discovery in GBM xenograft models

To comprehensively capture spatial transcriptional heterogeneity in vivo, we selected six well-characterized brain tumor-initiating cell (BTIC) lines established at our institution from surgical samples [5, 30, 31]. These derive from four patients spanning a range of clinical variables including sex, exposure to standard therapy, MGMT methylation status, and common GBM molecular drivers (Fig. 1b, Additional file 1: Table S1). In addition to the diversity of genetic drivers, our cohort further included BTIC lines from two patients (BT143 and BT238) that captured acquired phenotypic diversity during tumor evolution. In these patients, multiple lines were concurrently derived from anatomically distinct regions of the tumor as defined by pre-operative MRI, including the densely growing core (x), the contrast-enhancing highly vascularized tumor margin (y), and the highly diffuse leading edge of the tumor (z) that is typically outside of the surgical margin. Lines from each patient shared genomic drivers but maintained distinct growth phenotypes in vivo, with x lines growing densely relative to the diffusely infiltrating y and z lines (Fig. 1c; Additional File 2). Similar observations have been made for other lines similarly derived from edge versus core surgical samples [12] and in the context of adaptation to hypoxia [19]. The consistency of x/y/z growth patterns within and across cohorts indicated that predictable phenotypic adaptations can be acquired by GBM cells during tumor evolution, that these adaptations relate to spatial context (e.g., core vs edge), and that they are heritable, likely being fixed in the genome or epigenome [12, 19]. Our x/y/z BTIC lines thus offer a rare opportunity to investigate expression programs underlying dense versus diffuse growth. Finally, to capture invasion dynamics across tumor stages, our cohort includes early, mid, and late timepoints of growth based on known time to endpoint (ranging between 76 and 428 days across lines) (Fig. 1b, Additional file 1: Table S1).

We profiled 51,952 individual spatial transcriptomic RNAseq (stRNAseq) measurements (spots) using Visium from a total of 23 samples (Fig. 1c, “Methods”) spanning all lines and timepoints. In many cases, we were able to fit a complete coronal brain section diagonally within the capture area, ensuring profiling across the whole tumor and invasive front. In other cases, the injection side (i) and contralateral side (c) were mounted on separate capture areas (e.g., BT143y/z endpoint samples, Fig. 1c). Further, we also profiled sequential sections cut 30–40 µm apart but did not observe notable biological variability (data not shown). We used human-to-mouse transcriptome admixture to calculate the relative contribution of tumor cells to the transcriptional output of each spot. Admixture ranged from 0 in spots with 100% mouse cells to 1 in spots with 100% tumor cells (Fig. 1c–e, Additional file 2: Fig. S1a, b). The sensitivity with which we could distinguish mouse from human cells enabled us to focus our next analyses on regions of high (80–100%; D4), moderate-high (50–80%; D3), moderate-low (20–50%; D2), and low (5–20%; D1) tumor cell density, as well as spots of mouse brain without tumor (D0). We observed highly variable levels of tumor density among lines and timepoints across tumor regions defined in this way. BT143x stood out as the most densely growing line with many spots at > 80% tumor density by endpoint (Fig. 1e). BT161 and BT238 were the next-most densely growing, with many spots having 50–80% tumor cell density at endpoint. In contrast, BT134 endpoint tumors and all earlier timepoints across lines grew diffusely.

Gene expression program discovery in human and mouse cells

Since stRNAseq does not have single-cell resolution, we utilized unsupervised deconvolution to identify admixed cell types and states separately for human and mouse data. We used consensus non-negative matrix factorization (cNMF) [32] to identify robust transcriptional programs across all 23 samples together and quantify their relative usage within spots. Factorization yielded 15 human tumor cell programs (referred to as h1-h15) and 90 mouse brain and TME programs (referred to as m1-m90), in line with the greater transcriptional diversity of the mouse brain [33, 34]. To decipher if programs represented cell types, cell states, or anatomical structures, we (i) used established reference marker genes for normal mouse brain cell types [28, 35, 36], brain structures [37], TME-specific cell types [7, 28, 38], and human GBM states [4, 5, 17, 39,40,41,42,43,44], (ii) assessed program-specific enrichment of biological pathways, and (iii) compared spatial patterns of program usage to known mouse brain structures from the Allen Brain Atlas Common Coordinate Framework (CCFv3) (Additional file 1: Table S2a; Methods; Additional File 2).

Human tumor programs were most diverse in the core, with three to four discernable programs per spot in high-density regions and one to two programs present per spot in lower-density regions (Additional file 2: Fig. S1c). Mouse programs included 11 cell activities and 18 cell types resident in the normal brain, 19 TME-specific or enriched cell types and states, and 42 programs representing combinations of cell types that could not be further disambiguated at this level of factorization (Additional file 1: Table S2b). Comparison to the ABA confirmed that most of these combination programs corresponded to anatomic regions or structures of the mouse brain (Additional file 1: Table S2b; Additional file 2). We detected an average of five normal brain mouse programs per spot in regions without GBM (D0) and in the tumor leading edge (D1), with diversity decreasing at higher tumor densities (D2-D4; Additional file 2: Fig. S1c). In contrast, only one to two TME-related programs were observed per spot in tumor regions. Overall, across both human and mouse data, our approach was able to quantify the usage of up to eight programs per spot, showing far greater sensitivity than previously achieved with stRNAseq analyses (i.e., ~ 2 programs/spot) [16].

To ensure our workflow generated meaningful results, we evaluated a set of mouse programs expected to closely match previous literature for cell types and anatomic structures. A first example is program m25, which corresponds to ependymal cells based on the enrichment of marker genes (Additional file 1: Table S2d, 5b). Ependymal cells form a single-cell layer restricted to the lining of the ventricles, and indeed, m25 was localized to the ventricular lining with no background signal elsewhere (Fig. 1f–g). Usage values ranged between 0.3 and 0.5, indicating ependymal cells make up one-third to one-half of the signal in these spots, in line with known cell composition in the ventricular lining [45]. Analysis of marker genes and transcription factor (TF) activity further validated m25 as an ependymal program (Fig. 1h, i). Within the subventricular zone, we also identified progenitor cell programs at lower cellular frequency per spot (3–10%), including gliogenic progenitors concentrated in the dorsolateral ventricular region (m58), and two proliferating progenitor programs (m6, m83) located in the lateral ventricular lining and with diffuse spread into the brain parenchyma (Fig. 1f, g). As a second example, we could distinguish pericytes from endothelial cells (Additional file 2: Fig. S1e), even though these cell types always co-occur spatially within the vasculature, represent a minority of signals within a given Visium spot (max of 8% usage for the pericyte program), and have not previously been deconvoluted in spatial datasets [16, 17, 20, 28]. In a final example, we highlight the identification of cortical layers with high resolution (Additional file 2: Fig. S1f), including the outermost meninges (m4) and vascular leptomeningeal smooth muscle cells (m11). The cortical layers are arranged in a spatially overlapping manner that reflects the changing composition of cell types and states along the radial axis of the brain. This compositional gradient is quantitatively captured by the program usage values, recapitulating the level of spatial overlap between cortical programs (Additional file 2: Fig. S1f). Based on these observations, we conclude that our deconvolution and annotation approach is highly specific (enabling identification of unique cell types, states, and anatomic structures), is highly sensitive (capable of deconvoluting signals with low overall signals in the range of 3–10%), can identify both spatially coherent and diffuse programs, and provides interpretable usage values that quantify up to eight programs within each Visium spot.

Tumor programs span progenitor states to invasion phenotypes

Having established the sensitivity and specificity of our deconvolution approach, we characterized 15 de novo human tumor cell programs using a similar strategy (Fig. 2a–d, Additional file 2: Fig. S2a, b and Additional file 1: Table S2b). This included pathway-based and marker gene-based annotations, transcription factor (TF) activity-based similarity, and spatial co-occurrence. These analyses were used to label and assign each program to broader themes based on all information (Additional file 1: Table glycogenic S3a-e). The 15 human tumor programs are fully described in the context of known GBM transcriptional classes, pathways, and TME spatial neighborhoods in Supplemental Results.

We classified programs into six major themes, representing progenitor states, cell cycle, metabolism, astrocytic-like, oligodendrocytic-like, and invasion. Progenitor programs (h5_progenitor, h13_DNArepair) were most prevalent in high tumor cell density regions (D4) (Fig. 2e, Additional file 2: Fig. S2c, Additional file 1: Table S3f), along with one of the cell cycle programs (h7_telomere), indicating that cycling progenitors preferentially reside there. The OC-like programs (h12_OC1, h14_OC2) also showed enrichment in denser tumor regions (Fig. 2e), with h14_OC2 representing a genetic subclone within the BT143x cell line (Additional file 2: Fig. S2d-g, Supplemental Materials). Outside of the dense tumor core, a trio of programs formed a gradient of outward tumor expansion reflecting a phenotypic arrangement centering on regions of hypoxia (h9_hypoxia most centrally located), extending to regions of more diffuse tumor expansion without hypoxia (h2_AC), and finally to invasion into the normal brain (h11_invason) (Fig. 2f, g). This was reminiscent of hypoxia-centered cellular organization within human tumors [16, 46], indicating our xenograft cohort captures this facet of tumor biology and, importantly, extends the previous work by revealing a program of diffuse invasion (h11). The h11_invasion program scored highly for leading-edge gene sets (LE_IvyGAP; Fig. 2b), comprised a majority of D1–D2 spots at early to late timepoints, and was the only tumor program with significant over-representation in D1 regions that was common among multiple individual patients (Fig. 2g, Additional file 1: Table S3g). Finally, we noted that a few programs did not have cohort-wide usage (Fig. 2d). Instead, these were prevalent within individual patients suggesting that specific tumor genotypes could be linked to unique repertoires of tumor cell programs and that larger xenograft cohorts will likely reveal additional insights. In support of this, previous work establishing that EGFR plays a role in errant neovascularization [46] was in line with the enrichment of angiogenic tumor cell programs (h3, h10) in the tumor line BT134, which harbors high-level EGFR amplification (93 copies).

Finally, we gauged how prevalent xenograft tumor programs were in other datasets by comparing them with similarly identified programs in external cohorts of human tumors [5, 17, 30], brain tumor-initiating cell lines [30], and xenografts [30] (Additional file 2: Fig. S2i, Additional file 1: Table S3h). This comparison showed the highest correlations for progenitor, cell cycle, and metabolism programs among datasets. We saw less concordance for genotype-specific programs, indicating that these programs are less prevalent across GBM patients.

Tumor program association with survival and transcriptional subtypes

We next asked how the 15 human programs related to survival differences, anticipating that associations may be context-dependent with respect to the GBM transcriptional subtypes (classical, mesenchymal, proneural). For each tumor program, we first identified the genes most strongly contributing to program identity (i.e., top-scoring genes) (Fig. 3a, Additional file 2: Fig. S3, Additional file 1: Table S4a), then performed network analysis to select the subset acting as network hubs based on protein–protein interactions (Fig. 3a, Additional file 2: Fig. S3, Additional file 1: Table S4b). Hub genes were distinct among programs (Fig. 3b) and well aligned with previously described program themes (Fig. 2a). We performed survival analyses in the TCGA cohort using both top-scoring and hub genes and found that hub genes had a strong association with survival across all transcriptional subtypes (Fig. 3c, f, i, Additional file 1: Table S4e). This was true in both a classic survival analysis (comparing survival of patients ranked by gene set expression; Fig. 3d, g, j) and when comparing gene set enrichment in patients first stratified by survival outcome (Fig. 3e, h, k). Several programs emerged as robustly associated with survival in a subtype-specific manner. For instance, high expression of h1_metabolism and h12_OC1 was robustly associated with poor survival in classical tumors. The GSEA normalized enrichment scores (NES) of h1 and h12 gene sets were very highly correlated across the TCGA cohort (cor = 0.85; Additional file 1: Table S4f, g), indicating that OC-like cells enacting h1 metabolic activities (primarily cholesterol biosynthesis) operate as a unit, together influencing cell phenotypes related to survival. Mesenchymal tumors were stratified by cell cycle/epigenetic (h8_epigenetic) and invasion (h11_invasion) programs. Lower NES correlations between h8 and h11 (cor = 0.3; Additional file 1: Table S4f, g) indicated these programs independently influenced survival. Proneural tumors were also stratified by cell cycle (h7_telomere) and the three programs encompassing the hypoxia-to-invasion gradient (h9_hypoxic, h2_AC1, h11_invasion). NES correlation values among these programs in TCGA followed the same graded pattern of co-occurrence observed in xenografts (h9-h2-h11), strongly supporting their spatial and phenotypic relationship in human patients (Additional file 1: Table S4f, g)

Having established the prognostic value of the tumor programs, we more deeply investigated how these related to known GBM molecular transcriptional states (Neftel) [4]. This was specifically motivated by the weaker observed match of the strongly prognostic invasion programs (h2, h11) to NPC-like states (Fig. 2b), despite previous links between NPC-like states and GBM invasion [4, 6, 47]. Given that multiple states typically co-exist within a tumor and that state transitions can reflect adaptations to the TME [4, 6, 41], we anticipated correspondence between Neftel states and tumor density. As expected, each xenograft harbored all four states (Additional file 2: Fig. S4a-c), with NPC-like and MES states less abundant overall, but with significantly higher prevalence in D1–D2 regions relative to the tumor core (Additional file 2: Fig. S4a). We also found that tumor baseline (majority) states were evident in each xenograft sample, stable across timepoints and tumor density regions, independent of cNMF program usage, and therefore likely intrinsic to each line, as previously observed [6] (Additional file 2: Fig. S4c). By then calculating the rate of transition away from tumor baseline in spots with the usage of individual programs, we established that transitions were common across programs (Additional file 2: Fig. S4d). This supported the high plasticity of Neftel states among the tumor programs defined here and specifically for h11_invasion in regions of diffuse infiltration (D1) where dynamic transitions from an AC baseline to an NPC-like state predominated across all timepoints. This supported h11 as an invasion program largely orthogonal to the previously established NPC-like state.

Microenvironment programs encompass cell types and states

Of the 90 mouse programs, we focused on 19 that were either specific to or highly enriched in regions of the tumor, each with dynamic spatial and temporal kinetics. These TME programs were broadly categorized into 11 cell types and 8 cell activities based on marker gene, TF-activity and pathway annotations, and co-localization (Fig. 4a–d, Additional file 2: Fig. S5a-c, Additional file 1: Table S2b, 5a-d, 3e). Programs were labeled as cell types when marker gene-based annotations were unambiguous and strong; otherwise, labels thematically reflected highly enriched pathways indicating cell activities. All programs of the TME were detected in all lines (Fig. 4e). We describe the resulting in vivo spatial xenograft TME catalogue in Supplemental Results, covering astrocytes, vasculature, immune cells (microglia, monocytes, macrophages), and multiple activity programs.

TME programs were broadly organized along a spatial axis corresponding to tumor density. For instance, we observed that mature astrocytes (m74_Astro1) widespread and resident throughout the normal mouse brain were prevalent in the invasive tumor front but excluded from dense tumor regions (Fig. 4f, h). m49_AstroRNA, representing normal astrocytic activities (related to RNA processing, proliferation, glutamatergic synapses), was prevalent in astrocytes localized to the normal brain and tumor periphery (Supplementary Fig. 5d, Additional file 1: Table S5e). Within the context of the tumor, however, a dramatic state shift toward a reactive program (m43_AstroReac) was enacted across the full spectrum of tumor density and sustained at all timepoints of disease progression (Fig. 4f, Additional file 2: Fig. S5d, Additional file 1: Table S5e). Astrocytes in this reactive state were enriched for terms relating to proliferation, inflammatory response, and ECM organization and had a strong match to established reactive signatures [48] (Additional file 2: Fig. S5e). Altogether, these state change patterns supported phenotypic co-option of normal astrocytes along the advancing tumor front. We were surprised to also observe the tumor association of a regional astrocytic cell subtype or activity (m44_AstroHY), characterized by pathway enrichment of terms relating to precursor cell proliferation. This seemingly homeostatic program, highly used within hypothalamic regions at all timepoints (Fig. 4h), also showed early and sustained tumor enrichment (Additional file 2: Fig. S5d, Additional file 1: Table S5e), indicating that tumor-proximal astrocyte transitions to this state could play an important a role within the context of the GBM TME.

We observed similar regional specificity of endothelial programs, with normal vasculature in the normal brain and invasive front (m52_Endo1, m59_EndoLV) gradually replaced by tumor-enriched vascular programs (m86_Endo2, m30_EndoHyp) within the denser and more hypoxic regions of tumor [49, 50] (Fig. 4c, f). Microglial programs were also present diffusely throughout the normal brain (Fig. 4f, g), had an early response to sites of injury (injection tract), and persisted long-term within lower-density tumor areas (D2–D3) (Fig. 4f, Additional file 2: Fig. S5d, Additional file 1: Table S5e). Of these, m7_MG1 was more abundant overall, scoring strongly for response to injury and antigen processing and presentation, while m77_MG2 was involved in apoptotic cell clearance (Additional file 1: Table S5a), indicating these microglial subpopulations play different roles within the tumor. Two monocyte-derived macrophage (MDM) programs were rapidly recruited to early lesions and were otherwise absent from the normal brain (Fig. 4f, g). Although these two MDM programs had broadly similar spatial distributions, they enacted distinct activities possibly related to their micro-local spatial contexts—m40_MDM1 was enriched in terms relating to cell–matrix adhesion, chemotaxis, and migration, while m84_MDM2 scored highly for macrophage proliferation and phagocytosis (Fig. 4b, Additional file 1: Table S5a). Of the multiple immune cell activities identified (Supplemental Results), we highlight the cytotoxic program m57_Cytotoxic as primarily enacted by MG and MDMs (based on co-localization; Additional file 2: Fig. S5c, Additional file 1: Table S5d) and observed as the most prevalent immune activity in D4, indicating that cytotoxicity plays a more central role in denser regions (Fig. 4f, i).

We sought to compare this highly detailed catalogue of TME programs with programs identifiable in external cohorts. We used cNMF to similarly identify programs in five additional datasets, including single-cell and spatial data from non-tumor-bearing developing and adult mouse brain (Kleinman, Bayraktar) [28, 35], and single-cell data from CD45 + cells sorted from syngeneic (Movahedi) [38] and xenograft (Senger) [7] GBM models (Additional file 2: Fig. S5f, Additional file 1: Table S5f). After factorizing these datasets, we identified the best match between the TME programs and the resulting ensemble of cNMF solutions, observing general agreement along expected themes. For instance, astrocytic and endothelial programs had a higher correlation with programs in the normal brain datasets compared to CD45 + enriched datasets where non-immune cells were depleted. Interestingly, we saw better matches between m52_Endo1 and Visium rather than single-nuclei data (Bayraktar), indicating possible loss of normal vascular cells during single-cell sample processing. Brain-resident microglial programs were also more prevalent in the normal brain versus GBM datasets, whereas the opposite was true for MDM programs. Monocytes also showed specific enrichment in the GBM datasets, as expected from their recruitment during tumor progression. Overall, broad concordance among cohorts and data types indicated that the TME program catalogue defined here represents meaningful biological signals.

Tumor microenvironment crosstalk

Cellular communication is vital to the tumor ecosystem, with cell–cell contacts, tumor extracellular matrix (ECM) interactions, and secreted signaling all pivotal to tumor growth, adaptation, and invasion. To better understand how the glioma and non-malignant cells described here communicate across distinct niches, we quantified directional ligand receptor (LR) signaling using CellChat [51], expecting to capture distinct interactions in dense versus invasive regions reflective of cellular changes in TME composition. We surveyed all possible tumor–tumor, tumor–TME, and TME–TME interactions, stratified by density region (Fig. 5a, Additional file 2: Fig. S6a). In all, 63 pathways were involved in significant crosstalk, comprising 8 major groups based on the directional involvement of human versus TME ligands and receptors (LR1-LR8; Fig. 5b). Within the TME (LR1) for instance, we observed that CSF-CSF1R signaling in D1–D2 regions involved microglia as the primary signal-receiving cells, based on high scores for CSF1R in these programs (CSF1R gene rank m7 = 16, m77 = 22; Additional file 1: Table S6a). In LR2, significant Spp1-CD44 crosstalk took place between macrophages and reactive astrocytes at early timepoints and also between macrophages and tumor cells at later stages (Additional file 2: Fig. S6b-e, Additional file 1: Table S6a). In this case, the spatial data enabled the distinction of mouse versus human receivers based on spatial co-localization of mouse sender and receiver cells at early timepoints (Additional file 2: Fig. S6e). LR8 involved within-tumor signals related to invasion and stemness. In this group, CALCR signaling stood out based on the CALCRL receptor as the top-scoring gene in the h5_progenitor program (Additional file 1: Table S6a). CALCRL had previously been linked to glioma cell proliferation and has been negatively associated with glioma prognosis and promotion of angiogenesis [52]. Based on these associations, we speculate that the observed enrichment of the h5_progenitor program along with tumor-enriched vasculature programs in D4 could be related to vascular stem cell niche development and maintenance [53, 54]. Indeed, CALCRL is a marker of stemness and a promising candidate therapeutic target in other diseases [55].

By far, the most abundant group of pathways encompassed multidirectional signaling within and between tumor and TME (LR3), including MK and PTN pathways, AGRN, JAM, and NCAM (cell adhesion), tumorigenic NOTCH, PDGF, FGF, and angiogenic VEGF (Fig. 5a, b, Additional file 1: Table S6b, c). We noted that multiple LR3 pathways converged on the formation and maintenance of the ECM (tenascin, laminin, collagen, fibronectin). Laminin and tenascin ligands were primarily human (LAMB2, LAMA4, TNC, TNR), while fibronectin ligands were made by both mouse and human cells (Fig. 5d-f). Mouse (but not human) Fn1 predominated at early timepoints (Additional file 2: Fig. S6d) and in D1 regions (Fig. 6f), providing the majority of this tumor ECM component in areas of diffuse infiltration. Since vascular programs were strongly associated with Fn1 (m52, m30, m86), we conclude that a major invasion route in D1 is along vessels. Collagens were also derived from both tumor and TME components, with higher contribution from the tumor (Additional file 2: Fig. S6a-c, Additional file 1: Table S6c). The deposition of human collagens was spatially distinct, with COL6A1 dominant in D1–D3 and specifically associated with the h11_invasion program. In D4, COL9A2 and COL9A3 were more prevalent and associated with h1_metabolism and h7_telomere, indicating changing ECM deposition based on density and niche (Additional file 2: Fig. S6b, c). Mouse collagens Col4a1 and Col4a2 were highest in D4 and derived from the tumor-associated vasculature (with high program scores in m30_EndoHyp and m86_Vasc2). Altogether, cells of both tumor and TME differentially contributed to the tumor ECM, resulting in a dynamic structural scaffold underlying gradients of invasion and forming the basis for distinct tumor cell niches.

Within areas of diffuse infiltration (D1), we observed significant invasion-associated midkine signaling through human receptors PTPRZ1 and LRP1, both highly scoring genes in h11 and h2 programs. These receptors responded to midkine ligands from astrocytic sources (Mdk highly scoring in m74_Astro1 and m44_AstroHY). Notch signaling also stood out as a main source of mouse–human crosstalk in D1 (LR5; Fig. 5g-i, Additional file 2: Fig. S6d, Additional file 1: Table S6a). Notch1 is important in GBM invasion along white matter tracts [8] and for GBM cell survival within the perivascular niche [9]. Our data supports that in D1 regions, the human NOTCH1 receptor (highly scoring in h2 and h11) can bind to mouse ligands Dlk1 and Jag2. These ligands were associated with normal vasculature (Dlk1 in m52_Vasc1) and hypothalamic astrocytes (Jag2 in m44_AstroHYP), highlighting these cellular programs as distinct TME participants in the invasion signaling axis.

Molecular and structural contributors to invasion along distinct routes of tumor cell travel

Invasion involves tumor cell movement along white matter tracts, along perivascular routes, as well as directly through the brain parenchyma [2, 8, 9, 56]. We attempted to identify and characterize these routes through the association of mouse brain programs with tumor invasion, followed by CellChat [51] and differential expression analyses [57, 58]. We first assessed the over-representation of all 90 mouse brain programs in D1 spots, expecting to see enrichment of invasion-relevant cell types or anatomic regions (Fig. 6a, b). Indeed, white matter (WM)-related programs had the highest association with D1 tumor regions (Additional file 2: Fig. S7a, Additional file 1: Table S7a). We therefore selected D1 spots with high usage of mouse WM programs as representative of human tumor cells traveling along white matter tracts. Likewise, we observed that usage of caudoputamen programs (m1, m73) was also significantly enriched in D1 regions (Fig. 6a, b, Additional file 2: Fig. S7a). Further analysis showed that many caudoputamen spots also had usage of vascular programs (expected from the high diversity of programs identifiable per spot; Additional file 2: Fig. S1a); we therefore further stratified caudoputamen spots based on co-usage of the four vascular programs (Additional file 2: Fig. S7b). This strategy distinguished the D1 tumor cells traveling along each type of vasculature (perivascular routes), from those moving directly through the brain parenchyma (i.e., caudoputamen spots without vasculature). We observed that NPC-like and OPC-like states were more common in D1 tumor cells traveling along white matter tracts, while those within the brain parenchyma were skewed toward MES and AC-like states (Fig. 6c). Thus, although h11 and h2 invasion programs were highly plastic (Additional file 2: Fig. S4a-d), cell state decisions in the invasive front reflect adaptation to the local cellular context of specific routes.

We conducted another CellChat analysis on the resulting 6 groups of D1 spots: white matter (WM), parenchyma (Par), and perivascular routes (m52, m30, m86, m59) (Fig. 6d, Additional file 2: Fig. S7c-e). Collagen, fibronectin, and laminin interactions were significant and prevalent along all routes (Fig. 6d). Human integrin receptor interactions with various mouse collagens were largely stratified by vasculature programs, suggesting that tumor cells co-opted existing collagen scaffolds within the TME for movement. A notable exception involved human COL6A1 interactions in parenchyma and m52 regions, indicating that tumor cells specifically required and secreted this ECM component (Fig. 6e). Similarly, tumor cells not only bound to multiple mouse laminins (Fig. 6e; Additional file 1: Table S7b) but also contributed a select subset of ligands to the ECM (LAMA5 in the parenchyma and LAMB2 in the vasculature) (Fig. 6e). A single laminin was specific to the WM route (mouse Lama2) (Fig. 6e). This gene was highly ranked in program m71 (newly forming oligodendrocytes), indicating that this oligodendrocyte subtype contributes to the migration scaffold. Collectively, these results highlighted that TME-derived ECM could differentiate between invasion routes and, notably, that a subset of regionally specific components originated from the tumor cells themselves—indicating the importance of these molecules to glioma cell movement along distinct routes of travel.

CellChat analysis also revealed several signaling interactions specific to the parenchyma, including NRXN, CD99, PSAP, PTN, MK, CNTN, and NOTCH (Fig. 6d, Additional file 1: Table S7b). Of these, we highlight tumor (NRXN1) crosstalk with neuron-derived synaptic adhesion molecule neuroligins (Nlgn1, Nlgn2, Nlgn3) (Fig. 6d, e). This previously described mitogenic signaling axis co-opts excitatory neuronal activity toward tumor growth, suggesting that active neurons are playing a role in glioblastoma invasion through the parenchyma [14, 59]. In addition, we also further observed interactions between tumor cells via the NRCAM (neuronal cell adhesion molecule) receptor and the mouse Cntn1 ligand—a gene highly scoring in cholinergic and GABAergic neurons, suggesting that these neurons may play additional roles in invasion beyond excitatory activity (Fig. 6e, Additional file 1: Table S2d).

As an orthogonal approach to the identification of invasion-relevant genes, we performed differential expression analysis between spots with high versus low tumor density in each line (i.e., a program-agnostic approach; “Methods”) and shortlisted genes with recurrent upregulated expression in areas of low tumor density (Additional file 1: Table S7c). Moreover, to further distinguish between invasion routes, we performed differential expression analysis between WM and the other five groups of interest (parenchyma, m52, m30, m59, m86), shortlisting significant genes (Additional file 1: Table S7d, e). We observed high overlap of these differentially expressed genes with the top-scoring and hub genes of programs h2, h9, and h11 (and to a lesser extent the genotype-specific h6 and h10 programs) (Additional file 2: Fig. S7f).

Finally, we combined the evidence from all analyses focused on D1, including the full CellChat analysis (Fig. 5a), the invasion route-specific CellChat analysis, and the two differential expression analyses above. Intersecting the resulting genes with tumor programs revealed an overlap of 27 genes that were also top-scoring or network hub genes within the prognostic invasion-related programs h2 and h11 (Fig. 6f). Nine of these were hub genes in either h2 (PTN, PTPRZ1, LRP1, NRXN1, VEGFA) or h11 (SDC3, COL6A1, NCAM1, NOTCH2), while 2 genes were hubs in both (FN1, APP) (Fig. 6f). Moreover, in the TCGA cohort, three of these genes were independently associated with poor survival outcomes in proneural subtype tumors (TUBA1A, NRCAM, LGALS1; Fig. 6g). NRCAM is a neural adhesion molecule participating in cell proliferation, axon growth, and synapse formation during neural development and was previously observed to be overexpressed in brain tumors [60,61,62]. TUBA1A is a microtubule subunit with effector functions in neuronal migration and causally implicated in neurodevelopmental defects [63, 64]. LGALS1 (Galectin-1) has a role in GBM invasion through modulation of cell–cell and cell–matrix interactions [65] and also promotes an immunosuppressive TME [66]. Another two genes were independently and significantly associated with poor survival in TCGA mesenchymal tumors (COL6A1, FN1; Fig. 6h). Both were previously associated with increased deposition in the ECM by CD133 + glioma cancer stem cells relative to differentiated glioblastoma cells [67] and COL6A1 has been observed in perivascular and PAN regions in patient samples⁶⁸. Our results now place these genes in relation to each other and to other genes within modular networks that constitute prognostic invasion programs.

Discussion

In this study, we used global transcriptional profiling and unsupervised reference-free deconvolution to perform de novo discovery of cell types and states within the GBM ecosystem. We collated a comprehensive and in-depth catalogue of 15 tumor cell programs within the spatiotemporal context of 90 mouse brain cell types, activities, and anatomic structures. The xenograft platform was critical in providing the resolution needed to study the invasive front, where tumor cells were a small minority of the transcriptional output per spot. It has long been recognized that gliomas display histologically distinct patterns of growth relative to normal brain structures. As reported first by Scherer in 1938 [56], these include close interactions of glioma cells with neurons, close interactions with vessels, invasion along white matter tracts, and subpial accumulation. Our findings provide a glimpse into the underlying mechanisms of tumor-microenvironmental cell interactions that support those routes of migration. Importantly, while some of the genes highlighted here had been individually studied in previous work, our deconvolution approach was able to frame these within protein–protein interaction modules that relate expression programs to in vivo phenotypes. In particular, genes serving as program network hubs were highly prognostic, stratifying patients across transcriptional subtypes. Whether effective targeting of these prognostic programs could be achieved through perturbation of single or multiple hubs will require future functional testing in vivo to ensure that dependencies between invasion programs and invasion routes are faithfully maintained.

In addition to human tumor cell programs, we provide a spatiotemporal breakdown of the glioblastoma TME, encompassing multiple cell types and states. We anticipate that this catalogue will represent a tractable system for gauging the potential impact of rational therapies. For instance, both MDM programs (m40_MDM1, m84_MDM2) showed high program scores for Lilrb4a—a gene with established roles in promoting immunosuppression and the target of immunomodulatory therapies currently under development [69]. Similarly, Gpnmb was highly scoring in both MDM programs, has been implicated in proneural to mesenchymal state transitions in GBM [70], and is another promising TME-specific target in GBM. Our data suggest that targeting of these molecules could effectively impact MDMs but not microglial or monocyte programs, where these genes do not significantly contribute to program identity. More generally, coupling spatial profiling analysis as described here with in vivo preclinical assessment of targeted therapies toward either tumor or TME components would facilitate an ecosystem-level understanding of the immediate and long-term consequences of such perturbations. This would include the identification of compensatory programs and build toward the design of combination therapies with improved efficacy.

Our study was in part limited by the impact of cohort composition on program identification using unsupervised deconvolution. The cNMF programs identified here represented the most coherently co-varying signals in the data, and as such, the inclusion of more samples or of distinct biological conditions would likely increase the number of meaningful programs. For instance, the inclusion of xenografts treated with standard or rational therapies should reveal contextual tumor and TME adaptations to therapy. Understanding these relationships will be pivotal in designing and preclinical testing of more effective rational and combination therapies.

Conclusion

We present a compendium of gene expression programs that capture the landscape of tumor and TME interactions in diffuse glioblastoma. This study provides a basis for developing rational strategies that focus on targeting specific aspects of tumor-TME axis. We anticipate our analytic approach and results can facilitate an ecosystem-level understanding of the immediate and long-term consequences of such perturbations in preclinical models, including the identification of compensatory programs that can inform improved combinatorial therapies.

Methods

Intracranial BTIC models

Xenograft samples were generated from brain tumor-initiating cells (BTICs) established from patients with primary and recurrent GBM [30] surgical material available through our institution’s Clark Smith Tumor Biobank. Lines were maintained as previously described [31] prior to intracranial implantation into 6- to 8-week-old female SCID mice [7, 71]. Some BTIC lines were generated from patient surgical tissue collected from the center of the tumor (x), the highly vascular contrast-enhancing regions (y), or the infiltrating edge (z) (BT143x, BT143y, BT143z, BT238x, BT238z). Mice were housed in groups of three to five and maintained on a 12-h light/dark schedule with a temperature of 22°C ± 1°C and a relative humidity of 50 ± 5% and provided food and water ad libitum. All animal procedures were reviewed and approved by the University of Calgary Animal Care Committee (Animal Protocol #AC22-0053).

Visium library preparation and sequencing

The brain was removed and dissected into 2-mm coronal slabs. Fresh tissues were embedded in Tissue Tek OCT compound (Fisher Scientific 14–373-65) and snap-frozen in a chilled isopentane and dry ice bath; 10-µm cryosections were mounted on barcoded Visium slides (10 × Genomics), and libraries prepared with the 10 × Visium Spatial Gene Expression kit per manufacturer’s protocol. Briefly, sections were fixed in methanol, stained with H&E, and scanned on an EVOS FL Auto Imaging System (Thermo Scientific) using a 10 × objective. Permeabilization, reverse transcription, second strand synthesis, denaturation, and cDNA synthesis were performed as per protocol. Cycle number determination for cDNA amplification was performed on a BioRad Realtime qPCR system using KAPA SYBR FAST qPCR Master Mix (Roche, KK4600). cDNA QC and qualification was performed on an Agilent 2100 Bioanalyzer with Agilent High sensitivity DNA chips (Agilent 5067–4626). After enzymatic fragmentation and double size selection using SPRIselect reagent (Beckman Coulter, B23318), unique indexes and P5 and P7 Illumina primers were added to the libraries using Dual Index Kit TT Set A (PN-1000215). Libraries were sequenced on an Illumina NextSeq500/550 instrument using paired-end sequencing with the following parameters: 28 cycles for Read1 and 90 cycles for Read, 10–10 cycles for index i7 and index i5, loading concentration 1.8 pM on NextSeq 500/550 High Output Kit v2.5 150cycle (Illumina, 20024904).

Spatial transcriptomics data preprocessing

Illumina sequencing base call data (BCL) was converted to FASTQ files using bcl2fastq (SpaceRanger v1.3.1). Using 10 × Genomics SpaceRanger software (v1.3.1), the resulting FASTQ files were mapped to a hybrid genome reference sequence (GRCh38—mm10-2020-A) created by combining the human reference genome (GENCODE v32/Ensembl 98) and mouse reference (GENCODE vM23/Ensembl 98) genome. Data was aligned with STAR v2.7.2a and mapped to spatial coordinates using the spatial barcode information in SpaceRanger (default parameters). Additional details regarding species specificity of alignments are provided in Additional File 2.

Samples were aggregated using SpaceRanger Aggr. All reads outside the tissue region were removed in the SpaceRanger pipeline. The resulting filtered matrix output is used for subsequent analysis. This matrix consisted of human and mouse genes from 23 samples. The R package Seurat (v4.0.0) [57] was used to further process this, removing genes with expression in less than three spots. Spots with low number of detected genes (< 200) were excluded. Expression of human and mouse genes was normalized together by dividing expression in each spot by the total number of transcripts and multiplied by 10,000, followed by a natural-log transformation.

Publicly available scRNA-seq, stRNA-seq, and STARmap datasets

Multiple external datasets were downloaded and factorized to yield comparable gene expression programs based on cNMF. The Kleinman [33] scRNAseq data included five normal forebrains (developmental timepoints E12.5, E15.5, P0, P3, and P6), the hindbrain (E12.5), and pons (E15.5, P0, P3, P6). Data was downloaded as CellRanger outputs and cell annotations (GSE133531 [81]). Cells in each sample were filtered based on quality control metrics previously described, using the R package Seurat (v4.0.0). Libraries were scaled to 10,000 UMIs per cell and natural log normalized (Seurat::NormalizeData). The Bayraktar [28] data included matched snRNA-seq and stRNA-seq from adjacent brain sections from six adult mice. Data were downloaded from ArrayExpression, including Visium stRNA-seq (E-MTAB-1114 [82]; Space Ranger output) and snRNAseq (E-MATB-11115 [83]; Cell Ranger output and cell annotations). For Visium profiles, genes with expression in less than three spots were discarded. snRNA-seq data was subjected to quality control criteria described in the corresponding publication. Data were scaled to 10,000 UMIs per cell and natural log normalized (Seurat::NormalizeData). The Movahedi [36] study included CD45 + immune cells from orthotopic GL261 tumors in three adult mice. CellRanger outputs and cell annotations were downloaded from https://www.brainimmuneatlas.org [84]. Outlier cells and low-abundance genes were removed based on the workflow previously described. Raw counts were then scaled to 10,000 UMIs per cell and natural log normalized (Seurat::NormalizeData) (v4.0.0). The Senger [7] study included CD45 + immune cells from orthotopic xenograft models established from patient-derived BTICs. Cell ranger outputs and cell annotations were downloaded (GSE153487 [85]). Data was filtered using quality control metrics previously described and normalized, using R package Seurat (v4.0.0). The Pugh [5] study included scRNA-seq data from glioblastoma stem cells (GSC) (26 patient GSC cultures) and from tumors (7 patients). Raw and normalized gene expression matrices and cell annotations were downloaded from Broad Institute Single-Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP503 [86]). The Heiland [17, 30] study included Visium stRNA-seq profiles of 28 samples from 20 patients. Data from each sample was downloaded in the form of SPATA objects [87], converted to Seurat objects, and processed with guidelines described in the corresponding publication. The TFRI study included bulk RNA-seq data from 56 patients, including 44 tumor samples, 61 BTIC cell lines, and 13 xenograft samples, including both longitudinal and multiregional samples. Fastq files were downloaded (EGAS00001002709 [88]) and aligned (STAR v2.9.9a [66], using parameters: –runThreadN 16 –outSAMtype BAM SortedByCoordinate –quantMode TranscriptomeSAM GeneCounts –outFilterType BySJout –outFilterMatchNminOverLread 0 –outFilterScoreMinOverLread 0 –outSAMstrandField intronMotif –twopassMode Basic). Human reference genome (GENCODE v32/Ensembl 98) was used for mapping tumor and BTIC samples, while a hybrid genome reference (GRCh38—mm10-2020-A) was used for xenografts. The STARmap PLUS data from Hailing Shi et al. [67] encompasses in situ gene expression profiles of 1022 genes ~ 200 nm resolution, mapping to 1.09 million cells across 20 adult mouse brain and spinal cord tissue slices. By integrating with previously published scRNA-seq data and computational registration to the Allen Mouse Brain Common Coordinate Framework (CCFv3), the authors generate a comprehensive CNS spatial atlas of 230 of molecular cell types in 106 molecular tissue regions. In addition, imputed expression profiles of 11,844 genes across all cells in the cohort are provided, creating a transcriptome-wide spatial single-cell atlas. We utilize the STARmap PLUS data as ground truth to perform a systematic comparison of the spatial distribution of marker genes and cellular composition of structures in our xenografts. We downloaded processed STARmap PLUS profiles of a coronal tissue slice well3_5 from the Single-Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP1830/spatial-atlas-of-molecular-cell-types-and-aav-accessibility-across-the-whole-mouse-brain [89]). This includes the imputed gene expression matrix (11,844 genes × 47,607 cells) and metadata of coordinates, molecular cell type, and tissue regions per cell. In addition, gene markers for tissue regions and cell types were obtained from Supp. Tab 5 and Supp. Tab 6 from the reference study, respectively. We also employ the quantifications of cell types across tissue regions in the STARmap PLUS data to correlate gene expression patterns within mouse programs in our xenograft cohort profiled with the 10 × Visium technology (see Additional File 2).

Admixture calculation

We calculated the ratio of UMI counts from human genes relative to the total UMI count per spot in order to distinguish the tumor from non-neoplastic cells in xenografts. This ratio (human–mouse admixture) represented the transcriptional contribution of tumor cells relative to the total transcriptional output in each spot. Admixture was used to categorize spots based on tumor density: high (admixture 80–100%; D4), moderate (50–80%; D3), low (20–50%; D2), and very low (5–20%; D1). Tumor-free mouse brain was comprised of spots with admixture < 5% (D0).

Selection of over-dispersed (OD) genes

To select for features that are informative of identity and/or activity states in latent space, we selected genes significantly over-dispersed across all spots or samples within each dataset. For the newly generated xenograft stRNA-seq cohort and the reference sc/snRNA-seq and stRNA-seq datasets, we selected genes detected in more than 1% but less than 100% of spots or single cells that passed quality control thresholds as described above. Next, genes with higher-than-expected expression variance across the data were selected using a general additive model with a basis of 5 and an adjusted p-value cutoff of 1e-4. Importantly, in the case of the xenograft models, feature selection was performed individually for mouse and human genes across filtered spots. The calculations were based on the restrictCorpus() function from the R package STDeconvolve (v1.7.0) [26].

Matrix factorization and rank selection

We implemented consensus non-negative matrix factorization (cNMF v1.4) [32], a meta-analysis approach of matrix decomposition to independently infer gene expression programs. The command line version of cNMF was run separately for each dataset. In the xenograft models, the pipeline was run separately for mouse and human gene expression data. Run parameters were first prepared by providing the raw count matrix and a list of over-dispersed genes to be used for the factorization steps, all in tab-delimited text file formats to cnmf prepare. Prior to this step, spots with no expression of the OD genes were excluded. The number of factors (K) was set to range from 2 to 100 factors. The frobenius loss function was used. Other parameters included –n-iter 200, –total-workers 1, –seed 123456, –densify False. Default parameters were used in the next series of steps, cnmf factorize and cnmf combine, which factorize and merge results from 200 iterations. The final step, cnmf consensus, used to obtain consensus estimates of gene expression programs was performed for each value of K.

We selected an optimal rank based on the trade-off between stability and error, as well as the biological interpretability of the factors. This included the anticipated number of cell types and states within a dataset and manual review of spatial program profiles (in the case of stRNA-seq), as well as enrichment of previously annotated cell types in reference scRNA-seq cohorts. For the selected representative rank in each dataset, the program usage matrix was normalized such that usage values per spot/cell sum to 1 for downstream analyses. In the xenograft models, we selected a rank of 15 for human and 90 for mouse (brain + TME). To quantify the signal of human and mouse programs relative to the total signal, human program usages per spot were multiplied by the admixture ratio. Mouse program usages per spot were multiplied by 1-admixture ratio. The admix-scaled usage matrices were used for further analyses.

In reference datasets, we did not select an optimal rank, as these were used to assess the maximum concordance of the xenograft programs with programs found in each external cohort. Thus, we calculated all pairwise correlations between each xenograft program and each program (from all ranks) in the reference datasets to identify the best match (i.e., a reference program with maximum correlation to the xenograft cohort programs). Further, we calculated the proportion of patients in each external cohort that had usage of the maximally correlated program above a specific threshold (< 0.3, (0.3–0.4), (0.4–0.5), and > 0.5).

Annotation of gene expression programs

Gene expression programs identified in the xenograft models were annotated based on cell-type and cell-state gene signatures derived from multiple studies of the mouse brain, GBM, and the GBM microenvironment. The sources and the gene markers used are listed in Additional file 1: Table S2a. We used three methods: (1) calculation of marker gene scores, (2) pathway enrichment analysis using g:Profiler [68], and (3) assessment of spatial and molecular concordance with the Allen Brain Atlas Common Coordinate Framework (CCFv3).

1. Marker gene scores: we computed marker-gene scores for each program, defined as the rank weighted overlap of top 50 genes in the program with the gene set being queried. Specifically, for a given query gene set (q), with g number of genes, the intersection with the top 50 genes in a program (p) was derived. If n is the number of overlapping genes (gene₁, gene₂, …, gene_n), the marker gene score MGS_qp is then calculated as:

$${\mathrm{MGS}}_{qp}=\left[\mathrm n\;\mathrm x\;{\textstyle\sum_{i=0}^n}1/\mathrm{rank}\left(\mathrm{gene};\right)\right]/\log10\left(\mathrm g\right),$$

where the rank of a gene is obtained by rank ordering the genes in each cNMF program using the gene scores (i.e., the highest scoring gene in a program is rank 1) and g normalized for the magnitude of the gene set. This method resulted in unambiguous and strong matches to cell-type associated programs and in poor scores for programs representing cell states/activities and anatomic structures not well represented in the reference marker gene sets.

For instance, the astrocyte-like cell-state signature (AC_Neftel) from the Neftel et al. (2019) study comprises 39 markers. In h2_AC, 14 of the top 50 genes intersect this 39-gene signature. The genes and their ranks (in brackets) are AGT (6), CST3 (7), NDRG2 (9), SPARCL1 (13), EDNRB (18), PRCP (20), MT3 (21), GFAP (24), GATM (25), ATP1A2 (30), AQP4 (32), CLU (35), BCAN (38), and SLC1A3 (50). The sum of the inverse ranks (0. 872) is multiplied by the intersection size (14) and divided by the query size (log10(39)) to generate a score of 6.1. All variables and resulting scores for each marker list and program are available in Additional file 1: Table S3b.

2. Pathway enrichment analysis: to annotate cell activities, we used the top 1000 genes in each program to identify highly enriched pathways from various sources including GO:BP, GO:MF, KEGG, and REACTOME using g:Profiler. Terms with a size > 10 and < 2000 were included and adjusted p-value threshold was set to 0.05. Results were summarized as a heatmap of adjusted p-values of terms ordered by significance in each program.

3. Spatial and molecular concordance with the CCFv3: to annotate anatomic structures, a neuropathologist (JAC) manually reviewed the concordance of each program’s spatial usage against the template mouse brain (https://atlas.brain-map.org/). Briefly, we first manually extracted a slice image from the Allen CCFv3 that best corresponded to a Visium section (based on the layout of the major white matter tracts, ventricles, etc.). Each program’s spatial profile was then compared to the isocortical and subcortical areas present in the CCFv3, as well as to ventricle structures and fiber tracts—any matches were annotated and incorporated in the program name (Additional file 1: Table S2b).

We validated a subset of program-to-structure assignments made in this way using molecular markers. In the first strategy, we used the “Differential Search” function in the ABA to identify structure-specific genes with in situ hybridization (ISH) data for CCFv3 and intersected these with the genes in each cNMF structure program. We show high visual concordance between the ISH-based marker genes and the spatial cNMF programs and top-scoring genes (see Additional File 2). In the second strategy, we used a single-cell resolution STARmap spatial transcriptome dataset of mouse brain previously registered to the CCFv3 [35]. We intersected the marker genes of anatomic structures from the STARmap data with our cNMF structure programs to establish both molecular and spatial concordance (see Additional File 2 for details).

Spatial concordance of programs

To infer short-range spatial overlap between program usage among spots for each program, we calculated the proportion of spots that also contribute to a second query program. A usage threshold of > 0.1 was needed for a human tumor program to be considered present. Spots with mouse program usage of > 0.05 were considered to be positive for that program. The lower threshold used for mouse programs reflects the higher diversity of programs per spot in the normal brain relative to human tumor programs.

SCENIC

Active regulons in human and mouse programs were identified using the R package SCENIC (v1.1.2.1) [69] with default parameters. The matrix of z-scores of genes per program obtained from cNMF was used as input to GENIE3 to infer co-expression modules, where each module consisted of a transcription factor (TF) and its predicted targets based on co-expression. Next, using RcisTarget (v1.6.0) [69], each module (regulon) was pruned to include only targets for which the motif of the TF was enriched. Finally, AUCell (v1.8.0) was used to calculate regulon activity scores per program as the area under the curve (AUC).

CNV analysis

Copy number changes in individual spots were identified using the R package inferCNV [70] (v1.3.3). To obtain a robust signal, raw gene counts from each sample were subjected to a spatially aware smoothing method SPCS [71] using default parameters. As a normal reference dataset, we used a cortex section of the human brain from 10 × Genomics profiled with the Visium platform (https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Human_Brain_Section_1). The data were then passed to inferCNV for copy number variant inference, separately for each cell line and the following parameters were used: “denoise”, cutoff = 0.1, sd_amplifier = 2.5, scale_data = TRUE, analysis_mode = “subclusters,” window_length = 201, num_ref_groups = 15. All remaining parameters were set to default. Mitochondrial genes were excluded from this analysis. Individual CNV scores were averaged across clusters to visualize unique transcriptional tumor clones. To further assess how inferred CNVs may impact tumor biology, clones with variable CNVs across chromosome arms were identified by assigning genes to regions of gain or loss using inferCNV score cutoffs of 2.5 times the standard deviation below and above the mean. Spots harboring these clones were then used to assess spatial concordance with tumor programs.

GBM cell-state signatures

Module score and change across tumor density: gene signatures for GBM cell states (astrocyte-like (AC-like), oligodendrocyte-like (OPC-like), neural progenitor-like 1/2 (NPC-like), and mesenchymal-like 1/2 (MES-like)) were obtained from Neftel et al [4]. Spots with admixture ratio of ≥ 0.05 were selected for the analysis of tumor cell states. Spots were scored for cell states using the AddModuleScore() function in Seurat (v4.0.0). NPC1/2 and MES1/2 scores were averaged to represent a score for NPC-like and MES-like states, respectively. To analyze changes in GBM cell states as the density of tumor cells in the xenograft models increase (D1-D4), each spot was assigned to the cell state with the highest module score. Then significance of changes in the proportion of spots belonging to each state across the four groups was assessed using a Wilcoxon signed rank-sum test and visualized as boxplots.

Cell-state transition: to deduce plasticity among GBM cell states for each human transcriptional program, we calculated a transition score representing a change from the baseline of the sample to an end state observed in a given spot. First, a baseline state was established per sample. This was done by calculating the proportional overlap of highly contributing spots in a program (usage > 0.1) with spots assigned a GBM cell state based on their highest module score. The overlap metrics are assessed per sample and further by tumor density groups. The baseline state of a sample was then defined as that which is observed as the majority signal in end-point sample of patient (heatmaps and manual review). Next, the transition score for spots of a given baseline state at a given timepoint (early, mid, or late) and tumor-density range (D1–D4) was derived in two steps: (1) a matrix of proportional overlap between tumor programs and GBM cell states was created as described above and (2) for each program, the transition score was calculated as the relative contribution of a cell state over the baseline state and visualized as barplots, grouped by tumor stage and tumor density.

Cell–cell communication

Inference of ligand-receptor (LR) communication was based on the R package CellChat (v1.6.1) [47]. Given the expected cross-species communication in the xenograft models, ligand, and receptor information available in species-specific CellChatDB were used to create two additional databases, one enumerating human-to-mouse interactions and the other, mouse-to-human. To enable assessment of all directional interactions, the four databases (comprising tumor-to-tumor, TME-to-TME, tumor-to-TME and TME-to-tumor interactions) were re-built into a single database for further analysis. The log-normalized matrix consisting of mouse and human gene expression was used to create a CellChat object. Interactions were inferred across (1) D1–D4 groups and (2) white matter (WM), caudoputamen (parenchymal), m30, m52, m59, and m86 (vascular) mouse programs only within the D1 group. Spots within D1 (admixture: 0.05–0.2) were assigned to programs in (2) as follows: spots with usage > 0.1 in m19 or m71 were assigned to WM, > 0.1 in m1 or m73 to parenchymal and > 0.01 in m30, m52, m59, or m86 to the respective vascular programs. In case of spots with more than one of these programs, label assignment between the groups was prioritized as WM > vascular > parenchymal. In both analyses, overexpressed genes and interactions were identified using default parameters. Next, computeCommunProb() was used to derive significant (p.adj < 0.05) communications with 10% truncated mean for calculating average gene expression per spot. A data frame of all communications inferred was obtained using subsetCommunication() and used for downstream analysis. Chord diagrams of LR interactions per pathway were visualized using netVisual_chord_cell(). In (1), pathways were categorized into eight groups LR1-LR8 based on combinations of the four interaction categories (TME__TME, TME__Tumor, Tumor__TME, and Tumor__Tumor). Pathway enrichment analysis with GO:BP was performed by querying ligands and receptors of significant interactions in all pathways using g:Profiler. Only terms with a size > 10 and < 2000 were included and adjusted p-value threshold was set to 0.05. Results were summarized as a heatmap of adjusted p-values of terms ordered by significance in each LR group.

STRINGdb networks and selection of hub genes

Selection of hub genes

To build a network of protein–protein interactions per tumor program, we queried a set of “top-scoring” genes in STRING database (v2.6.5) [72] to obtain known and predicted interactions. Top-scoring genes of a program were derived by plotting the distribution of its gene-scores produced by cNMF, followed by the selection of genes that comprise the final component of Gaussian mixture model clusters. This was performed with default parameters of densityMclust() from the R package mclust (v6.0.0). Known + Predicted interactions obtained medium confidence levels from STRINGdb were plotted as a network with a Fruchterman–Reingold layout using the R package iGraph (v1.4.2). Nodes represented genes, and node size corresponded to the correlation between the expression of the gene and GEP usage across all spots. Finally, hub genes within the network were defined as nodes with a PageRank score in the 95th percentile. Hub genes selected using the PageRank algorithm were observed to be distributed over multiple “fast-greedy” communities within the networks.

Hub gene characterization

To functionally characterize hub genes identified in PPI networks of highly-scoring genes, we performed an enrichment analysis of GO:BP terms. Enrichment analysis was performed using the R package gProfiler2 (v0.2.2). The hub genes were queried using gost() with default parameters and converted to an enrichResult using the R package enrichplot (v1.14.2). For each program, terms with size > 5 and < 500 and p.adjust < 0.05 were selected and reduced to ~ 5 parent terms using iterative thresholds within go_reduce() from the rutils package (v0.99.2). The data was then visualized as a dot plot with the dot size representing the number of terms reduced per parent term.

Survival analysis

RNA-seq and clinical data for the TCGA-GBM cohort [73] were downloaded from http://xena.ucsc.edu [90]. For a given sample, normalized enrichment score (NES) for the gene set comprising hub genes of each tumor program was calculated using the R package fgsea (v1.20.0) with default parameters. Patients were stratified into three groups based on overall survival time—top and bottom 25%, 33%, and 50% and Kaplan–Meier curves were plotted for each group using R package survminer (v0.4.9) with, p-values determined by log-rank test. As a complementary method, a Student’s t test was also performed between NES of the stratified patients. Kaplan–Meier curves for individual genes were obtained from GEPIA2 [74].

Over-representation analysis

To assess the enrichment of a program across regions of varying tumor cell density, Pearson’s chi-square test was performed using a contingency table encompassing the number spots expressing a program (selected based on their usage) and the remainder of the spots in each category [Brain(D0), Tumor (D1-D4), and four tumor density regions (D1, D2, D3 and D4)]. “Brain” and “Tumor” categories were excluded while testing for the enrichment of tumor programs. A usage threshold of > 0.05 was used to select spots in each TME programs and > 0.1 for tumor and normal brain programs. The chi-squared residuals indicating the observed minus expected number of spots with the usage of a specific program per category (Brain(D0), Tumor (D1–D2), D1, D2, D3, and D4) are then plotted.

Differential expression analysis

To determine differentially expressed genes between regions of low and high tumor density with a program-agnostic approach, we used the R package ALDEx2 (v1.24.0) [52]. First, spots in each sample were categorized as low and high density (“lowHs” and “higHs”) by selecting those in the 15th and 75th percentile of the sample’s admixture values. Then, the ALDEx2 pipeline was implemented per sample using the aldex() function with the following parameters: mc.samples = 256, test = “t”, effect = TRUE, denom = “lvha”, paired.test = FALSE. Genes with BH adjusted p.value < 0.05 and absolute effect size > 1 were selected as differentially expressed in each category. Further, to obtain differentially expressed genes between WM and other groups (Par, m52, m30, m59, m86) within the D1 region, we used FindAllMarkers() with default parameters from R package Seurat (v4.0.0), following filtering (min.cells = 3, min.features = 5, customized nFeature_RNA and nCount_RNA filters per sample) and normalization with NormalizeData() of the selected spots. Genes with adjusted p.value < 0.05 and logfc.threshold > 0.25 were selected as differentially expressed.

Availability of data and materials

Space ranger and cNMF processed files of spatial transcriptomics data generated in this study have been deposited at Datadryad (https://doi.org/10.5061/dryad.wpzgmsbv6) [91]. FASTQ files have been deposited under BioProject ID PRJNA1137934 in the SRA repository [92]. The following publicly available data were obtained as follows: (1) scRNA-seq data of developing mice brain were obtained from GSE133531 [33, 81], (2) snRNA-seq and stRNA-seq data of adult mice brain were obtained from ArrayExpression under accession numbers E-MTAB-1114 and E-MATB-11115 [28, 82, 83], (3) scRNA-seq data from CD45 + immune cells from orthotopic tumors were obtained from https://www.brainimmuneatlas.org [36, 84], (4) scRNA-seq data of CD45 + immune cells from IL33 + and IL33- orthotopic xenograft models were obtained from GSE153487 [7, 85], (5) scRNA-seq data from glioblastoma stem cells (GSC) and GBM tumor cells were obtained Broad Institute Single-Cell Portal https://singlecell.broadinstitute.org/single_cell/study/SCP503 [5, 86], (6) stRNA-seq data of GBM patient samples were obtained from https://themilolab.github.io/SPATA2/articles/spata-v2-spata-data.html [17, 87], (7) raw RNA-seq data of GBM tumors, cell lines and xenografts were obtained from EGA under the accession number EGAS00001002709 [30, 88], and (8) TCGA-GBM RNA-seq data were obtained from http://xena.ucsc.edu [90]. Source data underlying main figures and analyses are included in Addition file 1. Code used to perform analyses presented in this manuscript are deposited in Datadryad (https://doi.org/10.5061/dryad.wpzgmsbv6) [91] and Github ( https://github.com/MorrissyLab/GBM_Xenograft_SpatialTranscriptomics) [93] and is released under an MIT license.

References

Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJB, et al. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. New England J Med. 2005;352(10):987–96.
CAS Google Scholar
Cuddapah VA, Robel S, Watkins S, Sontheimer H. A neurocentric perspective on glioma invasion. Nat Rev Neurosci. 2014;15(7):455–65.
CAS Google Scholar
Wang L, Jung J, Babikir H, Shamardani K, Jain S, Feng X, et al. A single-cell atlas of glioblastoma evolution under therapy reveals cell-intrinsic and cell-extrinsic therapeutic targets. Nat Cancer. 2022;3(12):1534–52.
CAS Google Scholar
Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell. 2019;178(4):835-849.e21.
CAS Google Scholar
Richards LM, Whitley OKN, MacLeod G, Cavalli FMG, Coutinho FJ, Jaramillo JE, et al. Gradient of developmental and injury response transcriptional states defines functional vulnerabilities underpinning glioblastoma heterogeneity. Nat Cancer. 2021;2(2):157–73.
CAS Google Scholar
Varn FS, Johnson KC, Martinek J, Huse JT, Nasrallah MP, Wesseling P, et al. Glioma progression is shaped by genetic evolution and microenvironment interactions. Cell. 2022;185(12):2184-2199.e16.
CAS Google Scholar
De Boeck A, Ahn BY, D’Mello C, Lun X, Menon SV, Alshehri MM, et al. Glioma-derived IL-33 orchestrates an inflammatory brain tumor microenvironment that accelerates glioma progression. Nat Commun. 2020;11(1):4997.
Google Scholar
Wang J, Xu SL, Duan JJ, Yi L, Guo YF, Shi Y, et al. Invasion of white matter tracts by glioma stem cells is regulated by a NOTCH1–SOX2 positive-feedback loop. Nat Neurosci. 2019;22(1):91–105.
CAS Google Scholar
Jung E, Osswald M, Ratliff M, Dogan H, Xie R, Weil S, et al. Tumor cell plasticity, heterogeneity, and resistance in crucial microenvironmental niches in glioma. Nat Commun. 2021;12(1):1014.
CAS Google Scholar
Sarkar S, Mirzaei R, Zemp FJ, Wei W, Senger DL, Robbins SM, et al. Activation of NOTCH signaling by Tenascin-C promotes growth of human brain tumor-initiating cells. Cancer Res. 2017;77(12):3231–43.
CAS Google Scholar
Comba A, Faisal SM, Dunn PJ, Argento AE, Hollon TC, Al-Holou WN, et al. Spatiotemporal analysis of glioma heterogeneity reveals COL1A1 as an actionable target to disrupt tumor progression. Nat Commun. 2022;13(1):3606.
CAS Google Scholar
Bastola S, Pavlyukov MS, Yamashita D, Ghosh S, Cho H, Kagaya N, et al. Glioma-initiating cells at tumor edge gain signals from tumor core cells to promote their malignancy. Nat Commun. 2020;11(1):4660.
CAS Google Scholar
Minata M, Audia A, Shi J, Lu S, Bernstock J, Pavlyukov MS, et al. Phenotypic plasticity of invasive edge glioma stem-like cells in response to ionizing radiation. Cell Rep. 2019;26(7):1893-1905.e7.
CAS Google Scholar
Venkatesh HS, Tam LT, Woo PJ, Lennon J, Nagaraja S, Gillespie SM, et al. Targeting neuronal activity-regulated neuroligin-3 dependency in high-grade glioma. Nature. 2017;549(7673):533–7.
Google Scholar
Bressan D, Battistoni G, Hannon GJ. The dawn of spatial omics. Science. 2023;381:eabq4964.
CAS Google Scholar
Greenwald AC, Darnell NG, Hoefflin R, Simkin D, Gonzalez-Castro LN, Mount C, et al. Integrative spatial analysis reveals a multi-layered organization of glioblastoma. bioRxiv. 2023. https://doi.org/10.1016/j.cell.2024.03.029.
Ravi VM, Will P, Kueckelhaus J, Sun N, Joseph K, Salié H, et al. Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma. Cancer Cell. 2022;40(6):639-655.e13.
CAS Google Scholar
Zheng Y, Carrillo-Perez F, Pizurica M, Heiland DH, Gevaert O. Spatial cellular architecture predicts prognosis in glioblastoma. Nat Commun. 2023;14(1):4122.
CAS Google Scholar
Gangoso E, Southgate B, Bradley L, Rus S, Galvez-Cancino F, McGivern N, et al. Glioblastomas acquire myeloid-affiliated transcriptional programs via epigenetic immunoediting to elicit immune evasion. Cell. 2021;184(9):2454-2470.e26.
CAS Google Scholar
Ren Y, Huang Z, Zhou L, Xiao P, Song J, He P, et al. Spatial transcriptomics reveals niche-specific enrichment and vulnerabilities of radial glial stem-like cells in malignant gliomas. Nat Commun. 2023;14(1):1028.
CAS Google Scholar
Andersson A, Bergenstråhle J, Asp M, Bergenstråhle L, Jurek A, Fernández Navarro J, et al. Single-cell and spatial transcriptomics enables probabilistic inference of cell type topography. Commun Biol. 2020;3(1):565.
Google Scholar
Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol. 2022;40(4):517-526.
Biancalani T, Scalia G, Buffoni L, Avasthi R, Lu Z, Sanger A, et al. Deep learning and alignment of spatially resolved single-cell transcriptomes with Tangram. Nat Methods. 2021;18(11):1352–62.
Google Scholar
Ma Y, Zhou X. Spatially informed cell-type deconvolution for spatial transcriptomics. Nat Biotechnol. 2022;40(9):1349–59.
CAS Google Scholar
Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H. SPOTlight: Seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes. Nucleic Acids Res. 2021;49(9):e50.
CAS Google Scholar
Miller BF, Huang F, Atta L, Sahoo A, Fan J. Reference-free cell type deconvolution of multi-cellular pixel-resolution spatially resolved transcriptomics data. Nat Commun. 2022;13(1):2339.
CAS Google Scholar
Rodriques SG, Stickels RR, Goeva A, Martin CA, Murray E, Vanderburg CR, et al. Slide-seq: a scalable technology for measuring genome-wide expression at high spatial resolution. Science (1979). 2019;363(6434):1463–7.
CAS Google Scholar
Kleshchevnikov V, Shmatko A, Dann E, Aivazidis A, King HW, Li T, et al. Cell 2location maps fine-grained cell types in spatial transcriptomics. Nat Biotechnol. 2022;40(5):661–71.
CAS Google Scholar
Coleman K, Hu J, Schroeder A, Lee EB, Li M. SpaDecon: cell-type deconvolution in spatial transcriptomics with semi-supervised learning. Commun Biol. 2023;6(1):378.
Google Scholar
Shen Y, Grisdale CJ, Islam SA, Bose P, Lever J, Zhao EY, et al. Comprehensive genomic profiling of glioblastoma tumors, BTICs, and xenografts reveals stability and adaptation to growth environments. Proc Natl Acad Sci U S A. 2019;116(38):19098–108.
CAS Google Scholar
Kelly JJP, Stechishin O, Chojnacki A, Lun X, Sun B, Senger DL, et al. Proliferation of human glioblastoma stem cells occurs independently of exogenous mitogens. Stem Cells. 2009;27(8):1722–33.
CAS Google Scholar
Kotliar D, Veres A, Nagy MA, Tabrizi S, Hodis E, Melton DA, et al. Identifying gene expression programs of cell-type identity and cellular activity with single-cell RNA-Seq. Elife. 2019;8:8.
Google Scholar
Zhang M, Pan X, Jung W, Halpern AR, Eichhorn SW, Lei Z, et al. Molecularly defined and spatially resolved cell atlas of the whole mouse brain. Nature. 2023;624:7991.
Google Scholar
Yao Z, van Velthoven CTJ, Kunst M, Zhang M, McMillen D, Lee C, et al. A high-resolution transcriptomic and spatial atlas of cell types in the whole mouse brain. Nature. 2023;624(7991):317–32.
CAS Google Scholar
Jessa S, Blanchet-Cohen A, Krug B, Vladoiu M, Coutelier M, Faury D, et al. Stalled developmental programs at the root of pediatric brain tumors. Nat Genet. 2019;51(12):1702–13.
CAS Google Scholar
Lein ES, Hawrylycz MJ, Ao N, Ayres M, Bensinger A, Bernard A, et al. Genome-wide atlas of gene expression in the adult mouse brain. Nature. 2007;445(7124):168–76.
CAS Google Scholar
Wang Q, Ding SL, Li Y, Royall J, Feng D, Lesnar P, et al. The allen mouse brain common coordinate framework: a 3D reference atlas. Cell. 2020;181(4):936-953.e20.
CAS Google Scholar
Pombo Antunes AR, Scheyltjens I, Lodi F, Messiaen J, Antoranz A, Duerinck J, et al. Single-cell profiling of myeloid cells in glioblastoma across species and disease stage reveals macrophage competition and specialization. Nat Neurosci. 2021;24(4):595–610.
CAS Google Scholar
Couturier CP, Ayyadhury S, Le PU, Nadaf J, Monlong J, Riva G, et al. Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy. Nat Commun. 2020;11(1):3406.
CAS Google Scholar
Puchalski RB, Shah N, Miller J, Dalley R, Nomura SR, Yoon JG, et al. An anatomic transcriptional atlas of human glioblastoma. Science (1979). 2018;360(6389):660–3.
CAS Google Scholar
Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science (1979). 2014;344(6190):1396–401.
CAS Google Scholar
Dusart P, Hallström BM, Renné T, Odeberg J, Uhlén M, Butler LM. A systems-based map of human brain cell-type enriched genes and malignancy-associated endothelial changes. Cell Rep. 2019;29(6):1690-1706.e4.
CAS Google Scholar
Darmanis S, Sloan SA, Croote D, Mignardi M, Chernikova S, Samghababi P, et al. Single-cell RNA-Seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep. 2017;21(5):1399–410.
CAS Google Scholar
Verhaak RGW, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010;17(1):98–110.
CAS Google Scholar
Doetsch F, García-Verdugo JM, Alvarez-Buylla A. Cellular composition and three-dimensional organization of the subventricular germinal zone in the adult mammalian brain. J Neurosci. 1997;17(13):5046–61.
CAS Google Scholar
Monteiro AR, Hill R, Pilkington GJ, Madureira PA. The role of hypoxia in glioblastoma invasion. Cells. 2017;6(4):45.
Google Scholar
Venkataramani V, Yang Y, Schubert MC, Reyhan E, Tetzlaff SK, Wißmann N, et al. Glioblastoma hijacks neuronal mechanisms for brain invasion. Cell. 2022;185(16):2899-2917.e31.
CAS Google Scholar
Das S, Li Z, Noori A, Hyman BT, Serrano-Pozo A. Meta-analysis of mouse transcriptomic studies supports a context-dependent astrocyte reaction in acute CNS injury versus neurodegeneration. J Neuroinflammation. 2020;17(1):227.
CAS Google Scholar
Hu X, Deng Q, Ma L, Li Q, Chen Y, Liao Y, et al. Meningeal lymphatic vessels regulate brain tumor drainage and immunity. Cell Res. 2020;30:3.
Google Scholar
Gragnano F, Sperlongano S, Golia E, Natale F, Bianchi R, Crisci M, et al. The role of von Willebrand factor in vascular inflammation: from pathogenesis to targeted therapy. Mediators Inflamm. 2017;2017:5620314.
Google Scholar
Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12(1):1088.
CAS Google Scholar
Angenendt L, Wöste M, Mikesch JH, Arteaga MF, Angenendt A, Sandmann S, et al. Calcitonin receptor-like (CALCRL) is a marker of stemness and an independent predictor of outcome in pediatric AML. Blood Adv. 2021;5(21):4413–21.
CAS Google Scholar
Heddleston JM, Li Z, McLendon RE, Hjelmeland AB, Rich JN. The hypoxic microenvironment maintains glioblastoma stem cells and promotes reprogramming towards a cancer stem cell phenotype. Cell Cycle. 2009;8(20):3274–84.
CAS Google Scholar
Li Z, Bao S, Wu Q, Wang H, Eyler C, Sathornsumetee S, et al. Hypoxia-inducible factors regulate tumorigenic capacity of glioma stem cells. Cancer Cell. 2009;15(6):501–13.
CAS Google Scholar
Gu S, Shu L, Zhou L, Wang Y, Xue H, Jin L, et al. Interfering with CALCRL expression inhibits glioma proliferation, promotes apoptosis, and predicts prognosis in low-grade gliomas. Ann Transl Med. 2022;10(23):1277.
CAS Google Scholar
Scherer HJ. Structural development in gliomas. Am J Cancer. 1938;34(3):333–51.
Google Scholar
Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573–87.
CAS Google Scholar
Fernandes AD, Reid JNS, Macklaim JM, McMurrough TA, Edgell DR, Gloor GB. Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome. 2014;2(1):15.
Google Scholar
Venkatesh HS, Johung TB, Caretti V, Noll A, Tang Y, Nagaraja S, et al. Neuronal activity promotes glioma growth through neuroligin-3 secretion. Cell. 2015;161(4):803–16.
CAS Google Scholar
Lustig M, Sakurai T, Grumet M. Nr-CAM promotes neurite outgrowth from peripheral ganglia by a mechanism involving axonin-1 as a neuronal receptor. Dev Biol. 1999;209(2):340–51.
CAS Google Scholar
Sehgal A, Boynton AL, Young RF, Vermeulen SS, Yonemura KS, Kohler EP, et al. Cell adhesion molecule Nr-CAM is over-expressed in human brain tumors. Int J Cancer. 1998;76(4):451–8.
CAS Google Scholar
Schuster A, Klein E, Neirinckx V, Knudsen AM, Fabian C, Hau AC, et al. AN1-type zinc finger protein 3 (ZFAND3) is a transcriptional regulator that drives Glioblastoma invasion. Nat Commun. 2020;11(1):6366.
CAS Google Scholar
Aiken J, Moore JK, Bates EA. TUBA1A mutations identified in lissencephaly patients dominantly disrupt neuronal migration and impair dynein activity. Hum Mol Genet. 2019;28(8):1227–43.
CAS Google Scholar
Belvindrah R, Natarajan K, Shabajee P, Bruel-Jungerman E, Bernard J, Goutierre M, et al. Mutation of the α-tubulin Tuba1a leads to straighter microtubules and perturbs neuronal migration. J Cell Biol. 2017;216(8):2443–61.
CAS Google Scholar
Camby I, Belot N, Lefranc F, Sadeghi N, De Launoit Y, Kaltner H, et al. Galectin-1 modulates human glioblastoma cell migration into the brain through modifications to the actin cytoskeleton and levels of expression of small GTPases. J Neuropathol Exp Neurol. 2002;61(7):585–96.
CAS Google Scholar
Van Woensel M, Mathivet T, Wauthoz N, Rosière R, Garg AD, Agostinis P, et al. Sensitization of glioblastoma tumor micro-environment to chemo- and immunotherapy by Galectin-1 intranasal knock-down strategy. Sci Rep. 2017;7(1):1217.
Google Scholar
Shevchenko V, Arnotskaya N, Pak O, Sharma A, Sharma HS, Khotimchenko Y, et al. Molecular determinants of the interaction between glioblastoma CD133+ cancer stem cells and the extracellular matrix. Int Rev Neurobiol. 2020;151:155–69.
CAS Google Scholar
Turtoi A, Blomme A, Bianchi E, Maris P, Vannozzi R, Naccarato AG, et al. Accessibilome of human glioblastoma: Collagen-VI-alpha-1 is a new target and a marker of poor outcome. J Proteome Res. 2014;13(12):5660–9.
CAS Google Scholar
Sharma N, Atolagbe OT, Ge Z, Allison JP. LILRB4 suppresses immunity in solid tumors and is a potential target for immunotherapy. J Exp Med. 2021;218(7):e20201811.
CAS Google Scholar
Xiong A, Zhang J, Chen Y, Zhang Y, Yang F. Integrated single-cell transcriptomic analyses reveal that GPNMB-high macrophages promote PN-MES transition and impede T cell activation in GBM. EBioMedicine. 2022;83:104239.
CAS Google Scholar
Lun X, Wells JC, Grinshtein N, King JC, Hao X, Dang NH, et al. Disulfiram when combined with copper enhances the therapeutic effects of temozolomide for the treatment of glioblastoma. Clin Cancer Res. 2016;22(15):3860–75.
CAS Google Scholar
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
CAS Google Scholar
Shi H, He Y, Zhou Y, Huang J, Maher K, Wang B, et al. Spatial atlas of the mouse central nervous system at molecular resolution. Nature. 2023;622(7983):552–61.
CAS Google Scholar
Kolberg L, Raudvere U, Kuzmin I, Adler P, Vilo J, Peterson H. G:Profiler-interoperable web service for functional enrichment analysis and gene identifier mapping (2023 update). Nucleic Acids Res. 2023;51(W1):W207–12.
CAS Google Scholar
Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14(11):1083–6.
CAS Google Scholar
Tickle T, Tirosh I, Georgescu C, Brown M, Haas B. inferCNV of the Trinity CTAT Project. 2019.
Liu Y, Wang T, Duggan B, Sharpnack M, Huang K, Zhang J, et al. SPCS: a spatial and pattern combined smoothing method for spatial transcriptomic expression. Brief Bioinform. 2022;23(3):bbac116.
Google Scholar
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43:43(D1).
Google Scholar
Weinstein JN, Collisson EA, Mills GB, Shaw KRM, Ozenberger BA, Ellrott K, et al. The cancer genome atlas pan-cancer analysis project. Nat Genet. 2013;45(10):1113–20.
Google Scholar
Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019;47:47(W1).
Google Scholar
Jessa S, Blanchet-Cohen A, Krug B, Vladoiu M et al. Single-cell atlas of the developing brain to investigate the cellular origins of pediatric brain tumors. Gene Expression Omnibus. 2019. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133531.
Kleshchevnikov V. 10X Visium spatial RNA-seq from adult mouse brain sections paired to single-nucleus RNA-seq. Array Expression. 2022. https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-11114
Kleshchevnikov V. Single-nucleus RNA-seq from adult mouse brain sections paired to 10X Visium spatial RNA-seq. Array Expression. 2022 https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-11114
Pombo Antunes AR, Scheyltjens I, Lodi F, Messiaen J et al. Single-cell profiling of myeloid cells in glioblastoma across species and disease stage reveals macrophage competition and specialization. Gene Expression Omnibus. 2020. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163120
De Boeck A, Ahn BY, D'Mello C, Lun X et al. Single cell profiling of cerebral immune cells in the Interleukin-33 (IL-33) driven tumor microenvironment. Gene Expression Omnibus. 2020. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE153487.
Laura M. Richards, Owen K. N. Whitley, Graham MacLeod, Florence M. G. Cavalli et al. Gradient of developmental and injury response transcriptional states define functional vulnerabilities underpinning glioblastoma heterogeneity. Single Cell Portal. 2021. https://singlecell.broadinstitute.org/single_cell/study/SCP503/gradient-of-developmental-and-injury-reponse-transcriptional-states-define-functional-vulnerabilities-underpinning-glioblastoma-heterogeneity#study-summary.
Ravi V, Will P, Kueckelhaus J, Sun N, Joseph K, et al. Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma. Dryad. 2022. https://doi.org/10.5061/dryad.h70rxwdmj.
Article Google Scholar
Shen Y, Grisdale CJ, Islan SA, Bose P, Lever J, et al. Comprehensive genomic profiling of matched glioblastoma tumours, cell-lines, and xenografts reveals genomic stability and adaptation to disparate growth environments. European Genome-Phenome Archive. 2019. https://ega-archive.org/studies/EGAS00001002709.
Shi H, He Y, Zhou Y, Huang J, Maher K, Wang B, et al. Spatial atlas of molecular cell types and AAV accessibility across the mouse central nervous system. Single Cell Portal. 2023. https://doi.org/10.1038/s41586-023-06569-5.
Goldman, M.J., Craft, B., Hastie, M. et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol (2020) UCSC Xena https://xenabrowser.net/datapages/?cohort=TCGA%20TARGET%20GTEx&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
Manoharan VT, Abdelkareem A, Brown S, Gillmor A, Hall C, Seo H, et al. Spatiotemporal modeling reveals high-resolution invasion states in glioblastoma. Dryad. 2024. https://doi.org/10.5061/dryad.wpzgmsbv6.
Article Google Scholar
Manoharan VT, Abdelkareem A, Brown S, Gillmor A, Hall C, Seo H, et al. 10x Visium spatial RNA-Seq of human GBM cell-line derived murine orthotopic xenograft. Sequence Read Archive, 2024. https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1137934
Manoharan VT, Abdelkareem A, Brown S, Gillmor A, Hall C, Seo H, et al. GBM xenograft – spatial transcriptomics. 2024.https://github.com/MorrissyLab/GBM_Xenograft_SpatialTranscriptomics

Download references

Funding

This research was supported by multiple grants and scholarships. ASM was supported by a Canadian Institutes of Health Research (CIHR) Operating Grant [grant number 438802]. ASM, FZ, DM, and JC were supported by The Alberta Cellular Therapy and Immune Oncology (ACTION) Initiative, funded by the Canadian Cancer Society (CCS) [grant number 2020–707161]. Funding for ACTION was also provided by generous community donors through the Alberta Children’s Hospital Foundation (ACHF). ASM holds a Canada Research Chair (CRC) Tier 2 in Precision Oncology. JC was also supported by a Kids Cancer Care Chair in Pediatric Oncology and philanthropic funding from Charbonneau Cancer Institute. DLS was supported by a Canadian Institutes of Health Research (CIHR) Operating Grant [grant number 419959]. VTM was supported by a Clark H Smith Brain Tumour Centre Graduate Scholarship. AA was supported by an Alberta Children's Hospital Research Institute Graduate Scholarship, an Alberta Innovates Graduate Student Scholarships for Data-Enabled Innovation, and a Clark Smith Brain Tumour Graduate Scholarship. AG was supported by an Alberta Graduate Excellence Scholarship (AGES), a University of Calgary Faculty of Medicine Graduate Council Scholarship, Alberta Innovates Graduate Student Scholarship, and Margaret Rosso Graduate Scholarship in Cancer Research. CH was supported by an Alberta Graduate Excellence Scholarship (AGES), and an Achievers in Medical Science Doctoral Scholarship. SG was supported by a Program for Undergraduate Research Experience (PURE) award. LM was supported by a Canadian Cancer Society Research Training Award—PhD in partnership with The Terry Fox Research Institute (CCS award #707977 /TFRI grant #1151–06). The funders had no role in the study design, data collection, and interpretation or the decision to submit the work for publication. Funding for open access charge: Canadian Institutes for Health Research.

Author information

Donna L. Senger, Jennifer A. Chan, and A. Sorana Morrissy are co-senior authors.

Authors and Affiliations

Department of Biochemistry and Molecular Biology, University of Calgary, Calgary, AB, Canada
Varsha Thoppey Manoharan, Aly Abdelkareem, Gurveer Gill, Samuel Brown, Aaron Gillmor, Courtney Hall, Heewon Seo, Kiran Narta, Sean Grewal, Franz Zemp & A. Sorana Morrissy
Charbonneau Cancer Institute, University of Calgary, Calgary, AB, Canada
Varsha Thoppey Manoharan, Aly Abdelkareem, Gurveer Gill, Samuel Brown, Aaron Gillmor, Courtney Hall, Heewon Seo, Kiran Narta, Sean Grewal, Ngoc Ha Dang, Bo Young Ahn, Kata Osz, Xueqing Lun, Laura Mah, Franz Zemp, Douglas Mahoney, Jennifer A. Chan & A. Sorana Morrissy
Department of Microbiology, Immunology and Infectious Diseases, Cumming School of Medicine, University of Calgary, Calgary, AB, Canada
Laura Mah & Douglas Mahoney
Gerald Bronfman Department of Oncology, McGill University, Montreal, QC, Canada
Donna L. Senger
Lady Davis Institute for Medical Research, Jewish General Hospital, Montreal, QC, Canada
Donna L. Senger
Department of Pathology and Laboratory Medicine, University of Calgary, Calgary, AB, Canada
Jennifer A. Chan
Alberta Children’s Hospital Research Institute, University of Calgary, Calgary, AB, Canada
Douglas Mahoney, Jennifer A. Chan & A. Sorana Morrissy

Authors

Varsha Thoppey Manoharan
View author publications
You can also search for this author in PubMed Google Scholar
Aly Abdelkareem
View author publications
You can also search for this author in PubMed Google Scholar
Gurveer Gill
View author publications
You can also search for this author in PubMed Google Scholar
Samuel Brown
View author publications
You can also search for this author in PubMed Google Scholar
Aaron Gillmor
View author publications
You can also search for this author in PubMed Google Scholar
Courtney Hall
View author publications
You can also search for this author in PubMed Google Scholar
Heewon Seo
View author publications
You can also search for this author in PubMed Google Scholar
Kiran Narta
View author publications
You can also search for this author in PubMed Google Scholar
Sean Grewal
View author publications
You can also search for this author in PubMed Google Scholar
Ngoc Ha Dang
View author publications
You can also search for this author in PubMed Google Scholar
Bo Young Ahn
View author publications
You can also search for this author in PubMed Google Scholar
Kata Osz
View author publications
You can also search for this author in PubMed Google Scholar
Xueqing Lun
View author publications
You can also search for this author in PubMed Google Scholar
Laura Mah
View author publications
You can also search for this author in PubMed Google Scholar
Franz Zemp
View author publications
You can also search for this author in PubMed Google Scholar
Douglas Mahoney
View author publications
You can also search for this author in PubMed Google Scholar
Donna L. Senger
View author publications
You can also search for this author in PubMed Google Scholar
Jennifer A. Chan
View author publications
You can also search for this author in PubMed Google Scholar
A. Sorana Morrissy
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

Varsha Thoppey Manoharan: conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing—original draft, and writing—review and editing. Aly Abdelkareem: data curation, formal analysis, and methodology. Gurveer Gill: data curation and formal analysis. Samuel Brown: data curation, methodology, investigation, and visualization. Aaron Gillmor: data curation and formal analysis. Courtney Hall: formal analysis. Heewon Seo: formal analysis. Kiran Narta: data curation and formal analysis. Sean Grewal: methodology. Ngoc Ha Dang: resources and methodology. Bo Young Ahn: resources and methodology. Kata Osz: resources and methodology. Xueching Lun: resources and methodology. Laura Mah: resources and methodology, investigation, and visualization. Franz Zemp: resources, supervision, and methodology. Douglas Mahoney: resources, supervision, and methodology. Donna Senger: conceptualization, resources, supervision, methodology, and writing—review and editing. Jennifer Chan: conceptualization, resources, supervision, methodology, and writing—review and editing. Sorana Morrissy: conceptualization, resources, supervision, funding acquisition, visualization, methodology, writing—original draft, and writing—review and editing.

Peer review information

Kevin Pang and Wenjing She were the primary editors of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Review history

The review history is available as Additional file 4.

Corresponding authors

Correspondence to Donna L. Senger, Jennifer A. Chan or A. Sorana Morrissy.

Ethics declarations

Ethics approval and consent to participate

Collection and use of human samples with participant consent was covered under research ethics board approval number HREBA.CC-16–0762 (The Clark Smith Tumour Biobank), approved by the Health Research Ethics Board of Alberta. The generation of patient-derived xenografts from the material was additionally approved by the Animal Care Committee of the University of Calgary under protocol AC21-0016. All experimental methods comply with the Helsinki Declaration.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Supplemental tables.

Additional file 2. Supplemental figures and materials.

Additional file 3. Descriptions of Additional file 1: Supplemental tables.

Additional file 4. Review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Manoharan, V.T., Abdelkareem, A., Gill, G. et al. Spatiotemporal modeling reveals high-resolution invasion states in glioblastoma. Genome Biol 25, 264 (2024). https://doi.org/10.1186/s13059-024-03407-3

Download citation

Received: 08 December 2023
Accepted: 29 September 2024
Published: 10 October 2024
DOI: https://doi.org/10.1186/s13059-024-03407-3

Spatiotemporal modeling reveals high-resolution invasion states in glioblastoma

Abstract

Background

Results

Conclusion

Background

Results

Spatial transcriptome profiling and unsupervised program discovery in GBM xenograft models

Gene expression program discovery in human and mouse cells

Tumor programs span progenitor states to invasion phenotypes

Tumor program association with survival and transcriptional subtypes

Microenvironment programs encompass cell types and states

Tumor microenvironment crosstalk

Molecular and structural contributors to invasion along distinct routes of tumor cell travel

Discussion

Conclusion

Methods

Intracranial BTIC models

Visium library preparation and sequencing

Spatial transcriptomics data preprocessing

Publicly available scRNA-seq, stRNA-seq, and STARmap datasets

Admixture calculation

Selection of over-dispersed (OD) genes

Matrix factorization and rank selection

Annotation of gene expression programs

Spatial concordance of programs

SCENIC

CNV analysis

GBM cell-state signatures

Cell–cell communication

STRINGdb networks and selection of hub genes

Selection of hub genes

Hub gene characterization

Survival analysis

Over-representation analysis

Differential expression analysis

Availability of data and materials

References

Funding

Author information

Authors and Affiliations

Contributions

Peer review information

Review history

Corresponding authors

Ethics declarations

Ethics approval and consent to participate

Competing interests

Additional information

Publisher’s Note

Supplementary Information

Additional file 1. Supplemental tables.

Additional file 2. Supplemental figures and materials.

Additional file 3. Descriptions of Additional file 1: Supplemental tables.

Additional file 4. Review history.

Rights and permissions

About this article

Cite this article

Share this article

Genome Biology

Contact us