Malignant Evaluation and Clinical Prognostic Values of m6A RNA Methylation Regulators in Glioblastoma

N6-methyladenosine (m6A) RNA methylation, the most common form of mRNA modification and regulated by the m6A RNA methylation regulators (“writers,” “erasers,” and “readers”), has been reported to be associated with the progression of the malignant tumor. However, its role in glioblastoma (GBM) has been poorly known. This study aimed to identify the expression, potential functions, and prognostic values of m6A RNA methylation regulators in GBM. Here, we revealed that the 13 central m6A RNA methylation regulators were firmly related to the clinical and molecular phenotype of GBM. Taking advantage of consensus cluster analysis, we obtained two categories of GBM samples and found malignancy-related processes of m6A methylation regulators and compounds that specifically targeted the malignant processes. Besides, we also obtained a list of genes with poor prognosis in GBM. Finally, we derived a risk-gene signature with three selected m6A RNA methylation regulators, which allowed us to extend the in-depth study and dichotomized the OS of patients with GBM into high- and low-risk subgroups. Notably, this risk-gene signature could be used as independent prognostic markers and accurate clinicopathological parameter predictors. In conclusion, m6A RNA methylation regulators are a type of vital participant in the malignant progression of GBM, with a critical potential in the prognostic stratification and treatment strategies of GBM.

The vital functions of RNA modification in processes of life have caught people's eyes in recent years. Substantial progress in regulating RNA transcription (18,19), the event of processing (20,21), splicing (5,22), RNA stabilities (23,24), and translation (25,26) was witnessed in m6A posttranscriptional modifications. However, to date, the functions of the majority of RNA modifications found in mRNAs need further exploration. Notably, the functional roles of m6A methylation in tumorigenesis, tumor differentiation (27), proliferation (28), and invasion (27) remain elusive.
GBM is the most common and devastating primary tumor in the brain. Even the combined surgical resection, radiation therapy, chemotherapy, and other therapies were broadly used, the recurrence of the patients with GBM is inevitable. Besides, the median survival of GBM patients is <15 months after a definite diagnosis (29)(30)(31). The m6A RNA methylation regulators were also reported to be associated with self-renewal, radio resistance, and tumorigenesis of GBM stem cells (32). However, there is no comprehensive investigation of the expression of m6A RNA methylation regulators in GBM.
In the current study, the 13 m6A RNA regulators, which have been widely reported, were systematically analyzed using the GBM RNA sequencing data from The Cancer Genome Atlas (TCGA) (n = 174) and Chinese Glioma Genome Atlas (CGGA) (n = 249) databases. Taking advantage of m6A RNA methylation regulator-based consensus clustering analysis, we demonstrated the malignant process and obtained a list of genes with poor prognosis in patients with GBM. Importantly, we further validated these genes in the CGGA database and identified potential drugs targeting the malignant process of GBM using the Connectivity Map (CMap) (33). Besides, the risk-gene signaturederived from m6A RNA methylation regulators might be used as a novel biomarker that could identify GBM patients' prognosis and predict the clinicopathological parameters of GBM.

Data Acquisition
The RNA-seq transcriptome data and corresponding clinicopathological parameters of GBM patients were obtained from the TCGA database (http://cancergenome.nih.gov/) and the CGGA database (http://www.cgga.org.cn). The RNA-seq transcriptome data of healthy human tissue was obtained from the Genotype-Tissue Expression (GTEx) database (http://commonfund.nih.gov/GTEx/). We combined GTEx and CGGA data, and then harmonized them using quantile normalization and svaseq-based batch effect removal (34). The clinicopathological parameters for the CGGA and TCGA datasets are summarized in Table S1.

Selection of m6A RNA Methylation Regulators
Thirteen widely recognized m6A RNA methylation regulators were retrieved from published literature. We then systematically compared the correlation between the expression of these m6A RNA methylation regulators and clinicopathological parameters in GBM patients.

Bioinformatic Analysis
To further explore the role of m6A RNA methylation regulators in GBM patients, we clustered the GBM patients into two clusters by using the R package ConsensusClusterPlus (35). Heatmaps were drawn based on the average linkage method and the Pearson distance measurement method. Principal Component Analysis (PCA) was carried out by an R package called PCA to observe the distribution of gene expression in two clusters. Differential analyses for each gene in the pre-classified samples performed using the limma package in R (36). Fold change (FC) > 2 and adjusted p-value (q-value) < 0.01 were set as the cutoff values to screen for differentially expressed genes (DEGs). Gene Ontology (GO) functional analyses and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to analyze the upregulated DEGs. The relationship between DEGs gene expression levels and overall patient survival time was illustrated by generating Kaplan-Meier plots. The correlation was tested using a log-rank test. Gene Set Enrichment Analysis (GSEA) was used to investigate the functions correlated with different clusters of GBM. |NES| > 1, adjusted p < 0.05, and FDR q < 0.25 were considered as statistically significant as described in a previous study (37).

Construction of Protein-Protein Interactions (PPI) Network
PPI among selected genes was analyzed using the STRING database (38) and reprocessing via Cytoscape software (39). For better visualization, the color of the node in the PPI network was applied to reflect the logFC value, and the size of the node was applied to indicate the number of source proteins with the target protein. Molecular COmplex Detection (MCODE) (version 1.4.3), which is clustered based on given network topology, was used to discover densely connected regions. Then, the most significant module was filtered out by MCODE from the PPI networks. The criteria for selection were as follows: MCODE scores > 5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2.

Construction of Gene-Signature
Univariate Cox regression analysis of the expression of 13 m6A RNA methylation regulators was conducted to determine the candidate genes associated with overall survival (OS). After that, an L1-penalized (LASSO) was performed to further identify the selected genes with independent prognostic value (40,41).
Finally, their regression coefficients were determined by the minimum criteria. The risk score for the signature was calculated accurately by the formula: where Coefi is the regression coefficient and xi is the expression of each selected gene. GBM patients were divided into lowand high-risk subgroups according to the median risk score. Kaplan-Meier plot was performed to compare the OS between two risk subgroups.

Identification of Potential Compounds Targeting the Malignancy-Related Pathways
CMap (updated in September 2017) (https://clue.io/), as the world's largest perturbation-driven gene expression dataset, was employed to search for candidate chemical compounds that might target GBM stemness related pathways (33). The compounds were discovered by interrogating the CMap database of signatures with a query (a list of DEGs relevant to biological features of interest). The final results involved a CMap connectivity score (from −1 to 1) that indicated the degree of specificity associated with our particular query. 300 DEGs (150 downregulated and 150 upregulated) were selected for query methodology. Noteworthily, the closer the connectivity score of a compound was to −1, the more likely it was to reverse the genetic pattern we are querying. Finally, the compounds with the absolute value of CMap connectivity score of 90 or higher were considered to be potential therapeutic agents for functional validation.

Statistical Analysis
Chi-square tests were used to compare the expression levels in GBM for age, gender, healthy samples, primary GBM and recurrence GBM, isocitrate dehydrogenase (IDH) status, and cytosine-phosphate-guanine island methylator phenotype (G-CIMP) status. One-way ANOVA was used to compare the distribution of the subtype of GBM (Classical, Mesenchymal, Neural, Proneural) (42). To evaluate the prediction accuracy of the risk score model, we performed a receiver operating characteristic (ROC) curve and calculated the area under the curve (AUC). Potential prognostic factors such as age (≤ 65 vs. > 65), gender (female vs. male), GBM subtype, and risk score (lowrisk vs. high-risk) were analyzed by Univariate and multivariate Cox hazard regression.

Expression Patterns of m6A RNA Methylation Regulators in GBM
According to the essential biological functions of methylation regulators in the development of GBM, we first analyzed the relationship between each m6A RNA methylation regulator and the clinical molecular phenotype of GBM. The expression level of individual m6A RNA methylation regulator and different types of samples was presented as a heatmap. The result strongly indicates that the expression of the majority of m6A RNA methylation regulators was associated with the occurrence of GBM ( Figure 1A). Importantly, the significant correlation between the occurrence of GBM and the expression levels of ALKBH5, METTL3, KIAA1429, HNRNPC, WATP, YTHDC2, YTHDF1, YTHDF2, and FTO were confirmed by the quantitative analysis of CGGA ( Figure 1B). Compared with the healthy samples, the expression of METTL3, HNRNPC, WTAP, KIAA1429, YTHDF2, and YTHDF1 was upregulated, while the expression of ALKBH5, YTHDC2, and FTO was downregulated in the GBM samples (Figures 1C,D). Correlation analysis was also employed to investigate the relationship between the expression level of m6A RNA methylation regulators and the different stages (primary tumor stage and recurrent tumor stage) of GBM. Among the 13 m6A RNA methylation regulators, only HNRNPC was significantly related to the cancer recurrence ( Figure 1E). Considering the dramatically imbalanced numbers of primary GBM (n = 156) and recurrent GBM (n = 13) in the TCGA database, the results of TCGA database analysis were not necessarily as accurate as those of CGGA, while the numbers of primary GBM (n = 140) and recurrent GBM (n = 109) were relatively balanced. Therefore, we analyzed the expression profile of the primary and recurrent GBM samples from the CGGA database and found that HNRNPC exhibited no correlation. Interestingly, as shown in Figure 1F, the expression of WTAP, ALKBH5, and METTL14 was significantly associated with the recurrence of GBM. These results suggested that WTAP, ALKBH5, and METTL14 were firmly related to the recurrence process of GBM (Figures 1G,H). We further explored the relationship between the expression of m6A RNA methylation regulators and GBM molecular subtypes. Notably, every expression of the majority of m6A RNA methylation regulators was associated with the subtype of GBM, except YTDHF2 and HNRNPC ( Figure S1A). We also investigated the relationship between IDH status, G-CIMP status, and expression levels of each m6A RNA methylation regulator in GBM. The results revealed that the expression levels of METTL14, KIAA1429, YTHDC1, ZC3H13, and FTO were significantly dysregulated between different IDH status in the TCGA dataset ( Figure S1B). As for the G-CIMP status, the expression of RBM15, YTHDF2, KIAA1429, and YTHDC1 exhibited a significant difference between G-CIMP+ and G-CIMP-( Figure S1C). We speculate that the change in the correlation of m6A RNA methylation regulators may be an internal characteristic that can reflect the external differences. From Figures 1I,J, different degrees of the relationship were observed between different m6A RNA methylation regulators. Most of the relationships between the regulators were positive correlations, especially YTHDC1, which contained the most active correlation with other regulators (Figures 1I,J).

Identification of Two Clusters of GBM Samples With Different Clinical Characteristics
Next, GBM samples with complete clinical parameters were selected for the subsequent consensus clustering analysis. From the view of the number of samples per group, an unbalanced  distribution was observed in the three groups when k = 3 ( Figure S2A). Hence, based on the expression similarity of the 13 regulators, k = 2 was the most optimum with clustering stability datasets increasing from k = 2-10 (Figures 2A-C and Figure S2). Then, GBM samples from the TCGA dataset were pre-classified into two groups (52 samples in one group labeled as RM1 and 106 samples in another group labeled as RM2) through consensus cluster analysis. The clinical features of the two groups are summarized in Table S2. The heatmap of cluster analysis showed that the 13 regulators could distinguish different samples, and the samples in the same cluster possessed a high correlation ( Figure 2C). Principal component analysis was performed to elucidate the difference in transcriptional profiles between the RM1 and RM2 subgroups. The results investigated that there was a clear distinction between these two subgroups ( Figure 2D). The survival curve according to Kaplan-Meier survival analysis for the clustered samples revealed a noticeable decrease in the OS in the RM2 subgroup compared with the RM1 subgroup, suggesting that the 13 methylation regulators could classify the GBM samples in prognostic level ( Figure 2E). We further found that the median survival of the RM1 group was 1.4 years, while the RM2 was only 1 year. In addition, the clinicopathological features of these two subgroups were compared. The RM1 subgroup was markedly correlated with younger age at diagnosis (P < 0.05), Neural or proneural subtypes (P < 0.001), and G-CIMP-status ( Figure 2F). The RM2 subgroup mainly contained GBM with older age at diagnosis, classical or mesenchymal subtypes, and G-CIMP+ status. Consistent with the report that classical and mesenchymal were more malignant compared to neural and proneural (42).

Functional Annotation of Classification Determined by Consensus Clustering Analysis
The above results indicate that the consensus clustering results were closely related to the degree of malignancy of GBM.
To better understand the mechanisms between the malignancy of GBM and the 13 m6A RNA methylation regulators, a total of 2,299 genes (599 genes were upregulated, and 1,700 genes were downregulated) were identified as DEGs by using differential analysis ( Figure 3A). To summarize the potential function of DEGs, we performed an annotation of the 599 significantly upregulated genes (onco role, Table S3) in the RM2 subgroup through GO function analysis and KEGG pathway analysis. The top 10 GO terms indicated that upregulated genes were enriched in malignancy-related processes, including neutrophil-mediated immunity, cell proliferation, cell junction, phagocytosis, and cell-substrate adhesion ( Figure 3B). KEGG pathway analysis top 10 terms testified that upregulated genes were related to the regulation of actin cytoskeleton pathway, focal adhesion pathway, proteoglycans in cancer pathway, and Fc gamma R-mediated phagocytosis pathway ( Figure 3C). Furthermore, GSEA suggested that the malignant hallmarks of tumors, including KRAS signaling (NES = 1.59, normalized P = 0.013), inflammatory response (NES = 1.68, normalized P = 0.052), myogenesis (NES = 2.02, normalized P < 0.001), and IL-6/JAK/ STAT3 signaling (NES = 2.0, normalized P < 0.001), were significantly associated with the RM2 subgroup (Figures 3D-G). All these results proved that the two categories derived from consensus clustering analysis were closely related to the malignancy of GBM.

Novel Candidate Compounds Targeting the Malignancy-Related Pathways and Biological Functions in GBM
Next, we sought to determine the potential compounds that target malignancy-related pathways and biological functions; DEGs based on consensus clustering were submitted to retrieve the CMap databases (33). The top 89 compounds capable of repressing the above gene expression of GBM were summarized in Table S4 and

Identification and Analysis of m6A-Related Genes With Poor Prognosis in GBM
To explore the significance of each upregulated gene for the survival time of GBM patients from the TCGA database, Kaplan-Meier survival curves were generated. Among the 599 upregulated DEGs in the RM2 subgroup, a total of 79 DEGs (Table S3) were able to predict the poor OS in the log-rank test (P < 0.05, representative figures were shown in Figure 5).
To better understand the interactions between the 79 genes with prognostic value and the 13 m6A RNA methylation regulators, we also analyzed the PPI among them using the STRING database. The network consists of four modules with a total of 88 nodes and 527 edges, which indicates that close interaction exists in this PPI network (Figure 6A). Among the 88 nodes, 54 central node genes (bold in Table S5) were selected with the filtering of degree > 10. The most significant 10 node degree genes were ITGAM, STAT3, SPI1, TNFRSF1B, MYO1F, SLC11A1, TCIRG1, RAP2B, FERMT3, and LCP1. The top two significant modules were selected by using the MCODE application for further analysis. To make the description more convenient, we named these modules "ITGAM module" and "RAP2B module, " respectively. Twenty-eight nodes and 165  edges were involved in ITGAM modules, with STAT3, SPI1, TNFRSF1B, SLC11A1, and FERMT3 being the remarkable nodes, as they had the most connections with other nodes in this module ( Figure 6B). In the RAP2B module, 8 edges involving 6 nodes were formed in this network ( Figure 6C). In addition, we predicted the function of the ITGAM module through GO analysis. It was related to the biological process of mRNA splicing via spliceosome, mRNA methylation, and oxidative single-stranded RNA demethylation. For instance, ITGAM was reported to play a critical role in invasive growth and angiogenesis in malignant gliomas (43). STAT3 methylation via STAT3 signaling could also promote tumorigenicity of GBM stem-like cells (44). The above results clearly demonstrate that m6A regulators participate in the critical malignant related biological regulatory network.
To find out whether the 79 OS-related DEGs found in the TCGA database were meaningful in the additional database, we further analyzed the expression profiles of 249 GBM cases from the CGGA database. Importantly, a total of 64 genes were validated to be significantly related to the poor prognosis, of  Table S3) were of particular interest, as they have not been previously reported for their prognostic value in GBM patients ( Figure S3).

Prognostic Value of m6A RNA Methylation Regulators
For the purpose of investigating the prognostic value of m6A RNA methylation regulators, univariate Cox regression analysis was performed on the expression profile data. Based on the information contained in these results, 4 of 13 genes exhibited a significant correlation with the prognosis. Among these four selected genes, HNRNPC, ALKBH5, and ZC3H13 were risky genes, with HR > 1, while FTO was a protective gene, with HR < 1 (Figure 7A).
Robust likelihood-based survival modeling and LASSO regression are widely used to screen prognostic genes in the context of high dimensional data and were therefore applied in our study. Compared with a single biomarker, integrating multiple biomarkers into one risk model could present a better prediction performance of the model. To remove the prediction errors and maintain the stability of the predictive prognosis, we specifically selected three genes (P < 0.05 and HR > 1) to develop the gene signature. Afterward, the above-selected genes with independent prognostic value, including NRNPC, ALKBH5, and ZC3H13, were screened as the candidate genes using LASSO regression. The regression coefficients based on the minimum criteria were used to assess the risk score for the GBM patients, and the coefficients of selected gene signatures were −0.014623, 0.017905, and −0.08661, respectively (Figures 7B,C).

Gene-Signature Showed Strong Associations With Clinical Features in GBM
To investigate the prognostic value of the risk gene signature in the TCGA database, GBM patients were dichotomized into lowand high-risk subgroups, based on the median risk score. We next sought to detect the correlation between the two risk subgroups and clinical features-a heatmap was designed and showed the expression of the three selected m6A RNA methylation regulators ( Figure 7D). Significant differences were observed in this heatmap between the high-and low-risk subgroups with respect to IDH1 status (P < 0.001), age (P < 0.001), molecular subtypes (P < 0.001), and RM1/2 subgroups (P < 0.001). To evaluate the predictive accuracy of the risk score model, we performed the ROC curve and calculated the AUC. The AUC was 0.701 in the 2-year ROC curve for the prognostic model ( Figure 7E). For molecular phenotypes, such as IDH1 and G-CIMP status, the risk model also performed relatively well, with AUC was 0.821 and 0.733, respectively (Figures 7F,G). A similar trend was observed in subgroup analyses for mesenchymal and proneural, with AUCs of 0.703 and 0.764, respectively (Figures 7H,I). Moreover, the predicting power of the risk score model was significantly increased in RM1/2 subgroup analysis, with an AUC of 0.887 ( Figure 7J). Notably, patients in the highrisk group exhibited significantly shorter survival time than those in the low-risk group (P < 0.05) (Figure 7H). Consistent with these findings, the patients with a high-risk score were also more sensitive to temozolomide chemotherapy, radiation therapy, and chemoradiation than low-risk score patients ( Figure S4).
We further performed univariate and multivariate Cox proportional hazard regression analyses for the TCGA dataset to determine whether the risk signature was an independent prognostic factor. By univariate Cox analysis, age (HR = 1.033, P < 0.001) and risk score (HR = 11.899, P < 0.001) were all correlated with the OS, while GBM subtype and gender were not (Figure 7L). A similar trend of risk score was also observed when including these factors in the multivariate Cox proportional hazard regression (Figure 7M). The results demonstrated that age and risk score were independent prognostic factors in the TCGA GBM dataset. According to our results, the independent prognostic value and excellent prediction accuracy of the gene signature derived from the 13 m6A RNA methylation regulators were identified.

Low Expression in Normal Brain Tissues of METTL3 and METTL14
Based on our results and the evidence in literature, METTL3 was overexpressed specifically in GBM and was significantly related to the occurrence of GBM (45)(46)(47). To comprehensively understand the function of METTL3, we retrieved the expression levels of healthy tissue and tumor tissue in different parts from the GTEx and GEPIA databases (48), respectively. We found that the expression values of METTL3 in the brain were lower than other tissues in the organism (Figure 8A). Notably, the expression level of METTL3 in most tumors was smaller than the corresponding healthy tissue, except for GBM ( Figure 8B). These results indicated that high expression of METTL3 might act as a driver of GBM and play a crucial role in GBM. Considering that METTL3 and METTL14 are both the most common and abundant mRNA modifications in eukaryotes, we also searched the expression profiles of METTL14 and found the same trend ( Figure S5). The above results provided evidence for METTL3 and METTL14 as proto-oncogenes of GBM.

DISCUSSION
In the current study, we systematically analyzed the expression of m6A RNA regulators with different clinicopathological parameters and revealed the potential values. In particular, by comparing the expression of 13 regulators in a large number of healthy tissues and primary and recurrent tumor tissues, we found that they were related to the occurrence and recurrence of GBM. Furthermore, the expression of m6A RNA methylation regulators was also associated with the GBM subtype, G-CIMP status, and IDH status. Besides, GBM samples were classified into two subgroups, RM1/2, through consensus cluster analysis based on the expression of the 13 regulators. The RM1/2 subgroup not only affected OS and clinical characteristics, but also closely related to malignancy-related processes, key signaling pathways, and GBM hallmarks. Taking advantage of CMap, we  also identified potential compounds targeting RNA methylation regulators in GBM. In addition, we obtained 79 genes with poor prognosis based on the RM1/2 subgroup by Kaplan-Meier analysis. Importantly, 64 genes with poor prognosis were validated in CGGA, a separate GBM database. Finally, we derived a prognostic gene signature, which dichotomized the OS of GBM patients into low-and high-risk subgroups and allowed us to extend the analysis. This risk gene-signature could be used as an independent prognostic marker and accurate clinicopathological parameters predictor. Glioma was divided into GBM and low-grade glioma (LGG). GBM, as the most destructive glioma (WHO: IV), possesses significantly different genomics, treatment methods, clinical manifestations, characteristics, and prognosis from LGG (WHO: I-III) (49)(50)(51)(52). Prior to this, the value of m6A methylation regulators in gliomas has been explored (53). Considering the comprehensive differences between LGG and GBM, we believed that this analysis was not sufficiently detailed and specific. Therefore, we specifically analyzed the specific value of these regulators in GBM. Similarly, we have all identified different hallmarks and pathways associated with malignancy. In particular, we further analyzed and obtained regulator-related specific targeted drugs and genes with poor prognosis, of which 37 genes have not been previously reported at the prognostic level. Especially, the PPI network between regulators and related genes was explored. Furthermore, we derived a prediction model that can predict the specific clinical characteristics and molecular phenotypes of GBM. Finally, we also provided evidence for the large transcriptome levels of METTL3 and METTL14 as cancer driver genes.
Among the m6A RNA methylation regulators, METTL3 or METTL14 is one of the most common and abundant mRNA modifications in eukaryotes. It has been reported that METTL3 or METTL14 inhibits the growth and selfrenewal of the GBM stem cells (32). ALKBH5 was reported to maintain tumorigenicity of GBM stem cells by sustaining FOXM1 expression and cell proliferation program (54), suggesting a crucial tumorigenic role. Most recently, it has been reported that FTO plays a carcinogenic role through the FTO/m6A/MYC/CEBPA signaling pathway in IDH mutant cancers, such as glioma and leukemia (19,55). The differences of involved genes among different tumor types give us a clue that altered the expression of key genes, which are sensitive to the function of m6A methylation regulators, can cause significant phenotype changes.
In this section, we rounded analyzed the expression of all m6A RNA methylation regulators in GBM at the occurrence and recurrence stages. Unlike our study, a previous trial showed that ALKBH5 was an oncogene to maintain tumorigenicity, while our study showed a significantly decreased trend in the GBM group compared with normal. However, unlike the previous trial, our study included a large number of clinical samples and was validated in two databases. This difference in the number of samples may account for the different results. Interestingly, the upward trend in ALKBH5 was significantly associated with tumor recurrence when we compared primary and recurrent tumors. It's worth mentioning that ALKBH5 belongs to the AlkB family of non-heme Fe(II)/a-ketoglutaratedependent dioxygenases, and the activity is iron-dependent (17). Given our results, we further speculate that iron metabolism is involved in GBM recurrence (56,57). However, this hypothesis needs to be more tested. Nevertheless, a tendency toward a lower expression of FTO was observed in GBM compared with normal tissues. Unlike ALKBH5, FTO was found to mediate the demethylation of m6Am instead of m6A preferentially. It seems that FTO and ALKBH5 mediate the demethylation of different methylation targets in GBM, which is worthy of future research. Based on this difference in demethylation targets and tendencies, we speculate that ALKBH5 and FTO have different functions and mechanisms in GBM and are worthy of further study.
Nearly all IDH-mutant GBMs harbored G-CIMP and patients carrying G-CIMP (G-CIMP+) have been confirmed to confer a better clinical outcome than those not carrying (G-CIMP-) (58). Collectively, we conclude that the expression of m6A RNA methylation regulators is closely associated with the occurrence, recurrence, IDH status, G-CIMP status, and molecular subtype of GBM. Moreover, these findings of the expression of each individual m6A methylation regulator can contribute to the development of new cancer therapies, as chemotherapy targeting m6A methylation is now at the forefront of cancer therapy.
We demonstrated that m6A RNA methylation regulators were also related to the biological processes, cellular component, and signaling pathways of GBM malignant progression. RNA m6A methylation is a nascent field as of yet, the significance of the above epigenetics marker in human cancer is just beginning to be appreciated. Although the m6A modification showed tissue-specific regulation and increased significantly throughout the brain development process (3), studies (59, 60) on the role of m6A modification in either brain lesions or brain cancers have only been reported sporadically (61,62). Several biological processes and signaling pathways have already been identified: tumor stem-like cell regulation, including maintenance, radio-resistance, and tumorigenesis; post-transcriptional modification, including in RNA transcript, RNA processing, RNA processing, RNA degradation, and RNA translation; FTO/m6A/MYC/CEBPA signaling pathways (19); JAK1/STAT5/C/EBP β pathways (63); and the IL-7/STAT5/SOCS pathways (64). This report provided the potential biological process and pathway between RNA m6A methylation and GBM-malignant progression, which represent a significant step toward developing therapeutic strategies to treat GBM by targeting m6A modification.
CMap can identify biomarkers for predicting specific drug reactions, mechanisms of treatment, and ways to overcome them (65)(66)(67). CMap analysis, which is based on a limited number of treated cell lines, accurately identified a number of compounds that have been shown to have an effect on m6A of other tumor types with specificity (33,(68)(69)(70). METTL3 has been reported to promote gastric cancer angiogenesis by secreting HDGF (71). This result verified the accuracy of our CMap-based drug prediction from the side. Based on these results, hence, we speculated that PDGFR tyrosine kinase receptor inhibitor, KIT inhibitor, and tubulin inhibitor could all be used as potential agents that specifically target m6A-related biological functions and pathways for subsequent research.
METTL3, served as a methyltransferase, has been reported to be essential for glioma stem-like cell maintenance and radioresistance (45). Our findings further confirm that METTL3 was a potential therapeutic target, and future research is expected to focus on studies that specifically target METTL3. Since the smallmolecule inhibitors of METTL3 have not yet been invented, future research should focus on this area (72).
This study identified and validated that 64 genes were associated with poor outcomes in GBM patients. Moreover, we were able to construct four PPI modules, all of which were related to critical GBM biological processes. Highly relevant nodes in the modules, including STAT3, SLC11A1, and ITGAM, have been reported to promote tumor proliferation, angiogenesis, migration, and invasiveness (73)(74)(75)(76)(77). Among the 64 genes validated, 27 (such as ALOX5, CAST, HS6ST1, ITGAM, PTPN6, SLC11A1, and SLC12A7) have been reported to be involved in the pathogenesis of GBM or critical in predicting OS. This suggests that our big data-based analyses using TCGA and CGGA cohorts harbor predictive value. Although the remaining 37 genes have not been previously reported to be associated with GBM prognosis, they can be used as a potential clinical prognostic indicator for GBM patients, which can facilitate clinicians to make more accurate diagnosis easily.
In this study, we attempted to introduce some concepts associated with the theory of the prognosis value of m6A RNA methylation regulators based on uncovered sets. METTL3 has been reported as a potential biomarker panel for prognostic prediction in colorectal carcinoma (78). The prognostic model of multiple m6A RNA methylation regulators for patients with GBM has not been developed. The GBM prognostic gene signature based on three selected m6A RNA methylation regulators was designed for the first time. As we observed, the risk score calculated by the correlation coefficient conferred the ability in prognosis and clinicopathological parameter prediction. Excitingly, Cox analysis results further confirmed the independent prognostic value of the risk score. Meanwhile, GBM patients with a high-risk score showed more sensitivity to temozolomide chemotherapy, radiation therapy, and chemoradiation than low-risk-score patients. These findings may deepen our understanding of m6A methylation regulators in prognosis level and tolerance to chemoradiotherapy.
To sum up, we attempted to identify the expressions, potential functions, and prognostic values of m6A RNA methylation regulators in GBM. Our study provides strategies for comprehensive analysis of cancer genomics based on consensus clustering analysis for systematic identification of specific m6Arelated targets and specific targeted drugs based on m6A RNA methylation regulators. The prognostic gene signature and genes with poor prognosis might contribute to the personalized prediction of GBM prognosis and serve as a potential biomarker reflecting GBM patients' response to therapies that specifically target m6A. Finally, further investigation of these genes could lead to novel insights into the potential association of m6A methylation regulators with GBM prognosis in a comprehensive manner.

AUTHOR CONTRIBUTIONS
JD, RX, LC, and SH conceived and designed the study and drafted the manuscript. JD and KH collected, analyzed, and interpreted the data. HJ, SMi, YB, and SMa participated in revising the manuscript. All authors have read and approved the final manuscript.