Submit an Article
Online ISSN: 1945-4589
Impact Journals, LLC. is the publisher of Aging: www.impactjournals.com.
Impact Journals is a member of the Wellcome Trust List of Compliant Publishers.
Impact Journals is a member of the Society for Scholarly Publishing.
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 15, Issue 10 pp 4051—4070

Integrating machine learning and bioinformatics analysis to m6A regulator-mediated methylation modification models for predicting glioblastoma patients’ prognosis and immunotherapy response

Chuanyu Li1, *, , Wangrui Liu1,2, *, , Chengming Liu3, *, , Qisheng Luo1, , Kunxiang Luo1, , Cuicui Wei4, , Xueyu Li1, , Jiancheng Qin1, , Chuanhua Zheng1, , Chuanliu Lan1, , Shiyin Wei1, , Rong Tan1, , Jiaxing Chen1, , Yuanbiao Chen1, , Huadong Huang1, &, , Gaolian Zhang5, , Haineng Huang1, &, , Xiangyu Wang6, ,

  • 1 Department of Neurosurgery, Affiliated Hospital of Youjiang Medical University for Nationalities, Baise 533000, Guangxi, China
  • 2 Department of Interventional Oncology, Renji Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai 200127, China
  • 3 Department of Neurosurgery, The First Affiliated Hospital of Soochow University, Suzhou 215006, Jiangsu, China
  • 4 Department of Outpatient, Affiliated Hospital of Youjiang Medical University for Nationalities, Baise 533000, Guangxi, China
  • 5 Department of Neurosurgery, The First Affiliated Hospital of Guangxi University of Chinese Medicine, Nanning 530023, Guangxi, China
  • 6 Department of Neurosurgery, The First Affiliated Hospital of Jinan University, Guangzhou 510632, Guangdong Province, China
* Equal contribution

Received: August 11, 2022       Accepted: November 30, 2022       Published: May 23, 2023      

https://doi.org/10.18632/aging.204495
How to Cite

Copyright: © 2023 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Epigenetic regulations of immune responses are essential for cancer development and growth. As a critical step, comprehensive and rigorous explorations of m6A methylation are important to determine its prognostic significance, tumor microenvironment (TME) infiltration characteristics and underlying relationship with glioblastoma (GBM).

Methods: To evaluate m6A modification patterns in GBM, we conducted unsupervised clustering to determine the expression levels of GBM-related m6A regulatory factors and performed differential analysis to obtain m6A-related genes. Consistent clustering was used to generate m6A regulators cluster A and B. Machine learning algorithms were implemented for identifying TME features and predicting the response of GBM patients receiving immunotherapy.

Results: It is found that the m6A regulatory factor significantly regulates the mutation of GBM and TME. Based on Europe, America, and China data, we established m6Ascore through the m6A model. The model accurately predicted the results of 1206 GBM patients from the discovery cohort. Additionally, a high m6A score was associated with poor prognoses. Significant TME features were found among the different m6A score groups, which demonstrated positive correlations with biological functions (i.e., EMT2) and immune checkpoints.

Conclusions: m6A modification was important to characterize the tumorigenesis and TME infiltration in GBM. The m6Ascore provided GBM patients with valuable and accurate prognosis and prediction of clinical response to various treatment modalities, which could be useful to guide patient treatments.

Introduction

Glioma is the highest incident primary brain tumor, accounting for ~80% of all adult brain cancer cases [1]. Glioblastoma multiforme (GBM) is classified as a grade IV cancer by the World Health Organization (WHO) and comprises 57.3% of all gliomas [2]. Its clinical presentation varies depending on its location and size and can be characterized by headaches, seizures, neurological dysfunctions, etc. [3]. Currently, little is known about its etiology [4]. Although multiple treatments have been developed for treating the patients, including surgery, chemotherapy, radiation and others, their survival rate remains far from satisfactory [5]. In addition, despite the promising efficacies of immunotherapy, most patients do not respond well [6, 7], urging the need for better treatments [8].

GBM develops in a highly complex tissue environment and depends on the continuous growth, invasion and metastasis of these environments, which is why it is constantly infiltrating and challenging to cure. Targeting the tumor microenvironment (TME) was shown to potentially decrease the risk of drug resistance and disease recurrence. TME has a variety of capabilities to induce the beneficial and unfavorable consequences of tumorigenesis, so how to destroy and promote the occurrence of GBM The TME is a challenging job [9–11].

N6-methyladenosine (M6A) methylation, a commonly observed modification on RNAs, is similar to DNA and histone modifications [12–15]. Its involvement has been reported in all the RNA stages life cycle, including regulation of RNA processes, translation, RNA degradation etc. [16–19]. Although recent discoveries indicated its association with cancer occurrence and development [17, 20, 21], its actual roles and underlying mechanisms remain unknown.

Recent literature reported that m6A was related to immunotherapy, which affects TME and its related immune cells by regulating targeted RNA. Therefore, m6A has become a potential target for immunotherapy. m6A may have different roles in different tumors, but the m6A detector has been shown to have broad significance in specific cancers [22–24]. In several cancers, including kidney, lung, gastric and other cancers [25–31], m6A-related signals have been identified as tumor immunity Phenotype and biomarker of anti-PD-1 immunotherapy response. These findings indicated that M6A modifications could mirror the TME status and anticipate immunotherapeutic effects in numerous cancers, not limited to specific cancer types. Another study found that m6A is inseparable from low-grade gliomas [32]. However, due to technical limitations, most studies focused on 1 or 2 m6A-related genes. TME is characterized by multiple highly coordinated components in a network, so there is no research on combining m6A and TME. In addition, deeper investigations on m6A methylation modification in GBM are needed to improve treatment outcomes.

In this study, we used 19 GBM-related m6A regulatory factors to analyze m6A modification patterns in GBM samples from the Chinese Glioma Genome Atlas (CGGA; http://www.cgga.org.cn/) and Cancer Genome Atlas (TCGA) databases, which were then validated in our own cohort. Next, we proposed a nomogram for quantifying the m6A modification patterns of GBM patients and predicting their potential immunotherapy response and prognosis. Altogether, our results showed that m6A modification promoted GBM progression and the potential clinical significance of our scoring system to guide the treatment and estimate the survival of GBM patients.

Materials and Methods

Data

The mRNA expression profile data and sample copy number variation (CNV) of Caucasian GBM samples are downloaded from the University of California, Santa Cruz (https://xenabrowser.net/datapages/). The expression profiles of Chinese GBM samples were obtained at the CGGA. In addition, were downloaded from The GEO database (https://www.ncbi.nlm.nih.gov/geo/) were queried to obtain the expression profiles of three groups of glioblastoma samples, namely GSE83300 [33], GSE74187 [34, 35]. For clinical information and data, the R package cgdsr and [36] are used for processing. Table 1 shows the retrieved information from the datasets. Perform consistency processing on the above data, including z-score [37] and batch correction [38, 39].

Table 1. The retrieved information from the datasets.

Data setClassificationNumber of samples
TCGAmRNA-seq15
AgilentG4502A_07_199
AgilentG4502A_07_2460
CGGACGGA.mRNAseq_325139
CGGA.mRNAseq_693249
CGGA.mRNA_array_30184
GEOGSE8330050
GSE7418760
GSE4337850

Unsupervised clustering

Expression data on 21 m6A regulators was extracted from 1,206 data in all data sets to determine m6A modification patterns. Some data set did not detect the expression of IGF2BP1 and METTL14, so the final expression extracted was 19 regulatory factors. These 19 m6A regulatory factors comprised seven writers (METTL3, RBM15, RBM15B, WTAP, KIAA1429, CBLL1 and ZC3H13), two erasers (ALKBH5 and FTO) and ten readers (YTHDC1, YTHDC2, YTHDF1), YTHDF2, NNPHHDF2MR, LRPPRC, ELAVL1), and based on their expression levels, various m6A modification patterns were determined and patients’ data were also analyzed [40].

Gene set variation analysis (GSVA) and single sample gene set enrichment analysis (ssGSEA)

GSVA enrichment analyses were conducted for determining the difference of m6A modification mode using c2.cp.kegg.v6.2 from the MSigDB database (https://www.gsea-msigdb.org/gsea/index.jsp) [41]. ssGSEA analyses were conducted to determine the ratio of 24 immune cells in GBM [42], the Wilcox test for comparing different samples, and Cox regression analysis for comparing their survival.

Differentially expressed genes (DEGs) between m6A clusters

The TCGA, CGGA and GEO datasets were grouped into two categories according to the 19 m6A gene expressions (Supplementary Tables 1, 2). Using the R “limma” package, we identified DEGs among the different m6A clusters [43] based on P<0.05 and difference multiple >1.5 times or <0.67 times.

m6A score calculation

Using the DEGs identified above, we implemented the random forest method to eliminate redundancy and perform survival analysis (P value <0.05). The DEGs were classified into two groups (positive or negative coefficients) through cox regressions. Refer to GGI scoring based on the equation to calculate m6Ascore [44].

m6Ascore=scale (∑X−∑Y)

Scale represents the standardization process, while X and Y represent the gene set expressions using a positive and negative Cox coefficient, respectively.

surv_cutpoint function was used for identifying the best threshold point (cutoff=-0.9884624) [45] to divide the samples into a high and low m6A score, performing correlation analysis and survival estimation.

Correlation between m6A score and other biological processes

GSVA analyses were conducted for quantifying the biological functions of the samples, and Pearson correlation analysis on the m6Ascore and ES scores of these biological functions to reveal underlying biological pathways involved in m6Ascore.

GISTIC (Q≤0.05, 95%CI) was computed to determine the common copy number change area in the samples [46].

The R “pRRophetic” package was used to obtain the IC50 of Cisplatin and Gemcitabine based on the expression profiles and to compare the difference in IC50 between the high and low m6A scores.

The TIDE tool (http://tide.dfci.harvard.edu/) was used to assess the clinical effect of immunotherapy. A higher TIDE score indicated poor immunotherapy efficacy and prognosis. Among five cancers with immune dysfunctions and rejection characteristics, only melanoma cases treated with anti-PD1 or anti-CTLA4 were available in public databases. Immune checkpoint treatment prognosis prediction in this analysis is completed by TIDE score.

Consent for publication

All authors agree to publish.

Results

Genetic variation of m6A regulatory factor of GBM in TCGA database

Here, we aimed to investigate the m6A regulators’ genetic backgrounds and variations in GBM. First, we selected all tumor samples in the TCGA database and analyzed the 21 m6A regulators, the eight writers, two erasers and 11 readers. First, we studied the mRNA levels of the m6A regulatory factors between TCGA-derived tumors and normal samples, which showed greater expressions in the tumor samples, including METL3 and WTAP (p<0.01) (Figure 1A). This shows that the role of m6A, in general, is to promote tumor growth.

Genetic variation of m6A regulator. (A) m6A expression of tumor and normal samples in TCGA; (B) Distribution of m6A gene mutations and different mutation types; (C) CNV incidence of m6A gene, blue indicates deletion and orange indicates amplification; (D) Position of m6A gene on chromosome; (E) PCA results of m6A gene in TCGA samples).

Figure 1. Genetic variation of m6A regulator. (A) m6A expression of tumor and normal samples in TCGA; (B) Distribution of m6A gene mutations and different mutation types; (C) CNV incidence of m6A gene, blue indicates deletion and orange indicates amplification; (D) Position of m6A gene on chromosome; (E) PCA results of m6A gene in TCGA samples).

Then, the CNV and somatic mutation rate of 20 m6A regulatory factors in GBM samples were investigated. HNRNPC has the highest mutation frequency, reaching 5% (Figure 1B). CNV mutations are generally changed in regulatory factors. The copy number of some genes is amplified. The frequency of CNV deletion is noticeable in genes such as HNRNPC, METL3 and ZC3H13 (Figure 1C and Supplementary Tables 3, 4). In addition, m6A regulated The position of the organ on chromosomes (Figure 1D). PCA analysis showed that tumors and normal samples could be distinguished (Figure 1E).

Unsupervised clustering of m6A genes in GBM samples

Due to the lack of IGF2BP1 and METTL14 expression levels in some data sets, we determined the consistency clustering of the m6A genes and m6A gene single-factor Cox regression using gene expression profile data of 19 m6A regulators and the survival data from the TCGA, CGGA and GEO data sets. Findings from the m6A regulatory networks indicated the associations between the m6A regulatory factors (Figure 2A) and that between regulatory factors and prognosis. The influence of m6A regulators on patients’ survival was demonstrated in Supplementary Table 5_celluar-types.tsv. Moreover, we found that m6A of both the same functional category and different functional categories showed a significant correlation (spearman); the statistical results are shown in Supplementary Table 5_celluar -interactions. tsv.

Unsupervised clustering of m6A genes in low-grade glioma samples. (A) Interaction between m6A genes. The size of the circle indicates the impact of each gene on survival prediction, and the larger the expression, the more relevant the prognosis. In the circle the green dots in the circle indicate prognostic protective factors, and the black dots in the circle indicate prognostic risk factors. The lines connecting genes show their interactions. The negative correlations are marked in blue and positive correlations in red. Gene clusters ABC are marked respectively in blue, red and brown; (B) Consistent clustering of m6A genes; (C) Kaplan-Meier curve showed no significant survival differences in two m6Aclusters; (D) GSVA enrichment analysis, showing the biological pathways with different m6Aclusters Activation state. Heat map is used to visualize these biological processes, red means activation and blue means inhibition; (E) Distribution of immune infiltration of 28 immune cells in two m6Aclusters; (F) Differential cell prognosis analysis.

Figure 2. Unsupervised clustering of m6A genes in low-grade glioma samples. (A) Interaction between m6A genes. The size of the circle indicates the impact of each gene on survival prediction, and the larger the expression, the more relevant the prognosis. In the circle the green dots in the circle indicate prognostic protective factors, and the black dots in the circle indicate prognostic risk factors. The lines connecting genes show their interactions. The negative correlations are marked in blue and positive correlations in red. Gene clusters ABC are marked respectively in blue, red and brown; (B) Consistent clustering of m6A genes; (C) Kaplan-Meier curve showed no significant survival differences in two m6Aclusters; (D) GSVA enrichment analysis, showing the biological pathways with different m6Aclusters Activation state. Heat map is used to visualize these biological processes, red means activation and blue means inhibition; (E) Distribution of immune infiltration of 28 immune cells in two m6Aclusters; (F) Differential cell prognosis analysis.

Our findings indicated that interactions between m6A regulators of different functional categories played essential roles for forming different m6A modifications in GBM. Then, the expressions of the regulators from TCGA, CGGA and GEO were obtained, and we utilized the R “ConsensusClusterPlus” package to conduct consistency clustering, based on which m6A clusters A and B were identified (Figure 2B) and further analyses showed that the prognosis between them was similar (P=0.23, Figure 2C).

Functional annotations and TME infiltration characterization of the m6A clusters

GSVA enrichment analysis was conducted to assess the differences between the biological behaviors of the regulatory factors in the two m6A modification subgroups. m6A cluster A was found to be significantly enriched in cell cycle pathways, while m6A cluster B in biosynthesis and oxidative phosphorylation of glycosphingolipids pathways (Figure 2D, Significant analysis results are shown in Supplementary Tables 6, 7).

Furthermore, ssGSEA analysis using the RNA-seq data was conducted to determine the proportions of 28 immune cells, such as B memory cells, activated Dendritic cells and M0 macrophages, in GBM samples (Figure 2E) and the findings showed significantly different distributions of immune cells abundance between the two subgroups. Figure 2F shows the results of single-factor Cox regression analysis of different proportions of immune cells between m6A clusters A and B (Analysis results are shown in Supplementary Table 8).

In the TCGA database, we found that the mutation difference of IDHI (chi-square test, P=0.1998181), TP53 (chi-square test, P=0.4284502) and EGFR (chi-square test, P=0.9538232) in subgroups m6AclusterA and m6AclusterB were not significant (Figure 3A, specific information is shown in Supplementary Table 9). No obvious differences were observed between their clinical characteristics (Figure 3B).

Comparative analysis between m6Acluster in the TCGA dataset. (A) The distribution of IDH1, EGFR, TP53 mutations in the 2 m6Aclusters; (B) The distribution of radiotherapy, gender, and age in m6Acluster; (C) The enrichment scores of different m6Acluster groups difference (**PPP

Figure 3. Comparative analysis between m6Acluster in the TCGA dataset. (A) The distribution of IDH1, EGFR, TP53 mutations in the 2 m6Aclusters; (B) The distribution of radiotherapy, gender, and age in m6Acluster; (C) The enrichment scores of different m6Acluster groups difference (**P<0.05, *** P<0.01, ****P<0.001).

GSVA analysis was carried out with gene set constructed by Mariathasan; there is a significant difference of enrichment scores between the m6Acluster groups (Figure 3C, analysis results are shown in Supplementary Table 10). Additionally, the m6A regulatory factor distributions in m6AclusterA and m6AclusterB are shown in Figure 4 (Analysis results are shown in Supplementary Table 11).

The expression of m6A regulatory factors in m6Acluster.

Figure 4. The expression of m6A regulatory factors in m6Acluster.

m6A-related DEGs and constructions of the m6Agenecluster

The biological behaviors of each m6Acluster were detected in 73 phenotype-related differential genes with a limma software package (Supplementary Table 12).

Furthermore, we implemented an unsupervised cluster analysis using the phenotype-related genes of m6A to classify patients into separate genomic subgroups according to different genomic subtypes based on their m6A-modified genomic phenotypes, m6A gene cluster A and m6A gene cluster B (Figure 5A and Supplementary Table 13). The findings showed the prognosis of m6A gene cluster A tumor was better than m6A gene cluster B, and part of the m6A regulatory factors expression levels were markedly greater compared to m6A gene cluster B (Figure 5C).

Comparison between m6Agenecluster. (A) Unsupervised clustering of m6A phenotype-related genes in low-grade glioma samples. The samples are divided into different genomic subtypes, called m6AgeneclusterA and m6AgeneclusterB; (B) Kaplan-Meier curve indicates that m6A modifies the genome table Type has an obvious relationship with overall survival rate; (C) Expression of 19 m6A genes in 2 gene clusters. The upper and lower ends of the box indicate the interquartile range of values. The line in the box indicates the median value, and the black dots indicate outliers. The T test is used to test the statistical differences between gene clusters (**PPP

Figure 5. Comparison between m6Agenecluster. (A) Unsupervised clustering of m6A phenotype-related genes in low-grade glioma samples. The samples are divided into different genomic subtypes, called m6AgeneclusterA and m6AgeneclusterB; (B) Kaplan-Meier curve indicates that m6A modifies the genome table Type has an obvious relationship with overall survival rate; (C) Expression of 19 m6A genes in 2 gene clusters. The upper and lower ends of the box indicate the interquartile range of values. The line in the box indicates the median value, and the black dots indicate outliers. The T test is used to test the statistical differences between gene clusters (**P<0.05, ***P<0.01, ****P<0.001).

Analysis of m6Ascore

The Random Forest method was implemented to remove the redundancy in the DEGs and identified significant genes associated with the classification (Supplementary Table 14). Cox regression was performed to establish the association between the significant genes and GBM patients’ survival. The genes’ coefficient values were used to separate them into two groups, and all samples were scored with the m6Ascore formula. Finally, based on the surv_cutpoint function of R package survminer, the optimal threshold of m6score was determined (cutoff=-0.9884624) as a standard to divide the samples into m6Ascorehigh group and m6Ascorelow group (Figure 6A and Supplementary Table 16). Our results indicated that the m6Ascorelow group had significantly better survival than the m6Ascorehigh group (P<0.0001), suggesting that the m6Ascore could be used to accurately characterize GBM patients’ prognoses (Figure 6B).

Establishment of m6Ascore. (A) Alluvial plot showing the changes of m6A cluster, gene cluster and m6Ascore; (B) Kaplan-Meier curve shows that m6Ascore high and low grouping has a significant relationship with overall survival rate; (C) Using Pearson analysis, the correlation between m6Ascore and known gene features in GBM. Negative correlation is marked in blue, and it is positively correlated with red. X in the figure indicates that the correlation is not significant, and the larger the circle, the more significant; (D) The distribution of the enrichment scores of known gene features in the m6Ascore high and low group samples in the TCGA+CGGA+GEOdata set (***PPE) The distribution of m6Ascore in m6Acluster (****PF) Distribution of m6Ascore in m6Agenecluster (****P

Figure 6. Establishment of m6Ascore. (A) Alluvial plot showing the changes of m6A cluster, gene cluster and m6Ascore; (B) Kaplan-Meier curve shows that m6Ascore high and low grouping has a significant relationship with overall survival rate; (C) Using Pearson analysis, the correlation between m6Ascore and known gene features in GBM. Negative correlation is marked in blue, and it is positively correlated with red. X in the figure indicates that the correlation is not significant, and the larger the circle, the more significant; (D) The distribution of the enrichment scores of known gene features in the m6Ascore high and low group samples in the TCGA+CGGA+GEOdata set (***P<0.01, ****P<0.001); (E) The distribution of m6Ascore in m6Acluster (****P<0.001); (F) Distribution of m6Ascore in m6Agenecluster (****P<0.001).

Correlation assessment between m6Ascores and known gene features indicated a significant positive correlation between m6Ascore and biological functions like EMT2 and immune checkpoint (Figure 6C and Supplementary Table 15). Wilcox test showed significant differences between m6Acluster and m6Agenecluster in m6ascore (Figure 6E, 6F). Further, the m6AclusterA and m6AgeneclusterA had significantly superior m6Ascore compared with the other groups.

Furthermore, in the TCGA cohort, a significant difference in m6A scores between the subgroups, such as IDH1 mutation status and TP53 mutations, were found.

Statuses, Wilcox test P value was 1.4e-09 and 0.0015) (Figure 7A, 7B). In addition, the m6Ascorelow subgroup had a significantly better prognosis than the m6Ascorehigh subgroup of GBM (P<0.0001; Figure 7C).

Comparative analysis and model verification of m6Ascore in TCGA dataset. (A, B) The distribution of m6Ascore in different subgroups; (C) there was significant difference in survival between m6Ascorehigh group and m6Ascorelow group in TCGA samples.

Figure 7. Comparative analysis and model verification of m6Ascore in TCGA dataset. (A, B) The distribution of m6Ascore in different subgroups; (C) there was significant difference in survival between m6Ascorehigh group and m6Ascorelow group in TCGA samples.

Differential molecular characteristics in m6Ascorehigh and m6Ascorelow group

Here, we compared the m6Ascorehigh with the m6Ascorelow groups in the TCGA dataset. Differences in somatic mutations were determined using the R “maftools” package. Significant alterations were observed in the frequency of TTN (m6Ascorehigh, 34%; m6Ascorelow, 28%), TP53 (m6Ascorehigh, 32%; m6Ascorelow, 62%) and MUC16 (m6Ascorehigh, 27%; m6Ascorelow, 28%) genes (Figure 8A, 8B). Figure 8C, 8D shows the distribution in CNV regions between the m6Ascorehigh and m6Ascorelow groups.

Analysis of molecular characteristics of m6Ascore high and low groups. (A, B) Distribution of gene mutations in samples of m6Ascore high and low groups; (C, D) The distribution of copy number amplification and deletion regions in the sample set of m6Ascore high and low groups.

Figure 8. Analysis of molecular characteristics of m6Ascore high and low groups. (A, B) Distribution of gene mutations in samples of m6Ascore high and low groups; (C, D) The distribution of copy number amplification and deletion regions in the sample set of m6Ascore high and low groups.

m6Ascore predicted chemotherapy and immunotherapy responses of GBM patients

In TCGA, CGGA and GEO datasets, we further estimated and compared the IC50 values of Cisplatin and Gemcitabine according to the expression levels of the m6Ascorehigh and m6Ascorelow groups using the R “pRRophetic” package, and found that a significantly greater IC50 value in the m6Ascorelow group compared with the m6Ascorehigh group, which suggest that patients in the m6Ascorehigh group had poorer drug resistance (Cisplatin: P<2.2e-16, Gemcitabine: P=8.1e-13; Figure 9A, 9B).

m6Ascore provided predictive outcomes for GBM patients receiving immunotherapies. (A, B) The difference between the IC50 values of Cisplatin and Gemcitabine in the samples of the high-risk group and the low-risk group; (C) The difference of TIDE score between samples of high-risk group and low-risk group; (D) Prediction of immunotherapy response results by using high and low groups of m6Ascore.

Figure 9. m6Ascore provided predictive outcomes for GBM patients receiving immunotherapies. (A, B) The difference between the IC50 values of Cisplatin and Gemcitabine in the samples of the high-risk group and the low-risk group; (C) The difference of TIDE score between samples of high-risk group and low-risk group; (D) Prediction of immunotherapy response results by using high and low groups of m6Ascore.

Meanwhile, in 154 samples with sequencing data (including 139 samples with chip data) from the TCGA database, a numerically greater TIDE score was observed for the m6Ascorehigh group compared with the m6Ascorelow group (Figure 9C), but the difference is not significant. The AUC of the ROC curve for the m6score in predicting immunotherapy response was 0.53.

Discussion

Multimodal therapies are the standard treatments for GBM, which include surgery resection, radiation therapy, and systemic treatment with temozolomide [47]. However, traditional treatment methods often have serious side effects, and the treatment effect is not very good, thereby cannot significantly improve patients’ survival [48]. We explored the molecular and clinical characteristics between different subgroups based on the 19 m6A genes through consistent clustering. Further, by dimensionality reduction in m6A-related DEGs between the subgroups, we determined their association with patient prognosis. Then m6Ascore of these genes were calculated, and the differences between the m6Ascorehigh group and m6Ascorelow group were compared.m6Ascore based classification was used for prediction of immunotherapy response. There exist some advantages of this study. Significant differences are shown in prognosis, clinical features or molecular characteristics between the groups of m6Agenecluster and m6Ascore. There are also some limitations. First, in the prediction of immunotherapy response, AUC of ROC curve is comparably low, only reached 0.53. Second, the difference of m6Acluster in prognosis and clinical characteristics were not significant. The clinical features of this part of glioblastoma did not find the status information of 1p19q. By consulting the data, this feature belongs to GBM.

In the past decades, many new concepts were investigated as an attempt to improve the treatments of GBM, including genetic testing, electromagnetic field therapy, and function-guided resection [2, 49]. However, the underlying causes of GBM malignant progression have not been explored, so there are still no biomarkers to accurately predict the prognosis and exacerbation of GBM patients [32, 50].

The TME and purity of tumor cells have important roles in tumor invasion and progression [51]. Therefore, by comprehensively analyzing the properties of TME and the cells it recruits in GBM, the immunophenotype of GBM at various stages can be identified, so as to find accurate biomarkers and discover new and effective therapeutic targets [52, 53].

Determining the roles and underlying mechanisms of m6A modification is a hot topic in cancer research and recent investigations showed that m6A regulators can regulate various tumor progression processes [11, 54]. Based on the limitations of previous related studies our group’s previous studies reported that the key markers and TME of LGG and GBM are not the same. In previous studies, the research group has found that m6Ascore accurately predicted LGG progression. Thus, in this present study, we carried out a precise analysis of GBM and found that the m6Acore could predict GBM patients’ prognosis and immunotherapy effects. These findings could be used as references to further determine the GBM’s etiology and assess more beneficial treatments.

Conclusions

Our findings showed that m6A modification played pivotal roles in the tumorigenesis and TME infiltration characterization of GBM using data from large cohorts. Our proposed m6Ascore accurately predicted patients’ prognosis and potential chemotherapy and immunotherapy responses, providing novel perspectives to deepen our understanding on the pathogenesis of GBM and identify potential targets to improve treatment outcomes.

Abbreviations

AUC: area under the curve; CGGA: Chinese Glioma Genome Atlas; CNV: Copy number variation; DEGs: Differential expressed genes; EMT: epithelial-mesenchymal transitions; GGI: Gene-gene interaction; GSVA: Gene set variation analysis; GBM: glioblastoma; OS: overall survival; ROC: Receiver operating characteristic; ssGSEA: single sample Gene Set Enrichment Analysis; TCGA: the Cancer Genome Atlas; TIDE: tumor immune dysfunction and exclusion; TME: tumor microenvironment.

Author Contributions

The work presented here was carried out with the collaboration of all authors. WR Liu and CY Li defined the research topic and discussed the analyses and their interpretation, and presentation. WR Liu, HD Huang and HN Huang developed the algorithm, found clinicopathological data and performed the statistical analyses. WR Liu and CM Liu drafted the manuscript, recorded the clinical data, analysed the m6A events and interpreted the results. CY Li, XY W and CM L participated in reviewing all clinical records and performed the associated data collection. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Funding

This study was supported by grants from: 1. 2020 Guangxi Zhuang Autonomous Region Health Committee self-funded scientific research project, project number 20201558, 2. In 2020, the general project of high-level talent scientific research project of the Affiliated Hospital of Youjiang Medical College for Nationalities (the young and middle-aged backbone talent project), contract number Y202011702, 3. 2021 Guangxi University’s young and middle-aged teachers’ basic research ability improvement project, project number 2021KY0542.

Editorial Note

&

This corresponding author has a verified history of publications using a personal email address for correspondence.

References

  • 1. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2012-2016. Neuro Oncol. 2019 (Suppl 5); 21:v1–100. https://doi.org/10.1093/neuonc/noz150 [PubMed]
  • 2. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD, Kleihues P, Ellison DW. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol. 2016; 131:803–20. https://doi.org/10.1007/s00401-016-1545-1 [PubMed]
  • 3. Davis ME. Glioblastoma: Overview of Disease and Treatment. Clin J Oncol Nurs. 2016; 20:S2–8. https://doi.org/10.1188/16.CJON.S1.2-8 [PubMed]
  • 4. Johnson DR, Fogh SE, Giannini C, Kaufmann TJ, Raghunathan A, Theodosopoulos PV, Clarke JL. Case-Based Review: newly diagnosed glioblastoma. Neurooncol Pract. 2015; 2:106–21. https://doi.org/10.1093/nop/npv020 [PubMed]
  • 5. Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJ, Belanger K, Brandes AA, Marosi C, Bogdahn U, Curschmann J, Jan zer RC, Ludwin SK, et al, and European Organisation for Research and Treatment of Cancer Brain Tumor and Radiotherapy Groups, and National Cancer Institute of Canada Clinical Trials Group. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med. 2005; 352:987–96. https://doi.org/10.1056/NEJMoa043330 [PubMed]
  • 6. Alarcón CR, Lee H, Goodarzi H, Halberg N, Tavazoie SF. N6-methyladenosine marks primary microRNAs for processing. Nature. 2015; 519:482–5. https://doi.org/10.1038/nature14281 [PubMed]
  • 7. Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, Jaffrey SR. m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature. 2016; 537:369–73. https://doi.org/10.1038/nature19342 [PubMed]
  • 8. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017; 18:31–42. https://doi.org/10.1038/nrm.2016.132 [PubMed]
  • 9. Wang S, Sun C, Li J, Zhang E, Ma Z, Xu W, Li H, Qiu M, Xu Y, Xia W, Xu L, Yin R. Roles of RNA methylation by means of N6-methyladenosine (m6A) in human cancers. Cancer Lett. 2017; 408:112–20. https://doi.org/10.1016/j.canlet.2017.08.030 [PubMed]
  • 10. Wang CX, Cui GS, Liu X, Xu K, Wang M, Zhang XX, Jiang LY, Li A, Yang Y, Lai WY, Sun BF, Jiang GB, Wang HL, et al. METTL3-mediated m6A modification is required for cerebellar development. PLoS Biol. 2018; 16:e2004880. https://doi.org/10.1371/journal.pbio.2004880 [PubMed]
  • 11. Lu X, Li C, Xu W, Wu Y, Wang J, Chen S, Zhang H, Huang H, Huang H, Liu W. Malignant Tumor Purity Reveals the Driven and Prognostic Role of CD3E in Low-Grade Glioma Microenvironment. Front Oncol. 2021; 11:676124. https://doi.org/10.3389/fonc.2021.676124 [PubMed]
  • 12. Pinello N, Sun S, Wong JJ. Aberrant expression of enzymes regulating m6A mRNA methylation: implication in cancer. Cancer Biol Med. 2018; 15:323–34. https://doi.org/10.20892/j.issn.2095-3941.2018.0365 [PubMed]
  • 13. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, Sun G, Lu Z, Huang Y, Yang CG, Riggs AD, He C, Shi Y. m6A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cell Rep. 2017; 18:2622–34. https://doi.org/10.1016/j.celrep.2017.02.059 [PubMed]
  • 14. Granier C, Karaki S, Roussel H, Badoual C, Tran T, Anson M, Fabre E, Oudard S, Tartour E. Immunothérapie des cancers : rationnel et avancées récentes [Cancer immunotherapy: Rational and recent breakthroughs]. Rev Med Interne. 2016; 37:694–700. https://doi.org/10.1016/j.revmed.2016.05.023 [PubMed]
  • 15. Luo X, Li H, Liang J, Zhao Q, Xie Y, Ren J, Zuo Z. RMVar: an updated database of functional variants involved in RNA modifications. Nucleic Acids Res. 2021; 49:D1405–12. https://doi.org/10.1093/nar/gkaa811 [PubMed]
  • 16. Helmy KY, Patel SA, Nahas GR, Rameshwar P. Cancer immunotherapy: accomplishments to date and future promise. Ther Deliv. 2013; 4:1307–20. https://doi.org/10.4155/tde.13.88 [PubMed]
  • 17. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, Powderly JD, Carvajal RD, Sosman JA, Atkins MB, Leming PD, Spigel DR, Antonia SJ, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012; 366:2443–54. https://doi.org/10.1056/NEJMoa1200690 [PubMed]
  • 18. Chen K, Song B, Tang Y, Wei Z, Xu Q, Su J, de Magalhães JP, Rigden DJ, Meng J. RMDisease: a database of genetic variants that affect RNA modifications, with implications for epitranscriptome pathogenesis. Nucleic Acids Res. 2021; 49:D1396–404. https://doi.org/10.1093/nar/gkaa790 [PubMed]
  • 19. Zheng J, Wang X, Qiu Y, Wang M, Yu H, Zhou Z, Wu Z, Jiang X. Identification of Critical m6A RNA Methylation Regulators with Prognostic Value in Lower-Grade Glioma. Biomed Res Int. 2021; 2021:9959212. https://doi.org/10.1155/2021/9959212 [PubMed]
  • 20. Tu Z, Wu L, Wang P, Hu Q, Tao C, Li K, Huang K, Zhu X. N6-Methylandenosine-Related lncRNAs Are Potential Biomarkers for Predicting the Overall Survival of Lower-Grade Glioma Patients. Front Cell Dev Biol. 2020; 8:642. https://doi.org/10.3389/fcell.2020.00642 [PubMed]
  • 21. Du J, Ji H, Ma S, Jin J, Mi S, Hou K, Dong J, Wang F, Zhang C, Li Y, Hu S. m6A regulator-mediated methylation modification patterns and characteristics of immunity and stemness in low-grade glioma. Brief Bioinform. 2021; 22:bbab013. https://doi.org/10.1093/bib/bbab013 [PubMed]
  • 22. Wood SL, Pernemalm M, Crosbie PA, Whetton AD. The role of the tumor-microenvironment in lung cancer-metastasis and its relationship to potential therapeutic targets. Cancer Treat Rev. 2014; 40:558–66. https://doi.org/10.1016/j.ctrv.2013.10.001 [PubMed]
  • 23. Wang Q, Li P, Wu W. A systematic analysis of immune genes and overall survival in cancer patients. BMC Cancer. 2019; 19:1225. https://doi.org/10.1186/s12885-019-6414-6 [PubMed]
  • 24. Quan C, Belaydi O, Hu J, Li H, Yu A, Liu P, Yi Z, Qiu D, Ren W, Ma H, Gong G, Ou Z, Chen M, et al. N6-Methyladenosine in Cancer Immunotherapy: An Undervalued Therapeutic Target. Front Immunol. 2021; 12:697026. https://doi.org/10.3389/fimmu.2021.697026 [PubMed]
  • 25. Zhong J, Liu Z, Cai C, Duan X, Deng T, Zeng G. m6A modification patterns and tumor immune landscape in clear cell renal carcinoma. J Immunother Cancer. 2021; 9:e001646. https://doi.org/10.1136/jitc-2020-001646 [PubMed]
  • 26. Wu X, Sheng H, Wang L, Xia P, Wang Y, Yu L, Lv W, Hu J. A five-m6A regulatory gene signature is a prognostic biomarker in lung adenocarcinoma patients. Aging (Albany NY). 2021; 13:10034–57. https://doi.org/10.18632/aging.202761 [PubMed]
  • 27. Mo P, Xie S, Cai W, Ruan J, Du Q, Ye J, Mao J. N6-methyladenosine (m6A) RNA methylation signature as a predictor of stomach adenocarcinoma outcomes and its association with immune checkpoint molecules. J Int Med Res. 2020; 48:300060520951405. https://doi.org/10.1177/0300060520951405 [PubMed]
  • 28. Zhao H, Xu Y, Xie Y, Zhang L, Gao M, Li S, Wang F. m6A Regulators Is Differently Expressed and Correlated With Immune Response of Esophageal Cancer. Front Cell Dev Biol. 2021; 9:650023. https://doi.org/10.3389/fcell.2021.650023 [PubMed]
  • 29. Li J, Wang W, Zhou Y, Liu L, Zhang G, Guan K, Cui X, Liu X, Huang M, Cui G, Sun R. m6A Regulator-Associated Modification Patterns and Immune Infiltration of the Tumor Microenvironment in Hepatocarcinoma. Front Cell Dev Biol. 2021; 9:687756. https://doi.org/10.3389/fcell.2021.687756 [PubMed]
  • 30. Zhao S, Zhang Y, Lu X, Ding H, Han B, Song X, Miao H, Cui X, Wei S, Liu W, Chen S, Wang J. CDC20 regulates the cell proliferation and radiosensitivity of P53 mutant HCC cells through the Bcl-2/Bax pathway. Int J Biol Sci. 2021; 17:3608–21. https://doi.org/10.7150/ijbs.64003 [PubMed]
  • 31. Liu W, Xu W, Chen Y, Gu L, Sun X, Qu Y, Zhang H, Liu X, Huang H. Elevated double-strand break repair protein RAD50 predicts poor prognosis in hepatitis B virus-related hepatocellular carcinoma: A study based on Chinese high-risk cohorts. J Cancer. 2020; 11:5941–52. https://doi.org/10.7150/jca.46703 [PubMed]
  • 32. Liu WR, Li CY, Xu WH, Liu XJ, Tang HD, Huang HN. Genome-wide analyses of the prognosis-related mRNA alternative splicing landscape and novel splicing factors based on large-scale low grade glioma cohort. Aging (Albany NY). 2020; 12:13684–700. https://doi.org/10.18632/aging.103491 [PubMed]
  • 33. Feng L, Qian H, Yu X, Liu K, Xiao T, Zhang C, Kuang M, Cheng S, Li X, Wan J, Zhang K. Heterogeneity of tumor-infiltrating lymphocytes ascribed to local immune status rather than neoantigens by multi-omics analysis of glioblastoma multiforme. Sci Rep. 2017; 7:6968. https://doi.org/10.1038/s41598-017-05538-z [PubMed]
  • 34. Yu X, Feng L, Liu D, Zhang L, Wu B, Jiang W, Han Z, Cheng S. Quantitative proteomics reveals the novel co-expression signatures in early brain development for prognosis of glioblastoma multiforme. Oncotarget. 2016; 7:14161–71. https://doi.org/10.18632/oncotarget.7416 [PubMed]
  • 35. Kawaguchi A, Yajima N, Tsuchiya N, Homma J, Sano M, Natsumeda M, Takahashi H, Fujii Y, Kakuma T, Yamanaka R. Gene expression signature-based prognostic risk score in patients with glioblastoma. Cancer Sci. 2013; 104:1205–10. https://doi.org/10.1111/cas.12214 [PubMed]
  • 36. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016; 44:e71. https://doi.org/10.1093/nar/gkv1507 [PubMed]
  • 37. Curtis AE, Smith TA, Ziganshin BA, Elefteriades JA. The Mystery of the Z-Score. Aorta (Stamford). 2016; 4:124–30. https://doi.org/10.12945/j.aorta.2016.16.014 [PubMed]
  • 38. Hart T, Komori HK, LaMere S, Podshivalova K, Salomon DR. Finding the active genes in deep RNA-seq gene expression studies. BMC Genomics. 2013; 14:778. https://doi.org/10.1186/1471-2164-14-778 [PubMed]
  • 39. Hancks DC, Kazazian HH Jr. SVA retrotransposons: Evolution and genetic instability. Semin Cancer Biol. 2010; 20:234–45. https://doi.org/10.1016/j.semcancer.2010.04.001 [PubMed]
  • 40. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
  • 41. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011; 27:1739–40. https://doi.org/10.1093/bioinformatics/btr260 [PubMed]
  • 42. Xiao B, Liu L, Li A, Xiang C, Wang P, Li H, Xiao T. Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front Oncol. 2020; 10:607622. https://doi.org/10.3389/fonc.2020.607622 [PubMed]
  • 43. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 44. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, Nordgren H, Farmer P, Praz V, Haibe-Kains B, Desmedt C, Larsimont D, Cardoso F, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006; 98:262–72. https://doi.org/10.1093/jnci/djj052 [PubMed]
  • 45. Sun L, Ke X, Wang D, Yin H, Jin B, Xu H, Du S, Xu Y, Zhao H, Lu X, Sang X, Zhong S, Yang H, Mao Y. Prognostic Value of the Albumin-to-γ-glutamyltransferase Ratio for Gallbladder Cancer Patients and Establishing a Nomogram for Overall Survival. J Cancer. 2021; 12:4172–82. https://doi.org/10.7150/jca.49242 [PubMed]
  • 46. Lawrence MS, Stojanov P, Polak P, Kryukov GV, Cibulskis K, Sivachenko A, Carter SL, Stewart C, Mermel CH, Roberts SA, Kiezun A, Hammerman PS, McKenna A, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013; 499:214–8. https://doi.org/10.1038/nature12213 [PubMed]
  • 47. Karachi A, Dastmalchi F, Mitchell DA, Rahman M. Temozolomide for immunomodulation in the treatment of glioblastoma. Neuro Oncol. 2018; 20:1566–72. https://doi.org/10.1093/neuonc/noy072 [PubMed]
  • 48. Sengupta S, Marrinan J, Frishman C, Sampath P. Impact of temozolomide on immune response during malignant glioma chemotherapy. Clin Dev Immunol. 2012; 2012:831090. https://doi.org/10.1155/2012/831090 [PubMed]
  • 49. Duffau H. Lessons from brain mapping in surgery for low-grade glioma: insights into associations between tumour and brain plasticity. Lancet Neurol. 2005; 4:476–86. https://doi.org/10.1016/S1474-4422(05)70140-X [PubMed]
  • 50. Liu W, Xu W, Li C, Xu J, Huang K, Hu R, Huang H, Liu X. Network pharmacological systems study of Huang-Lian-Tang in the treatment of glioblastoma multiforme. Oncol Lett. 2021; 21:18. https://doi.org/10.3892/ol.2020.12279 [PubMed]
  • 51. Zhang C, Cheng W, Ren X, Wang Z, Liu X, Li G, Han S, Jiang T, Wu A. Tumor Purity as an Underlying Key Factor in Glioma. Clin Cancer Res. 2017; 23:6279–91. https://doi.org/10.1158/1078-0432.CCR-16-2598 [PubMed]
  • 52. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
  • 53. Rhee JK, Jung YC, Kim KR, Yoo J, Kim J, Lee YJ, Ko YH, Lee HH, Cho BC, Kim TM. Impact of Tumor Purity on Immune Gene Expression and Clustering Analyses across Multiple Cancer Types. Cancer Immunol Res. 2018; 6:87–97. https://doi.org/10.1158/2326-6066.CIR-17-0201 [PubMed]
  • 54. Liu W, Li C, Wu Y, Xu W, Chen S, Zhang H, Huang H, Zhao S, Wang J. Integrating m6A Regulators-Mediated Methylation Modification Models and Tumor Immune Microenvironment Characterization in Caucasian and Chinese Low-Grade Gliomas. Front Cell Dev Biol. 2021; 9:725764. https://doi.org/10.3389/fcell.2021.725764 [PubMed]
Download PDF
Corresponding Authors

Huadong Huang
dongdongh@126.com
https://orcid.org/0009-0004-8490-7669

Haineng Huang
passpsycho@163.com
https://orcid.org/0000-0002-4320-4623

Xiangyu Wang
wang_xy123@126.com
https://orcid.org/0000-0001-5103-433X

Keywords
m6A prognosis consistent clustering metastasis glioblastoma
Paper Sections