The methods of this study are described in Figure 1 . To determine differentially expressed and methylated genes, we first extracted mRNA expression and DNA methylation profiles of 259 glioma samples with WHO grade I–III from TCGA database. Patients were divided into two subgroups according to the status of IDH mutation and 1p19q integrity. The clinical and pathological characteristics between subtypes are presented in Table S1 . Hierarchical bi-clustering was performed for IDH^mut/1p19q^codel samples (n = 165) and IDH^wt/1p19q^non-codel samples (n = 94). As a result, 137 candidates including 74 downregulated genes and 63 upregulated genes were selected ( Figure 2A and Table S2 ). Subsequently, the MethyMix method was used to filtrate DMDGs. A total of 433 DMDGs including 318 hypomethylated genes and 115 hypermethylated genes were determined, among which the adjusted p-value was less than 0.05 between the hyper- and hypomethylation groups and the correlation coefficient was less than −0.3 between DNA methylation and gene expression ( Figure 2B and Table S3 ). Afterwards, gene ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to elucidate DMDG functional property (p < 0.05). The results demonstrated that multiple inflammation and tumor progress-related GO terms and signaling pathways were significantly enriched in IDH^wt/1p19q^non-codel gliomas ( Figures S1A, B and Tables S4 , S5 ). The Venn diagram of DMDGs and DEGs revealed 31 DME genes, which were simultaneously hypomethylated and upregulated at the transcriptional level or hypermethylated and downregulated ( Figure 2C ).

To further narrow down the scope of the candidate DME genes, the least absolute shrinkage and selection operator (LASSO) regression was performed to select the most suitable predictive variables. To this end, six candidate DME genes (named risk gene cluster) from a further LASSO Cox regression model were selected based on minimum lambda with 10-fold cross-validation. These genes included the following: DNA Damage Inducible Transcript 4 Like (DDIT4L), Epithelial Membrane Protein 3 (EMP3), Mesenchyme Homeobox 2 (MEOX2), Ovarian Cancer Immunoreactive Antigen Domain Containing 2 (OCIAD2), Transforming Growth Factor Beta 2 (TGFB2), and Tumor Necrosis Factor Receptor Superfamily Member 12A (TNFRSF12A) ( Figures 3A, B ). The transcriptional risk score predictive model was developed by adding the mRNA expression level and relevant coefficient of each gene in the LASSO regression as follows: transcriptional risk score = 0.0350970 × DDIT4L mRNA expression + 0.1368395 × EMP3 mRNA expression + 0.0974575 × MEOX2 mRNA expression + 0.0723336 × OCAID2 mRNA expression + 0.0738469 × TGFB2 mRNA expression + 0.2045352 × TNFRSF12A mRNA expression. Positive coefficients of all genes in the LASSO regression suggested that mRNA high expression levels were correlated with poor OS in LGG patients and the Kaplan–Meier (K–M) analysis was performed to confirm the relationship between transcriptional risk score, risk gene cluster expressions, and OS. The OS of high transcriptional risk score or mRNA high expression group was significantly shorter ( Figures 3C–I ). Of note, the X-tile method was utilized to distinguish the optimal cutoff value. Additionally, principal component analysis (PCA) was performed to assess the distinguished accuracy based on the DME genes. Compared with the thirty-one DME gene expression levels, the contributing rate of the first principal component was observably promoted to 76.2% using the risk gene cluster expression levels ( Figures S2A, B ). Despite the fact that the contributing rate of the first principal component was also ascending based on DNA methylation levels or the combination of methylation and expression using risk gene cluster, the clinical feasibility was inconvenient ( Figures S2C, D ). Therefore, the transcriptional risk score model depending on the expression of risk gene cluster was adopted for further analysis. As shown in the risk factor association diagram ( Figure 3J ), the blue dots in the figure represented the surviving LGG patients while the red dots represented death, and the corresponding risk gene cluster expression profiles were visualized as a heatmap. The dotted line indicated that the optimal cutoff value of transcriptional risk score, with which all LGG patients were divided into two groups including 176 low transcriptional risk score samples and 83 high transcriptional risk score samples. The results showed that along with the increasing of the transcriptional risk score, the number of deaths gradually increased as well as the mRNA expression levels of the risk gene cluster, demonstrating that the patients in the high transcriptional risk score group exhibited more severe survival and higher risk of death. To further elucidate signaling pathways underlying our risk score model, we perform the Gene Set Enrichment Analysis (GSEA) in the two groups. As shown in Figure S3 , a wide range of signaling pathways related to the tumor immunological process were enriched in the high transcriptional risk score gliomas, indicating that signals from immune cells or its presence within the tumor might be crucial factors that affect progression and recurrence of glioma.

To confirm the consistency between methylation level and gene expression of the risk gene cluster in LGG patients, a difference analysis and a correlation analysis were performed. The results showed that the hypermethylation status and downregulation level of the risk gene cluster were coincidently enriched in the IDH^mut/1p19q^codel group ( Figures S4 , S5 ) and significant negative correlation between methylation and mRNA expression could be observed ( Figure S6 ). Figure S7 also demonstrated the relative hypermethylation in the IDH^mut/1p19^codel samples and the relative hypomethylation in IDH^wt/1p19q^non-codel samples. Taken together, these findings indicated a consistent tendency between DNA methylation and the expression of risk gene cluster and implied a potential prediction model for gliomas.

To verify that the expression of the risk gene cluster was induced by corresponding DNA methylation alterations rather than by copy number alterations (CNAs) or mutations, the cBioPortal database and COSMIC database were used to investigate the CNA or mutation levels of the risk gene cluster. As shown in Figure S8 , the risk gene cluster was not detected in either the top 20 CAN genes or mutation genes; moreover, the genetic mutation ratio of these genes was less than 3% in all glioma samples. These results indicated that the transcriptional regulation of the risk gene cluster was driven by DNA methylation alterations.

To provide further insight into CpG methylation, the preprocessed CpG site methylation value of the risk gene cluster was selected. The correlation coefficients between CpG sites and risk gene cluster expression were calculated ( Table S6 ). Subsequently, 53 CpG sites were initially filtrated based on univariate Cox regression analysis among 79 CpG sites ( Figure S9A ), and the relevant genomic information was presented using a heatmap ( Figure S9B ), which also presented a significant difference between the high and low transcriptional risk score groups ( Table S7 ). The results indicated extensive hypomethylation located mainly in the CpG island, CpG shelf, CpG shore, and open sea of CpG sites in the high transcriptional risk score group, which might contribute to the upregulation of the corresponding risk gene cluster. Moreover, the correlation between 53 CpG sites and de-/methyltransferase derived from the Molecular Signatures Database (MSigDB) with GO_DEMETHYLATION and GO_METHYLATION was tested ( Figure S9C ). |r| > 0.7 and an adjusted p-value < 0.05 were set as cutoff criteria for further filtration of the de-/methyltransferase, of which the annotation was diagrammatized to present positive and negative de-/methyltransferase using a Sankey diagram ( Figures S9D, E ). Moreover, the protein–protein interaction (PPI) was performed by GeneMANIA ( Figure S9F ). The results indicated that de-/methyltransferase could participate in the regulation of the CpG sites of the risk gene cluster. Figure S9G demonstrated that 53 CpG sites were mainly distributed in promoter regions (60.4%), which included 1500 bp upstream of the transcriptional start site (TSS 1500) (20.8%), TSS 200 (13.2%), the 5’-untranslated region (5’UTR) (9.4%), and the first exon (1stExon) (17.0%). Interestingly, the 48 (90.57%) of hypomethylation CpG sites were also mainly distributed in the promoter region (63.8%) among 53 CpG sites. Therefore, we can speculate that the upregulation of risk gene cluster expression could appear due to the hypomethylation of CpG sites in the promoter regions of the corresponding genes under the effect of de-/methyltransferase.

To further contrast with transcriptional risk score in the accurate prediction of patient outcomes, LASSO regression was used to narrow down the candidate CpG sites. As a result, 3 CpG sites, namely, cg03208951, cg23344780, and cg23545105, were selected based on minimum lambda with 10-fold cross-validation ( Figures S10A, B ). The CpG risk score predictive model was developed by adding the product of the CpG methylation level and relevant coefficient of CpG site in the LASSO regression as follows: CpG risk score = (−1.0703883 × cg03208951) + (−1.7594301 × cg23344780) + (−0.4950028 × cg23545105). Afterwards, the 6-gene transcriptional risk score model and 3-CpG methylation risk score model were further assessed by concordance index (C-index) and time-dependent receiver operating characteristic (ROC) analysis. The results demonstrated that the prediction accuracy of about 3 years of transcriptional risk score was higher than CpG risk score and the long-term prediction accuracy was almost identical ( Figures S10C–F and Table S8 ). Although the long-term prediction accuracy was not statistically different, the prognosis of LGG patients was different and some IDH^wt/1p19q^non-codel LGG patients experienced rapid recurrence in the short term (8). Therefore, the transcriptional risk score model was selected as the most efficient prediction methods and was used to establish the nomogram.

We further investigate the predictive accuracies of the transcriptional risk score indicating the pathological subtypes, therapy reaction, and patient survival by using TCGA database. The results showed that the transcriptional risk score was markedly elevated in IDH^wt/1p19q^non-codel gliomas ( Figure 4A ). Consistently, an increased transcriptional risk score could also be observed in WHO grade III gliomas compared to those with lower WHO grade ( Figure 4B ). Interestingly, we also found that the transcriptional risk score was significantly increased in anaplastic astrocytoma (AA), which is considered to be more malignant and undergoes transition to glioblastoma more frequently, compared to other pathological subgroups of LGG ( Figure 4C ). Moreover, the results of ROC analysis indicated that our transcriptional risk score had encouraging sensitivity and specificity for distinguishing IDH/1p19q subtypes, WHO grades, and particular pathology subtypes, especially for discrimination of astrocytoma from oligodendroglioma or mixed glioma (MG) ( Figure 4D ). However, no significant difference could be observed between MG and anaplastic oligodendroglioma (AO) ( Figures 4C, D ). Next, the predictive efficiency of our transcriptional risk score model in primary or long-term treatment gliomas was investigated. According to our data, the transcriptional risk score was significantly increased in advanced gliomas as opposed to the stable-remission ones ( Figure 4E ). Also, the ROC analysis indicated that our transcriptional risk score was significantly more efficient for treatment prediction than WHO grade or pathological classification in short-term outcomes, but the AUC had no significant statistical difference between the transcriptional risk score and IDH/1p19q ( Figure 4F ). The transcriptional risk score was also significantly increased in advanced gliomas as opposed to the stable-remission ones in long-term outcomes ( Figure 4G ). Moreover, the ROC analysis indicated that our transcriptional risk score was significantly more efficient for treatment prediction than WHO grade or IDH/1p19q in long-term outcomes, but the AUC had no significant statistical difference between transcriptional risk score and pathology ( Figure 4H ). Finally, we used ROC curve to demonstrate that the transcriptional risk score model was more efficient than the calculation based on the single indicator to patient survival ( Figure 4I ). Taken together, our transcriptional risk score model showed high sensitivity and specificity and could be used as a reliable prognostic prediction model in glioma.

It should be noted that this transcription risk score can only be evaluated after surgical resection and therefore cannot be used for pre-surgical evaluation of malignancies of LGG. Thus, in addition to the transcription risk model, a non-invasive pre-operation quantification method should be established. To this end, we used a radiomics method to explore the related features with transcription risk score and further estimate the risk score level for LGG patients in pre-surgery ( Figure 5A ). Preprocessed contrast-enhanced MR images of 85 patients with pathological diagnosis and continuous follow-up were used to identify the most correlated radiological features. All of the 107 radiological signatures with intra-class correlation coefficient >0.80 were enrolled to establish a radiomics risk score model. The predictive model was established by adding the product of the radiomics feature value and relevant coefficient of each radiomics features in the LASSO regression based on minimum lambda with 10-fold cross-validation ( Figures 5B, C ). Finally, 13 radiomics features were selected and the radiomics risk score was calculated as follows: Radiomics risk score = 0.1967277 × Voxel Volume + 0.0086076 × Mesh Volume + (−0.1419602 × Sphericity) + 0.1840789 × Maximum 2D Diameter Column + (−0.5034319 × Large Dependence High Gray Level Emphasis) + 0.6774184 × Inverse Difference Moment Normalized + (−0.0740082 × Inverse Variance) + (−0.6442577 × Cluster Prominence) + (−0.0165486 × Skewness) + 0.0521741 × Gray Level Non Uniformity_GLSZM + (-0.0176127 × Large Area High Gray Level Emphasis) + 0.3468764 × Zone Entropy + 0.1971442 × Strength. We then verified a suitable calibration using the calibration curve analysis. The solid straight line (the 45-degree line) showed an ideal prediction radiomics model, and the broken lines represented the observed radiomics model, in which a closer fit to the dashed line means a better prediction model, and the result showed a satisfying consequence of this model and indicated that the radiomics risk score had a more favorable fitting to the transcriptional risk score ( Figure 5D ). Importantly, the radiomics risk score showed a statistically significant negative correlation with the transcriptional risk score model ( Figures 5E, F ). Moreover, ROC analysis also indicated that the radiomics risk score model exhibited remarkably improved sensitivity and specificity compared to the usage of the single radiomics feature ( Figure 5G ). Collectively, we analyzed 85 MR post-contrast T1-weighted images of LGG patients and identified 13 transcription risk score-specific radiomic signatures, and these results demonstrated that the radiomics-dependent model could be used as a dependable method for pre-operational assessment of the transcriptional risk score.

To explore the consistency of the risk gene cluster in clinical samples, immunohistochemistry (IHC) staining was performed using 61 glioma samples to address the expression level of the risk gene cluster. The results indicated that the risk gene cluster was significantly enriched in the IDH^wt/1p19q^non-codel group compared with their corresponding counterparts ( Figures 6A, B ). Moreover, the K–M analysis indicated that the higher immunohistochemical score (IHS) of the risk gene cluster correlated with poor prognosis for glioma patients, which is consistent with the results obtained from TCGA database ( Figure 6C ). The quantitative reverse transcription PCR (qRT-PCR) analysis was also performed to quantify the expression levels of the risk gene cluster by using 37 glioma samples with complete radiographic data and survival data (22 IDH^mut/1p19q^codel samples and 15 IDH^wt/1p19q^non-codel samples) and a non-tumor tissue derived from epilepsy patient used as a control. Consistently, the results showed that the risk gene cluster was significantly elevated in the IDH^wt/1p19q^non-codel group ( Figures 7A–F ). Similar to the results from IHC, K–M analysis also demonstrated prolonged OS for the patients with lower expression of these risk genes ( Figures 7G–M ). To further validate the reliability of our model, the transcription risk score model and radiomics model were calculated using clinical samples. The transcriptional risk score exhibited satisfying AUCs of 1-, 2- and 3-year OS (0.717, 0.802, and 0.923, Figure 7N ) based on qRT-PCR analysis, and the radiomics score also presented appropriate AUC (0.706) based on pre-surgical MRIs ( Figure 7O ). Of note, the cutoff values were recalculated by X-title. Taken together, these results indicated that our transcriptome model and radiomics model showed reasonably good reliability in our clinical cohort.

Based on our previous results, the transcriptional risk score was the most appropriate and accurate method compared to the others. Therefore, we adopted transcriptional risk score combined with clinical indicators to establish a novel nomogram for pre-surgical assessment of patient survival and therapy reaction. To this end, univariate Cox regression followed by multivariate Cox regression were performed to identify the most significant independent risk/protective factors. As a result, transcriptional risk score presented the most significant hazard ratio (HR) (HR = 2.94, 95% CI: 1.60–5.42, p < 0.01); in addition, patient age was also confirmed to be an independent risk factor with an HR of 2.55 (95% CI: 1.76–3.71, p < 0.001) ( Table 1 ). Therefore, patient age was enrolled in addition to the transcriptional risk score to establish the prediction nomogram.

Given that the age is well known to affect the methylation status of genes (26), the interaction between age and transcriptional risk score was verified via the interaction test. The result showed that the interaction was statistically significant, indicating that patient age might affect the expression of the risk gene cluster ( Table 2 , p for interaction = 0.033). Therefore, we performed further stratified analyses to eliminate this ambiguous association ( Table 2 ). Two hundred and fifty-nine patients were divided into 4 subgroups according to quartile categories of age, and the transcriptional risk score was divided into 3 subgroups according to tertile categories of risk score: Q1 (17–37 years old), Q2 (38–48 years old), Q3 (49–58 years old), and Q4 (59–87 years old), and low risk (Q1, transcriptional risk score: 0.3376742–0.6905997), median risk (Q2, transcriptional risk score: 0.6943530–1.6467377), and high risk (Q3, transcriptional risk score: 1.6655889–3.4480996); the median of each subgroup was used for statistical comparison ( Table 2 ). Significant differences were observed in all age subgroups (total HR: 2.371, 7.279, 5.285, and 2.078, p for trend: 0.038, <0.001, <0.001, and 0.002). The results demonstrated that the mortality risk of LGG patients was gradually elevated along with the increase of transcriptional risk score in each age subgroup. In particular, we found that the HR of transcriptional risk score showed an inverted U-shaped distribution along with the increase of age with the peak value appearing in the Q2 subgroup ( Figure S11A ). Collectively, these data suggested that enrichment of the risk gene cluster implied the highest risk of death of LGG in 35–45 years.

Next, we used variance inflation factor (VIF) to test the collinearity, which leads to some weaknesses such as unstable parameter estimation, unreliable models, and weak predictive ability. Given the result of the interaction test between age and transcriptional risk score, we defined the interaction term: age × risk score (A.R.). The VIFs of age, risk score, and A.R. were 3.918, 20.954, and 28.425, respectively, which indicated that the collinearity could exist between risk score and A.R. However, after mean-centering (each independent variable minus the corresponding average), the VIFs of age, risk score, and A.R. were 1.119, 1.231, and 1.106, respectively, pointing out the nonessential collinearity (27). Afterwards, two models were respectively constructed to assess whether adding A.R. can increase the performance of the model. As shown in Table 3 , there were no significant increase in terms of log-likelihood ratio (LLR), C-index, Akaike information criterion (AIC), Bayesian information criterion (BIC), and AUC of time-dependent ROC. Compared with model 2, the net reclassification improvement (NRI) and integrated discrimination improvement (IDI) of model 1 also did not improve ( Table 3 , all p > 0.05). Therefore, the age and risk score were included in the nomogram based on simplicity and efficiency.

The forest plot presented the age and risk score as independent risk factors ( Figure 8A ). The Schoenfeld residual test showed that all of the variables met equally proportional hazards (PH) assumption ( Figure 8B ) and there were no outliers based on the Deviance residual test ( Figure 8C ). The Martingale residuals demonstrated the linear relationship between age and transcriptional risk score with the logit transformation value of the hazard and the restricted cubic spline (RCS) analysis also verified ( Figures S11B–E ). Considering all the previously mentioned significant predictive factors, we established a comprehensive nomogram including age and transcriptional risk score ( Figure 8D ). We have also calculated the uncorrected and corrected C-index, which were 0.873 and 0.870, respectively ( Table 3 ). The calibration curves of 1 year, 3 years, and 5 years indicated a suitable calibration efficiency while a closer fitness to the dashed line indicates a better prediction performance ( Figures 8E–G ). The decision curve analysis (DCA) was used to assess the clinical applicability of nomogram and a net benefit for diverse prediction models at different threshold probabilities by adding the benefits and minimizing the harms. As demonstrated by the favorable probability, the comprehensive nomogram showed better net benefit than age and risk score ( Figures 8H–J ). Moreover, the time-dependent ROC curves verified that the prediction performance of the nomogram was gradually elevated along with the increase in time and also was better compared to the single index ( Figures 8K–M ). To create an intuitive application, the nomogram was further developed into a web version and could be dynamically operated online: https://drw576223193.shinyapps.io/Nomo/. Thus, the comprehensive nomogram was established according to the multiple prognostic factors that surpassed each single factor taken alone. The nomogram could help clinicians make more accurate assessment of patient prognosis.

To confirm the reliability of the comprehensive nomogram, the gene expression profiles were extracted from the CGGA dataset and then used for further model validation. Consistently, the expression level of the risk gene cluster was increased in IDH^wt/1p19q^non-codel samples and correlated with more severe prognosis of glioma patients ( Figures S12 , S13 ). Furthermore, the calibration curves of the comprehensive nomogram for the possibility of 1-, 3-, and 5-year OS displayed obvious concordance between the predicted results and observations ( Figures S14A–C ). In addition, the uncorrected and corrected C-index were 0.837 and 0.831, respectively, indicating that the comprehensive nomogram had an appropriate discrimination in the CGGA cohort. Similar to the results from the TCGA cohort, the ROC analysis demonstrated that our nomogram exhibited gratifying sensitivity and specificity on prognostic prediction with the AUCs of 0.845, 0.900, and 0.883 for the 1-, 3-, and 5-year survival, respectively ( Figure S14D ). Collectively, the results demonstrated that the comprehensive nomogram model validated both the training and the validation cohorts well.