Integrative spatial analysis reveals a multi-layered organization of glioblastoma

Abstract

Glioma contains malignant cells in diverse states. Here, we combine spatial transcriptomics, spatial proteomics, and computational approaches to define glioma cellular states and uncover their organization. We find three prominent modes of organization. First, gliomas are composed of small local environments, each typically enriched with one major cellular state. Second, specific pairs of states preferentially reside in proximity across multiple scales. This pairing of states is consistent across tumors. Third, these pairwise interactions collectively define a global architecture composed of five layers. Hypoxia appears to drive the layers, as it is associated with a long-range organization that includes all cancer cell states. Accordingly, tumor regions distant from any hypoxic/necrotic foci and tumors that lack hypoxia such as low-grade IDH-mutant glioma are less organized. In summary, we provide a conceptual framework for the organization of cellular states in glioma, highlighting hypoxia as a long-range tissue organizer.

Keywords: glioblastoma; glioma; hypoxia; intratumor heterogeneity; spatial proteomics; spatial transcriptomics.

Figure 1.. Experimental design and spot classification per sample.

(A) Experimental design and patient cohort. Fresh frozen tissue sections from GBM (n=13) and IDH-mutant gliomas (n=6) were profiled by 10X Visium. CODEX was performed on 12 near-adjacent GBM tissue sections with a panel of 40 antibodies. Four tumors were spatially annotated by the neurosurgeon during navigation-guided surgery. This GBM Visium cohort was combined with an external GBM Visium cohort for joint analysis. Scheme created with BioRender.com. (B) Copy number aberrations (CNAs) were inferred by average relative expression in sliding windows of 150 analyzed genes after sorting genes by their chromosomal location. ZH1019, from which we profiled both an infiltrating sample and a T1-constrast enhancing sample, is shown here as an example. Rows correspond to spots arranged by malignancy level as inferred from CNA; columns correspond to genes arranged by chromosomal position. Annotation bars correspond to region from which the spot was derived (T1 or infiltrating) and malignancy level. (C) Spatial maps of ZH1019 T1 contrast-enhancing and ZH1019 infiltrating samples with spots annotated by malignancy level as described in (B). (D) Per sample Leiden (left) and NMF clustering (right) of ZH916bulk for vascular and hypoxia clusters (all other clusters are shown in gray).

Figure 2.. Deriving and annotating spatial metaprograms (MPs).

(A) Scheme of metaprogram (MP) generation approach created with BioRender.com. Each sample is clustered individually by Leiden and NMF. All Leiden and all NMF programs across all samples are clustered by their gene overlap. Each cluster is collapsed to a consensus MP by selecting for the most recurrent genes across programs within the cluster. (B) Similarity matrix based on gene identity overlap (quantified by Jaccard index), for all programs derived from NMF and Leiden clusters. Programs are annotated by sample and region identity. Cluster numbers correspond to the table in (C). (C) Table of metaprogram names and selected genes corresponding to clusters numbered in (B). Malignant metaprograms are depicted in bold. (D) Enrichments of spatial metaprograms (rows) with gene-sets (columns) previously defined from studies indicated at the top.^,,,– Enrichment calculated by hypergeometric test (-log10 FDR-adjusted P-values). (E) For each sample (represented by one bar), proportions of spots assigned to each MP are shown. Four tumors with multiple sections are at the right and are grouped by tumor. Below are spatial MP maps of several ZH881 tissue sections, demonstrating differences in composition and spatial organization between sections isolated from the same tumor.

Figure 3.. Spatial profiling of gliomas by CODEX.

(A) Schematic workflow of CODEX experiment, image processing and computational analysis. Created with BioRender.com. (B) Protein markers profiled by CODEX. (C) Relative protein expression per cell type (or state) by z-score. Columns represent proteins. Each row represents the indicated cell type/state in one sample. Number of rows differ between cell types/states because some samples lack specific cell types/states. (D) CODEX staining of representative examples: Top row shows images with the indicated markers. Bottom row shows the corresponding nuclear segmentation masks with cytoplasmic expansion by 3 μm. Masks are colored by cell state. Scale bar = 20μM. (E) Voronoi diagram showing cell type/state abundances across all samples. Dashed lines separate malignant from non-malignant and immune cells. (F Spatial maps of cell types/states annotations for (i) CODEX single cell data, (ii) CODEX pseudospots and (iii) Near-adjacent section of ZH1019 infiltrating following alignment between Visium and CODEX with STalign (Methods), median Pearson correlation between CODEX and Visium =0.664, p-value = 0.0105 by t-test). Visium spots annotated by their highest scoring MP and CODEX pseudospots are annotated by the most abundant state within each spot. (G) Cumulative distribution plot showing the frequency of pseudospots in which the highest cell type/state, the second highest, or both cover a higher proportion of the cells than the value specified on the x-axis. (H) Exemplary CODEX image overlaid with pseudospot grid (white) and nuclear segmentation masks colored by cell type/state.

Figure 4.. Spatial distribution of cell states and organization zones.

(A) Scheme of ZH881_T1 illustrating spatial coherence score calculation depicting high spatial coherence score for a grouped cell state and a low spatial coherence score for a scattered cell state. (B) MP spatial coherence by sample. Standard deviation is shown in error bars. (C) MP mean spatial coherence across all samples; standard deviation is shown in error bars. Two-way ANOVA to compare variance across samples (p=0.000003) vs. across states (p=0.23). (D) Left: Spatial maps of two samples from the same tumor profiled by Visium showing multiple organization zones. Right: Heatmap of the relative abundance of the different organizational patterns (rows) in each sample (columns). Annotation bars shows per sample MES-Hyp abundance (yellow) and assignment to glioma type (green). (E) Left: Scheme describing the assignment to organizational zones. Spatial coherence and spot malignancy level (inferred from CNA) were used to assign spots to Structured-Malignant (Struct-Malig), Structured-Normal (Struct-Norm), Disorganized-Malignant (Disorg-Malig) and Disorganized-Normal (Disorg-Norm) in Visium data. Right: Structured and disorganized regions are also found in CODEX data. H&E of Struct-Malig zone shows microvascular proliferation (MVP) but no pseudopalisades or necrosis in a region that is structured at the CODEX cell state level. Upper scale bar=100μM; lower scale bar=50μM. (F) Boxplot of spatial coherence score of malignant MPs across regions (windows) with high abundance (>10%) of MES-Hyp vs. low abundance of MES-Hyp. P-values are shown on the plot (by Bonferroni adjusted t-test). Boxes indicate the median and the 1^st and 3^rd quartiles. The upper and lower whiskers extend to the maximal and minimal values. (G) Fraction of structured regions based on H&E staining (MVP and PAN) vs. fraction of structured regions based on CODEX spatial coherence (p=0.0019, by Wilcoxon rank sum test), boxes as defined in (F).

Figure 5.. Spatial associations between states across scale.

(A) Scheme depicting the three measures of spatial relationships between MPs across scales (from low to high resolution): (i) Regional composition (across sliding windows of r=2–15), (ii) Adjacency, and (iii) Colocalization (Methods). The example shown highlights coupling between MES-Ast and MES-Hyp across scales of resolution. (B) Nested circle plots depicting spatial relationship strength across scales of resolution. Top: Schematic of nested circle plot. Bottom: Example of associations with nested circle plots (Visium) adjacent to CODEX images, scale bar indicated in the images. Yellow line shows the distance from hypoxia (NDRG1) to blood vessels (CD31) ~160μm. Note: MES-Hyp/Vasc association is NA (white) at the colocalization and adjacency levels due to too low abundance of vascular spots within hypoxic regions to perform the calculation. (C) Left: Line plot showing the number of significant interactions across different analyses in structured vs. disorganized regions. Right: Venn diagram of the number of significant interactions across the regional composition analysis (see criteria for ‘consensus interactions’ in Methods). (D) CODEX images highlighting selected pairs of coupled cell types from (C). Left: AC/OPC coupling (scale bar=20μM). Right: Spatial relationships in the vascular niche (scale bar=40μM). While Vasc and immune cells are immediately adjacent, Vasc/MES coupling occurs over larger areas. (E) Summary heatmap of state-state associations across scales in structured samples covering most consensus interactions and other pairs of interest (see also Data S1). Each column represents the summary of a different spatial relationship measure (Visium data) across all GBM samples, with the final two columns corresponding to CODEX colocalization and regional composition, respectively. Dots are colored by mean scaled relationship strength (coupling score) (Methods) and dot size corresponds to the fraction of samples in which the relationship is significant by Fisher’s exact test. (F) Top: Network graphs representing couplings between normal brain cell types or Vasc with the malignant states, Bottom: Couplings of immune cells with malignant states. In each graph, the central node is a non-malignant cell type, and the surrounding nodes represent cancer states. Edge strength represents the mean scaled relationship strength across all measures (Visium). Only consensus interactions are shown.

Figure 6.. A layered model of GBM spatial organization

(A) Top: Network graph with nodes representing cell types/states and edges representing recurrent interactions (mean scaled coupling score >0.35 across samples and significantly coupled in > 20% of samples for at least two measures of spatial relationships). Edge color represents mean coupling score across levels of resolution, for structured regions of Visium data. Edges with dashed lines represent connections with cell types coming from CODEX (T-cells, B-cells). Bottom: Scheme showing gradients of hypoxia, hypoxia response and infiltration in alignment with the layers. (B) CODEX image showing the indicated cell types and markers, representing layers 1 – 4. Layer 5 was added as scheme depicting normal brain. Created with BioRender.com. (C) Spatial maps for sample ZH881_1A: (i) Visium sample annotated by IVY Gap histological feature-associated transcriptional programs and state layers (ii) CODEX sample annotated by state layers and neuropathologist annotations of H&E staining from the CODEX sample. (D) Spatial maps comparison for Visium samples annotated by IVY Gap histological feature-associated transcriptional programs (left) vs. annotated by our layers (right) (E) Stacked bar plot showing MP composition per Ivy GAP histological feature transcriptional program annotation (left) and CODEX MP composition per histological annotations (right). Bars below 0 represent non-malignant MPs.