An N6-methyladenosine regulation- and mRNAsi-related prognostic index reveals the distinct immune microenvironment and immunotherapy responses in lower-grade glioma
BMC Bioinformatics volume 24, Article number: 225 (2023)
N6-methyladenosine (m6A) modification is involved in tumorigenesis and progression as well as closely correlated with stem cell differentiation and pluripotency. Moreover, tumor progression includes the acquisition of stemness characteristics and accumulating loss of differentiation phenotype. Therefore, we integrated m6A modification and stemness indicator mRNAsi to classify patients and predict prognosis for LGG.
We performed consensus clustering, weighted gene co-expression network analysis, and least absolute shrinkage and selection operator Cox regression analysis to identify an m6A regulation- and mRNAsi-related prognostic index (MRMRPI). Based on this prognostic index, we also explored the differences in immune microenvironments between high- and low-risk populations. Next, immunotherapy responses were also predicted. Moreover, single-cell RNA sequencing data was further used to verify the expression of these genes in MRMRPI. At last, the tumor-promoting and tumor-associated macrophage polarization roles of TIMP1 in LGG were validated by in vitro experiments.
Ten genes (DGCR10, CYP2E1, CSMD3, HOXB3, CABP4, AVIL, PTCRA, TIMP1, CLEC18A, and SAMD9) were identified to construct the MRMRPI, which was able to successfully classify patients into high- and low-risk group. Significant differences in prognosis, immune microenvironment, and immunotherapy responses were found between distinct groups. A nomogram integrating the MRMRPI and other prognostic factors were also developed to accurately predict prognosis. Moreover, in vitro experiments illustrated that inhibition of TIMP1 could inhibit the proliferation, migration, and invasion of LGG cells and also inhibit the polarization of tumor-associated macrophages.
These findings provide novel insights into understanding the interactions of m6A methylation regulation and tumor stemness on LGG development and contribute to guiding more precise immunotherapy strategies.
Gliomas are heterogeneous neuroepithelial neoplasms deriving from the glial cells, which are the most common and lethal malignant brain tumor  in the central nervous system (CNS) . Classically, they can be classified into grades I–IV based on their histopathological features. Diffuse grade II and III gliomas, including astrocytomas, oligodendrogliomas, and oligoastrocytomas, are considered lower-grade gliomas (LGG), and grade IV glioma is called glioblastoma multiforme (GBM) . Due to the high infiltration and aggressiveness, LGGs are stubbornly resistant to first-line therapies referring to maximum neurosurgical removal followed by adjuvant radiotherapy and chemotherapy, ultimately, these tumors inevitably progress into a higher grade or experience relapse . Nevertheless, because of the obvious intratumoral heterogeneity, the patients affected by LGG still have different biological and clinical characteristics, and their response to active therapy varies from person to person. Therefore, the median survival is showing an extreme range from 5.6 to 13.3 years . Aimed at classifying patients more accurately, the genotypic features consisting of isocitrate dehydrogenase (IDH) mutation and 1p/19q co-deletion status were integrated into the traditional classification . However, it provided valuable but insufficient and imprecise risk stratification and prognosis prediction, especially for genetically heterogeneous groups. Thus, it is urgent to uncover novel biomarkers to develop risk stratification and provide a new perspective for the personalized management of patients.
Recently, N6-methyladenosine (m6A) modification has gained increasing attention. It is the most prevalent and abundant form of modification in eukaryotic mRNA and is a dynamic reversible process regulated by a methyltransferase complex involved in binding proteins, methyltransferases, demethylases, also named “readers”, “writers”, and “erasers” . Not only does m6A modification acts a vital role in mRNA metabolism ranging from RNA stability, splicing, export, intracellular distribution and translation , but also affects multiple biological processes such as regulating cell cycle and differentiation, and maintenance of circadian rhythm . Additionally, the disorder of m6A regulators leads to weakened self-renewal capacity, developmental defects, dysregulation of cell proliferation, and cell death . What is more, m6A modification is involved in many complex diseases, especially tumorigenesis and progression . Alternatively, m6A methylation is also closely correlated with stem cell differentiation and pluripotency . Excitedly, Tathiane M. Malta and colleagues  used machine learning to perform a multi-platform comprehensive analysis of 33 tumor types in The Cancer Genome Atlas (TCGA) database to obtain stemness indices that can quantify tumor stemness. It has also been demonstrated that higher mRNAsi values are accompanied by greater tumor dedifferentiation and more active biological processes related to cancer stem cells. Considering tumor progression includes the acquisition of stemness characteristics and accumulating loss of differentiation phenotype. Therefore, we integrated m6A methylation and stemness indicator mRNAsi to conduct a comprehensive analysis.
In the present study, first, we find genes related to m6A regulation and stemness index mRNAsi through consensus clustering and weighted gene co-expression network analysis (WGCNA). Next, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to identify an m6A regulation and mRNAsi-related prognostic index (MRMRPI), which could classify patients into different prognosis groups. Subsequently, we further explored the differences in immune status and immune microenvironment (such as immune cells, and immune pathways) between high- and low-risk populations. Additionally, the immunotherapy responses were predicted based on this MRMRPI. Then, a nomogram based on the MRMRPI was established to quantitatively predict prognosis. Finally, the biological function of one gene (TIMP1) in the MRMRPI was validated by in vitro experiments. It may be valuable to give novel insights into personalized management and fighting against LGGs by combining m6A and mRNAsi for providing promising prognostic targets.
Identifying m6A methylation modification patterns of LGG
To explore whether the m6A methylation modification plays a vital role in LGG, consensus clustering was performed based on the expression of 21 regulators in TCGA and CGGA datasets. These regulators included 8 writers (CBLL1, KIAA1429, METTL3, METTL14, WTAP, ZC3H13, RBM15, RBM15B), 11 readers (ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, LRPPRC, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2), and 2 erasers (ALKBH5, FTO). It should be noted that IGF2BP1 is excluded from the CGGA analysis because of the lack of available sequencing data. In total, 481 samples in TCGA and 404 samples in CGGA with complete survival information were enrolled in the study. In TCGA cohort, patients were divided into three groups when k = 3 (Additional file 1: Fig. S1A-C), and there were significant differences in patient survival between different groups (Additional file 1: Fig. S1D, log-rank test). Meanwhile, similar results can be obtained in the CGGA cohort (Additional file 1: Fig. S1E-H).
Identifying MRGs and screening modules related to mRNAsi by WGCNA
As shown in Fig. 1A, 3136 MRGs were identified. To further screen MRGs related to mRNAsi, WGCNA was performed to establish a co-expression network. We select β = 4 as the soft threshold to ensure a scale-free network after outlier samples are removed. Eventually, these MRGs were clustered into six modules, each with no less than 30 genes (Fig. 1B) and the genes non-clustered were assigned to the gray module. Next, the correlations between these modules and each phenotype were calculated separately as shown in Fig. 1C. Among these modules, the brown (positive correlation) and turquoise (negative correlation) modules attracted our attention due to the highest correlation with mRNAsi. At the same time, for mRNAsi, the gene significance of these two modules is the most remarkable (Fig. 1D). Besides, we plotted to scatter diagrams based on the GS and MM of each gene in the two modules (Fig. 1E), the high correlations between them revealed the importance of module genes and their close correlation with mRNAsi. Therefore, a total of 260 genes were selected, including 58 brown module genes and 202 turquoise module genes for subsequent analysis, and they were significantly enriched in immune-related signaling pathways (Fig. 1F).
Development and validation of MRMRPI
To construct a robust MRMRPI, TCGA acted as the training set and CGGA as the validation set. First of all, the candidate genes in the training set were subjected to univariate Cox regression to determine the top 50 genes ranked by P value, which were further submitted to LASSO Cox regression analysis (Additional file 2: Fig. S2A). Finally, 10 prognostic genes were chosen for the risk-scoring system, and their univariate Cox results were visualized using a forest plot (Additional file 2: Fig. S2B). Consistently, the K–M curves with the median expression level of these genes as the thresholds also demonstrated their prognostic significance (Additional file 3: Fig. S3A). And they also showed the same prognostic trend in the CGGA data set (Additional file 3: Fig. S3B). The calculating formula of this risk-scoring system was equal to a line combination of the expression value of these genes and the optimal coefficients (Additional file 2: Fig. S2C) derived from LASSO regression. Based on this system, patients ranked by their risk score (RS) were distinguished into high- and low-risk groups according to the median RS value. Besides, the patients with high risk showed a significantly shortened survival time (Additional file 2: Fig. S2D). Simultaneously, aimed at validating the robustness of MRMRPI, the same RS formula and similar procedures were conducted in the validation set. 404 patients were classified into different groups, as expected, the one with low risk showed a significantly prolonged survival time (Additional file 2: Fig. S2E). Moreover, the ROC curves also revealed the good predictive performance of this MRMRPI in TCGA (Additional file 2: Fig. S2F) and CGGA (Additional file 2: Fig. S2G) cohorts. Finally, the DCA has also demonstrated the superiority of the MRMRPI in training (Additional file 2: Fig. S2H) and validation (Additional file 2: Fig. S2I) cohorts, which would make the clinical application more convincing.
Subgroup analysis of immune landscape and relevant biological processes based on the MRMRPI
The result that screened module genes were significantly enriched in immune-related signaling pathways (Fig. 1F) such as cell adhesion molecules, Toll-like receptor signaling pathway, and so on, which indicated that the tumor immune microenvironment probably acted an important role in the LGG progress. Therefore, we further investigate the differences in the immune microenvironment between the high and low-risk groups of the MRMRPI. First of all, GSEA was performed in TCGA cohort, and results showed that there were 41 immune-related gene sets (Additional file 4: Fig. S4A) enriched in the high-risk subgroup but no one in the low-risk subgroup. The top 10 terms ranked by FDR were further displayed in Additional file 4: Fig. S4B.
Afterward, we employed the ESTIMATE algorithm to calculate immune and stromal scores, and the MCP-counter algorithm to estimate the abundance of infiltrating immune cells (T cells, CD8+ T cells, cytotoxic lymphocytes, B lineage, NK cells, Monocytic lineage, Myeloid dendritic cells, Neutrophils), stromal cells (Endothelial cells, Fibroblasts). As the results showed (Fig. 2A, E), there was an obvious positive correlation between RS and infiltrating cells as well as various infiltrating cell subpopulations except for the statistically non-significant relationship between Endothelial cells and Monocytic lineage in TCGA set (Fig. 2A). Also, the abundance of each cell subpopulation was all significantly higher infiltrated in the high-risk group compared with the low-risk group (Fig. 2B, F). Consistently, the high-risk group showed significantly higher immune and stromal scores than the low-risk group (Student’s t test, P < 0.0001) in both TCGA and CGGA cohorts, which were respectively exhibited on the box plots of Fig. 2C, D, G, H. Simultaneously, we further explored the correlation between RS and the expression of several common immune checkpoints. Pearson’s correlation analysis indicated significantly positive correlations between RS and immune checkpoints including GZMB, CD27, CTLA-4, ICOS, LAG-3, OX40, PD-1, and PD-L2. The results originated from TCGA and CGGA cohorts were visualized in Additional file 5: Fig. S5A-H and S5I-P respectively.
Next, to have a better understanding of the immune landscape and relevant biological processes between subgroups, we further annotated a series of gene signatures. The results showed that immune pathways related to immune cell recruitment, antigen processing and presentation, immune suppression, cytotoxicity, inflammation, and adaptive and innate immunity were obviously activated in the high-risk group both in TCGA (Fig. 3) and CGGA (Additional file 6: Fig. S6) cohorts. Consistently, functional annotation demonstrated that gene sets including immune activation, stromal activation, and DNA damage repair were remarkably enhanced in the high-risk group. This phenomenon was observed in both TCGA (Fig. 4A–D) and CGGA (Fig. 4E–H) cohorts.
Chemotherapeutic drug sensitivity and immunotherapy response prediction
Chemotherapy is an important strategy for postoperative treatment; thus, we analyzed 138 chemotherapeutic drugs' sensitivity, and the results with statistically different IC50s between groups are shown in Fig. 5. Additionally, previous studies demonstrated that the tumor immune microenvironment plays important roles in tumorigenesis and immunotherapy [12, 13]. Considering that MRMRPI greatly affected the immune landscape and relevant biological processes, we inferred that the MRMRPI could be used to predict the responses to immunotherapy for LGGs patients. Eventually, the results that patients in the low-risk group responded significantly better than those in the high-risk group in both TCGA (Fig. 6A) and CGGA cohorts (Fig. 6B), which confirmed our speculation.
Construction and assessment of the nomogram
To construct a nomogram that could quantitatively predict the survival probability of LGG patients, the independent prognostic risk factors were screened by univariate Cox analysis followed by multivariate Cox analysis in TCGA cohort. The RS and clinical characteristics including age, gender, grade, mRNAsi, and IDH_1p/19q status were incorporated into the above screening process. However, only the RS and age demonstrated independent prognostic values (Fig. 7A), which were finally integrated into the construction of the nomogram (Fig. 7B). The high C-index of 0.848 showed a good discrimination ability, and the calibrations presented strong coherence for the predicted and actual probabilities of 1-, 3-, and 5-year survival (Fig. 7C). Moreover, the ROC analyses also suggested an excellent predictive ability for sensitivity and specificity (Fig. 7D) with its 1-, 3-, and 5-year predicting AUC equal to 0.907, 0.893, and 0.806 respectively.
Analysis of gene expressions at the single-cell level
To further analyze the expression of the gene signature at a single cell level, we performed scRNA-seq analysis in LGG. Totally, 7 clusters of cells were visualized by the UMAP dimensionality reduction algorithm (Fig. 8A). Moreover, 5 major cell types (neoplastic cells, oligodendrocyte precursor cells (OPCs), macrophages, astrocytes, and oligodendrocytes) were identified from the eight LGG samples (Fig. 8B). As Fig. 8C illustrated, genes like CSMD3, and CABP4 were expressed in most cell types, whereas TIMP1 was mostly expressed in macrophages and astrocytes. HOBX3 is almost only expressed in astrocytes. Gene like PTCRA was expressed very low in specific cell types.
TIMP1 promotes the proliferation, migration and invasion of LGG cells and macrophage polarization
Among the ten genes that make up the MRMRPI, the biological functions of TIMP1 in LGG remained unknown. To validate the tumor-promoting role of TIMP1 in the MRMRPI, siRNA was used to inhibited the expression of TIMP1 in HS683 and SHG44 cells (Fig. 9A). The CCK-8 assays and colony formation assays showed that the inhibition of TIMP1 significantly inhibited the proliferation ability of LGG cells (Figs. 9B, C). Meanwhile, the results of wound healing and Transwell assays demonstrated that knockdown of the expression of TIMP1 remarkably inhibited the migration and invasive capacity of LGG cells (Fig. 9D–E). These results illustrated that TIMP1 played a tumor-promoting role in LGG cells, and its risk role in the MRMRPI was validated.
In both HS683 and SHG44 cell lines, immunofluorescence staining analysis revealed that the fluorescence intensity ratio of CD163 and CD68 was considerably lower in the si-TIMP1 group than in the control group (Additional file 7: Fig. S7A). The statistical analyses were showed in Fig. S7B. These results showed that knocking down TIMP1 in glioma cells inhibited tumor-associated macrophage polarization.
Mounting evidence suggested that m6A modification serves a crucial role in cancer and non-cancerous diseases, aberrant m6A RNA methylation could affect inflammation, innate immunity as well as the response to antitumor therapy, and the interaction between regulatory factors can promote or suppress the progression of tumors . As most studies concentrate on single or several regulators rather than comprehensive and overall analysis, to fill this insufficiency, we paid attention to MRGs based on distinct m6A modification patterns that demonstrated taking effect in tumors such as gastric cancer , colon cancer , hepatocellular carcinoma  as well as gliomas . On the other hand, as well known that undifferentiated malignancies are responsible for tumor recurrence and anti-tumor resistance, and contribute to disease progression and poor prognosis. Recently, Malta et al. utilized a machine-learning algorithm  to develop a stemness index mRNAsi, which describes the similarity between stem and tumor cells according to gene expression data and possesses satisfactory adaptability and quantification of stemness. Besides, it was observed that mRNAsi was tightly associated with survival in several tumors including LGGs in pan-cancer cohorts. Hence, we combined MRGs and mRNAsi in the present study to develop a robust prognostic index for LGGs. As far as we know, it is the first study integrating MRGs and mRNAsi, which provides a novel insight into a more systematic understanding of LGG progression, facilitating more effective looking for therapeutic targets and enhancing more personalized management of LGGs.
In the present study, the m6A modification patterns classified by Consensus Clustering instead of just m6A regulators attracted our attention. Subsequently, MRGs were obtained from differential expressed analyses between various patterns, which could help us more comprehensively understand the role of m6A in LGG. Next, the WGCNA algorithm was used to integrate MRGs and mRNAsi, and two module genes were identified. Finally, there were ten genes incorporated into the development of the MRMRPI determined by univariate Cox regression and LASSO regression analyses. Among these genes, DGCR10, CYP2E1, and CSMD3 were identified as protective genes with their high expression closely associated with prolonged survival. In contrast, HOXB3, CABP4, AVIL, PTCRA, TIMP1, CLEC18A, and SAMD9 were confirmed as risk-associated genes. In previous publications, DGCR10 (also known as DGCR5) acted as a tumor-suppressive factor in glioma, consistent with our findings, which upregulated significantly inhibited glioma cell proliferation, migration, and invasion, whereas promoted apoptosis . CYP2E1 genetic polymorphism affects the susceptibility of multiple tumors [20, 21], while CSMD3 mutation is significantly linked to prognosis in some cancers such as non-small cell lung cancer , esophageal squamous cell carcinoma , but their specific mechanism still needed further exploring. Regarding the risk-associated genes, previous studies demonstrated that HOXB3 promotes cell proliferation and invasion in glioblastoma , overexpressing a cytoskeleton regulator AVIL accelerate cell proliferation and migration, enables fibroblasts to transform into immortalized astrocytes to drive tumorigenesis of glioblastoma , but their roles in LGGs is still ambiguous. TIMP1 showed a positive correlation with glioma malignancy  and its expression level was suggested as an independent predictor of glioblastoma survival . However, there was still a lack of experiments to validate the biological function of TIMP1. Herein, we performed in vitro experiment and the results showed that TIMP1 promotes the proliferation, migration and invasion of LGG cells. Moreover, the scRNA-seq analysis showed that TIMP1 is mostly expressed in macrophages, which maybe regulate the biological function of tumor-associated macrophages in the microenvironment of LGG. More importantly, in vitro experiments, we further confirmed the polarizing effect of TIMP1 on tumor-associated macrophages. Currently, there was almost no report about CABP4, PTCRA, CLEC18A, and SAMD9 in LGGs, which remains to be fully explored.
Then, we compared the immune status, as well as the distribution of immune scores, stromal scores, immune cells, immune pathways and relevant biological processes between high- and low-risk groups. Stromal scores in the high-risk group were higher than that in the low-risk group as expected, which was consistent with the reports that tumor stromal strongly facilitated the growth, progression, differentiation, and metastasis of tumor cells by nourishing tumor parenchyma . Stromal-related signatures activated in the high-risk group also confirmed such a viewpoint. Interestingly, the immune scores, as well as infiltrated immune cells, were significantly higher in the high-risk rather than low-risk group. This may be related to the formation of an immune-excluded microenvironment. As reported that although such a microenvironment was abundant in immune cells, they stayed in the matrix surrounding tumor cell nests instead of penetrating parenchyma [15, 29]. The matrix may penetrate the tumor itself or be limited to the tumor envelope, which makes immune cells seem to be truly inside the tumor. Besides, the results that the immune checkpoints expressed much higher in the high-risk group are in accord with this assumption, which could contribute to the state of immunosuppression and lead to a poor prognosis . Furthermore, antigen and inflammation persist within the tumor microenvironment, which will lead to T-cell exhaustion and dysfunction . This explains the activation of inflammatory and immune pathways in the high-risk group rather than the low-risk group. In addition, the immune system needs to maintain immune homeostasis, which contributes to avoiding potential tissue damage and autoimmunity as well as generates a successful immune defense . Our GSEA results presented significantly enriched immune-related GO terms in the high-risk group, reflecting an activated immune system in this population. However, the activated immune system of high-risk patients generated an adverse rather than beneficial effect on the prognosis of LGGs, which may be caused by the imbalance between suppressed and activated responses. This also revealed its microenvironment dominated by chronic inflammation and immunosuppression. Therefore, it might be a valuable immunotherapy strategy that maintains an equilibrium between amplifying and suppressing the immune response, which was supported by the results that better responses to immunotherapy of patients in the low-risk group.
Immunotherapy has ushered in a new era of cancer treatment, however, increasing reports [32, 33] suggest that cancer stem cells (CSCs) may potentially play an important role in treatment resistance and have been suggested to accelerate the progression and recurrence of tumors. It has been shown that in gliomas, CSCs evade immune clearance by activating regulatory T (Treg) cells, inactivating dendritic and natural killer cells , suppressing T cell proliferation and recruiting infiltration of type 2 macrophages (M2) , which results in local or systemic immunosuppression . Alternatively, tumor-associated macrophages (TAMs) secrete cytokines to enhance the self-renewal  of CSCs as well as to stimulate invasion and drug resistance  of CSCs. Therefore, strategies developed to target the stemness of tumor cells may identify new therapeutic opportunities for glioma treatment. We also established a quantitative prognostic prediction nomogram that combined the MRMRPI and independent clinical prognostic factors determined by univariate and multivariate Cox analyses. According to this nomogram, clinicians can classify patients into distinct risk stratification more accurately and provide more scientifically personalized management of LGG patients.
In conclusion, as far as we know, it was the first time to construct an MRMRPI in LGG that integrated MRGs and mRNAsi, classifying patients into distinct risk stratification. Importantly, patients in different risk stratifications have completely different prognoses and immune microenvironment. These findings provide novel insights into the understanding interactions of m6A methylation regulation and tumor stemness on LGG development and contribute to guiding more precise immunotherapy strategies.
Materials and methods
Data collection and processing
The LGGs mRNA-Seq data was extracted from TCGA and Chinese Glioma Genome Atlas (CGGA) databases. The normalized former by log2(x + 1) transformed, along with relevant clinical information were downloaded from the UCSC Xena website (https://xena.ucsc.edu/), while the latter and corresponding clinical data were downloaded from the CGGA official website (http://www.cgga.org.cn/index.jsp). Only samples with complete survival data and survival time greater than 30 days were selected for subsequent analysis. At the same time, the sequencing data in CGGA is also log2(x + 1) transformed. Besides, the stemness index mRNAsi corresponding to the sample in TCGA-LGG was obtained from the previous research .
Consensus clustering and differential expression analysis
To investigate whether the m6A regulation patterns are suitable for LGG, we mined the 21 regulators involving 8 writers, 11 readers, and 2 erasers if possible from a current study  to identify subtypes by Consensus Clustering analysis  in TCGA and CGGA datasets. 1000 times repetitions and Pearson correlation were conducted using the “ConsensuClusterPlus” package.
To screen m6A-related genes (MRGs), the “limma” R package was utilized in TCGA dataset to determine differential expression genes (DEGs) between different subtypes classified by the consensus clustering of m6A regulators expression . Adjusted P value < 0.01 was defined as the significant criteria for identifying DEGs, and the overlaps are considered to be MRGs.
Weighted gene co-expression network analysis (WGCNA)
The overlap of DEGs was used to establish a co-expression network [40, 41] in TCGA cohort and identify the modules most related to mRNAsi. Primarily, a gene expression similarity matrix was constructed based on Pearson correlation. Next, a proper soft-thresholding power β = 4 was selected to establish a signed weighted adjacency matrix, which was subsequently translated into a topological overlap matrix (TOM). Finally, the average linkage hierarchy was clustered with the parameter height = 0.25 and gene modules were identified. Besides, the module and clinical traits correlations were calculated by module eigengene (ME) representing the first principal component of each module. Next, gene significance (GS) the correlation between genes and clinical traits as well as module membership (MM) the correlation between module genes and gene expression profiles, which were both calculated. The high correlation between GS and MM means the importance of the genes in the module and the close correlation with the clinical trait.
To explore the underlying biological functions of the chosen module genes, the Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis [42,43,44,45,46] was completed by using the “clusterProfiler” package. The significant cut-off value was defined as P < 0.05.
Development and validation of the MRMRPI
After WGCNA filtering, the module genes related to m6A and mRNAsi were subjected to further analysis. Aimed at constructing an MRMRPI, of which the development was in TCGA set, while validation in CGGA set. First, univariate Cox regression was employed in the training set to screen the top 50 genes ranked by P value. Then, to identify the optimal genes with prognostic value and corresponding coefficients, the least absolute shrinkage and selection operator (LASSO) Cox analysis was performed just like the present studies [3, 47, 48]. Moreover, the MRMRPI was validated in CGGA dataset by using the same coefficients. Patients were classified into high- and low-risk subgroups based on the median cutoff value, which was proven by Kaplan–Meier (K–M) analysis with the log-rank test. Besides, time-dependent receiver operating characteristic (ROC) curves and decision curve analysis (DCA) were applied to evaluate the prognosis predicting performance.
Functional and pathway enrichment analyses
To investigate the immune status variations between different risk groups, we performed GSEA with program gsea-3.0.jar, and the immune-related gene ontology gene sets were regarded as the reference gene sets, which were obtained from the Molecular Signatures Database (http://software.broadinstitute.org/gsea/msigdb/). The enrichment items are consideredsignificant only when the nominal P value < 0.01 and the false discovery rate (FDR < 0.25). We also performed a GSVA (gene set variation analysis) algorithm to identify the differences in biological processes between distinct groups. The gene sets included c2.cp.kegg.v7.4.symbols downloaded from the Molecular Signatures Database and other relevant biological processes retrieved from the supplementary material in a previous study .
ESTIMATE algorithm and MCP-counter
ESTIMATE algorithm  was used to calculate the immune score and stromal score that could represent all immune cells and stromal cells respectively . The Microenvironment Cell Populations-counter (MCP-counter) method was applied to further estimate the infiltrating cell types. It is a validated method that enables the reliable quantification of the abundance of 8 immune and 2 matrix populations in the transcriptome of malignant tissues, and the abundance of these cells allows an inter-sample comparison .
Prediction of chemotherapeutic drug sensitivity and immunotherapy responses
To assess the predictive capability of MRMRPI in chemotherapeutic agent sensitivity, we calculated the half-maximal inhibitory concentrations (IC50) of 138 chemotherapeutic components by uing the "pRRophetic" package [53, 54] and compared the differences between the low- and high-risk groups. To further predict the response of immune checkpoint blockade (ICB) therapy, the ImmuCellAI platform (http://bioinfo.life.hust.edu.cn/ImmuCellAI/#!/analysis), a powerful tool  with high accuracy, was employed.
Nomogram construction and validation
To strengthen the predictive performance of MRMRPI and quantitatively predict prognosis, we developed a nomogram  in TCGA dataset, which integrates RS derived from MRMRPI and other important prognostic biomarkers. Univariate and multivariate Cox regression analyses were employed to determine independent prognostic biomarkers incorporated into the nomogram. Additionally, the concordance index(C-index), calibration curves, and ROC curves were used to evaluate the nomogram.
Single-cell RNA sequencing
According to a previous study, single-cell RNA sequencing (scRNA-seq) analysis was performed . We downloaded a glioma dataset (GSE89567) in the GEO database and extracted the single-cell data expression matrix of LGG. The R package Seurat was employed to analyze the single-cell data. First, a Seurat object to store the data matrix was created by the“CreateSeuratObject” function. Then, quality control was performed to discard features and cells that do not meet the basic standards: (1) genes detected in more than 3 cells; (2) cells with more than 200 total detected genes; (3) cells with less than 5% of mitochondrial genes. Next, the “NormalizeData” function was used to normalize the data, and the function“FindVariableGenes” was applied to identify 2000 highly variable genes. Then, the principal component analysis (PCA) was performed. Afterward, a uniform manifold approximation and projection (UMAP) algorithm was used for further visualization. Finally, the “scCATCH” R package combined with manual annotation was used to annotate the cell types and “FeaturePlot” was used to visualize expressions.
Cell culture and quantitative real-time PCR
LGG cell lines (HS683 and SHG44) were obtained from Xiangya School of Medicine, Central South University, Changsha, China. HS683 and SHG44 cells were cultured in high-glucose DEME (Gibco) supplemented with 10% fetal bovine serum. The Small interfering RNAs (siRNAs) against the TIMP1 gene were synthesized by RiboBio Corporation (Guangzhou, China). We used Lipofectamine 2000 transfection reagent (Invitrogen) for the siRNA transfection according to the manufacturer's protocol. The siRNA of TIMP1 (sense: CCACCUUAUACCAGCGUUATT, antisense: UAACGCUGGUAUAAGGUGGTT). The TRIzol lysis method was utilized to extract total RNA from cells. The Thermo Scientific RevertAid First Strand cDNA Synthesis Kit was used to synthesize cDNAs. The mRNA level of TIMP1 was detected by quantitative real-time PCR (qRT-PCR). The 2-ΔΔCt method was used to calculate the mRNA expression levels. The qRT-PCR primers were synthesized by Sangon Biotech (Shanghai, China), and the sequences were as follows: for TIMP1, the forward primer was 5ʹ-CTTCTGCAATTCCGACCTCGT-3ʹ and the reverse primer was 5ʹ-ACGCTGGTATAAGGTGGTCTG-3ʹ for GAPDH, the forward primer was 5ʹ-CATTGACCTCAACTACATGGTT-3ʹ and the reverse primer was 5ʹ-CCATTGATGACAAGCTTCCC-3ʹ.
Wound healing and Transwell assays
The wound healing and Transwell assays were performed by previously described methods .
Cell colony formation assay
After transfected with TIMP1 or control siRNAs, about 1000 cells/well were plated into 6-well plates and cultured for two weeks to allow colony formation. The colonies were fixed with 4% paraformaldehyde and stained with 0.01% crystal violet. Then we judged the cell growth ability according to the colony numbers.
Cell proliferation assay
After HS683 and SHG44 cells were transfected with TIMP1 or control siRNAs, the Cell Counting Kit-8 (CCK-8) assay (Vazyme, Nanjing, China) was conducted to monitor cell proliferation ability. HS683 and SHG44 cells (1.5 × 103 cells/well) were seeded into 96-well plates. Then, 10 μL of CCK-8 reagent was added to each well and incubated for 2 h at 37 °C. The Optical Density (450 nm) was determined on 0, 24, 48, and 72 h.
Macrophage differentiation and co-culture system
Macrophages (M0) were induced by THP-1 monocytes. Briefly, the THP-1 cells were seeded into the 6-well plate at 1 × 106 cells/ml. Then, 100 ng/ml PMA (Phorbol-12-myristate-13 acetate) was added for 48 h to obtain macrophages (M0). To establish a co-culture system, HS683 and SHG44 cells were seeded on top of the culture inserts, and macrophages (M0) were seeded in a 24-well plate. Subsequently, the macrophages (M0) were harvested for further analysis after 72 h.
Immunofluorescence staining was conducted to observe the expression level of CD68 and CD163. After discarding the culture medium, macrophages growing in a 24-well plate were fixed with 4% paraformaldehyde at room temperature for 10 min and then permeabilized by 0.3% Triton X-100 for 20 min. Next, at room temperature, 5% BSA (Bovine Serum Albumin) was used to block the unspecific binding sites for 2 h. The cells were incubated overnight with the primary CD68 (1:100; mouse; Proteintech 66,231–2-Ig) and CD163 (1:100; rabbit; Proteintech 16,646–1-AP) antibodies. Then, at room temperature, the slides were incubated in Alexa Fluor 568-conjugated donkey anti-mouse secondary antibody (1:500, Invitrogen) and Alexa Fluor 488-conjugated donkey anti-rabbit secondary antibody (1:500, Invitrogen) for 1 h. DAPI (1:500, Sigma, United States) was used to label the nuclei.
Statistical analyses were carried out using the R software (version 4.0.0) and GraphPad Prism (version 8.0). The log-rank test was employed in the Kaplan–Meier survival analysis. Pearson’s correlation analyses were conducted to calculate the correlation between the two groups. Student's t test was used in the two-group comparisons. A P value less than 0.05 was considered statistically significant.
Central nervous system
The Cancer Genome Atlas
Weighted gene co-expression network analysis
Least absolute shrinkage and selection operator
M6A regulation and mRNAsi related prognostic index
Chinese Glioma Genome Atlas
Differential expression genes
Topological overlap matrix
Kyoto encyclopedia of genes and genomes
Receiver operating characteristic
Gene set enrichment analysis
False discovery rate
Microenvironment cell populations-counter
Guichet PO, Masliantsev K, Tachon G, Petropoulos C, Godet J, Larrieu D, Milin S, Wager M, Karayan-Tapon L. Fatal correlation between YAP1 expression and glioma aggressiveness: clinical and molecular evidence. J Pathol. 2018. https://doi.org/10.1002/path.5133.
Liang T, Zhou X, Li P, You G, Wang F, Wang P, Feng E. DZIP3 is a key factor to stratify IDH1 wild-type lower-grade gliomas. Aging. 2020. https://doi.org/10.18632/aging.103817.
Yin W, Jiang X, Tan J, Xin Z, Zhou Q, Zhan C, Fu X, Wu Z, Guo Y, Jiang Z, et al. Development and validation of a tumor mutation burden-related immune prognostic model for lower-grade glioma. Front Oncol. 2020;10:1409. https://doi.org/10.3389/fonc.2020.01409.
Richardson LG, Nieman LT, Stemmer-Rachamimov AO, Zheng XS, Stafford K, Nagashima H, Miller JJ, Kiyokawa J, Ting DT, Wakimoto H, et al. IDH-mutant gliomas harbor fewer regulatory T cells in humans and mice. Oncoimmunology. 2020;9(1):1806662. https://doi.org/10.1080/2162402x.2020.1806662.
Lombardi G, Barresi V, Castellano A, Tabouret E, Pasqualetti F, Salvalaggio A, Cerretti G, Caccese M, Padovan M, Zagonel V, et al. Clinical Management of Diffuse Low-Grade Gliomas. Cancers. 2020. https://doi.org/10.3390/cancers12103008.
He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of N6-methyladenosine and its role in cancer. Mol Cancer. 2019;18(1):176. https://doi.org/10.1186/s12943-019-1109-9.
Zhou J, Wang J, Hong B, Ma K, Xie H, Li L, Zhang K, Zhou B, Cai L, Gong K. Gene signatures and prognostic values of m6A regulators in clear cell renal cell carcinoma—a retrospective study using TCGA database. Aging. 2019;11(6):1633–47. https://doi.org/10.18632/aging.101856.
Meyer KD, Jaffrey SR. Rethinking m(6)A readers, writers, and erasers. Annu Rev Cell Dev Biol. 2017;33:319–42. https://doi.org/10.1146/annurev-cellbio-100616-060758.
Su Y, Huang J, Hu J. m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gastric cancer. Front Oncol. 2019;9:1038. https://doi.org/10.3389/fonc.2019.01038.
Zhao X, Cui L. Development and validation of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res. 2019;9(10):2156–69.
Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, Kamińska B, Huelsken J, Omberg L, Gevaert O, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell. 2018;173(2):338-354.e315. https://doi.org/10.1016/j.cell.2018.03.034.
Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, Li J, Li F, Tan HB. Immune cells within the tumor microenvironment: biological functions and roles in cancer immunotherapy. Cancer Lett. 2020;470:126–33. https://doi.org/10.1016/j.canlet.2019.11.009.
Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol. 2016;27(8):1482–92. https://doi.org/10.1093/annonc/mdw168.
Pan J, Xu L, Pan H. Development and validation of an m6A RNA methylation regulator-based signature for prognostic prediction in cervical squamous cell carcinoma. Front Oncol. 2020;10:1444. https://doi.org/10.3389/fonc.2020.01444.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53. https://doi.org/10.1186/s12943-020-01170-0.
Chong W, Shang L, Liu J, Fang Z, Du F, Wu H, Liu Y, Wang Z, Chen Y, Jia S, et al. m(6)A regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer. Theranostics. 2021;11(5):2201–17. https://doi.org/10.7150/thno.52717.
Shen X, Hu B, Xu J, Qin W, Fu Y, Wang S, Dong Q, Qin L. The m6A methylation landscape stratifies hepatocellular carcinoma into 3 subtypes with distinct metabolic characteristics. Cancer Biol Med. 2020;17(4):937–52. https://doi.org/10.20892/j.issn.2095-3941.2020.0402.
Xu S, Tang L, Dai G, Luo C, Liu Z. Expression of m6A regulators correlated with immune microenvironment predicts therapeutic efficacy and prognosis in gliomas. Front Cell Dev Biol. 2020;8:594112. https://doi.org/10.3389/fcell.2020.594112.
He Z, Long J, Yang C, Gong B, Cheng M, Wang Q, Tang J. LncRNA DGCR5 plays a tumor-suppressive role in glioma via the miR-21/Smad7 and miR-23a/PTEN axes. Aging. 2020;12(20):20285–307. https://doi.org/10.18632/aging.103800.
Shen ZT, Wu XH, Li B, Shen JS, Wang Z, Li J, Zhu XX. CYP2E1 Rsa Ι/Pst Ι polymorphism and lung cancer susceptibility: a meta-analysis involving 10,947 subjects. J Cell Mol Med. 2015;19(9):2136–42. https://doi.org/10.1111/jcmm.12579.
Pellé L, Cipollini M, Tremmel R, Romei C, Figlioli G, Gemignani F, Melaiu O, De Santi C, Barone E, Elisei R, et al. Association between CYP2E1 polymorphisms and risk of differentiated thyroid carcinoma. Arch Toxicol. 2016;90(12):3099–109. https://doi.org/10.1007/s00204-016-1660-8.
La Fleur L, Falk-Sörqvist E, Smeds P, Berglund A, Sundström M, Mattsson JS, Brandén E, Koyi H, Isaksson J, Brunnström H, et al. Mutation patterns in a population-based non-small cell lung cancer cohort and prognostic impact of concomitant mutations in KRAS and TP53 or STK11. Lung Cancer (Amsterdam, Netherlands). 2019;130:50–8. https://doi.org/10.1016/j.lungcan.2019.01.003.
Deng J, Chen H, Zhou D, Zhang J, Chen Y, Liu Q, Ai D, Zhu H, Chu L, Ren W, et al. Comparative genomic analysis of esophageal squamous cell carcinoma between Asian and Caucasian patient populations. Nat Commun. 2017;8(1):1533. https://doi.org/10.1038/s41467-017-01730-x.
Xu K, Qiu C, Pei H, Mehmood MA, Wang H, Li L, Xia Q. Homeobox B3 promotes tumor cell proliferation and invasion in glioblastoma. Oncol Lett. 2018;15(3):3712–8. https://doi.org/10.3892/ol.2018.7750.
Xie Z, Janczyk P, Zhang Y, Liu A, Shi X, Singh S, Facemire L, Kubow K, Li Z, Jia Y, et al. A cytoskeleton regulator AVIL drives tumorigenesis in glioblastoma. Nat Commun. 2020;11(1):3457. https://doi.org/10.1038/s41467-020-17279-1.
Groft LL, Muzik H, Rewcastle NB, Johnston RN, Knäuper V, Lafleur MA, Forsyth PA, Edwards DR. Differential expression and localization of TIMP-1 and TIMP-4 in human gliomas. Br J Cancer. 2001;85(1):55–63. https://doi.org/10.1054/bjoc.2001.1854.
Crocker M, Ashley S, Giddings I, Petrik V, Hardcastle A, Aherne W, Pearson A, Bell BA, Zacharoulis S, Papadopoulos MC. Serum angiogenic profile of patients with glioblastoma identifies distinct tumor subtypes and shows that TIMP-1 is a prognostic factor. Neuro Oncol. 2011;13(1):99–108. https://doi.org/10.1093/neuonc/noq170.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144(5):646–74. https://doi.org/10.1016/j.cell.2011.02.013.
Joyce JA, Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science. 2015;348(6230):74–80. https://doi.org/10.1126/science.aaa6204.
Ghouzlani A, Rafii S, Karkouri M, Lakhdar A, Badou A. The promising IgSF11 immune checkpoint is highly expressed in advanced human gliomas and associates to poor prognosis. Front Oncol. 2020;10:608609. https://doi.org/10.3389/fonc.2020.608609.
Schneider AK, Chevalier MF, Derré L. The multifaceted immune regulation of bladder cancer. Nat Rev Urol. 2019;16(10):613–30. https://doi.org/10.1038/s41585-019-0226-y.
Wang C, Li Y, Jia L, Kim JK, Li J, Deng P, Zhang W, Krebsbach PH, Wang CY. CD276 expression enables squamous cell carcinoma stem cells to evade immune surveillance. Cell Stem Cell. 2021;28(9):1597-1613.e1597. https://doi.org/10.1016/j.stem.2021.04.011.
Jia L, Zhang W, Wang CY. BMI1 inhibition eliminates residual cancer stem cells after PD1 blockade and activates antitumor immunity to prevent metastasis and relapse. Cell Stem Cell. 2020;27(2):238-253.e236. https://doi.org/10.1016/j.stem.2020.06.022.
Galassi C, Musella M, Manduca N, Maccafeo E, Sistigu A. The immune privilege of cancer stem cells: a key to understanding tumor immune escape and therapy failure. Cells. 2021. https://doi.org/10.3390/cells10092361.
Guo X, Zhao Y, Yan H, Yang Y, Shen S, Dai X, Ji X, Ji F, Gong XG, Li L, et al. Single tumor-initiating cells evade immune clearance by recruiting type II macrophages. Genes Dev. 2017;31(3):247–59. https://doi.org/10.1101/gad.294348.116.
Saha D, Martuza RL, Rabkin SD. Macrophage polarization contributes to glioblastoma eradication by combination immunovirotherapy and immune checkpoint blockade. Cancer Cell. 2017;32(2):253-267.e255. https://doi.org/10.1016/j.ccell.2017.07.006.
Mortezaee K, Majidpoor J. Roles for macrophage-polarizing interleukins in cancer immunity and immunotherapy. Cell Oncol (Dordr). 2022;45(3):333–53. https://doi.org/10.1007/s13402-022-00667-8.
Taniguchi S, Elhance A, Van Duzer A, Kumar S, Leitenberger JJ, Oshimori N. Tumor-initiating cells establish an IL-33-TGF-β niche signaling loop to promote cancer progression. Science. 2020. https://doi.org/10.1126/science.aay1813.
Chi H, Jiang P, Xu K, Zhao Y, Song B, Peng G, He B, Liu X, Xia Z, Tian G. A novel anoikis-related gene signature predicts prognosis in patients with head and neck squamous cell carcinoma and reveals immune infiltration. Front Genet. 2022;13:984273. https://doi.org/10.3389/fgene.2022.984273.
Liang D, Xue J, Geng L, Zhou L, Lv B, Zeng Q, Xiong K, Zhou H, Xie D, Zhang F, et al. Cellular and molecular landscape of mammalian sinoatrial node revealed by single-cell RNA sequencing. Nat Commun. 2021;12(1):287. https://doi.org/10.1038/s41467-020-20448-x.
Wu CL, Dicks A, Steward N, Tang R, Katz DB, Choi YR, Guilak F. Single cell transcriptomic analysis of human pluripotent stem cell chondrogenesis. Nat Commun. 2021;12(1):362. https://doi.org/10.1038/s41467-020-20598-y.
Du R, Huang C, Liu K, Li X, Dong Z. Targeting AURKA in cancer: molecular mechanisms and opportunities for cancer therapy. Mol Cancer. 2021;20(1):15. https://doi.org/10.1186/s12943-020-01305-3.
Chi H, Peng G, Yang J, Zhang J, Song G, Xie X, Strohmer DF, Lai G, Zhao S, Wang R, et al. Machine learning to construct sphingolipid metabolism genes signature to characterize the immune landscape and prognosis of patients with uveal melanoma. Front Endocrinol (Lausanne). 2022;13:1056310. https://doi.org/10.3389/fendo.2022.1056310.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51. https://doi.org/10.1002/pro.3715.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587-d592. https://doi.org/10.1093/nar/gkac963.
Tang G, Yin W. Development of an immune infiltration-related prognostic scoring system based on the genomic landscape analysis of glioblastoma multiforme. Front Oncol. 2020;10:154. https://doi.org/10.3389/fonc.2020.00154.
Peng G, Chi H, Gao X, Zhang J, Song G, Xie X, Su K, Song B, Yang J, Gu T, et al. Identification and validation of neurotrophic factor-related genes signature in HNSCC to predict survival and immune landscapes. Front Genet. 2022;13:1010044. https://doi.org/10.3389/fgene.2022.1010044.
Cao R, Ma B, Wang G, Xiong Y, Tian Y, Yuan L. Characterization of hypoxia response patterns identified prognosis and immunotherapy response in bladder cancer. Mol Ther Oncolytics. 2021;22:277–93. https://doi.org/10.1016/j.omto.2021.06.011.
Zhang X, Zhao H, Shi X, Jia X, Yang Y. Identification and validation of an immune-related gene signature predictive of overall survival in colon cancer. Aging. 2020;12(24):26095–120. https://doi.org/10.18632/aging.202317.
Yang Z, Wei X, Pan Y, Xu J, Si Y, Min Z, Yu B. A new risk factor indicator for papillary thyroid cancer based on immune infiltration. Cell Death Dis. 2021;12(1):51. https://doi.org/10.1038/s41419-020-03294-z.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17(1):218. https://doi.org/10.1186/s13059-016-1070-5.
Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9(9):e107468. https://doi.org/10.1371/journal.pone.0107468.
Chi H, Xie X, Yan Y, Peng G, Strohmer DF, Lai G, Zhao S, Xia Z, Tian G. Natural killer cell-related prognosis signature characterizes immune landscape and predicts prognosis of HNSCC. Front Immunol. 2022;13:1018685. https://doi.org/10.3389/fimmu.2022.1018685.
Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, Guo AY. ImmuCellAI: a unique method for comprehensive t-cell subsets abundance prediction and its application in cancer immunotherapy. Adv Sci (Weinh). 2020;7(7):1902880. https://doi.org/10.1002/advs.201902880.
Weiser MR, Hsu M, Bauer PS, Chapman WC Jr, González IA, Chatterjee D, Lingam D, Mutch MG, Keshinro A, Shia J, et al. Clinical calculator based on molecular and clinicopathologic characteristics predicts recurrence following resection of stage I-III colon cancer. J Clin Oncol. 2021. https://doi.org/10.1200/jco.20.02553.
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, Hao Y, Stoeckius M, Smibert P, Satija R. Comprehensive integration of single-cell data. Cell. 2019;177(7):1888-1902.e1821. https://doi.org/10.1016/j.cell.2019.05.031.
Peng M, Zheng Z, Chen S, Fang L, Feng R, Zhang L, Tang Q, Liu X. Sensitization of non-small cell lung cancer cells to gefitinib and reversal of epithelial-mesenchymal transition by aloe-emodin via PI3K/Akt/TWIS1 signal blockage. Front Oncol. 2022;12:908031. https://doi.org/10.3389/fonc.2022.908031.
We sincerely acknowledge the public databases CGGA and TCGA.
This study was not funding-supported.
Ethics approval and consent to participate
Consent for publication
The authors declare no conflicts of interest in this manuscript.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The consensus clustering of m6A regulators could classify LGG patients into three groups in TCGA and CGGA glioma datasets. (A) Consensus clustering matrix of 481 samples from TCGA dataset for k = 3. (B) Relative change in area under the cumulative distribution function (CDF) curves according to different k values (TCGA). (C) Principal component analysis (PCA) based on the expression of m6A regulators showed distinct groups of glioma patients (TCGA). (D) Survival analysis of patients in different groups in TCGA cohort. (E) Consensus clustering matrix of 404 samples from the CGGA dataset for k = 3. (F) Relative change in area under the CDF curve according to different k values (CGGA). (G) PCA based on the expression of m6A regulators showed distinct groups of glioma patients (CGGA). (H) Survival analysis of patients in different groups in the CGGA cohort.
Construction and validation of the m6A regulation and mRNAsi-related prognostic index (MRMRPI). (A) The 10 genes were selected by least absolute shrinkage and selection operator (LASSO) Cox analysis in TCGA dataset. (B) Forest plot of the univariate Cox results of the 10 genes. (C) Coefficient values for each gene in the LASSO Cox analysis. Risk scores, living status, and Kaplan-Meier curves in the training (D) and validation cohorts (E). Time-dependent ROC curve analysis of the MRMRPI in the training (F) and validation (G) cohorts (H).
Kaplan-Meier curves of the 10 prognostic genes in TCGA (A) and CGGA (B) datasets.
GSEA showed the immune-related GO terms between low- and high-risk groups. (A) A total of 41 immune-related GO terms were significantly enriched in the high-risk group. (B) The visualization of the top 10 enrichments in the high-risk group.
The correlation analyses between risk score (RS) and immune checkpoints in TCGA (A-H) and CGGA dataset (I-P).
Activation of several immune pathways in the high-risk groups in the CGGA cohort. These pathways are involved in immune cell recruitment, antigen presentation and processing, innate immunity, immune suppression, cytotoxicity, inflammation, and adaptive immunity. Green font represents the gene overexpressed in the low-risk group, while red represents the gene overexpressed in the high-risk group. Statistical test: Wilcoxon. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.
TIMP1 promotes macrophage differentiation toward M2 in vitro. (A) The expression of CD68 and CD163 in macrophages treated differently detected by immunofluorescence. (B) Statistical analysis of different groups (** Represents p < 0.01).
About this article
Cite this article
Tang, G., Peng, J., Huo, L. et al. An N6-methyladenosine regulation- and mRNAsi-related prognostic index reveals the distinct immune microenvironment and immunotherapy responses in lower-grade glioma. BMC Bioinformatics 24, 225 (2023). https://doi.org/10.1186/s12859-023-05328-7
- Lower-grade glioma (LGG)
- N6-methyladenosine (m6A) regulation
- Prognostic index
- Immune microenvironment
- Immunotherapy responses