Integrative analyses reveal biological function and prognostic role of m7G methylation regulators in high-grade glioma

Patients and samples

The training and validation datasets were provided by the Chinese Glioma Genome Atlas (CGGA, cgga.org.cn, The CGGA database contains clinical and sequencing data of over 2,000 brain tumor samples from Chinese cohorts, and is equipped with a user-friendly web application for data storage and exploration.) and the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). Prior to performing further analyses, the normalization of both datasets took place, and samples lacking follow-up or clinical data were eliminated. Gliomas classified as WHO III and IV were of high grade [16]. After conducting a thorough analysis of the literature, we discovered 29 genes responsible for the regulation of m7G RNA methylation (Supplementary Table 1). As the information was obtained from publicly accessible databases, there was no need for ethical approval.

Bioinformatic analyses

GSEA and GSVA analyses

GSEA v4.2.3 was utilized to perform the gene set enrichment analysis (GSEA) [19]. Gene-set variation analysis (GSVA) was conducted using the R package called ‘GSVA’ [20].

Development and verification of a predictive model

To build and validate a model, as well as determine the final number of genes, we employed the LASSO regression, which is known as the least absolute shrinkage and selection operator [21]. Multivariate Cox regression modeling was subsequently conducted. The risk score for each sample was calculated using the following formula: coef1 multiplied by geneExpression1, plus coef (2) multiplied by gene Expression (2), and so on, up to coef (n) multiplied by gene Expression [22, 23]. Based on the median risk score value, the glioma patients were categorized into either a high-risk or low-risk group. Survival curves of these groups were compared using Kaplan-Meier (KM) analyses, and their ability to predict was assessed using a receiver operating characteristic curve (ROC). To establish a risk classification, the utilization of principal component analysis (PCA) was implemented. An individual risk assessment nomograph was generated. Furthermore, a comparison was made between the high-risk and low-risk groups in terms of the tumor microenvironment and levels of immune infiltration. The correlation between signature genes and drug IC50 values was evaluated.

Statistical analysis

Differences in categorical variables were examined using a Chi-squared test, whereas the t-test was employed to compare means. Multiple groups were compared using one-way ANOVA, and Dunnett's method was used for making multiple comparisons. To examine the associations between two continuous variables, a Spearman correlation analysis was conducted. The hazard ratio (HR) and its 95% confidence intervals (CIs) were calculated using both univariate and multivariate Cox regression. Statistical analyses were conducted using R version 4.0.1, and statistical significance was reported with a significance level of p < 0.05.

Data availability statement

The public dataset used in this study can be found in CGGA (http://www.cgga.org.cn/) and TCGA (https://portal.gdc.cancer.gov/).

Characteristics at the molecular level of regulators of m7G in high-grade glioma

The entire study flow was depicted in Figure 1. High-grade glioma exhibited high expression levels of twelve genes involved in m7G regulation (METTL1, WDR4, NSUN2, DCPS, NUDT4B, CYFIP1, NCBP1, NCBP2, EIF3D, EIF4A1, LSM1, and SNUPN), while eight genes (NUDT10, NUDT11, NUDT3, EIF4E1B, EIF4E3, LARP1, EIF4G3, and NCBP2L) showed low expression levels. The remaining genes did not display significant differences between tumor and normal samples (Figure 2A, 2B). The levels of gene alterations in m7G genes were below 1%, encompassing CYFIP1, GEMIN5, EIF4G3, IFIT5, NSUN2, AGO2, and NCBP2 (as shown in Figure 2C). The frequency of gain was higher than loss for eight genes, while the frequency of gain for 21 genes was lower (Figure 2D). The distribution of these genes occurred on chromosomes 2, 3, 5, 6, 8, 9, 10, 11, 12, 15, 17, 21, 22, and X can be observed in Figure 2E.

Figure 1. The flow chart of the study.

Figure 2. The landscape of m7G regulators in glioma. (A) Expression levels of m7G regulators between tumor and normal samples. (B) PCA of m7G regulators between tumor and normal samples. (C) The gene alterations of m7G regulators in advanced glioma. (D, E) Copy number variation frequencies and chromosomal location of m7G regulators in high-grade glioma.

m7G clusters and their characteristics in high-grade glioma

Using 29 m7G genes, we conducted a consensus matrix analysis and acquired three m7G clusters (Cluster 1, Cluster 2, and Cluster 3, as shown in Figure 3A). According to the individual political action committee (PAC), three evident distributions were observed (as shown in Figure 3B). Additionally, the overall survival (OS) displayed notable disparities, as depicted by the Kaplan-Meier curve (Figure 3C). According to the GSVA results, cluster 1 showed significant enrichment in oxidative phosphorylation, arachidonic acid metabolism, xenobiotic metabolism by cytochrome P450, and cytochrome P450 drug metabolism. Cluster 2, on the other hand, exhibited enrichment in ubiquitin-mediated proteolysis, progesterone-mediated oocyte maturation, oocyte meiosis, adherents’ junction, ERBB signaling pathway, and TGF beta signaling pathway. Lastly, cluster 3 demonstrated enrichment in lysosome, GPI anchor biosynthesis, protein export, pyrimidine metabolism, spliceosome, lysosome, and glycolate and dicarboxylate metabolism (Figure 3D–3F).

Figure 3. Clustering of advanced glioma based on m7G regulators. (A) Consensus clustering identified three subtypes. (B) PCA showed three distinct clusters. (C) Kaplan-Meier analysis indicated significant differences in overall survival among three clusters. (D–F) GSVA indicated the enriched pathways in three m7G clusters. (G–J) The differences in ESTIMAT score, stromal score, immune score, and tumor purity across three m7G clusters. (K, L) Comparisons of immune cell infiltrations level function among three m7G clusters. (M) Comparisons of Immune checkpoint-related genes expression among three m7G clusters.

Immune status in three m7G clusters

Initially, we assessed the tumor microenvironment in three m7G clusters. Significant variations were observed in the ESTIMATE, stromal, immune scores, and tumor purity. Cluster 1 exhibited elevated ESTIMATE, stromal, and immune scores in comparison to cluster 2 (Figure 3G–3I). Additionally, a similar trend was observed in tumor purity (Figure 3J). Moreover, notable variations in immune function exist among the three clusters. Cluster A exhibited elevated levels of APC co-stimulation, CCR, check-point, cytolytic activity, inflammation-promoting, T cell co-stimulation and inhibition, as well as type I IFN response. Cluster C exhibited greater APC co-inhibition, T cell co-inhibition, and type II IFN response, with cluster B ranking second (Figure 3K). The levels of infiltration by immune cells exhibited a comparable pattern across the three clusters (Figure 3L). Cluster 1 exhibited elevated levels of aDCx, mastocytes, NK cells, pDCs, Tfh, and Th1 cells. Cluster 2 exhibited elevated levels of CD8+ T cells and neutrophils. Cluster 3 exhibited elevated levels of B cells, dendritic cells (DCs), immature dendritic cells (iDCs), T helper cells, Th2 cells, and regulatory T cells (Treg). After assessing the genes related to immune check points, we noticed significant variations among the three clusters for all genes. Cluster 2 exhibited the greatest level of CD274 expression, while cluster 1 followed closely. On the other hand, cluster 3 displayed the lowest expression level. Additionally, the remaining genes demonstrated comparable patterns (Figure 3M).

Detecting gene clusters by comparing the overlap of differentially expressed genes within m7G clusters

Among the three m7G clusters, we identified genes that were expressed differently (DGEs), resulting in a total of 440 DGEs (Figure 4A). The GO function enrichment analysis showed that these differentially expressed genes (DGEs) were enriched in the detoxification of copper ions, stress response to copper ions, cellular response to copper ions, response to zinc ions, cellular response to zinc ions, and detoxification of inorganic compounds (Figure 4B). According to the KEGG pathway, these genes were found to be enriched in protein processing in the endoplasmic reticulum, proteasome, focal adhesion, oxidative phosphorylation, phagosome, human papillomavirus infection, and cell cycle, as shown in Figure 4C.

Figure 4. Identification of gene clusters-based DEGs among three m7G clusters. (A) Venn diagram indicating 440 DEGs among three m7G clusters. (B, C) GO and KEGG enrichment analyses based on 440 DEGs. (D) Consensus matrix identified the number of gene clustering. (E) PCA indicated three distinct gene clusters. (F) Kaplan-Meier analysis indicated significant differences in overall survival among three gene clusters. (G) Correlations of expression profiling, m7G clusters and gene clusters and clinical features. (H) Ggalluvial analysis indicated changes from m7G clusters to survival outcomes. (I–L) Differences in stromal, immune, estimate, and tumor purity. (M, N) Comparisons of immune cell infiltrations and immune functions. (O) Comparisons of immune checkpoint-related genes among three gene clusters.

By utilizing the 440 DGEs, we successfully conducted the clustering process and acquired three gene clusters (as depicted in Figure 4D). Figure 4E showed that the PCA revealed three distinct components. According to Figure 4F, the Kaplan-Meier curve showed that gene cluster 2 had the lowest overall survival (OS), followed by gene cluster 3, in comparison to gene cluster 1. Prognosis outcomes, m7G clusters, age, and 1p19q codeletion status exhibited distinct distributions in the clinical characteristics, as depicted in Figure 4G, 4H.

In addition, we assessed the tumor microenvironment and levels of immune infiltration. No significant variations were observed in ESTIMATE, stromal, and tumor purity across gene clusters 1, 2, and 3. However, a notable distinction was found in the immune score, with gene cluster 1 exhibiting a higher immune score compared to gene cluster 2 (Figure 4I–4L). Three gene clusters also exhibited notable variations in immune function and immune cells. Cluster 1 exhibited elevated levels of aDCs, mastocytes, natural killer cells, pDCs, Tfh, Th1 lymphocytes, and tumor-infiltrating lymphocytes. In Figure 4M, gene cluster 2 exhibited elevated levels of CD8+ T cells and neutrophils, while gene cluster 3 showed increased levels of B cells, DCs, iDCs, T helper cells, Th2 cells, and Treg. The highest check-point, cytolytic activity, T cell co-stimulation, and type I IFN response were observed in gene cluster 1, while gene cluster 3 had a slightly lower level, and gene cluster 2 exhibited the lowest immune function. In Figure 4N, gene cluster 3 exhibited the greatest APC co-inhibition, CCR, T cell co-inhibition, and type II IFN response, while gene cluster 2 followed closely behind. On the other hand, gene cluster 1 displayed the lowest levels. All genes exhibited significant differences in the immune check-related gene, as depicted in Figure 4O.

Creation and verification of a forecasting model relying on m7G regulators

Initially, we conducted the univariate Cox regression analysis and discovered 17 m7G genes that are linked to prognosis. In high-grade glioma (Figure 5A), EIF4E1B and NUDT11 acted as protective factors for OS, while the remaining factors posed a risk for Spatterwares, we conducted the LASSO regression analysis and discovered 13 genes that were incorporated into the ultimate prognostic model (AGO2, CYFIP1, DCP2, EIF4E1B, EIF4G3, GEMIN5, METTL1, NCBP1, NUDT11, NUDT16, SNUPN, WDR4, LARP1, (Figure 5B, 5C, Supplementary Table 2). Within the CGGA training cohort, we computed the risk score and subsequently categorized the participants into groups of high and low risk. According to the Kaplan-Meier plot, the high-risk category exhibited a less favorable prognosis compared to the low-risk category (Figure 5D, 5E). Additionally, PCA demonstrated distinct distributions for two different individuals (Figure 5F). Significant and comparable findings were also observed in the TCGA validation cohort (Figure 5G–5I).

Figure 5. Development and validation of m7G regulator-based prognostic signatures. (A) Forest plot of univariate cox regression in CGGA. (B, C) Parameter selection tuning by cross-validation using LASSO regression. (D, E) Kaplan-Meier curves of risk groups and distribution of risk score and patients in CGGA. (F) Display of two components by PCA in CGGA. (G) Kaplan-Meier curves of two risk groups and distribution of risk score and patients in TCGA. (H, I) Display of two components by PCA in TCGA. (J, K) GO and KEGG enrichment analyses based on 440 DEGs between high- and low-risk groups.

In order to examine the disparity in function, we acquired the differentially expressed genes (DEGs) between the high-risk and low-risk groups (Supplementary Table 3). The enrichment analysis using the GO function revealed that the high-risk groups showed enrichment in regionalization, regulation of the process that specifies patterns in the postsynaptic membrane potential, GABA-ergic synapse, GABA-A receptor complex, and activity (Figure 5J). Additionally, the KEGG analysis indicated that the high-risk group exhibited enrichment in GABAergic synapse, retrograde endocannabinoid signaling, interaction between neuroactive ligand receptors, and interaction between cytokines and cytokine receptors (Figure 5K).

The immune landscape and tumor microenvironment of various risk categories

The findings from our study showed that the high-risk group exhibited elevated ESTAMTE score, stromal score, and immune score, whereas the tumor purity was lower in the high-risk group (Figure 6A–6D). These results imply the presence of distinct tumor microenvironments in the high-risk and low-risk groups. In the high-risk group, there was a notable increase in the immune function (APC co-inhibition, APC c0-stimulation, CCR, Check-point, cytolytic activity, HLA, inflammation promoting, MHC class I, Para inflammation, T cell co-inhibition and stimulation, type I and II IFN response, Figure 6E). Likewise, the high-risk group exhibited elevated levels of immune cells such as aDCs, B cells, CD8+T cells, iDCs, macrophages, neutrophils, pDCs, T helper cells, Th1 and Th2 cells, and TIL Treg, compared to the low-risk group (Figure 6F). The correlation analysis showed that the risk score had a positive correlation with T cells gamma delta, M2 and M0 macrophages, and neutrophils (Figure 6G–6J). Conversely, the risk score exhibited an inverse correlation with B cells memory, B cells naïve, NK cell activated, and monocytes (as shown in Figure 6K–6N). The immune-regulated checkpoint genes showed a significant increase in the high-risk category (Figure 6O), suggesting a heightened immune status within this group.

Figure 6. Immune infiltration analysis of two risk groups. (A–D) Boxplot showing the comparative analysis of Estimation, stromal, immune scores, and tumor purity between two risk groups. (E, F) Comparisons of immune-related cells and functions of two risk groups. (G–N) Scatter plot showing risk score’s association with regulatory T cells gamma delta, Macrophages M2, M0, Neutrophils, B cells memory and naïve, NK cells activated, and monocytes. (O) Comparisons of immune checkpoint-related genes between high- and low-risk groups.

Independent analysis and risk assessment

To confirm the autonomy of the risk score in predicting outcomes, both univariate and multivariate cox regression analyses were conducted. The risk score in the training group showed a significant association with overall survival (univariate hazard ratio (HR) of 3.021, 95% confidence interval (CI) 2.453–3.722, P < 0.001; multivariate HR of 1.836, 95% CI 1.425–2.366, P < 0.001, as depicted in Figure 7A, 7B). Prognosis in high-grade glioma was also influenced by factors such as PRS categories, histological characteristics, patient age, chemotherapy treatment, presence of IDH mutation, and 1p19q status. Prognosis was not linked to the risk score in the validating group (Figure 7C, 7D). However, the age appears to be a separate prognostic factor in high-grade glioma. Figure 7E–7G displayed a calibration plot showing a precise match between the observed overall survival (OS) and the predicted probability from the nomogram for 1-year, 3-year, and 5-year OS.

Figure 7. Independent prognosis analysis of risk score. (A, B) Univariate and multivariate cox forest plot of the risk score in CGGA. (C, D) Univariate and multivariate cox forest plot of the risk score in TCGA. (E–G) Calibration plots of the nomogram to predict OS over one, three, and five years in the CGGA. (H) Nomograph predicting one, three, and five-year OS probabilities as per the m7G prognostic signature.

Using risk score, age, IDH mutation status, chemotherapy status, 1p19q codeletion status, histology, and PRS type, we developed a personalized risk evaluation tool. According to this tool, the projected overall survival rates for high-grade glioma at 1-year, 3-year, and 5-year were 0.423, 0.069, and 0.017 respectively (Figure 7H).

Chemotherapy sensitivity

In order to investigate the small molecular compound, we conducted a comparison of the IC50 levels of certain compounds between the high-risk and low-risk groups. In Figure 8A–8L, it was discovered that glioma cells may exhibit resistance to AZD7762, Paclitaxel, Nilotinib, JNK inhibitor VIII, JNK.9L, GSK269962A, Gefitinib, FTI.277, Doctaxel, Cyclopamine, Camptothecin, and Bicalutamide. The other potential compounds could be seen in the Supplementary Table 4.

Figure 8. Effect of m7G risk score on chemotherapy sensitivity. (A–L) AZD7762, Paclitaxel, Nilotinib, JNK Inhibitor VIII, JNK.9L, GSK269962A, Gefitinib, FTI.277, Docetaxel, Cyclopamine, Camptothecin, Bicalutamide.