mRNA expression data with clinical information of 506 and 596 LGG patients were collected from TCGA (training dataset) and CGGA (validation dataset), respectively. The corresponding clinical data include: gender, age, histology, subtype, IDH1 mutational status, MGMT promoter status, 1p/19q co-deletion status and survival information. Somatic mutation data were obtained from the TCGA portal. Considering that the GLASS dataset (https://www.glass-consortium.org) contained only the minimal clinical information (n = 330) [25], the GLASS dataset only served as an additional validation dataset to validate a portion of results. A total of 150 genes (Supplementary Table 1) associated with DDR were retrieved from Molecular Signature Database V7.0 (MSigDB) [26] (http://www.broad.mit.edu/gsea/msigdb/, HALLMARK_DNA_REPAIR).

The mRNA expression data for samples were pre-processed according to the following steps. Samples with missing clinical information and those with an overall survival (OS) of less than 30 days were excluded, and 477 LGG samples from TCGA were eventually included in this study for follow-up analysis (Supplementary Table 1). In addition, the GLASS dataset was used as an independent validation dataset, in which samples from TCGA were excluded.

The DDR-related genes shared among the TCGA, CGGA, and GLASS datasets were retained for further analysis. Consequently, the final dataset comprised expression profiles of these common DDR-related genes derived from the aforementioned three datasets.

The CGGA mRNA sequencing data included both mRNAseq_693 (Illumina HiSeq platform) and mRNAseq_325 (Illumina HiSeq 2000 and 2500 platforms). Therefore, the batch effects between these two datasets were removed using the “sva” R package [27] and a total of 594 LGG patients in CGGA were retrieved (Supplementary Table 1).

In order to identify gene signature and construct risk score, the R package (http://cran.r-project.org) “survival” and “survminer” were used to perform univariate Cox proportional hazard regression on 477 LGGs in TCGA dataset to select DDR-related genes that are highly correlated with prognosis (p 0.05 was selected as the threshold). Then “randomForest” R package was performed to ranked the importance of the selected DDR-related genes related to the prognosis by Random Survival Forest algorithm (nrep and nstep respectively represent the number of Monte Carlo simulation iterations and the number of steps forward, with values of 100 and 5). Genes with the relatively importance value $>$0.2 were subsequently incorporated into the multivariate stepwise cox regression analysis to established a prognostic signature. In three datasets, the Cox regression coefficient and prognostic gene expression were used to calculate the risk score for each patient. Finally, the risk model containing 6 genes was established, and the final formula was as follows:

Risk score$={\sum }_{i=1}^N(\beta \text{i}\times \text{Expi})$, where $\beta$i represented the coefficient of gene i, and Expi represented the expression value of gene i.

According to the risk score formula, samples from the training and validation datasets were divided into high-risk and low-risk groups based on the median risk scores. Kaplan-Meier (KM) curves were used to estimate the survival divergence between different risk groups. To determine the independent predictive power of each variable, we conducted multivariate and univariate regression analyses using R “survival” package and visualized by the “forestplot” package in R. Additionally, we retrieved the gene protein expression level from the Human Protein Atlas (HPA) database to evaluate the expression of prognostic genes.

A nomogram is a graphical representation of a mathematical model that estimates the outcome or probability of an event based on multiple variables or inputs, serving as a visual tool to simplify complex calculations. Nomograms facilitate clinicians’ decision-making processes, enabling them to tailor treatment strategies for individual patients according to their unique characteristics and risk factors. Univariate and multivariate Cox regression were conducted to identify the independent prognostic risk factors for nomogram construction. Using the “rms” packages, calibration curves were drawn to assess the OS probability for LGG patients at 1-, 3-, and 5-years. In addition, nomogram accuracy was also determined by concordance index (C index).

Gene set enrichment analysis (GSEA) was performed using the GSEA V3 software (Broad Institute, Inc., Massachusetts Institute of Technology, and Regents of the University of California) [28] to identify gene sets that were statistically different between different risk groups in LGG samples of TCGA and CGGA mRNA expression profiles. GSEA V3 software was used to select the MSigDB “c2.cp.kegg.v7.1.symbols.gmt” gene set as a reference to screen out significantly enriched pathways with NOM p-value and false discovery rate (FDR) less than 0.05 and 0.25, respectively.

We consulted the recently updated Connectivity Map (CMap, http://cmap.topcoder.com/) to screen for potential small molecule compounds targeting DDR-related genes in the hope of achieving new breakthroughs in LGG therapy. “limma” R package [29] was applied to identify differential expression genes (DEGs) with FDR 0.05 and $|$Fold change (FC)$|$ $\ge$2 between different risk score groups. The up- and down-regulated genes were uploaded to the CMap (https://clue.io) database [30]. Small molecule compounds with an absolute value of enrichment greater than or equal to 0.7 and p 0.05 were considered as candidates.

“ConsensusClusterPlus” [31] R package was performed for consensus clustering of the prognostic DDR-related genes. Similarity distance among samples was calculated using Euclidean distance and clustering was conducted using K-means. Resampling analysis was conducted to estimate the optimal cluster number based on cumulative distribution function (CDF) and consensus matrix. And the clustering results determined by the above methods were also validated in GLASS dataset. The principal component analysis (PCA) was applied to evaluate the distribution of the molecular subtypes in TCGA and GLASS datasets via the “princomp” package. The differences in OS among different clusters in TCGA and GLASS datasets were present using KM curves. We further estimated single nucleotide polymorphism (SNP) distribution of each cluster. The somatic mutation profiles of LGG were firstly obtained from TCGA data portal and were divided into distinct clusters for analysis according to consensus clustering. The somatic mutation data of each cluster were summarized, analyzed, and visualized based on the Mutation Annotation format (MAF) file using the “maftool” [32] R package, and the driver genes were identified. The copy number alterations (CNA) data in each cluster were also generated by GISTIC2.0 from GDAC Firehose (https://gdac.broadinstitute.org) [33].

782 marker genes for predicting the abundance of 28 tumor-infiltrating immune cells (TIICs) were obtained from Broad Institute and Bindea et al. [34]. According to the characteristics of expressed genes, gene set variation analysis (GSVA) software “gsva” [35] R package was used for single-sample gene set enrichment analysis (ssGSEA), and the enrichment scores of 28 different immune cell types in each LGG sample were quantitatively analyzed. The levels and functions of immune infiltration were compared using Kruskal-Wallis test, and heatmaps and boxplots were used to visualize the results. In addition, the immune score and matrix score were calculated through the ESTIMATE algorithm [36].

The chemotherapy response of LGG samples in TCGA and GLASS were predicted using the Genomics of Drug Sensitivity in Cancer (GDSC) database [37]. GDSC is a large-scale research project that aims to identify genomic biomarkers associated with the response of cancer cells to various anti-cancer drugs. We assessed the sensitivity of eight chemotherapy drugs (Cisplatin, Erlotinib, Methotrexate, Vincristine, Carmustine, TMZ, Rapamycin and Doxorubicin) in gliomas. By using the “pRRophetic” [38] R package we estimated LGG’s half maximal inhibitory concentration (IC50) by ridge regression, and 10-fold cross-validation was performed to determine the accuracy.

To construct a DDR-related score, the “limma” [29] R package was applied to filter DEGs among DDR clusters (adjusted p 0.05 and $|$log2FC$|$ $>$1). Boruta algorithm was used for feature selection based on DEGs, and then feature genes were analyzed by univariate Cox regression. The Boruta algorithm is a feature selection method used in machine learning and data analysis to identify the most relevant and informative features in a dataset. It is based on the Random Forest algorithm, a powerful ensemble learning technique that constructs multiple decision trees to make predictions or classifications. Further analysis was performed on the genes with the significant prognosis. We then performed PCA to establish DDR related gene signature. Both the principal component 1 and 2 were selected to serve as signature scores. We then define the DDR score similar to genomic grade index [39, 40]: DDR score = Ʃ (PC1i + PC2i), where i represent the expression of the prognostic genes. We subsequently validated the prognostic value of DDR score in the disease special survival (DSS), progression free survival (PFS) and GLASS dataset, respectively. In addition, we systematically searched the gene expression profiles of DDR therapy, and a advanced urothelial cancer with intervention of atezolizumab (IMvigor210 cohort) was included in our analysis [41].

Statistical analysis and tests were implemented through R software version 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria). One-way ANOVA with Tukey’s test (q-value 0.05, ǀlog2FCǀ $>$1) was implemented to determine DEGs between different risk groups or molecular subgroups by employing the “limma” R package. Wilcoxon rank sum test was applied for continuous data of two groups, while Kruskal-Wallis test was used for continuous data of more than two groups. And the Chi-Squared Test was used to the categorical data. The optimal cut-off point of each subgroup in each dataset was determined by the “surv_cutpoint” function in the “survminer” R package. Multivariate and univariate Cox regression analyses were utilized to assess association between OS and clinicopathological characteristics and risk scores. Survival curves among different clusters were generated using KM curves with log-rank tests. To evaluate the predictive efficacy of the risk score, Receiver Operating Characteristic (ROC) and area under the curve (AUC) were delineated using the “pROC” R package. p-value was corrected by Bonferroni’s test. Statistical significance was defined as two-sided p 0.05.

An LGG prognostic risk model was constructed based on DDR gene in training set. Firstly, the relationship between 150 DDR-related genes and survival was analyzed by univariate Cox regression and 87 DDR-related genes that were associated with LGG survival (p 0.05) were identified (Supplementary Table 2). Further, random forest algorithm (“randomForestSRC” R software package) was used for feature selection and the importance of survival-related DDR genes was ranked. Finally, the genes with relative importance greater than 0.2 were selected as the final signature (as shown in Supplementary Fig. 1A,B). We subsequently performed multivariate stepwise cox regression analysis on the seven genes identified by the random forest to established the prognostic signature, and the importance and relative importance of the Hazard ratios (HRs), coefficients, confidence intervals and out-of-bag estimates of the six genes are shown in Table 1. The gene expression levels used for calculating the Risk score were absolute values, and the unit of each gene expression level was transcripts per million (TPM). The risk score formula was as follows:

Risk Score = TARBP2 $\times$.512574212 + STX3 $\times$ 0.300994954 + POLD3 $\times$ 0.748353899 + SDCBp $\times$ 0.433562531 + POLR3GL $\times$ (-)1.10017873 + POLL $\times$ (-)1.114708732

As shown in Fig. 1, there were more deaths in TCGA cohort with a high-risk score than low-risk score, while the heatmap (Fig. 1A) showed that changes in TARBP2, POLR3GL, POLL, STX3, POLD3 and SDCBP gene expression increased the risk value in different samples. KM curve analysis results indicated that the prognosis of patients in different risk groups was significantly distinct (p 0.00001), and patients in the low-risk group had a better prognosis than those in the high-risk group (Fig. 1B). We also evaluated the accuracy of risk score to predict prognosis. We performed the R software package “SurvivalROC” to conduct ROC analysis on the prognostic classification of risk score and evaluate the accuracy of risk signature. The AUC value of the 1-, 3- and 5-year in TCGA dataset were 0.901, 0.832 and 0.771 (Fig. 1C).

Fig. 1.

Construction of prognostic gene signature in the Cancer Genome Atlas (TCGA) dataset. (A) Risk score, survival time, survival state and expression of the 6 genes in TCGA training dataset. The panel consists of three rows: top row showed the risk score distribution for the high- and low-risk score group, color transition from green to red indicates the increasing risk score of corresponding expression level of DNA damage repair (DDR) genes from low to high; middle row represents the low-grade gliomas (LGG) patients’ distribution and survival status; the bottom row shows that the heatmap of 6 prognostic DDR-related genes expression. (B) Kaplan–Meier curves of prognostic models for high- (n = 238) and low-risk (n = 239) subgroups of patients with LGG in TCGA datasets. (C) Receiver Operating Characteristic (ROC) curve and area under the curve (AUC) of the 6-gene signature classification.

To validate the survival risk score, we calculated the risk score of patients in CGGA and GLASS validation datasets with the same regression coefficient and repeated the above analysis, and as expected, the validation datasets also achieved consistent results (Supplementary Fig. 2A–C and Supplementary Fig. 3A,B). We also compared our signature with other previously reported risk models [42, 43, 44, 45], the result showed that our signature had excellent AUC value when compared to the other signature, indicating that our risk model is the optimal (Supplementary Fig. 4). In addition, we also evaluated the protein expression level through the HPA database and found that TARBP2 and POLD3 were highly expressed in the tumor tissue, while POLR3GL and SDCBP showed low expression in normal tissue (Supplementary Fig. 5).

The association between risk score and major clinical characteristics in LGG was investigated. We found a significant divergence of risk score between different IDH1, MGMT promoter methylation, 1p/19q deletion, and pathological grades in TCGA and CGGA datasets. Patients with IDH1-wildtype (WT), 1p/19q non-codel, G3 grade and unmethylated MGMT promoter in LGG had significantly higher risk scores than those with IDH1-Mutant (Fig. 2C), 1p/19q codel (Fig. 2D), G2 grade (Fig. 2F) and methylated MGMT promoter (Fig. 2E) in TCGA and CGGA datasets (Supplementary Fig. 6C–F). The above results were similar to our previous studies [46].

Fig. 2.

Comparisons of risk score plot between age (A), gender (B), IDH (C), 1p/19q (D), MGMT (E), and grade (F) in TCGA dataset. IDH, isocitrate dehydrogenase; MGMT, O6-methylguanine-DNA methyltransferase.

However, differences in risk scores were not consistent across age stratification and gender in both datasets. In the TCGA dataset, patients over 50 years old had significantly higher risk scores than those aged 50 years or younger (Fig. 2A), and there was no difference between male and female patients (Fig. 2B). However, in CGGA dataset, patients of different ages (Supplementary Fig. 6A) and genders (Supplementary Fig. 6B) did not have significantly different risk scores. Possibly, this is due to differences in basic clinical characteristics between the two datasets, such as bias or ethnicity.

In addition, we validated the prognostic differences between high- and low-risk group in stratified cohorts of different clinical characteristics (e.g., IDH1-WT/-Mutant, MGMT promoter methylated/unmethylated and 1p/19q codel/non-codel cohorts). Results were consistent with expectations. In the stratified cohorts of TCGA (Fig. 3) and CGGA (Supplementary Fig. 7) datasets with different clinical characteristics, except for the IDH1-WT (Fig. 3F) group of TCGA, grade G2 (Supplementary Fig. 7K) and IDH1-WT (Supplementary Fig. 7F) of CGGA, the high-risk groups showed poor prognosis. A possible reason for this due to the small sample size of the IDH1-WT group in TCGA and CGGA.

Fig. 3.

Kaplan-Meier survival curves for the high- and low-risk groups stratified by clinical factors in TCGA dataset. (A,B) Age. (C,D) Gender. (E,F) IDH. (G,H) 1p/19q. (I,J) MGMT. (K,L) Grade.

Clinical characteristics (age, gender, risk score, grade, IDH1, MGMT promoter methylated, 1p/19q codeletion) in TCGA and CGGA datasets were analyzed by conducting univariate and multivariate Cox regression. As shown in the forest plot (Fig. 4), the risk score constructed by the 6 DDR-related genes was an independent prognostic risk factor in TCGA (p 0.05). The independence was also validated in CGGA dataset (Supplementary Fig. 8).

Fig. 4.

DDR signature as an independent prognostic factor shows strong prognostic power. (A) Univariate and (B) multivariate Cox regression analysis of clinical characteristics and DDR signature for survival of LGG in TCGA training dataset.

To further improve the accuracy of prognosis, a nomogram integrating risk score and independent clinical characteristics was constructed to visualize the prognostic value of different survival-related variables in TCGA and CGGA datasets (Fig. 5A and Supplementary Fig. 9A). In addition, to assess the performance of nomogram, a calibration curve analysis was conducted (“rms” R package) and we observed that the prediction curve was close to the ideal curve (Fig. 5B–D and Supplementary Fig. 9B–D), indicating that the nomogram had excellent prediction function. We discovered that the nomogram showed a high accuracy in the prognosis of 1-, 3- and 5- year (AUC value $>$0.8) (Fig. 5E and Supplementary Fig. 9E). Interestingly, the prediction accuracy of nomogram was also higher than the other clinical factors (Fig. 5F and Supplementary Fig. 9F).

Fig. 5.

Construction of nomogram. (A) A nomogram constructed based on independent prognostic factors. The calibration plot for internal validation of the nomogram in 1- (B), 3- (C) and 5-years (D) OS. (E) ROC curve and AUC of the Nomogram (F) Comparisons of the ROC curve among the clinical factors and nomogram. Abbreviations: OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve (AUC). * p 0.05, *** p 0.001.

GSEA provides valuable insights into the biological processes and pathways that are potentially dysregulated in the context of the studied conditions, helping researchers understand the underlying molecular mechanisms and guiding the development of targeted therapies or interventions. GSEA analysis of different risk groups in the TCGA cohort revealed different biological processes. Significant candidate pathways were selected according to the criteria with FDR 0.25 and nominal p-value 0.05. GSEA showed that cancer-related signaling pathways including focal adhesion, glioma, oocyte meiosis and apoptosis were significantly enriched in high-risk patients, while none of these pathways were observed in low-risk patients (Supplementary Fig. 10).

In order to identify candidate compounds that may potentially influence the expression of DDR prognostic genes for the treatment of LGG, DEGs between different risk groups (191 up-regulated genes and 81 down-regulated genes, “limma” R package, Supplementary Table 1) were uploaded to the CMap database. As a result, 35 compounds with 25 mechanisms of action (MoA) were identified (Fig. 6). Results indicated that multiple compounds shared the same MOA. For example, six compounds (scriptaid, ISOX, vorinostat, belinostat, THM-I-94 and HC-toxin) shared the MOA of histone deacetylase (HDAC) inhibitor; four compounds including amsacrine, irinotecan, teniposide and SN-38 shared the MoA of Topoisomerase inhibitor; another two compounds (NU-7441and KU-0060648) shared the MoA of DNA dependent protein kinase inhibitor. In addition, a total of five compounds (TW-37, butein, everolimus, NVP-TAE684 and XMD-892) shared the following five tumor-related MoAs, respectively: B cell lymphoma 2 (BCL) inhibitor, epidermal growth factor receptor (EGFR) inhibitor, mammalian target of rapamycin (MTOR) inhibitor, anaplastic lymphoma kinase (ALK) inhibitor and mitogen-activated protein (MAP) kinase inhibitor. The present research identified potential small molecule compounds that target DDR.

Fig. 6.

Identification of potential molecular compounds and its MoAs in LGG. Heatmap showed each compound from the CMap that shares a MoA (rows), sorted by descending number of compounds with a shared MoA. The above compounds have an absolute value of enrichment score $\ge$0.7 and might be capable of targeting the DDR-related signature. Abbreviations: MoA, mechanisms of action; CMap, Connectivity Map; BCL, B cell lymphoma; HDAC, histone deacetylase; EGFR, epidermal growth factor receptor; MTOR, mammalian target of rapamycin; ALK, anaplastic lymphoma kinase; MAP, mitogen-activated protein; FLT3, Fms-like tyrosine kinase 3; HIF, Hypoxia-Inducible Factor; IKK, IkappaB kinase; KIT, tyrosine kinase.

We calculated the enrichment levels of 28 immune-related cells for each risk group using ssGSEA algorithm. We found that most immune-related cells have a significant divergence between different risk groups in TCGA and CGGA dataset, such as Memory B cell, Regulatory T cell, Eosinophil, Macrophage, Natural killer T cell. Interestingly, the enrichment level was significantly higher in high-risk group. (Fig. 7 and Supplementary Fig. 11).

Fig. 7.

The landscape of 28 immune cell infiltration level in high-risk (n = 238) and low-risk (n = 239) groups in TCGA dataset. ns, not significant; * p 0.05; ** p 0.01; *** p 0.001.

Consensus clustering was conducted for clustering analysis of 477 samples and 150 DDR-related genes expression profiles in TCGA RNA-seq dataset. All LGG samples were divided into K clusters (K = 2–9) using the “ConsesusClusterPlus” R package. Based on the CDF curve, 477 samples were divided into two clusters after a thorough consideration (Fig. 8A–C). OS significantly differed between the two clusters based on the KM curve analysis (Fig. 8D). PCA results showed that the two clusters were independent of each other (Fig. 8E).

Fig. 8.

DDR-related genes could distinguish LGG patients in TCGA with different clinical and molecular features. (A) Consensus clustering CDF for k = 2 to k = 9. (B) Delta area plot exhibits relative change in area under CDF curve when at a certain k value compared with k-1. (C) Consensus clustering matrix heatmap plots of 477 samples from TCGA datasets for k = 2. (D) KM analysis of patients among two clusters,the prognosis of Cluster A (n = 361) is better than Cluster B (n = 116). (E) PCA analysis of the DDR-related genes expression when k = 2. Abbreviations: CDF, cumulative distribution function; KM, Kaplan-Meier.

The above analyses were repeated in GLASS dataset and the results obtained were shown in Supplementary Fig. 12. Heatmap of the two clusters (the top 100 variably expressed genes) showed that a significant divergence existed in the IDH1 status, 1p/19q co-deletion status and MGMT methylation status (Fig. 9). Furthermore, we compared the risk score differences among the two different clusters in TCGA. Interestingly, as shown in Supplementary Fig. 13A, Cluster B with the worse prognosis had the higher risk scores. Conversely, Cluster A with the better prognosis had the lower risk scores. Therefore, we identified two LGG clusters based on DDR gene expression. Clustering results showed that there was a wide range of heterogeneity in LGG, and DDR-related genes were closely associated with LGG malignancy grade and prognosis.

Fig. 9.

The landscape of DDR-related genes among the two subtypes (Cluster A, n = 361; Cluster B, n = 116). Heatmap of two clusters defined by the top 100 variable expression genes.

We further showed the mutant landscape in different clusters of LGG patients. By analyzing the somatic mutation data in TCGA dataset, we evaluated the somatic variation distribution of driver genes among LGG subtypes, to clarify whether there are differences in the frequency of somatic mutation and observe the different patterns of mutations among LGG clusters. LGG driver genes were obtained by Maftools. The analysis of TCGA mutation annotation files showed that there were obvious mutation characteristics between Cluster A and Cluster B. Among them, IDH1, TP53 and ATRX were the top three mutation genes and Cluster A corresponded to more mutation frequency than Cluster B (Supplementary Fig. 14A,B). Our cluster-specific CNAs analysis revealed that chromosome deletions and amplifications, for example, deletion of 9p21.3 were significantly enriched in Cluster B. The 9p21.3 deletions is associated with the poor survival outcome in LGG [47] (Supplementary Fig. 15).

The scores of different clusters were compared, as shown in Fig. 10A–C, Cluster B with the worst prognosis in TCGA dataset had the higher immune and stromal scores compared to Cluster A, while tumor purity was lower than the other cluster (ANOVA, p 0.001). Cluster A with the better prognosis had the opposite result, with the lower immune and stromal scores and the higher tumor purity.

Fig. 10.

The landscape of immune infiltration and TME characteristics in two subtypes. Violin diagram showed the differences in stromal score (A), immune score (B) and tumor purity (C) among 2 clusters. (D) Unsupervised clustering of LGG patients from the TCGA cohort using ssGSEA scores which represent the 28 cell infiltration. The cluster, age, survival status, gender, grade, IDH1 status, 1p/19q status and MGMT status are shown as patient annotations in the lower panel.

The correlation between LGG subtypes and TME features was further explored in our present study. Results of ssGSEA showed that in LGG patients of TCGA cohort, the enrich scores of activated B cells (Ba), eosinophils and monocytes in Cluster B patients with worse prognosis were highly expressed compared to Cluster A. Fig. 10D shows the high expression level of immune cells in Cluster A.

Currently, surgery and chemotherapy are the most common treatments for cancer. Therefore, we sought to evaluate the differences in estimated IC50 levels between the clusters for eight chemotherapeutic agents, including Cisplatin, Erlotinib, Methotrexate, Vincristine, Carmustine, TMZ, Rapamycin and Doxorubicin. Ridge regression was used to train the prediction model on the dataset of GDSC cell lines, and 10-fold cross-validation was performed to evaluate prediction accuracy. The IC50 value of each LGG patient was then calculated according to the prediction model of the eight chemotherapeutic drugs. Results showed that except for Erlotinib (Wilcoxon rank test, p = 0.483 (Fig. 11D), the IC50 sensitivity of the seven chemotherapeutic agents was statistically different in the two clusters. Interestingly, the IC50 estimates for the seven chemotherapeutic agents in Cluster B were significantly lower than those in the other cluster, including Carmustine (Wilcoxon rank test, p = 1.927 $\times$ 10${}^{-5}$) (Fig. 11A), Cisplatin (Wilcoxon rank test, p = 1.985 $\times$ 10${}^{-11}$) (Fig. 11B), Doxorubicin (Wilcoxon rank test, p = 1.857 $\times$ 10${}^{-6}$) (Fig. 11C), TMZ (Wilcoxon rank test, p = 1.157 $\times$ 10${}^{-6}$) (Fig. 11E), Methotrexate (Wilcoxon rank test, p = 0.0492) (Fig. 11F), Vincristine (Wilcoxon rank test, p = 5.504 $\times$ 10${}^{-6}$) (Fig. 11G) and Rapamycin (Wilcoxon rank test, p = 1.587 $\times$ 10${}^{-7}$) (Fig. 11H), indicating that Cluster B was more sensitive to the above seven chemotherapeutic agents.

Fig. 11.

Cluster B is more sensitive to most chemotherapies. Box plots depicted the differences in the estimated IC50 levels of (A) Carmustine; (B) Cisplatin; (C) Doxorubicin; (D) Erlotinib; (E) Temozolomide; (F) Methotrexate; (G) Vincristine and (H) Rapamycin between two clusters.

A total of 574 DEGs were identified between the two clusters (Supplementary Fig. 16A,B). The GO enrichment analysis exhibited that these DEGs mainly enriched in chromosome segregation, nuclear division, and DNA replication (Supplementary Fig. 16C). The KEGG pathway enrichment analysis result was shown as Supplementary Fig. 16D.

We first performed the univariate cox regression analysis and selected the genes with p-value 0.05. We subsequently performed feature selection for significant DEGs exploiting Boruta algorithm. 95 genes from the feature selection were eventually applied to construct the DDR score based on the formula: DDRscore = Ʃ (PC1i + PC2i). The KM curve analysis result revealed that low DDR score corresponded to a better survival outcome, while high DDR score was associated with poor prognosis in TCGA dataset (Fig. 12A) and GLASS dataset (Supplementary Fig. 17). Moreover, the DDR score in Cluster B were significantly increased compared to Cluster A (Wilcoxon rank test p 2.2 $\times$ 10${}^{-16}$) (Supplementary Fig. 13B). The ROC analysis revealed that the DDR score had high AUC value in 1-, 3- and 5-years, suggesting that the DDR score could be used as potential biomarker to the prognosis of LGG (Fig. 12B). In addition, we also retrieved the progression free survival (PFS) and disease special survival (DSS) clinical info of LGG to validate the prognostic value of DDR score. As showed in Fig. 12C,E, the high DDR score group had a significantly better survival outcome compared to the low DDR score group in the PFS and DSS. The high AUC value also demonstrated that DDR score had a good prognostic value in PFS and DSS, respectively (Fig. 12D,F).

Fig. 12.

Construction and validation of the DDR score. (A) KM survival curve analysis in the overall survival (A), progression free survival (C) and disease special survival (DSS) (E). ROC analysis was applied to evaluate the specificity and sensitivity of the DDR score in overall survival (B), progression free survival (D) and DSS (F). Abbreviations: KM, Kaplan-Meier; OS, overall survival; DSS, disease special survival.

We further evaluated the response to immune checkpoint blockade (ICB) of DDR score in the IMvigor210 cohort (bladder cancer). Of note, survival outcome was significantly different from LGG, with low DDR scores associated with poorer survival outcomes (Supplementary Fig. 18).