Hub Genes and Key Pathway Identication in Wilms Tumor Based on Bioinformatics Analysis

Wilms tumor (WT) is a childhood kidney cancer with unknown etiology. Gene expression analysis has become very essential in WT. Thus, we performed an integrated analysis of gene expression data to identify new molecular mechanisms and key functional genes in WT. Gene expression (GSE60850) dataset was downloaded from Gene Expression Omnibus. Differentially expressed genes (DEGs) were identied using limma. Pathway and Gene Ontology (GO) enrichment analyses were performed for the DEGs by ToppGene database. Then, protein–protein interaction (PPI) networks and modules were established by the Mentha database and PEWCC1, and visualized by Cytoscape software. Target gene miRNA regulatory network and target gene - TF regulatory network were established by the Network Analyst database and visualized by Cytoscape software. Finally, survival analysis, expression analysis, stage analysis, mutation analysis, immunohistochemical (IHC) analysis, receiver operating characteristic (ROC), reverse transcription polymerase chain reaction (RT-PCR) and immune inltration analysis of hub genes was performed. We identied 988 DEGs ultimately including 502 up regulated genes and 486 down regulated genes. Pathway and GO enrichment analysis revealed that DEGs were mainly enriched in D-myo-inositol (3,4,5,6)-tetrakisphosphate biosynthesis, platelet activation, cholesterol biosynthesis III, and complement, coagulation cascades, embryo development, cell surface, DNA-binding transcription factor activity, carboxylic acid metabolic process, extracellular space and signaling receptor binding. FN1, AURKA, TRIM41, NFKBIA, TXNDC5, SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 were ltrated as the hub genes. These identied DEGs and hub genes facilitate our knowledge of the underlying molecular mechanism of WT and have the potential to be used as diagnostic and prognostic biomarkers or therapeutic targets for WT. ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins, pathways in cancer, innate immune system and ATP binding. cytokine signaling in immune system, signaling receptor binding, molecular Wnt


Introduction
Wilms tumor (WT) is the rare diagnosed pediatric tumor worldwide and is named as nephroblastoma [1].
WT is form of kidney cancer that mostly advances in children under age under 10 years [2]. Because of routine early screening and recent advances in treatment techniques, long-term survival rates have upgraded [3]. However, in developing countries, most WT patients are diagnosed at an end stage, with poor prognosis [4]. Therefore, further studies should still be emphasized for the early diagnoses, prognosis and targeted therapy of WT.
Genetic aberrations and its related pathways have been reported to be signi cant factors contributing to the progression of WT. Genes such as IGF2 [5], WT1 [6], RASSF1A [7], PAF1 [8], and DROSHA and DICER1 [9] as well as signaling pathways such as WNT/β-catenin pathway [10], IGF signaling pathway [11], S1P/S1P1 signaling pathway [12], PTEN/PI3K/AKT signaling pathway [13] and VEGF-C/VEGFR-2 signaling pathway [14] were responsible for pathogenesis of WT. Despite improvement and progress in WT diagnosis, prognosis and treatment, the underlying WT molecular mechanisms are not entirely clear and novel diagnosis, prognosis and treatment options are still needed for more effective control of WT development.
Gene expression pro le analysis is a high-throughput method for detecting messenger RNA expression in various cancer tissues or cell samples. By analyzing the different gene expression between cancer patients and normal controls, an improved understanding of the molecular mechanism of a various tumors can be obtained, facilitating the identi cation of the potential key genes and pathways for diagnostics markers, prognostics markers and targeted therapy [15][16].
The current study aimed to explore the molecular pathogenesis of WT by a computational bioinformatics analysis of gene expression. Gene expression data from the Gene Expression Omnibus (GEO) database was extracted, and differentially expressed genes (DEGs) between WT and normal samples were identi ed. The possible functions of the DEGs were predicted using pathway and gene ontology (GO) enrichment analysis. Furthermore, protein-protein interaction (PPI) networks were constructed using mentha PPI database, and visualized and module analysis was conducted using Cytoscape software to search for essential hub genes that may be associated in the progression of WT. Dysregulation of microRNAs (miRNAs) and transcription factors (TFs) have been indicated to be associated with the pathogenesis of WT, the WT speci c regulatory networks of target gene and miRNA, and target gene and TFs were constructed. Validation of the hub genes was performed to screen genes with prognostic and diagnostics signi cance in WT.

Microarray data
Human gene expression microarray data of WT samples (n = 36) and normal samples (n = 36) were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) with an accession ID of GSE60850. The platform of GSE60850 is GPL19130 Breakthrough Human 17K 2.1.2.

Data preprocessing
The raw data in GSE60850 were preprocessed using limma [17], an R software package and it implemented background correcting, quantile normalization and expression calculation automatically.
Then, probe values were translated to gene-symbol values based on message associated in microarray platform, and probes without proper gene-symbols were excluded.

Differentially Expressed Genes
Based on the gene expression microarray data, DEGs between WT samples and normal samples were identi ed using limma [17], an R software package. The corresponding p-values were calculated using ttest provided by limma. The genes met the criteria of p-value<0.05 and |log2 fold change (FC)|≥1.22 for up regulated genes and |log2 fold change (FC)|≥ -1.39 for down regulated genes were de ned as signi cant DEGs between the two groups.
Construction of target gene -TF regulatory network Experimentally-validated target genes and their TFs were screened in one TF database ChEA database (http://amp.pharm.mssm.edu/lib/chea.jsp) [44]. TFs that have a regulatory relationship with the target genes in the constructed network were identi ed. The NetworkAnalyst (https://www.networkanalyst.ca/) [41] online tool was used to predict TF-regulating genes in the network. Cytoscape (version: 3.7.2) [34], an open-source platform for visualizing complex networks, was used to visualize target genes -TF regulatory network.

Validation of hub genes
The survival probability study was implemented using Kaplan-Meier method to compare overall survival curves between high and low expression gene groups UALCAN (https://ualcan.path.uab.edu/index.html) online dataset [45], which is a user-friendly, interactive web resource for the analysis of cancer transcriptome data.. P<0.05 was considered to indicate a statistically signi cant difference. The expression analysis and stage analysis of hub genes were analyzed using UALCAN online dataset [45]. The mutation frequencies of up and down hub genes were inquired in cBioportal online database (http://www.cbioportal.org/) [46]. In addition, up and down regulated hub genes were further validated for their prognostic values (immunohistochemical (IHC) analysis in normal and cancer tissue) using The Cancer Genome Atlas database (https://www.proteinatlas.org/) [47]. Receiver operating characteristic (ROC) analyses are generally used to check out the conduct of disease diagnosis and prognosis. The area under the curve (AUC) was used to demonstrate the accuracy of an individual gene for predicting recurrence using R package"pROC" [48]. Reverse transcription polymerase chain reaction (RT-PCR) was carried out for validation of up and down regulated hub genes. Total RNA was extracted from the WT tissue sample and normal kidney tissue samples using TRI Reagent® (Sigma, USA) according to the manufacturer's protocol. A RNA was reverse transcribed into cDNA using FastQuant RT kit (with gDNase; Tiangen Biotech Co., Ltd.), according to the manufacturer's protocol. The primer sequences (Genewiz, Inc.) used for RT-PCR are listed in Table 1. The mRNA expression levels of hub genes were measured by Real time-PCR using the QuantStudio 7 Flex real-time PCR system (Thermo Fisher Scienti c, Waltham, MA, USA) . The following reaction conditions were used for RT-PCR: Initial denaturation at 95˚C for 3 min followed by 40 cycles of denaturation at 95˚C for 10 sec and annealing and elongation at 60˚C for 30 sec.
The relative expression levels of up and down regulated hub genes were determined using the 2 -ΔΔCt method [49] and normalized to the internal reference gene, β-actin. Immune in ltration analysis was performed using the TIMER (https://cistrome.shinyapps.io/timer/) [50] is a RNA-Seq expression pro ling database from The Cancer Genome Atlas (TCGA) portal for up and down regulated hub genes, which is used to check the immune in ltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) across WT.

Identi cation of DEGs
After data, including 36 WT samples and 36 normal samples, was downloaded from GEO database and preprocessed. The results before and after normalization are shown Fig. 1A and Fig. 1B. 988 DEGs, including 486 up genes and 502 down genes were identi ed using limma packages on the basis of the cut off criteria (P<0.05 and |log2 fold change (FC)| > 1.39 for up regulated genes and |log2 fold change (FC)| < -1.22 for down regulated genes) in WT samples compared with normal samples ( Table 1). The volcano plot showed the up regulated and down regulated genes in dataset GSE60850 is shown in Fig. 2. The details of up and down regulated gene expression heat map are shown in Fig. 3 and Fig. 4.

Pathway enrichment analysis of DEGs
Pathway enrichment analysis of the DEGs (up and down regulated genes) was performed using ToppGene. Pathways were identi ed for the up regulated genes, including the cholesterol biosynthesis III (via desmosterol), superpathway of methionine degradation, complement and coagulation cascades, ECM-receptor interaction, FOXA1 transcription factor network, direct p53 effectors, hemostasis, extracellular matrix organization, phenylalanine tyrosine and tryptophan biosynthesis, tyrosine metabolism, ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins, genes encoding enzymes and their regulators involved in the remodeling of the extracellular matrix, plasminogen activating cascade, blood coagulation, altered lipoprotein metabolic, gluconeogenesis pathway, phenylalanine and tyrosine metabolism. Similarly, pathways were identi ed for the up regulated genes including the D-myo-inositol (3,4,5,6)-tetrakisphosphate biosynthesis, 1D-myoinositol hexakisphosphate biosynthesis V (from Ins(1,3,4)P3), platelet activation, protein digestion and absorption, endothelins, alpha-synucleinsignaling, extracellular matrix organization, degradation of the extracellular matrix, MAP kinase kinase activity, glycolysis, gluconeogenesis, ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins, ensemble of genes encoding core extracellular matrix including ECM glycoproteins, collagens and proteoglycans, Wnt signaling pathway, integrin signalling pathway, activinsignalin, parkinson disease, quinapril pathway and diltiazem pathway. The detailed results of pathway enrichment analysis for up and down regulated genes are presented in Table 2 and Table 3.

GO enrichment analysis of DEGs
All up and down regulated genes were uploaded to the ToppGene software to identify GO function. GO enrichment analysis results for up and down regulated genes are presented in Table 4 and Table 5. For biological processes (BP), the top GO terms of up and down regulated genes were signi cantly enriched in carboxylic acid metabolic process, oxoacid metabolic process, embryo development and animal organ morphogenesis, were included. For cell component (CC), top GO terms of up and down regulated genes were signi cantly enriched in cell surface, endoplasmic reticulum, neuron projection and neuron part. For molecular function (MF), the top GO terms of up and down regulated genes were signi cantly enriched in signaling receptor binding, identical protein binding, DNA-binding transcription factor activity and calcium ion binding.

PPI network construction and module analysis
The Mentha PPI database was used to construct PPI networks. The PPI network of the up regulated genes is illustrated in Fig. 5 with 7649 nodes and 17236 edges. The topology analysis (higest node degree distribution, betweenness centrality, stress centrality, closeness centrality and lowest clustring coe cient) for up regulated genes showed that ESR1, FN1, AURKA, SMURF1, PDK1, NANOG, SLC25A5, NUDT21, KCNQ3, ADM, CEL, CXCL3 and GABRA5 were the hub genes (Table. 6) and statistical results in scatter plot for node degree distribution, betweenness centrality, stress centrality, closeness centrality and clustring coe cient are shown in Fig. 6A -6E. These identi ed hug genes were enriched in neuron part, ECM-receptor interaction, metabolism of proteins, negative regulation of response to stimulus, carboxylic acid metabolic process, response to oxygen-containing compound, programmed cell death, identical protein binding, cell surface, signaling receptor binding, metabolic pathways, ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins, and regulation of response to stress. Similarly, PPI network of the down regulated genes is illustrated in Fig. 7 with 7691 nodes and 16050 edges. The topology analysis (higest node degree distribution, betweenness centrality, stress centrality, closeness centrality and lowest clustring coe cient) for down regulated genes showed that VCAM1, DDIT4L, TCF4, PLK1, RB1, MEOX2, SYK, PLXDC1, TCF7L2, MAPK10, MAGI1 and MRPL15 were the hub genes (Table. 6) and statistical results in scatter plot for node degree distribution, betweenness centrality, stress centrality, closeness centrality and clustring coe cient are shown in Fig. 8A -8E. These identi ed hug genes were enriched in cell adhesion molecules (CAMs), regulation of Wnt-mediated beta catenin signaling and target gene transcription, FoxO family signaling, regulation of retinoblastoma protein, embryo development, animal organ morphogenesis, ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins, pathways in cancer, innate immune system and ATP binding.
Based on the hub genes (up and down regulated) from the PPI module, pathway and GO terms for further analysis. We chose the four most signi cant modules (up regulated genes) for further analysis (Fig.9). Hub genes in these PPI modules were mainly enriched in the ECM-receptor interaction, metabolism of proteins, negative regulation of response to stimulus, programmed cell death, response to endogenous stimulus, neuron part, protein-containing complex binding, focal adhesion, proteoglycans in cancer, regulation of cell differentiation, signaling receptor binding, enzyme regulator activity, carboxylic acid metabolic process, regulation of response to stress, metabolism of amino acids and derivatives, metabolism of lipids and lipoproteins, cell motility and enzyme binding. Finally, we chose the four most signi cant modules (down regulated genes) for further analysis (Fig.10). Module 17 consisted of 145 nodes and 186 edges, module 24 consisted of 122 nodes and 188 edges, module 34 consisted of 100 nodes and 117 edges, and module 40 consisted of 93 nodes and 95 edges. Hub genes in these PPI modules were mainly enriched in the pathways including DNA-binding transcription factor activity, transcription regulatory region DNA binding, sequence-speci c DNA binding, pathways in cancer, extracellular matrix organization, hemostasis, innate immune system, PDGFR-beta signaling pathway, cytokine signaling in immune system, signaling receptor binding, molecular function regulator, Wnt signaling pathway, embryo development, neurogenesis, regulation of cell differentiation, positive regulation of multicellular organismal process and cell surface.

Construction of target gene -miRNA regulatory network
Based on the interaction information of target genes and miRNAs in corresponding miRNA databases, the integrated regulatory network of target genes (up and down regulated) and relevant miRNAs were constructed ( Fig. 11 and Fig. 12). We found that up regulated target genes such as CCND1 can be targeted by 197 miRNAs (ex, hsa-mir-2392), SCD can be targeted by 167 miRNAs (ex, hsa-mir-1269a), PTP4A1 can be targeted by 132 miRNAs (ex, hsa-mir-6731-5p), LDLR can be targeted by 123 miRNAs (ex, hsa-mir-4295) and RRM2 can be targeted by 102 miRNAs (ex, hsa-mir-4458) are listed in Table 7. These identi ed target genes were enriched in focal adhesion, PPAR signaling pathway, cell motility, organic substance catabolic process and superpathway of purine nucleotide salvage. Similarly, we found that down regulated target genes such as ZNF703 can be targeted by 115 miRNAs (ex, hsa-mir-3938), ENPP5 can be targeted by 114 miRNAs (ex, hsa-mir-4768-3p), MYLIP can be targeted by 113 miRNAs (ex, hsa-mir-552-5p), ENAH can be targeted by 92 miRNAs (ex, hsa-mir-4282) and ZBTB20 can be targeted by 85 miRNAs (ex, hsa-mir-4282) are listed in Table 7. These identi ed target genes were enriched in regulation of multicellular organismal development, integral component of plasma membrane, adaptive immune system, axon guidance and positive regulation of multicellular organismal process.

Construction of target gene -TF regulatory network
Based on the interaction information of target genes and TFs in corresponding TF database, the integrated regulatory network of target genes (up and down regulated) and relevant TFs were constructed ( Fig. 13 and Fig. 14). We found that up regulated target genes such as MAGEC2 can be targeted by 207 TFs (ex, SOX2), TSPAN7 can be targeted by 173 TFs (ex, MYC), ESR1 can be targeted by 172 TFs (ex, HNF4A), PCSK6 can be targeted by 166 TFs (ex, EGR1) and LDLR can be targeted by 145 TFs (ex, TP63) are listed in Table 8. These identi ed target genes were enriched in organic substance catabolic process, intrinsic component of plasma membrane, neuron part, golgi apparatus and identical protein binding.Similarly, we found that down regulated target genes such as PLEKHO1 can be targeted by 184 TFs (ex, SOX2), CACHD1 can be targeted by 151 TFs (ex, AR), CASD1 can be targeted by 139 TFs (ex, NANOG), GLIS3 can be targeted by 132 TFs (ex, STAT3) and AFF3 can be targeted by 130 TFs (ex, TP53) are listed in Table 8. These identi ed target genes were enriched in cell projection part, regulation of transcription by RNA polymerase II and DNA-binding transcription factor activity.

Validation of hub genes
Te overall survival rates of patients with high expression of UCHL1, FN1, AURKA, TRIM41 and TXNDC5 were all signi cantly lower than those of patients with low/medium expression (Fig. 15), while overall survival rates of patients with low expression of SIN3A, MAGI1, GPRASP2, FXYD6 and NFKBIA were all signi cantly lower than those of patients with high expression (Fig. 16). The box plots (expression analysis) showed that the expression levels of FN1, AURKA, TRIM41, NFKBIA and TXNDC5 were signi cantly higher in primary tumor than those in the normal kidney for WT patients from TCGA ( Fig.  17A -17E), while the expression levels of SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 were signi cantly lower in primary tumor than those in the normal kidney for WT patients from TCGA ( Fig. 17F -17J). The box plot suggested (stage analysis) that the high expression level of FN1, AURKA, TRIM41, NFKBIA and TXNDC5 show signi cant distance in different pathological stages in KT compared to normal ( Fig. 18A -18E), while low expression level of SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 show signi cant distance in different pathological stages in KT compared to normal ( Fig. 18F -18J). Up and down regulated hub genes' alteration statuses in TCGA WT patients were analyzed using the CbioPortal database. FN1 altered (2%), and missense mutation, truncating mutation, ampli cation and deep dilation were the main type. AURKA altered (0%). TRIM41 altered (8%), and missense mutation and ampli cation were the main type. NFKBIA altered (0.3%), and ampli cation was the main type. TXNDC5 altered (0.7%), and missense mutation and truncating mutation were the main type. SIN3A altered (0.3%), and missense mutation was the main type. MAGI1 altered (2.8%), and inframe mutation, ampli cation and deep dilation was the main type. GPRASP2 altered (2%), and truncating mutation and ampli cation were the main type. UCHL1 altered (0%). FXYD6 altered (0.3%), and ampli cation was the main type. The frequencies of alteration of each hub gene are shown in Fig. 19. The Human Protein Atlas database, which indicated the expression level of FN1, AURKA, TRIM41, NFKBIA and TXNDC5 were higher in WT tissue compared to normal kidney tissues (Fig. 20A-20E), while expression level of SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 were lower in WT tissue compared to normal kidney tissues (Fig. 20F-20J). The ROC curve de ned an optimal threshold to predict the recurrence risk of WT, and the AUC values of were signi cantly lower in WT tissues compared with normal kidney tissues ( Fig. 22F -22J). The Immune in ltration analysis of up and down hub genes from the TIMER was investigated using TCGA database. The results demonstrated that the higher expression level of FN1, AURKA, TRIM41, NFKBIA and TXNDC5 were all negatively associated with tumor purity (Fig. 23A -23E), while lower expression level of SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 positively associated with tumor purity (Fig. 23F -23J).
PPI network was constructed and analyzed for up regulated genes. AURKA was important for pathogenesis WT [259]. SMURF1 was responsible for invasion of breast cancer cells [260], but this gene may be liable for invasion of WT cells. NUDT21 was involved in proliferation of glioblastoma cells [261], but this gene may be associated with proliferation WT cells. Our study found that NANOG (nanoghomeobox), SLC25A5 and KCNQ3 are up regulated in WT and has potential as a novel diagnostic and prognostic biomarker, and therapeutic target. Similarly, PPI network was constructed and analyzed for down regulated genes. PLK1 was associated with proliferation of kidney cancer cells [262], but this gene may be liable for proliferation of WT cells. Low expression of MAGI1 was linked with progression of kidney cancer [263], but decrease expression of this gene may be responsible for pathogenesis of WT.
Our study found that DDIT4L and MRPL15 are down regulated in WT and has potential as a novel diagnostic and prognostic biomarker, and therapeutic target.
Module analysis was performed for up regulated genes. Genes such as IGF2BP1 [264] and PIR (Pirin (iron-binding nuclear protein)) [265] were linked with invasion of various cancer cells types, bur these genes may be involved in invasion of WT cells. Over expression of CCND1 was involved in pathogenesis of breast cancer [266], but high expression of this gene may be linked with progression of WT. Our study found that APRT, HBZ, EIF2S1, CUL7 and TKT are up regulated in WT and has potential as a novel diagnostic and prognostic biomarker, and therapeutic target. Similarly, module analysis was performed for down regulated genes. FANCC (fanconianemia, complementation group C) was important for advancement of WT [267].
Target gene -miRNA network was constructed and analyzed for up regulated genes. PTP4A1 was important for invasion of breast cancer cells [268], but this gene may be linked with invasion of WT cells. High expression RRM2 of was involved in advancement of cervical cancer [269], but elevated expression of this gene may be associated with development of WT. Similarly, target gene -miRNA network was constructed and analyzed for down regulated genes. ZNF703 was liable for invasion of colorectal cancer cells [270], but this gene may be responsible for invasion of WT cells.
Target gene -TF network was constructed and analyzed for up regulated genes. MAGEC2 was linked with invasion of breast cancer cells [271], but this gene may be involved in invasion of WT cells. Similarly, target gene -TF network was constructed and analyzed for down regulated genes. Methylation inactivation of tumor suppressor PLEKHO1 was responsible for advancement of gastric cancer [272], but loss of this gene may be important for pathogenesis of WT. Our study found that CACHD1 and CASD1 are down regulated in WT and has potential as a novel diagnostic and prognostic biomarker, and therapeutic target.
In the current investigation, the DEGs between WT and normal tissue samples in the GSE60850 dataset were determined, and the up and down regulated hub genes among the DEGs were demonstrated to be associated with the prognosis and diagonsis of patients with WT. Furthermore, FN1, AURKA, TRIM41, NFKBIA, TXNDC5, SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 were identi ed as possible candidate biomarkers for patients with WT. High FN1, AURKA, TRIM41, NFKBIA, TXNDC5 mRNA expression levels and low SIN3A, MAGI1, GPRASP2, UCHL1 and FXYD6 mRNA expression levels were validated by TCGA database, human protein atlas database and subsequent ROC analysis and RT-qPCR analysis, which may preliminarily discover the pathophysiological role of these hub genes in WT at the molecular level.
In conclusion, 988 DEGs and 10 hub genes were identi ed as potential diagnostic or prognostic biomarkers of WT. The current investigation identi ed several genes which had not been already associated with WT and implemented evidence that these genes were associated with this disease.
Encourage examines are recommended to authenticate these results and to more precisely analyze the associations between these genes and WT. Overall, the current investigation highlights possibly new targets for more individualized treatment of patients with WT. Declarations Acknowledgement I thank Richard Dafydd Williams, UCL Institute of Child Health, Developmental Biology and Cancer, 30 Guilford Street, London, United Kingdom, very much, the author who deposited their microarray dataset, GSE60850, into the public GEO database.

Con ict of interest
The authors declare that they have no con ict of interest.

Ethical approval
This article does not contain any studies with human participants or animals performed by any of the authors.

Informed consent
No informed consent because this study does not contain human or animals participants.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.   Page 38/43 Volcano plot of differentially expressed genes. Genes with a signi cant change of more than two-fold were selected. Green dot signi cant up regulated genes and red dot signi cant down regulated genes. Modules in PPI network. The red nodes denote the down regulated genes.