An Integrated Immune-Related Bioinformatics Analysis in Glioma: Prognostic Signature’s Identification and Multi-Omics Mechanisms’ Exploration

As the traditional treatment for glioma, the most common central nervous system malignancy with poor prognosis, the efficacy of high-intensity surgery combined with radiotherapy and chemotherapy is not satisfactory. The development of individualized scientific treatment strategy urgently requires the guidance of signature with clinical predictive value. In this study, five prognosis-related differentially expressed immune-related genes (PR-DE-IRGs) (CCNA2, HMGB2, CASP3, APOBEC3C, and BMP2) highly associated with glioma were identified for a prognostic model through weighted gene co-expression network analysis, univariate Cox and lasso regression. Kaplan-Meier survival curves, receiver operating characteristic curves and other methods have shown that the model has good performance in predicting the glioma patients’ prognosis. Further combined nomogram provided better predictive performance. The signature’s guiding value in clinical treatment has also been verified by multiple analysis results. We also constructed a comprehensive competing endogenous RNA (ceRNA) regulatory network based on the protective factor BMP2 to further explore its potential role in glioma progression. Numerous immune-related biological functions and pathways were enriched in a high-risk population. Further multi-omics integrative analysis revealed a strong correlation between tumor immunosuppressive environment/IDH1 mutation and signature, suggesting that their cooperation plays an important role in glioma progression.

Introduction

Gliomas are the most common primary intracranial brain tumors, accounting for 81% of malignant brain tumors (Ostrom et al., 2014). As the most malignant and aggressive form of brain tumors, gliomas can cause significant death and morbidity (Ostrom et al., 2014; Ludwig and Kornblum, 2017). Glioma is composed of a variety of malignant cells and non-malignant cells, which can develop in the special environment of the tumor microenvironment (TME), and the tumor evolution of glioma is related to the immune changes in this microenvironment (Wang Q. et al., 2017; Tan et al., 2020). At present, the treatment of glioma is mainly surgery, radiotherapy, chemotherapy, immunotherapy and targeted therapy (Lin et al., 2017; Verburg and de Witt Hamer, 2021).

Intra-tumor heterogeneity, as the main factor affecting the therapeutic effect, has brought enormous scope for the improvement of these therapeutic approaches (Van Meir et al., 2010). It is reported that TME with immunosuppressive properties can help cancer cells evade immune detection, thus leading to cancer progression (Jackson et al., 2011). Studies reveal that TME not only plays a vital role in tumor initiation, progression, and migration, but also affects the generation of therapeutic resistance and malignancy (Huang et al., 2020). TME in human glioblastoma exhibits considerable immune cell infiltration, and the disproportion of immune cells in TME may play an essential role in gliomas (Kong et al., 2020; Qiu et al., 2020). However, due to the strong immunosuppressive microenvironment of gliomas, immunotherapy strategies exhibit a very limited effect on gliomas (Locarno et al., 2020; Xu et al., 2020). Mounting evidence shows that the isocitrate dehydrogenase (IDH) mutation is crucial for the alterations in tumor immunological microenvironment, as indicated by suppression of tumor-infiltrating lymphocytes, natural killer cells and cytotoxic T cells (Bunse et al., 2018; Ren et al., 2019). Moreover, IDH mutations cause neomorphic enzymatic activity that would result in the production of the oncometabolite 2-hydroxyglutarate (2-HG), which can then directly affect the TME (Leca et al., 2021). In gliomas, IDH mutation correlates with decreased PD1/PD-L1 expression (Buege et al., 2018; Mu et al., 2018), and specific inhibitors of mutated IDH may improve the efficacy of immunotherapy in patients with IDH mutated gliomas (Kohanbash et al., 2017). Meanwhile, inhibiting 2-HG production may enhance a host’s ability to immunotherapy response (Leca et al., 2021).

In 2011, Salmena et al. proposed a hypothesis that the crosstalk among messenger RNA (mRNA), transcribed pseudogenes and long non-coding RNA (lncRNA) based on microRNA response elements (MRE) formed a network to regulate RNA transcripts (Salmena et al., 2011). Theoretically, any transcript containing MRE can act as a potential competing endogenous RNA (ceRNA), including mRNAs, lncRNAs, pseudogene RNAs and circular RNAs (circRNAs) (Karreth and Pandolfi, 2013; Qi et al., 2015). CeRNA is reported to be involved in biological processes and plays an important role in disease pathogenesis, such as ovarian cancer (Braga et al., 2020), gastric cancer (Yang et al., 2018), and human colon adenocarcinoma (Zhang et al., 2018). Many lncRNAs play significant regulatory roles in the progression of glioma and can be used as prognostic biomarkers (Langfelder and Horvath, 2008; Liu Z. et al., 2020; Mu et al., 2020; Zhu et al., 2020).

Weighted gene co-expression network analysis (WGCNA) is a new biological method that can describe the connectivity of modules within a comprehensive network and correlate the modules with external sample traits (Langfelder and Horvath, 2008). At present, WGCNA has been successfully applied to the research of numerous cancers, such as breast cancer, non-small-cell lung cancer and ovarian cancer (Niemira et al., 2019; Yin et al., 2020; Su et al., 2021). WGCNA provides an effective way to screen genes that play an important role in tumors. This study aims to screen out prognosis-related differentially expressed immune-related genes (PR-DE-IRGs) that are highly associated with gliomas to construct a prognostic model. Not only should it have excellent prognostic performance, but also abundant clinical application value. Moreover, considering the important role of TME, mutation, cell stemness, ceRNA regulatory network in tumor progression and treatment, we also hope to explore their corresponding potential biological processes in gliomas.

Materials and Methods

Collection of Glioma Samples and Identification of Differentially Expressed Immune-Related Genes

The flowchart shows the RNA sequencing and clinical data sources used in this study (Figure 1): The Cancer Genome Atlas (TCGA; cancergenome.nih.gov/), Gene Expression Omnibus (GEO; ncbi.nlm.nih.gov/geo/) and China Glioma Genome Atlas (CGGA; cgga.org.cn/) databases. While the CGGA database shared the cohort of 1,018 glioma samples, the TCGA database provided the cohort of 703 samples (698 glioma and 5 adjacent normal tissues). The GEO database covered 4 cohorts, including the GSE108474 cohort (148 astrocytoma, 228 glioblastoma multiforme, 67 oligodendroglioma, and 28 normal brain tissues), GSE4290 cohort (26 astrocytoma, 81 glioblastoma, 50 oligodendroglioma, and 23 normal brain tissues), GSE4412 cohort (85 glioma samples) and GSE43378 cohort (50 glioma samples). In addition, the ImmPort (immport.org/home) and InnateDB (innatedb.ca/) databases provided us with a gene list of 6196 immune-related genes (IRGs).

FIGURE 1

FIGURE 1. The flowchart of the whole study.

Based on this list, we extracted the RNA sequencing data of 2365, 1879, 1879, 1638, 1879, and 1700 IRGs from the TCGA, GSE108474, GSE4290, GSE4412, GSE43378, and CGGA cohorts, respectively. To identify differentially expressed immune-related genes (DE-IRGs) between tumor and normal tissues from TCGA, we set |log2 fold change| (|log2FC|) >1 and false discovery rate (FDR) < 0.05 as the filter condition. Similarly, FDR <0.05 was used as a new filter condition to identify DE-IRGs from the GSE108474 and GSE4290 cohorts, respectively. Finally, the R package “venn” was used to visualize the overlapping process of DE-IRGs lists and IRGs lists.

Identification of Differentially Expressed Immune-Related Genes Highly Associated With Glioma Based on Weighted Gene Co-Expression Network Analysis

The RNA sequencing values of 163 common DE-IRGs were extracted from the TCGA cohort to construct the co-expression network among them. After clustering the TCGA samples for eliminating free samples, we used the function pickSoftThreshold to select the best soft power β = 4 to construct the best scale-free network. Based on the formula:

$a_{ij}=|S_{ij/\beta }$

(a_ij: adjacency matrix between gene i and gene j, S_ij: similarity matrix which is done by Pearson correlation of all gene pairs, β: softpower value), we created an adjacency matrix, and converted it to a topological overlap matrix (TOM) and a corresponding dissimilarity (1-TOM) (Yip and Horvath, 2007). Then, we took 1-TOM as the distance to cluster the genes, and built a dynamic pruning tree to identify the gene modules (Yip and Horvath, 2007). Finally, after merging similar modules with 75% similarity, we identified 3 modules. Likewise, we chose the best soft power β = 4 and β = 5 to identify two modules from the GSE108474 and GSE4290 cohorts, respectively. The turquoise modules of the GSE4290, GSE108474, and TCGA cohorts all exhibited the strongest positive correlation with tumor status. The genes of these three modules were extracted separately to obtain the common DE-IRGs highly associated with glioma.

Acquisition of Prognosis-Related Differentially Expressed Immune-Related Genes and Model Construction

Samples with complete overall survival (OS) and RNA sequencing data from 6 cohorts were extracted for subsequent analysis, respectively. To obtain corresponding PR-DE-IRGs in the TCGA, GSE4412, GSE43378 and CGGA cohorts, we performed univariate Cox analysis with a cutoff value of p < 0.05. At the same time, we also performed Kaplan-Meier survival analysis on the common PR-DE-IRGs in these four cohorts to explore the relationship between their expression and OS.

668 samples from TCGA cohort were randomly matched to the training set (n = 334) and test set (n = 334) on average. Lasso regression analysis can screen out highly relevant PR-DE-IRGs from the training set, thereby minimizing the risk of overfitting the screening features, and achieving the purpose of accurately predicting the patients’ prognosis (Lu et al., 2022). The optimal penalty parameter (λ) obtained by the minimum 10-fold cross validation was used to determine five PR-DE-IRGs and corresponding coefficients for constructing the prognostic model. We calculated each sample’s risk score from four cohorts by the formula:

$Risk\mathrm{s}core={\displaystyle \sum \left(PR-DE-IRGs expression values\times corresponding coefficient\right)}$

Verification of Model’s Predictive Ability

CCGA set, GEO’s GSE4412 and GSE43378 sets, and TGGA’s training, test and whole sets were used for the verification of model’s predictive ability. We used the median risk score of samples from each cohort as a cutoff point to classify samples into high-risk and low-risk groups. Kaplan-Meier survival analysis was used to validate the model’s ability in differentiating the glioma patients’ prognosis. We also plotted the receiver operating characteristic (ROC) curves to evaluate the accuracy of the model in predicting prognosis. Univariate and multivariate Cox regression analyses further identified the role of risk score as an independent predictor of prognosis. In addition, we also performed principal component analysis (PCA) and t-distributed random neighborhood embedding (t-SNE) to evaluate the model’s ability to discretize samples through the R package “Rtsne".

Gene Set Enrichment Analysis

To enrich for potential biological functions and pathways involved in different risk groups, we ran Gene Set Enrichment Analysis (GSEA) based on the R package “clusterProfiler”. In this process, we used “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” as reference gene sets and set nominal p value < 0.05 as the filter condition.

Tumor Immunosuppressive Environment and Immune Infiltration Type Analyses

First, the overall stromal and immune cell scores for each TCGA sample were calculated using the ESTIMATE algorithm (sourceforge.net/projects/estimateproject/). We compared their differences between different risk groups, and analyzed their correlations with risk score through Spearman. Next, based on the single-sample gene set enrichment analysis (ssGSEA) of the R packages “GSEABase” and “gsva”, we quantified the scores of 16 immune cells and 13 immune functions in each sample. After visualizing their distribution across all samples using a heatmap, similar differential expression and correlation analyses were again applied to explore their relationship with the model.

Six types of immune infiltration were identified in human tumors, which corresponded from tumor promoting to tumor suppressing respectively, namely C1 (wound healing), C2 (INF-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-b dominant) (Lin et al., 2021). The two-way ANOVA was used to explore association between risk score and different types of immune infiltration. To compare the distribution of C3, C4, and C5 subtypes between different risk groups, we ran a chi-square test.

Cell Stemness Analysis of Glioma

The one-class logistic regression machine learning algorithm (OCLR) provided targeted training for stem cells (embryonic stem cells; induced pluripotent stem cells) and their differentiated ectoderm, mesoderm and endoderm progenitor cells (Malta et al., 2018). Malta et al. calculated stemness index (mDNAsi) for each TCGA sample by OCLR, ranging from low (zero) to high (one) (Malta et al., 2018). The mDNAsi data from 564 glioma samples from this study were used for our analysis. Next, we explored the mDNAsi difference between different risk groups/clinical subgroups as well as the correlation between mDNAsi and risk score.

Mutation Analysis

After counting the somatic gene mutation data of gliomas from TCGA, we obtained the top 30 genes with the highest mutation frequency. The R package “GenVisR” was used to visualize the statistical results of these 30 genes’ mutation in different risk groups. A similar method was applied to the mutated DE-IRGs highly associated with glioma. While analyzing the association between tumor mutation burden (TMB) and model, we also explored the effect of TMB on prognosis. Based on the mutation status of IDH1, the samples were divided into wild and mutant groups. In addition to the differences in the expression of 5 PR-DE-IRGs and CD274 between the two groups, the differences in the risk score, survival probability, 16 immune cells and 13 immune functions were compared.

Construction of a Comprehensive Regulatory Network Composed of Interconnected ceRNAs

To further explore the possible mechanisms involved in glioma progression from the perspective of ceRNA, we selected BMP2 as the target gene to construct a comprehensive regulatory network composed of interconnected ceRNAs. After annotating the miRNA expression data from TCGA using mature miRNA file from mirbase database (mirbase.org/), we obtained a miRNA expression matrix for 535 samples (530 glioma and 5 adjacent normal tissues). The starBase database (starbase.sysu.edu.cn/) contains multiple target gene prediction programs (PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan). We selected miRNAs that appeared more than 2 times in all prediction programs as candidate miRNAs for BMP2. During this period, Cytoscape (v3.7.2) was used to map the co-expression network of all predicted miRNAs and BMP2. Based on the correlation analyses between these miRNAs and BMP2 expression (filter condition: correlation coefficient < −0.4 and p < 0.001), we further screened out the differentially expressed miRNAs between glioma and normal tissues (filter condition: |log2FC| >1 and p < 0.05). Ultimately, only 3 miRNAs (hsa-miR-365a-3p, hsa-let-7e-5p and hsa-miR-98-5p) showed significant survival differences in the further Kaplan-Meier survival analyses (filter condition: p < 0.05).

Next, starBase was used to predict candidate lncRNAs that may bind to hsa-miR-365a-3p, hsa-let-7e-5p and hsa-miR-98-5p, respectively. Similarly, the correlation analyses between these lncRNAs and hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p expression (filter condition: correlation coefficient < −0.2 and p < 0.001), as well as their difference analyses between glioma and normal tissues (filter condition: p < 0.05) were run separately. Next, 13/13/11 lncRNAs that might bind to hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p were used for further screening. Cytoscape was run again to map the comprehensive regulatory network composed of interconnected ceRNAs. Finally, 4/2/4 lncRNAs upstream of hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p were identified from the following three analysis steps: 1) The correlation analysis between 13/13/11 lncRNAs and hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p expression (filter condition: coefficient < −0.28 and p < 0.001); 2) The correlation analysis between 13/13/11 lncRNAs and BMP2 expression (filter condition: correlation coefficient < −0.4 and p < 0.001); 3) The differential analysis of 13/13/11 lncRNAs expression between glioma and normal tissues (filter condition: |log2FC| >1 and p < 0.05). Kaplan-Meier survival analysis was again used to explore the influence of these lncRNAs on prognosis. We highlighted these lncRNAs with yellow in the comprehensive regulatory network.

To further predict the regulatory network constructed by circRNAs, miRNAs, and BMP2, we also downloaded the circRNAs expression data of the GSE165926 set (12 gliomas and 4 normal brain tissues) and the miRNAs expression data of the GSE138764 set (33 astrocytoma and 9 normal brain tissues). Based on the TCGA and GSE138764 sets, respectively, the differentially expressed miRNAs (DE-miRNAs) between glioma and normal brain tissues were obtained (filter condition: |log2FC| >2 and p < 0.05). Prognosis related differentially expressed miRNAs (PR-DE-miRNAs) were obtained by Kaplan-Meier survival analysis based on TCGA data (filter condition: p < 0.001). We re-identified the common miRNAs (hsa-miR-129-5p and hsa-miR-381-3p) from the PR-DE-miRNAs of the TCGA set, the DE-miRNAs of the GSE138764 set, and the miRNAs upstream of BMP2 predicted by starBase as candidate miRNAs that may regulate BMP2 expression. StarBase was again used to predict circRNAs that might bind to hsa-miR-129-5p and hsa-miR-381-3p, respectively. A similar approach was used to obtain differentially expressed circRNAs (DE-circRNAs) from the GSE165926 set. We got the common circRNAs from DE-circRNAs of the GSE165926 set and the circRNAs predicted by starBase as candidate circRNAs that might bind to hsa-miR-129-5p/hsa-miR-381-3p, respectively. Finally, only the up-regulated hsa_circ_0004662 and hsa_circ_0007548 in gliomas upstream of hsa-miR-129-5p satisfy the regulation of ceRNA. Cytoscape plots the corresponding regulatory network.

Clinical Application of Prognostic Model

We obtained the immunophenotypic score (IPS) of glioma patient from the Cancer Immunology Atlas (TCIA) database (tcia.at/home). Previous studies showed that immunogenicity increases with increasing IPS score (Liu et al., 2020b). By analyzing the gene expression of the four cell types that determine immunogenicity (effector cells, immunosuppressive cells, major histocompatibility complex molecules and immunomodulators), the IPS of the sample was obtained (Liu et al., 2020b). Spearman correlation analysis was used to analyze the correlation between 4 types of IPSs and 5 PR-DE-IRGs expression/risk score. We used the violin plot to show the difference of Tumor Immune Dysfunction and Exclusion (TIDE), microsatellite instability (MSI), T cell dysfunction and exclusion score between high and low risk groups. TIDE scores were calculated online at the TIDE website (tide.dfci.harvard.edu/).

We accessed the NCI-60 database containing 60 different cancer cell lines from 9 different types of tumors through the CellMiner interface (discover.nci.nih.gov/cellminer). The Pearson correlation analysis was run to explore the relationship between 5 PR-DE-IRGs expression and 263 drugs approved by the food and drug administration or clinical trials.

Establishment and Verification of a Combined Nomogram

We used the R package “regplot” to combine the risk score with two clinical prognostic factors (grade, age) to establish a combined nomogram, which can more accurately predict the survival probability of glioma patients at 1-, 2-, and 3- years. To better prove the accuracy and effectiveness of nomogram, we constructed internal calibration curves and ROC curves based on the training, test, and whole sets.

Immune Checkpoint Inhibitors/N6-Methyladenine/Multidrug-Resistance Related Gene Expression Analysis

As the most common covalent modification of RNA at the posttranscriptional level, N6-methyladenine (m6A) mRNA modification plays a key role in gliomas through various mechanisms and associates with clinicopathological features and prognosis of gliomas patients, showing a great clinical significance (Zhang et al., 2020; Zhang et al., 2021). To explore the correlation between risk score and Immune Checkpoint Inhibitors (ICIs)/m6A related gene expression, we again used Spearman correlation analysis. In addition, we also compared the expression differences of ICIs/m6A related genes between high and low risk groups. Studies have shown that the multidrug resistance-associated protein (MRP or ABCC) expression may be related to the inherent multidrug resistance in human gliomas (Mohri et al., 2000), so we respectively analyzed the correlation between the ABCC1/ABCC3 expression and the risk score, as well as their expression differences between different risk groups.

Stratified Analyses and Immunohistochemical Staining Images Verification

We used a heatmap to show the distribution of clinical characteristics across all TCGA samples from different risk groups. To visualize the differences in risk score between different subgroups for each clinical characteristic, we drew boxplots. Kaplan-Meier survival analysis was used to evaluate the model’s predictive ability in each subgroup with different clinical characteristics. We obtained the immunohistochemical (IHC) staining results from the human protein atlas database (proteinatlas.org/) to verify the differences in protein levels of 5 modelled genes between normal tissues and glioma tissues. Finally, only 4 PR-DE-IRGs (HMGB2, CCNA2, CASP3, and APOBEC3C) IHC staining images were obtained.

Statistical Method

We used Student’s t-test or Wilcoxon signed-rank test to compare the differences between continuous variables, while the differences between categorical variables were compared by Chi-square test or Fisher’s exact test (Fan et al., 2021). Univariate cox regression analysis was used to select PR-DE-IRGs. Lasso and multiple Cox regression were used to construct prognostic model and nomogram. Univariate and multivariate Cox regression analyses were also used to identify independent prognostic factors. Spearman or Pearson correlation analyses were used to analyze the correlation among variables. The above analyses were run in the R software (version 4.0.3), Perl and SPSS Statistics 22.

Results

Extraction of Common Differentially Expressed Immune-Related Genes

We extracted 455 DE-IRGs (189 genes were down-regulated and 266 genes were up-regulated in glioma) from the TCGA cohort (Supplementary Figure S1) and 968 DE-IRGs (492 genes were up-regulated and 476 genes were down-regulated in glioma) from the GSE108474 cohort (Supplementary Figure S2). In the same way, we also identified 986 DE-IRGs from the GSE4290 cohort (Supplementary Figure S3). Finally, by Venn diagram based on the 6 cohorts, we obtained 163 common DE-IRGs (Figure 2A).

FIGURE 2

FIGURE 2. Extraction of PR-DE-IRGs. (A) Venn diagram showing the overlapping process of DE-IRGs lists from three sets and IRGs lists from three other sets. (B) Venn diagram showing the overlapping process of genes from TCGA turquoise module, GSE4290 turquoise module and GSE108474 turquoise module. (C) Heatmap reflecting the expression levels of 15 PR-DE-IRGs in TCGA samples. (D) Forest plot showing the univariate Cox regression analysis results of 15 PR-DE-IRGs in TCGA.

Identification of Differentially Expressed Immune-Related Genes Highly Associated With Glioma Based on Weighted Gene Co-Expression Network Analysis

The merged dynamic pruning tree in Figures 3A,C,E showed the results of genes and corresponding modules matched in the TCGA, GSE4290 and GSE108474 cohorts, respectively. Then we identified 3 modules from the TCGA cohort (Figure 3B), 2 modules from the GSE4290 cohort (Figure 3D), and 2 modules from the GSE108474 cohort (Figure 3F). The heatmaps (Figures 3B,D,F), reflecting the correlation between modules and tumor status, showed that the turquoise module from the GSE4290 cohort had the strongest correlation with tumor status (|correlation coefficient| = 0.64, p < 0.001, Figure 3D). Finally, we got 41 common DE-IRGs highly associated with glioma using the Venn diagram (Figure 2B).

FIGURE 3

FIGURE 3. Merged dynamic pruning trees and heatmaps of TCGA, GSE4290, and GSE108474 cohorts. (A,C,E) Dynamic pruning trees merging similar gene modules into one gene module. (B,D,F) Heatmaps showing correlation between gene modules and tumor status. The color label on the left is followed by the corresponding module. The two squares after the gene modules shows the correlation coefficients and p values between the gene modules and normal/tumor status, respectively.

Construction of Prognostic Model Based on Prognosis-Related Differentially Expressed Immune-Related Genes

668 samples from the TCGA cohort, 85 samples from the GSE4412 cohort, 50 samples from the GSE43378 cohort, and 983 samples from the CGGA cohort were used in the subsequent analyses. Supplementary Tables S1,2 showed the statistical results of these samples’ clinical characteristics. During the overlapping process of TCGA cohort’s 41 PR-DE-IRGs, GSE4412 cohort’s 17 PR-DE-IRGs, GSE43378 cohort’s 25 PR-DE-IRGs, and CGGA cohort’s 40 PR-DE-IRGs, we obtained 15 common PR-DE-IRGs (Supplementary Figure S4A). The forest plots showed the univariate Cox analyses’ results of 15 common PR-DE-IRGs from 4 cohorts (TCGA cohort: Figure 2D; CGGA cohort: Supplementary Figure S4B; GSE4412 cohort: Supplementary Figure S4C; GSE43378 Cohort: Supplementary Figure S4D). Except for BMP2 and PRKX (HR < 1), the remaining 13 genes (HR > 1) were identified as prognostic risk factors in all cohorts. This conclusion was supported by further Kaplan-Meier survival analysis results (TCGA cohort: Figure 4; CGGA cohort: Supplementary Figure S5; GSE4412 cohort: Supplementary Figure S6; GSE43378 Cohort: Supplementary Figure S7). In addition, we also used a heatmap to visualize the expression levels of these 15 PR-DE-IRGs in the TCGA cohort (Figure 2C).

FIGURE 4

FIGURE 4. Survival curves of 15 PR-DE-IRGs from TCGA cohort.

Finally, lasso regression analysis determined 5 PR-DE-IRGs and corresponding coefficients (Supplementary Table S3) to build the model according to the optimal penalty parameter (λ) (Supplementary Figures S8A,B).

Verification of Model’s Predictive Ability

With survival point plots, we observed more deceased samples in the high-risk group (Figures 5A–C and Figures 6A–C). Kaplan-Meier survival curves further showed the lower survival probability for high-risk group samples (Figures 5D–F and Figures 6D–F). To better verify the model’s predictive performance, we drew the ROC curves. In the ROC curves of three TCGA’s sets, all area under curve (AUC) values were greater than 0.85 (Figures 5G–I), showing our model’s excellent performance. This conclusion was also supported by the ROC curves of the three external cohorts (all AUC >0.7) (Figures 6G–I). The risk score was identified as a prognostic risk factor by the univariate Cox regression analysis (Figures 5J–L and Figures 6J–L). After adjusting for various clinical factors, the multivariate Cox regression analysis further determined the risk score as an independent prognostic risk factor (Figures 5M–O and Figures 6M–O). Finally, both PCA and t-SNE dot plots showed that the risk grouping separated the samples well (Supplementary Figure S9).

FIGURE 5

FIGURE 5. Validation of model’s predictive performance using the data of the training, test, and whole set from TCGA. (A–C) Risk curves and survival point plots with increasing risk score. (D–F) Survival curves. (G–I) ROC curves of OS at 1-, 2-, 3-years. AUC >0.5 is considered to have a predictive value. (J–O) Forest plots reflecting the results of univariate and multivariate cox regression analyses.

FIGURE 6

FIGURE 6. Validation of model’s predictive performance using the data of the CGGA, GSE4412, and GSE43378 cohorts. (A–C) Risk curves and survival point plots with increasing risk score. (D–F) Survival curves. (G–I) ROC curves of OS at 1-, 2-, 3-years. AUC >0.5 is considered to have a predictive value. (J–O) Forest plots reflecting the results of univariate and multivariate cox regression analyses.

Gene Set Enrichment Analysis Enrichment Analysis and Tumor Immunosuppressive Environment

In the low-risk group, we found GO analysis results (neurotransmitter secretion, neurotransmitter transport, regulation of neurotransmitter levels, regulation of postsynaptic membrane potential, and regulation of synaptic plasticity) (Figure 7B) and KEGG analysis results (calcium signaling pathway, cardiac muscle contraction, long-term potentiation, and neuroactive ligand-receptor interaction) (Figure 7D) were all closely related to the nervous system. In the high risk group, both GO analysis results (CD8-positive alpha-beta T cell activation, interferon gamma mediated signaling pathway, interferon gamma production, positive regulation of B cell activation, regulation of T cell activation, T cell activation involved in immune response, T cell mediated immunity and toll-like receptor 2 signaling pathway) (Figure 7A) and KEGG analysis results (cytokine-cytokine receptor interaction, JAK-STAT signaling pathway, natural killer cell mediated cytotoxicity, NOD-like receptor signaling pathway, p53 signaling pathway and toll-like receptor signaling pathway) (Figure 7C) contained numerous immune-related biological processes. Because GSEA exhibited a strong link between high-risk group and immunity, we further analyzed the relationship between the model and the immune microenvironment. Figures 7E,F showed a positive correlation between risk score and overall immune cell score, and higher overall immune cell score in the high-risk group. Similar results were observed in the overall stromal cell score (Figures 7G,H). After visualizing 16 immune cells and 13 immune functions scores for each sample (Figure 7I), we further analyzed their correlation with risk score, and their differences between high-risk and low-risk groups. The results showed that the remaining immune cells and functions scores were positively correlated with risk score except for natural killer (NK) and mast cells (Figure 7J). The above results were consistent with the differential analysis results presented in Figure 7K.

FIGURE 7

FIGURE 7. GSEA and tumor immunosuppressive environment analysis. (A,B) GO enrichment analysis results of high-risk group and low-risk group. (C,D) KEGG enrichment analysis results of high-risk group and low-risk group. (E,G) Correlation analysis results between risk score and overall immune score/stromal score. (F,H) Difference analysis results of overall immune score/stromal score between the high-risk and low-risk groups. (I) Heatmap visualizing 16 immune cell and 13 immune function scores for each sample. (J) Correlation analysis results between risk score and 16 immune cells/13 immune functions. (K) Difference analysis results of 16 immune cells/13 immune functions between the high-risk and low-risk groups. The symbols on the right and top of the graph represent different p values, respectively. ns: no significance; *p < 0.05; **p < 0.01; ***p < 0.001.

Immune Infiltration Type and Cell Stemness Analyses

Figure 8A showed the distribution of C3, C4, and C5 between the high-risk and low-risk groups. We observed the most C5 subtypes in the low-risk group and the most C4 subtypes in the high-risk group. The results in Figure 8B further verified the above situation. Cell stemness analysis was of great significance to the epigenetic characteristics of glioma patients. We observed a significant correlation between the mDNAsi and risk score (Figure 8C), and higher mDNAsi in the high-risk group (Figure 8D). Finally, we observed higher mDNAsi in the age >60 group (Figure 8E), deceased group (Figure 8F) and G3 group (Figure 8H). But, there was no significant difference in mDNAsi between different gender subgroups (Figure 8G).

FIGURE 8

FIGURE 8. Immune infiltration type and cell stemness analyses. (A) Table showing the distribution of immune infiltrating subtypes (C3, C4, and C5) between the different risk groups. (B) Risk score difference among the 3 immune infiltrating subtypes. (C) Correlation analysis results between mDNAsi and risk score. (D) Difference analysis results of mDNAsi between high-risk and low-risk groups. (E–H) Difference analysis results of mDNAsi between different clinical feature subgroups.

Mutation Analysis

From Figure 9A, we found that 32 of the 41 DE-IRGs highly associated with glioma were mutated. The mutation frequency of IDH1 was the highest in different risk groups (Figures 9B,C). We observed that TMB was positively correlated with risk score, which was supported by further differential analysis (Figures 9D,E). The survival analysis in Figure 9F showed that the high TMB group samples have significantly lower survival probabilities (p < 0.001). Lower risk score and higher survival probability were observed in the IDH1 mutant group (Figures 9G,H). Except for BMP2, the expression levels of other genes (CD274, APOBEC3C, CASP3, CCNA2 and HMGB) in the wild-type group were higher (Figures 9I,K–O). Figure 9J showed that NK cells infiltrated more abundantly in the mutant group, while the others infiltrated more abundantly in the wild-type group.

FIGURE 9

FIGURE 9. Mutation analysis. (A) Waterfall plot reflecting the statistical results of mutations in 32 mutated DE-IRGs highly associated with gliomas. (B,C) Waterfall plots visualizing the statistical results of these 30 genes’ mutation in different risk groups. The right panel of the waterfall plot showing the mutation frequency. The different colors at the bottom of the waterfall plot showing different mutation types and clinical characteristics types. The histogram above the waterfall plot showing the TMB statistical results for each sample. (D) Correlation analysis results between risk score and TMB. (E) TMB difference between the high-risk and low-risk groups. (F) Difference in survival probability between high and low TMB groups. (G) Risk score difference between IDH1 mutant and wild-type group. (H) Difference in survival probability between IDH1 mutant and wild-type groups. (I,K–O) Expression difference of CD274/APOBEC3C/BMP2/CASP3/CCNA2/HMGB2 between IDH1 mutant and wild-type groups, respectively. (J) Difference in 16 immune cells and 13 immune functions between the IDH1 mutant and wild-type groups.

Construction of a Comprehensive Regulatory Network Composed of Interconnected ceRNAs

After reasonable prediction, we obtained 52 upstream miRNAs that may bind to BMP2. Figure 10J showed the regulatory network between these 52 miRNAs and BMP2. Supplementary Table S4 showed all predicted miRNAs’ correlation and differential analysis results. Figures 10A–I visualized the results of the correlation, difference and Kaplan-Meier survival analyses for 3 miRNAs (hsa-miR-365a-3p, hsa-let-7e-5p, and hsa-miR-98-5p). Therefore, these three miRNAs may become the most potential miRNAs regulating BMP2 expression in glioma.

FIGURE 10

FIGURE 10. Construction of a comprehensive regulatory network composed of interconnected ceRNAs. (A,C,E) Correlation analysis results between BMP2 and hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p. (B,D,F) Difference analysis results of miR-365a-3p/let-7e-5p/miR-98-5p between the glioma and the normal groups. (G–I) Kaplan-Meier survival curves of hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p. (J) Regulatory network between 52 predicted miRNAs and BMP2. (K–M) CeRNA regulatory networks based on lncRNAs, hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p and BMP2.

Based on the starBase database, we successfully obtained 112/101/112 possible upstream lncRNAs of miR-365a-3p/let-7e-5p/miR-98-5p, respectively. Figures 10K–M showed the comprehensive regulatory network composed of interconnected ceRNAs (13 lncRNAs, hsa-miR-365a-3p and BMP2; 13 lncRNAs, hsa-let-7e-5p and BMP2; 11 lncRNAs, hsa-miR-98- 5p and BMP2). Then, 4 lncRNAs (TUG1, LINC00689, NNT-AS1 and ZEB1-AS1), 2 lncRNAs (SNHG16 and AL590666.2), and 4 lncRNAs (PXN-AS1, ZEB1-AS1, TUG1 And NNT-AS1) were screened out separately through the analysis of next three steps. Supplementary Figures S10A–J showed the results of the correlation and difference analyses for these 10 lncRNAs. Supplementary Figure S10K showed the Kaplan-Meier survival curves of these seven lncRNAs.

We obtained 2 common miRNAs from 59 PR-DE-miRNAs of TCGA set, 92 DE-miRNAs of GSE138764 set, and 52 miRNAs of starBase (Supplementary Figure S11A). Supplementary Figures S11B–E showed the results of difference and survival analyses corresponding to these two miRNAs. 11/1 common circRNAs were obtained from 269 DE-circRNAs of the GSE165926 set and 6618/3647 circRNAs upstream of hsa-miR-129-5p/hsa-miR-381-3p predicted by starBase, respectively (Supplementary Figures 11F,G). Supplementary Figure S11H visualizes the differential expression results of these 12 circRNAs. The mechanisms by which miRNAs regulate target gene expression suggest that there should be negative correlations between miRNAs and BMP2, and between miRNAs and circRNAs. Only hsa_circ_0004662 and hsa_circ_0007548, which were up-regulated in gliomas, satisfied this regulation. Finally, the regulatory network composed of hsa_circ_0004662/hsa_circ_0007548, hsa-miR-129-5p and BMP2 was shown in Supplementary Figure S11I.

Clinical Application of Prognostic Model

IPS is used to assess which patients are more inclined to respond to ICIs therapy. A higher IPS evaluation value implies a better therapeutic effect of ICIs (Liu et al., 2020b). IPS showed broad correlation with risk score, CASP3, APOBEC3C, HMGB2, and BMP2 expression (Figures 11A,B). Such results suggested that BMP2, CASP3, APOBEC3C, HMGB2 and risk score can be used to predict the efficacy of ICIs. TIDE score was negatively correlated with the therapeutic effect of ICIs (Liu et al., 2020b). Figures 11C,F showed that there was no significant difference in TIDE score and T cell dysfunction score between the high-risk and low-risk groups. But, the T cell exclusion and MSI scores in the low-risk group were higher (Figures 11D,E).

FIGURE 11

FIGURE 11. Model’s guiding value in ICIs therapy. (A) Correlation analysis results between 4 kinds of IPSs and 5 PR-DE-IRGs. (B) Correlation analysis results between 4 kinds of IPSs and risk score. (C–F) Difference in TIDE, MSI, T cell dysfunction and exclusion scores between high-risk and low-risk groups. ns: meaningless; *p < 0.05; **p < 0.01; ***p < 0.001.

The correlation analysis results of the 5 drugs most closely related to 5 PR-DE-IRGs were visualized, respectively (Figures 12A–E). These results implied the guiding value of the five PR-DE-IRGs in these drugs’ clinical application. Unexpectedly, the expression of CASP3 was positively correlated with the sensitivity of carmustine, a drug recommended by the latest National Comprehensive Cancer Network (NCCN) for the treatment of glioma (Figure 12A).

FIGURE 12

FIGURE 12. Correlation analysis between 5 PR-DE-IRGs and the 5 drugs approved by the food and drug administration (FDA). (A) CASP3. (B) APOBEC3C. (C) BMP2. (D) CCNA2. (E) HMGB2.

Finally, we built a nomogram based on age, grade and risk score to predict the glioma patients’ survival probability (Figure 13A). The internal calibration curves showed that the predicted results were basically consistent with the actual results (Figures 13B–D). The ROC curves further confirmed the excellent predictive performance of the nomogram (AUC >0.8) (Figures 13E–G).

FIGURE 13

FIGURE 13. Establishment and verification of a combined nomogram. (A) Nomogram based on age, grade, and risk group. (B–D) Internal calibration curves for verifying model’s prediction accuracy. (E–G) ROC curves for evaluating model’s prediction performance.

Immune Checkpoint Inhibitors/N6-Methyladenine/Multidrug-Resistance Related Gene Expression Analysis

Except for VTCN1, CD200, TNFSF9, HHLA2, and ADORA2A, which were negatively correlated with the risk score, the expressions of the other ICIs-related genes were all positively correlated with the risk score (Figure 14A). We also observed that half of the m6A-related genes’ expression were positively correlated with risk score, while the other half were negatively associated with the risk score (Figure 14A). The results of the above ICIs-related and m6A-related genes were verified in the further difference analysis (Figures 14B,C). We observed a positive correlation between ABCC1/ABCC3 expression and risk score as well as their higher expression in the high-risk group, further supporting the potential value of the model in drug resistance prediction (Figures 14D–G).

FIGURE 14

FIGURE 14. Correlation analysis between risk score and ICIs-related genes/m6a-related genes/multidrug resistance related genes expression, and comparison of them between different risk groups. (A,B) ICIs-related genes. (A,C) M6a-related genes. (D,E) ABCC1. (F,G) ABCC3.

Clinical Stratification Analysis Based on Prognostic Model

The heat map showed the distribution of clinical characteristics across all TCGA samples from different risk groups (Figure 15A). Subsequently, we further explored the relationship between survival status (Figure 15B)/age (Figure 15C)/gender (Figure 15D)/grade (Figure 15E) and risk score. The risk score of patients in the deceased group, age >60 group and G3 group were higher (Figures 15B,C,E). Figure 15F showed that the model still maintained the excellent performance of discriminating prognosis in all clinical characteristics subgroups. Figures 15G–J showed the IHC staining images reflecting the protein expression levels of the four PR-DE-IRGs (HMGB2, CCNA2, CASP3, and APOBEC3C) in the model, respectively. We only observed higher expression of three genes (HMGB2, CCNA2, and CASP3) in the glioma tissues (Figures 15G–I), consistent with the conclusions from the previous analysis based on the TCGA cohort. The above results further verified the model’s stability.

FIGURE 15

FIGURE 15. Clinical stratification analysis and IHC staining images verification. (A) Heatmap reflecting the clinical characteristics of each sample in TCGA. no asterisk: no significance; ***p < 0.001. (B–E) Risk score difference of between different clinical subgroups. (F) Kaplan-Meier survival curves reflecting the model’s ability to discriminate prognosis in each clinical subgroup. (G–J) IHC staining images reflecting the protein expression levels of HMGB2, CCNA2, CASP3, and APOBEC3C in glioma and normal brain tissue.

Immunohistochemical Staining Images Verification

Figures 15G–J showed the IHC staining images reflecting the protein expression levels of the four PR-DE-IRGs (HMGB2, CCNA2, CASP3 and APOBEC3C) in the model, respectively. We only observed higher expression of three genes (HMGB2, CCNA2 and CASP3) in the glioma tissues (Figures 15G–I), consistent with the conclusions from the previous analysis based on the TCGA cohort. The above results further verified the model’s stability.

Discussion

This study extracted multiple datasets from TCGA, GEO, and CGGA databases for analysis. WGCNA identified 41 DE-IRGs highly associated with gliomas. The 15 PR-DE-IRGs screened out by prognostic analysis were used in lasso regression analysis. Finally, five PR-DE-IRGs highly associated with glioma were used to construct a prognostic model. Kaplan-Meier analysis, multivariate Cox regression and ROC curves validated the excellent performance of risk score as an independent predictor in predicting prognosis. In addition to external cohorts, the data from clinical subgroups also supported this conclusion. The IHC staining images of HMGB2, CCNA2, and CASP3 demonstrated their differences in protein expression between glioma and normal brain tissues, further supporting the stability of the model. Through correlation analysis, we determined the excellent value of the model in guiding immunotherapy and chemotherapy. In addition, integrating the analysis of Muti-Omics data, we also found many results closely related to the progression and prognosis of glioma in terms of cell stemness, ceRNA regulatory axis, mutation and tumor immunosuppressive microenvironment.

Combined with previous studies, 5 genes are closely related to numerous cancers and immunity (Kwon et al., 2010; Galluzzi et al., 2016; Gan et al., 2018; Bernard et al., 2019; Constantin et al., 2021). Cyclin A2 (CCNA2) is a cyclin gene that acts as a biomarker for various tumors (Hydbring et al., 2016; Gan et al., 2018; Guo et al., 2020; Wang et al., 2020). Cancer cells, such as colorectal cancer, esophageal carcinoma, and ovarian cancer, can be inhibited from growing and progressing through the cell cycle by the silent CCNA2 (Gan et al., 2018; Guo et al., 2020; Wang et al., 2020). In particular, the study of Xi et al. showed that the invasion and metastasis of gliomas could be inhibited by reducing the expression of CCNA2 protein (Xi et al., 2019). HMGB2 belongs to the high mobility group box (HMGB) protein family, which can participate in the innate immune response of mammalian nucleic acid mediators (Rao et al., 2013; Morinaga et al., 2021). In addition, HMGB2 has been proved to be an independent influencing factor of the patient’s prognosis by its protein or mRNA level in hepatocellular carcinoma, epithelial ovarian cancer and glioblastoma multiforme (Ouellet et al., 2006; Kwon et al., 2010; Wu et al., 2013). Caspase 3 (CASP3) belongs to the caspase protease family, whose activation can cause the cleavage of many important functional proteins in cells, thereby achieving the purpose of apoptosis (Eckhart et al., 2005; Galluzzi et al., 2016). Furthermore, Bernard et al.’s study showed that the inhibition of CASP3 can reduce the development of spontaneous tumors and make cells sensitive to chemotherapeutics (Bernard et al., 2019). Apolipoprotein B mRNA-editing catalytic polypeptide-like 3C (APOBEC3C), a member of the APOBEC3 subfamily, is implicated in HIV-1 restriction, with an unclear relationship to cancer (Refsland et al., 2010; Constantin et al., 2021). As a member of the transforming growth factor-β super-family, enhanced expression of bone morphogenetic protein 2 (BMP2) can check the growth of colorectal cancer cells, enhance apoptosis, and reduce tumor development in the body (Vishnubalaji et al., 2016; Zhou et al., 2016). The conclusions of the above studies all support the accuracy of our analysis results.

A growing body of research has revealed that ceRNA regulatory networks are closely associated with the pathogenesis of glioma. The members that emerged in our ceRNA network were all closely related to tumors, and a large proportion was closely related to gliomas. The hsa-let-7e-5p was found to be a tumor suppressor to inhibit the progression of head and neck squamous cell carcinoma by targeting CCR7 expression (Wang et al., 2019) and a potential prognosis marker for rectal carcinoma with liver metastases (Chen et al., 2018). Many previous studies have also demonstrated that hsa-miR-365 functions as a tumor-suppressor in various cancers (Nie et al., 2012; Wang Y. et al., 2017). Li et al. revealed that ZEB1-AS1 interacts with hsa-miR-365a-3p and inhibits hsa-miR-365a-3p function (Li et al., 2019). Interestingly, both hsa-let-7e-5p and hsa-miR-98-5p belong to a tumor suppressor family of miRNA, the let-7 family, which is down-regulated in many tumor types and is closely related to tumorigenesis (Lu et al., 2005; Croce, 2009; Blandino et al., 2014). The expression level of let-7 is positively correlated with the malignant grade of gliomas, indicating that let-7 is likely to be a crucial inhibitory factor in the progression of gliomas (Lee et al., 2011). The abnormal expression of hsa-miR-365a-3p has been reported in many tumors, and it is also a tumor suppressor gene in gastric cancer (Guo et al., 2013; Hamada et al., 2014). Hsa-miR-129-5p, a tumor suppressor, is down-regulated in gliomas and inhibits glioma proliferation, migration and development by targeting TGIF2, WNT5A, DNMT3A, HOXC10, FNDC3B (Xu et al., 2017; Gu et al., 2018; Zeng et al., 2018; Liu et al., 2020a). The relationship between BMP2 and hsa-miR-129-5p has been reported in many diseases, such as hepatocellular carcinoma (Liu et al., 2021) and intervertebral disc degeneration (Yang and Sun, 2019). Liu et al. proved that BMP2 is the target gene of hsa-miR-129-5p and is post-transcriptional regulated by miR-129-5p (Liu et al., 2021). All these researches support our analytics results.

Previous studies have shown that lncRNAs are involved in glioma progression and prognosis (Kiang et al., 2015; Peng et al., 2018; Xi et al., 2018). The aberrant expression of taurine up-regulated gene 1 (TUG1) is associated with cell proliferation, migration, cell cycle changes, apoptosis and drug resistance in numerous tumors (Khalil et al., 2009). In glioma, TUG1 functions as a tumor suppressor by promoting cell apoptosis via activating caspase-3 and caspase-9 mediated intrinsic pathways and inhibiting Bcl-2 mediated anti-apoptotic pathways (Li et al., 2016). Additionally, small nucleolar RNA host gene 16 (SNHG16) suppressed the expression of p21, caspase-3 and caspase-9, while promoting cyclin-D1 and cyclin-B1 expression to inhibit apoptosis in glioma cells (Zhou et al., 2020). CircRNAs can recruit other RNA species, and affect the transcriptional silencing, translation and/or decay of specific mRNAs through the binding of miRNAs (Patop et al., 2019). Han et al. observed a high expression of hsa_circ_0004662 in hepatocellular carcinoma, suggesting its important role in the occurrence and development of hepatocellular carcinoma (Han et al., 2022). The crosstalk based on shared MREs in RNA transcript forms a complex network. The expression of lncRNAs and circRNAs alters the levels of miRNAs, which in turn affects the expression of their target mRNAs (Yang et al., 2014; Panda, 2018). Theoretically, the expression level of lncRNA/circRNA is negatively correlated with that of miRNA, and our results are highly consistent with this (Gawronski et al., 2018; Goodall and Wickramasinghe, 2021). Our results suggest that the down-regulation of miRNAs as suppressors in ceRNAs leads to up-regulation of the protective factor BMP2 and the corresponding lncRNAs and circRNAs. Compared with studying the isolated ceRNA axis, the analysis of a more extensive interconnected ceRNA network may be closer to physiological conditions and is conducive to more profound insight into ceRNA-mediated gene regulation. In our study, the mRNAs, miRNAs, and lncRNAs in the network showed a significant correlation with the prognosis of glioma, suggesting that miRNAs, circRNAs and lncRNAs may affect mRNAs expression in a combined way involving multiple pathways, thus affecting the progression and prognosis of the tumor.

Except for activated dendritic cells (aDCs), mast cells, and NK cells, the remaining 13 immune cells were more abundant in the high-risk group, which further supported the observation in the GSEA enrichment analysis that the high-risk group covered more immune-related biological functions and pathways. Neutrophils and macrophages play a vital role in regulating the immune response to inflammation in cancer (Galdiero et al., 2013). The presence of neutrophils at melanoma ulcer sites is strongly associated with the cell proliferation at these sites, which is associated with poor prognosis (Antonio et al., 2015). When macrophages are exposed to T helper 2 (Th2) cytokines (such as IL-4 and IL-13), they polarize into M2 macrophages and promote tumor growth, while tumor-associated macrophages (TAMs) are mainly characterized by M2 macrophages (Galdiero et al., 2013).

Hambardzumyan et al. found that TAMs interact with tumor cells, providing a conducive microenvironment that allows the tumor to escape from immune detection, thereby promoting the proliferation and metastasis of gliomas (Hambardzumyan et al., 2016). Plasmacytoid dendritic cells (pDCs) are responsible for creating an immunosuppressive microenvironment in various tumors (Vermi et al., 2011). Peritoneal pDCs infiltration may represent an immune pathogenic pathogen microenvironment and can be used to predict poor prognosis in patients undergoing therapeutic profiling for intrahepatic cholangiocarcinoma (Hu et al., 2020). The above analyses suggest that the poor prognosis of high-risk groups may be closely related to the tumor immunosuppressive environment generated by these immune cells. According to previous reports, IDH1 mutation is common in many malignant tumors, such as glioma, acute myeloid leukemia, thyroid cancer and chondroma (Murugan et al., 2010; Amary et al., 2011; Losman et al., 2013; Bai et al., 2016). Studies have shown that mutation in the IDH1 gene, especially the R132H mutation, can promote NK cell recruitment through CX3CL1/CX3CR1 chemotherapy and are associated with a better prognosis in gliomas (Ren et al., 2019). This finding suggests that NK cells can be enriched in the IDH1 mutant group, suggesting a good prognosis for glioma, which coincides with our results. Furthermore, this may explain why NK cells were negatively correlated with risk score in our study. BMP2 gene was expressed more in the IDH1 mutant group. Blood circulating NK cells express type I and type II BMP receptors, BMP-2 and BMP-6 ligands, which mediate the signaling of the BMP family members (Robson et al., 2014). The results of the study by Robson et al. showed that the inhibition of BMP signal could effectively inhibit the effect or function of NK cells, providing a new idea for immunotherapy to kill tumor cells (Robson et al., 2014). Through the above studies, we speculate that the better prognosis of patients with IDH1 mutation is due to the fact that BMP signaling pathway can stimulate the effect or function of NK cells.

As the results show, our model has excellent capabilities in predicting the prognosis of patients with glioma and guiding clinical treatment. However, the research is not perfect and still has numerous limitations. For example, the difference in protein expression of APOBEC3C is still not supported by IHC images obtained based on experiments. Moreover, this study did not combine basic experiments to verify the results of the study. Still, it only used the results of others to explain and further speculate on our results rationally. Even if limited by these shortcomings, the rich mechanism discussion results may provide new directions for follow-up research.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

XF and ML designed the research; XF prepared the figures and drafted the manuscript; XF analyzed the data; XF, YZ, LZ, JH, YF and TZ participated in the writing of the manuscript. All authors have read and approved the final manuscript.

Funding

The authors sincerely thank all the participants and the support from the Research Fund Project of Jiangxi Provincial Department of Education [180017].

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank Prof. Xine Fan from the Department of English Teaching and Research, Ministry of Humanities and Social Sciences of Jiangxi Medical College, and Huijie Liu from the Affiliated Stomatological Hospital of Nanchang University for assistance with language editing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.889629/full#supplementary-material

Supplementary Figure S1 | Difference analysis results of IRGs from the TCGA cohort. (A) Heatmap showing the expression levels of the top 100 DE-IRGs in glioma and normal tissues. (B) Volcano plot showing the expression status of all IRGs. The red/black/green dots represent up-regulated genes/genes with no significant difference/down-regulated genes, respectively.

Supplementary Figure S2 | Difference analysis results of IRGs from the GSE108474 cohort. (A) Heatmap showing the expression levels of the top 100 DE-IRGs in glioma and normal tissues. (B) Volcano plot showing the expression status of all IRGs. The red/ black/green dots represent up-regulated genes/genes with no significant difference/down-regulated genes, respectively.

Supplementary Figure S3 | Difference analysis results of IRGs from the GSE4290 cohort. (A) Heatmap showing the expression levels of the top 100 DE-IRGs in glioma and normal tissues. (B) Volcano plot showing the expression status of all IRGs. The red/black/green dots represent up-regulated genes/genes with no significant difference/down-regulated genes, respectively.

Supplementary Figure S4 | PR-DE-IRGs shared by the TCGA, CGGA, GSE4412, and GSE43378 cohorts. (A) Venn diagram showing the overlapping process of PR-DE-IRGs from 4 cohorts. (B–D) Forest plot showing the univariate Cox regression analysis results of 15 PR-DE-IRGs from CGGA, GSE4412, and GSE43378 cohorts, respectively.

Supplementary Figure S5 | Survival curves of 15 PR-DE-IRGs from the CGGA cohort.

Supplementary Figure S6 | Survival curves of 15 PR-DE-IRGs from the GSE4412 cohort.

Supplementary Figure S7 | Survival curves of 15 PR-DE-IRGs from the GSE43378 cohort.

Supplementary Figure S8 | The process of determining 5 PR-DE-IRGs and corresponding coefficients based on the optimal penalty parameter (λ) via Lasso regression. (A) The penalty coefficient was utilized to minimize the mean square error of the model. (B) LASSO regression model profile for the 5 selected PR-DE-IRGs.

Supplementary Figure S9 | PCA and tSNE dot plots based on the training, testing, whole, CGGA, GSE4412, and GSE43378 sets.

Supplementary Figure S10 | Corresponding analysis results of 10 screened lncRNAs. (A–D) Correlation and difference analysis results of TUG1, LINC00689, NNT-AS1, and ZEB1-AS1 upstream of hsa-miR-365a-3p. (E,F) Correlation and difference analysis results of SNHG16 and AL590666.2 upstream of hsa-let-7e-5p. (G–J) Correlation and difference analysis results of PXN-AS1, ZEB1-AS1, TUG1, and NNT-AS1 upstream of hsa-miR-98-5p. (K) Kaplan-Meier survival curves of TUG1, LINC00689, NNT-AS1, ZEB1-AS1, SNHG16, AL590666.2, and PXN-AS1.

Supplementary Figure S11 | Construction of the regulatory network consisting of circRNAs, miRNAs, and BMP2. (A) Obtaining common miRNAs by Venn diagram. (B,C) Differential expression results of hsa-miR-129-5p and hsa-miR-381-3p based on GSE138764 set and TCGA set, respectively. (D,E) Kaplan-Meier survival curves for hsa-miR-129-5p and hsa-miR-381-3p based on the TCGA set, respectively. (F,G) Obtaining common circRNAs from 269 DE-circRNAs of the GSE165926 set and 6618/3647 circRNAs upstream of hsa-miR-129-5p/hsa-miR-381-3p predicted by starBase, respectively. (H) Differential expression results of 12 common circRNAs based on GSE165926 set. (I) The regulatory network consisting of hsa_circ_0004662/hsa_circ_0007548, hsa-miR-129-5p and BMP2. *p < 0.05; **p < 0.01; ***p < 0.001.

References

Amary, M. F., Bacsi, K., Maggiani, F., Damato, S., Halai, D., Berisha, F., et al. (2011). IDH1 and IDH2 Mutations Are Frequent Events in central Chondrosarcoma and central and Periosteal Chondromas but Not in Other Mesenchymal Tumours. J. Pathol. 224, 334–343. doi:10.1002/path.2913

PubMed Abstract | CrossRef Full Text | Google Scholar

Antonio, N., Bønnelykke‐Behrndtz, M. L., Ward, L. C., Collin, J., Christensen, I. J., Steiniche, T., et al. (2015). The Wound Inflammatory Response Exacerbates Growth of Pre‐neoplastic Cells and Progression to Cancer. EMBO J. 34, 2219–2236. doi:10.15252/embj.201490147

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, H., Harmancı, A. S., Erson-Omay, E. Z., Li, J., Coşkun, S., Simon, M., et al. (2016). Integrated Genomic Characterization of IDH1-Mutant Glioma Malignant Progression. Nat. Genet. 48, 59–66. doi:10.1038/ng.3457

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernard, A., Chevrier, S., Beltjens, F., Dosset, M., Viltard, E., Lagrange, A., et al. (2019). Cleaved Caspase-3 Transcriptionally Regulates Angiogenesis-Promoting Chemotherapy Resistance. Cancer Res. 79, 5958–5970. doi:10.1158/0008-5472.CAN-19-0840

PubMed Abstract | CrossRef Full Text | Google Scholar

Blandino, G., Fazi, F., Donzelli, S., Kedmi, M., Sas-Chen, A., Muti, P., et al. (2014). Tumor Suppressor microRNAs: A Novel Non-coding Alliance against Cancer. FEBS Lett. 588, 2639–2652. doi:10.1016/j.febslet.2014.03.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Braga, E. A., Fridman, M. V., Moscovtsev, A. A., Filippova, E. A., Dmitriev, A. A., and Kushlinskii, N. E. (2020). LncRNAs in Ovarian Cancer Progression, Metastasis, and Main Pathways: ceRNA and Alternative Mechanisms. Int. J. Mol. Sci. 21, 8855. doi:10.3390/ijms21228855

PubMed Abstract | CrossRef Full Text | Google Scholar

Buege, M., DiPippo, A., and DiNardo, C. (2018). Evolving Treatment Strategies for Elderly Leukemia Patients with IDH Mutations. Cancers 10, 187. doi:10.3390/cancers10060187

PubMed Abstract | CrossRef Full Text | Google Scholar

Bunse, L., Pusch, S., Bunse, T., Sahm, F., Sanghvi, K., Friedrich, M., et al. (2018). Suppression of Antitumor T Cell Immunity by the Oncometabolite (R)-2-hydroxyglutarate. Nat. Med. 24, 1192–1203. doi:10.1038/s41591-018-0095-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Lin, G., Yao, Y., Chen, J., Shui, H., Yang, Q., et al. (2018). MicroRNA Hsa-Let-7e-5p as a Potential Prognosis Marker for Rectal Carcinoma with Liver Metastases. Oncol. Lett. 15, 6913–6924. doi:10.3892/ol.2018.8181

PubMed Abstract | CrossRef Full Text | Google Scholar

Constantin, D., Dubuis, G., Conde‐Rubio, M. D. C., and Widmann, C. (2021). APOBEC3C, a Nucleolar Protein Induced by Genotoxins, is Excluded from DNA Damage Sites. FEBS J. 289, 808–831. doi:10.1111/febs.16202

PubMed Abstract | CrossRef Full Text | Google Scholar

Croce, C. M. (2009). Causes and Consequences of microRNA Dysregulation in Cancer. Nat. Rev. Genet. 10, 704–714. doi:10.1038/nrg2634

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckhart, L., Ballaun, C., Uthman, A., Kittel, C., Stichenwirth, M., Buchberger, M., et al. (2005). Identification and Characterization of a Novel Mammalian Caspase with Proapoptotic Activity. J. Biol. Chem. 280, 35077–35080. doi:10.1074/jbc.C500282200

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, X., Ou, Y., Liu, H., Zhan, L., Zhu, X., Cheng, M., et al. (2021). A Ferroptosis-Related Prognostic Signature Based on Antitumor Immunity and Tumor Protein P53 Mutation Exploration for Guiding Treatment in Patients with Head and Neck Squamous Cell Carcinoma. Front. Genet. 12, 732211. doi:10.3389/fgene.2021.732211

PubMed Abstract | CrossRef Full Text | Google Scholar

Galdiero, M. R., Bonavita, E., Barajon, I., Garlanda, C., Mantovani, A., and Jaillon, S. (2013). Tumor Associated Macrophages and Neutrophils in Cancer. Immunobiology 218, 1402–1410. doi:10.1016/j.imbio.2013.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Galluzzi, L., López-Soto, A., Kumar, S., and Kroemer, G. (2016). Caspases Connect Cell-Death Signaling to Organismal Homeostasis. Immunity 44, 221–231. doi:10.1016/j.immuni.2016.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Gan, Y., Li, Y., Li, T., Shu, G., and Yin, G. (2018). CCNA2 Acts as a Novel Biomarker in Regulating the Growth and Apoptosis of Colorectal Cancer. Cancer Manag. Res. 10, 5113–5124. doi:10.2147/CMAR.S176833

PubMed Abstract | CrossRef Full Text | Google Scholar

Gawronski, A. R., Uhl, M., Zhang, Y., Lin, Y.-Y., Niknafs, Y. S., Ramnarine, V. R., et al. (2018). MechRNA: Prediction of lncRNA Mechanisms from RNA-RNA and RNA-Protein Interactions. Bioinformatics 34, 3101–3110. doi:10.1093/bioinformatics/bty208

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodall, G. J., and Wickramasinghe, V. O. (2021). RNA in Cancer. Nat. Rev. Cancer 21, 22–36. doi:10.1038/s41568-020-00306-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, X., Gong, H., Shen, L., and Gu, Q. (2018). MicroRNA-129-5p Inhibits Human Glioma Cell Proliferation and Induces Cell Cycle Arrest by Directly Targeting DNMT3A. Am. J. Transl. Res. 10, 2834–2847.

PubMed Abstract | Google Scholar

Guo, S.-L., Ye, H., Teng, Y., Wang, Y.-L., Yang, G., Li, X.-B., et al. (2013). Akt-p53-miR-365-cyclin D1/cdc25A axis Contributes to Gastric Tumorigenesis Induced by PTEN Deficiency. Nat. Commun. 4, 2544. doi:10.1038/ncomms3544

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, F., Zhang, K., Li, M., Cui, L., Liu, G., Yan, Y., et al. (2020). miR-508-3p Suppresses the Development of Ovarian Carcinoma by Targeting CCNA2 and MMP7. Int. J. Oncol. 57, 264–276. doi:10.3892/ijo.2020.5055

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamada, S., Masamune, A., Miura, S., Satoh, K., and Shimosegawa, T. (2014). MiR-365 Induces Gemcitabine Resistance in Pancreatic Cancer Cells by Targeting the Adaptor Protein SHC1 and Pro-apoptotic Regulator BAX. Cell Signal. 26, 179–185. doi:10.1016/j.cellsig.2013.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambardzumyan, D., Gutmann, D. H., and Kettenmann, H. (2016). The Role of Microglia and Macrophages in Glioma Maintenance and Progression. Nat. Neurosci. 19, 20–27. doi:10.1038/nn.4185

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, L., Wang, M., Yang, Y., Xu, H., Wei, L., and Huang, X. (2022). Detection of Prognostic Biomarkers for Hepatocellular Carcinoma through CircRNA-Associated CeRNA Analysis. J. Clin. Transl. Hepatol. 10, 80–89. doi:10.14218/JCTH.2020.00144

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Z.-Q., Zhou, Z.-J., Luo, C.-B., Xin, H.-Y., Li, J., Yu, S.-Y., et al. (2020). Peritumoral Plasmacytoid Dendritic Cells Predict a Poor Prognosis for Intrahepatic Cholangiocarcinoma after Curative Resection. Cancer Cel Int. 20, 582. doi:10.1186/s12935-020-01676-z

CrossRef Full Text | Google Scholar

Huang, S., Song, Z., Zhang, T., He, X., Huang, K., Zhang, Q., et al. (2020). Identification of Immune Cell Infiltration and Immune-Related Genes in the Tumor Microenvironment of Glioblastomas. Front. Immunol. 11, 585034. doi:10.3389/fimmu.2020.585034

PubMed Abstract | CrossRef Full Text | Google Scholar

Hydbring, P., Malumbres, M., and Sicinski, P. (2016). Non-canonical Functions of Cell Cycle Cyclins and Cyclin-dependent Kinases. Nat. Rev. Mol. Cel Biol. 17, 280–292. doi:10.1038/nrm.2016.27

CrossRef Full Text | Google Scholar

Jackson, C., Ruzevick, J., Phallen, J., Belcaid, Z., and Lim, M. (2011). Challenges in Immunotherapy Presented by the Glioblastoma Multiforme Microenvironment. Clin. Develop. Immunol. 2011, 1–20. doi:10.1155/2011/732413

CrossRef Full Text | Google Scholar

Karreth, F. A., and Pandolfi, P. P. (2013). ceRNA Cross-Talk in Cancer: When Ce-Bling Rivalries Go Awry. Cancer Discov. 3, 1113–1121. doi:10.1158/2159-8290.CD-13-0202

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalil, A. M., Guttman, M., Huarte, M., Garber, M., Raj, A., Rivea Morales, D., et al. (2009). Many Human Large Intergenic Noncoding RNAs Associate with Chromatin-Modifying Complexes and Affect Gene Expression. Proc. Natl. Acad. Sci. U.S.A. 106, 11667–11672. doi:10.1073/pnas.0904715106

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiang, K., Zhang, X.-Q., and Leung, G. (2015). Long Non-coding RNAs: The Key Players in Glioma Pathogenesis. Cancers 7, 1406–1424. doi:10.3390/cancers7030843

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohanbash, G., Carrera, D. A., Shrivastav, S., Ahn, B. J., Jahan, N., Mazor, T., et al. (2017). Isocitrate Dehydrogenase Mutations Suppress STAT1 and CD8+ T Cell Accumulation in Gliomas. J. Clin. Invest. 127, 1425–1437. doi:10.1172/JCI90644

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, Y., Feng, Z.-C., Zhang, Y.-L., Liu, X.-F., Ma, Y., Zhao, Z.-M., et al. (2020). Identification of Immune-Related Genes Contributing to the Development of Glioblastoma Using Weighted Gene Co-expression Network Analysis. Front. Immunol. 11, 1281. doi:10.3389/fimmu.2020.01281

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, J.-H., Kim, J., Park, J. Y., Hong, S. M., Park, C. W., Hong, S. J., et al. (2010). Overexpression of High-Mobility Group Box 2 is Associated with Tumor Aggressiveness and Prognosis of Hepatocellular Carcinoma. Clin. Cancer Res. 16, 5511–5521. doi:10.1158/1078-0432.CCR-10-0825

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinform. 9, 559. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Leca, J., Fortin, J., and Mak, T. W. (2021). Illuminating the Cross-Talk between Tumor Metabolism and Immunity in IDH-Mutated Cancers. Curr. Opin. Biotechnol. 68, 181–185. doi:10.1016/j.copbio.2020.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S.-T., Chu, K., Oh, H.-J., Im, W.-S., Lim, J.-Y., Kim, S.-K., et al. (2011). Let-7 microRNA Inhibits the Proliferation of Human Glioblastoma Cells. J. Neurooncol. 102, 19–24. doi:10.1007/s11060-010-0286-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Zhang, M., An, G., and Ma, Q. (2016). LncRNA TUG1 Acts as a Tumor Suppressor in Human Glioma by Promoting Cell Apoptosis. Exp. Biol. Med. (Maywood) 241, 644–649. doi:10.1177/1535370215622708

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Guan, H., Liu, Y., and Gan, X. (2019). LncRNA ZEB1-AS1 Reduces Liver Cancer Cell Proliferation by Targeting miR-365a-3p. Exp. Ther. Med. 17, 3539–3547. doi:10.3892/etm.2019.7358

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L., Cai, J., and Jiang, C. (2017). Recent Advances in Targeted Therapy for Glioma. Curr. Med. Chem. 24, 1365–1381. doi:10.2174/0929867323666161223150242

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z., Xu, Q., Miao, D., and Yu, F. (2021). An Inflammatory Response-Related Gene Signature Can Impact the Immune Status and Predict the Prognosis of Hepatocellular Carcinoma. Front. Oncol. 11, 644416. doi:10.3389/fonc.2021.644416

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Cheng, C., Jiao, J., Huang, W., Huang, J., Sun, J., et al. (2020a). MircoRNA-129-5p Suppresses the Development of Glioma by Targeting HOXC10. Pathol. Res. Pract. 216, 152868. doi:10.1016/j.prp.2020.152868

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Meng, H., Nie, S., Sun, Y., Jiang, P., Li, S., et al. (2020b). Identification of a Prognostic Signature of Epithelial Ovarian Cancer Based on Tumor Immune Microenvironment Exploration. Genomics 112, 4827–4841. doi:10.1016/j.ygeno.2020.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Wang, X., Yang, G., Zhong, C., Zhang, R., Ye, J., et al. (2020c). Construction of lncRNA‐associated ceRNA Networks to Identify Prognostic lncRNA Biomarkers for Glioblastoma. J. Cel. Biochem. 121, 3502–3515. doi:10.1002/jcb.29625

CrossRef Full Text | Google Scholar

Liu, Z., Sun, J., Wang, X., and Cao, Z. (2021). MicroRNA-129-5p Promotes Proliferation and Metastasis of Hepatocellular Carcinoma by Regulating the BMP2 Gene. Exp. Ther. Med. 21, 257. doi:10.3892/etm.2021.9688

PubMed Abstract | CrossRef Full Text | Google Scholar

Locarno, C. V., Simonelli, M., Carenza, C., Capucetti, A., Stanzani, E., Lorenzi, E., et al. (2020). Role of Myeloid Cells in the Immunosuppressive Microenvironment in Gliomas. Immunobiology 225, 151853. doi:10.1016/j.imbio.2019.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Losman, J.-A., Looper, R. E., Koivunen, P., Lee, S., Schneider, R. K., McMahon, C., et al. (2013). (R)-2-hydroxyglutarate is Sufficient to Promote Leukemogenesis and its Effects are Reversible. Science 339, 1621–1625. doi:10.1126/science.1231677

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, J., Getz, G., Miska, E. A., Alvarez-Saavedra, E., Lamb, J., Peck, D., et al. (2005). MicroRNA Expression Profiles Classify Human Cancers. Nature 435, 834–838. doi:10.1038/nature03702

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, M., Li, J., Fan, X., Xie, F., Fan, J., and Xiong, Y. (2022). Novel Immune-Related Ferroptosis Signature in Esophageal Cancer: An Informatics Exploration of Biological Processes Related to the TMEM161B-AS1/hsa-miR-27a-3p/GCH1 Regulatory Network. Front. Genet. 13, 829384. doi:10.3389/fgene.2022.829384

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludwig, K., and Kornblum, H. I. (2017). Molecular Markers in Glioma. J. Neurooncol. 134, 505–512. doi:10.1007/s11060-017-2379-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 173, 338. doi:10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohri, M., Nitta, H., and Yamashita, J. (2000). Expression of Multidrug Resistance-Associated Protein (MRP) in Human Gliomas. J. Neurooncol. 49, 105–115. doi:10.1023/a:1026528926482

PubMed Abstract | CrossRef Full Text | Google Scholar

Morinaga, H., Muta, Y., Tanaka, T., Tanabe, M., Hamaguchi, Y., and Yanase, T. (2021). High-mobility Group Box 2 Protein is Essential for the Early Phase of Adipogenesis. Biochem. Biophys. Res. Commun. 557, 97–103. doi:10.1016/j.bbrc.2021.03.149

PubMed Abstract | CrossRef Full Text | Google Scholar

Mu, L., Long, Y., Yang, C., Jin, L., Tao, H., Ge, H., et al. (2018). The IDH1 Mutation-Induced Oncometabolite, 2-Hydroxyglutarate, May Affect DNA Methylation and Expression of PD-L1 in Gliomas. Front. Mol. Neurosci. 11, 82. doi:10.3389/fnmol.2018.00082

PubMed Abstract | CrossRef Full Text | Google Scholar

Mu, M., Niu, W., Zhang, X., Hu, S., and Niu, C. (2020). LncRNA BCYRN1 Inhibits Glioma Tumorigenesis by Competitively Binding with miR-619-5p to Regulate CUEDC2 Expression and the PTEN/AKT/p21 Pathway. Oncogene 39, 6879–6892. doi:10.1038/s41388-020-01466-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Murugan, A. K., Bojdani, E., and Xing, M. (2010). Identification and Functional Characterization of Isocitrate Dehydrogenase 1 (IDH1) Mutations in Thyroid Cancer. Biochem. Biophys. Res. Commun. 393, 555–559. doi:10.1016/j.bbrc.2010.02.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Nie, J., Liu, L., Zheng, W., Chen, L., Wu, X., Xu, Y., et al. (2012). microRNA-365, Down-Regulated in colon Cancer, Inhibits Cell Cycle Progression and Promotes Apoptosis of colon Cancer Cells by Probably Targeting Cyclin D1 and Bcl-2. Carcinogenesis 33, 220–225. doi:10.1093/carcin/bgr245

PubMed Abstract | CrossRef Full Text | Google Scholar

Niemira, M., Collin, F., Szalkowska, A., Bielska, A., Chwialkowska, K., Reszec, J., et al. (2019). Molecular Signature of Subtypes of Non-small-cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-expression Network Analysis (WGCNA). Cancers 12, 37. doi:10.3390/cancers12010037

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostrom, Q. T., Bauchet, L., Davis, F. G., Deltour, I., Fisher, J. L., Langer, C. E., et al. (2014). The Epidemiology of Glioma in Adults: A “State of the Science” Review. Neuro Oncol. 16, 896–913. doi:10.1093/neuonc/nou087

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouellet, V., Page, C. L., Guyot, M.-C., Lussier, C., Tonin, P. N., Provencher, D. M., et al. (2006). SET Complex in Serous Epithelial Ovarian Cancer. Int. J. Cancer 119, 2119–2126. doi:10.1002/ijc.22054

PubMed Abstract | CrossRef Full Text | Google Scholar

Panda, A. C. (2018). Circular RNAs Act as miRNA Sponges. Adv. Exp. Med. Biol. 1087, 67–79. doi:10.1007/978-981-13-1426-1_6

PubMed Abstract | CrossRef Full Text | Google Scholar

Patop, I. L., Wüst, S., and Kadener, S. (2019). Past, Present, and Future of Circ RNA S. EMBO J. 38, e100836. doi:10.15252/embj.2018100836

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, Z., Liu, C., and Wu, M. (2018). New Insights into Long Noncoding RNAs and Their Roles in Glioma. Mol. Cancer 17, 61. doi:10.1186/s12943-018-0812-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, X., Zhang, D.-H., Wu, N., Xiao, J.-H., Wang, X., and Ma, W. (2015). ceRNA in Cancer: Possible Functions and Clinical Implications. J. Med. Genet. 52, 710–718. doi:10.1136/jmedgenet-2015-103334

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, H., Li, Y., Cheng, S., Li, J., He, C., and Li, J. (2020). A Prognostic Microenvironment-Related Immune Signature via ESTIMATE (PROMISE Model) Predicts Overall Survival of Patients with Glioma. Front. Oncol. 10, 580263. doi:10.3389/fonc.2020.580263

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, Y., Su, J., Yang, C., Peng, L., Feng, X., and Li, Q. (2013). Characterizations of Two Grass Carp Ctenopharyngodon Idella HMGB2 Genes and Potential Roles in Innate Immunity. Develop. Comp. Immunol. 41, 164–177. doi:10.1016/j.dci.2013.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Refsland, E. W., Stenglein, M. D., Shindo, K., Albin, J. S., Brown, W. L., and Harris, R. S. (2010). Quantitative Profiling of the Full APOBEC3 mRNA Repertoire in Lymphocytes and Tissues: Implications for HIV-1 Restriction. Nucleic Acids Res. 38, 4274–4284. doi:10.1093/nar/gkq174

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, F., Zhao, Q., Huang, L., Zheng, Y., Li, L., He, Q., et al. (2019). The R132H Mutation in IDH 1 Promotes the Recruitment of NK Cells through CX 3 CL 1/CX 3 CR 1 Chemotaxis and is Correlated with a Better Prognosis in Gliomas. Immunol. Cel Biol. 97, 457–469. doi:10.1111/imcb.12225

PubMed Abstract | CrossRef Full Text | Google Scholar

Robson, N. C., Hidalgo, L., McAlpine, T., Wei, H., Martínez, V. G., Entrena, A., et al. (2014). Optimal Effector Functions in Human Natural Killer Cells Rely upon Autocrine Bone Morphogenetic Protein Signaling. Cancer Res. 74, 5019–5031. doi:10.1158/0008-5472.CAN-13-2845

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA Hypothesis: the Rosetta Stone of a Hidden RNA Language? Cell 146, 353–358. doi:10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, R., Jin, C., Zhou, L., Cao, Y., Kuang, M., Li, L., et al. (2021). Construction of a ceRNA Network of Hub Genes Affecting Immune Infiltration in Ovarian Cancer Identified by WGCNA. BMC Cancer 21, 970. doi:10.1186/s12885-021-08711-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Y. Q., Li, Y. T., Yan, T. F., Xu, Y., Liu, B. H., Yang, J. A., et al. (2020). Six Immune Associated Genes Construct Prognostic Model Evaluate Low-Grade Glioma. Front. Immunol. 11, 606164. doi:10.3389/fimmu.2020.606164

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Meir, E. G., Hadjipanayis, C. G., Norden, A. D., Shu, H. K., Wen, P. Y., and Olson, J. J. (2010). Exciting New Advances in Neuro-Oncology: the Avenue to a Cure for Malignant Glioma. CA Cancer J. Clin. 60, 166–193. doi:10.3322/caac.20069

PubMed Abstract | CrossRef Full Text | Google Scholar

Verburg, N., and de Witt Hamer, P. C. (2021). State-of-the-art Imaging for Glioma Surgery. Neurosurg. Rev. 44, 1331–1343. doi:10.1007/s10143-020-01337-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Vermi, W., Soncini, M., Melocchi, L., Sozzani, S., and Facchetti, F. (2011). Plasmacytoid Dendritic Cells and Cancer. J. Leukoc. Biol. 90, 681–690. doi:10.1189/jlb.0411190

PubMed Abstract | CrossRef Full Text | Google Scholar

Vishnubalaji, R., Yue, S., Alfayez, M., Kassem, M., Liu, F.-F., Aldahmash, A., et al. (2016). Bone Morphogenetic Protein 2 (BMP2) Induces Growth Suppression and Enhances Chemosensitivity of Human colon Cancer Cells. Cancer Cel Int 16, 77. doi:10.1186/s12935-016-0355-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Hu, B., Hu, X., Kim, H., Squatrito, M., Scarpace, L., et al. (2017a). Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell 32, 42–56. doi:10.1016/j.ccell.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Xu, C., Wang, Y., and Zhang, X. (2017b). MicroRNA-365 Inhibits Ovarian Cancer Progression by Targeting Wnt5a. Am. J. Cancer Res. 7, 1096–1106.

PubMed Abstract | Google Scholar

Wang, S., Jin, S., Liu, M.-D., Pang, P., Wu, H., Qi, Z.-Z., et al. (2019). Hsa-let-7e-5p Inhibits the Proliferation and Metastasis of Head and Neck Squamous Cell Carcinoma Cells by Targeting Chemokine Receptor 7. J. Cancer 10, 1941–1948. doi:10.7150/jca.29536

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Fu, L., Wei, D., Wang, B., Zhang, C., Zhu, T., et al. (2020). MiR-29c-3p Suppresses the Migration, Invasion and Cell Cycle in Esophageal Carcinoma via CCNA2/p53 Axis. Front. Bioeng. Biotechnol. 8, 75. doi:10.3389/fbioe.2020.00075

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z. B., Cai, L., Lin, S. J., Xiong, Z. K., Lu, J. L., Mao, Y., et al. (2013). High-mobility Group Box 2 is Associated with Prognosis of Glioblastoma by Promoting Cell Viability, Invasion, and Chemotherapeutic Resistance. Neuro. Oncol. 15, 1264–1275. doi:10.1093/neuonc/not078

PubMed Abstract | CrossRef Full Text | Google Scholar

Xi, J., Sun, Q., Ma, L., and Kang, J. (2018). Long Non-coding RNAs in Glioma Progression. Cancer Lett. 419, 203–209. doi:10.1016/j.canlet.2018.01.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Xi, X., Chu, Y., Liu, N., Wang, Q., Yin, Z., Lu, Y., et al. (2019). Joint Bioinformatics Analysis of Underlying Potential Functions of Hsa-Let-7b-5p and Core Genes in Human Glioma. J. Transl. Med. 17, 129. doi:10.1186/s12967-019-1882-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, H., Hu, Y., and Qiu, W. (2017). Potential Mechanisms of microRNA-129-5p in Inhibiting Cell Processes Including Viability, Proliferation, Migration and Invasiveness of Glioblastoma Cells U87 through Targeting FNDC3B. Biomed. Pharmacother. 87, 405–411. doi:10.1016/j.biopha.2016.12.100

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, S., Tang, L., Li, X., Fan, F., and Liu, Z. (2020). Immunotherapy for Glioma: Current Management and Future Application. Cancer Lett. 476, 1–12. doi:10.1016/j.canlet.2020.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., and Sun, P. (2019). Downregulation of microRNA‐129‐5p Increases the Risk of Intervertebral Disc Degeneration by Promoting the Apoptosis of Nucleus Pulposus Cells via Targeting BMP2. J. Cel. Biochem. 120, 19684–19690. doi:10.1002/jcb.29274

CrossRef Full Text | Google Scholar

Yang, G., Lu, X., and Yuan, L. (2014). LncRNA: A Link between RNA and Cancer. Biochim. Biophys. Acta Gene Regul. Mech. 1839, 1097–1109. doi:10.1016/j.bbagrm.2014.08.012

CrossRef Full Text | Google Scholar

Yang, X.-Z., Cheng, T.-T., He, Q.-J., Lei, Z.-Y., Chi, J., Tang, Z., et al. (2018). LINC01133 as ceRNA Inhibits Gastric Cancer Progression by Sponging miR-106a-3p to Regulate APC Expression and the Wnt/β-Catenin Pathway. Mol. Cancer 17, 126. doi:10.1186/s12943-018-0874-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, X., Wang, P., Yang, T., Li, G., Teng, X., Huang, W., et al. (2020). Identification of Key Modules and Genes Associated with Breast Cancer Prognosis Using WGCNA and ceRNA Network Analysis. Aging 13, 2519–2538. doi:10.18632/aging.202285

PubMed Abstract | CrossRef Full Text | Google Scholar

Yip, A. M., and Horvath, S. (2007). Gene Network Interconnectedness and the Generalized Topological Overlap Measure. BMC Bioinformatics 8, 22. doi:10.1186/1471-2105-8-22

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, A., Yin, J., Li, Y., Li, R., Wang, Z., Zhou, X., et al. (2018). miR-129-5p Targets Wnt5a to Block PKC/ERK/NF-κB and JNK Pathways in Glioblastoma. Cell Death Dis. 9, 394. doi:10.1038/s41419-018-0343-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Qian, W., Wang, S., Ji, D., Wang, Q., Li, J., et al. (2018). Analysis of lncRNA-Associated ceRNA Network Reveals Potential lncRNA Biomarkers in Human Colon Adenocarcinoma. Cell. Physiol. Biochem. 49, 1778–1791. doi:10.1159/000493623

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Geng, X., Li, Q., Xu, J., Tan, Y., Xiao, M., et al. (2020). m6A Modification in RNA: Biogenesis, Functions and Roles in Gliomas. J. Exp. Clin. Cancer Res. 39, 192. doi:10.1186/s13046-020-01706-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Geng, X., Xu, J., Li, Q., Hao, L., Zeng, Z., et al. (2021). Identification and Characterization of N6‐methyladenosine Modification of circRNAs in Glioblastoma. J. Cel. Mol. Med. 25, 7204–7217. doi:10.1111/jcmm.16750

CrossRef Full Text | Google Scholar

Zhou, N., Li, Q., Lin, X., Hu, N., Liao, J.-Y., Lin, L.-B., et al. (2016). BMP2 Induces Chondrogenic Differentiation, Osteogenic Differentiation and Endochondral Ossification in Stem Cells. Cell Tissue Res. 366, 101–111. doi:10.1007/s00441-016-2403-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X.-Y., Liu, H., Ding, Z.-B., Xi, H.-P., and Wang, G.-W. (2020). lncRNA SNHG16 Exerts Oncogenic Functions in Promoting Proliferation of Glioma through Suppressing P21. Pathol. Oncol. Res. 26, 1021–1028. doi:10.1007/s12253-019-00648-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., Gu, W., and Yu, C. (2020). MATN1‐AS1 Promotes Glioma Progression by Functioning as ceRNA of miR‐200b/c/429 to Regulate CHD1 Expression. Cell Prolif. 53, e12700. doi:10.1111/cpr.12700

PubMed Abstract | CrossRef Full Text | Google Scholar