Identification of Core Genes and Screening of Potential Targets in Glioblastoma Multiforme by Integrated Bioinformatic Analysis

Glioblastoma multiforme is the most common primary intracranial malignancy, but its etiology and pathogenesis are still unclear. With the deepening of human genome research, the research of glioma subtype screening based on core molecules has become more in-depth. In the present study, we screened out differentially expressed genes (DEGs) through reanalyzing the glioblastoma multiforme (GBM) datasets GSE90598 from the Gene Expression Omnibus (GEO), the GBM dataset TCGA-GBM and the low-grade glioma (LGG) dataset TCGA-LGG from the Cancer Genome Atlas (TCGA). A total of 150 intersecting DEGs were found, of which 48 were upregulated and 102 were downregulated. These DEGs from GSE90598 dataset were enriched using the overrepresentation method, and multiple enriched gene ontology (GO) function terms were significantly correlated with neural cell signal transduction. DEGs between GBM and LGG were analyzed by gene set enrichment analysis (GSEA), and the significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways involved in synapse signaling and oxytocin signaling pathways. Then, a protein-protein interaction (PPI) network was constructed to assess the interaction of proteins encoded by the DEGs. MCODE identified 2 modules from the PPI network. The 11 genes with the highest degrees in module 1 were designated as core molecules, namely, GABRD, KCNC1, KCNA1, SYT1, CACNG3, OPALIN, CD163, HPCAL4, ANK3, KIF5A, and MS4A6A, which were mainly enriched in ionic signaling-related pathways. Survival analysis of the GSE83300 dataset verified the significant relationship between expression levels of the 11 core genes and survival. Finally, the core molecules of GBM and the DrugBank database were assessed by a hypergeometric test to identify 10 drugs included tetrachlorodecaoxide related to cancer and neuropsychiatric diseases. Further studies are required to explore these core genes for their potentiality in diagnosis, prognosis, and targeted therapy and explain the relationship among ionic signaling-related pathways, neuropsychiatric diseases and neurological tumors.


INTRODUCTION
Neuroepithelial tumors are collectively referred to as gliomas, which represent approximately 30% of all primary brain and other central nervous system (CNS) tumors and are the most common intracranial malignancies. Malignant gliomas grow rapidly and spread to brain tissue, leading to distinct surrounding brain edema and obvious focal symptoms. The annual incidence rate is more than three case per 100,000 people. Based on histological criteria, gliomas are generally graded by cell activity and aggressiveness on a scale of I to IV. They are divided into low-grade gliomas (WHO I-II; LGGs) and high-grade gliomas (WHO grade III-IV; HGGs).
LGGs are well-differentiated gliomas. Although these tumors are not biologically benign, their prognosis is relatively good. HGGs, including glioblastoma multiforme (GBM), are poorly differentiated gliomas. Such tumors are malignant tumors with a poor prognosis. GBM is the most commonly occurring type of glioma, with a 5-year survival rate of approximately 5.6% (1). However, this classification system has high interobserver variability, and the survival rate varies greatly between patients with different disease grades (2). Advances in the molecular characterization of gliomas have rebuilt our understanding of their biological process and led to the development of new glioma classification systems (3). The approach of searching efficient biomarkers to improve the characteristics of clinically relevant subgroups may help address this dilemma.
In the past two decades, there have been many studies on establishing specific molecular classifications for gliomas. For example, low-grade gliomas without isocitrate dehydrogenase (IDH) mutations exhibit similar molecular and clinical characteristics to glioblastomas (4). However, patients with IDH mutant gliomas have a better prognosis, and gliomas with low-grade histology tend to develop slowly (5). The methylation status of the O(6)-methylguanine-DNA methyltransferase (MGMT) gene is not only helpful in determining the grade and prognosis of glioma but also an important predictive biomarker of the benefit from alkylating agent therapy in glioblastoma (6,7). With the development of high-throughput technologies, such as microarrays, next-generation sequencing and single-cell sequencing, combined with the use of multidimensional analysis, the understanding of the biological processes of glioma has become increasingly deeper. An increasing number of microarray and sequencing data are being deposited in public databases such as the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and the Cancer Genome Atlas (TCGA, https://cancergenome. nih.gov/), which is convenient for data downloading and reanalysis. These accumulated massive datasets have facilitated the development of classifications of gliomas, especially glioblastomas, to facilitate diagnosis, prognosis prediction, and treatment management.
In this study, we downloaded data for a large number of glioma samples deposited in the GEO and TCGA databases and then analyzed the differentially expressed genes (DEGs) between gliomas (comparing the different grades of gliomas) and normal controls with a bioinformatics approach. Additionally, other methods, including gene ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) pathway analyses, construction of a protein-protein interaction (PPI) network and identification of putative molecules, were used to analyze these data and to determine key genes and pathways. The specific flow of data preparing, processing and analysis is shown in Supplemental Figure 1. The aim of this study was to identify hub genes and core pathways and to provide potential candidate biomarkers for the diagnosis, survival prediction and targeted treatment of gliomas.

Glioma Data Sets
Two sets of GBM microarray data, GSE90598 and GSE83300, were downloaded from GEO database. The data preprocessing method was as follows: the downloaded data were in the format of log2-transformed quantile-normalized signal intensity. First, the probe was mapped to the gene, and the empty probe was removed. The maximum value was taken as the expression level of the gene when multiple probes corresponded to the same gene. We downloaded the TCGA-GBM and TCGA-LGG datasets from TCGA. The data types downloaded included read counts and fragments per kilobase of transcript per million mapped reads upper quartile (FPKM-UQ). The baseline characteristics of patients with LGG (n=511) and patients with GBM (n=154) from TCGA datasets are summarized in Supplemental Table 1.
There is no significant difference between the two groups in gender. Compared with the LGG group, the GBM group had a lower proportion of patients older than 60 years old, and a lower IDH mutation and 1p/19q codeletion rate.

Differentially Expressed Gene Screening
Differential analysis was performed on the GSE90598 dataset using the "limma" (8) package in R. The differential genes were filtered with Benjamini's corrected P value <1e-10-2 and a fold change of >2-fold as the thresholds. The "DESeq2" (9) package in R was used to analyze and filter the differentially expressed genes between the TCGA-GBM and TCGA-LGG datasets by selecting Benjamini's corrected P value <1e-10-2 and fold change more than 2 as thresholds. Subsequently, the intersection set was obtained and used as the final set of differential genes. The results were displayed via the "VennDiagram" (https://CRAN.Rproject.org/package=VennDiagram) package of R. Based on the expression of DEGs in the GSE90598 dataset, cluster analysis was performed. Sample types included healthy human brain tissue samples, GBM patient tissue samples and astrocyte cell line samples.

Gene Ontology and Pathway Enrichment Analysis of Differentially Expressed Genes
To describe the biological characteristics and functional annotations of DEGs, the "overrepresentation" method of the "clusterProfiler" (10) package in R was used to assess GO function and KEGG pathway enrichment, respectively. The top 10 resulting pathways sorted by p value are displayed.
The "GSEA" method was used to identify the enriched KEGG pathways based on the expression profiles between the LGG and GBM subtypes.

Protein-Protein Interaction Network Construction and Module Analysis
A protein interaction analysis of the DEGs was performed to obtain the PPI network in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database via the following website: https://string-db.org/ (11). The genes with the highest degrees were designated as the core genes. The PPI network results contained a total of 153 links. The resulting file was imported into Cytoscape (12) for editing and mapping. Two modules were obtained through the MCODE plug-in in Cytoscape. The genes contained in the main module (module 1) were assessed by KEGG pathway enrichment analysis to obtain the main functional pathways of the module.

Survival Analysis of Putative Hub Genes
Through differential expression analysis and PPI network analysis, 11 core genes of pleomorphic glioblastoma were finally determined. Using the "survival" (https://CRAN.Rproject.org/package=survival) package in R, the patients were divided into high and low expression groups by the quartiles of the expression of each core gene. Survival analysis was performed to assess the association of gene expression levels with survival time. In addition, to combine the expression values of multiple genes for survival prediction, the method described below was used: the expression value of each gene in each sample was standardized (Xij, where i represents the gene and j represents the sample). First, a z-score matrix was created by subtracting the mean between samples and dividing by the standard deviation to compare between genes: Then, in each sample, the average value of the upregulated gene z-score in the core genes was subtracted from the average value of the downregulated gene z-score to obtain the final expression level characteristic index sj: Finally, patients were divided into a high expression group, a low expression group and an intermediate expression group according to the core gene expression characteristic index (sj) of each sample. The R package "survival" was used to analyze patient survival.
Verification was performed with the GSE83300 dataset. Based on the expression levels of the 11 core genes, the patients were divided into a high expression group and a low expression group, and survival analysis was performed to generate survival curves.

Drug Screening for Glioma Treatment
Based on the expression data of the core genes of GBM, combined with the information of all drug molecular targets in the DrugBank (13) database, the significance analysis was performed using hypergeometric tests to find the drugs that could have a therapeutic effect on gliomas. Through differential expression analysis and PPI network analysis, k core genes were selected among the total N genes. Among them, n genes were related to a drug, and M genes had a probability distribution that obeyed the hypergeometric distribution. According to the significance of the P value, the relevance of the selected core genes to the drug was judged. The formula is as follows:

Data Preprocessing and Differential Gene Identification
Two GBM microarray datasets (GSE83300 and GSE90598) were downloaded from the GEO database. The GSE90598 dataset contained 16 GBM and 9 normal samples, and the GSE83300 dataset included 50 GBM patient samples, as shown in Table 1.
The downloaded data were preprocessed to log2-transformed quantile-normalized signal intensity format. The TCGA-GBM and TCGA-LGG datasets were obtained from the TCGA database, as shown in Table 2. The data types downloaded included read counts and FPKM-UQ.
The "limma" package in R was used to perform differential expression analysis on the GSE90598 dataset (Supplemental Table 2), selecting Benjamini's corrected P value <1e-10-2 and fold change of more than 2 as the thresholds. The screening results showed that the 258 DEGs included 124 upregulated genes and 134 downregulated genes. Then, we used the "DESeq2" package in R to analyze the DEGs between the TCGA-GBM and TCGA-LGG datasets (Supplemental Table  3), selecting Benjamini's corrected P value <1e-10-2 and fold change of more than 2 as the threshold, and obtained 7152 DEGs. Among these DEGs, 2797 genes were upregulated and 4355 genes were downregulated. After intersecting the results for each dataset, 150 DEGs were obtained, of which 48 were upregulated and 102 were downregulated. The results were  displayed with the "VennDiagram" package in R using Venn diagrams ( Figure 1). The results of cluster analysis of the DEGs from the GSE90598 dataset are shown in Figure 2. The top 30 DEGs were selected and displayed in the heatmap. Sample types included healthy human brain tissue, GBM tissue and astrocyte cell line samples.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analyses of Differentially Expressed Genes
To describe the biological characteristics and functional annotations of the DEGs, the overrepresentation method of the "clusterProfiler" package in R was used to perform GO and KEGG enrichment analyses. The top 10 pathways were sorted by p value and displayed ( Figure 3). The detailed enrichment results are shown in Tables 3, 4, and Supplemental Table 4. In the biological process category of the GO database, the DEGs were mainly enriched in signal pathways such as modulation of chemical synaptic transmission (GO: 0050804), regulation of transsynaptic signaling (GO: 0099177) and positive regulation of transporter activity (GO: 0032411). In the KEGG pathway enrichment results, the DEGrelated pathways mainly involved GABAergic synapse (hsa04727), morphine addiction (hsa05032) and the oxytocin signaling pathway (hsa04921).  KEGG pathway enrichment analysis of DEGs between the LGG and GBM subtypes was performed using the GSEA method. The enrichment pathway results are shown in Figure  4, Table 5 and Supplemental Table 5. Most of the pathways with different expression levels between LGG and GBM were related to signal transduction.

Module Screening From the Protein-Protein Interaction Network
To obtain the protein interaction network results, further analysis of key gene sets was performed with the protein interaction information in the STRING database. From this analysis, 11 high degree (>= 7) genes were selected as the core gene set, including GABRD, KCNC1, KCNA1, SYT1, CACNG3, OPALIN, CD163, HPCAL4, ANK3, KIF5A, and MS4A6A. The PPI network revealed a total of 153 links. The resulting file was input into Cytoscape for further analysis and mapping ( Figure  5). Through the MCODE plug-in in Cytoscape, two main modules were obtained. Module 2 contained only four genes, and the links were weak. Therefore, GO term functional enrichment and KEGG pathway enrichment analysis were performed on the gene set in module 1, and the main functional pathways of the module are shown in Tables 6, 7, and Supplemental Table 6. In order to further analyze whether these core genes still have significant differences in different histological subtypes of glioma. We performed differential analysis of the transcriptome expression data between LGGs and GBMs from 320 cases of astrocytoma in TCGA datasets, and found that among the 11 core genes selected above, 10 genes: GABRD, KCNA1, SYT1, CACNG3, OPALIN, CD163, HPCAL4, ANK3, KIF5A, and MS4A6A were differentially expressed (Supplemental Figure 2A). We also evaluated whether these  core genes have significant differences in different molecular subtypes through performing differential analysis of the transcriptome expression data between LGGs and GBMs without IDH mutation or without 1p/19q deletion, and we found that between IDH wild-type TCGA-LGG and TCGA-GBM samples, 8 genes (GABRD, KCNA1, SYT1, CACNG3, OPALIN, HPCAL4, and KIF5A) showed differential expression patterns (Supplemental Figure 2B). In addition, there were

FIGURE 4 | KEGG pathway analysis of DEGs between LGG and GBM subtypes. (A)
The KEGG pathway enrichment analysis of differential genes between LGG and GBM subtypes; (B) GSEA analysis of DEGs significantly expressed between LGG and GBM subtypes using the KEGG pathway datasets.

Survival Analysis of Putative Hub Genes and Verification
Through differential gene analysis and MCODE module analysis, 11 core GBM genes were ultimately determined. According to the expression levels of core genes, population quartiles were used as thresholds, and patients with GBM were divided into core gene high expression and low expression groups. Survival analysis was performed, and the survival curves are shown in Figure 6A. Among them, two genes with a significant correlation between gene expression and survival (p value < 0.05) included GABRD and SYT1. However, the relationship between other genes and survival was not very significant. For further analysis, the expression levels of the 11 core genes were combined, and patients were divided into high expression, low expression, and intermediate expression groups according to the overall expression level of the core genes, and survival analysis was performed. The results are shown in Figure 6B. When comparing core gene expression and survival, the p value was 0.11, and the correlation was not significant. Subsequently, the survival analysis results were verified in the GSE83300 dataset. By combining the expression levels of core genes, the feature index sj was calculated, and 50 samples were divided into high-, low-and medium-expression groups by quartiles. The p value when comparing core gene expression and survival was 0.039. The expression levels of the 11 screened core genes in the validation group were significantly correlated with patient survival. We analyzed the protein expression level of the 11 core genes in clinical specimens from The Human Protein Atlas. Protein expression level of survival-related gene SYT1 was shown in Supplemental Figure 3. It had positive strong expression in GBM cerebral cortex tissues and negative weak expression in LGG cerebral cortex tissues.

Core Gene Related Drug Targets Screening
The drug related data were downloaded from the DrugBank database and significance analysis was performed by hypergeometric test. The significantly (p value < 0.01) related drugs are shown in Table 8 and Supplemental Table 7 and include tetrachlorodecaoxide, ketamine, cinolazepam, quazepam, fludiazepam, clotiazepam, adinazolam, prazepam, estazolam, and halazepam. Among them, tetrachlorodecaoxide is related to the treatment of cancer, and the other drugs are related to the treatment of brain-related diseases such as anxiety and insomnia.

DISCUSSION
A large number of current studies indicate that the occurrence and development mechanism of glioma may be related to the interaction of polygenic genetic factors with environmental factors (14-16). As early as 2002, some researchers developed a classification system based on DNA microarray gene expression data. According to gene expression characteristics, childhood medulloblastoma can be distinguished from other brain tumors, including primitive neuroectodermal tumors (PNETs), atypical teratoid/rhabdoid tumors (AT/RTs) and malignant gliomas (17). Subsequent studies have found that although all belong to the same grade of astrocytoma, patients with different molecular expression patterns have different clinical progression features (18). Furthermore, researchers found that molecular subclasses are closely related to the prognosis, disease progression and treatment of glioma (19)(20)(21). With the development of nextgeneration sequencing technology, neurooncologists have a better understanding of gliomas, especially glioblastomas. IDH1 mutations were found in 12% of glioblastoma patients by second-generation sequencing. Such mutations often occur in young patients and patients with secondary glioblastomas and are associated with favorable overall survival (22). Based on the work of the TCGA network, multidimensional genomic data have been integrated to classify GBM into proneural, neural, classical, and mesenchymal subtypes based on the aberrations and gene expression of EGFR, NF1, and PDGFRA/IDH1. Further research revealed that classical subtype patients benefit the most from aggressive therapy, while proneural subtype patients do no benefit (23). Since then, the clinical importance of tumor protein p53 (TP53) abnormalities, deletions involving chromosomes 1p and 19q, MGMT promoter methylation status, abnormalities in the PTEN tumor suppressor gene and the BRAF oncogene, and IDH mutations in gliomas have become better illuminated (5,24,25). This refined classification of glioma into different molecular entities may lead to the use of different treatment regimens (26,27).
In the present study, we obtained the gene expression profile from the GEO and TCGA databases and then identified the DEGs in glioma using bioinformatics analysis. Functional enrichment analysis showed that these DEGs were mainly involved in signal transduction and similar functional patterns in glioma. Most of these DEGs were mainly related to signal transmission pathways such as modulation of chemical synaptic transmission, regulation of transsynaptic signaling and positive regulation of transporter activity. In addition, the main enriched pathways involving the DEGs were the GABAergic synapse, morphine addiction and oxytocin signaling pathways. Chemical synaptic transmission is one of two main modalities of synaptic transmission. The released neurotransmitter transfers information from one cell to an adjacent cell, which is required for communication between neurons (28). A recent study showed that synaptic transmission plays a role in the pathogenesis of neurological diseases such as epilepsy, and its chemical modulators may become potential therapeutic targets (29). Transsynaptic signals between the preand postsynapse are essential for synapse assembly. The dysregulation of transsynaptic signaling could cause severe synapse loss and impair many facets of organization transsynaptically and cell autonomously (30). At the same time, the enriched KEGG pathways, such as the GABAergic synapse and oxytocin signaling pathways, are also mainly related to signal transduction and are involved in neuroanatomy and pathophysiological processes of nervous system diseases (31)(32)(33). Reports on signal transduction have mainly focused on neurodevelopment and nervous system diseases such as epilepsy and neurodegenerative diseases, and few reports have mentioned a relation to neurological tumors. However, a recent study showed that antitumor therapy could contribute to seizure control and that antiepileptic drugs might have beneficial effects on tumors, which provides new ideas for future tumor research (34).
The STRING database is an open database used to assess proteins and their functional interactions for a full understanding of biological phenomena (11). In the present study, this PPI database was utilized to construct interactive networks by screening DEGs. Through online analysis, the core gene set containing 11 hub genes, GABRD, KCNC1, KCNA1, SYT1, CACNG3, OPALIN, CD163, HPCAL4, ANK3, KIF5A, and MS4A6A, was identified. Among the core genes, the expression levels of GABRD and SYT1 were closely related to the survival time of glioma patients. GABRD (gamma-aminobutyric acid type    A receptor subunit delta) is a member of the g-aminobutyric acid A (GABAA) receptors, a class of transmembrane ligand-gated chloride channels expressed in the mammalian brain. dGABAA receptors have been reported to participate in certain DGdependent memory behaviors and facilitate neurogenesis. Gabrd (-/-) mice exhibited impaired recognition memory and contextual discrimination because the migration, maturation, and dendritic development of adult-born neurons were impaired (35). In adult IDH wild-type (WT) diffuse LGG, GABRD expression was independently associated with overall survival (OS) status and showed a moderate negative correlation with tumor-infiltrating macrophages (TIMs) (36). Our analysis showed that among GBM patients, the prognosis of the high GABRD expression group was worse than that of the low expression group. SYT1, synaptotagmin 1, is a calcium-binding protein that triggers fast Ca(2+)-dependent transmitter release in response to membrane depolarization (37). SYT1 not only participates in various biological processes, such as lipid transport (38) and the pathophysiology of neurological diseases (39) but also may play a crucial role in the pathogenesis and treatments of various tumors (40)(41)(42). Nord H et al. analyzed the detailed genetic profiling of a set of 78 glioblastomas and found that eight genes, such as SYT1, might be considered novel candidate oncogenes (43). A later study further analyzed more glioblastoma and normal brain tissue samples from the TCGA and GEO databases, and the results showed that SYT1 is a core gene among 552 differentially expressed genes (44). Our findings are similar to those of the two previous studies. KCNC1 encodes KV3.1, a subunit of the KV3 voltage-gated potassium ion channels. A recurrent de novo mutation, c.959G>A (p.Arg320His), in KCNC1 was identified as a new major cause for progressive myoclonus epilepsies (PMEs) (45). The voltagegated K+ channel gene KCNA1 (KV1.1) is also closely related to epilepsy, as confirmed by animal experiments (46). CACNG3, a member of calcium channels, was evaluated as a susceptibility locus by linkage and association analysis in childhood absence epilepsy (CAE) (47). TMEM10/OPALIN, which encodes a novel transmembrane glycoprotein, could be used as a specific marker for myelinating oligodendrocytes and perhaps for the evaluation of myelination diseases (48). Immunohistochemistry experiments revealed that the macrophage scavenger receptor CD163 could contribute to neuropathological changes in "high inflammation" schizophrenia brains (49). The higher expression level of hippocalcin like 4 (HPCAL4), a neural calcium sensor, was detected in the lateral nucleus of the amygdala (LA) with stimuli associated with danger through fear conditioning (50). Ankyrin-3 (ANK3), encoding the adaptor protein Ankyrin-G (AnkG), has been indicated as an important factor in the pathophysiology of schizophrenia and bipolar disorder (51). Kinesin family member 5A (KIF5A) has been reported to regulate neuronal surface expression of GABA(A)Rs via an interaction with GABA(A)Rassociated protein (GABARAP), indicating that KIF5A could be involved in inhibitory neural transmission regulation related to epilepsy (52). Another hub gene, membrane spanning 4-domain A6A (MS4A6A), has been identified as a new Alzheimer's disease (AD) susceptibility locus by genome-wide association studies (53). With almost no exceptions, all the above-mentioned hub genes are involved in the process of signal transduction in the nervous system and are related to neurodevelopmental, neurodegenerative and neuropsychiatric disorders. In addition, GABRD and SYT1 are implicated in gliomagenesis, and all other core genes, with the exception of OPALIN, have been reported to play a role in tumorrelated diseases, including glioma. As our investigation has shown, abnormal signal transduction is the important property of glioma, and further in-depth understanding of the molecular mechanisms holds promise for the development of cancer drugs that target this process. Interestingly, the core gene-related drugs were identified by combined significance analysis of the expression of core genes and the data from the DrugBank database. The top 10 candidate drugs filtered by p value included tetrachlorodecaoxide, ketamine, cinolazepam, quazepam, fludiazepam, clotiazepam, adinazolam, prazepam, estazolam, and halazepam. Among them, tetrachlorodecaoxide (TCDO, WF10), a chlorite-based nontoxic compound, is implicated in oncotherapy, and the other drugs are involved in the therapeutic management of brain-related diseases such as anxiety and insomnia. As early as 30 years ago, TCDO was shown to improve the oxygenation status of solid tumors and then sensitize tumors to treatments, such as radiotherapy (54). Animal experiment results showed that TCDO could act as a therapeutic agent in acute radiation syndrome and that it significantly reduced the carcinogenic risk in rats after exposure to ionizing radiation (55). Subsequent researchers designed randomized controlled trials in an attempt to evaluate the effectiveness of TCDO in treating oral mucositis or late hemorrhagic radiation cystitis in patients receiving radiochemotherapy, and the results showed that TCDO significantly reduced the recurrence of objective hematuria but was not found to be effective in treating mucositis (56,57). Although the other identified drugs have not been reported to play therapeutic roles in human tumors, particularly gliomas, our analysis provided a starting point for understanding the commonalities between nervous system tumors and neurodevelopmental or neurodegenerative diseases. In addition, studies have indicated that antitumor therapy can contribute to seizure control and that antiepileptic drugs might have beneficial effects on tumors (34), which also provides new ideas for the indepth study of the common molecular mechanisms and shareable therapeutic regimens of both neuropsychiatric diseases and neurological tumors. In summary, our investigation, which is based on a large number of cases and bioinformatics methods, can help researchers reexamine the data from a new perspective. Overall, the prognostic significance and subtype-based expression of hub genes were evaluated in a large cohort and verified experimentally, indicating that these hub genes from the bioinformatics analysis assume on oncogenic role in the progression of glioma.

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
JY contributed to the publication search, data extraction, draft writing, and conception and design. QY contributed to the quality assessment, conception and design, statistical analysis, and editing. All authors contributed to the article and approved the submitted version.

FUNDING
The study was supported by grant from the National Natural Science Foundation of China (81700600 to QY).

ACKNOWLEDGMENTS
The results in this research were based upon data from the Gene Expression Omnibus and The Cancer Genome Atlas established by the NCI and NHGRI. Information about GEO and TCGA and the investigators and institutions that constitute the GEO and TCGA research network can be found at http://www.ncbi.nlm. nih.gov/geo/ and http://cancergenome.nih.gov/.
Supplementary Figure 3 | The expression of SYT1 in cerebral cortex tissue of GBM patient and negative weak expression in cerebral cortex tissue of LGG patient. Images were taken from the Human Protein Atlas (http://www.proteinatlas.org) online database. Four samples of LGG and six samples of GBM were analyzed. The ID numbers of patients with LGG were 122, 164, 2858, and 2890. The ID numbers of patients with GBM were 3, 1587, 1608, 1644, 2811, and 2849. The antibody used in the staining was from Sigma-Aldrich (HPA008394, 1: 90). Antigen retrieval did by HIER pH6.