Pharmacogenomic profiling reveals molecular features of chemotherapy resistance in IDH wild-type primary glioblastoma

Abstract

Background: Although temozolomide (TMZ) has been used as a standard adjuvant chemotherapeutic agent for primary glioblastoma (GBM), treating isocitrate dehydrogenase wild-type (IDH-wt) cases remains challenging due to intrinsic and acquired drug resistance. Therefore, elucidation of the molecular mechanisms of TMZ resistance is critical for its precision application.

Methods: We stratified 69 primary IDH-wt GBM patients into TMZ-resistant (n = 29) and sensitive (n = 40) groups, using TMZ screening of the corresponding patient-derived glioma stem-like cells (GSCs). Genomic and transcriptomic features were then examined to identify TMZ-associated molecular alterations. Subsequently, we developed a machine learning (ML) model to predict TMZ response from combined signatures. Moreover, TMZ response in multisector samples (52 tumor sectors from 18 cases) was evaluated to validate findings and investigate the impact of intra-tumoral heterogeneity on TMZ efficacy.

Results: In vitro TMZ sensitivity of patient-derived GSCs classified patients into groups with different survival outcomes (P = 1.12e-4 for progression-free survival (PFS) and 3.63e-4 for overall survival (OS)). Moreover, we found that elevated gene expression of EGR4, PAPPA, LRRC3, and ANXA3 was associated to intrinsic TMZ resistance. In addition, other features such as 5-aminolevulinic acid negative, mesenchymal/proneural expression subtypes, and hypermutation phenomena were prone to promote TMZ resistance. In contrast, concurrent copy-number-alteration in PTEN, EGFR, and CDKN2A/B was more frequent in TMZ-sensitive samples (Fisher's exact P = 0.0102), subsequently consolidated by multi-sector sequencing analyses. Integrating all features, we trained a ML tool to segregate TMZ-resistant and sensitive groups. Notably, our method segregated IDH-wt GBM patients from The Cancer Genome Atlas (TCGA) into two groups with divergent survival outcomes (P = 4.58e-4 for PFS and 3.66e-4 for OS). Furthermore, we showed a highly heterogeneous TMZ-response pattern within each GBM patient using in vitro TMZ screening and genomic characterization of multisector GSCs. Lastly, the prediction model that evaluates the TMZ efficacy for primary IDH-wt GBMs was developed into a webserver for public usage ( http://www.wang-lab-hkust.com:3838/TMZEP ).

Conclusions: We identified molecular characteristics associated to TMZ sensitivity, and illustrate the potential clinical value of a ML model trained from pharmacogenomic profiling of patient-derived GSC against IDH-wt GBMs.

Keywords: Cancer genomics; Glioblastoma; Intra-tumoral heterogeneity; Machine learning; Pharmacogenomics; Temozolomide.

Fig. 1

In vitro TMZ screening of patient-derived glioblastoma (GBM) stem cell (GSC) predicts GBM prognosis. a TMZ sensitivity determination pipeline for patient-derived GSC. GR: growth rate; AUC: area under the conventional drug response curve; GSC: Glioblastoma Stem Cell; TMZ: temozolomide. b GR50 and AUC z-score comparison between MGMT methylated and un-methylated samples. Wilcoxon rank sum tests were performed for p-values (* P < 0.05). c, d Comparison of progression-free survival (c) and overall survival (d) between the TMZ-resistant (red) and TMZ-sensitive (blue) patients. The cohort was determined by TMZ screening in panel a. Sen, TMZ-sensitive; Res, TMZ-resistant. P-values were calculated by logrank test

Fig. 2

Somatic mutational landscape of the main cohort. TMZ sensitivity was determined based on in vitro TMZ screening. The single nucleic variants (SNV) including point mutations, short insertion/deletions, and copy number amplifications and deletions of selected GBM driver genes were included. Alteration frequencies are shown on the right side. When available, copy number alteration results were inferred from the methods in the following order; WES (highest priority), Gliomascan, and RNAseq. EGFRvIII was identified from RNA-seq data. Moderate: missense or inframe deletion; high: frameshift, stop gained, splice donor or splice acceptor; M, methylated; UM, unmethylated; N/A, not available; C, classical; P, proneural; M, mesenchymal; GS, Gliomascan; WES, whole exome sequencing

Fig. 3

Elucidation of expression markers associated with the resistance to temozolomide in patients with IDH-wt primary glioblastoma (GBM). a TMZ-resistant expression marker identification using RNA sequencing (n = 34). b Progression-free survival (upper panel) and overall survival (bottom panel) of IDH-wt, TMZ-treated primary GBM from the TCGA. High risk, z-score of gene expression > 2 in at least one of the TMZ-resistant marker genes (n = 13); others, the rest (n = 83). c Comparison of the expression level of the TMZ-Resistant expression markers in the initial and recurrent paired samples from the longitudinal sequencing cohort with 40 IDH-wt, TMZ-treated primary GBMs. Each gray line connects the gene expression level in one initial and recurrent pair. Wilcoxon rank sum tests were performed for p-values (* P < 0.05, ** P < 0.01). qnorm, quantile normalized

Fig. 4

Machine learning from the combined genomic and expression features predicts patient prognosis. a Bubble plot showing the trends of features in terms of TMZ-resistant and TMZ-sensitive. The bubble size indicates P-value, the color and location of the bubble indicate the log2 of TMZ-resistant ratio/TMZ-Sensitive ratio value. If the log2 TMZ-Resistant ratio/TMZ-Sensitive ratio value is positive, the bubble is colored in red, and if negative, it is colored in blue. Copy number gain and loss were not counted in this plot. del: deletion, amp: amplification, exp: expression, subtype: GBM subtype; CL, classical; PN, proneural; MES, mesenchymal; M, methylated; UM, unmethylated. P values on gene expression and MGMT fusion bubbles are by t-test, the rest are by Fisher’s exact test. b ROC curve in the training set (n = 69). All features include 25 features shown in a. P < 0.01, using features that are P < 0.01 in a (gene expression of ANXA3, PAPPA, EGR4; AUC: 0,81); P < 0.001, using features that are P < 0.001 in a (ANXA3 expression; AUC: 0.77); MGMT, only using MGMT promoter status as prediction feature. c Sankey diagram showing confusion matrix of resistant and sensitive samples in the training dataset. Sen, TMZ-sensitive; Res, TMZ-resistant. d, e Survival curves of TCGA IDH-wt, TMZ-treated primary GBM samples which the TMZ response has been predicted by machine learning. P-values were calculated by logrank test. f, g Survival curves of mesenchymal TCGA samples separated by predicted TMZ response. P-values were computed by log-rank test

Fig. 5

TMZ screening in multi-sector samples underscores intra-tumor heterogeneity of drug response. a Experimental design for screening multi-sector samples. Sen, TMZ-sensitive; Res, TMZ-resistant. b Molecular features of multi-sector samples. I, initial tumor; R, recurrent tumor; WT, wild-type; mut, mutant; M, methylated; UM, unmethylated; N/A, not available; C, classical; P, proneural; M, mesenchymal. c, d Phylogenetic trees of somatic mutation evolution in multi-sector samples from c patient M13 and d patient M14. The length of the branch is relative to the number of mutations. Dashed lines indicate a relatively larger number of mutations that cannot be scaled for visualization. Indicated alterations are GBM driver alterations and RNA expression of TMZ-resistant markers. amp, amplification; del, deletion; higher_exp, higher transcriptomic expression compared to other sample/samples. Blue, TMZ-sensitive; red, TMZ-resistant. e Comparison of Concurrent CNAs in PTEN, EGFR, and CDKN2A/B in the main cohort. P value by Fisher’s exact test. f Overall survival difference in patients with multi-sector samples identified as S, all sensitive (M1~M4); H, heterogeneous (M11~M18); R, all resistant (M5~M10). P-values calculated by logrank test. g Detection rate of TMZ heterogeneity by the number of multi-sector samples. h Relative distribution of TMZ response by the number of multi-sector samples

Fig. 6

Summarizing scheme for the known and newly identified molecular features associated to TMZ response in GBM. References supporting the associations are shown next to the arrows [, , –45]. Features with gray dotted lines are the proposed association from this study. Small arrows pointing upward inside the bubbles indicate activated signaling or highly expressed