Gene expression profiles and survival data for patients with LGG were downloaded from the CGGA (http://www.cgga.org.cn/) database. Overall, 495 samples from TCGA and 172 samples from the CGGA were analyzed (Table 1). Perl and R in R software were utilized to evaluate the samples.

The expression level of a total of 26 m6A RNA methylation regulatory genes in LGG samples and normal samples in TCGA and CGTA datasets were detected by the R package “limma” and visualized by the R package “pheatmap.” The correlations between 26 m6A RNA methylation regulatory genes were detected by the R package “corrplot.”

A total of 172 patients from the CGGA database were randomly assigned in a 1:1 ratio to a training set and a validation set. The hazard ratio of OS of 26 m6A regulatory genes in the internal training set was calculated using the univariate Cox proportional hazard regression model, and genes with p < 0.05 were considered statistically significant and included in subsequent analyses. Then, m6A-related genes that were identified as significant in the univariate analysis in the internal training set were included in the penalized Cox regression model with the least absolute shrinkage and selection operator (LASSO) Cox regression model for 10-fold cross validations to select the most significant genes. Finally, an m6A-related prognostic gene signature was constructed based on a linear combination of the regression coefficient derived from the LASSO Cox regression model coefficients multiplied by the mRNA expression level (Tibshirani,1997).

The risk score was calculated as follows: ∑ C o e f i   ×   E x p i , where Coefi is the coefficient of prognostic biomarker and Expi is the expression of the corresponding prognostic biomarker. According to the median value of the risk score, LGG patients in the CGGA cohort were divided into low-risk and high-risk groups. The overall Kaplan–Meier survival curves of low-risk and high-risk patients were analyzed by the R package “survival” and “survminer.”

The receiver operating characteristic curve (ROC) analysis was used to assess the accuracy of the diagnostic gene signature. The R software package “pROC” was used to generate an ROC curve (Robin et al., 2011). ROC area under curve (AUC) was calculated to evaluate the performance of each signature, and AUC>0.7 was considered as a diagnostic gene signature.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis led to identification of enrichment pathways, which help to determine the biological pathways to which the identified genes belonged (Yu et al., 2012). The KEGG analysis is performed by the R package “clusterProfiler.”

Gene Set Enrichment Analysis v2.2.2 (http://www.broadinstitute.org/gsea) was used to investigate the biological differences among patients with different expression patterns of the 5-gene signature. In addition, C2:CP KEGG gene sets from the Molecular Signatures database (MSigDB) were used as the reference gene sets.

The gene signatures were used to construct a nomogram by the “rms” R package, and the accuracy of the nomogram was evaluated by the calibration curve (1,000 bootstrap resamples).

Continuous variables were analyzed using Student’s t-tests or non-parametric tests. Categorical variables were analyzed using chi-squared tests or Fisher’s exact tests. The R package copy was used for univariate and multivariate analyses. All data were analyzed using SPSS or R software (http://www.r-project.org/). p < 0.05 was considered statistically significant.

A total of 26 m6A RNA methylation regulator genes were included in this study to evaluate the methylation status of tumor tissues. We constructed a profile containing 26 m6A regulatory genes in patients with LGG. The m6A regulatory genes and corresponding clinical parameters are shown in Figure 2A. As shown by the heatmap, the expression levels of the m6A regulatory genes differed among the patients. To better understand the relationships among m6A regulatory genes in LGG, we performed a regression analysis of the m6A regulatory genes (Figure 2B). The results showed that YTHDF2 was highly correlated with RBM15. In addition, there was a strong correlation between VIRMA and YTHDF3. These results indicated that these pairs of genes may be involved in the same biological functions. To determine which m6A regulatory genes had prognostic value, we conducted a univariate COX analysis. This analysis identified 16 of the 26 methylated genes as prognostic indicators (Figure 2C).

A total of 172 patients in the CGGA cohort were randomly divided in a 1:1 ratio into an internal training cohort and a validation cohort. The 16 m6A regulatory genes with prognostic value were included in the LASSO Cox model in the internal training CGGA cohort to screen for characteristic variables (Figures 3A,B). As shown in Figure 2A, the minimum lambda value was reached with inclusion of 5 genes. Therefore, we selected 5 genes to construct a prognostic model. Finally, we constructed the following formula for prognosis of OS in patients with LGG: (formula= (0.033 ×HNRNPC) + (0.237 × IGF2BP2) + (0.260 × IGF2BP3) − (0.271 × LRPPRC) + (0.713 × YTHDF2)). According to the cutoff value (2.03) obtained using the survminer package, the patients in the training cohort were divided into high- and low-PRI groups. Patients with high risk had more events and a worse prognosis (Figure 3C). Application of the cutoff value to the internal validation CCGA cohorts and external validation TCGA cohorts produced the same result as that observed in the training group (Figures 3D,E). To evaluate the accuracy of this 5-gene signature for determination of prognosis, we generated an ROC curve. In the internal training cohort, the accuracy of the 5-gene signature was investigated at 2, 3, and 5 years, which resulted in AUC values of 0.917, 0.889, and 0.874, respectively (Figure 3F). In the internal validation cohort, the AUC values were 0.787, 0.768, and 0.834, respectively (Figure 3G). In the external validation TCGA cohort, the AUC values at 2, 3, and 5 years were 0.783, 0.721, and 0.701, respectively (Figure 3H). In summary, the 5-gene profile was an excellent prognostic indicator for patients with LGG.

Based on the LASSO Cox model, the prognostic risk score for each patient in the CGGA cohort was calculated according to the factor value and expression level. According to the median risk score, patients were divided into a high-risk group and a low-risk group (Figures 4A–C). The high-risk group had higher levels of IGF2BP2, IGF2BP3, HNRNPC, and YTHDC2 and lower levels of LRPPRC. More patients survived in the low-risk group than in the high-risk group. Further correlation analysis found a negative correlation between LRPPRC and IGF2BP2 and between LRPPRC and IGF2BP2 (Figure 4D). We further evaluated the stability of 5-gene signature in different groups. The results showed that the 5-gene signature was an excellent prognostic indicator regardless of gender (Figures 4E,F), grade (Figures 4G,H), IDH mutation status (Figures 4I,J), and the co-mutation state of chromosome 1p and 19q (Figures 4K,L).

To elucidate the differences in biological characteristics between the high- and low-risk groups using the 5-gene signature, we performed Spearman’s correlation analysis and selected the first 1,000 genes. Then, the clusterProfiler package in R software was used to perform KEGG enrichment analysis. These genes are significantly enriched in the focal adhesion and cell cycle pathways (Figure 5A). Specifically enriched genes for each KEGG term are summarized in Figures 5C,D. Next, we divided the patients into high- and low-risk groups according to their risk scores for the GSEA analysis, which showed enrichment of the 5-gene signature in cell cycle pathways (Figure 5B), which indicated that the cell cycle pathway may be a critical factor in poor prognosis in patients with LGG. In addition, expression levels of genes in the T-cell receptor signaling pathway were abnormally high in patients with LGG, which indicated that immune status may have differed with level of risk. Therefore, we further explored the expression levels of immune checkpoint (IC) genes in the different patient risk groups, which may provide additional information regarding personalization of treatment.

Previous studies have shown that the expression levels of immune checkpoint (IC) genes are associated with immunotherapy efficacy. We compared the expression patterns of IC genes (PD-1, PD-L1, and CTLA-4) in different risk groups of patients with LGG. As shown in Figure 6A, patients in the high-risk group had higher expression levels of IC genes. We further evaluated whether our research model for IC genes was generalizable to patients with similar expression levels of IC genes. As shown in Figure 6B, the overall survival of patients with low risk scores and high PD-1 expression was significantly better than that of patients with high risk scores and high PD-1 expression. Furthermore, the survival of patients with low risk scores and low PD-1 expression was longer than that of patients with high risk scores and low PD-1 expression. Similar results were observed for PD-L1 and CTLA-4 (Figures 6C,D). Stratification of patients with low risk scores according to IC gene expression showed that IC gene expression was significantly associated with survival of patients with low risk scores. However, there was no differences in survival of patients with high risk scores when stratified based on IC gene expression. In addition, patients with low risk scores and low IC gene expression tended to have much higher survival rates than patients in the other three groups. These results suggested that our 5-gene model may be a predictor of ICI immunotherapy efficacy.

To provide patients with more accurate prognostic predictions, we incorporated the 5-gene signature and clinical parameters that had prognostic value in the univariate analysis and performed a multivariate analysis in the CGGA cohort. The results showed that the 5-gene signature, WHO grade, gender, and X1p19q codeletion status were stable predictors (Supplement Tables S1,2). These results were used to construct a nomogram to predict the prognosis of patients (Figure 7A). The calibration curve shows that the nomogram had stable predictive values at 3 and 5 years (Figures 7B,C).