Integrated molecular and multiparametric MRI mapping of high-grade glioma identifies regional biologic signatures

Abstract

Sampling restrictions have hindered the comprehensive study of invasive non-enhancing (NE) high-grade glioma (HGG) cell populations driving tumor progression. Here, we present an integrated multi-omic analysis of spatially matched molecular and multi-parametric magnetic resonance imaging (MRI) profiling across 313 multi-regional tumor biopsies, including 111 from the NE, across 68 HGG patients. Whole exome and RNA sequencing uncover unique genomic alterations to unresectable invasive NE tumor, including subclonal events, which inform genomic models predictive of geographic evolution. Infiltrative NE tumor is alternatively enriched with tumor cells exhibiting neuronal or glycolytic/plurimetabolic cellular states, two principal transcriptomic pathway-based glioma subtypes, which respectively demonstrate abundant private mutations or enrichment in immune cell signatures. These NE phenotypes are non-invasively identified through normalized K2 imaging signatures, which discern cell size heterogeneity on dynamic susceptibility contrast (DSC)-MRI. NE tumor populations predicted to display increased cellular proliferation by mean diffusivity (MD) MRI metrics are uniquely associated with EGFR amplification and CDKN2A homozygous deletion. The biophysical mapping of infiltrative HGG potentially enables the clinical recognition of tumor subpopulations with aggressive molecular signatures driving tumor progression, thereby informing precision medicine targeting.

Fig. 1. Multiregional biopsy and MRI-based tumor sampling from a cohort of glioma patients.

a MRI contrast enhancement-based sampling of glioma specimens. b Circos plot indicating the molecular assay and MRI annotation of multiregional samples. c Tumor purity has been inferred from WES (available for 302 samples) and compared between contrast-enhancing (CE) and contrast non-enhancing (NE) samples within IDH wild-type group (left; MRI annotated IDH wild-type samples = 253), within IDH-mutant group (middle; MRI annotated IDH-mutant samples = 40), and between IDH wild-type and IDH-mutant samples (right). The middle line corresponds to the median; the lower and upper lines show the first and third quartiles. Difference of the purity mean among groups was assessed using the two-sided Mann–Whitney–Wilcox test. Source data are provided as a Source data file. d Schematic of imaging features extracted for this study and their phenotypic correlates. Features are separated into conventional and advanced MRI features.

Fig. 2. Somatic genetic mutations, copy number alterations, and correlates of imaging features in IDH-mutant glioma.

a For each IDH-mutant tumor (listed on horizontal axis), genetic variants have been annotated as private (exclusively occurring in one sample), shared (occurring in two or more samples, but not in all samples) and truncal (occurring in all samples) and reported as a percentage of the total number of somatic variations. The proportion of mutation types was significantly different between CE and NE (two-sided Fisher’s exact test p = 3.43e−67). A schematic example of multiregional tumor evolution is represented as a tree in which truncal, shared, and private branches are distinguished. b Overview of somatic alterations in IDH-mutant samples grouped by patient. Mutation load and tumor purity are reported in the top barplot and heatmap track, respectively. Clinical annotation and gene expression classification are indicated in the bottom tracks. Gene alteration frequency in the patient cohort is indicated as percentage on the left, known with driver mutations highlighted in red. c MEM model derived estimated marginal mean of T2W in IDH-mutant vs. IDH wild-type samples in the NE. Error bars show 95% confidence interval. Two-sided t-test with Tukey correction (n = 86). d MEM model derived estimated marginal mean of EPI + C in IDH-mutant vs. IDH wild-type samples in the NE. Error bars show 95% confidence interval. Two-sided t-test with Tukey correction (n = 89). c, d Data are presented as mean values +/−SD. a, c, d Source data are provided as a Source data file.

Fig. 3. Molecular alteration landscape of IDH wild-type glioma.

a The spectrum of somatic genetic alterations occurring in multiregional samples (n = 255) from IDH wild-type glioma patients (n = 48) indicated the frequency of truncal, shared, and private events within each single patient (left bars), within all patients (right bar, Total), and within contrast-enhancing and contrast non-enhancing samples (right bars, CE, and NE, respectively). The proportion of mutation types was significantly different between CE and NE (two-sided Fisher’s exact test p = 4.85e−43). Source data are provided as a Source data file. b Somatic genetic mutations and copy number alterations occurring in IDH wild-type glioma (260 samples from 53 patients). Samples are grouped by patients. Mutation load and tumor purity are reported in the top barplot and heatmap track, respectively. Clinical annotation and gene expression classification are indicated in the bottom tracks. Gene alteration frequency in patient cohort is indicated as percentage on the left.

Fig. 4. Mutual genetic alteration profiles of EGFR and NF1 in IDH wild-type glioma.

a EGFR and NF1 somatic genetic alterations occurring in multiregional samples from IDH wild-type glioma. Samples are grouped by patient and clinical annotations are indicated in the bottom tracks. EGFR (red box) and NF1 (green box) alterations are mutually exclusive in 98.7% samples (152 out of 154, two-sided Fisher’s Exact test p = 1.72e−07). The mutual exclusivity between EGFR and NF1 alterations was also observed within the same patient (orange box). Co-occurring genetic alterations have been identified only in 2 samples from 2 patients (blue box). b Evolutionary model of glioblastoma from patient P065 that included 5 contrast-enhancing samples harboring mutual exclusive EGFR and NF1 alterations. NF1 truncating mutation specifically occurred in two samples with wild-type EGFR locus; conversely, EGFR was amplified in the other two NF1 wild-type samples. Tumor purity is indicated for each sample (percentages displayed). VAF, variant allele frequency. Number of truncal, shared, and private alterations are indicated.

Fig. 5. Spatial heterogeneity of EGFR alteration and EGFR-associated imaging phenotypes.

a Three-dimensional visualization of EGFR alterations and their associated MRI features. b The molecular evolution of glioblastoma from patient P129 inferred from the occurrence of genetic alterations as truncal, shared, and private events across the multiregional specimens (n = 11, 2 contrast-enhancing and 9 contrast non-enhancing samples). The length of branches in the evolutionary tree (left panel) is proportional to the number of occurring alterations. Truncal driver alterations (CDKN2A deletion and PTEN frameshift mutation), and non-truncal multiple EGFR alterations have been reported along the evolutionary tree and in the oncoprint (right panel). Mutation load, tumor purity, MRI contrast enhancement annotation, and gene expression classification have been reported. c The percent variance attributed to each fixed term (y-axis) in MEMs for each imaging variable (x-axis) separated by region. *p-value < 0.05. EGFR*CDKN2A indicates the interaction of EGFR and CDKN2A (n = 221 biopsy samples). Statistical test: ANOVA. d MEM model derived estimated marginal mean of MD for EGFR and CDKN2A genotypes in the NE. Intermediate genotypes (Int. Gen.) denote genotypes that are not double wild type or EGFR amplified and CDKN2A homozygous loss. Error bars show 95% confidence interval (n = 74 biopsy samples). Statistical test: two-sided t-test with Tukey correction. Data are presented as mean values +/−SD. c, d Source data are provided as a Source data file.

Fig. 6. Molecular tumor evolution in IDH wild-type glioma.

a Alterations exclusive to the NE region are displayed from a subset of IDH wild-type patients. Alterations specific to the NE region were annotated with gene ontology terms and grouped based on ontology similarity. b Evolutionary trajectories of genetic alterations have been predicted by comparing the multiregional molecular profiles in a subset of IDH wild-type glioma patients (n = 34) with more than one truncal driver event identified. Four evolutionary models of IDH wild-type glioma have been proposed from the supervised clustering of repeated initiating trajectories. The number of trajectories observed within each cluster is reported (light green boxes); the number of clonal (red) and subclonal (blue) events is indicated for each alteration.

Fig. 7. Conventional MRI, transcriptomic, and genotypic characterization of NE region phenotypes.

a Unsupervised hierarchical clustering of the 158 multiregional glioma samples; rows are the 2826 most variable genes. b Pie charts show the frequencies of pathway-based classifications of MRI CE samples (top) and MRI NE samples of unsupervised cluster 1 and 2, respectively (bottom). c Correlations between the genotypic and euclidean distance of paired samples with the same pathway-based classification. d Genetic distance of CE samples to NE samples of each pathway-based classification (n = 94 biopsy samples). Boxplots represent data minimum, 25th percentile, 50th percentile, 75th percentile, and maximum. The p-values are indicated above each comparison in the figure. e Private, shared, and truncal alterations in individual samples in the NE region classified as each pathway-based subtype (from left to right: glycolytic/plurimetabolic, mitochondrial, neuronal, and proliferative/progenitor), with the average of private, shared, and truncal mutations for each pathway-based classification displayed to the right. f The proportion of truncal mutations vs. non-truncal (private and shared) mutations in samples of NEU subtype was significantly different than the proportion of truncal mutations vs. non-truncal mutations in the other subtypes (one-tailed Fisher’s exact test p = 6.13e−27). g Box and whisker plots show the absolute number of total (left) and private (right) mutations in each pathway-based classification and the distribution of mutational burden across samples (n = 51 biopsy samples). Boxplots represent data 25th percentile, 50th percentile, and 75th percentile. The upper whisker extends from the upper hinge to the largest value no further than 1.5X IQR (inter-quartile range, or distance between the first and third quartiles) and the lower whisker extends from the lower hinge to the smallest value at most 1.5X IQR. c, d, e, g Source data are provided as a Source data file.

Fig. 8. Associations between imaging variables and pathway-based signatures with subsequent phenotypic modeling.

a Transcriptomic pathway enrichment analysis for samples binned as high vs low for each image feature. Gene Set Enrichment Analysis (GSEA, Kolmogorov–Smirnov-like test as implemented in clusterProfiler1) of supervised differential analysis (Mann–Whitney–Wilcox test) among samples labeled as high vs low for each image feature. Size of the dot represents the adjusted p-value of significant enriched GO BP terms (Benjamini-Hochberg adjustment); color of the dots represents the Normalized Enrichment Score (NES) of the terms. b, c Scatter plots showing significant correlations between T1 + C and NEU (b) or PPR (c) signatures across all samples. d, e Scatter plots showing significant correlations between nK2 and NEU (d) or GPM (e) in NE samples only (indicated by blue center on data points). b–e Statistical test: two sample correlation t-test. f, g Average R2*(t) curves for regions corresponding to a NEU (f) and GPM (g) sample. R2*(t) is used to derive nK2. h, i 3D renderings demonstrating conditions illustrative of homogenous uniform cell sizes (h) with 0% coefficient of variation and heterogenous mixed cell sizes (i) with 6.5% coefficient of variation. b–i Source data are provided as a Source data file.