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

In vitro screening using patient-derived GSCs reflects personalized TMZ efficacy

To evaluate GBM’s response to TMZ, we performed in vitro TMZ cytotoxicity assays in short-term (6 days) cultured patient-derived GSCs (n = 69, main cohort) obtained from surgically resected IDH-wt primary GBM specimens. Since conventional metrics such as the effective concentration at 50% (IC50) or maximum inhibition % (Emax) highly depends on cell division rate obscuring accurate sensitivity prediction, we adopted GR inhibition metrics, which are independent of division number and therefore superior to conventional metrics for assessing the effects of drugs in fast dividing cells [34]. We calculated GR50 values for each sample, and for those with infinite GR50 values, we measured conventional AUC values. By calculating Z-scores for GR50 and AUC, we divided our samples into TMZ-Sensitive and TMZ-Resistant groups (Fig. 1a). As expected, MGMT promoter methylation was observed to be related to Z-scores of GR50 and AUC values (Fig. 1b, Wilcoxon rank sum test P = 0.018) [35].

Fig. 1figure 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

Strikingly, TMZ-Resistant (n = 29) and TMZ-Sensitive (n = 40) groups defined from in vitro sensitivity were highly predictive of survival outcomes for patients who were under a TMZ-based treatment regimen (Fig. 1c and d; PFS, P = 1.12e−4; OS, P = 3.63e−4; by log-rank test). Notably, the above-defined in vitro sensitivity surpasses the well-known MGMT promoter methylation status in predicting the patient prognosis (Additional file 2: Fig. S4a). Additionally, a Cox-regression multivariate survival analysis considering age, gender, extent of resection, and MGMT promoter methylation revealed that in vitro TMZ sensitivity and the extent of resection were independent factors associated with PFS and OS, while MGMT promoter methylation being related to in vitro sensitivity (Additional file 2: Fig. S3b, Fisher’s exact test P = 0.0156) was marginally significant (Table 1). Collectively, these data reflect the reliability of our preclinical TMZ testing system for assessing clinical response to TMZ in patients newly diagnosed with IDH1-wt GBM.

Table 1 Cox-regression analysis of survival outcomes in the 69 cases screened with TMZ in vitroGenomic analysis reveals somatic mutational landscape of TMZ-resistant and sensitive groups

To identify genetic factors contributing to TMZ response, we explored somatic genomic alterations in the TMZ-resistant and sensitive groups in our main cohort. WES and/or GliomaSCAN on 57 tissue specimens (with matched blood controls) and RNA-seq on 34 tissue specimens were either newly performed or downloaded from previous publications [11, 20] (Additional file 1: Table S1). Somatic SNVs and short insertions/deletions were identified by SAVI2 [18] (Additional file 1: Table S2). A sample was labeled as hypermutated if the total number of somatic mutations was over 350 by WES. CNAs were calculated from WES, GliomaSCAN, or were predicted from RNA-seq by CNAPE [26] (Additional file 1: Table S2-3, Additional file 2: Fig. S5-S6, and Additional file 3: Supplementary Methods). Variants with VAF over 5% and CNAs in previously reported GBM driver genes, together with EGFRvIII (Additional file 1: Table S4) and the expression-based GBM subtyping were shown in Fig. 2. Overall, no significant genomic difference was observed between the responder and non-responder groups. Yet, mesenchymal/proneural subtype and somatic mutations in genes including NF1, NF2, and PTEN were more often observed in TMZ-resistant samples, while PIK3R1 somatic mutations were slightly more frequent in TMZ-sensitive samples (Fig. 2). These observations indicate that GBM’s response to TMZ might be determined by the combination of multiple factors but not by single ones.

Fig. 2figure 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

Transcriptomic sequencing reveals marker genes of TMZ resistance

To identify marker genes of TMZ response, we first separated conserved and non-conserved gene expression between GSCs and the initial tumor tissue through a Gaussian mixture model (Additional file 2: Fig. S2). The conserved genes were then used to perform differential gene expression analysis between tissue RNA-seq data of TMZ-resistant (n = 12) and sensitive (n = 22) samples. Principal component analysis on these 34 samples showed slight clustering by GBM subtype (proneural vs. mesenchymal/classical) but by no other factors including age, gender, and MGMT promoter methylation status (Additional file 2: Fig. S3a, b). We identified four genes (EGR4, PAPPA, LRRC3, and ANXA3) significantly up-regulated in the TMZ-resistant group (Fig. 3a, Additional file 2: Fig. S3c, log2 fold change > 2.5, adjusted P < 0.01). To explore the prognostic value of the TMZ-resistant marker genes, we extracted 96 RNA-seq available TMZ-treated IDH-wt primary GBM patients from the TCGA dataset and classified them into high-risk and low-risk groups based on the expression of these four genes. Notably, the high-risk group had significantly worse PFS (P = 1.59e−03 by log-rank test) and OS (P = 3.46e−03 by log-rank test, Fig. 3b), compared to that of the low-risk group.

Fig. 3figure 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

To further investigate the expression change of these genes before and after TMZ treatment, we integrated a total of 40 paired RNA-seq data of initial and matched TMZ-treated recurrent IDH-wt GBM samples [18, 19]. As shown in Fig. 3c, the expression level of TMZ-resistant markers increased in the recurrent samples compared to the initial, suggesting that the TMZ-resistant marker-expressed cell population survived TMZ treatment and expanded in the recurrent GBM.

A machine learning (ML) approach for integrating key features to predict TMZ response of IDH1-wt GBM

Figure 4a presents the overall relevance of the genomic, transcriptomic, and other features on TMZ response. Along with the expression of four TMZ-resistant markers, MGMT expression, MGMT promoter methylation status, hypermutation status, GBM subtype, somatic mutations, and CNAs identified from the main cohort, we added 5-aminolevulinic acid (5-ALA) tendency [36] as another feature (Additional file 1: Table S1). In order to integrate these features for patient evaluation, we constructed an XGBoost classifier to identify the TMZ response of a patient as TMZ-resistant or TMZ-sensitive. Among the 30 features shown in Fig. 4a, 25 features were used to train the machine learning model in the main cohort, excluding NF2 mutation, hypermutation, 5-ALA positive and 5-ALA negative which were not available in the TCGA testing cohort (Additional file 1: Table S5). Compared with MGMT promoter status as the only feature, adding other features provided more information for recognition of TMZ non-responders (Fig. 4b). Notably, the top five informative features from the model were the expression level of ANXA3 and LRRC3, proneural subtype, ERG4 and MGMT expression (Additional file 2: Fig. S7a). In addition, incorporating the four expression markers together with other features achieved a stronger discrimination power compared to the presence of just an individual marker (Additional file 2: Fig. S7b). Within the training cohort (main cohort), a prediction of 88.4% (61 out of 69) of the samples matched the in vitro TMZ-response (Fig. 4c). We then tested our model in an independent cohort with 262 IDH-wt, TMZ-treated primary GBM patients from TCGA (inclusive of the 96 RNA-seq available patients from Fig. 3b). Importantly, patients predicted to be TMZ-resistant by the classifier had significantly worse PFS (Fig. 4d, P = 4.58e−04 by log-rank test) and OS (Fig. 4e, P = 3.66e−04 by log-rank test) validating the power of our model to predict prognostic outcome in patients treated by TMZ. Moreover, we investigated the survival difference across four subtypes (classical, proneural, neural, and mesenchymal) in the TCGA cohort and no significant segmentation was observed for PFS (P = 0.531 by multivariate log-rank test) and OS (P = 0.412 by multivariate log-rank test) in Additional file 2: Fig. S8a, which is expected and compatible with the observations previously reported by the TCGA group. We further correlated the four GBM subtypes with the TMZ response predicted from our machine learning model. Notably, the mesenchymal subtype is associated with TMZ resistance (P = 0.039 by Fisher exact test, Additional file 2: Fig. S8b). Among the mesenchymal cases, the resistant group demonstrated worse PFS (P = 1.98e−02 by log-rank test, Fig. 4f) and OS (P = 1.26e−04 by log-rank test, Fig. 4g), compared to that of the sensitive group, highlighting the value of our model to unveil new responders/non-responders within the subtypes. In addition, when compared to using only MGMT promoter methylation status in the TCGA dataset (n = 203), our model integrating multiple features provided a better way to segregate patients with different outcomes of both PFS and OS (Additional file 2: Fig. S9). Within MGMT methylated group, our model identified a limited number of high-risk resistant cases with worse PFS and OS (Additional file 2: Fig. S9c). Furthermore, to facilitate the use of our model, we designed a freely accessible website named TMZep that provides the function for evaluating potential TMZ response for GBM patients (http://www.wang-lab-hkust.com:3838/TMZEP) [37]. Users can input patient’s data on part or all of the 25 features to the website, which will evaluate the potential TMZ treatment response of the corresponding GBM patient.

Fig. 4figure 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

Multi-sector TMZ screening underlines intratumoral heterogeneity in drug responsiveness

Intratumoral heterogeneity (ITH) is a key factor that causes therapeutic resistance and recurrence in GBM [38, 39]. To reveal the impact of ITH on TMZ treatment, we compiled 52 GBM tumor tissue specimens from 18 patients, where each patient had 2 to 4 multi-sector samples taken from the same patient (Multi-sector cohort, Additional file 1: Table S1). Of the Multi-sector cohort, 15 tissues and patients overlapped with the main cohort. The three additional patients were patients with recurrent GBM, IDH mutant GBM, or IDH-wt primary GBM without TMZ treatment. We performed in vitro TMZ screening of the multi-sector GSCs (Fig. 5a) followed by WES in 19 samples and RNA-seq in 26 samples. Interestingly, almost half of the patients (8/18) carried both TMZ-resistant and TMZ-sensitive tumor samples (Fig. 5b and Additional file 2: Fig. S10). We termed these patients as TMZ-ITH, which harbored heterogeneous GSCs within one tumor in terms of in vitro TMZ treatment response. We confirmed several TMZ-associated factors identified earlier in this study by comparing the molecular signatures of multi-sectors. In particular, the TMZ-resistance markers were upregulated in the resistant sectors of M13 and M14 (Fig. 5c, d, Additional file 2: Fig. S11a-b). Meanwhile, a combination of PTEN loss, EGFR gain, and deeper deletion of CDKN2A/B was observed specifically in the sensitive sectors of these two patients (Fig. 5c, d, and Additional file 2: Fig. S11c). Motivated by this observation, we checked the concurrent CNAs in PTEN, EGFR, and CDKN2A/B back in our main cohort and found that it was significantly more frequent in TMZ-sensitive samples (Fisher’s exact P = 0.0102, Fig. 5e), while each individual factor did not have statistical significance.

Fig. 5figure 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

Back to the ITH analysis, eight TMZ-ITH patients had comparable survival time with patients harboring only TMZ-resistant sectors and significantly worse survival time than patients with only TMZ-sensitive sectors (OS, P = 0.027; PFS, P = 0.015; by log-rank test), indicating that although the TMZ treatment achieved the particular effect by eliminating a sensitive group of tumor cells, the resistant GSCs might quickly lead to tumor relapse (Fig. 5f and Additional file 2: Fig. S12). This observation underscores the importance of careful consideration of ITH via multi-sector evaluation before treatment delivery.

Since the number of sectors analyzed from a tumor may influence the possibility of TMZ-ITH detection, we evaluated the optimal number of sectors to observe TMZ-ITH in a patient. We demonstrated that when two sectors were taken from a tumor, the TMZ-ITH detection rate was around 30% (17/56), followed by 43% (12/28) with three sectors, and 50% (3/6) with four sectors (Fig. 5g). Interestingly, while the TMZ-ITH detection increases with multi-sector number, purely resistant groups decreased but not the purely sensitive patients, underscoring the existence of good responders of TMZ treatment (Fig. 5h).

留言 (0)

沒有登入
gif