Prognostic value of m6A regulators and the nomogram construction in glioma patients.

Although N6-methyladenosine (m6A) has been implicated in various biological functions in human cancers, its role in predicting the prognosis of glioma remains unclear. In this study, the transcriptome expression profiles and the clinical data of 961 patients were derived from the Chinese Glioma Genome Atlas (CGGA). We comprehensively evaluated the association between the expression of m6A regulators and the prognosis of glioma and established a 3-gene (YTHDF2, FTO, and ALKBH5) risk signature using least absolute shrinkage and selection operator (LASSO) analysis. Patients with a high-risk signature had significantly adverse prognoses. Gene set enrichment analysis (GSEA) analysis revealed that the G2M checkpoint, MTORC1 signaling, epithelial mesenchymal transition, and PI3K-AKT-mTOR signaling were significantly enriched in the high-risk group. Univariate and multivariate Cox regression analyses confirmed the independent prognostic value of this risk signature. We then constructed a nomogram for individualized prediction of overall survival (OS) by integrating clinicopathological features (age, World Health Organization [WHO] grade), treatment information (radiotherapy, temozolomide therapy), and m6A risk signature. The calibration curves showed excellent agreement between the predicted and actual probabilities for the 1-, 3-, and 5-year OS, with a C-index of 0.780 in the training cohort and 0.717 in the validation cohort. Altogether, our study elucidated the important role of m6A regulators in glioma prognosis, which is valuable for the selection of therapeutic methods and clinical management of patients with glioma.

1. Introduction

Gliomas are the most common and lethal tumors of the central nervous system (CNS), representing approximately 25.1% of all primary CNS tumors and 80.8% of malignant tumors according to the 2020 CBTRUS report.[ Assessment of prognosis in patients with glioma is complicated and affected by multiple factors, including clinical features, treatment strategies, histological types, and molecular alterations. The World Health Organization (WHO) classification of gliomas (2007) primarily relies on tumor histology,[ which results in substantial variation in patient survival within the same WHO grade. With the advancement of molecular research in gliomas, some significant biomarkers have been discovered, such as the isocitrate dehydrogenase mutation, 1p/19q codeletion, and MGMT methylation.[ In 2016, the WHO revised the classification system of gliomas by integrating morphology with molecular alterations,[ which deepened our insight into glioma progression and prognosis. Even though great progress has been made in the prognosis assessment of gliomas, there remains much room for improvement. Therefore, more effective biomarkers and prognostic prediction models are urgently required.

N6-methyladenosine (m6A), which has been implicated in various biological functions, is 1 of the most common forms of ribonucleic acid (RNA) modification in eukaryotic mRNAs. The extent of m6A methylation is mediated by multiple m6A regulators, including methyltransferases (“writers”), demethylases (“erasers”), and binding proteins (“readers”).[ “Writers” are important catalytic enzymes that stimulate methylation at the N6-position of adenosine.[ These proteins do not exist in isolation, but form complexes to execute their biological functions. “Erasers” are responsible for removing m6A from RNA.[ “Readers” recognize and bind to m6A sites, resulting in different destinies of target mRNA.[ Accumulating evidence shows that m6A plays a dual role in cancer pathogenesis and progression.[ Unlike miRNA regulation, which usually causes gene silencing, m6A can promote or inhibit the expression of tumor-related genes at the post-transcriptional regulatory level, thereby acting as an oncogene or tumor suppressor role.[ Although the regulatory mechanisms of m6A have been revealed, its role, especially in the prognostic prediction of glioma, remains unclear.

In this study, we explored the potential application of m6A regulators in predicting the OS of glioma patients. A risk signature based on the key prognostic m6A regulators was established using LASSO analysis. Subgroup and Cox regression analyses were performed to confirm the independent prognostic value of the risk signature. Finally, a nomogram was built for the individualized prediction of OS in patients with glioma.

2. Material and Methods

2.1. Data extraction

The transcriptome expression profiles (fragments per kilobase of transcript per million mapped reads normalized) and the corresponding clinicopathological data were downloaded from the Chinese Glioma Genome Atlas (CGGA). Patients with missing WHO grades, missing OS values, or OS < 30 days were excluded to reduce statistical bias. Finally, 305 and 656 patients from the CGGA _325 and CGGA _693 datasets were defined as the training and validation cohorts, respectively. The patient characteristics in the 2 cohorts are summarized in Table 1. The CGGA database is free of patient identification information. Ethics committee approval was not applicable to this study.

Table 1

Clinicopathological characteristics of glioma patients from CGGA database.

	Training cohort	Validation cohort
No. of patients	305	656
Age
Mean (SD)	43.3 (12.0)	43.4 (12.4)
Gender (%)
Female	115 (37.7)	282 (43.0)
Male	190 (62.3)	374 (57.0)
Grade (%)
WHO II	97 (31.8)	172 (26.2)
WHO III	73 (23.9)	248 (37.8)
WHO IV	135 (44.3)	236 (36.0)
Radiotherapy (%)
Untreated	60 (19.7)	130 (19.8)
Treated	236 (77.4)	501 (76.4)
NA	9 (3.0)	25 (3.8)
TMZ therapy (%)
Untreated	107 (35.1)	155 (23.6)
Treated	186 (61.0)	480 (73.2)
NA	12 (3.9)	21 (3.2)
OS [yr (%)]
≤1	93 (30.5)	172 (26.2)
1~3	85 (27.9)	205 (31.2)
3~5	31 (10.2)	104 (15.9)
>5	96 (31.5)	175 (26.7)

CGGA = Chinese Glioma Genome Atlas, OS = overall survival, TMZ = temozolomide; WHO = World Health Organization.

2.2. M6a regulators

Based on previously published literature,[ a total of 20 m6A regulators containing 8 writers (METTL3, METTL14, METTL16, WTAP, KIAA1429, ZC3H13, RBM15, and RBM15B), 10 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, and RBMX), and 2 erasers (FTO and ALKBH5) were collected in this study.

2.3. Consensus clustering

Consensus clustering was performed using the “ConsensusClusterPlus” R package to categorize patients with gliomas into different prognostic subgroups (k). Euclidean distance, calculated by the expression of m6A regulators, was used to measure sample similarity.

2.4. LASSO analysis

The prognostic m6A risk signature was constructed using least absolute shrinkage and selection operator (LASSO) analysis in the training cohort. A 10-fold cross-validation for tuning parameter selection was applied to the LASSO model to avoid overfitting. The largest lambda, whose error was within 1 standard error of the optimum (lambda.1se), was selected to obtain parsimonious models as recommended.[ The coefficients obtained from the LASSO regression algorithm were used to yield the following risk-score equation: . According to this equation, the risk score for each patient was calculated separately for the training and validation cohorts. Subsequently, the patients were divided into high-risk and low-risk signature groups, and the median risk score was set as the cutoff point.

2.5. GSEA analysis

To explore the potential molecular mechanisms between high-risk and low-risk signature groups, we performed gene set enrichment analysis (GSEA) to identify enriched terms. Hallmark gene sets were chosen from the Molecular Signature Database as the reference gene sets. Normalized enrichment score (|NES|) > 1, nominal P value < .05 and FDR < 0.25 were set as the cutoff criteria.

2.6. Independent prognostic analysis

Univariate and multivariate Cox regression analyses were conducted to determine the independent prognostic value of the m6A risk signature and other clinical features, including age, sex, grade, radiotherapy, and TMZ therapy. Backward stepwise regression was used in multivariate analysis to obtain the best-fit model with the lowest Akaike information criterion.[ Hazard ratios (HR) with 95% confidence intervals (CI) were calculated to estimate the risk of death. Variables with P values < .05 both in the univariate and multivariate analyses were considered independent prognostic variables.

2.7. Construction and validation of nomogram

Variables in the best-fit model were included in nomogram construction. The nomogram for predicting OS was established by using the “survival” and “rms” package in R (version 4.0.3).

The nomogram was subjected to 1000 bootstrap resamples for internal validation in the training cohort and external validation in the validation cohort. The predictive performance was assessed by calculating the concordance index (C-index), which ranges from 0.5 1.0. A C-index that is closer to 1.0 represents more accurate prognostic predictions. Calibration curves for the 1-, 3-, and 5-year OS were obtained by comparing the predicted probabilities from the nomogram with the actual probabilities. A perfect calibration curve falls along the 45-degree line.

2.8. Statistical analysis

Statistical analysis and visualization were performed using R version 4.0.3. The Kruskal–Wallis test was used to compare mRNA expression in the WHO grade subgroups. Subgroup analysis was performed to evaluate the stability of the risk signature. Kaplan–Meier analysis and log-rank tests were used to construct survival curves and determine the significance between different groups. P < .05 was used as the significant threshold.

3. Results

3.1. Expression of m6a regulators in glioma patients

To explore the biological role of m6A regulators in gliomas, we systematically compared the differences in gene expression between WHO grade subgroups in the training cohort. Fourteen of the 20 m6A regulators showed significant differences in gene expression (Fig. 1, P < .05). The results of the corresponding Kaplan–Meier survival analysis stratified by WHO grade are presented in Figure S1, Supplemental Digital Content 1, http://links.lww.com/MD/H347. This finding indicates that m6A regulators have the potential to affect the survival outcomes of glioma patients.

Figure 1.

Expression of m6A regulators in gliomas with different WHO grades. m6A = N6-methyladenosine, WHO = World Health Organization.

3.2. Different expression patterns of m6a regulators

To further investigate the influence of m6A regulators on the survival prognosis of patients with glioma, we performed consensus clustering based on the expression levels of m6A regulators in the training cohort. A heatmap of the consensus matrix for each k and the cumulative distribution function were used to assess the cluster stability and determine the optimal k. As shown in Figure S2, Supplemental Digital Content 2, http://links.lww.com/MD/H348 k = 2 was identified with optimal clustering stability from k = 2 to 9. Patients were clustered into 2 subgroups: cluster1 and cluster2. Kaplan–Meier analysis revealed that patients in cluster1 had a markedly worse overall survival (OS) than that of those in cluster2 (Fig. 2). This result revealed a close connection between the different expression patterns of m6A regulators and survival prognosis in glioma patients.

Figure 2.

Survival outcome in glioma patients with different expression patterns of m6A regulators. (A) Heatmap of the consensus matrix for k = 2. (B) Kaplan–Meier curves of overall survival in different clusters. m6A = N6-methyladenosine.

3.3. Identification of prognosis-related genes for m6a regulators

Considering the multicollinearity between m6A regulators, we conducted LASSO analysis to screen the key prognosis-related m6A regulators. Three m6A regulators were identified: YTHDF2, FTO, and ALKBH5 (Fig. 3A and B). The risk score based on YTHDF2, FTO, and ALKBH5 was calculated using the formula: Risk score = 0.940 × YTHDF2 expression - 0.255 × FTO expression + 0.259 × ALKBH5 expression. The patients were divided into high-risk and low-risk signature groups according to the median value of the risk score. Patients with a high-risk signature had a significantly adverse prognosis in both the training and validation cohorts (Fig. 3C and D). Multivariate Cox regression was performed on YTHDF2, FTO, and ALKBH5.The high expression of YTHDF2 and ALKBH5 (HR > 1, P < .001) and low expression of FTO (HR < 1, P < .001) were closely related to poor prognosis in patients with glioma (Fig. 3E). These results substantiate the reliability of YTHDF2, FTO, and ALKBH5 in glioma prognosis prediction.

Figure 3.

(A) LASSO coefficients of the 20 m6A regulators. (B) Partial likelihood deviance in the LASSO model. (C and D) Kaplan–Meier curves for OS in training cohort and validation cohort between high-risk and low-risk signature. (E) Multivariate Cox regression of 3 key prognostic m6A regulators. LASSO = least absolute shrinkage and selection operator, m6A = N6-methyladenosine, OS = overall survival.

3.4. Regulatory mechanisms in patients with different risk signature

To elucidate the regulatory mechanisms contributing to the differences in prognosis among patients with different risk signatures, we performed GSEA analysis. The results revealed that gene sets including G2M checkpoint (NES = 1.851, P = .009, FDR = 0.028), MTORC1 signaling (NES = 1.842, P = .004, FDR = 0.027), epithelial mesenchymal transition (NES = 1.720, P = .03, FDR = 0.054), and PI3K-AKT-mTOR signaling (NES = 1.707, P = .009, FDR = 0.054), were upregulated and significantly enriched in patients with high-risk signatures (Fig. 4). The activation of these cancer-promoting signaling pathways might be the reason why the high-risk signature based on YTHDF2, FTO, and ALKBH5 portends a poor survival prognosis.

Figure 4.

GSEA analysis among patients with different risk signature. GSEA = gene set enrichment analysis, NES = normalized enrichment score.

3.5. Independent prognostic analysis for risk signature

To confirm the independent prognostic value of the risk signature based on YTHDF2, FTO, and ALKBH5 in OS, we first conducted a subgroup Cox regression analysis. The patients were stratified by age, sex, grade, radiotherapy, and Temozolomide (TMZ) therapy status. As shown in Figure 5A, we found that the high-risk signature was a risk factor in all other subgroups (HR > 1, P < .001), except for the WHO II subgroup. For patients in WHO II, HR > 1 indicates that a high-risk signature is a risk factor, but there is no statistical difference in the P value. This phenomenon may be related to the disproportion in the number of patients with high- and low-risk signatures (Fig. 5B). In the WHO II subgroup, 92.78% of the population had a low-risk signature, and only 7.22% of the population had high-risk signatures. Therefore, the prognostic effect of the risk signature in the WHO II population requires further study with a larger sample size.

Figure 5.

(A) Subgroup univariate analysis for risk signature. (B) The distribution of patients with high-risk and low-risk signature in the different WHO subgroups. WHO = World Health Organization.

The results of univariate analysis for each variable are shown in Figure 6A. Sex was excluded from the final multivariate analysis after backward stepwise regression (Fig. 6B). Risk signature, WHO grade, radiotherapy, and TMZ therapy were significantly associated with the OS in both the univariate and multivariate analyses. Regarding age, there was a significant difference (P < .001) in univariate analysis, but no significant difference (P = .10) in multivariate analysis, which could be attributed to the correlation between age and risk signature. As shown in Figure 6C, patients in the high-risk group were significantly older than those in the low-risk group (P = .001). Although age lost its significance in multivariate analysis, the prediction performance and fitting ability of the model were weakened (Akaike information criterion value would increase) after excluding age. Given the borderline significance of age and its clinical relevance for predicting OS,[ age was also incorporated in the subsequent construction of the nomogram. TMZ is an important treatment option for glioma. Clinically, patients with higher WHO grades and more risk factors are more likely to be recommended TMZ therapy. As shown in Figure 6D, only 38.04% of WHO II patients received TMZ therapy, whereas this proportion was 75.71% and 74.81% in WHO III and IV patients, respectively. This phenomenon resulted in the HR of TMZ therapy being >1 in univariate analysis. After adjusting for other prognostic factors in the multivariate analysis, the HR value was <1, indicating that receiving TMZ chemotherapy was a protective factor. Univariate analysis only included TMZ therapy status but did not consider the influence of other variables, which masked the protective effect of TMZ.

Figure 6.

Independent prognostic analysis for risk signature. (A and B) Univariate and multivariate Cox regression analyses of the association between clinical features and OS. (C) The correlation between age and risk signature. (D) The distribution of TMZ therapy status in the different WHO subgroups. OS = overall survival, TMZ = temozolomide, WHO = World Health Organization.

Altogether, these results confirmed that the risk signature based on YTHDF2, FTO, and ALKBH5 was an independent prognostic factor for patients with glioma.

3.6. Construction and validation of nomogram

Age, WHO grade, radiotherapy, TMZ therapy, and risk signature were incorporated into the construction of a nomogram to predict the 1-, 3-, and 5-year OS probability of glioma patients. The corresponding value for each factor was assigned a score on a point scale. The total score was calculated by adding the scores for each factor. We can easily obtain the estimated probability of the OS by projecting the total score to the bottom scale (Fig. 7A).

Figure 7.

Construction and validation of nomogram. (A) The prognostic nomogram for predicting the 1-, 3-, and 5-year OS of glioma patients. (B and C) The calibration curves for predictions of OS in training cohort and validation cohort. OS = overall survival.

The predictive performance of the nomogram was evaluated by calculating its C-index. The C-index values in the training and validation cohorts were 0.780 (95% CI: 0.754–0.806) and 0.717 (95% CI: 0.692–0.742), respectively. The calibration curves showed excellent agreement between the predicted and actual probabilities in the training and validation cohorts for 1-, 3-, and 5-year OS (Fig. 7B and C).

In summary, the nomogram constructed in this study has good predictive value and provides clinicians with a more accurate and practical prediction tool, which is conducive to individualized and accurate assessment of OS in glioma patients.

4. Discussion

Accurate prognostic evaluation is of great significance for the management and follow-up of patients with glioma. Some classical molecular markers are closely associated with prognosis, such as isocitrate dehydrogenase mutations, 1p/19q co-deletion, EGFR, and others.[ With advances in epigenetic research, some potential novel prognostic markers have emerged in recent years, including m6A. In this study, we demonstrated the independent prognostic value of m6A regulators for survival outcomes in glioma patients and proposed possible regulatory mechanisms. Moreover, an m6A regulator-based nomogram was constructed, providing a visualization of the predicted results for 1-, 3-, and 5-year OS of glioma patients.

Initially, we comprehensively analyzed the expression of 20 m6A regulators in the different WHO grade subgroups. Multiple m6A regulators, such as ALKBH5, KIAA1429, RBM15, YTHDF1, and YTHDF2, showed a remarkably incremental trend with the progression of glioma. However, the expression of FTO gradually decreased with the progression of glioma. These findings demonstrate that m6A regulators have the potential to influence the prognosis of glioma patients, which was further proven by the results of consensus clustering. Three key prognosis-related m6A regulators were identified by LASSO analysis, including YTHDF2, FTO, and ALKBH5. As an m6A reader, YTHDF2 plays a critical role in regulating mRNA stability by inducing degradation of transcripts via interaction with the CCR4-NOT and ribonuclease P/MRP complexes.[ YTHDF2 mediates mRNA degradation of the tumor suppressors LHPP and NKX3-1 in prostate cancer, which results in the activation of PI3K-AKT-mTOR signaling and promotes cell proliferation and migration.[ Notably, YTHDF2 can promote the cancer stem cell enrichment and enhance tumor growth in liver cancer and glioblastoma.[ ALKBH5, an m6A demethylase, belongs to the AlkB family of non-heme Fe(II)/α-ketoglutarate-dependent dioxygenases, but its role in human cancers is controversial.[ ALKBH5 serves as a tumor suppressor in pancreatic cancer and hepatocellular carcinoma by inhibiting the proliferation and invasion of cancer cells.[ In our study, the expression of ALKBH5 was significantly upregulated in WHO grades III and IV, indicating an adverse prognosis and oncogenic role in glioma. ALKBH5 can induce epithelial-to-mesenchymal transition, promoting the migration and invasion of uveal melanoma cells.[ In addition, ALKBH5 can maintain the tumorigenicity of glioblastoma stem cells and promote their proliferation and invasion capability.[ Moreover, Kowalski-Chauvel et al reported that ALKBH5 promotes the radio resistance of glioblastoma stem cells by controlling homologous repair.[ FTO, which encodes an alpha-ketoglutarate-dependent dioxygenase, was the first demethylase to remove m6A.[ FTO knockout in GC-1 cells triggers G2/M and cell cycle arrest.[ Prolonged G2M arrest can inhibit cell proliferation and induce apoptosis. FTO has been reported to inhibit aggression in gliomas by interacting with FOXO3a.[ FTO has also been reported to function as a tumor suppressor in ovarian and thyroid cancer.[ Although some of the findings mentioned above have been reported in other types of cancer, these studies still provide support for our results of prognosis analysis among YTHDF2, FTO, and ALKBH5.

A nomogram is an effective tool for graphic representation of the model fit and has been widely used in clinical study.[ Several studies have identified gene signatures associated with glioma prognosis. Wang et al identified 7 risk genes and established a nomogram to predict the survival of patients with low-grade glioma.[ Liu et al proposed a prognostic model based on 5 risk pseudogenes.[ Compared with previous studies, the nomogram in this study showed the following strengths. First, we incorporated treatment information into the prognosis prediction, including radiotherapy and TMZ therapy status. Second, we identified 3 key prognostic genes, which make it easier for clinicians to evaluate gene expression compared with those of over 5 genes. Third, the nomogram showed satisfactory performance in predicting the OS of patients with glioma, and the prediction reliability was further strengthened in the validation cohort. Our study has several limitations. First, we did not conduct a prospective multicenter clinical trial to validate the m6A risk signature. Second, experimental studies should be performed to elucidate the potential molecular mechanisms. Third, surgical treatment is an important factor that affects the prognosis of patients with glioma. As there was no specific information on surgical treatment in the CGGA database, we did not incorporate surgical information into the nomogram.

In conclusion, the current study systematically analyzed the expression of m6A regulators and established a 3-gene risk signature. Clinicopathological features, treatment information, and m6A risk signature were integrated into the prognostic nomogram of glioma patients, which is valuable for individualized prognostic assessment and clinical management of patients with glioma. Our study provides new insights into the mechanisms of m6A regulation and therapeutic targets in glioma. However, further experimental and clinical studies are still required.

Acknowledgments

The authors would like to thank the CGGA database for the data availability.

Author contributions

Conceptualization: Pengdi Liu, Jianqiang Qu.

Data curation: Chengwen Ma.

Methodology: Pengdi Liu, Xianxia Yan.

Supervision: Jianqiang Qu.

Validation: Junxiang Gu, Fuyu Tian.

Visualization: Pengdi Liu, Chengwen Ma.

Writing – original draft: Pengdi Liu, Chengwen Ma.

Writing – review & editing: Jianqiang Qu, Xianxia Yan.