Whole transcriptome sequencing data and corresponding clinical information (histology, gender, age, grade, therapeutic approaches, and survival information) were downloaded from the CGGA database (http://www.cgga.org.cn) as the training set. IDH mutations were detected using pyrosequencing and performed on a Pyro-Mark Q96 ID System (Qiagen, Valencia, CA, USA) as previously described (10). Primers used for PCR amplification are listed as follows: IDH1 5′-GCTTGTGAGTGGATGGGTAAAAC-3′ and 5′-Biotin-TTGCCAACATGACTTACTTGATC-3′; IDH2 5′-ATCCTGGGGGGGACTGTCTT-3′ and 5′-Biotin-CTCTCCACCCTGGCCTACCT-3′. The primer sequences used for pyrosequencing are 5′-TGGATGGGTAAAACCT-3′ for IDH1 and 5′-AGCCCATCACCATTG-3′ for IDH2. The status of chromosome 1p/19q was calculated as reported (11). A total of 362 LGG patients with clinical characteristics, transcriptome sequencing, and molecular data (ATRX, 1p/19q, IDH) were downloaded from The Cancer Genome Atlas (TCGA) as the validation set (http://cancergenome.nih.gov/). Whole-genome mRNA expression microarray data and clinical information from GSE16011 and Repository of Molecular Brain Neoplasia Data (Rembrandt) were obtained for analyses as well. Moreover, the Genotype-Tissue Expression (GTEx) transcriptome data from GEPIA2 was obtained to compare the expression between normal and tumor tissues.

The known and putative α-KG-dependent dioxygenases were obtained as previously reported, including 10 DNA/RNA-modifying enzymes (TET1, TET2, TET3, ABH1, ABH2, ABH3, ABH4, ABH5, ABH6, and FTO), 29 JmjC domain-containing enzymes (KDM2A, KDM2B, KDM3A, KDM3B, KDM4A, KDM4B, KDM4C, KDM4D, KDM5A, KDM5B, KDM5C, KDM5D, KDM6A, KDM6B, KDM7A, KDM8, HR, JARID2, JHDM1C, JMJD1C, JMJD4, JMJD6, JMJD7, JMJD8, MINA, NO66, PHF2, PHF8, and UTY), 15 proline/lysine hydroxylases (EGLN1, EGLN2, EGLN3, P4HA1, P4HA2, P4HA3, P4HB, P4HTM, PLOD1, PLOD2, PLOD3, LEPRE1, LEPREL1, LEPREL2, and TMLHE), and 11 other hydroxylases (ASPH, ASPHD1, ASPHD2, BBOX1, FIH1, HSPBAP1, OGFOD1, OGFOD2, PAHX-AP1, PHYH, and PHYHD1) (12).

To identify differential genes between IDH wild type (IDH^WT) and IDH mutant (IDH^MUT) patients, Student’s t-test was applied, and genes with corresponding p< 0.05 were considered as candidates. The diagnostic performance of the IDH-related gene signature was measured using generalized linear models (GLMs) embedded in Matlab software. The GLM generalizes linear regression by allowing the linear model to be linked to the response variable and by allowing the magnitude of the variance of each measurement to be a function of its predicted value. A receiver operating characteristic (ROC) curve with the calculated area under the curve (AUC) was delineated based on the screened genes.

According to Harrell’s guideline, the number of samples should exceed at least 10 times of included variables in a multivariate analysis. To address this issue, we performed the Lasso Cox regression model for dimensionality reduction (13). The selected genes were used in developing a linear combination weighted by their respective coefficients generated by the Lasso Cox model. The risk score for the OS time of each individual was calculated as follows:

Next, the patients in the training dataset were classified into high-risk and low-risk groups using the median risk score as the cutoff point. The same coefficients for each gene and median risk score cutoff were applied to TCGA, GSE16011, and Rembrandt datasets for validation.

To classify the latent molecular alterations between high- and low-risk patients, we screened the genes significantly associated with the evaluation of risk scores. The gene with Pearson’s correlation coefficient greater than 0.5 and a p-value less than 0.05 were subjected to Gene Ontology (GO) analysis using the online Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov/) (14). Heatmaps were constructed using Gene Cluster and Gene Tree View software. Top significant biological processes or cellular components were figured with R programming language (http://cran.r-project.org). Meanwhile, we explored the differences between risk groups utilizing gene set enrichment analysis (GSEA). The annotated gene sets were obtained from the Molecular Signatures Database v5.1 (MSigDB) (http://www.broad.mit.edu/gsea/msigdb/).

To identify therapeutic approaches that can be used to specifically treat LGG patients with a worse risk score, genetic and pharmacological profiles were obtained from the Catalogue of Somatic Mutations in Cancer (COSMIC, http://cancer.sanger.ac.uk/cell_lines) (15). Pearson’s correlation algorithm was utilized to identify the correlation between the IC₅₀ value of anti-cancer drugs and gene expression data.

The significant differences between the two groups were estimated using Student’s t-test. The chi-square test and Fisher’s exact test were used to compare the frequencies between groups. The Kaplan–Meier survival analysis was implemented to estimate the survival distributions. Cox regression was used to determine the prognostic value of each variable with OS in LGG patients. Nomograms were depicted based on the results of a backward step-down selection process with the Akaike information criterion (AIC) by using the package of rms in R (13, 16). The performance of the nomogram was measured by the concordance index (C-index), which was indicative of the accuracy of the prognostic prediction (16). Calibration curves were drawn for each dataset for predicted 1-, 2-, and 3-year survival. All differences were considered statistically significant at the level of two-sided p< 0.05. Statistics were carried out using the SPSS 13.0 statistical software and GraphPad Prism 6.0.

To evaluate the relationship between IDH mutation and enzyme-encoding genes in LGG patients, we first screened the differentially expressed genes between the IDH^WT and IDH^MUT groups. A total of 65 known and putative α-KG-dependent dioxygenases were utilized for analyses (12). After genes with no significance were eliminated, an enzyme-related signature of 41 genes was established (Figure 1A). Among which, 26 were upregulated and 15 were downregulated in IDH^MUT LGGs compared to the IDH^WT counterpart. Then, we applied the GLM algorithm to access the diagnostic value of these genes in discriminating IDH phenotype. In the CGGA cohort, these 41 genes could distinguish LGG patients as IDH^WT and IDH^MUT with an AUC of 0.973 (Figure 1B), while this signature could stratify samples with an AUC of 0.992 in TCGA cohort (Figure 1C).

Previous studies highlight the importance of IDH mutation in prognostic prediction for LGG patients. To determine the prognostic significance of these differential genes, the Lasso Cox algorithm was used to pick the most potent prognostic characteristics from the CGGA dataset (Figures 2A, B). The genes with a non-zero coefficient in the Lasso Cox regression model were TET1 (−0.418), TMLHE (−0.179), ASPH (−0.142), EGLN3 (−0.031), EGLN1 (−0.013), PLOD3 (0.004), P4HB (0.006), KDM5A (0.080), LEPRE1 (0.235), and PHYH (0.559) (Figure 2C). These 10 genes presented a similar correlation in both the CGGA and TCGA cohorts (Figure 2D) and were significantly dysregulated between different IDH phenotypes, WHO grades, and chromosome 1p/19q status (Figures 2E-G). EGLN1, KDM5A, P4HB, and PLOD3 were statistically significant between normal and LGG tissues (Figure S1). Moreover, these genes could serve as predictors for OS of LGG patients (Figures S2, S3). To access the general efficiency and effectiveness of the 10 genes comprised signature in prognostic prediction, we calculated the risk score of each individual based on the gene expression levels and corresponding regression coefficients. Patients were classified into a high-risk group (n = 86) and a low-risk group (n = 86) based on the median risk value as the cutoff. Notably, patients with lower risk score survived longer than those in the high-risk group (hazard ratio [HR] = 0.0983, 95% confidence interval [CI]: 0.055–0.176, p< 0.0001; Figure 3A). Furthermore, this signature-based risk score could divide WHO II and III patients into two groups with the distinct OS using the same cutoff (Figures 3B, C). Similarly, in different subtypes of LGG, this signature could be used for the prediction of prognosis (Figures 3D-F). Cox regression analysis revealed that this risk model was associated with the survival of patients, independent of the clinical parameters (HR = 2.642, 95% CI: 1.612–4.332, p< 0.001; Table 1).

The prognostic value of the α-KG-dependent enzymes was validated using TCGA RNA sequencing data. The same formula was applied to calculate the risk score of each sample in the validating datasets, and the median risk score was used as a grouping criterion. In TCGA cohort, the survival curve showed that LGG patients with high-risk scores had a shorter OS than low-risk ones (p< 0.0001; Figure S4). Cox regression analysis revealed that the signature-based risk score could serve as an independent predictive indicator for LGG patients (p = 0.007; Table 1). To further confirm the robustness of this risk signature, microarray data from GSE16011 and Rembrandt datasets were analyzed. Consistently, the OS of LGG patients in GSE16011 and Rembrandt dataset could be stratified distinctly (Figures S4, S5), while the multivariate Cox model revealed that this risk signature could serve as an independent indicator for LGG patients (GSE16011, HR = 1.965, 95% CI: 1.475–2.618, p< 0.001; Rembrandt, HR = 2.890, 95% CI: 1.789–4.672, p< 0.001; Table S1).

Considering the significance of the risk score in prognostic prediction, we attempted to explore the latent molecular alterations. GO analysis of the high-risk positively related genes revealed that several malignant biological processes, including angiogenesis, cell adhesion, inflammatory response, and extracellular matrix organization, were significantly enriched (Figures 4A, B). Similarly, the GO results of TCGA, GSE16011, and Rembrandt suggested that the biological processes were mainly in relation to cell adhesion, angiogenesis, inflammatory response, apoptotic process, and cell proliferation (Figure S6). Afterward, the differences in clinical and pathological characteristics between high- and low-risk populations were compared. There are more WHO grade 3 gliomas, lower IDH mutation frequency, and lower 1p/19q codeletion ratio in the high-risk group in both the CGGA and TCGA datasets (Table 2). In GSE16011 sets, the low-risk group showed a higher IDH mutation frequency, while high-risk patients in Rembrandt sets were diagnosed with more WHO grade 3 glioma (Figure S7). GSEA identified that significant hallmarks involved in the high-risk group were inflammatory response, epithelial–mesenchymal transition, glycolysis, and angiogenesis (Figure 4C).

The prognostic nomogram that integrated all significant factors derived from AIC selection was depicted. In the CGGA cohort, the integrated nomogram with an enzyme-based risk signature could increase the C-index from 0.803 to 0.856 (Figure 5A). The calibration curves for the probability of survival at 1, 2, or 3 years showed an optimal agreement between the estimation and actual observation (Figure 5B). Moreover, this nomogram was independently validated in TCGA cohort. The C-index for validating the nomogram is 0.754, which showed a good calibration as well (Figure 5C).

To explore the therapeutic strategies for the high-risk population, the pharmacogenomic interaction dates of LGG cell lines were analyzed. With the exacerbation of risk score, LGG cell lines presented a higher IC₅₀ value for a spectrum of anti-cancer drugs, such as tozasertib, paclitaxel, and sorafenib. On the contrary, PI3K/Akt targeted drugs, including pictilisib (PI3K), MK-2206 (Akt1, Akt2), NU7441 (DNAPK, PI3K), AZD6482 (PI3K), and idelalisib (PI3K), effectively inhibited cell viability in vitro, which is attributed to the hyperactivation of PI3K/Akt signaling pathway in high-risk tumors (Figure S8). Interestingly, the IC₅₀ of JQ-1, a BET bromodomain inhibitor, was negatively correlated with the enzyme-based risk score (Figure S8D), suggesting a therapeutic implication for epigenetic inhibition in the treatment of this population.