Integrated Transcriptome Profiling Identifies Prognostic Hub Genes as Therapeutic Targets of Glioblastoma: Evidenced by Bioinformatics Analysis.

Glioblastoma (GBM) is the most devastating and frequent type of primary brain tumor with high morbidity and mortality. Despite the use of surgical resection followed by radio- and chemotherapy as standard therapy, the progression of GBM remains dismal with a median overall survival of <15 months. GBM embodies a populace of cancer stem cells (GSCs) that is associated with tumor initiation, invasion, therapeutic resistance, and post-treatment reoccurrence. However, understanding the potential mechanisms of stemness and their candidate biomarkers remains limited. Hence in this investigation, we aimed to illuminate potential candidate hub genes and key pathways associated with the pathogenesis of GSC in the development of GBM. The integrated analysis discovered differentially expressed genes (DEGs) between the brain cancer tissues (GBM and GSC) and normal brain tissues. Multiple approaches, including gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were employed to functionally annotate the DEGs and visualize them through the R program. The significant hub genes were identified through the protein-protein interaction network, Venn diagram analysis, and survival analysis. We observed that the upregulated DEGs were prominently involved in the ECM-receptor interaction pathway. The downregulated genes were mainly associated with the axon guidance pathway. Five significant hub genes (CTNNB1, ITGB1, TNC, EGFR, and SHOX2) were screened out through multiple analyses. GO and KEGG analyses of hub genes uncovered that these genes were primarily enriched in disease-associated pathways such as the inhibition of apoptosis and the DNA damage repair mechanism, activation of the cell cycle, EMT (epithelial-mesenchymal transition), hormone AR (androgen receptor), hormone ER (estrogen receptor), PI3K/AKT (phosphatidylinositol 3-kinase and AKT), RTK (receptor tyrosine kinase), and TSC/mTOR (tuberous sclerosis complex and mammalian target of rapamycin). Consequently, the epigenetic regulatory network disclosed that hub genes played a vital role in the progression of GBM. Finally, candidate drugs were predicted that can be used as possible drugs to treat GBM patients. Overall, our investigation offered five hub genes (CTNNB1, ITGB1, TNC, EGFR, and SHOX2) that could be used as precise diagnostic and prognostic candidate biomarkers of GBM and might be used as personalized therapeutic targets to obstruct gliomagenesis.

Introduction

Glioblastoma (GBM) is the most devastating and frequent type of primary brain tumor with high morbidity and mortality. Despite treatment regimens that include surgical resection with radiation and concomitant adjuvant chemotherapy, the median survival time for patients with GBM is 12–15 months, with survival rates of 25% and 10% after 2 and 5 years, respectively.[1−4] Accumulating evidence in recent years shows that GBM consists of the subpopulation of cells displaying various stem cell-like properties including long-term self-renewal with the capacity to generate phenotypically diverse hierarchical neoplastic progeny and stromal cells referred as glioma stem cells (GSCs).[5] GSCs also can recapitulate the essential phenotypes of the original tumor, such as tumor cell heterogeneity, invasiveness, and vascularity promoting resistance to chemotherapy and radiotherapy.[6,7] Thus, neoplastic cells displaying stem-like phenotypes are currently believed to be the main barriers for successful treatment of GBM that associated inexplicably in tumor growth and recurrence after therapy.[8−10]

The biggest challenge in glioma is to monitor the diagnosis and prognosis process. The typical way of disease monitoring in patients is radiographic utilizing computed tomography, magnetic resonance imaging, or positron emission tomography, which is entirely dependent on the experience of the neurosurgeon.[11] Unfortunately, the detection of anatomic changes in the face of progressing or regressing disease is often nonspecific and slow to change by these imaging modalities. As a result, additional surgery is required for definitive tissue diagnosis, or inappropriately waiting for radiographic findings to change as the disease progresses leads to more patient deaths.[12,13] Thus, there is a dire need for more sensitive and specific tumor biomarkers for the early detection of gliomas. On the other hand, literature surveys demonstrate that gene expression profiling affords an unbiased technique to classify tumors and its correlation is more satisfactory than tumor histology, with prognosis.[14,15] Various specific tumor-related biomarkers have been found to better understand the molecular pathogenesis of GBM due to the development of next-generation sequencing technology.[16,17] Some molecular signatures including IDH mutation (favorable prognoses, secondary GBM),[18] copy number variations of chromosomes 9 and 10 and of the telomerase reverse transcriptase promoter gene (TERTp), B-Raf gene (BRAF) and H3F3A mutation, O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation,[19] and 1p/19q co-deletion (chemosensitivity)11 have been clinically used for the detection of gliomas.[19,20] However, the molecular biomarkers were limited in specific classification and precise outcomes in the prediction of GBM.[11,21] Therefore, there is an urgent need to identify novel molecular biomarkers and new therapeutic strategies to improve the diagnostics and prognosis of GBM and elucidate the mechanisms of GBM or increase overall patient survival.

Currently, many chemotherapeutic molecules have been approved by FDA as a standard of care (SOC) for the treatment of GBM.[22] However, only five drugs and one device including lomustine (CCNU), carmustine (BCNU), temozolomide (TMZ), bevacizumab (BVZ), carmustine wafer implants (BCNU wafers), and optune device (TTFields) are mostly used to treat glioma and its symptoms. In comparison to other SOC, only TTFields has shown a marginal increase in the overall survival (20.5 vs 15.6 months) and progression-free survival at 6 months (PFS6) (56% vs. 37%) in GBM patients.[23,24] However, TTFields is the latest addition to the list of FDA-approved drugs, which has not yet been acknowledged worldwide as a part of SOC,[24] emphasizing the significance of the emerging drug molecule.

The combined bioinformatics analysis of transcriptomic data such as microarray and high-throughput sequencing (NGS) data has gradually become a powerful tool for the identification of potential candidate biomarkers and related pathways. In this study, we strive to identify overlapping molecular biomarkers and new therapeutic strategies for GBM and GSC and investigate their potential clinical roles and molecular mechanisms in the pathogenesis of GBM. To address this, we performed an integrated analysis of GBM and GSC samples utilizing one RNA-seq data and one microarray data from the Gene Expression Omnibus (GEO) database and TCGA_GBM from The Cancer Genome Atlas (TCGA) to discover novel overlapped molecular signatures. The identified overlapped differentially expressed genes (DEGs) from all datasets were assimilated with a protein–protein interaction network to ascertain hub genes primarily associated with gliomagenesis. Further, the gene-level expression of hub genes was validated through GEPIA2 and two microarray data sets, and protein-level expression was predicted through “The Human Protein Atlas”. Additionally, the impact of the expression on the overall survival was evaluated to detect a diagnostic or prognostic model that might be helpful for the precise diagnosis and prognostic estimation of GBM. Finally, an miRNA–hub gene network was constructed to explore the epigenetic regulation, and candidate small molecular drugs were examined for GBM.

Results and Discussion

Gene Intersection between Differentially Expressed Genes

The distribution of the altered expression of genes in normal and abnormal samples is often displayed in certain diseases. The significantly deregulated genes in the glioblastoma (GBM) and glioma stem cells (GSC) in comparison with normal stem cells (NSC) were detected using a t-test of the LIMMA package in R and applying cutoff criteria as FC (FDR < 0.05, |log2 FC| > 1). The analysis of dataset GSE119834 between GSC and NSC revealed 596 upregulated and 1152 downregulated genes. A total of 2679 upregulated and 1751 downregulated genes were identified for the GSE119834 dataset by comparing GBM to NSC. The TCGA_GBM dataset analysis disclosed that a total of 6599 genes were deregulated, which included 3482 upregulated and 3117 downregulated genes. The GSE41031 dataset analysis provided 8193 overexpressed and 1110 underexpressed genes. The volcano plot for each analysis is shown in Figure S1. Venn diagram analysis was employed to detect common genes within four analyses, which disclosed a total of 44 intersecting genes including 32 markedly upregulated genes (Figure a) and 12 significantly downregulated genes (Figure b). Table lists the intersecting genes and their symbol.

Figure 1

Venn diagram represents 44 statistically significant intersecting DEGs among four different sample analyses with adjusted P < 0.05 and |log2 FC| > 1 as the cutoff criterion for their significance. (a) Upregulated genes and (b) downregulated genes.

Table 1

List of 44 Intersecting Differentially Expressed Gene (DEGs) with Cutoff Criteria as FC (FDR < 0.05, |log2 FC| > 1)

intersecting genes	gene symbol
upregulated DEGs	ARID5A; METTL7B; PCDH12; HEY1; HOXA11; ZNF107; TNFRSF19; DENND2A; ITGB8; TRIM14; SHOX2; SOX2; STK32A; ITPRIPL1; PDLIM3; C1orf94; TRIM47; GSC; PCDHB9; HERC5; HOXA6; KCNE4; B3GNT5; SLC27A3; HOXC8; OAS3; FREM2; KIAA0040; TNC, HOXD9; NMD3; CXCL16
downregulated DEGs	SORBS2; NEFM; RAB3B; PTPN3; L1CAM; INA; ST6GALNAC5; SLIT2; TMEM200A; SFRP1; NTN4; MAGI1

Functional Annotation for 44 Intersecting DEGs through GO and KEGG Analysis

The biological functions of a set of genes/proteins can be extracted through the enrichment analysis method. To acquire the biological importance of the intersecting genes in the development of GBM and GSC, DAVID was employed. The enriched GO and KEGG pathway terms were considered significant when a term possessed a p-value of less than 0.05. Top enriched GO and KEGG pathway terms are shown in Figure and listed in Table with the details of the gene count and symbols. The biological process (BP) of the GO analysis result indicated that 32 upregulated intersecting genes were markedly enriched in the embryonic skeletal system morphogenesis, cartilage development involved in endochondral bone morphogenesis, innate immune response, regulation of transcription from an RNA polymerase II promoter, embryonic forelimb morphogenesis, etc. (Figure a); 12 downregulated genes were involved in ureteric bud development, cell differentiation, cellular response to heparin, and axon guidance (Figure b). Variation in the molecular function (MF) revealed that 32 upregulaed genes mostly played role in sequence-specific DNA binding (Figure a); downregulated genes were primarily enriched in the structural constituent of the cytoskeleton and laminin-1 binding (Figure b). Changes in the cellular component (CC) of upregulated genes were significantly annotated in chromatin and transcription factor complex; downregulated genes were found in neurofilament, plasma membrane, and cell surface (Figure ). After deep analysis of 32 upregulated intersecting DEGs through GO and KEGG annotation, we found that BP, CC, and MF of three genes, namely KIAA0040, C1orf94, and ITPRIPL1, were unclear. It has been noticed that protein binding is the MF term for C1orf94; however, BP and CC terms were not enriched. Only the CC term, i.e., an integral component of the membrane was enriched for KIAA0040, while the integral component of the membrane and protein binding was annotated for ITPRIPL1 as CC and MF, respectively (Table S2). More interestingly, we found that KIAA0040 and C1orf94 were two uncharacterized proteins whose biological function is still imprecise. Subsequently, it has been observed that many of the genes were understudied such as B3GNT5, METTL7B, PCDHB9, ZNF107, and some others listed in Table S3 with less than 20 publications. Furthermore, the KEGG pathway analysis uncovered that three upregulated genes such as ITGA1 (integrin subunit alpha 8), TNC (tenascin C), and FREM2 (FRAS1 related extracellular matrix 2) were prominently involved in the ECM–receptor interaction pathway (Figures c and 3) and three downregulated genes, i.e., NTN4 (netrin 4), SLIT2 (slit guidance ligand 2), and L1CAM (L1 cell adhesion molecule) significantly enriched the axon guidance pathway (Figures c and 4).

Figure 2

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for the intersecting significant DEGs. (a) Enriched biological function by GO analysis for upregulated genes associated with GBM, (b) enriched biological function by GO analysis for downregulated genes involved in GBM progression, and (c) enriched pathway by KEGG analysis for both up- and downregulated genes related to GBM.

Table 2

GO Enrichment and KEGG Pathway Analysis of Intersecting DEGs Participating in GBM Progression with the Gene Count and Symbol

GO_term	count	p-value	genes
Upregulated Genes
BP:GO:0009952∼anterior/posterior pattern specification	5	1.00 × 10–5	HEY1, HOXA6, HOXD9, HOXC8, HOXA11
MF:GO:1990837∼sequence-specific double-stranded DNA binding	7	2.10 × 10–4	HEY1, GSC, SHOX2, HOXA6, HOXD9, HOXC8, HOXA11
BP:GO:0035115∼embryonic forelimb morphogenesis	3	9.70 × 10–4	SHOX2, HOXD9, HOXA11
BP:GO:0048704∼embryonic skeletal system morphogenesis	3	0.0016	GSC, HOXA6, HOXD9
BP:GO:0000122∼negative regulation of transcription from the RNA polymerase II promoter	7	0.0021	SOX2, HEY1, GSC, SHOX2, ARID5A, HOXD9, HOXC8
BP:GO:0006357∼regulation of transcription from the RNA polymerase II promoter	9	0.0023	HEY1, GSC, SHOX2, ARID5A, ZNF107, HOXA6, HOXD9, HOXC8, HOXA11
MF:GO:0000981∼RNA polymerase II transcription factor activity, sequence-specific DNA binding	8	0.0034	SOX2, HEY1, GSC, SHOX2, HOXA6, HOXD9, HOXC8, HOXA11
CC:GO:0000785∼chromatin	7	0.0035	SOX2, HEY1, GSC, SHOX2, HOXA6, HOXD9, HOXC8
CC:GO:0005667∼transcription factor complex	4	0.0046	SOX2, GSC, ARID5A, HOXA11
BP:GO:0007155∼cell adhesion	5	0.0077	TNC, ITGB8, PCDH12, PCDHB9, FREM2
BP:GO:0060351∼cartilage development involved in endochondral bone morphogenesis	2	0.01	SHOX2, HOXA11
MF:GO:0000978∼RNA polymerase II core promoter proximal region sequence-specific DNA binding	7	0.01	SOX2, HEY1, GSC, ZNF107, HOXA6, HOXD9, HOXA11
BP:GO:0001649∼osteoblast differentiation	3	0.012	SOX2, SHOX2, TNC
BP:GO:0060272∼embryonic skeletal joint morphogenesis	2	0.014	SHOX2, HOXA11
BP:GO:0009954∼proximal/distal pattern formation	2	0.036	HOXD9, HOXA11
BP:GO:0045087∼innate immune response	4	0.055	HERC5, OAS3, ARID5A, TRIM14
Downregulated Genes
MF:GO:0005200∼structural constituent of the cytoskeleton	3	0.0016	NEFM, SORBS2, INA
BP:GO:0071504∼cellular response to heparin	2	0.0025	SFRP1, SLIT2
MF:GO:0043237∼laminin-1 binding	2	0.0035	NTN4, SLIT2
CC:GO:0005883∼neurofilament	2	0.0054	NEFM, INA
BP:GO:0007411∼axon guidance	3	0.0061	NTN4, SLIT2, L1CAM
BP:GO:0001657∼ureteric bud development	2	0.018	SFRP1, SLIT2
CC:GO:0005886∼plasma membrane	7	0.031	MAGI1, RAB3B, SFRP1, NTN4, SORBS2, L1CAM, PTPN3
CC:GO:0009986∼cell surface	3	0.043	SFRP1, SLIT2, L1CAM
BP:GO:0030154∼cell differentiation	3	0.044	SFRP1, SLIT2, INA
KEGG Pathway
UP_hsa04512:ECM–receptor interaction	3	0.007	TNC, ITGB8, FREM2
DOWN_hsa04360:axon guidance	3	0.0048	NTN4, SLIT2, L1CAM

Figure 3

Extracellular matrix (ECM)–receptor interaction pathway generated using the pathview R package. The overexpressed genes in GBM are highlighted in red color.

Figure 4

Schematic representation of the axon guidance pathway produced by the pathview R package. The lower expressed genes in GBM are highlighted in green color.

The extracellular matrix (ECM) serves diverse functions in tissue and organ morphogenesis including the establishment and maintenance of the cell and tissue structure and function.[25] The ECM communicates with cells by specific interactions, mediated by transmembrane molecules, primarily integrins and also possibly proteoglycans or other cell-surface-associated components (Figure ).[26] These interactions can directly or indirectly regulate several cellular processes such as adhesion, migration, differentiation, proliferation, and apoptosis.[27] Subsequently, integrins comprise noncovalently bound alpha- and beta-subunits that acts as mechanoreceptors to establish a force-transmitting physical link between the ECM and the cytoskeleton.[28] Further, the remodulation of the ECM is a crucial mechanism to monitor certain biological processes including morphogenesis, bone remodeling, angiogenesis, and wound repair. However, the dysregulation of ECM dynamics contributes to pathological conditions such as osteoarthritis, inflammatory diseases, tissue fibrosis, and invasive cancer.[29] TNC, ITGB8, and FREM2 were found to be highly expressed and involved in ECM–receptor interactions in our analysis. It has been demonstrated in various research articles that TNC (tenascin C), an extracellular matrix molecule, is overexpressed in many types of cancers such as breast cancer, colorectal cancer, lung cancer, and GBM.[30] TNC acts as a driver for tumor cell survival, proliferation, invasion, stemness, and metastasis in different signaling pathways. In human macrophage cells, TNC stimulates proinflammatory factors such as TNFα via a toll-like receptor 4 and STAT3-dependent mechanism. In GBM, TNC plays a vital role in immunomodulation by interfering in the antitumor function of brain innate immune cells.[31] ITGB8 and FREM2 are overexpressed in some cancers including breast, ovarian, lung, prostate, bladder cancer, and GBM and play a major role in tumor cell metastasis, less adhesion, invasion, and growth.[32] More interestingly, ITGB8 functions as a key driver of transforming growth factor b (TGF-b) to impede the antitumor immunity leading to the tumor growth.[33]

PPI Network Construction and Hub Gene Identification

The importance of protein and biological modules can be simplified through protein–protein interactive relationship network analysis. Hence, the STRING database was adopted to predict the interactive potential relationship between the identified intersecting DEGs with a confidence score of ≥0.04 at the protein level. Subsequently, the interactive relationship network was imported and reconstructed through the Cytoscape plugin. As revealed in Figure a, 54 nodes (genes) and 85 edges (interactions) were established in the protein interactive network. Furthermore, the most relevant hub genes for the progression of GBM and GSC were identified through the cytoHubba plugin of Cytoscape. Owing to the heterogeneity of the biological network, it makes sense to employ many methods for capturing important proteins.[34] Four different algorithms, namely degree of connectivity, DMNC, MNC, and MCC were utilized to extract the hub genes (Figure b; Figure S2). As a result, a group of 12 overlapping genes was discriminated as common DEGs from these four methods (Figure c). Table displays the list of 12 overlapping genes, their descriptions, and scores obtained from different algorithms. Subsequently, functional enrichment analysis of 12 hub genes revealed that most of the genes are involved in various important biological processes, which trigger tumorigenesis such as osteoblast differentiation, positive regulation of transcription from an RNA polymerase II promoter, positive regulation of cell proliferation, cell fate specification, negative regulation of the apoptotic process, etc. Cellular component terms, i.e., cell–cell adherens junction, focal adhesion, cell surface, extracellular space, perinuclear region of cytoplasm, and membrane raft were enriched for 12 overlapped hub genes. Molecular functions of the 12 common hub genes were annotated as estrogen receptor binding, cadherin binding involved in cell–cell adhesion, integrin binding, enzyme binding, and RNA polymerase II activating transcription factor binding. According to KEGG analysis, the highly enriched pathways for 12 overlapping hub genes included signaling pathways governing stem cell pluripotency, cancer pathways, proteoglycans in cancer, Rap1 signaling pathway, and basal cell carcinoma, among others. (Figure d). All of the significant detailed information such as GO and KEGG terms, gene counts, and names are listed in Table S4. Overall, from GO and KEGG enrichment analyses, we collected evidence that our identified hub genes participated in the stemness of tumor cells and various cancer progression activities. These data signposted that these candidate genes have a prime influence on the progression of glioma and can be considered prognostic biomarkers.

Figure 5

Visual representation of the protein–protein interaction (PPI) network of overlapped DEGs and the hub genes. (a) PPI network constructed by overlapping genes. The thickness of the connection line indicates the level of betweenness within proteins. Color gradients from cyan to pink of nodes represent the change of log FC. Green color denotes the coexpression nodes added by the STRING database. (b) Top 15 hub genes screened by cytoHubba plugin of Cytoscape based on the MNC score where a higher score is represented in red color and a lower score in yellow color, (c) Venn diagram representing the overlapping hub genes within these four methods, and (d) GO and KEGG enrichment analyses of 12 hub genes.

Table 3

List of 12 Overlapping DEG Names with Their Four Different (Degree, DMNC, MNC, and MCC) cytoHubba Algorithm Scores

gene symbol	gene description	degree score	DMNC score	MNC score	MCC score
CTNNB1	catenin beta 1	14	0.408739	13	445
EGFR	epidermal growth factor receptor	13	0.424202	11	382
SHH	sonic hedgehog signaling molecule	12	0.492074	11	439
SOX2	SRY-box transcription factor 2	11	0.44117	11	412
BMP4	bone morphogenetic protein 4	11	0.44117	11	320
ITGB1	integrin subunit beta 1	10	0.419005	10	260
SRC	SRC proto-oncogene, nonreceptor tyrosine kinase	10	0.399052	10	148
ISL1	ISL LIM homeobox 1	8	0.408202	8	60
HEY1	Hes-related family BHLH transcription factor with YRPW motif 1	6	0.453784	5	27
SFRP1	secreted frizzled-related protein 1	5	0.648263	5	120
TNC	tenascin C	4	0.473661	4	12
SHOX2	short stature homeobox 2	4	0.378929	4	8

Validation of mRNA Expression of Overlapped Hub Genes in GBM

To explore the prognostic potential of 12 overlapped hub genes (CTNNB1, EGFR, SHH, SOX2, BMP4, ITGB1, SRC, ISL1, HEY1, SFRP1, TNC, SHOX2), we first visualized the expression levels of hub genes in GBM compared with LGG and normal glial cells using GEPIA2 server and two data sets, namely GSE4290 and GSE15824. The validated result of the expression levels of 12 overlapped hub genes from the GEPIA2 server and GSE4290 data set is presented in Figure and Figure S3. Both analyses possessed similar results as most of the hub genes were significantly overexpressed with a p-value of <0.05 in both GBM and LGG except for CTNNB1 and BMP4. As shown in Figure a, from the GEPIA2 server, CTNNB1 was markedly upregulated in both GBM and LGG, while GSE4290 data sets revealed that the expression of CTNNB1 was not significant (Figure b). Similar to CTNNB1, BMP4 exhibited high expression in GBM and LGG through the GEPIA2 server, while the GSE4290 data set displayed downregulation. Then, we revalidated both CTNNB1 and BMP4 gene expression levels with another data set GSE15824 (data not shown). The result uncovered that CTNNB1 was significantly upregulated in GBM with a p-value of 0.00094; however, its expression in LGG was not significant with a p-value of 0.12. The expression level of BMP4 was shown to be downregulated in both microarray data sets, while in TCGA_GBM-and-LGG data sets, it was found to be overexpressed but not significant (Figure S3). Furthermore, we noticed that sufficient data for the expression level of ISL1 were not available in TCGA_GBM-and-LGG data sets. Then, we analyzed ISL1 expression through two data sets (GSE4290 and GSE15824) and found that ISL1 was upregulated in both data sets but the expression level was not significant. From the result, we concluded that HEY1, ITGB1, SOX2, TNC, EGFR, and SHOX2 (only in GBM) were significantly overexpressed in TCGA_GBM-and-LGG data sets, as evidenced by the GEPIA2 server and cross-validation via GSE4290 and GSE15824 data set analysis. Also, we hypothesized that the expression level of CTNNB1 may differ from patient to patient or depend on some external factors for its regulation.

Figure 6

Representative box plots for the validated result of six overlapping hub gene expression levels in GBM and LGG in comparison with N. (a) Validation result of gene expression levels through the GEPIA2 server. TCGA_GBM tissues (n = 163) and TCGA_LGG tissues (n = 518) marked in red color and noncancerous tissues marked in gray (n = 207); a red asterisk indicates a p-value <0.01. (b) Cross-validation of hub gene expression levels in glioma and normal samples based on GSE4290. Two-tailed Student t-tests were used to assess the statistical significance of differences. These results indicate that our findings are reliable with P < .05. GBM: glioblastoma, LGG: low-grade glioma, and N: normal.

After examining the expression level of 12 overlapped hub genes, the protein levels were investigated through immunohistochemical analysis using the HPA database. Remarkably, the protein levels of SOX2 (Figure a), ITGB1 (Figure b), EGFR (Figure c), and TNC (Figure d) were not expressed in normal tissues, while medium and high expression levels of these genes were perceived in glioma tissues (Figure a–d). Additionally, the low-level protein expressions of SRC and CTNNB1 were observed in normal glial tissues, whereas the protein-level expressions of these genes were medium in glioma tissues (Figure e,f). Furthermore, the expressions of SFRP1, BMP4, HEY1, and ISL1 genes were not detected in protein-level expression supporting the expressional validation through data sets except HEY1 (Figure g–i). Moreover, the gene-level expression of HEY1 was noted high in glioma tissues, while its protein-level expression was not detected in glioma tissues. The protein-level expression of HEY1 exhibited a conflicting trend to gene-level expression, which needs more investigation, or it could be postulated that post-translation modification may be playing a major role in HEY1 protein-level expression. The protein-level expression data were unavailable for SHOX2 and SHH genes. Due to insufficient data for ISL1 in TCGA_GBM-and-LGG data sets, ISL1 was not considered for our further studies. Additionally, we were interested in the localization and protein-level expression of all 32 genes and found that most of the genes were found to be medium or highly expressed in glioma cells in comparison to normal glial cells (Table 5). In summary, the transcriptional and translational expression levels of six overlapped hub genes were upregulated in patients with glioma. Further, these 12 genes were carried forward to revalidate the association between their expression level and the overall survival of glioma patients through survival analysis.

Figure 7

Immunohistochemistry images of 10 hub genes in normal brain tissues (cerebral cortex) and glioma cancer specimens derived from the Human Protein Atlas (HPA) database. (a) SOX2, (b) ITGB1, (c) EGFR, (d)TNC, (e) SRC, (f) CTNNB1, (g) SFRP1, (h) BMP4, (i) HEY1, and (j) ISL1.

Survival Analysis of Overlapped Hub Genes

To get a clear vision of the relationship between the overexpression of overlapped hub genes and the overall survival of glioma patients, Kaplan–Meier (K–M) survival analysis was performed on all TCGA glioma cohorts. The result disclosed that most of the genes were significantly correlated with the worst survival of glioma patients with a p-value < 0.05 (Figure ). Additionally, the higher expression of ITGB1 (HR = 2.4; P = 7.7 × 10–6), SHOX2 (HR = 3.2; P = 8.4 × 10–10), CTNNB1 (HR = 1.9; P = 1.8 × 10–6), EGFR (HR = 1.3; P = 0.023), and TNC (HR = 2.2; P = 2.4 × 10–5) mRNAs was markedly associated with worse survival in GBM and LGG patients. This suggested that higher expression levels of these five genes at diagnosis can be considered as unfavorable prognostic genes, which can shorten the overall survival of glioma patients. The lower expression of BMP4, HEY1, SFRP1, SHH, SOX2, and SRC genes was significantly correlated with worse prognosis and lower overall survival in glioma patients. Most importantly, significant high-level expression of SOX2 has been noticed at both gene levels (Figure c) and also in protein levels (Figure a). In contrast, the K–M survival analysis of the TCGA data set suggested that the lower expression of SOX2 was associated with worse survival in glioma patients. Meanwhile, for further confirmation, the CGGA database was also used to evaluate the prognostic role of 11 overlapped genes. Figure S4 reveals that the high-level ITGB1 (p < 0.0001), SHOX2 (p < 0.0001), CTNNB1 (p < 0.01), and TNC (p < 0.0001) mRNAs in glioma patients were significantly interrelated with shorter overall survival, which was analogous to the TCGA data sets except for EGFR with a p-value of 0.38 (Figure S4). The lower expression of BMP4, HEY1, SFRP1, SHH, and SRC genes was markedly correlated with the poorer overall survival of glioma patients except for SOX2, which was in agreement with the result of the TCGA data set. The expression level of SOX2 differed between TCGA and CGGA data sets but not significantly. Overall, these outcomes validated the prognostic value and the relationships between the 11 overlapped hub genes in gliomagenesis.

Figure 8

Kaplan–Meier plotters and log-rank tests for the prognostic value of the hub genes in GBM and LGG generated through the GEPIEA2 server. The dashed lines are upper confidence interval and the lower confidence interval.

The overall observations such as survival analysis, gene-level, and protein-level expression analysis confirmed that overexpression levels of ITGB1, SHOX2, CTNNB1, TNC, and EGFR at diagnosis can be considered as unfavorable prognostic genes, whose expression can reduce the overall survival of glioma patients. The high level of BMP4, HEY1, SFRP1, SHH, and SRC mRNAs at diagnosis can be reflected as favorable prognostic genes, whose presence significantly led to higher overall survival for the patients. For our further studies, we have considered five unfavorable prognostic hub genes, namely ITGB1, SHOX2, CTNNB1, TNC, and EGFR. These five unfavorable genes mostly participated in pathways of inhibition of apoptosis, the DNA damage repair mechanism, activation of the cell cycle, EMT (epithelial–mesenchymal transition), hormone AR (androgen receptor), hormone ER (estrogen receptor), PI3K/AKT (phosphatidylinositol 3-kinase and AKT), RTK (receptor tyrosine kinase), and TSC/mTOR (tuberous sclerosis complex and mammalian target of rapamycin), as shown in Figure S5.

Growing evidence also confirmed our results that integrin subunit beta 1 (ITGB1) is widely expressed in several malignant cancers, including breast cancers,[35] colorectal cancer,[36] and lung cancer.[37] The overexpression of ITGB1 in cancer is mostly associated with tumorigenicity, tumor cell metastasis, proliferation, invasion, and chemotherapeutic resistance, suggesting it to be a prognostic biomarker.[38] Additionally, β1 integrins govern stem cell renewal by stimulating cell division and determining the axis of polarity for asymmetric cell progression.[39]

Previous studies have demonstrated that the elevated expression of the SRY-box transcription factor 2 (SHOX2) gene was significantly associated with cancer pathogenesis driving tumorigenesis, cancer cell development, metastasis, drug resistance, and stemness and recurrence.[40,41] Furthermore, SHOX2 has also been designated as a powerful biomarker for the early diagnosis of various cancers including gliomas,[42] lung cancer,[43] hepatocellular carcinoma,[44] malignant pleural effusion (MPE),[45] bladder cancer,[40] colorectal cancer (CRC),[46] and biliary tract carcinoma (BTC).[47] The expression of SHOX2 induces cancer stem cell-like characteristics, migration, invasion, and recurrence of bladder cancer, and the knockdown of SHOX2 reduced the colony formation and tightened the cellular morphology.[40] The detection of methylated SHOX2 showed more sensitivity and specificity in diagnosing lung cancer at an early stage rather than other commonly used methods.[48]

Catenin β 1 (CTNNB1) gene encodes a beta-catenin protein that primarily participated in cell-to-cell adhesion to maintain tissue architecture and polarity. Beta-catenin also plays a vital role in the Wnt signaling pathway to promote cell growth and division, which determines the specialized functions of cells.[49] However, mutant CTNNB1 has been identified in various cancers including colorectal, liver, thyroid, ovarian, endometrial, and skin cancers, as well as brain tumors.[50,51] Studies suggested that mutation in the CTNNB1 gene interrupts the processing of beta-catenin, leading to the consequent nuclear accumulation of beta-catenin, which traffics into the nucleus and stimulates the unchecked cell division and proliferation, permitting the development and stemness of cancer cells.[52] Additionally, the detection of CTNNB1 circulating tumor DNA at diagnosis can severe as a putative biomarker for assessing hepatoblastoma patient prognosis.[53] The involvement of CTNNB1 in the Wnt signaling pathway facilitates GBM formation, and stemness and leads to drug resistance. Targeting CTNNB1 in the Wnt pathway can inhibit the stemness and/or malignant cellular phenotypes of glioma considered as a therapeutic target for the treatment of GBM.[54,55]

An epidermal growth factor receptor (EGFR) is a transmembrane tyrosine kinase that spans the cell membrane and gets activated by different types of ligands. The activated EGFR stimulates autophosphorylation of tyrosine residues, eliciting the activation of downstream pathways, which leads to cell growth, proliferation, DNA synthesis, and survival.[56] The alterations, mutations, amplification, and overexpression of EGFR have been perceived in cancer tissues such as lung, colon/rectum, pancreas, head and neck, gastrointestinal tract, breast, and brain.[57] The altered and elevated expression of EGFR induces multiple aspects of tumorigenesis including angiogenesis, cell growth, proliferation, invasiveness, and resistance to chemo- and radiotherapy.[58] Furthermore, it has been demonstrated by Lee et al. that EGFR also regulates CTNNB1 expression, localization, and stability and their coactivation might be a malignancy marker of oral cancer.[59]

Tenascin C (TNC) is an extracellular matrix (ECM) glycoprotein that is overexpressed during organogenesis, particularly in the development of the central nervous system, in migration of neural crest cells, and epithelial–mesenchymal transition (EMT), but its distributions typically remain limited in adult tissues.[60] However, increased deposition of TNC has been noticed in different types of malignancies, for instance, melanoma, bladder, breast, uterus (both the cervix and body), lung, tongue, colon, stomach, larynx, skin, and urinary tract cancer.[61,62] The elevated expression of TNC has also been observed in high-grade gliomas that positively correlates with the invasiveness of glioma cells and is suspected to be a potential biomarker of earlier detection of GBM and GSC considered as a better therapeutic target of glioblastoma.[63,64] TNC promotes tumor progression by triggering cell proliferation and differentiation, migration, and angiogenesis and participates in the development of CSC through the HH signaling pathway.[65]

Hub Gene Expression and miRNA Regulatory Network

To further investigate the regulatory relationship between the five hub genes (ITGB1, SHOX2, CTNNB1, TNC, and EGFR) and miRNAs, the miRSystem was employed to predict the targeted miRNAs of hub genes which assimilate seven well-known miRNA target gene prediction programs such as DIANA, miRanda, miRBridge, PicTar, PITA, rna22, and TargetScan.[66] The miRNAs predicted by at least three databases were considered as the targeted miRNAs of hub genes (Table S6). The correlation network associated with the development of GBM was constructed by Cytoscape software. The coexpression network consists of five hub genes (ITGB1, SHOX2, CTNNB1, TNC, and EGFR) and 92 miRNAs (Figure ). SHOX2, CTNNB1, and TNC were associated with most of the miRNAs considered as hub genes. The miRNAs-hub gene coexpression network revealed that has-miR-330-3p and hsa-miR-578 modulate SHOX2 and CTNNB1 mRNAs and hsa-miR-141-3p regulates the expression of CTNNB1 and EGFR. Furthermore, hsa-miR-198 and hsa-miR-495-3p were predicted to be involved in the expression of SHOX2 and TNC.

Figure 9

miRNA–hub gene regulatory network of GBM. The green parallelogram node represented the hub genes. The purple V-shape nodes stand the miRNAs. The edge between two nodes indicates the interaction between hub genes and miRNAs.

A previous report stated that hsa-miR-330-3p promotes tumor invasion proliferation, migration, epithelial–mesenchymal transition (EMT), and survival in multiple cancers, highlighting the role of hsa-miR-330-3p in carcinogenesis.[67] The coexpression of hsa-miR-330-3p with CTNNB1 and SHOX2 may be positively correlated with the tumorigenesis promoting tumor invasion and metastasis. The activity of hsa-miR-578 suppresses the tumor cell growth, migration, metastasis, and glycolysis and boosted cell apoptosis by targeting zinc finger RNA binding protein in breast cancer and by targeting a proliferation-inducing ligand (APRIL) in hepatocellular carcinoma.[68,69] Hsa-miR-578 controls angiogenesis by targeting significant modulators of Focal adhesion, Vascular Endothelial Growth Factor (VEGF), and Hypoxia-Inducible Factor-1 (HIF-1) signaling pathways.[70] The previous studies indicated that the expression of SHOX2 and CTNNB1 might be negatively regulated by hsa-miR-578. Additionally, miR-141-3p inhibited cell growth, proliferation, and migration by targeting EGFR and CTNNB1 in osteosarcoma, renal interstitial fibrosis, and CRC and is considered as a key negative regulator of EGFR and CTNNB1.[71,72] In addition, has-miR-198 suppressed cell mobility and invasion of CRC by repressing TNC.[73] This miRNA–hub gene interaction network provided a basic powerful molecular mechanism to understand the epigenetic regulation of GBM.

Hub Genes and Drug Interaction

Five hub genes were explored for the drug–gene interaction through the DGIdb database. A total of 174 drugs were possibly used to treat GBM by targeting five hub genes, reported and screened through DGIdb database (Figure ). The final list consists of only FDA-approved drugs that were compiled and selected for three promising targets i.e., CTNNB1, ITGB1, and EGFR. The three FDA-approved drugs such as imatinib, temsirolimus, and trametinib were identified as common drugs regulating the function of EGFR and CTNNB1, as displayed in Figure and Table . Imatinib was the first kinase inhibitor approved by the US Food and Drug Administration (FDA) in 2001 to treat chronic myelogenous leukemia (CML) and acute lymphoblastic leukemia with the Philadelphia chromosome.[74,75] Temsirolimus, an inhibitor of mammalian target of rapamycin (mTOR), was used to treat nonsmall-cell lung cancer (NSCLC),[76] refractory solid tumors,[77,78] and renal cell carcinoma.[79] Temsirolimus is mostly used in combination with other inhibitors targeting EGFR to treat cancer types with high expression of EGFR.[80,81] Trametinib is used to treat lung adenocarcinoma patients with an EGFR mutant.[82,83]Table S7 highlights that most of the drugs (62/74) might target EGFR to inhibit its function and some of the drugs, i.e., verteporfin and rindopepimut, were in clinical trials for treating recurrent GBM targeting EGFR.[84,85] Among the listed drugs, osimertinib was deliberated as a potential drug of EGFR to treat cancer patients including GBM by inhibiting the MAPK pathway.[86] Recently, FDA approved several drugs targeting EGFR or the EGFR mutant, for instance, mobocertinib[87] and amivantamab-vmjw,[88] to treat nonsmall-cell lung cancer (NSCLC). Furthermore, a downstream network of EGFR was constructed using the STITCH database to investigate possible effects instigated by EGFR inhibition. The model disclosed that inhibition of EGFR might have possible downstream stimulation on ubiquitin C (UBC) and heat shock protein 90 kDa alpha (cytosolic), class A member 1 (HSP90AA1) (Figure S6). According to the drug–EGFR interaction network, erlotinib plays a vital role in the regulation of EGFR and HSP90AA1. As per our best knowledge, gefitinib,[89] afatinib,[90] erlotinib,[91] Lapatinib,[92] vandetanib,[93] and dacomitinib[94] have been taken for phase-II clinical trials against EGFR to treat GBM, recurrent GBM, and glioma patients. Overall, these data might provide novel insight for targeted therapy in GBM patients.

Figure 10

Candidate drugs targeting the hub gene interaction network. The blue ellipses stand for hub genes and pink diamonds represent drug molecules. All reported drug molecules have been utilized to draw this diagram.

Table 4

Three Common FDA-Approved Candidate Drug Molecules Targeting EGFR and CTNNB1

no.	gene	drug	approved	scores	reference (PubMed ID)
1	CTNNB1, EGFR	trametinib	FDA	1	26343583
					28783719
					26036643
					26582713
					27312529
2	CTNNB1, EGFR	temsirolimus	FDA	1	27016228
					24470557
3	CTNNB1, EGFR	imatinib	FDA	1	26861905
					28762371
					22323597

Several worthy limitations were there in monitoring this investigation. First, we primarily investigated the potential role and functions of hub genes without deeply evaluating the other DEGs and most importantly two uncharacterized genes. In the future, extensive studies are required considering these aspects. Second, TCGA data and two microarray data sets were utilized to validate the expression level of hub genes. As we know that the expression of genes depends on the external and internal environmental factors of a person, hence, for Indian patients, further experimental studies are needed to endorse the above findings. Lastly, the clinical information of GBM patients was not extensively studied owing to the inaccessibility of data. Despite this, our deep investigation provided novel finding to treat GBM patients. Our integrated bioinformatics examination might deliver more precise outcomes in comparison to a single data set analysis. Furthermore, the therapeutic targets and drugs from this exploration are significant and novel and could be used as personalized therapy for GBM patients. Moreover, the miRNA–hub gene interaction network might have revealed the significance of epigenetic regulation of gliomagenesis.

Conclusions

Through integrated analysis and protein–protein interaction network analysis, five potential candidate hub genes (CTNNB1, ITGB1, TNC, EGFR, and SHOX2) were identified, which were positively correlated with the stemness of GBM and negatively correlated with the overall survival of patients. These five unfavorable genes might serve as a potential biomarker for the precise diagnosis, prognosis, and targeted therapy of GBM patients. Additionally, a group of drug molecules were scrutinized which could be exploited for the treatment of GBM patients. Consequently, we believe that clinical validation of our findings with GBM patients would deliver new insights toward designing more potential therapeutic strategies and personalized therapies that may reduce the mortality rates. In our future perspective, we are aiming to design/discover/repurpose drug molecules targeting SHOX2 or TNC and characterize identified uncharacterized proteins to understand their role in gliomagenesis.

Materials and Methodology

Raw Transcriptomic Data Set Collection

Publicly available mRNA expression data sets were retrieved from the TCGA data portal[95] and Gene Expression Omnibus (GEO) database.[96] From TCGA projects, level 3 RNA-seq data set and clinical information of GBM were downloaded to ascertain differentially expressed genes (DEGs). Three microarray data sets such as GSE41031,[97] GSE4290,[98] and GSE15824[99] were obtained from the GEO database for our study, where GSE41031 was utilized for the identification of DEGs, and the gene expression profiles of GSE4290 and GSE15824 were exploited for external validation and confirmation of the hub gene expression level. Additionally, one expression profiling by high-throughput sequencing (NGS) data set, i.e., GSE119834,[100] was acquired from GEO to screen DEGs. The platform used and the sample sizes of selected data sets are summarized in Table S1.

Data Processing and Identification of DEGs

For microarray data sets, all probe ids were converted to the corresponding gene symbol using platform annotation packages. Duplicate genes were removed considering the highest expression values of all probes as representative of the final expression level.[101] For the RNA-seq data set (GSE119834), raw reads were mapped to the reference genome GRCh38. latest using HiSat2.2.0 (hierarchical indexing for spliced alignment of transcripts 2),[102] and SAMtools was used to sort and remove duplicates.[103] The quantification of gene expression was performed using featureCounts.[104] Meanwhile, the linear model microarray analysis (LIMMA) R package was applied for quantile normalization, log2 conversion, and to screen differential expressed genes (DEGs) between GBM, GSC, and noncancerous samples.[105,106] The “Empirical Bayes” model was computed for statistical analysis that generates the statistical significance value (p-value) for each DE gene.[107,108] The genes were filtered with log2 FC > 1.0 and FDR < 0.05 reflected as upregulated or overexpressed in GSC and GBM samples, while those with log2 FC < −1.0 and FDR < 0.05 were considered downregulated or underexpressed.[109−111] A web-based tool “InteractiVenn” was employed to identify intersecting up- and downregulated genes in four types of data analysis.[112]

Functional Enrichment Analysis of DEGs

The enrichment of gene ontology (GO) terms for overlapped DEGs was carried out using the Database for Annotation, Visualization, and Integrated Discovery (DAVIDv6.8).[113−115] The GO analysis mainly includes three individual modules such as biological process (BP), cellular component (CC), and molecular function (MF).[116] The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a comprehensive database that contains information about genomes, biological processes, illnesses, pharmaceuticals, and chemicals.[117,118] The KEGG analysis was executed to identify the biological pathways of DEGs. The Benjamini–Hochberg FDR < 0.05 is set as the threshold for indicating the statistically significant genes.[119,120]

Protein–Protein Interaction (PPI) Network Construction, Module Analysis, and Identification of Hub Genes

Protein–protein interaction (PPI) networks can provide novel insights into the protein function and may assist in classifying key genes[121,122] and pivotal gene modules responsible for the development of glioma stem cells in GBM at an interaction level. In this study, the translated DEGs were employed to construct a protein–protein interaction (PPI) network using the Search Tool for the Retrieval of Interacting Genes (STRINGv11) database.[123,124] The STRING database provides an integrative and critical assessment of interactions between proteins, including prediction and experimental interaction data.[125] The PPIs of translated DEGs with confidence scores greater than or equal to 0.4 were retrieved and imported to Cytoscapev3.9.0 software to reconstruct the network.[126,127] Subsequently, the hub genes were identified by cytoHubba using four methods such as DEGREE, maximum clique centrality (MCC), maximum neighborhood component (MNC), and density of maximum neighborhood component (DMNC).[128−130]

Validation of Hub Gene Expression

The expression levels of hub genes were validated using the Gene Expression Profiling Interactive Analysis 2 (GEPIA2[131]) server and two microarray data sets, i.e., GSE15824 and GSE4290.[132] GEPIA2 is a web-based server that analyzes the expression of tumors and normal samples from the Cancer Genome Atlas (TCGA). GSE4290 consists of 76 LGG samples (includes 26 astrocytomas and 50 oligodendrogliomas), 81 GBM samples, and 23 nontumor samples.[132] GSE15824 contains 30 brain tumor samples (12 primary glioblastomas (GBM), 3 secondary glioblastomas (GBM-2), 8 astrocytomas (Astro) and 7 oligodendrogliomas (Oligo)) and 5 normal brain tissues.[99] Thus, the expression and correlation of the hub genes can be validated in GBM tissues and normal tissues. Further, the basic expression level of hub genes in the cerebral cortex was predicted using the online tool “The Human Protein Atlas”.[133−135]

Overall Survival Analysis of Hub Genes

Survival analysis of hub genes was performed using GEPIA2[131] and CGGA databases.[136−139] The Chinese Glioma Genome Atlas (CGGA) database is a user-friendly, interactive web application for data storage and analysis, which explores over 2000 brain tumor samples from Chinese cohorts. The web resource also can be used to validate the gene expression level and their impact on patients’ survival time.[140] Kaplan–Meier plotter online tool in GEPIA2 was used to determine the overall survival of identified hub genes based on the TCGA database.[141]

miRNA–Hub Gene Coexpression Network

The miRNAs targeting hub genes were predicted by the target prediction database miRSystem,[66,142] which integrates seven well-known miRNA target gene prediction programs: DIANA, miRanda, miRBridge, PicTar, PITA, rna22, and TargetScan.[66,143] The miRNAs identified by at least three predictors were considered as the targeted miRNAs of hub genes. The correlation between miRNA and hub genes was visualized by constructing a coexpression network using Cytoscape.[143,144] In the network, parallelogram nodes denoted the mRNAs, while the V-shape nodes represented the miRNAs.

Drug–Hub Gene Interaction Analysis

Drugs were identified through the Drug Gene Interaction Database (DGIdbv3.0.2;[145]) for hub genes, which can serve as promising targets for GBM.[146,147] DGIdb is a user-friendly drug–gene interaction and the druggable genome data mining database, which mined the data from over 30 trusted sources such as ChEMBL, DrugBank, Ensembl, NCBI Entrez, PharmGKB, PubChem, clinical trial databases, and the literature in NCBI PubMed.[148] Drugs were selected, which are supported by more than one database or have PUBMED references. Finally, the identified drugs, which are only approved by FDA were considered for further study. Further, Search Tool for Interacting Chemicals (STITCHv5.0;[149]) was used to visualize the hub genes and identify the drug interaction network. STITCH integrated disparate data sources for over 430000 chemicals.[150−152]