Patients

The TCGA database was downloaded from UCSC Xena (https://xenabrowser.net/), which included 530 samples of LGG gene expression RNA sequencing (RNAseq) from the IlluminaHiSeq_RNASeqV2. 524 samples of LGG miRNA mature strand expression RNAseq from the IlluminaHiSeq_miRNASeq, relevant clinical information and 511 samples of LGG somatic mutation. The Chinese Glioma Genome Atlas (CGGA) database was downloaded from the CGGA website (http://www.cgga.org.cn/), which included 198 samples of GBMLGG (including 107 LGG and 91 GBM) gene expression miRNA from the V2.0 miRNA Expression BeadChip (Illumina). LGGs with missing overall survival (OS) values were excluded in order to reduce statistical bias during the analysis.

Selection of m⁶A genes and m⁶A-related miRNAs

We obtained the transcriptional profiles of miRNAs and m⁶A genes from the TCGA database. A total of 23 m⁶A regulators were included for study in this research, including writers (CBLL1, KIAA1429, METTL14, METTL3, RBM15, RBM15B, WTAP and ZC3H13), erasers (ALKBH5 and FTO), and readers (ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC, YTHDC1, YTHDC2, YTHDF1, YTHDF2 and YTHDF3) [27]. LASSO regression analysis was applied using the R package glmnet, and 13 m⁶A genes were obtained to be significantly correlated with the prognosis of LGGs. Then, we screened m⁶A-related miRNAs by Pearson’s correlation analysis, and 434 m⁶A-related miRNAs were identified. The criteria for screening by this procedure was |Pearson R| > 0.3.

Immunofluorescence and immunohistochemistry

Immunofluorescence and immunohistochemistry of m⁶A-related regulatory proteins in cerebral cortex and LGGs were obtained from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/).

Risk signature construction and validation

The TCGA database was randomly divided into a training dataset and a validation dataset. Baseline characterization showed no statistically significant differences in molecular and clinical characteristics between these two cohorts (P > 0.05) (Table 1). Next, using univariate Cox regression analysis, we screened prognosis-related miRNAs from the 434 m⁶A-related miRNAs in the TCGA database (P < 0.05). Using the R package glmnet for LASSO regression, we screened 22 m⁶A-related miRNAs that were significantly associated with OS in LGGs. Then, 22 m⁶A-related miRNAs were identified using multivariate Cox regression analysis to identify six m⁶A-related miRNAs with independent prognostic value. Finally, those six m⁶A-related miRNAs were used to construct a risk model. The risk scores were calculated as follows. Risk score = coef(miRNA1) × expr(miRNA1) + coef(miRNA2) × expr(miRNA2) + ...... + coef(miRNAn) × expr(miRNAn), where coef(miRNAn) is the coefficient of prognosis-related miRNAs and expr(miRNAn) is the expression of miRNAs. The classification of low- and high-risk groups was defined by the mean of the risk scores.

Bioinformatic analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) is used to carry out enrichment analysis. ESTIMATE was used to assess the percentage of stromal and immune cells in the tumor. The proportions of various immune cells in the tumor were evaluated using CIBERSORT. The R package maftools and the tumor immune dysfunction and exclusion (TIDE) algorithm were respectively used to analyze mutation data and immunotherapy response prediction. Single sample gene set enrichment analysis (GSVA; ssGSVA) was performed using GSVA to obtain a functional enrichment score for each sample [28].

Principal-component analysis and Kaplan-Meier survival analysis

High-dimensional data of 23 m⁶A genes and 6 m⁶A-associated miRNAs for risk modeling were efficiently downscaled, model identified and grouped for visualization using principal-component analysis (PCA). Kaplan-Meier survival analysis was used to assess the prognosis between the low- and high-risk groups.

Exploration of promising clinical therapeutic compounds targeting m⁶A-related miRNA risk model

To identify potential compounds for clinical treatment of LGGs, the half-maximal inhibitory concentration (IC₅₀) of listed compounds were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) website. The R package oncoPredict was used to predict the IC₅₀ of the above compounds in patients with LGG.

Independence of the m⁶A-related miRNA model

Univariate and multivariate Cox regression analyses were performed to identify independent variables for risk model and other clinicopathological characteristics of LGG patients. Simultaneously, we also comprehensively considered whether there is interaction and correlation between variables via R package corrgram and car.

Nomogram model construction and validation

The predictive performance of the nomogram and other variables [risk score, age, WHO grade, Karnofsky performance status (KPS), and histological type] for the 1-, 3-, and 5-year OS was established. The receiver operating characteristic (ROC) curve is used to evaluate the discrimination ability of the model. A network-based interactive dynamic nomogram application was constructed using Shiny version 1.6.0.

Selection of m⁶A-related miRNAs in LGGs

We screened prognostic m⁶A genes from 23 m⁶A genes in the TCGA database using LASSO regression analysis (Figure 1A, 1B). We defined m⁶A-related miRNAs as miRNAs associated with one or more of the 13 m⁶A genes (|Pearson R| > 0.3). Sankey plots were used to visualize the m⁶A-miRNA co-expression network (Figure 1C), and 434 m⁶A-related miRNAs were discerned as m⁶A-related miRNAs. Finally, immunofluorescence results revealed that m⁶A-related regulatory proteins were mainly expressed in the nucleus and cytoplasm of glioma cells. Immunohistochemical results demonstrated that WTAP, RBM15, IGF2BP2, HNRNPA2B1, and ELAVL1 were more expressed in LGG than in cerebral cortex, while YTHDC1 and LRPPRC were more expressed in cerebral cortex than in LGG (Figure 2, Figure S1).

Display full size

Figure 1. 

Identification of prognostic m⁶A genes in LGG patients. (A, B) LASSO regression analysis using 10-fold cross-validation; (C) Sankey plots of m⁶A genes and m⁶A-related miRNAs

Display full size

Figure 2. 

Immunofluorescence and immunohistochemistry of m⁶A-related regulatory proteins in cerebral cortex and low-grade gliomas

Risk signature construction and validation

We screened m⁶A-related prognostic miRNAs from 434 m⁶A-associated miRNAs in the TCGA database using univariate Cox regression analysis. Those m⁶A-related miRNAs that were significantly associated with prognosis were screened in the TCGA database (P < 0.05). We further obtained 22 m⁶A-related miRNAs that were significantly associated with OS through the LASSO regression analysis (Figure 3A). Next, in the TCGA database, six m⁶A-related miRNAs with independent prognostic value were selected based on clinical information and used to construct a risk model (Figure 3B). The correlation analysis revealed that the six m⁶A-related miRNAs were significantly correlated with m⁶A regulators (Figure S2A). Then, we found that MIMAT0004951 (MIR887) has m⁶A modification sites (Figure S2B) based on these six miRNA transcripts through online website sequence-based RNA adenosine methylation site predictor (SRAMP; a mammalian m⁶A sites predictor).

To assess the prognostic ability of the risk model, the LGG sample was divided into low- and high-risk groups based on the mean of the risk scores in the training and validation cohorts. The heatmap depicts the expression of six m⁶A-related miRNAs involved in the construction of the risk model in the training cohort (Figure 3C) and the validation cohort (Figure 3D). Kaplan-Meier survival analysis of both the training and validation cohorts showed that patients with LGG in the high-risk group had worse OS than those in the low-risk group (Figure 3E, 3F). To further assess the prognostic predictive ability of the risk model, we included disease-specific survival (DSS) and progression-free interval (PFI) metrics to assess the differences between low- and high-risk LGG patients. The results demonstrated that LGGs in the high-risk group had worse DSS and PFI than low-risk group. The above results suggest that the risk model of m⁶A-related miRNAs can effectively assess the prognosis of LGGs (Figure 3E, 3F). A total of 523 LGG patients, who were included in the study, were randomly divided into a training cohort (n = 393) and a validation cohort (n = 130). Detailed characteristics for these two cohorts showed homogeneity in these cohorts (Table 1).

Display full size

Figure 3. 

Risk model for LGG patients based on m⁶A-related miRNAs. (A) LASSO regression analysis using 10-fold cross-validation; (B) univariate and multivariate Cox regression analysis showed six independent prognostic miRNAs; (C, D) heat plots displayed the expression of six m⁶A-related miRNAs in the training and validation cohorts; (E, F) Kaplan-Meier analysis of the prognosis in the low-risk and high-risk groups. HR: hazard ratio; CI: confidence interval; ^* statistical significance

PCA of m⁶A related miRNA risk signature model

PCA was performed to assess the differences between the low- and high-risk groups based on the transcriptional profiles of 23 m⁶A genes and 6 m⁶A-related miRNAs involved in the construction of the risk model (Figure S2C, 2D). The results of PCA indicated that 6 m⁶A-related miRNAs had better distinction than 23 m⁶A genes in terms of the distribution of the low- and high-risk groups. Next, the ssGSEA results showed compared with low-risk group, LGGs in high-risk group had a higher m⁶A writer and reader scores, but a lower eraser score (Figure S2E).

We analyzed the prognostic differences between the low- and high-risk groups in the TCGA and CGGA database stratified by clinicopathological characteristics. First of all, compared with age ≤ 40, KPS ≥ 70, primary tumor, WHO grade 2, oligodendroglioma, and IDH mutant groups, the risk score was higher in age > 40, KPS < 70, recurrent tumor, WHO grade 3, astrocytoma, and IDH wildtype groups (Figure 4A–G; Figure S3A–G). Then, according to the subgroups stratified by age, gender, KPS, sample type, WHO grade, histological type, and IDH status, the OS of the high-risk group continued to be inferior to that of the low-risk group (Figure 4H–N; Figure S3H–N).

Display full size

Figure 4. 

Kaplan-Meier analysis of survival differences stratified by clinicopathological variables between the low- and high-risk groups in the TCGA database. (A–G) Risk score differences stratified by clinicopathological variables between the low- and high-risk groups; (H–N) Kaplan-Meier analysis of survival differences stratified by clinicopathological variables between the low- and high-risk groups

Correlation between m⁶A related miRNA risk signature and immune microenvironment of LGGs

We found that the risk model of m⁶A-related miRNAs was mainly correlated with immune-associated biological processes by performing gene ontology (GO) enrichment analysis of differently expressed genes (DEGs) between the low- and high-risk groups (Figure 5A, 5B). Moreover, the ESTIMATE results suggested that compared with low-risk group, LGGs in high-risk group had higher immune and stromal scores, while the purity of tumors was lowered (Figure 5C, 5D).

Display full size

Figure 5. 

Enrichment analysis and ESTIMATE scores in the low- and high-risk groups of LGGs. (A, B) Enrichment analysis of DEGs in the low- and high-risk groups of LGGs in the TCGA database. The top 20 items were illustrated in the bubblechart; (C) the heat map displayed the distribution features of ESTIMATE scores between the low- and high-risk groups of LGGs in the TCGA database; (D) the scatter plot depicted the results of the ESTIMATE scores between the low- and high-risk groups

We next explored the relationship between risk models of m⁶A-related miRNAs and tumor-infiltrating immune cells (TICs). Immune infiltration analysis revealed higher proportion of M1 macrophage, regulatory T (Treg) cell, M0 macrophage, and M2 macrophage in the high-risk group was than that in the low-risk group (Figure 6A, 6B). We also founded that the expression of exhausted immune response markers was markedly increased in the high-risk group compared to the low-risk group (Figure 6C). Since M1 macrophages exert pro-inflammatory and anti-tumor functions, while M2 macrophages cells have anti-inflammatory and pro-tumor functions, we further analyzed the difference in secretion of immune factors between M1 and M2 macrophages in the low- and high-risk groups. The results showed that the anti-inflammatory cytokines transforming growth factor beta (TGFB), interleukin 10 (IL10) and indoleamine 2,3-dioxygenase (IDO) secreted by M2 macrophages were significantly higher in the high-risk group than in the low-risk group. The pro-inflammatory cytokines tumor necrosis factor (TNF), IL23, and IL12 secreted by M1 macrophages in the high-risk group were not different from those in the low-risk group or were lower than those in the low-risk group (Figure 6D). The above results indicate that M1 macrophages in the high-risk group occurred disability.

Display full size

Figure 6. 

Proportion of TICs and exhausted immune response markers in the low- and high-risk groups of LGGs. (A) Bar chart showed the percentage of various TICs in the samples from the low-risk and high-risk groups; (B, C) bar chart demonstrated the difference in the expression of various immune cells and exhausted immune response markers between the low-risk and high-risk groups; (D) expression of M1- and M2-associated genes differences between the low- and high-risk groups in the TCGA database; (E) TIDE prediction difference in the high- and low-risk patients, and accumulative bar diagram showed that the LGGs with low-risk group (49.67%, 150/302) were more likely to respond to immune checkpoint blockade (ICB) immunotherapy than the LGGs with high-risk group (37.1%, 82/221). ^* P < 0.05, ^** P < 0.01, ^*** P < 0.001, ^**** P < 0.0001

In addition, the results of TIDE discovered that LGGs in the low-risk group had a higher response rate to immunotherapy than in the high-risk group. The above results suggest that m⁶A-related miRNA-based risk scores may be a potential predictor of immunotherapy response (Figure 6E).

Then, we analyzed and summarized the mutation data. The waterfall plot shows the top 20 driver genes with the highest frequency of non-silent mutations between the low- and high-risk groups (Figure 7A, 7B). At the same time, the mutation frequency of IDH1 in the low-risk group (92%) was obviously higher than that in the high-risk group (56%). Therefore, we evaluated whether a risk model based on m⁶A-related miRNAs could better predict OS than IDH mutation status. The results showed that LGGs with IDH mutations and IDH wildtype in the high-risk group (defined as IDH_Mut high risk and IDH_wt high risk, respectively) had worse OS than those with IDH mutations and IDH wildtype in the low-risk group (defined as IDH_Mut low risk and IDH_wt low risk, respectively) (Figure 7C, 7D). Interestingly, when OS was less than 2,660 days, LGGs with IDH mutant in the high-risk group (IDH_Mut high risk) had better survival outcomes than those with IDH wildtype in the low-risk group (IDH_wt low risk). However, when OS was greater than 2,869 days (in the TCGA database) or 1,314 days (in the CGGA database), LGGs with IDH mutant in the high-risk group (IDH_Mut high risk) had worse survival outcomes than those with IDH wildtype in the low-risk group (IDH_wt low risk) (Figure 7C, 7D). Therefore, the above results suggest that risk models for m⁶A-related miRNAs are likely to have greater prognostic reference value than IDH mutation status.

Display full size

Figure 7. 

Genetic alteration analyses between the low- and high-risk groups of LGGs. Waterfall plot indicated gene mutations in the (A) low-risk and (B) high-risk groups; (C, D) prognostic evaluation of LGGs survival by IDH mutation status and m⁶A-related miRNA risk score in TCGA and CGGA databases

Exploration of promising clinical therapeutic compounds targeting m⁶A-related miRNA risk model

To identify potential drugs for clinical treatment of LGGs, we used the GDSC database to calculate the IC₅₀ for each sample using the oncoPredict algorithm. A total of 104 compounds were screened for significant differences in IC₅₀ between the low- and high-risk groups. The results of comparing the IC₅₀ of the compounds between the two groups revealed that more compounds were more sensitive in the low- risk group compared to the high-risk group (Figure S4A). The difference between the top 5 compounds in the low- and high-risk groups of LGGs (Figure S4B–F).

Prognostic and clinical characterization of m⁶A-related miRNA risk model

Univariate and multivariate Cox regression analyses revealed risk score, age, WHO grade, and histological type were independent prognostic variables for LGGs (Figure S5A). Simultaneously, the results of the ROC curve analysis indicated that the area under the curve (AUC) of the risk score was always greater than other clinicopathological variables at 1, 3 and 5 years of OS. The above results suggest that the risk model based on six m⁶A-related miRNAs is more reliable for the prognostic evaluation of LGGs (Figure S5B). Moreover, we comprehensively considered whether there is interaction and correlation between variables (Table S1, S2), and found that the risk score and IDH mutation status have a significant correlation (Pearson r = –0.659).

Construction and evaluation of the prognostic dynamic nomogram

The nomogram was used to predict the incidence of OS at 1, 3, and 5 years, the risk scores of the nomogram model showed dominant predictive power by comparison with other clinicopathological characteristics (Figure 8A). The ROC curve for the nomogram and its AUC displayed in Figure 8B. The final Cox model incorporated five independent predictors (risk score, age, WHO grade, KPS, and histological type) and was developed as a simple-to-use dynamic nomogram that is available online (https://fengyuan1993.shinyapps.io/LGG_m6A_Immunotherapy/) as screenshotted in Figure 8C, 8D.

Display full size

Figure 8. 

Construction and evaluation of a prognostic nomogram. (A) The nomogram predicts the probability of the 1-, 2-, and 3-year OS; (B) ROC curve of the nomogram predicts the probability of the 1-, 2-, and 3-year OS; (C, D) online dynamic nomogram accessible at https://fengyuan1993.shinyapps.io/LGG_m6A_Immunotherapy/