Screening the Significant Hub Genes by Comparing Tumor Cells, Normoxic and Hypoxic Glioblastoma Stem-like Cell Lines Using Co-Expression Analysis in Glioblastoma

Glioblastoma multiforme (GBM) is categorized by rapid malignant cellular growth in the central nervous system (CNS) tumors. It is one of the most prevailing primary brain tumors, particularly in human male adults. Even though the combination therapy comprises surgery, chemotherapy, and adjuvant therapies, the survival rate is on average 14.6 months. Glioma stem cells (GSCs) have key roles in tumorigenesis, progression, and counteracting chemotherapy and radiotherapy. In our study, firstly, the gene expression dataset GSE45117 was retrieved and differentially expressed genes (DEGs) were spotted. The co-expression network analysis was employed on DEGs to find the significant modules. The most significant module resulting from co-expression analysis was the turquoise module. The turquoise module related to the tumor cells, hypoxia, normoxic treatments of glioblastoma tumor (GBT), and GSCs were screened. Sixty-one common genes in the turquoise module were selected generated through the co-expression analysis and protein–protein interaction (PPI) network. Moreover, the GO and KEGG pathway enrichment results were studied. Twenty common hub genes were screened by the NetworkAnalyst web instrument constructed on the PPI network through the STRING database. After survival analysis via the Kaplan–Meier (KM) plotter from The Cancer Genome Atlas (TCGA) database, we identified the five most significant hub genes strongly related to the progression of GBM. We further observed these five most significant hub genes also up-regulated in another GBM gene expression dataset. The protein–protein interaction (PPI) network of the turquoise module genes was constructed and a KEGG pathway enrichments study of the turquoise module genes was performed. The VEGF signaling pathway was emphasized because of the strong link with GBM. A gene–disease association network was further constructed to demonstrate the information of the progression of GBM and other related brain neoplasms. All hub genes assessed through this study would be potential markers for the prognosis and diagnosis of GBM.


Introduction
One of the most prevailing and highly deadly heterogeneous forms of brain tumors is glioblastoma multiforme (GBM) or grade-IV glioma [1]. The diagnosis of GBM patients is very challenging, and the patient survival rate is 12-15 months even with combinational therapies [2]. The low efficiency of all therapeutic methods including surgery, chemotherapy, and radiotherapy [3] demands pointed to new therapeutic targets for GBM in recent years.
GBM is an extremely heterogeneous tumor at the pathological and cellular level [4,5]. Gene expression and cell proliferation levels also highly differ in GBM. Glioma stem cells The GEOquery package in Bioconductor is used to analyze the GSE45117 dataset [15]. The list of other packages is Biobase, biomaRT, and gplots of R studio [15][16][17]. The Benjamini-Hochberg technique is used to correct multiple testing and calculate the adjusted p-value to avoid Type I errors. A hypergeometric model was performed for both the down and up-regulated DEGs in GO enrichment in categories and KEGG pathway studies [18,19]. Moreover, adjusting the statistical tests locally is done by calculation of a false discovery rate (FDR) [20,21]. A workflow of the data analysis step by step is drawn in Figure 1. tumor (GBT), normoxic glioblastoma stem-like cell lines (GSN), normoxic glioblastoma stem-like cell line, exposed to hypoxia for 48 h (GSN_H), hypoxic glioblastoma stem-like cell line (GSH), and hypoxic glioblastoma stem-like cell line, exposed to normoxia for 48 h (GSH_N).
The GEOquery package in Bioconductor is used to analyze the GSE45117 dataset [15]. The list of other packages is Biobase, biomaRT, and gplots of R studio [15][16][17]. The Benjamini-Hochberg technique is used to correct multiple testing and calculate the adjusted p-value to avoid Type I errors. A hypergeometric model was performed for both the down and up-regulated DEGs in GO enrichment in categories and KEGG pathway studies [18,19]. Moreover, adjusting the statistical tests locally is done by calculation of a false discovery rate (FDR) [20,21]. A workflow of the data analysis step by step is drawn in Figure 1.

Gene Expression Analysis
The study was done in R version 3.6.3. The GBM dataset with low quality and low reads was excluded, yet the rest of the expression set was transformed to a base-2 logarithmic scale. Moreover, gene expression levels were normalized by averaging the treatments before conducting analyses. A general assessment of statistical implementation can be obtained by clustering samples utilizing the correlation metric. Dendrograms based on the correlation metric are useful for identifying outlying samples [22]. Samples presenting atypical distribution of noisy intensities might be an extensive issue. This can be balanced by utilizing non-normalized data to create a box plot of the log intensities, before using absolute signal intensities, which warrants a more even representation of data [23,24].

Statistics and Differentially Expressed Genes
The dataset was retrieved utilizing the GEOquery package in Bioconductor [25]. The statistical significance of p-value < 0.05 and |log2(FC)| > 0.5 was set to determine DEGs

Gene Expression Analysis
The study was done in R version 3.6.3. The GBM dataset with low quality and low reads was excluded, yet the rest of the expression set was transformed to a base-2 logarithmic scale. Moreover, gene expression levels were normalized by averaging the treatments before conducting analyses. A general assessment of statistical implementation can be obtained by clustering samples utilizing the correlation metric. Dendrograms based on the correlation metric are useful for identifying outlying samples [22]. Samples presenting atypical distribution of noisy intensities might be an extensive issue. This can be balanced by utilizing non-normalized data to create a box plot of the log intensities, before using absolute signal intensities, which warrants a more even representation of data [23,24].

Statistics and Differentially Expressed Genes
The dataset was retrieved utilizing the GEOquery package in Bioconductor [25]. The statistical significance of p-value < 0.05 and |log2(FC)| > 0.5 was set to determine DEGs between each treatment category using student t-test for additional review. The study used the heatmap.2 function in the gplots package to generate heatmap plots of DEGs [17,26]. Moreover, the expression values were normalized for each data point in each case of expression data using a log(FC) transformation as follows: where case ij represents the expression value of the gene i of the jth case and mean i is the average expression value of the gene i in the control samples [27].

Weighted Co-Expression Analysis
The union of DEG sets with the corresponding expression values within the treatments was utilized to detect the scale-independent gene modules of co-expression and strongly linked genes generated by WGCNA [28,29]. The topological overlap matrix (TOM) was estimated by adjacency conversion and picked the number (1-TOM) as a measure of distance for detecting genes and modules via hierarchical clustering [30]. Moreover, the parameter blockSize set to 30 and TOMType were assigned to nothing. The soft threshold power β was fixed to 7, the lowermost power founded on the scale-independent topology to construct a weighted gene network.

The PPI Network
A most significant module was selected, and the common hub genes were outlined by treatment weights and module connectivity, which resulted from co-expression analysis. Furthermore, to pick hub seeds, all the candidate genes in the significant module were uploaded by their corresponding average gene expression values to the NetworkAnalyst platform to build the PPI network utilizing STRING Interactome [31]. Additionally, we utilized a zero-order network tool, particularly one that enables to keep hub proteins that interact with each other directly. The proteins were determined with the degree connectivity > 2 (total edges/total nodes) as the hub genes in the PPI network to implement the co-expression network [32].

GO and KEGG Enrichments of the Pathways
Before studying annotations of the DEGs in the most significant module, affy IDs were matched with corresponding gene symbols using the Biomart package [16]. Consequently, gene ontology annotations regarding BP, CC, and MF using DAVID 6.8 and KEGG pathways were studied [19,33]. All of the annotations and subsequent hub genes were cautiously evaluated and separated based on the features of their biological and molecular significances.

Common Hub Gene Survival Analysis
The common hub genes portion in the stemness of GBM was examined via the positively correlated genes in TCGA of the gene expression dataset utilizing the UALCAN database [34]. The gene expression levels of hub genes at significance (p-value < 0.05) are studied. For deep analysis and validation, survival analysis of GBM patients and the significance of survival effect is measured by log-rank test [35]. The GEPIA2 multiple gene comparison tool [36] was utilized to pair TCGA normal and GTEx data of most significant hub genes.

Analysis of a Separate Glioma Gene Expression Dataset
This study retrieved another human glioblastoma and glioma stem cells from the NIH Gene Expression Omnibus (GEO) by typing in the search box the word "glioma" in the GEO database. The GSE124145 gene expression dataset [37] includes total RNAs from three human glioblastoma multiforme tissues (hGBM), six human glioma stem cells (GSCs), and three glioma cell lines from direct tumor resection of a 54-year-old female patient. The DEGs of the GSE124145 gene expression dataset were studied separately. In particular, the expression levels of the most significant hub DEGs from the GSE45117 were focused on to validate the significance of them referencing the GSE124145 dataset. This study further constructed a gene-disease network to confirm the association of the common hub genes with GBM and related genetic diseases or disorders via NetworkAnalysist. The network association information was collected from the DisGeNET database [38], which is a broad platform combining knowledge on human disease-associated genes and variants. Figure 2A illustrates the GSE45117 dataset examining hierarchical clustering that is beneficial for picturing clusters or groups of samples and comparative adjacencies. Figure 2B demonstrates the boxplot of the GSE45117 gene expression dataset within each treatment. The expression values were used to confer a paired comparison of volcano plots and designed heatmaps ( Figures S1 and S2). three human glioblastoma multiforme tissues (hGBM), six human glioma stem cells (GSCs), and three glioma cell lines from direct tumor resection of a 54-year-old female patient. The DEGs of the GSE124145 gene expression dataset were studied separately. In particular, the expression levels of the most significant hub DEGs from the GSE45117 were focused on to validate the significance of them referencing the GSE124145 dataset.

The Gene-Disease Association Network of the Common Hub Genes
This study further constructed a gene-disease network to confirm the association of the common hub genes with GBM and related genetic diseases or disorders via Net-workAnalysist. The network association information was collected from the DisGeNET database [38], which is a broad platform combining knowledge on human disease-associated genes and variants. Figure 2A illustrates the GSE45117 dataset examining hierarchical clustering that is beneficial for picturing clusters or groups of samples and comparative adjacencies.

DEGs of Paired Tumor Subtypes
Following data preprocessing and quality evaluation, this study revealed the expression values from the 12 samples in GSE45117. A set of 2447 DEGs (1133 down-regulated and 1313 up-regulated) in GBT, GSN, GSN_H, GSH, and GSH_N clinical features, under the threshold of p-value < 0.05 and |log2(FC)| > 0.5 (Table 1) were screened for the consequent analyses. The heatmaps and volcano plots of DEGs are provided in Figure S1.

DEGs of Paired Tumor Subtypes
Following data preprocessing and quality evaluation, this study revealed the expression values from the 12 samples in GSE45117. A set of 2447 DEGs (1133 down-regulated and 1313 up-regulated) in GBT, GSN, GSN_H, GSH, and GSH_N clinical features, under the threshold of p-value < 0.05 and |log 2 (FC)| > 0.5 (Table 1) were screened for the consequent analyses. The heatmaps and volcano plots of DEGs are provided in Figure S1.

Co-Expression Analysis
Employing the average linkage-clustering technique ( Figure 3A) based on samples of DEGs of the GSE45117 dataset and to offer a scale-independent network, β = 7 (unscaled R 2 = 0.90) power was chosen ( Figure 3B). The significant five modules were observed ( Figure 3C) using two techniques to assess the relationship among modules and progression of the GBM treatments. Module hierarchical clustering analysis ( Figure 3C) revealed the turquoise module had a greater correlation with treatments than different modules ( Figure 3D). Moreover, the turquoise module was positively correlated with GBT cells but was not correlated with GSH_N and GSH (white shading), and was negatively correlated with GSN_H and GSN (blue shading) ( Figure 4). Hence, the analysis noted the turquoise module of progression in treatments as the most correlated based on the two methods.

Co-Expression Analysis
Employing the average linkage-clustering technique ( Figure 3A) based on samples of DEGs of the GSE45117 dataset and to offer a scale-independent network, β = 7 (unscaled R 2 = 0.90) power was chosen ( Figure 3B). The significant five modules were observed ( Figure 3C) using two techniques to assess the relationship among modules and progression of the GBM treatments. Module hierarchical clustering analysis ( Figure 3C) revealed the turquoise module had a greater correlation with treatments than different modules ( Figure 3D). Moreover, the turquoise module was positively correlated with GBT cells but was not correlated with GSH_N and GSH (white shading), and was negatively correlated with GSN_H and GSN (blue shading) ( Figure 4). Hence, the analysis noted the turquoise module of progression in treatments as the most correlated based on the two methods.

The Hub Module: GO and KEGG Pathway Enrichments
In Table 2 and Figure 5, the BP, CC, and MF of the GO annotations on DAVID are listed for 61 common hub genes in the turquoise module at a significance level (p-value < 0. 05) (See Table S1 for the entire list).

The Hub Module: GO and KEGG Pathway Enrichments
In Table 2 and Figure 5, the BP, CC, and MF of the GO annotations on DAVID are listed for 61 common hub genes in the turquoise module at a significance level (p-value < 0. 05) (See Table S1 for the entire list).

Screening Common Hub Genes and PPI Network
To achieve the analysis of protein-protein interactions in the turquoise module (1025 DEGs), the PPI network was built on the common hub genes list. Through the PPI net- This shows the significant enrichments of most candidate hub genes in BP terms are GO:0006955~immune response, GO:0006952~defense response, GO:0050900~leukocyte migration, GO:0006954~inflammatory response, and GO:0002682~regulation of immune system process. The significant enrichment GO terms in CC are revealed as GO:0005887~integral component of plasma membrane, GO:0031226~intrinsic component of plasma membrane, GO:0005615~extracellular space, GO:0044421~extracellular region part, and GO:0031988 membrane-bounded vesicle. Lastly, the significant enrichment of the hub genes in MF contains GO:0005102~receptor binding, GO:0032403~protein complex binding, GO:0005178 integrin binding, GO:0004872~receptor activity, and GO:0060089~molecular transducer activity. KEGG signaling pathway enrichment study reported that the hub genes were significantly enriched in hsa05150:staphylococcus aureus infection, hsa05144:Malaria, hsa04380:Osteoclast differentiation, hsa04064:NF-kappa B signaling pathway, and hsa05134:Legionellosis.

Screening Common Hub Genes and PPI Network
To achieve the analysis of protein-protein interactions in the turquoise module (1025 DEGs), the PPI network was built on the common hub genes list. Through the PPI network, 123 nodes and 219 edges were selected ( Figure 6) A total of 61 seeds (proteins) which were highly connected with the turquoise module were screened as common hub genes as a result of PPI network construction.
Twenty primary nodes (proteins) with the greatest degrees in the gene expression of GSE45117 dataset common hub genes were picked and are listed in Table 3. These are SRC, SYK, EGFR, LYN, RAC2, SORBS2, IL1R1, PLCG2, S100A8, CCL8, PIK3CG, RHOG, CD44, TLR4, RHOU, EGR1, DAB2, KDR, ITGB2, HCK. Hub genes in Figure 6 following a gradual color change from green to light green and brown to light brown represent expression intensity. The size of nodes represents fold change (FC). The turquoise module was linked to the connectivity degree in the co-expression network. In Figure 6, a negative correlation in the brown nodes and a positive correlation in the green nodes can be seen.
The PPI network of common hub gene KEGG enrichment analysis picks the VEGF signaling pathway as the most significant (p-value < 0.05) pathway ( Table 4). The VEGF signaling pathway demonstrates an important role in many cancers is one of the leading angiogenic regulator pathways involving GBM through hyperactivation. It is also of concern in various biomarkers of tumorigenic progression such as proliferation and survival.

Authentication of the Five Most Significant Hub Genes
The top 20 common hub genes (Table 4) of the co-expression network were confirmed at expression levels and overall survival (OS) in TCGA datasets. After survival analysis using the GSE45117 and TCGA dataset, as illustrated in Figure 7 and Table 5, we identified the five most significant hub genes (IL1R, SORBS2, S100A8, CCL8, DAB2) strongly linked to the progression of GBM. Figure 8A demonstrates a histogram of the five most significant hub gene expression levels by treatments (GBT, GSN, GSN_H, GSH, and GSH_N). The multiple gene comparison analysis of hub genes using TCGA normal and GTEx datasets is shown in Figure 8B.
All five hub genes have shown a decrease in expression with GSN_H and GSH_N treatments. A dramatic drop in the expression of IL1R1 and S100A8 with GSN_H and GSH_N treatments is further reported in Table 5. The relation to gene expression study outcomes shows details about relative expression levels of IL1R1 and S100A8 in normal versus GB tumor samples as demonstrated in Figure 8C,D.

Authentication of the Five Most Significant Hub Genes
The top 20 common hub genes (Table 4) of the co-expression network were confirmed at expression levels and overall survival (OS) in TCGA datasets. After survival analysis using the GSE45117 and TCGA dataset, as illustrated in Figure 7 and Table 5, we identified the five most significant hub genes (IL1R, SORBS2, S100A8, CCL8, DAB2) strongly linked to the progression of GBM. Figure 8A demonstrates a histogram of the Figure 6. The protein-protein interactions in the turquoise module. The gradual color change from green to light green and brown to light brown represents expression intensity. The turquoise module was akin to the connectivity degree in the co-expression network, correlated negatively in the brown and correlated positively in the green nodes. The perimeters of the nodes were analogous to the fold change (FC).   five most significant hub gene expression levels by treatments (GBT, GSN, GSN_H, GSH, and GSH_N). The multiple gene comparison analysis of hub genes using TCGA normal and GTEx datasets is shown in Figure 8B. Orange lines represent the low expression of the most significant hub genes, whereas green lines represent high expression. All five hub genes have shown a decrease in expression with GSN_H and GSH_N treatments. A dramatic drop in the expression of IL1R1 and S100A8 with GSN_H and GSH_N treatments is further reported in Table 5. The relation to gene expression study outcomes shows details about relative expression levels of IL1R1 and S100A8 in normal versus GB tumor samples as demonstrated in Figure 8C,D.

Validation of Five Most Significant Hub Genes
The expression level of IL1R1, SORBS2, S100A8, CCL8, and DAB2 was focused on GSE124145. The outcomes of the analysis to validate the significance are demonstrated in Table 6; it was noted that all of the most significant hub genes (IL1R1, SORBS2, S100A8, CCL8, and DAB2) were expressed significantly upregulated in both GSE124145 as well.

Validation of Five Most Significant Hub Genes
The expression level of IL1R1, SORBS2, S100A8, CCL8, and DAB2 was focused on GSE124145. The outcomes of the analysis to validate the significance are demonstrated in Table 6; it was noted that all of the most significant hub genes (IL1R1, SORBS2, S100A8, CCL8, and DAB2) were expressed significantly upregulated in both GSE124145 as well. We further constructed a gene-disease association network as shown in Figure 9. The most significant hub genes were observed in the intersection of glioblastoma, giant glioblastoma, status epilepticus, and head and neck Neoplasms.  Figure 9. A gene-disease association network of the most common hub genes of co-expression.

Discussion
GBM is the most prevalent, destructive, and fatal brain cancer. Present treatment decisions including surgery, adjuvant therapy, and chemotherapy cannot fully treat the disease, because the tumor is highly defiant to these treatments. GSCs have the self-renewal capacity and are accountable for the tumor resistance in treating GBM.

Discussion
GBM is the most prevalent, destructive, and fatal brain cancer. Present treatment decisions including surgery, adjuvant therapy, and chemotherapy cannot fully treat the disease, because the tumor is highly defiant to these treatments. GSCs have the self-renewal capacity and are accountable for the tumor resistance in treating GBM.
This study identified the DEGs in the GBT-GSN, GBT-GSN_H, GBT-GSH, GBT-GSH_N, GSN-GSN_H, GSN-GSH, GSN-GSH_N, GSN_H-GSH, and GSN_H-GSH_N treatment groups. Prior to the identification of DEGs, normalization is done based on a method used by Chung and Lee (2021). In this study, we performed WGCNA for mutual DEGs derived from the GEO database and reconstructed gene co-expression networks. First, we applied the WGCNA approach to DEGs (Table 1) of the GSE45117 dataset to evaluate the gene expression profile differences including GSCs, and a human GBT sample. As a data analysis approach or a gene filtering (screening) technique to identify groups (modules) of favorably related proteins, the WGCNA is an R software for weighted co-expression study and can be utilized [28,29]. Subsequently, GO and KEGG pathway enrichments were implemented on the turquoise module (p-value < 0.05), which resulted in the most significant module ( Table 2).
The GO pathway analysis is revealed in biological process (BP) terms GO:0006955 immune response, GO:0006952~defense response, GO:0050900~leukocyte migration, GO:0006954~inflammatory response, and GO:0002682~regulation of immune system process. The significant enrichment GO terms in the cellular component (CC) are revealed as GO:0005887~integral component of plasma membrane, GO:0031226~intrinsic component of plasma membrane, GO:0005615~extracellular space, GO:0044421~extracellular region part, and GO:0031988~membrane-bounded vesicle. Lastly, the significant enrichment of the hub genes in molecular function (MF) contains GO:0005102~receptor binding, GO:0032403~protein complex binding, GO:0005178~integrin binding, GO:0004872~receptor activity, and GO:0060089~molecular transducer activity. The KEGG signaling pathway study reported that the hub genes were significantly enriched in hsa05150:staphylococcus aureus infection, hsa05144:Malaria, hsa04380:Osteoclast differentiation, hsa04064:NF-kappa B signaling pathway, and hsa05134:Legionellosis (Table 2 and Figure 5).
Twenty primary nodes ( Figure 6) with the top degrees in the gene expression levels of DEGs are presented in Table 3. These genes can be listed as SRC, SYK, EGFR, LYN, RAC2, SORBS2, IL1R1, PLCG2, S100A8, CCL8, PIK3CG, RHOG, CD44, TLR4, RHOU, EGR1, DAB2, KDR, ITGB2, and HCK, and so-called common hub genes (Table 4 and Figure 7). The VEGF signaling pathway was screened as the most significant (p-value < 0.05) pathway ( Table 4). The VEGF signaling pathway demonstrates an important role in many cancers involving GBM through hyperactivation and is of concern in various biomarkers of tumorigenic progression such as proliferation and survival. Furthermore, the VEGF signaling pathway is one of the leading angiogenic regulator pathways in these tumors [39,40]. Other KEGG pathway enrichments of the top 20 hub genes have resulted in hsa05205:Proteoglycans in cancer, hsa04666:Fc gamma R-mediated phagocytosis, hsa04664:Fc epsilon RI signaling pathway, hsa04662:B cell receptor signaling pathway, and hsa04064:NF-kappa B signaling pathway.
The survival and expression analyses of the common hub genes pick the five most significant hub genes such as Interleukin 1 Receptor Type 1 (IL1R1), Sorbin and SH3 Domain Containing 2 (SORBS2), S100 calcium-binding protein A8 (S100A8), C-C Motif Chemokine Ligand 8 (CCL8), and DAB Adaptor Protein 2 (DAB2). The hub genes are powerfully connected with the development of GBM and they might be useful as potential therapeutic agents as shown in Figure 9. A dramatic drop in the expression of IL1R1 and S100A8 with GSN_H and GSH_N treatments is further reported in Table 5. The association to gene expression analysis results proposes details about comparative expression levels of IL1R1 and S100A8 in normal versus GB tumor as shown in Figure 8C,D. Interleukin-1 signaling is established as an appealing and key therapeutic target for the controlling of glioblastoma-related cerebral edemas [41]. The role of the IL-1 gene family in glioblastoma linked angiogenesis and tumor development has been reported in several studies [41][42][43]. Therapeutically, a knockdown of the IL-1R1 might evaluate inhibition of IL-1 signaling as a novel therapy for GBM [44,45]. In a recent study, SORBS2 in TCGA GBM cohorts has been reported among other genes as possibly being linked with inferior consequences and PDE1C silencing down-regulated their expression [46], consequently proving to be promising concerning patient survival. Furthermore, WFDC1 is an instance of SORBS2bound transcripts which is mediated by SORBS2 and a key metastasis suppressor [47]. Previous research has reported that WFDC1 expression was considerably downregulated in mesenchymal cells in brain cancer [48,49].
Recently, S100A8 has been reported as a prospective indicator with prognostic and diagnostic value in GBM [50]. Gielen et al. (2016) proposed that glioma patients have enlarged quantities of intracellular S100A8/9 compared with healthy controls. Glioma patients further have boosted S100A8/9 serum quantities correlated with amplified arginase bustle in serum. S100A8/9 can be expressed by various myloid cells and tumor cells in glioma, where it can promote tumor cell growth and migration [51].
In a recent laboratory investigation [52], the data uncovered that CCL8 is a tumorassociated macrophage element that resolves penetration and GBM stemness and has resistance to therapies. Moreover, it is reported in another study that CCL8 stimulates the development of tumor cells in the glioma microenvironment [52]. Our study also verified targeting CCL8 offers a new prospect for GBM treatment.
The usual treatment protocol for GBM involves surgical removal of the tumor at a maximal and healthy level, radiation therapy, and temozolomide (TMZ) chemotherapy which is broadly used for silencing GBM. It is shown that the loss of DAB2IP (DOC-2/DAB2 interacting protein) is liable for TMZ resistance in GBM through autophagy-related protein 9b (ATG9B). In a fresh study that listed four subtypes of GBM, DAB2 is reported as one of the fifteen selected genes that belong to the classical (CL) subtype. S100A4 was found in the CL subtype of GBM [53].
To validate our results we further analyzed another clinical GBM gene expression dataset of GSE124145. All these genes were significantly up-regulated in both of the GBM datasets. Thus, targeting these five most significant hub genes (IL1R1, SORBS2, S100A8, CCL8, and DAB2) may offer insightful strategies for GBM treatment. To confirm the association with GBM, we constructed a gene-disease association network as shown in Figure 9. While we validated the results in the TCGA dataset, the accuracy of the results requires molecular and cellular experiments. Thus, by utilizing a sequence of bioinformatics investigation, this current study demonstrated the five most significant hub genes which may be tangled in the diagnosis and prognosis and efficient concerning the characterization of GBM and treatment options. A limitation of this study is the lack of experimental validation to confirm our results. The essential pathways enriched in the candidate hub genes were cell migration, cell motility, localization of cell, locomotion, and leukocyte migration. These results would significantly offer to uncover the progression of GBM.

Conclusions
Here, we focused on one of the GBM gene expression datasets in the GEO database. In this study, IL1R1, SORBS2, S100A8, CCL8, and DAB2 were filtered as the most significant hub genes for the prospective molecular, metabolism, functional studies in GBM. We further validated the five most significant hub genes' significance level using a similar clinical GBM gene expression dataset. Moreover, a gene-disease association network was constructed to confirm the impact of these five hub genes in GBM. These hub genes can be offered to the candidate markers of future research for therapeutic targets in GBM. The rest of the analysis in this study would help to explore the causes of the gliomas, in particular GBM, underlying biological, cellular, and functional events.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13030518/s1, Figure S1: Volcano plots of each paired treatments of DEGs of GSE45117 gene expression dataset; Figure S2: Heatmap plots of each paired treatments of DEGs of GSE45117 gene expression dataset.