An m1A/m6A/m7G/m5C regulator-mediated methylation modification pattern and Landscape of immune microenvironment infiltration characterization in Lower-Grade Glioma cohorts from three continents based on machine learning

doi:10.21203/rs.3.rs-2371294/v1

Download PDF

Research Article

An m1A/m6A/m7G/m5C regulator-mediated methylation modification pattern and Landscape of immune microenvironment infiltration characterization in Lower-Grade Glioma cohorts from three continents based on machine learning

https://doi.org/10.21203/rs.3.rs-2371294/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Although many studies have highlighted RNA modification processes such as N1-methyladenosine (m1A), N6-methyladenosine (m6A), N7-methylguanosine (m7G), and 5-methylcytosine methylation (m5C)’s role in the prognosis of patients suffering from different cancers, their prospective involvement in lower-grade gliomas (LGG) has not yet been outlined.

Methods: This work aims to assess the 64 genes related to m1A/m6A/m7G/m5C modification. Based on the expression of methylation-related regulators (MRRs), unsupervised clustering was conducted to identify molecular subtypes. The m1A/m6A/m7G/m5C modification patterns, tumor microenvironment (TME) cell infiltration features, and correlation with immune infiltration markers were assessed. Additionally, the first stage of MMR screening was conducted using univariate Cox analysis, and the prognostic model for the m1A/m6A/m7G/m5C risk score was constructed using different machine learning algorithms analysis.

Results: The m1A/m6A/m7G/m5C risk model, including five genes illustrated better prognostic ability for LGG in both the training and validation datasets, wherein the patients were classified into the low and high-risk groups. The LGG patients who were categorized into the high-risk groups displayed poor prognoses. In addition, the role played by five genes at the protein expression level was confirmed using immunohistochemical sections in the HPA database. Finally, functional analysis revealed the richness of pathways and biological processes related to MRR regulation and immune function.

Conclusion: An m1A/m6A/m7G/m5C-related risk model was developed and validated in this study to offer valuable new insights into the role played by m1A/m6A/m7G/m5C modification patterns in predicting the prognosis of LGG patients from three continents and developing better and improved treatment strategies for LGG.

m1A/m6A/m7G/m5C

RNA modification

lower-grade glioma

prognostic signature

tumor microenvironment

Gliomas are a primary type of primary malignant brain tumors [1], while lower-grade gliomas (LGG) are a type of heterogeneous neuroepithelial tumors originating in the central nervous system (CNS)-supporting glial cells [2]. LGG accounts for 6.4% of all the primary CNS tumors affecting adults, and the estimated incidence of astrocytoma and oligodendroglioma in the United States were 0.51/100000 and 0.25/100000, respectively [1]. Classically, the World Health Organization (WHO) categorized gliomas into four grades according to their histopathological characteristics. The WHO grade I and II gliomas were considered to be low-grade. The advancements in targeted sequencing and histological techniques have led to the identification of many molecular markers associated with tumorigenesis and development, malignant regression, treatment, and prognosis of gliomas, as well as some new tumor types that gradually emerged in clinical practice [3–5]. These findings led to a paradigm shift in glioma classification. The latest WHO CNS tumor classification released in 2021 focuses more on the role of molecular diagnosis in CNS tumor classification. LGG is no longer considered a lower histological grade but is reported to be more dependent on its MGMT promoter methylation state, isocitrate dehydrogenase (IDH) state, Histone 3 (H3) state, 1p/19q joint deletion, etc. Because of the convergence of several traits, WHO grade II and III gliomas are now classed as LGG [6]. Patients with IDH-mutant astrocytomas (9.3 years) and glioblastomas (3.6 years) have longer median overall survival (OS) duration than patients with IDH-wild type glioblastomas (1.9 years) and astrocytomas (1.2 years) [7], suggesting that the prognosis of LGG patients is better than that of glioblastoma patients, but LGG patients can still progress to higher grade gliomas [8].

Recurrence following surgical resection, radiotherapy, and chemotherapy is still quite common since the pathophysiological process and particular mechanisms of gliomas are still not fully understood. Additionally, there are few treatment options available, and the overall treatment outcome is not satisfactory. The prognosis for patients with LGG is still dismal, despite significant advancements in new anticancer medications that have increased the survival rate of patients with various cancers. Since the advent of the alkylating agent, temozolomide, roughly 15 years ago, there have been very few new pharmacologic additions for early or recurring treatments for gliomas [9]. Therefore, it is imperative to develop new biomarkers that may be used to anticipate an LGG patient's prognosis, understand the molecular mechanisms responsible for the onset and progression of LGG, and identify possible therapeutic targets.

Technology advancements like single-cell RNA sequencing (scRNA-seq) have made it possible to fully understand the considerable intra-tumoral heterogeneity present in gliomas [10], including genetic heterogeneity, epigenetic heterogeneity, and environmental heterogeneity. The diverse malignant cell phenotype offers several methods for coping with stress and treatment resistance that contribute to a highly-resilient form of the disease. To improve the prognosis of LGG patients, there is a higher need to understand glioma heterogeneity, particularly epigenetic modifications [11]. Traditional epigenetic studies have focused on DNA methylation, chromatin remodeling, histone modification, and non-coding RNA [12]. The epigenetic modification of RNA has received a lot of attention recently owing to its involvement in controlling a range of biological functions, including the onset and development of cancer. It has emerged as the most important post-transcriptional regulator in gene expression programs, wherein a majority of the programs include RNA base methylation [13, 14]. It may influence gene expression through RNA metabolism, shearing, stability, and translation. Furthermore, mounting data suggest that the dynamic RNA modification pathway is improperly controlled in human cancer and would provide a good target for cancer therapy [15]. N1-methyladenosine (m1A), N6-methyladenosine (m6A), N7-methylguanosine (m7G), and 5-Methylcytidine (m5C) are among the most prevalent RNA modifications [16]. The M6A modification is the most prevalent and extensively researched internal eukaryotic mRNA and non-coding RNA form of modification [17]. It has been suggested that some m6A components (writers, readers, and erasers) may serve as therapeutic targets for cancer [14]. The development of the tumor immune microenvironment (TME) and the differentiation of glioma stem cells (GSC) may depend primarily on M6A modification [18]. One of the earliest identified m6A demethylases, FTO, causes IDH mutation via the FTO/MYC/CEBPA signaling pathway, which leads to cancer [19]. M1A is primarily noted in tRNA and rRNA. M1A retains the tertiary structure of tRNA and rRNA and affects the translation process by combining with RNA containing m1A. YTHDC1, YTHDF1, YTHDF2, and YTHDF3 act as readers, and ALKBH1 and ALKBH3 are seen to demethylate m1A mRNA modification [20]. Studies have revealed that M1A methylation affects the development of the immunological microenvironment, and also influences the prognosis of patients suffering from ovarian cancer [21] and gliomas [22]. Prokaryotes and eukaryotes both contain a common modification of tRNA called m7G that is mediated by methyltransferase-like proteins such as methyltransferase-like protein-1 (METTL1) and WD repeat domain 4 (WDR4). Increasing numbers of studies have highlighted the need to investigate the involvement of m7G in cancer inhibition [23]. A study has demonstrated that the NOL1/NOP2/SUN domain (NSUN) family proteins (NSUN1, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, and NSUN7) and the DNA methyltransferase (DNMT) homologue DNMT2 catalyze the modification of the m5C RNA [20]. According to Wang et al. [24] m5C may modify PKM2 mRNA in the HIF-1 / ALYREF /PKM2 axis, which can increase glucose metabolism in bladder cancer cells. This presents a novel and promising therapeutic target for treating bladder cancer. Recently, a few researchers highlighted the essential role of Methylation-related regulators (MRRs) in the progression and metastasis of malignant tumors. A proper understanding of the different forms of RNA modifications occurring in malignant tumor cells would allow for their accurate prognosis and precise treatment [25]. Very few studies have investigated the specific role of MRRs in LGG, and the general mechanism of m1A/m6A/m7G/m5C in LGG is still unknown. As a result, the relationship between MRRs and LGG cannot be fully summarized. A thorough investigation of the m1A/m6A/m7G/m5C-related genes in LGG is therefore required to develop an accurate prognosis model and identify potential prognostic biomarkers.

In this study, the RNA sequence data and clinical data from 975 human LGG patients were extracted using The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases. Afterward, large-scale bioinformatics analysis was conducted to investigate the possible role of m1A/m6A/m7G/m5C-related genes in LGG. Thereafter, the m1A/m6A/m7G/m5C risk score model was constructed, and the developed model was validated using the CGGA data set. Machine learning techniques were employed to determine the diagnostic capacities of the genes used for constructing risk models in predicting tumorigenesis. Additionally, the functions of genes associated with m1A/m6A/m7G/m5C in the TME of LGG, immune cell infiltration, clinical traits, and chemosensitivity were also investigated. These findings led to the conclusion that the m1A/m6A/m7G/m5C score could be utilized as a promising biomarker to anticipate the clinical outcomes of LGG patients and presents a new approach for treating LGG.

2.1 Data acquisition and pre-processing of Lower-Grade Gliomas

In this study, the RNA-Seq data from LGG samples with WHO classification II-III were extracted from TCGA and CGGA databases and included for further analysis. On the other hand, data with repeated sequencing, OS time < 1 day, and non-primary LGG were excluded. Finally, 481 samples from the TCGA-LGG cohort, 162 samples from the CGGA-325 cohort, and 332 samples from the CGGA-693 cohort were included in this study. After TPM transformation and log (x + 1) normalization, the samples from the above databases were background corrected, normalized, and the expressions were estimated using the combat function of the sva software. Finally, a meta-cohort was generated, which contained 975 LGG samples. In addition, the scRNA-seq data for GSE138794 were downloaded from the GEO database, keeping the LGG samples for GSM4119535, GSM4119536, and GSM4119538. All MRRs related to M1A, M6A, M7G, and M5C were identified from previous researches [26, 27]. Finally, 64 genes were annotated in the meta-cohort (Supplementary File 1). The immunohistochemical data related to ALKBH3, DCP2, EIF3D, IGF2BP2, and NSUN7 were derived from the HPA database and collected from the normal cerebral cortex and LGG tumor samples, respectively. The whole exon sequencing data of the TCGA-LGG cohort were also retrieved from the TCGA database, while the MutSigCV algorithm was employed for selecting the oncogenes with a mutation frequency higher than the background. Furthermore, the Maftools package was used for displaying the mutation map of oncogenes.

2.2 Unsupervised Clustering of m1A/ m6A/m7G/m5C modification pattern

The prognostic significance of each MRR in the meta-cohort was determined by means of univariate cox regression analysis. Based on predictive MRRs (p < 0.05) and DEGs (p < 0.05), unsupervised clustering analysis and principal component analysis (PCA) were conducted to ascertain if each subtype was relatively independent of other subtypes. The number of clusters was calculated using the "consensusClusterPlus" tool in R software. The stability of the subtypes was verified using a log-rank test on 100 replicates with a pltem = 0.8.

2.3 Enrichment analysis for distinct m1A/ m6A/m7G/m5C modification patterns

The differentially expressed genes (DEGs) across different subtypes were examined using the limma program. After DEGs were identified, the gene pathway was annotated using the clusterProfiler package for Gene Ontology (GO) and the Kyoto Encyclopedia of Gene and Genome (KEGG) (adj. p༜0.05, |logFC|༞1). Significant enrichment pathways were defined as P-value < 0.05 and q-value < 0.05. Additionally, GSVA was also employed to assess the differences (variations) in biological pathways across subtypes, whereas c2.cp.kegg.v7.0.symbol and h.all.v7.4.symbols were utilized as reference gene sets, and FDR < 0.05 was considered as the threshold.

2.4 Identification and verification of MMRs risk scores

The CGA-LGG cohort was used for modeling and the CGGA-693 and CGGA-325 cohorts were used for validation. Firstly, the univariate cox regression model was employed for screening the MRRs in the modeling cohort (p༜0.05), followed by the use of the Least absolute shrinkage and selection operator (LASSO) model to remove the redundant genes in MRRs. Finally, the risk scoring formula in the multivariate Cox regression analysis was used for assessing the gene expression value and integration coefficient. The CGGA-693 and CGGA-325 cohorts were classified into high and low-risk groups as per the median risk scores in the TCGA-LGG cohort. The time-dependent Receiver operating characteristic (ROC) curve was employed for comparing the predictive accuracy of the risk scores using conventional clinicopathological parameters.

2.5 Evaluation of Tumor Microenvironment (TME) and Immune Cells Infiltration in lower-grade gliomas

A variety of algorithms, including ssGSEA, TIMER, CIBERSORT, QUANTISEQ, MCP-counter, XCELL, and EPIC were used in the study to determine the number of immune cells in different samples. Immune score and stromal score were derived using the ESTIMATE algorithm to reflect the condition of the TME. In addition, the IMvigor-210 cohort (treated with immune checkpoints) was downloaded from previous literature to generate risk scores and verify the effectiveness of risk scores in the immunotherapy dataset.

2.6 Single-Cell Sequencing analysis

The original data were converted into Seurat objects, retaining the genes that were expressed in a minimum of 3 cells and the cells that expressed at least 200 genes. A harmony algorithm was used to combine different datasets. The downstream analysis determined the cell clustering profiles based on the corrected Harmony embeddings and FindNeighbors functions. The t-sne function was used to visualize the cell cluster and a singleR package was employed for the purpose of annotating the cell types.

2.7 Machine Learning Pipelines

The normal sample group consisted of data from 103 normal cortical samples derived from the GTEx project, whereas the disease occurrence group included 481 samples from the TCGA-LGG cohort. The input matrix consisted of RNASeq data in TPM format from the TCGA and GTEx datasets that were uniformly processed by the Toil process (https://xenabrowser.net/datapages/). Machine learning channels were used to incorporate the genes used in risk model construction to determine the diagnostic capacity of risk genes in tumorigenesis prediction. Various machine learning models (Bayes, DT, FDA, GBM, NNET, RF, SVM, and Logistic) were constructed using the default settings. The superiorities and drawbacks of various models were compared using the AUC values, and the contributions of various genes in the above models were also compared.

2.8 Drug sensitivity analysis

Chemotherapeutic medications were retrieved from the Genome of Drug Sensitivity in Cancer (GDSC) database, and the IC50 value was obtained using the prophetic package in the R program.

2.9 Molecular docking

Molecular docking with AutoDock Vina 1.1.2 was used to virtually screen the m1A/m6A/m7G/m5C-related Hub genes and identify the most likely optimal ligand. The component structures were derived from Pubchem. The following optimization settings were imported into the SYBYL-X program during the molecular optimization process: energy gradient was confined to 0.005 kcal/(mol-A), the Tripos force field was employed, the Gasteiger-Hückel charge was chosen, the maximum iteration factor was set to 10,000, and all other parameters were left at their default settings. The target structure was downloaded from RCSB (https://www1.rcsb.org). The Surflex-Dock module of SYBYL-X was used to remove the protein's crystalline water, process the terminal residues, and extract, hydrogenate, and optimize the energy of the ligands to construct a binding pocket, followed by molecular docking (Table 1).

Table 1

Information on the four targets(ALKBH3, IGF2BP2, DCP2, EIF3D).
Target	PDB ID	X	Y	Z
ALKBH3	4RFR	-2.36871	2.381238	3.126286
IGF2BP2	6ROL	27.17724	-11.726	-18.9552
DCP2	5J3Y	28.67744	-26.5806	-32.5834
EIF3D	5K4D	20.19526	55.83009	-7.53488

2.10 Statistical Analysis

R software was used to perform all statistical analyses. Pairwise comparisons were conducted using the Wilcoxon test. The OS rates across all the groups were compared using the Kaplan-Meier (KM) analysis with the log-rank test. Statistics were deemed significant at P < 0.05.

3.1. Modification patterns mediated by m1A/m6A/m7G/m5C-regulators in lower-grade gliomas

PCA analysis revealed that batch effects were more effectively eliminated when CGGA-325, CGGA-693, and TCGA-LGG were included in a Meta cohort and the COMBAT technique was used (Fig. 1A-B). Subsequently, a univariate cox regression analysis of 64 MRRs in the Meta cohort was performed, and the results indicated that 42 MRRs exhibited a prognostic indication value, with YTHDF2 being probably the strongest risk factor and FMR1 being the best protective factor in the LGG cohort (Fig. 1C). Patients were subsequently classified based on the expression of prognostic MRRs using the ConsensusClusterPlus package, and the best classification was achieved when the K value was set as 2. The heat map shows the distribution of clinical data of patients having different molecular subtypes in the CCGA-LGG (Fig. 1D) and TCGA-LGG cohorts (Fig. 1E), respectively. Interestingly, the distribution of cluster A and cluster B was seen to be relatively consistent across the cohorts and the majority of MRRs were significantly upregulated in cluster B.

3.2. Exploring the m1A/m6A/m7G/m5C modification pattern regulatory genes:

To explore the indicative effect of different clusters on biological functions and prognosis, the data were analyzed using several bioinformatic processes. Cluster A in the Meta cohort consisted of 774 cases, whereas Cluster B included 201 cases. The prognostic analysis indicated that people with Cluster B molecular subtypes displayed a significant survival disadvantage (Fig. 2A). PCA analysis further showed that both the molecular subtypes had a comparatively discrete nature (Fig. 2B). Comparing the expression of HLA molecules and ICI mRNA in various molecular subtypes revealed that cluster B showed higher mRNA expression levels for most HLA and ICI (Fig. 2C-D). In addition, ssGSEA analysis revealed the TME status in several molecular subtypes, and it was noted that cluster B was likely to display hot tumor characteristics, such as elevated levels of activated NK cells, macrophages, and activated T cells, along with a very active TME (Fig. 2E). Finally, the ESTIMATE algorithm also confirmed the hypothesis that the overall immune score was significantly upregulated in cluster B (Fig. 2F).

3.3. GSVA enrichment analysis

Subsequently, GSVA enrichment analysis was carried out based on different gene sets, and the results indicated that cluster B presented significant enrichment in KEGG pathways associated with HRD such as cell cycle pathway and DNA mismatch repair compared to cluster A (Fig. 3A). Furthermore, it was also noted that cluster B was significantly activated in the classic cancer-related pathways (Fig. 3B). This illustrates the primary cause of the various prognoses for the molecular subtypes that were identified using MRRs. In addition, a differential analysis between the two molecular subtypes was carried out. After analyzing the 48 identified DEGs (Fig. 3C), it was noted that the primary enrichment pathways of the aforementioned genes may be linked to biological processes such as focal adhesion, ECM-receptor interaction, etc. (Fig. 3D). These findings may offer evidence that could support different prognosis profiles across various subtypes and variations in TME. Thereafter, patients were classified depending on the expression levels of different prognostic DEGs, and the best classification was achieved when the K value = 3. Three different regulatory subtypes were finally identified (Fig. 3E), of which geneCluster B displayed the worst prognosis while geneCluster A presented the best prognosis (Fig. 3F). The heatmap analysis showed that geneCluster B had higher expression of the related MRRs and distribution of clinical characteristics across different regulatory subtypes (Fig. 3G).

3.4. Construction of m1A/m6A/m7G/m5C Regulated Gene-Related Risk Model

In this study, TCGA-LGG was used as a training cohort with the OS period as the outcome. In the first step, the MRRs were screened using a univariate Cox regression model in the training pair-up (p < 0.05), and 34 MRRs were prognostically indicated (Fig. 4A). After redundant genes were eliminated using the LASSO model, multivariate Cox regression analysis was conducted to determine the coefficients of every gene (Fig. 4B-C). Finally, a signature of five genes was established. Each risk score for every patient was determined using the following formula: riskscore = (-0.4783 × ALKBH3 expression) + (0.9376 × NSUN7 expression) + (0.4419 × IGF2BP2 expression) + (1.5324 × DCP2 expression) + (-0.6637 × EIF3D expression). Based on the median values of the risk scores in the TCGA dataset, two risk groups for the LGG patients were established, i.e., the low-risk group and the high-risk group. After analyzing the data, it was noted that NSUN7, IGF2BP2, and DCP2 genes were expressed at high levels in the high-risk patients in the TCGA-LGG cohort, whereas EIF3D and ALKBH3 were underexpressed (Fig. 4D). The KM curve indicated a significant decrease in the survival probability of the high-risk patients in the TCGA-LGG cohort (Fig. 4E). Furthermore, the ROC curves also presented the survival prediction values over 1 year (AUC = 0.886), 3 years (AUC = 0.854), and 5 years (AUC = 0.775) (Fig. 4F). Surprisingly, a significant agreement was noted between the gene expression levels and patient prognosis in the CGGA-325 Fig. 4G-I) and CGGA-693 datasets (Fig. 4J-L). Lastly, the OS scatter plots demonstrated that patients with a higher risk score in various cohorts were more likely to exhibit a fatal outcome (Fig. 4M).

3.5. Landscape of m1A/m6A/m7G/m5C Regulators in LGG

To further characterize the five risk genes involved in the risk model: ALKBH3, NSUN7, IGF2BP2, DCP2, and EIF3D, the location of MRRs on the chromosomes was determined (Fig. 5A). The network of five MRRs outlined the combined gene interactions and the significance of those interactions for patients' genetic prognosis (Fig. 5B), with ALKBH3 and EIF3D having protective value and NSUN7, IGF2BP2, and DCP2 having risk predictive value (Fig. 5C). The normal samples from the GTEx database were combined and compared, and the results indicated that the ALKBH3, DCP2, and EIF3D genes were significantly upregulated in tumor samples, IGF2BP2 was significantly downregulated, whereas NSUN7 expression did not show any significant differences (Fig. 5D). Additionally, CNV modifications were common in MRRs, where copy number amplification was common in IGF2BP2, whereas CNV deletion frequency was primarily noted in EIF3D (Fig. 5E). However, the above five genes were less frequently mutated, with only EIF3D and IGF2BP2 detected in 0.2% of the TCGA-LGG cohort, while no significant mutations were seen in the remaining MRRs (Fig. 5F). To further examine the cell type distribution of the five MRRs, three LGG single-cell samples were integrated (Fig. 6A-B) and subsequently four cell types (i.e., astrocytes, endothelial cells, neurons, and macrophages) were annotated using singleR software (Fig. 6C). The cellular distribution of different MRRs was subsequently demonstrated in t-sne plots (Fig. 6D-H), which showed that ALKBH3/EIF3D/DCP2 were more evenly distributed in different cell types, NSUN7 was mostly distributed in glial cells, while IGF2BP2 was expressed at low levels in a majority of cell types.

3.6. m1A/m6A/m7G/m5C Regulators as potential biomarkers for LGG

The input matrix was created using the processed RNA-seq data in TPM format from TCGA and GTEx databases to incorporate genes that were involved in the construction of risk models into the machine-learning pipeline. This helped in identifying the diagnostic power of risk genes in predicting tumorigenesis. Multiple machine learning models were built using default parameters: Bayes, DT, FDA, GBM, NNET, RF, SVM, Logistic, and the contributions of different MRRs in DT, FDA, GBM, NNET, RF, and Logistic were compared.

The results highlighted the significance of DCP2 in the three models, i.e., DT, FDA, and GBM; EIF3D showed the highest importance in RF, and Logistic; while ALKBH3 presented the highest importance in the NNET model (Fig. 7A). The AUC values were used to compare the advantages and disadvantages of various models, and the results showed that the five MRRs in the RF model had the best diagnostic ability for identifying the tumor samples with AUC = 0.938313 (Fig. 7B). Finally, the hypothesis proposed in this study was further confirmed at the protein expression level utilizing immunohistochemical sections from the HPA database (Fig. 7C-G).

3.7. Prognostic value of m1A/m6A/m7G/m5C risk scores in individual LGG patients

The Sankey diagram showed a substantial association of risk with molecular subtypes. Most patients in cluster A, geneCluster A, and low-risk groups had good survival rates (Fig. 8A), whereas patients who displayed a better prognosis molecular subtype had a lower risk score (Fig. 8B-C). Correlation and difference analysis of hallmark-related pathways and classic cancer-related pathways were significantly activated with increasing risk scores, and boxplots showed differences in activated pathways across risk groups (Fig. 8D-E). In addition, the combined tumor mutational burden (TMB) status allowed further stratification of patients for survival risk, in which patients with high TMB-high risk scores displayed the worst prognosis (Fig. 8F) and high-risk patients presented higher TMB scores (Fig. 8G). Finally, the results of the full exon sequencing data showed that there was no significant difference in Top mutations between both the risk groups, i.e., IDH1, TP53, and ATRX (Fig. 8H-I).

3.8. Identifying the indicative role of m1A/m6A/m7G/m5C scores in TME

The proportion of lymphocyte-depleted subtypes was considerably higher in the high-risk patients in comparison to the pre-existing immune subtypes (Fig. 9A), and the risk scores were seen to be varied across molecular subtypes (Fig. 9B). In addition, the risk formula could distinguish the survival risk in the immunotherapy cohort (Fig. 9C), and the results indicated that the immunization score increased with an increasing risk score in the ESTIMATE analysis (Fig. 9D). Finally, different algorithms were simultaneously used, i.e., TIMER, CIBERSORT, MCP-counter, QUANTISEQ, XCELL, and EPIC to estimate the immune cell infiltration fraction in different samples. The differential analysis heat map showed that the high-risk group had a very active TME (Fig. 9E), whereas the correlation analysis showed that the infiltration fraction of killer immune cells increased significantly with an increasing risk score (Fig. 9F). Even though most studies suggested that LGG is resistant to chemotherapy, this study investigated whether the m1A/m6A/m7G/m5C score-based grouping might be indicative of conventional cytotoxic drugs. Hence, the prophetic package was used to calculate the IC50 values for different drugs. Chemotherapeutic medications (AS601245, Acadesine, Axitinib, Veliparib dihydrochloride, Rucaparib phosphate, Navitoclax, A-443654, A-770041, Saracatinib AZ628) were found to be more effective in the high-m1A/m6A/m7G/m5C score group (Figure S1 A-J).

3.9. Investigating the best-fitting compounds on m1A/m6A/m7G/m5C scores hub genes

From a drug perspective, various drugs have varying affinities for important targets. It has been noted that CID 9549184, CID 11676786, and CID 10172943 have strong binding abilities to every target, but each drug has a slightly higher affinity for IGF2BP2 and DCP2. CID 11676786 has a strong affinity and binding energy of -9.9 kcal/mol for IGF2BP2, whereas CID 9549184 has a strong affinity (Table 2).

Table 2

An Affinity of different drugs for key targets.
Name	ALKBH3	IGF2BP2	DCP2	EIF3D
CID_10109823	-7.5	-8.5	-8.5	-8.8
CID_10172943	-9	-9	-9.7	-9.8
Saracatinib	-8.5	-7.4	-8.3	-9
CID_11676786	-9.1	-9.9	-9.1	-9.6
Veliparib_dihydrochloride	-6.5	-7.5	-7.3	-7.9
Acadesine	-7.6	-6.9	-6.4	-7
Navitoclax	-8.2	-9.2	-9.9	-8.2
Axitinib	-8.2	-8.5	-9.2	-9
CID_9549184	-9.6	-9.8	-10.9	-8.8
Rucaparib_phosphate	-8.2	-8.6	-8.7	-8.5

The results of a Pymol visualization research showed that CID_9549184 bound to a cavity on the ALKBH3 protein surface by forming hydrogen bonds with two amino acid residues (SER129 and TYR76) within the binding pocket. On the other hand, DCP2 could bind to the ASN72, PRO74, and LYS218 amino acids, whereas cid_10172943 could form hydrogen bonds to the EIF3D amino acid residues of ARG271 and LYS282. CID_9549184 formed hydrogen bonds with IGF2BP2 amino acid residues of ASN50, GLU506, HIS564 (Fig. 10A-D).

Growing evidence suggests that RNA modifications such as m1A/m6A/m7G/m5C play a crucial role in both the physiological and pathological processes, particularly in the onset, advancement, migration, and invasion of cancer [28]. Inflammation, innate immunity, and anticancer effects are influenced by the interaction between the MRRs and RNA modifications [29]. Immunotherapy has displayed many positive outcomes in the treatment of several malignancies. A randomized, multi-center, open-label Phase 3 trial was conducted, where the researchers observed that Nivolumab, a programmed cell death (PD)-1 inhibitor, demonstrated a positive effect when it was administered to patients with adenocarcinoma of the gastroesophageal junction, stomach, and esophagus while maintaining an acceptable safety profile [30]. Cloughesy et al. [31] carried out a study where they investigated the PD-1/L1 inhibitors for glioma and observed that neoadjuvant injection of PD-1 blockade increases both the local and systemic antitumor immune responses and could represent a more effective strategy for treating the uniformly deadly brain tumors. The TME of LGG, however, necessitates more research. Furthermore, the majority of previous studies have focused on a single regulatory molecule, and it is yet unclear how the combined activities of different m1A/m6A/m7G/m5C regulators affect the overall characteristics. Additional research is required to better understand the glioma-related m1A/m6A/m7G/m5C gene modification patterns, prognosis predictive roles, and potential effects on the TME.

In this study, the researchers explored the m1A/m6A/m7G/m5C modification patterns based on 64 regulators of MMRs and identified new molecular subtypes based on m1A/m6A/m7G/m5C-related regulators. These two subtypes showed different molecular characteristics, with cluster B showing a significant survival disadvantage and cluster B showing active m1A/m6A/m7G/m5C modifications and enrichment mainly in cell cycle-related pathways, DNA mismatch repair, etc. Then, functional enrichment analysis showed that the prognostic signature was significantly linked to the immune cells and immune-related pathways, which is consistent with poor prognosis, again demonstrating the main reason for the differentiation of different prognosis based on m1A/m6A/m7G/m5C-related regulators.

A novel m1A/m6A/m7G/m5C risk score model was developed in this study, which consisted of five genes (ALKBH3, NSUN7, IGF2BP2, DCP2, EIF3D) based on the m1A/m6A/m7G/m5C modification pattern. This was seen to be significantly associated with LGG prognosis, and a strong relationship was noted between risk groups and molecular subtypes. Several studies indicated that high TMB was associated with clinical benefits and prolonged OS duration in patients [32, 33]. Therefore, TMB could be employed as a predictive biomarker for LGG patients, which facilitated the clinical decision-making process. The results in this study showed that the m1A/m6A/m7G/m5C scores were positively and significantly correlated with TMB, i.e., the low-m1A/m6A/m7G/m5C score group showed lower TMB scores, whereas the Low-TMB group represented a better prognostic outcome when low-m1A/m6A/m7G/m5C scores were combined with low-TMB. This validated the superiority of using the m1A/m6A/m7G/m5C score in the immunotherapy of LGG patients. Additionally, it was noted that the high m1A/m6A/m7G/m5C score group was more sensitive to chemotherapeutic agents (AS601245, Acadesine, Rucaparib phosphate, Veliparib dihydrochloride, Navitoclax, A-770041, A-443654, Saracatinib, AZ628, Axitinib), which could improve the clinical efficacy of LGG.

Several researchers have found that the m1A/m6A/m7G/m5C-related regulatory hub genes are essential for the progression and metastasis of tumors. ALKBH3 belongs to the alpha-ketoglutarate-dependent dioxygenase AlkB homolog family and can eliminate alkyl adducts from nucleobases using the oxidative dealkylation process [34]. The sequencing and functional studies showed that ATP5D, a significant subunit of adenosine-5'-triphosphate synthase, was involved in m1A demethylase ALKBH3-regulated glycolysis in malignant cells [35]. It was noted that the m1A demethylase ALKBH3 regulates glycolysis in malignant cells in a demethylation activity-based manner. Kuang et al. found that ALKBH3 depletion enhanced Aurora A mRNA degradation and impeded its translation. They also showed that ALKBH3 morphants displayed cilia abnormalities that were greatly reversed by wild-type ALKBH3. Thereafter, they concluded that the ALKBH3-dependent m1A demethylation process was involved in Aurora A mRNA regulation [36]. It was also shown that ASCC3-mediated breakdown of colliding ribosomes allows the demethylation of aberrant m1A and m3C by ALKBH3. This presented the first evidence for selective purification of aberrant mRNA methyl bases and the RNA-mediated effects of chemical alkylating agents commonly used in clinical practice require further investigation [37]. Additionally, You et al. detected the presence of 1, N6-dimethyladenosine (m1,6A), a new adenosine bimethylation modification, in the tRNAs of living organisms. The results of this study showed that m1,6A was located at position 58 in the tRNA in the mammalian cells and tissues. TRMT6/61A can catalyze the formation of m1,6A in tRNAs, and ALKBH3 can demethylate m1,6A. The identification of m1, 6A broadens the range of RNA modifications and might indicate a novel mechanism for tRNA modification-mediated gene regulation [38]. EIF3D plays a crucial role in the early assembly stage of the translation machinery by acting as a molecular link between the small ribosomal subunit (the 40S) and eukaryotic translation initiation factor 3 (eIF3) [39]. According to Lamper, when human cells were under metabolic stress, the non-classical 5' cap-binding protein EIF3D was activated. The members of the mammalian target of the rapamycin (mTOR) pathway, whose directed translational adaptation is essential for cell survival under chronic glucose deprivation, are among the factors necessary for glucose homeostasis that are enriched in the genetic program controlled by EIF3d [40]. Cieśla et al. revealed that the core spliceosome component, i.e., SF3A3 regulates the transformation in an EIF3D-dependent manner via RNA stem-loop in a translation-based program, which fuels the tumorigenic potential of MYC in vivo and predicts the molecular and phenotypic characteristics of aggressive human breast cancer [41]. Colon cancer [42], acute myeloid leukemia [43] and other cancers have all been linked to the expression of EIF3D, which enhances protein synthesis. The human insulin-like growth factor 2 (IGF2) mRNA-binding protein family (IMPs/IGF2BPs) is engaged in a variety of biological processes, such as tumorigenesis, stemness, and development [44]. IGF2BP plays a crucial role as an m6A reader in post-transcriptional gene regulation and cancer biology [45]. Yao et al. revealed that circEZH2 interacts with IGF2BP2 and prevents its ubiquitination-dependent degradation, thereby promoting CREB1 mRNA stability in colorectal cancer (CRC) and thus aggravating the CRC progression [46]. Aguilo et al. noted that the NSUN7 deletion leads to target gene depletion of peroxisome proliferator-activated receptor gamma coactivator-1α (PGC-1α) in the hepatocyte model system. They also observed that the interaction between NSUN7 and PGC-1α highlighted the enrichment of m5C-modified eRNAs on enhancers of specific target genes [47]. Additionally, an earlier study suggested a significant association between NSUN7 immunoreactivity and the OS rate in patients with Ewing sarcoma, suggesting that NSUN7 could serve as an independent prognostic marker [48]. DCP22 belongs to the Nudix hydrolase family, which is an mRNA-decapping enzyme that displays a clear specificity for m7G cap RNA and interacts with other mRNA attenuators [49]. It has been shown that DCP2 acts as a miR-4293 target and that its expression levels were suppressed by miR-4293. The interaction between miRNAs and lncRNAs via DCP2 helped in regulating the pathogenesis of cancer [50]. In summary, ALKBH3, NSUN7, IGF2BP2, DCP2, and EIF3D genes, included in the MMR risk score model, were associated with m1A/m6A/m7G/m5C. This study showed that ALKBH3 and EIF3D genes were associated with an improved prognosis of LGG patients and displayed a protective value, whereas the NSUN7, IGF2BP2, and DCP2 genes exhibited a risk predictive value. It is still unknown how all the above genes interact during m1A/m6A/m7G/m5C modifications.

The functional processes of the LGG-related genes m1A/m6A/m7G/m5C have not yet been studied. ALKBH3, NSUN7, IGF2BP2, DCP2, and EIF3D are seen to be the major m1A/m6A/m7G/m5C-related genes in LGG, and the results in this study showed that m1A/m6A/m7G/m5C modifications were significantly associated with LGG. Of these, ALKBH3 and EIF3D exhibited a protective role, whereas NSUN7, IGF2BP2, and DCP2 showed a risk predictive value. The RNA modification pattern in LGG patients can be evaluated using the m1A/m6A/m7G/m5C score, which can also be used to elaborate on the characteristics of the TME cell invasion. The results of this study have increased our understanding of the LGG immunophenotype and improved the effectiveness of clinical treatments. The m1A/m6A/m7G/m5C score may therefore be used as an independent predictive biomarker to direct LGG clinical treatment, and the inherent immune regulatory mechanisms in patients with active m1A/m6A/m7G/m5C changes may support the functional analyses and clinical trials. However, the inherent limitations of genetic heterogeneity could have an impact on the findings and future clinical applications. Additionally, this study primarily focused on the genetic expression of hypothesized genes using a bioinformatic approach, but the protein expression needs to be validated. In conclusion, more research is needed to determine how the specific processes used by the m1A/m6A/m7G/m5C-related genes, such as ALKBH3, NSUN7, IGF2BP2, DCP2, and EIF3D, influence the onset, progression, and patient prognosis of LGG.

In conclusion, the results of this study showed that the genes related to m1A/m6A/m7G/m5C were significantly associated with LGG. LGG is significantly influenced by M1A/m6A/m7G/m5C modifications, which also have considerable potential for clinical diagnosis and treatment. Post-transcriptional editing of defective m1A/m6A/m7G/m5C loci in combination with chemotherapy or immunotherapy may improve therapeutic outcomes in the future because of its role in chemoresistance and TME reconfiguration. The investigation into the M1A/m6A/m7G/m5C modification is a research area that needs to be investigated further. However, based on its significance in tumor growth and value in clinicopathology, M1A/m6A/m7G/m5C modifications offer promising strategies for the therapeutic management and therapy of LGG.

m1A: N1-methyladenosine; m6A: N6-methyladenosine; m7G: N7-methylguanosine; m5C: 5-methylcytosine methylation; LGG: Lower-grade glioma; TCGA: The Cancer Genome Atlas; CGGA: The Chinese Glioma Genome Atlas; MRRs: methylation-related regulators; TME: Tumor microenvironment; HPA: Human Protein Atlas; CNS: central nervous system; WHO: the World Health Organization; IDH: isocitrate dehydrogenase; H3: Histone 3; OS: overall survival; scRNA-seq: single-cell RNA sequencing; GSC: glioma stem cells; METTL1: methyltransferase-like protein-1; WDR4: WD repeat domain 4; DNMT: DNA methyltransferase; PCA: principal component analysis; GSVA: Gene set variation analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; HRD: Homologous Recombination Deficiency; DEGs: differentially expressed genes; LASSO: Least Absolute Shrinkage and Selection Operator; GTEx: Genotype-Tissue Expression; GDSC: Genome of Drug Sensitivity in Cancer; CNV: copy number variation; TMB: Tumor Mutational Burden; PD-1: Programmed cell death protein 1; PD-L1:Programmed cell death protein ligand 1.

Availability of data and material

Data and download URLs involved in this study had been described in detail in the materials and methods section. The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Code availability

Bioinformatics analysis of this study was based on R 4.0.3, and the involved R packages were described in detail in the materials and methods section. We uploaded the R code in the supplementary file. All codes were written and reviewed by AM, respectively. AM had the ultimate power of interpretation and ownership of the code.

Ethics approval and consent to paticipate

All datasets used in the present study were downloaded from public databases, including TCGA, GEO and CGGA database. These public databases allow researchers to download and analyze public datasets for scientific purposes and thus no ethical approval nor informed consent was required. The current research follows the TCGA, GEO and CGGA data access policies and publication guidelines. All methods/protocols were performed in accordance with the relevant guidelines and regulations.

Consent to participate

All authors voluntarily participated in this study.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Funding

This work was supported by the Xinjiang Uygur Autonomous Region Regional Collaborative Innovation Project (Shanghai Cooperation Organization Science and Technology Partnership Program and International Science and Technology Cooperation Program) (grant number: 2022E01068) (grant number: 2021E00131), Xinjiang Natural Science Foundation Surface Project (grant number: 2021D01A02), (grant number: 2020D01C235),Xinjiang Uygur Autonomous Region "Young Scientific and Technological Talents - Rural Revitalization" Project (grant number: WJWY-XCZX202208),

Author Contributions

AM & YB: Conceptualization, Methodology, Validation, Investigation, Supervision, Visualization, Writing - original draft, Writing- Reviewing. MT, AA, XC, YP, YA, DA, LJ, YW, ZW, and ZF participated in the coordination of data acquisition and data analysis and reviewed the manuscript. All codes were written and reviewed by AM and YB.

Acknowledgments

We are grateful to the contributors to the public databases used in this study.

Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2012-2016. Neuro-oncology. 2019;21.
Brat DJ, Verhaak RGW, Aldape KD, Yung WKA, Salama SR, Cooper LAD, et al. Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N Engl J Med. 2015;372:2481-98.
Eckel-Passow JE, Lachance DH, Molinaro AM, Walsh KM, Decker PA, Sicotte H, et al. Glioma Groups Based on 1p/19q, IDH, and TERT Promoter Mutations in Tumors. N Engl J Med. 2015;372:2499-508.
Xu J, Liu F, Li Y, Shen L. A 1p/19q Codeletion-Associated Immune Signature for Predicting Lower Grade Glioma Prognosis. Cell Mol Neurobiol. 2022;42:709-22.
Chen CCL, Deshmukh S, Jessa S, Hadjadj D, Lisi V, Andrade AF, et al. Histone H3.3G34-Mutant Interneuron Progenitors Co-opt PDGFRA for Gliomagenesis. Cell. 2020;183.
Louis D, Perry A, Wesseling P, Brat D, Cree I, Figarella-Branger D, et al. The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro-oncology. 2021;23:1231-51.
Pekmezci M, Rice T, Molinaro AM, Walsh KM, Decker PA, Hansen H, et al. Adult infiltrating gliomas with WHO 2016 integrated diagnosis: additional prognostic roles of ATRX and TERT. Acta Neuropathol. 2017;133:1001-16.
Dhawan S, Patil CG, Chen C, Venteicher AS. Early versus delayed postoperative radiotherapy for treatment of low-grade gliomas. Cochrane Database Syst Rev. 2020;1:CD009229.
Jiang T, Nam D, Ram Z, Poon W, Wang J, Boldbaatar D, et al. Clinical practice guidelines for the management of adult diffuse gliomas. Cancer letters. 2021;499:60-72.
Tirosh I, Suvà ML. Dissecting human gliomas by single-cell RNA sequencing. Neuro-oncology. 2018;20:37-43.
Nicholson JG, Fine HA. Diffuse Glioma Heterogeneity and Its Therapeutic Implications. Cancer Discov. 2021;11:575-90.
Li T, Chen X, Gu M, Deng A, Qian C. Identification of the subtypes of gastric cancer based on DNA methylation and the prediction of prognosis. Clin Epigenetics. 2020;12:161.
Nombela P, Miguel-López B, Blanco S. The role of mA, mC and Ψ RNA modifications in cancer: Novel therapeutic opportunities. Mol Cancer. 2021;20:18.
Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, et al. Cross-talk of four types of RNA modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer. 2021;20:29.
Huang W, Kong F, Li R, Chen X, Wang K. Emerging Roles of mA RNA Methylation Regulators in Gynecological Cancer. Front Oncol. 2022;12:827956.
Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer. 2020;20:303-22.
Xu Y, Ye S, Zhang N, Zheng S, Liu H, Zhou K, et al. The FTO/miR-181b-3p/ARL5B signaling pathway regulates cell migration and invasion in breast cancer. Cancer communications (London, England). 2020;40:484-500.
Tang B, Yang Y, Kang M, Wang Y, Wang Y, Bi Y, et al. mA demethylase ALKBH5 inhibits pancreatic cancer tumorigenesis by decreasing WIF-1 RNA methylation and mediating Wnt signaling. Mol Cancer. 2020;19:3.
Su R, Dong L, Li C, Nachtergaele S, Wunderlich M, Qing Y, et al. R-2HG Exhibits Anti-tumor Activity by Targeting FTO/mA/MYC/CEBPA Signaling. Cell. 2018;172.
Shi H, Chai P, Jia R, Fan X. Novel insight into the regulatory roles of diverse RNA modifications: Re-defining the bridge between transcription and translation. Mol Cancer. 2020;19:78.
Liu J, Chen C, Wang Y, Qian C, Wei J, Xing Y, et al. Comprehensive of N1-Methyladenosine Modifications Patterns and Immunological Characteristics in Ovarian Cancer. Front Immunol. 2021;12:746647.
Mao M, Chu Q, Lou Y, Lv P, Wang L-J. RNA N1-methyladenosine regulator-mediated methylation modification patterns and heterogeneous signatures in glioma. Front Immunol. 2022;13:948630.
Chen Y, Lin H, Miao L, He J. Role of N7-methylguanosine (m<sup>7</sup>G) in cancer. Trends Cell Biol. 2022;32:819-24.
Wang J-Z, Zhu W, Han J, Yang X, Zhou R, Lu H-C, et al. The role of the HIF-1α/ALYREF/PKM2 axis in glycolysis and tumorigenesis of bladder cancer. Cancer Communications (London, England). 2021;41:560-75.
Han X, Wang M, Zhao Y-L, Yang Y, Yang Y-G. RNA methylations in human cancers. Semin Cancer Biol. 2021;75.
Li D, Li K, Zhang W, Yang K-W, Mu D-A, Jiang G-J, et al. The m6A/m5C/m1A Regulated Gene Signature Predicts the Prognosis and Correlates With the Immune Status of Hepatocellular Carcinoma. Front Immunol. 2022;13:918140.
Dong Y, Li Y, Yao Y, Song Q. A novel defined m7G regulator signature to investigate the association between molecular characterization and clinical significance in lung adenocarcinoma. Front Oncol. 2022;12:897323.
Michalak EM, Burr ML, Bannister AJ, Dawson MA. The roles of DNA, RNA and histone methylation in ageing and cancer. Nat Rev Mol Cell Biol. 2019;20:573-89.
Xie S, Chen W, Chen K, Chang Y, Yang F, Lin A, et al. Emerging roles of RNA methylation in gastrointestinal cancers. Cancer Cell Int. 2020;20:585.
Janjigian YY, Shitara K, Moehler M, Garrido M, Salman P, Shen L, et al. First-line nivolumab plus chemotherapy versus chemotherapy alone for advanced gastric, gastro-oesophageal junction, and oesophageal adenocarcinoma (CheckMate 649): a randomised, open-label, phase 3 trial. Lancet. 2021;398:27-40.
Cloughesy TF, Mochizuki AY, Orpilla JR, Hugo W, Lee AH, Davidson TB, et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med. 2019;25:477-86.
Zhu J, Zhang T, Li J, Lin J, Liang W, Huang W, et al. Association Between Tumor Mutation Burden (TMB) and Outcomes of Cancer Patients Treated With PD-1/PD-L1 Inhibitions: A Meta-Analysis. Front Pharmacol. 2019;10:673.
Owonikoko TK, Park K, Govindan R, Ready N, Reck M, Peters S, et al. Nivolumab and Ipilimumab as Maintenance Therapy in Extensive-Disease Small-Cell Lung Cancer: CheckMate 451. J Clin Oncol. 2021;39:1349-59.
Tsao N, Brickner JR, Rodell R, Ganguly A, Wood M, Oyeniran C, et al. Aberrant RNA methylation triggers recruitment of an alkylation repair complex. Mol Cell. 2021;81.
Wu Y, Chen Z, Xie G, Zhang H, Wang Z, Zhou J, et al. RNA mA methylation regulates glycolysis of cancer cells through modulating ATP5D. Proc Natl Acad Sci U S A. 2022;119:e2119038119.
Kuang W, Jin H, Yang F, Chen X, Liu J, Li T, et al. ALKBH3-dependent mA demethylation of Aurora A mRNA inhibits ciliogenesis. Cell Discov. 2022;8:25.
Wollen KL, Hagen L, Vågbø CB, Rabe R, Iveland TS, Aas PA, et al. ALKBH3 partner ASCC3 mediates P-body formation and selective clearance of MMS-induced 1-methyladenosine and 3-methylcytosine from mRNA. J Transl Med. 2021;19:287.
You X-J, Zhang S, Chen J-J, Tang F, He J, Wang J, et al. Formation and removal of 1,N6-dimethyladenosine in mammalian transfer RNA. Nucleic Acids Res. 2022;50:9858-72.
Herrmannová A, Prilepskaja T, Wagner S, Šikrová D, Zeman J, Poncová K, et al. Adapted formaldehyde gradient cross-linking protocol implicates human eIF3d and eIF3c, k and l subunits in the 43S and 48S pre-initiation complex assembly, respectively. Nucleic Acids Res. 2020;48:1969-84.
Lamper AM, Fleming RH, Ladd KM, Lee ASY. A phosphorylation-regulated eIF3d translation switch mediates cellular adaptation to metabolic stress. Science. 2020;370:853-6.
Cieśla M, Ngoc PCT, Cordero E, Martinez ÁS, Morsing M, Muthukumar S, et al. Oncogenic translation directs spliceosome dynamics revealing an integral role for SF3A3 in breast cancer. Mol Cell. 2021;81.
Yu X, Zheng Ba, Chai R. Lentivirus-mediated knockdown of eukaryotic translation initiation factor 3 subunit D inhibits proliferation of HCT116 colon cancer cells. Biosci Rep. 2014;34:e00161.
Liu G-Z, Liu J-Z, Li X-Q, Zhang L, Li S-J, Xiao T-W, et al. Knockdown of eukaryotic translation initiation factor 3 subunit D (eIF3D) inhibits proliferation of acute myeloid leukemia cells. Mol Cell Biochem. 2018;438:191-8.
Dai N. The Diverse Functions of IMP2/IGF2BP2 in Metabolism. Trends Endocrinol Metab. 2020;31:670-9.
Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, et al. Recognition of RNA N-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat Cell Biol. 2018;20:285-95.
Yao B, Zhang Q, Yang Z, An F, Nie H, Wang H, et al. CircEZH2/miR-133b/IGF2BP2 aggravates colorectal cancer progression via enhancing the stability of mA-modified CREB1 mRNA. Mol Cancer. 2022;21:140.
Aguilo F, Li S, Balasubramaniyan N, Sancho A, Benko S, Zhang F, et al. Deposition of 5-Methylcytosine on Enhancer RNAs Enables the Coactivator Function of PGC-1α. Cell Rep. 2016;14:479-92.
Chen Y, Su H, Su Y, Zhang Y, Lin Y, Haglund F. Identification of an RNA-Binding-Protein-Based Prognostic Model for Ewing Sarcoma. Cancers (Basel). 2021;13.
Vidya E, Duchaine TF. Eukaryotic mRNA Decapping Activation. Front Genet. 2022;13:832547.
Zhang Q, Yan Y-F, Lv Q, Li Y-J, Wang R-R, Sun G-B, et al. miR-4293 upregulates lncRNA WFDC21P by suppressing mRNA-decapping enzyme 2 to promote lung carcinoma proliferation. Cell Death Dis. 2021;12:735.

No competing interests reported.

FigureS1.tif
Figure S1: A-J. The half-maximal Inhibitory concentration (IC50) of 10 common chemotherapeutic drugs (AS601245, Acadesine, Rucaparib phosphate, Veliparib dihydrochloride, Navitoclax, A-770041, A-443654, Saracatinib, AZ628, and Axitinib).
SupplementaryFiles1.xlsx

Download PDF

Version 1

posted

You are reading this latest preprint version

An m1A/m6A/m7G/m5C regulator-mediated methylation modification pattern and Landscape of immune microenvironment infiltration characterization in Lower-Grade Glioma cohorts from three continents based on machine learning

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods

2.1 Data acquisition and pre-processing of Lower-Grade Gliomas

2.2 Unsupervised Clustering of m1A/ m6A/m7G/m5C modification pattern

2.3 Enrichment analysis for distinct m1A/ m6A/m7G/m5C modification patterns

2.4 Identification and verification of MMRs risk scores

2.5 Evaluation of Tumor Microenvironment (TME) and Immune Cells Infiltration in lower-grade gliomas

2.6 Single-Cell Sequencing analysis

2.7 Machine Learning Pipelines

2.8 Drug sensitivity analysis

2.9 Molecular docking

2.10 Statistical Analysis

3. Results

3.1. Modification patterns mediated by m1A/m6A/m7G/m5C-regulators in lower-grade gliomas

3.2. Exploring the m1A/m6A/m7G/m5C modification pattern regulatory genes:

3.3. GSVA enrichment analysis

3.4. Construction of m1A/m6A/m7G/m5C Regulated Gene-Related Risk Model

3.5. Landscape of m1A/m6A/m7G/m5C Regulators in LGG

3.6. m1A/m6A/m7G/m5C Regulators as potential biomarkers for LGG

3.7. Prognostic value of m1A/m6A/m7G/m5C risk scores in individual LGG patients

3.8. Identifying the indicative role of m1A/m6A/m7G/m5C scores in TME

3.9. Investigating the best-fitting compounds on m1A/m6A/m7G/m5C scores hub genes

4. Discussion

5. Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1