Identication of a Four-miRNA Risk Score Model and Hub Genes as a Novel Potential Prognostic Biomarker in Patients with Glioma

Background Novel biomarkers for malignant glioma are urgently needed due to the poor prognosis of this cancer. Numerous discoveries have demonstrated that microRNAs (miRNAs) participate in the occurrence and progression of glioma. The objective of this study is to explore new biomarkers related to glioma prognosis by bioinformatic analysis. Methods The miRNA, mRNA expression proles and relevant clinical information were extracted from TCGA database. Differently expressed miRNAs (DEMs) and differently expressed mRNAs/genes (DEGs) were screened, and a risk score comprising four miRNAs was constructed. We evaluated and tested the model accuracy through the ROC curve and a survival status analysis. Four-miRNA nomography were performed to quantitatively and visually evaluate the model. The potential biological functions of those network genes were performed through KEGG pathways and Gene Ontology (GO) terms analysis. Furthermore, survival analysis and the clinical information was used to verify whether the hub genes are associated with prognosis. The CGGA and GEPIA database were used to veried the accuracy of these analysis. To investigate the potential functions of those hub genes in glioma, we performed GSEA and GSVA. Results The results of univariate Cox regression showed that twenty miRNAs were signicantly associated with overall survival of glioma. Multivariate Cox regression analyses were performed to screen four miRNAs (hsa-miR-744-5p, hsa-miR-1249-3p, hsa-miR-3677-3p and hsa-miR-10b-3p). The estimated area under the receiver operating characteristic (ROC) curve (AUC) were 0.902, 0.896 and 0.907 for groups train, test and entire, respectively, indicating that the model had a signicant ability to predict survival in glioma patients. Furthermore, univariate and multivariate Cox regression were performed using relevant clinical factors associated with glioma, showed that the model could act as an independent prognostic factor. CGGA and GEPIA database were used to veried the accuracy of these hub genes. GSEA and GSVA revealed that the lower expressions of these hub genes were most closely associated with vascularization and tumor proliferation.

radiotherapy, the overall survival rates for patients continue to be dismal. Therefore, it is important to develop new reliable diagnostic and prognostic biomarkers for glioma.
MiRNAs are essential for maintaining homeostasis in health and disease by nely tuning the expression of genes [4,5]. Low expression or overexpression of some miRNAs (tumour suppressors or oncogenic) also cause disruption in delicate equilibrium which are involved in several aggressive cancers [6].
Accumulating evidence suggests that various miRNAs contribute to cancer initiation, progression, metastasis, diagnosis, prognosis and have potential applications in cancer therapy in numerous malignancies [7,8]. Several miRNAs are identi ed as being closely associated with cancers by binding oncogenes or tumour suppresser genes in tumour cells [9]. A growing body of research has suggested that some miRNAs are related to the oncogenesis, development and prognosis of glioma [10,11]. MiR-34a-5p has been found to inhibit glioma cell migration, proliferation and potentiate temozolomideinduced cytotoxicity by targeting high mobility group AT-hook 2 (HMGA2) [12]. It has been reported that miR-520e inhibits Wnt/β-catenin signaling by blinding to broblast growth factor 19 during glioma cells invasion and proliferation [13]. MiR-195, miR-15b and miR-103 have been shown to be tumors suppressors in glioma cells by downregulating SALL4 expression and suppressing glioma cell proliferation and invasion [14]. Therefore, studying the role of miRNA involved in the initiation and progression of gliomas is essential for making miRNA therapy for gliomas, and screening miRNA target genes for tumors is of vital importance during tumour research. Numerous miRNA-related bioinformatics technologies, such as microarray and RNA sequencing, chip sequencing, computational prediction, have been developed to validate miRNA targets, and miRNA-target interaction databases have been established.
In this study, we constructed a four-miRNA risk score model that could effectively predict the survival prognosis of glioma patients. Analysing the functional pathways and biological processes affected by these miRNAs and screening the targeted hub genes may help to elucidate the molecular mechanism of glioma and identify therapeutic targets.

Data acquisition and processing
The miRNA, mRNA expression pro les and relevant clinical information were extracted from The Cancer Genome Atlas (TCGA) website (https://www.cancer.gov). The clinical information of hub genes was downloaded from the CGGA database (Chinese Glioma Genome Atlas) which we downloaded RNA-seq data of 1018 Glioma samples as a validation set. We described a brief owchart to illustrate our study design in Fig. 1.

Identi cation Of Differently Expressed Mirnas And Mrnas
EdgeR package of R language 3.6.1 was used to screen differently expressed miRNAs (DEMs) and differently expressed mRNAs/genes (DEGs) with |log2FC|>1 and corrected p-values (FDR) < 0.05. We selected patients who had a survival time ≥ 30 days and conducted an analysis of mRNA integrated with miRNA expression pro les.

Establishment Of A Prognostic Mirna Risk Score Model
The R software "caret" package was used to classify patients into two groups (the train group and the test group). Univariate and multivariate Cox regression analyses were performed to establish the prognostic miRNA risk score model through the "Coxph" function and R language "direction = both" survival package [15]. By calculating the sum of each miRNA multiplied by its coe cient, a risk score comprising four miRNAs was constructed that was capable of predicting the prognosis of glioma.
Depending on the median of the risk score, the patients were randomized to two groups: low-risk group and high-risk group. We evaluated and tested the model accuracy in the three groups: test group, training group and entire group. The accuracy of a prediction model for glioma using the receiver operating characteristic (ROC) curve has been undertaken by using "survival ROC" package [16]. A survival status analysis and four-miRNA nomography were performed to quantitatively and visually evaluate the model using the R language "survival" and "rms" packages. Subsequently, univariate and multivariate Cox regression analysis were performed to evaluate the independent prognostic risk factor between the risk score and other clinical parameters including age and gender of patients, and stage of disease. The validity of the risk score and clinical variables were further veri ed by using ROC curves.

Network Between Mirnas And Degs And Their Functions
Data analysis using the largest and most recognized resources (Target Scan, miRTarBase and miRDB) to reveal overlapping target by Perl language. The associations between the miRNA and the target genes of microRNA was visualized using Cytoscape 3.6.1. The potential biological functions of those network genes were analysed using the "cluster Pro ler" and "org.Hs.eg.db" R packages for Kyoto Encyclopedia of Genes and Genomes (KEGG) database pathways and Gene Ontology (GO) terms. An adjusted p-value < 0.05 and q value < 0.05 were considered to be signi cant [17].

Identi cation Of Hub Genes Related To Prognostic Survival
A protein-protein interaction (PPI) network of the targeted genes was obtained from the STRING database (https://string-db.org/) [18], parameter con dence was set to 0.400. Subsequently, we identi ed the top 8 hub genes by using CytoHubba plug-in of Cytoscape 3.6.1. Furthermore, Kaplan-Meier survival analysis was used to verify whether the hub genes are associated with prognostic survival (log-rank test, p-values < 0.05).

Hub Genes Veri cation And Analysis
CGGA and GEPIA database were used to veri ed the accuracy of these hub genes. The clinical information of cancer histology (Fig. 8), grade (Fig. 9A) and the isocitrate dehydrogenase (IDH) mutation status (Fig. 9B) were signi cant correlation with these hub genes expression. GEPIA database analysis indicated that those hub genes were signi cantly lower expression in cancer tissues, however, the MAPK1 was not signi cant expression between the normal and cancer tissues (Fig. 10). The survival prognosis analysis result of CGGA database showed that the higher these hub genes expression, the better prognosis patients had (Fig. 11). Therefore, those genes may play important role in the glioma development. Combined with all these analyses supposed that these gene were tumor suppressor genes. Statistical analysis R version 3.6.1 was used for statistical analyses, and a p-value < 0.05 was considered to be signi cant ( * P < 0.05, ** P < 0.01).

Identi cation of DEMs and DEGs
There were 267 differently expressed DEMs, including 134 upregulated DEMs and 133 downregulated DEMs. We also screened 1601 DEGs, 759 of which were upregulated and 842 of which were downregulated, through the TCGA database.
Grouping of samples and construction of a prognostic miRNA risk score model We screened the entire group of patients (N = 478) who had more than 30 days of follow-up time and combined it with clinical information. Then, the entire group was randomly classi ed into two groups as train group (N = 240) and test group (N = 238). The survival random forest and univariate Cox regression results showed that twenty miRNAs were established to be signi cantly correlated with overall survival in the training group with p-value < 0.05 (Table 1). Ten miRNAs were subsequently performed for further analysis according to the multivariate Cox regression analysis. Moreover, Kaplan-Meier analysis revealed that these miRNAs were associated with the overall survival time of patients with p-value < 0.05 (Fig. 2). Nevertheless, there were contradictions between prognosis survival analysis and these miRNAs exhibited expression in tumors, such as hsa-miR-10b-5p, hsa-miR-21-3p, hsa-miR-135a-3p, hsa-miR-181c-5p, hsa-miR-221-3p and hsa-miR-548ba. Finally, four miRNAs (hsa-miR-744-5p, hsa-miR-1249-3p, hsa-miR-3677-3p and hsa-miR-10b-3p) were selected by using multivariate Cox regression analysis (Table 1). A prognostic model was built based on calculating the sum of each miRNA multiplied by its corresponding coe cient: miRNA risk score = (0.57 × expression of hsa-miR-3677-3p) + (0.66 × expression of hsa-miR-10b-3p) -(0.52 × expression of hsa-miR-1249-3p) -(0.79 × expression of hsa-miR-744-5p).
Veri ed the validity and quantitative of the miRNA risk score model Patients were randomized to two groups: high-risk group and low-risk group. Analysis of the Kaplan-Meier survival curves indicated that the overall survival was signi cantly higher in the low-risk group compared with that of the high-risk group in the test group (p = 7.435E − 04), entire group (p = 1.558E − 10) and train group (p = 8.603E − 09; Fig. 3A-C). The test group and entire group were tested, and the accuracy and effectiveness of the four-miRNA risk score model were veri ed. The estimated area under the receiver operating characteristic (ROC) curve (AUC) of the test group, entire group and train group were 0.896, 0.907 and 0.902, respectively, which showed that the prediction model has a better predictive ability ( Fig. 3D-F). Then, we built a survival status analysis for brief quantitative and visual analysis of the model (Fig. 3G-I). In addition, univariate and multivariate Cox regression were performed using relevant clinical factors associated with glioma, showed that the four-miRNA risk score model could act as an independent prognostic factor ( Table 2). Univariate Cox regression analysis showed that the risk score was markedly related to overall survival in patients (hazard ratio HR = 1.010, con dence interval 95% CI = 1.007-1.013, p < 0.01). Multivariate Cox regression analysis revealed that the risk score was independent factor of overall survival, considering other clinical factors (hazard ratio HR = 1.008, con dence interval 95% CI = 1.004-1.012, p < 0.01; Table 2), such as clinical grade and age. The multivariate ROC curve indicated that risk score (0.908) had a high predictive ability (Fig. 4A). Finally, we built a nomography for quantitative and visual analysis of the four-miRNA risk score model (Fig. 4B).

Network Between Mirnas And Degs And Functional Enrichment Analysis
There were 195 overlapping genes identi ed between miRNA target genes and DEGs. The network was built on those genes (Fig. 5A) and screened the eight hub genes (Fig. 5B). We predicted the potential biological functions of those genes through KEGG signaling pathway and GO enrichment analysis. The GO and KEGG pathway analysis results are listed in Tables 3 and 4, respectively. A p-value < 0.05 was used as a cut-off point. The top terms of GO annotation in biological process (BP), molecular function (MF) and cellular component (CC) are displayed. The enriched GO functions for target genes of the miRNAs were primarily associated with neural axonogenesis and synapses (Fig. 6A-F). The KEGG pathway analysis indicated that cancer-related pathways were clearly activated, including the cAMP signaling pathway, oocyte meiosis pathway and MAPK signaling pathway (Fig. 6G-I).

Survival Prognosis Analysis Of The Hub Genes
According the node degrees, eight hub genes were nally screened out by the cytoHubba plug-in, including SNAP25, MAPK1, CAMK2A, SYT1, NRXN1, CALM1, ATP2B2, and CALM3 in Table 5. Kaplan-Meier survival analysis indicated that the eight top hub genes were positively associated with the survival prognosis of glioma (Fig. 7).

Potential Functions Of These Hub Genes
To further investigate the potential functions of those hub genes in glioma, we performed GSEA and GSVA. As shown in Fig. 12, the high expression of these hub genes were most enriched in "long term potentiation", "neuroactive ligand receptor interaction" and "phosphatidylinositol signaling system" pathways. These pathways were closely associated with signal transduction and synaptic repair. The lower expressions of these hub genes were most enriched in "p53 signaling pathway", "antigen processing and presentation" and "complement and coagulation cascades" pathways. These gene sets with the highest enrichment scores were most closely associated with vascularization and tumor proliferation. Furthermore, GSVA veri ed that these gene sets were signi cantly upregulated in the lower expression (Fig. 13).

Discussion
Glioma is a highly malignant brain tumour that has a high risk of recurrence and mortality, and a poor prognosis [19]. Thus, identifying more sensitive, speci c and timely biomarkers for diagnosis, as well as investigating the mechanisms of development and progression in glioma, is becoming increasingly urgent [20]. A growing body of research has demonstrated that distinct miRNA expression pro les have been identi ed as a crucial factor responsible for the initiation and progression of various human malignant cancers [21,22]. The differential expression of miRNA targeted genes in tumor tissues may be triggered by several factors, including positional effects of cancer-related genome regions, epigenetic regulatory mechanisms, and alterations in the miRNA processing machinery [23]. MiRNAs are involved in tumor invasion by targeting many critical protein-coding genes [24]. Recent studies have highlighted that miRNAs play important roles in cancer development, indicating that miRNAs have strong potential to become clinically useful diagnostic markers for early cancer detection [25,26]. The multiple miRNA risk score is superior to a single miRNA in terms of statistically robust analysis. Although many studies have reported many prognostic markers implicated in glioma prognosis at present, using miRNA risk scores and sample groupings to validate the model has not been supported by previous studies.
Meanwhile, the network between these miRNA target genes and the DEGs was described, functional analysis and pathway enrichment analysis were performed on the network genes. GO enrichment analysis illustrated that those target genes were primarily participated in the regulation of cell morphogenesis and neuron projection development, modulation of chemical synaptic transmission, synaptic membrane, glutamatergic synapse, neuron to neuron synapse, actin binding, protein serine/threonine kinase activity and channel regulator activity. The enriched pathways of target genes of included cAMP signalling pathway, oocyte meiosis pathway and long-term potentiation pathway and so on. All of these studies showed that these target genes played active role in promoting tumour maintenance to some degree. Thus, the four miRNAs have great signi cance in the regulation of tumour signaling pathways.
Finally, we screened eight hub genes (i.e., SNAP25, MAPK1, CAMK2A, SYT1, NRXN1, CALM1, ATP2B2, and CALM3). The result of Kaplan-Meier indicated that high expression of SNAP25, MAPK1, CAMK2A, SYT1, NRXN1, CALM1, ATP2B2, and CALM3 was associated with favourable prognosis in glioma patients. Therefore, we supposed that these genes were tumor suppressor genes. To further validated the relevance of our hub genes to clinical outcomes, we performed the CGGA and GEPIA database. In this study, the result showed that the glioma histology was signi cant correlation with these hub genes expression. It was well known that the higher the tumor grade, the worse the survival prognosis. Those hub genes were lower expression in higher grade glioma. Many studies had reported that IDH mutant gliomas had better overall survival [30][31][32]. In this study, these hub genes were also higher expression in IDH mutant glioma.
The GEPIA database analysis indicated that those hub genes were signi cantly lower expression in cancer tissues. The survival prognosis analysis result of CGGA database also showed that the higher these hub genes expression, the better prognosis patients had. GSEA and GSVA also revealed that the lower expressions of these hub genes were most closely associated with vascularization and tumor proliferation. These results further veri ed our hypothesis and illustrated the reliability of our research that these hub genes might play more important roles in the glioma development and act as tumor suppressor genes. Synaptosomal-associated protein 25 (SNAP25), a key hub gene in the PPI network that could regulate calcium dynamics and neuronal plasticity, was related to the favorable survival of patients. Many studies have reported that lower expression of SNAP25 contributes to the occurrence and progression of tumours, and the expression of SNAP25 involves in the restoration of dendrites which may lead to terminal differentiation with subsequent loss of tumorigenicity [33]. Therefore, these miRNAs may affect the survival prognosis of glioma patients by regulating these hub genes.

Conclusions
We constructed a four-miRNA prognosis risk score model that had good reliability and strong predictive ability that could represent as a novel prognostic marker in glioma. Further we screened 8 hub genes which acted as the tumor suppressor genes which were regulated by the four miRNAs. The present study analysed the miRNA and target mRNA data of glioma through multiple tumor databases using a bioinformatic approach. To verify the effectiveness and feasibility of this proposed model, other databases and numbers of experiments should be employed in future studies, and the model may provide a promising therapeutic target for glioma treatment.
Declarations SS L and ZZ M wrote the manuscript. XY H, JJ W and DF D analyzed the data and created the images. YC L and JG X contributed to the modi cation of the manuscript. All authors read and approved the nal manuscript.

Funding
This work was nancially supported by the National Key R&D Program of China (2018YFC2001600).
Availability of data and material The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Ethics approval and consent to participate Not applicable

Consent for publication
Not applicable.

Competing interests
The authors have no con icts of interest to declare.  Flowchart of the bioinformatic analysis.

Figure 2
Ten miRNAs associated with overall survival in glioma patients.  (B) Nomography of the four-miRNA risk score model.