Bioinformatics Analysis of Key Genes and Pathways for Medulloblastoma as a Therapeutic Target

Introduction: One of the major challenges in cancer treatment is the lack of specific and accurate treatment in cancer. Data analysis can help to understand the underlying molecular mechanism that leads to better treatment. Increasing availability and reliability of DNA microarray data leads to increase the use of these data in a variety of cancers. This study aimed at applying and evaluating microarray data analyzing, identification of important pathways and gene network for medulloblastoma patients to improve treatment approaches especially target therapy. Methods: In the current study, Microarray gene expression data (GSE50161) were extracted from Geo datasets and then analyzed by the affylmGUI package to predict and investigate upregulated and downregulated genes in medulloblastoma. Then, the important pathways were determined by using software and gene enrichment analyses. Pathways visualization and network analyses were performed by Cytoscape. Results: A total number of 249 differentially expressed genes (DEGs) were identified in medulloblastoma compared to normal samples. Cell cycle, p53, and FoxO signaling pathways were indicated in medulloblastoma, and CDK1, CCNB1, CDK2, and WEE1 were identified as some of the important genes in the medulloblastoma. Conclusion: Identification of critical and specific pathway in any disease, in our case medulloblastoma, can lead us to better clinical management and accurate treatment and target therapy.


Introduction
Medulloblastoma (MB) tumors classify as a embryonal tumors of central nervous system which occurs most frequently in children with the rate of 15-20% (Bloom et al., 1969;Northcott et al., 2012). They are commonly found in children between the ages of three and eight, with a higher occurrence in males (Rutkowski et al., 2010). Medulloblastoma patients indicate a low rate of survival and higher mortality (Armstrong et al., 2009;Edelstein et al., 2011). The exact etiology of medulloblastoma has not yet been determined in most cases, however, it seems that its pathogenesis factors include the following: Genetic factors, Head injury, N-nitroso compounds, electromagnetic field exposure and Race (Gilbertson, 2004). Remarkable genomic heterogeneity among medulloblastomas patients have been demonstrated by recent genome-wide profiling studies that provide evidence for the existence of molecular subclasses within medulloblastoma that shows different respond to treatment, in fact paying attention to these subclasses may lead to a better outcome in medulloblastoma treatment Taylor et al., 2012). Recently, methods for treatment of patients with medulloblastoma are consist of the common methods using for any other cancers treatment, such as radiotherapy, chemotherapy, and surgery (Van Dyk et al., 1977;Cumberlin et al., 1979). But all of these methods do not provide targeted and accurate treatment for medulloblastoma also, these methods are followed by many other Side effects for patients. For these reasons, target therapy is considered as a today's challenge and issue in aim to finding better treatment outcome for cancer. In this type of treatment, molecular targets are commonly used, such as genes, signaling pathways, and proteins (Kesari et al., 2005;Movafagh et al., 2007;He et al., 2017). Thus, identification and do a right action for treatment of patients in medulloblastoma can be effective and help to better treatment and a high survival rate (Rossi et al., 2008;Huse and Holland, 2010).
A medical researcher is interested in Bioinformatics methods and statistical analysis to predict and determine the underlying pathways and mechanisms and achieve distinct treatment for many diseases (Yeh et al., 2006;Huang et al., 2008a;Elbers et al., 2009). These biological mechanisms help physicians and drug designers to overcome health care problems and decrease the mortality level. In this regard, data analysis and its applications can be more effective, if used properly. Therefore the application of data analysis on medulloblastoma data sets, do the appropriate test and extract useful output can develop therapeutic methods.
One of the recent techniques in which many studies using is microarray data analysis, Gene set enrichment analysis (GSEA) and pathway analysis to recognize set of the most important genes, their ontology, and pathways in any disease pathogenesis, which can be used as a target for target therapy (Hong et al., 2009;Peng et al., 2010). Using these new technologies in medulloblastoma may play a pivotal role in the treatment processes and an effective therapeutic strategy.
Several Bioinformatics tools are proposed for microarray data analysis, GSEA and pathway analysis. Gene set enrichment analysis (GSEA) (also functional enrichment analysis) is a method to identify classes of genes or proteins that are over-represented in a large set of genes or proteins and may have an association with disease phenotypes. There are many tools for Functional GSEA study, including Broad Institute and FunRich (Pathan et al., 2015;Das and Bhattacharya, 2017).
Pathway analysis tools offer intuitive bioinformatics tools for the visualization, interpretation, and analysis of pathway in order to support basic research, genome analysis, modeling, systems biology, and education. There are several databases for a pathway analyzing including, PANTHER, Reactome, and KEGG (Zhong et al., 2010;Dabaghian et al., 2015).
Cytoscape is a powerful and popular software platform for both visualizing complex networks and integrating these with any type of attribute data. Its additional features are available as plugins. Plugins perform network and molecular profiling analyses, new layouts, additional file format support and connection with databases and searching in large networks. It customizes network data showing by using powerful visual styles. Cytoscape can analyze large network of approximately more than 100 thousand nodes by its rendering engine (Shannon et al., 2003;Cline et al., 2007;Kohl et al., 2011).
Regarding the importance and critical role of medulloblastoma and its treatment as well as the understanding the effective role of pathways, it seems essential to investigate the specific pathway and use the appropriate bioinformatics tool. This current study aimed at using bioinformatics tools to predict the important pathways in medulloblastoma that can hold a value for target based therapeutic means.

Microarray gene expression data analyses
Microarray gene expression data (GSE50161) of 22 medulloblastoma patients and 13 healthy individuals were extracted from NCBI Gene Expression Omnibus (GEO) database, which based on GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix Inc, Denver, USA) and deposited by Donson et al., (2014) and Griesinger et al., (2013). This data was submitted in September 2013 and updated in January 2018. Gene expression profiles were generated from surgical tumor and normal brain samples (n=35) using Affymetrix HG-U133plus2 chips (Platform GPL570). This data was analyzed by the affylmGUI package in R program and created unnormalized data, then these data were normalized by Robust Multiarray Average (RMA) algorithm using the Affy package in Bioconductor (Wettenhall et al., 2006). |logFC| ≥1 and a P value <.05 were considered as threshold values. In this step, the normalized box plot was drowned to investigate samples differences and MA plot was used to approve sample normalization. Next, all normalized and appropriate data were analyzed by the linear method, that all result was demonstrated through the volcano plot and differential expression (logfc). Through this analyses, differentially expressed genes (DEGs) were determined and then genes were sorted for their p-value (a P-value of <0.05 was considered statistically significant) and log fold change (logfc greater than 4 and less than -4 was considered) to determine the specific genes in medulloblastoma tumors.

Construction of Protein-Protein Interaction (PPI) network
STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database compose of experimental data, computational prediction method, and using several classification systems including GO, KEGG, and Pfam to show protein-protein interaction network. The STRING visualizes protein networks and shows the biological relationships of these gene products. In the present research, protein interaction analysis was performed by using Cytoscape (STRING tool kit) to understand DEGs network (p<0.05 and coefficient=0.7 were significant) (Shannon et al., 2003;Szklarczyk et al., 2016).

Gene enrichment analysis
FunRich is a stand-alone software tool that is used for functional enrichment of the gene set. Also, FunRich (Version 3.1.3) can leverage the Gene Ontology (GO) to determine the biological processes, protein domains and cellular components represented in the gene profile. In this study, DEGs that their protein product was associated with each other, were selected and analyzed to investigate the function of DEGs. A p-value of < 0.05 was considered significant (Pathan et al., 2015).

Pathway analysis
The Database for Annotation, Visualization and Integrated Discovery (DAVID) is a bioinformatics tool to investigate the biological meaning behind a list of genes. BiNGO, a Cytoscape plugin, provides specific data to facilitate the interception of biological meaning of experimental data. For this step, we used data from STRING analysis, so that we select the genes that their medulloblastoma. The p-value<0.05 was as a threshold value.

Protein network and gene enrichment analyses
Interaction analyses were shown 238 nodes, 51 edges, 0.429 average node degree and 0.106 Avg. Local clustering coefficient for DEGs in medulloblastoma (PPI enrichment p-value= 1.0e-16) (Figure 3). Approximately 53 proteins were found that they had interacted together.
FunRich arranged the DEGs that interact with each other according to STRING analysis, into two functional categories including, biological processes and cellular component. The parameters for interception of biological processes were indicated in Figure 4. The blue chart for a percentage of genes, the yellow chart for p-value=0.05 reference and the red chart for showing the p-value were determined. Overall, biological processes were included in the mitotic cell cycle, DNA replication, cell cycle checkpoint and mitotic M-M/G1 phases, in which 73.1% and 57.7% belong to the mitotic cell cycle (p<0.001) and DNA replication (p<0.001) (Figure 4). The gene expressions involved in mitotic cell cycle and DNA replication processes might play a significant role in the initiation and progression of medulloblastoma.
The 82.1% percent of the cellular component of DEGs was the nucleus (p-value=0.005) ( Figure 5). 51.9 percentages of the genes had a coiled-coil domain which is found in oncoprotein and transcriptional factors (p-value=0.217). A significant domain in DEGs protein products are associated with each other. Finally, their related pathways were indicated and visualized by using the DAVID database and BiNGO tool kit. In all analyses, a P-value of <0.05 was considered statistically significant (Movafagh et al., 2011;Huang et al., 2008b;Huang et al., 2008a).

Predictions of miRNA-target genes and their enrichment analysis
The regulatory elements miRNA of DEGs was extracted from validated miRNAs-target databases including miRTarBase and TarBase. Furthermore, the molecular functions of the target genes of miRNA were identified by using the FunRich software (Pathan et al., 2015;Khazaei Koohpar et al., 2015).

Differentially Expressed Gene (DEGs) in medulloblastoma
All data were normalized and had suitable quality and value according to the q-q plot and box plot (Figure 1). Volcano plot was also designed to show expression changes ( Figure 2). After data correction and normalization, a total of 21,051 genes was identified in 22 medulloblastoma samples.
Based on the analysis, a total of 249 DEG compose of 145 genes had been expressed highly and about 104 genes had been shown to decrease expression in  Mad3-BUB1 (p-value=0.003), which is crucial for preventing cell and play an important role in mitotic checkpoint.

Pathway enrichment analyses
The significantly enriched pathways, by DEGs, included the cell cycle pathway, p53 signaling pathway, FoxO signaling pathway. Furthermore, cell cycle checkpoint and G2/M Checkpoints were enriched in the cell cycle pathway (Table 1) ( Figure 6).

miRNAs of DEGs and their enrichment analysis
Most of the DEGs were regulated by hsa-miRNA-125a-5p, has-miRNA-125b-5p, hsa-miRNA-19a-3p and hsa-miRNA-19b-3p (Supplementary Table 1). The hsa-miRNA-125a-5p influenced DEGs played a role in transcription regulation. The hsa-miRNA-19a-3p targeted DEGs influenced cell proliferation and metastasis. The hsa-miRNA-19a-3p influenced lymph-node metastasis, liver metastasis and advanced stages in colorectal cancer. Gene enrichment analysis indicated that significant biological pathway for DEGs targeted against miRNAs were included, proteoglycan syndecan-mediated signaling, glypican pathway, beta 1 integrin cell surface interaction, VEGFR signaling and Erb-B receptor signaling network. Their significantly molecular functions were such as transcription factor activity, transcription regulator activity and transporter activity.

Discussion
According to data analyses, it was determined that this data had potential value for this study. In the present study, by using a combined approach of microarray data analysis-bioinformatics tools, we understand a total number of 249 DEGs, including 104 downregulated genes and 145 upregulated genes in medulloblastoma compared to normal brain samples. Our findings revealed that CDK1, CCNB1, CDK2, and WEE1 had strong interaction in the PPI network. Moreover, the cell cycle, p53 signaling pathway, and the FoxO signaling pathway were significantly enriched by DEGs.
Previous studies indicated that CDK1 plays an important role in neuroblastoma, an embryonal tumor,   (Ringer et al., 2011). In addition, another study shown siRNA, MK-1775, could mediate knockdown of WEE1, a kinase known to participate in the G2/M cell cycle and checkpoint and DNA replication during S phase, to reduce medulloblastoma cell proliferation (Hirai et al., 2009). Furthermore, Matheson et al., (2016. indicated that inhibition of WEE1 by AZD1775 can induce apoptosis in medulloblastoma tumor cells. WEE1 plays a critical role in cell proliferation and growth and increased tumor progression and metastasis in medulloblastoma. The other experiment determined that WEE1 can be   a potential therapeutic target for medulloblastoma (Moreira et al., 2016). Several studies have shown that both expressions of MYC and CCNB1 can be a strong prognosis biomarker to predict relapse in patients with medulloblastoma (Tamayo et al., 2011). Guessous et al., (2008) suggested that CDK2 kinase activity was involved in neuronal differentiation, cell proliferation and cell cycle progression in medulloblastoma. Therefore, our findings are consistent with previous results and introduce WEE1, CDK1, CDK2, and CCNB1 as key genes that associated with the development of medulloblastoma.
In the present study, it has been shown that significant amount of DEGs are associated with cell cycle and p53 signaling pathways. The BUB1 is a part of the cell cycle pathway and can regulate APC by controlling cell cycle checkpoint in medulloblastoma. Furthermore, miR-631 was downregulated in medulloblastoma, its downregulation may lead to cancer cell proliferation and invasion by regulating its target gene VEGFA in MB. Besides, FoxO is one of the important pathways in medulloblastoma, and absence of FoxO genes can induce cell cycle and cell proliferation. Srivastava et al., (2009) proposed that FOXO1, one of the FoxO genes, phosphorylated by Akt that results in cytoplasmic sequestration of FOXO1 and inhibition of cell death in the medulloblastoma cancer cell. This data shows that FOXO1 can play a role in medulloblastoma cell growth. Evidence showed that there is a cross-talk with p53 signaling pathways, therefore, FoxO can also play an inhibitory role in cell cycle and apoptosis. Since the cell cycle plays an important and mediatory role in the development of medulloblastoma (Dimitrova and Arcaro, 2015). Therefore, our results are consistent with those previous studies, and the cell cycle and p53 signaling pathway have a pivotal role in medulloblastoma progression and development.
In addition, based on this study, hsa-miRNA-125a-5p, h a s -m i R N A -1 2 5 b -5 p , h s a -m i R N A -1 9 a -3 p , and hsa-miRNA-19b-3p were frequently seen among DEGs that enriched in proteoglycan syndecan-mediated signaling, glypican pathway, beta 1 integrin cell surface interaction, VEGFR signaling, and Erb-B receptor signaling network. Proteoglycan syndecan-mediated signaling pathways play a role in cancer cell proliferation, invasion, and metastasis. Park et al., (2002) indicated that anti-sense of syndecan-2 can induce cell cycle arrest in G0/G1 with enhanced expression of p21, p27, and p53 in colon cancer. Some studies have shown that p53 mutation can upregulate VEGF and VEGFR expression (Wheler et al., 2016). Moreover, our studies are in line with previous research and approve them for medulloblastoma. Therefore, these miRNAs can arrest the cell cycle by inhibiting critical genes that play a role in essential pathways in the development of medulloblastoma.
In conclusion, WEE1, CDK1, CDK2, and CCNB1 may be a critical gene for medulloblastoma. In addition, cell cycle and p53 signaling pathway play a pivotal role in the invasion and progression of medulloblastoma from normal samples. Therefore, these pathways and genes can be a potential therapeutic target for target therapy in medulloblastoma. We also introduce several miRNAs that can use to target these critical genes and pathways.

Statement conflict of Interest
There is no conflict of interest.