The expression profiling data (FPKM) of patients with glioblastoma multiforme (GBM) and lower-grade glioma (LGG) were downloaded from TCGA GDC official website (https://portal.gdc.cancer.gov/), and the FPKM was then converted to the TPM value. The clinical data included age, sex, and survival prognosis, and after deleting the missing information, 638 tumor tissues and five normal tissues were obtained. The copy number variation (CNV) of glioma patients was also downloaded, and the RCircos package (Zhang et al., 2013) was used to map the genetic copy number variation in 23 pairs of chromosomes. After selecting “Masked Somatic Mutation” as the somatic mutation data, we used the maftools package (Mayakonda et al., 2018) to visualize somatic mutations and obtain the tumor mutation burden (TMB) of each patient. In addition, the sequencing results of 970 glioma patients were downloaded from the CGGA database, and the sequencing results of 2,642 normal brain tissues were downloaded from the GTEx database and converted into TPM values. Finally, 2,647 normal brain tissues and 1,608 glioma tissues were sequenced.

To analyze the expression of m6A-related genes in glioma, we first analyzed the differential expression of m6A-related genes between glioma and normal tissues, the correlation of gene expression, and its influence on the prognosis of patients with glioma. Subsequently, using the expression profiling of both TCGA-glioma data and CGGA-glioma data, the m6A-related genes were incorporated into the risk model, and the least absolute shrinkage and selection operator (LASSO) algorithm was used to perform dimensionality reduction analysis and obtain prognosis-related genes. The normalized gene expression value weighted by the penalty coefficient obtained by LASSO Cox analysis established a risk score formula, and the median value of the risk score was used to divide the patients into high- and low-risk groups. risk S c o r e   = ∑ i C o e f f i c i e n t   ( h u b   g e n e i ) ∗ m R N A   E x p r e s s i o n   ( h u b   g e n e i )  

Gene ontology (GO) analysis is a common method for large-scale functional enrichment research, including biological processes (BP), molecular functions (MF), and cellular components (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database that stores information on genomes, biological pathways, diseases, and drugs. The R clusterProfiler package (Wu et al., 2021) was used to perform GO annotation and KEGG pathway enrichment analysis for the signature genes. The critical value of FDR <0.05 is considered statistically significant.

To study the differences in biological processes based on the gene expression profiling data of glioma patients, we performed gene set enrichment analysis (GSEA). GSEA calculates the statistical difference between two biological states in a specific gene set (Canzler and Hackermüller, 2020) and is usually used to estimate the changes in pathway and biological process activity in the dataset. The “h.all.v7.2. symbols.gmt” gene set was downloaded from the MSigDB database (Liberzon et al., 2015) for the analysis. Statistical significance was set at an adjusted p-value of less than 0.05.

We used the GSVA method (Ferreira et al., 2021) to analyze the correlation between different groups and biological processes. Mariathasan et al. (2018) constructed a set of genes to store those related to certain biological processes, including 1) immune checkpoints, 2) antigen processing, 3) CD8⁺ T cell characteristics, 4) epithelial–mesenchymal transition (EMT) markers, including EMT1, EMT2, and EMT3, 5) angiogenesis characteristics, 6) pan-fibroblast TGF-β response characteristics (Pan-FTBRS), 7) WNT characteristics, 8) DNA damage repair, 9) mismatch repair, (10) nucleotide excision repair, 11) DNA replication, and 12) antigen processing and presentation. The gene sets were downloaded according to different biological characteristics to calculate the enrichment scores corresponding to the patients and compare the differences between the two groups.

The R limma package (Ritchie et al., 2015) was used to analyze the differentially expressed genes (DEGs) between the high-and low-risk glioma patient groups to determine the genes associated with the m6A risk model. The DEGs were defined as those with an absolute value of log(fold change) > 1.5 and an FDR <0.01. The tumor was divided into different groups based on the Euclidean distance, and named as genecluster using the hierarchical clustering method, where the R ConsensuClusterPlus package (Wilkerson and Hayes, 2010) was used to determine the number of clusters in the dataset, and repeated 1,000 times to ensure the stability of classification. At the same time, based on the expression profile of specific genes, they were divided into two groups: signature genes A and B.

First, according to the DEG value, unsupervised clustering classified the patients in TCGA. According to the Boruta algorithm, the m6A Signature-A and B gene clusters are reduced in dimensionality, and PC1 is extracted as a score using the PCA algorithm. Finally, we applied a method similar to that of the gene expression grade index to define the ICI score of each patient: m 6 A   s c o r e   =   ∑ i P C 1 A − ∑ i P C 1 B

To further quantify the proportion of different immune cells in the glioma sample, we used a single-sample GSEA algorithm to distinguish human immune cell phenotypes in the tumor immune microenvironment (TME) with high sensitivity and specificity. The ssGSEA algorithm was used to quantify the relative content of tumor-infiltrating immune cells in patients with glioma (Minjie Wang et al., 2021). The algorithm identified 28 genes for marking different tumor-infiltrating immune cell types through the research of Bindea et al. (2013). The gene set contained various human immune cell subtypes, including CD8⁺ T cells, dendritic cells, macrophages, and regulatory T cells. The enrichment score calculated by ssGSEA analysis with the R GSVA package (Minjie Wang et al., 2021) can be used to represent the infiltration level of each immune cell type in each sample.

The ESTIMATE package (Yoshihara et al., 2013) is used to evaluate the immune activity of the tumors. ESTIMATE analysis quantifies the immune activity (immune infiltration level) in the tumor sample based on its gene expression profile and obtains an immune score for each tumor sample. This can be used to compare the differences in immune infiltration between the high- and low-risk groups of GLIOMA patients.

To analyze the copy number changes in different risk score groups of TCGA-glioma patients, we used R’s TCGAbiolinks package to download the masked copy number segment data of patients. GISTIC 2.0 analysis was performed on the downloaded CNV fragments using GenePattern5. During the GISTIC 2.0 analysis, the default settings were used, except for a few parameters (for example, the confidence level was 0.99; the X chromosome was not excluded before the analysis). Finally, we used R’s Maftools package to visualize the analysis results.

The Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) is a public database for molecular cancer therapy and mutation exploration (Yang et al., 2013). We used R’s pRRophetic package (Geeleher et al., 2014) to download cell line gene mutation data and IC₅₀ values of different anti-cancer drugs. We then analyzed the correlation between patients with high- and low-risk scores and the sensitivity to different anti-cancer drugs.

Univariate and multivariate Cox analyses were used to analyze the risk score combined with the patient’s clinicopathological characteristics to predict the overall survival (OS) to prove that the risk score combined with clinicopathological characteristics can evaluate the prognosis for the individual patient. Subsequently, the risk scoring model combined with clinicopathological characteristics was selected to construct a clinical prediction nomogram. Harrell’s consistency index (C-index) was measured to quantify the discrimination performance. A calibration curve was generated to evaluate the performance of the nomogram by comparing its predicted value with the actual survival rate.

To validate m6A gene expression, we performed immunohistochemistry (IHC) on surgical human glioma samples. Using available antibodies, we selected three genes of interest: IGF2BP3 (14642-1-AP, Proteintech, China), RBM15B (22249-1-AP, Proteintech, China), and RBM15 (10587-1-AP, Proteintech, China). The comparison was performed between the glioma sample and the para-tumor area, performed under the approval of local medical ethics (No.2020SQ119) in Shanghai General Hospital, affiliated to Shanghai Jiaotong University.

All data processing and analyses were performed using the R software (version 3.6.2). To compare two groups of continuous variables, the statistical significance of normally distributed variables was estimated using the independent Student’s t-test, and the differences between non-normally distributed variables were analyzed using the Mann–Whitney U test (i.e., Wilcoxon rank sum test). Chi-square or Fisher’s exact tests were used to compare and analyze the statistical differences between the two groups of categorical variables. Pearson correlation analysis calculated the correlation coefficients between different gene sets. The survival package in R was used for survival analysis. The Kaplan–Meier survival curve was used to show the survival difference. The log-rank test was employed to determine the significance of different survival times between the two groups. Univariate and multivariate Cox analyses were used to determine independent prognostic factors. All statistical p-values were two-sided, and statistical significance was set at p < 0.05.

To analyze the overall expression of m6A-related genes in patients with glioma, we analyzed genomic mutations and mRNA expression, including gene expression levels, single nucleotide polymorphisms, and copy number variations. First, we conducted a comprehensive analysis of the expression in gliomas and normal brain tissues in TCGA, CGGA, and GTEx databases and used the de-batch method. PCA results showed that the characteristics of m6A-related genes differed between glioma and normal brain tissues (Figure 1A). Subsequently, the differential analysis showed that, between glioma and normal brain tissue, a variety of m6A-related genes were differentially expressed, including METTL14, METTL16, ZC3H13, YTHDC1, YTHDC2, YTHDF2 (Zhang et al., 2017; Dixit et al., 2021; Fang et al., 2021)etc. (Figure 1B).

The results of single nucleotide polymorphism (SNP) analysis showed that among GBM samples, 12 had single nucleotide mutations in m6A-related genes, among which the mutation rate of the ZC3HI3 gene was the highest. In contrast, in the LGG samples, 16 had single nucleotide mutations in m6A-related genes, and the mutation rate of METTL3 was the highest (Figure 1C). Moreover, studies on the frequency of CNV changes have shown that m6A-related gene changes in CNV levels in glioma patients are common, and most of them are concentrated on copy number loss (Figure 1D).

The heat map resulting from Pearson’s analysis revealed a positive correlation between m6A-related gene expression and glioma tissue (TCGA dataset) (Figure 2A). The detailed number of each coefficient is displayed in Supplementary data. We further analyzed the influence of m6A-related genes on the prognosis of patients with glioma in TCGA and CGGA databases. The gene network depicts the interaction of m6A-related genes in glioma and their impact on the overall survival of glioma patients (Figure 2B).

Subsequently, we integrated the expression of m6A-related genes to construct a risk-scoring system to quantify the impact of m6A-related genes on the prognosis of each glioma patient. First, m6A-related genes were included in the LASSO Cox analysis, and 15 genes with the best prognostic value were obtained (Figures 2C,D). Based on the penalty coefficients of important characteristic genes calculated by LASSO Cox analysis, the gene expression and the corresponding coefficients were multiplied, and the final risk score of each sample was calculated. The distribution of risk scores, survival status, and expression patterns of the feature genes is shown in Figure 2E. Kaplan-Meier analysis showed that the overall survival (OS) of patients with high-risk scores was relatively poor (log-rank p < 0.001; Figure 2F). At the same time, the differential analysis results showed significant differences in the expression of m6A-related genes between the high- and low-risk models (Figure 2G).

Based on the median value of the m6A risk score of glioma patients, we placed the patients into the high- or low-risk group and assessed the changes in biological function between the two groups. The GSVA method was used to determine the enrichment scores of these patients, and heat maps were used to show the relevant signaling pathways and analyze their variations in the two groups (Figure 3A). In addition, the results showed that there were significant differences in the enrichment of certain related biological pathways, such as CD8⁺ T cell effectors, immune checkpoints, EMT pathways, angiogenesis, and others (p < 0.05; Figure 3B).

The limma package was used to analyze DEGs between different risk models, and 443 genes were obtained to determine the potential biological characteristics of different m6A-related phenotypes. Subsequently, based on the expression of DEGs, unsupervised clustering was used to divide glioma patients into three subtypes: genecluster-A, -B, and -C (Figure 4A). At the same time, tSNE analysis showed certain differences in gene expression levels between genes A, B, and C (Figure 4B). The heat map shows the gene expression characteristics of the three genotypes (Figure 4C). Meanwhile, the survival analysis results showed that there were significant differences in the prognosis of patients with the three different genotypes, among which the patients in the Genecluster-A group had the worst prognosis (log-rank p < 0.001; Figure 4D).

According to the expression and correlation of DEGs in these groups, the genes were divided into m6A signature-A and m6A signature-B. There were 268 m6A signature-A and 51 m6A signature-B gene sets. To explore the differences in biological functions between the two groups, we conducted a functional enrichment analysis. KEGG enrichment analysis showed that signature genes A and B showed different, unique biological processes (Tables 1, 2). m6A gene-A is involved in the synaptic vesicle cycle, insulin secretion, nicotine addiction, and GABAergic synapse pathways (Figure 4E), while the gene set overexpressing signature gene-B mainly manifests as the relaxin signaling pathway, AGE-RAGE signaling pathway in diabetic complications, ECM-receptor interaction, and protein digestion and absorption pathways (Figure 4F).

We constructed a new prognostic-related risk-scoring system to better predict the impact of m6A features on patient prognosis. According to the expression of m6A signature genes A and B in glioma patients, principal component analysis was used to calculate the corresponding PCA1 of each patient, and the corresponding m6A score was obtained by subtraction and named m6A group. Similarly, based on the median score of the prognostic models, patients were divided into high- and low-risk groups. The Sankey diagram shows the correspondence between the gene cluster corresponding to each glioma patient, prognostic model of the m6A group, and patient survival status (Figure 5A). At the same time, the results of the survival analysis showed that the prognostic score model could predict well the OS of glioma patients (log-rank p < 0.001; Figure 5B).

Subsequently, we analyzed the biological effects of the high and low m6A groups. GSEA showed that hallmark oxidative phosphorylation, adipogenesis, hedgehog signaling and MYC target V1 were significantly enriched in glioma patients in the low m6A group (Table 3; Figures 5C–H).

Subsequently, we evaluated the effect of the m6A risk score on the genetic variation in glioma patients, including SNPs and copy number variations (CNVs). The single nucleotide mutation analysis of driver genes in common tumorigenesis showed that their SNP levels differed between the high and low groups (Figures 6A,B). At the same time, the overall analysis showed that there was no significant difference in the tumor mutation burden (TMB) between the high and low m6A groups of GBM patients (p = 0.73; Figure 6C), while there was a significant difference in the LGG patient group (p < 0.001; Figure 6D). Moreover, research on the frequency of CNV changes showed that in patients in the high m6A group, the CNV changes were mainly in the deletion of gene copy numbers. In contrast, in patients in the low m6A group, they reflected gene amplification (Figure 6E).

Next, we evaluated the effect of the m6A risk score on the overall immune characteristics and the different levels of immune cell infiltration in glioma patients. The results showed that compared with the low-risk group, the immune-related and stromal-related scores of patients in the high-risk group were significantly increased (p < 0.001, Figure 7A). At the same time, we further used the ssGSEA algorithm to evaluate the infiltration level of 28 different immune cells (Figure 7B). Differential analysis showed that the infiltration levels of multiple immune subgroups were significantly different between the high- and low-risk groups (Figure 7C), including CD8+T cells, activated memory CD4 +T cells, follicular helper T cells, and M1 macrophages. Further analysis showed that the expression levels of multiple HLA family genes and immunotherapy-related targets were significantly different between the high and low m6A groups (Figures 7D,E).

To analyze glioma patient sensitivity to small molecule drugs based on the m6A risk score, we downloaded the cell line gene mutation data and IC₅₀ values of different anti-cancer drugs from the GDSC database. Based on the responsiveness of the cell lines to 138 different chemotherapeutics and small molecule anti-cancer drugs, the IC₅₀ values of glioma patients for different drugs were predicted. The results showed significant differences between patients with high and low m6A risk scores (p < 0.001; Figure 8), especially Nutlin.3a (p53 activator) (Awan et al., 2021), EHT. 1864 (Rac GTPase inhibitor), (Onesto et al., 2008), and BIRB. 0,796 (pan p38MAPK inhibitor) (Tang et al., 2018).

Subsequently, we assessed the impact of the m6A risk score on the prognosis of patients with glioma. Univariate and multivariate Cox analyses showed that the m6A risk score was an independent risk factor for glioma patient prognosis prediction (Table 4; Figure 9A). The m6A group was combined with different clinicopathological characteristics to construct a nomogram to predict the OS of the patient (Figure 9B). We used the C-Index to calculate the discriminative ability of the nomogram, which showed a high degree of discrimination (0.717 (0.701–0.733)). Moreover, the calibration curve shows that by comparing the one-, two-, and three-year OS estimates, the actual values observed in patients are in agreement (Figures 9C–E).

To validate the expression of m6A genes in glioma, we performed IHC of patient samples. We found that the expression of IGF2BP3 and RBM15B (Figure 10A,B) was increased in the glioma area compared to that in the para-tumor area. In contrast, the expression of RBM15 was obviously reduced in the tumor area (Figure 10C). These findings are consistent with the LASSO Cox model (Table 5).