Topographic mapping of the glioblastoma proteome reveals a triple-axis model of intra-tumoral heterogeneity

Abstract

Glioblastoma is an aggressive form of brain cancer with well-established patterns of intra-tumoral heterogeneity implicated in treatment resistance and progression. While regional and single cell transcriptomic variations of glioblastoma have been recently resolved, downstream phenotype-level proteomic programs have yet to be assigned across glioblastoma's hallmark histomorphologic niches. Here, we leverage mass spectrometry to spatially align abundance levels of 4,794 proteins to distinct histologic patterns across 20 patients and propose diverse molecular programs operational within these regional tumor compartments. Using machine learning, we overlay concordant transcriptional information, and define two distinct proteogenomic programs, MYC- and KRAS-axis hereon, that cooperate with hypoxia to produce a tri-dimensional model of intra-tumoral heterogeneity. Moreover, we highlight differential drug sensitivities and relative chemoresistance in glioblastoma cell lines with enhanced KRAS programs. Importantly, these pharmacological differences are less pronounced in transcriptional glioblastoma subgroups suggesting that this model may provide insights for targeting heterogeneity and overcoming therapy resistance.

Fig. 1. Schematic overview detailing data generation, analysis and presentation.

a Hematoxylin and Eosin (H&E) images detailing the anatomical niches within GBM: leading edge (LE), infiltrating tumor (IT), cellular tumor (CT), palisading cells around necrosis (PAN), and microvascular proliferations (MVP). b Collection of clinical data for a cohort of 20 GBMs. c FFPE tissue was sectioned and stained with H&E to identify histomorphological features. These were then excised using LCM and prepared for LC-MS/MS-based proteomics. d Bioinformatic pipeline including multi-dimensional scaling, differential expression matrix analysis, and gene set enrichment analysis (GSEA). e Selected region-specific biomarkers were validated within an external cohort of patients by immunohistochemistry (IHC). Clinical significance of regions specific markers was approximated using stratified RNA levels from TCGA. f Data was deposited in an accessible online data portal to allow for interactive and real-time exploration and comparison of proteomic profiles within GBM: (https://www.brainproteinatlas.org/dash/apps/GPA).

Fig. 2. Gene expression in anatomic features.

a Annotated H&E section of GBM highlighting areas of palisading cells around necrosis (PAN), microvascular proliferations (MVP), infiltrating tumor (IT), and cellular tumor (CT) regions. b, c Multidimensional scaling of all proteins using principal component analysis highlight (b) region- and (c) patient-specific protein spatial distributions. Note the reduction of inter-patient variability of regional signatures compared to whole (bulk) tissue profiling. Overall, samples derived from the same anatomic feature, regardless of patient origin, were more similar to each other. d Pearson correlation within and between tumors based on hierarchical clustering of all proteins (n = 78 samples). Anatomical features are a strong driver of proteomic intra-tumoral heterogeneity in GBM. e Gene Set Enrichment Analysis (GSEA) is based on all samples and their comparisons against other sample types. Normalized enrichment score (NES) is derived from the GSEA output and accounts for differences in gene set size and in correlations between gene sets in the expression dataset. Source data are provided as a Source Data file.

Fig. 3. Region-specific biomarker identification and validation.

a Spearman Rank correlation between histomorphological transcriptomic (Ivy GAP) and our own proteomic profiles, leading-edge (LE), infiltrating tumor (IT), cellular tumor (CT), palisading cells around necrosis (PAN), and microvascular proliferations (MVP). Only genes that were identified in both datasets were included in the analysis. b Differential expression matrix analysis based on genes identified based on all samples based on a two-sided t-test at p < 0.01 and 76 degrees of freedom. Values along the diagonal are the number of genes that are differentially abundant against all other regions. c Differential expression matrix comparison by Venn diagram across transcriptionally (Ivy GAP) and proteomically significant biomarkers. Only genes that were identified in both datasets were included. d Regional enrichment of AKAP12 by boxplot highlights increased expression within PAN (n = 154). Data are presented as median values ±  IQR and min/max values (whiskers). e H&E image of GBM highlighting necrosis, PAN, and CT. f Region-specific validation within an external cohort IHC demonstrates high staining of AKAP12 within PAN. g Regional enrichment of CD276 using boxplot highlights increased expression within MVP (n = 154). Data are presented as median values ±  IQR and min/max values (whiskers). Region-specific validation within an external cohort by IHC demonstrates (h) high staining of CD276 within MVP but (i) no staining of vasculature in LE. j Regional enrichment of PTPRZ1 by boxplot highlights increased expression within IT (n = 154). Data are presented as median values ±  IQR and min/max values (whiskers). k H&E image of GBM highlighting LE and IT. l Region-specific validation within an external cohort by IHC demonstrates high staining of PTPRZ1 within tumor cells of the IT. Source data are provided as a Source Data file.

Fig. 4. Single-sample gene signature analysis of highly pure GBM tumour fractions reveal two biologically distinct clusters.

a Two clusters of samples, independent of histomorphologic niches, are revealed by hierarchical clustering of samples using the ssGSEA scores from 64 selected gene signatures. Teal represents palisading cells around necrosis (PAN) samples while orange represents cellular tumor (CT) samples. b Functional drivers associated with each cluster using psSubPathways software package. c Inverse correlation between the KRAS targets and MYC targets signatures in the protein dataset generated in this study (n = 34 samples); additional null hypothesis testing was performed by permuted linear model fit analysis (lmPerm R package, p = 0.024). The 95% confidence interval is shown as grey areas. d 3D scatter plot of ssGSEA values of KRAS targets, MYC targets, and hypoxia gene signatures using proteomics (this study). e Transcriptomics (Ivy Gap) data from histomorphologically-defined samples. f Functional profiling by GSEA of Ivy Gap CT samples for the groups of samples labeled as MYC targets-high (n = 9) vs KRAS targets-high (n = 5) while presenting low activation of the hypoxia signature. g ROC analysis of clustering based on KRAS_targets, MYC_targets, or HYPOXIA gene signatures compared to that using Verhaak et al. mesenchymal, proneural and classical gene sets using cut-off optimization with the Youden method. h Distribution of KRAS_targets and MYC_targets signature enrichments in GBM at the single-cell level (Richards et al.) where (i) Mitotic spindle enrichment is quantified across the three detected populations. j Cell migration process enrichment across the three detected populations. Source data are provided as a Source Data file.

Fig. 5. High levels of the “KRAS targets” signature is associated with clinical aggressiveness in GBM patients.

a Schematic depicting non-tumour tissue content influencing overall molecular signatures and how computational methods that allow for niche-specific inference may lead to robust actionable insights. b A synthetic dataset was built from pathologically defined leading edge (LE), cellular tumor (CT), palisading cells around necrosis (PAN), and microvascular proliferations (MVP) tissue samples by generating computationally simulated mixtures; this dataset was used to train XGBoost models to infer the KRAS targets, MYC targets and hypoxia signatures in a CT-niche specific manner. c Association between KRAS and MYC target signatures in the RNA Ivy Gap dataset (n = 54). Null hypothesis testing was performed by permuted linear model fit analysis. The 95% confidence interval is shown as grey areas. d CT-niche-specific inference in the TCGA-GBM dataset (n = 171). Null hypothesis testing was performed by permuted linear model fit analysis. The 95% confidence interval is shown as grey areas. e CT-niche-specific inference in the CPTAC-GBM dataset (n = 100). Null hypothesis testing was performed by permuted linear model fit analysis. The 95% confidence interval is shown as grey areas. f Kaplan-Meier survival curves for samples split into “high” and “low” KRAS targets activity using ssGSEA of CT-niche specific inference in a proteomics dataset (CPTAC) and (g) at the RNA level (TCGA-GBM). Source data are provided as a Source Data file.

Fig. 6. Pharmacological profiling reveal axis-specific drug vulnerabilities and resistance.

a Volcano plots showing pharmacotranscriptomic comparisons of drug sensitivities across cells lines ranking high and low along with the MYC-, KRAS- and hypoxia axis. b Schematic overview of pharmacological profiling of hypoxic and normoxic cell populations to highlight potential differences in drug sensitivities along this axis (left). c Western blots analysis highlights increased CA9 expression with decreasing oxygen levels indicating a downstream response to hypoxia. Relative expression of palisading cells around necrosis (PAN) associated marker AKAP12 by western blot reveals increased expression in hypoxic conditions. Vinculin serves as a loading control. Densitometry analysis of AKAP12 western blots highlights an average 2-fold increase between 21% oxygen and 0.2% oxygen conditions across 4 different cell lines. d Differential viability of hypoxic versus normoxic cell populations (y-axis) upon treatment with 188 compounds and averaged over two replicates, overall cell viability is the average of the cell populations against the reference population (x-axis) and targets with the greatest differential viability are colored by the target as specified in the top legends. Blue dots represent compounds with the minimal differential response. d Relative viability effects on GSC proliferation in a kinome screen under differential oxygen concentrations. e Spheroid images upon treatment of GSK1398477, a PDGFRB inhibitor, in hypoxic and normoxic conditions. f Model of the heterogeneous co-existence of GBM populations driven by the KRAS or MYC target genes and signatures defined on the hypoxia pathway axis. Source data are provided as a Source Data file.