Bioinformatics analysis and identification of genes and pathways involved in patients with Wilms tumor.

Background: Wilms tumor is the most common childhood kidney malignant tumor. However, the genes and signaling pathways associated with the disease remain incompletely understood.
Methods: GSE66405, GSE73209, and GSE11151 were collected from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) were detected using R software. A protein-protein interaction (PPI) network was constructed using the STRING database, and the clustering modules and hub genes were analyzed with the Cytoscape software. Genes functional enrichment analyses were performed using the package "clusterProfiler" in R software, and the gene set enrichment analysis (GSEA) analysis was performed using GSEA v4.1.0 software.
Results: Respectively, 3,092, 620, and 3,567 DEGs were screened in GSE66405, GSE73209, and GSE11151, with a total of 474 common DEGs detected in three expression profiles. For the common DEGs, the top 30 significant results of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses were presented. Furthermore, five modules were found as the most related modules to Wilms tumor. GO term and KEGG pathway enrichment analyses of the genes in all the modules identified 10 GO terms and 5 KEGG pathways as significantly enriched. The top 10 hub DEGs of the PPI network were ALB, CDH1, EGF, AQP2, REN, SLC2A2, SPP1, UMOD, NPHS2, and FOXM1, with ALB identified as the highest degree. GSEA results showed 11 pathways were correlated with ALB expression in GSE66405 and 10 pathways were related to the expression of the ALB gene in GSE73209. Conclusions: Our study revealed robust gene signatures in Wilms tumor. Dysregulations of the signaling pathways were associated with the development and progression of the Wilms tumor, and 10 hub genes may play important roles in its diagnosis and therapy. 2022 Translational Cancer Research. All rights reserved.

Introduction

Wilms tumor, or nephroblastoma, is an embryonic malignant tumor in children which accounts for almost 6% of all pediatric cancers (1,2). The tumor was first reported in 1814 and its detailed pathological description was made by Wilms in 1899. Wilms tumor is mainly seen in children under 5 years of age, while 90% of new cases occur in those less than 3 years of age (3). Although the current cure rate approaches 85% the early stage of the disease lack specific symptoms, which often delays both diagnosis and optimal treatment (4,5). Moreover, about 25% of Wilms tumor patients will face chronic health problems for many years after diagnosis and treatment (6,7). Therefore, it is of great significance to explore the molecular mechanism of Wilms tumor, such as genes and pathways, to improve its diagnosis, prediction, and treatment. In this study we sought to detail specific biomarkers for Wilms tumor by identifying DEGs between it and normal tissues.

Gene chipping, or gene profiling, is an important technology platform in the field of life sciences which has been used for more than 10 years. It is an effective means of screening differentially expressed genes (DEGs), which has the advantages of high throughput and fast measurement (8). Gene chipping has been widely used in gene expression analysis, gene discovery, gene mutation and polymorphism analysis, genome library mapping, disease diagnosis and prediction, drug screening, gene sequencing, and other fields. The development of high throughput microarray hybridization and sequencing technology has generated a great deal of genetic data, and biological information has rapidly expanded into the ocean of data. To adapt to this kind of high throughput gene expression data and the increasing need for data sharing, a variety of databases have been generated, including the national biological information technology center (NCBI) and the high throughput Gene Expression Omnibus (GEO) database. These public databases store the world’s largest high throughput molecular abundance data. Many valuable clues for new research have been discovered through integrating and re-analyzing these databases, and in recent years, many microarray data analysis studies have been published. For example, Zhang et al. (9) used bioinformatics methods to detect the gene expression profile data GSE13601 and revealed the potential mechanism of tongue squamous cell carcinoma, while Long et al. (10) reported that by analyzing the gene expression profile data GSE19804, the DEGs and hub genes enhanced the cognitive of the occurrence and development of lung cancer. The GEO and TCGA databases were used by Zhang et al. (11) to identify key genes and microRNAs in patients with Wilms tumor and Li et al. (12) to screen immune-related prognostic genes in clear cell renal tumors. At present, there is no reliable biomarker for Wilms tumor and the molecular mechanism of its occurrence and development have not been fully clarified.

In this study, the microarray profiles GSE66450, GSE73209, and GSE11151 were selected and collected from the GEO database and analyzed using bioinformatics methods. The purpose of this study was to find common DEGs, analyze them with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) methods, and construct a protein-protein interaction (PPI) network. These data sets may help to improve the understanding of Wilms tumor and have a significant impact on its diagnosis and treatment. The purpose of this study is to find out potential biomarkers for the treatment of Wilms tumor. The combination of multiple chips can greatly increase the accuracy of screening hub genes. We present the following article in accordance with the STREGA reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1847/rc).

Methods

Microarray data information and DEGs identification

The Wilms tumor gene expression profile of GSE66405, GSE73209, and GSE11151 were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), and the raw expression data were downloaded as a series matrix file. The platform information of GSE66405 microarray data was as follows: GPL17077, Agilent-039494 Sure Print G3 Human GE v2 8×60K Microarray 039381 (probe name version) (submission date: 02 Mar 2016) (13). The platform information of GSE73209 microarray data was as follows: GPL10558, Illumina HumanHT-12 V4.0 expression bead chip (Submission date: 18 Sep 2016) (14), and for GSE11151 microarray data was GPL570 (HG-U133_plus_2) Affymetrix Human Genome U133 plus 2.0 Array (submission date: 11 Apr 2008) (15-17). The Affy installation package (library “affy” in R, https://www.r-project.org/) of R software was used for data correction and normalization of each chip. Applying the “limma” package (18) (https://bioconductor.org/packages/limma/) in the R statistics environment, DEGs were then extracted from GSE66405, GSE73209, and GSE11151, respectively. FDR ≤0.05 and |log2 fold change (FC)| ≥1 were set as the thresholds for DEGs. Finally, the “Robust Rank Aggregate” package (19) in R language was applied to identify common DEGs among those from the expression profile data sets. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Construction of PPI network

The online database STRING (available online: http://string-db.org) (20) was employed to analyze the PPI network, and after obtaining the common DEGs with a combined score >0.04 were extracted.

Module analysis and hub genes selection

Subsequently, we used the Molecular Complex Detection (MCODE) tool with a cutoff MCODE score >5 in Cytoscape software (version 3.7.2) analysis of the PPI network. In addition, according to the Closeness method, identify hub genes using the CytoHubba plugin, where genes of the top 10 degrees were considered hub genes in the PPI network.

Functional enrichment analysis

For further functional analysis, the package “clusterProfiler” in R software was used to conduct KEGG and GO enrichment analysis on the common DEGs. Moreover, KEGG and GO pathway analyses of modules genes were conducted with the “clusterprofiler” package in R. For pathway analysis, P<0.05 was seemed as the significant gene functionals, and the results were visualized using the “ggplot2” and “GOplot” packages in R software.

Gene set enrichment analysis (GSEA)

GSEA v4.01.0 software was downloaded from the GSEA home website (http://www.gsea-msigdb.org/gsea/index.jsp) and the software implemented in a Java environment. GSEA was run using c2.cp.kegg.v.7.2.symbols.gmt (Curated) gene sets database, with 1,000 number of permutations, and GSEA FDR <0.05 was considered statistically significant. Multiple GSEA results were visualized using the R software.

Statistical analysis

All statistical analyses were completed in R software. The DEGs were analyzed via the limma package. The Wilcoxon non-parametric test was used for compare differences in two groups, and the correlation analyses were conducted using Spearman’s method. A P value less than 0.05 was considered statistically significant.

Results

Identification of DEGs in Wilms tumor

The GSE66405 data included 28 Wilms tumor and 4 normal tissues, the GSE73209 data had 32 Wilms tumor and 6 normal tissues, and the GSE11151 data included 4 Wilms tissues and 5 normal tissues (). Before data analysis, the R bioconductor “affy” package was used for the normalization of Affymetrix data (). The identification of DEGs of microarray data after normalization was performed via “Limma” package in R. From the expression profile datasets GSE66405, GSE73209, and GSE11151, 3,092 (1,270 up- and 1,822 down-regulated genes), 620 (129 up- and 491 down-regulated genes), and 3,567 (1,698 up- and 1,869 down-regulated genes) DEGs were extracted, respectively. The heatmaps and volcano plots of the DEGs were generated using R language (), with its “Robust Rank Aggreg” package used to analyze the common DEGs from the three profile datasets. A total of 474 common DEGs were identified, including 82 up- and 392 down-regulated genes in Wilms tumor tissues, compared to normal tissues. The package in R language was used to draw the heat map of the top 20 up- and top 20 down-regulated DEGs from the common DEGs ().

Table 1

Characteristics of the individual studies

GSE ID	Platform	Tumor	Normal	Year
GSE66405	GPL17077	28	4	2016
GSE73209	GPL10558	32	6	2016
GSE11151	GPL570	4	5	2008

Figure 1

Normalization of GSE data. (A) GSE66405 data before normalization; (B) GSE73209 data before normalization; (C) GSE11151 data before normalization; (D) GSE66405 data after normalization; (E) GSE73209 data after normalization; (F) GSE11151 data after normalization. The Y-axis stands for the relative mRNA level and the X-axis represents individual sample.

Figure 2

Volcano plots and heat maps of DEGs. (A-C) Volcano plots of DEGs in expression profiling data GSE66405, GSE73209, and GSE11151, respectively. Red and green dots represent up-regulated and down-regulated DEGs respectively, and black dots represent non DEGs. (D-F) Heat maps of DEGs in GSE66405, GSE73209, and GSE11151, respectively. DEGs, differentially expressed genes.

Figure 3

Heat map of 474 common DEGs among Wilms tumor compared with normal controls. DEGs, differentially expressed genes.

PPI network analysis results

To further highlight connections between the common DEGs, the STRING database was used. A total of 474 common DEGs (82 up- and 392 down-regulated genes) were filtered into the DEGs PPI network complex ().

Figure 4

PPI network of 474 common DEGs. PPI, protein-protein interaction; DEGs, differentially expressed genes.

Modules and hub genes analysis

Using the MCODE app in Cytoscape software with a score ≥5, the five most significant modules were aggregated and extracted from the PPI network. Module 1 contained 23 nodes and 244 edges with a score of 22.182; module 2 contained 11 nodes and 50 edges with an MCODE score of 10.000; module 3 contained 11 nodes and 45 edges with an MCODE score of 9.000; module 4 contained eight nodes and 26 edges with a score of 7.429; and module 5 contained 26 nodes and 66 edges with a score of 5.280, respectively (). Subsequently, with the using of cytoHubba we identify the top 10 hub genes of the PPI network, which included ALB, CDH1, EGF, AQP2, REN, SLC2A2, SPP1, UMOD, NPHS2, and FOXM1 (). By the final step, with a degree of 75, ALB was identified as the highest degree protein in the network ().

Figure 5

The top five significant clustering modules that were identified in the PPI networks. (A) Module 1; (B) module 2; (C) module 3; (D) module 4; (E) module 5. PPI, protein-protein interaction.

Figure 6

Hub genes screened from PPI network. (A) The top 10 hub genes explored by CytoHubba; (B) with a degree of 75, ALB was identified as the highest degree in the PPI network. The X-axis represented the edge degree of the hub genes, and the Y-axis represented the hub genes. PPI, protein-protein interaction.

GO and KEGG analysis in Wilms tumor

To acquire a more comprehensive and in-depth understanding of the selected DEGs, KEGG pathways and GO pathways enrichment were analyzed in R software. GO enrichment analysis was classified into three aspects: biological process (BP), molecular function (MF), and cellular component (CC). The top 30 KEGG and GO pathway enrichments are shown (). Moreover, to further investigate the biological functions of the DEGs in all the modules, KEGG and GO analysis were performed. DEGs in the modules were enriched in five KEGG pathways, including mineral absorption, glutathione metabolism, protein digestion and absorption, aldosterone-regulated sodium reabsorption, and platinum drug resistance (). GO enrichment analysis of the genes in modules showed 10 GO terms (GO:0006882, GO:0098754, GO:0055069, GO:0071280, GO:0010043, GO:0010273, GO:1990169, GO:0061687, GO:0097501, GO:0000280) represented significantly ().

Figure 7

KEGG and GO pathway enrichment analysis of common DEGs. (A) The top 30 enriched KEGG pathways are presented; (B) the top 30 GO terms associated with biological process are listed. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; DEGs, differentially expressed genes; BP, biological process; MF, molecular function; CC, cellular component.

Figure 8

Signaling pathway enrichment analysis of DEGs in five modules. (A) Significantly enriched five KEGG pathways of genes in modules; (B) GO enrichment of genes in modules. FC, fold change; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

GSEA enrichment analysis

After determining ALB was the most related protein in the PPI network, we used the GSEA enrichment analysis method to detect the signal pathway related to the ALB gene in the expression profile data. In GSE66405, GSEA results showed apoptosis, the chemokine signaling pathway, FC epsilon RI signaling pathway, FC gamma R mediated phagocytosis, Mark signaling pathway, Neurotrophin signaling pathway, Nod-like receptor signaling pathway, Notch signaling pathway, Rig I like receptor signaling pathway, Toll-like receptor signaling pathway, and VEGF signaling pathway were correlated with ALB expression (). In GSE73209, GSEA results showed allograft rejection, the B cell receptor signaling pathway, Cell adhesion molecules cams, Chemokine signaling pathway, FC epsilon RI signaling pathway, Glycolysis gluconeogenesis, GRAFT versus HOST disease, Regulation of autophagy, RIG I like receptor signaling pathway, and T cell receptor signaling pathway were correlated with ALB expression ().

Figure 9

GSEA enrichment analysis between ALB gene high expression data and ALB low expression data in expression profiling data. (A) GSE66405; (B) GSE73209. KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis.

Discussion

In the present study, three Wilms tumor series matrix files, GSE66405, GSE73209, and GSE11151, were collected from the GEO online database. Moreover, all data were normalized before been analyzed into DEGs. Based on the “limma” package in R language, we extracted 3,092, 620, and 3,567 DEGs from the datasets GSE66405, GSE73209, and GSE11151, respectively. Using RobustRankAggregate packages in R language, we further detected 474 common DEGs between Wilms tumor tissue and normal tissue. As a comprehensive database that stores large amounts of genetic data, GEO has been widely used for bioinformatics analysis, and recently, massive genetic studies of cancer have been published (21-23). Several studies have reported the analysis of Wilms tumor using the GEO database, including Zhang et al. (11) who showed the key genes, miRNA, and mRNA-miRNA regulatory networks may contribute to further comprehending of the molecular mechanisms underlying the progress of the disease. The bioinformatics analysis from the GEO database by Zong et al. showed miR-30d was low-expressed in Wilms tumor and could induce apoptosis and inhibit proliferation, invasion, and migration by mediating Sox4 (24). However, there are few articles in Wilms tumor research utilizing multiple expression profiling data for analysis. Therefore, our study has some notable strengths, by analyzing the three GEO expression profiles.

After obtaining the 474 common DEGs, GO and KEGG analysis were performed via the “clusterProfiler” package, while STRING database analysis revealed the PPI network. Utilizing the MCODE app in Cytoscape software, five modules that might serve an important role in the development of Wilms tumor were detected. Studies have demonstrated the TGF-β/Smad signaling pathway plays a key role in Wilms tumor progression (25), while Li et al. showed S1P/S1P1 signaling stimulates cell migration and invasion in the disease (26). Therefore, our KEGG and GO results suggest these signaling pathways may serve a crucial role in Wilms tumor progression.

Previous studies have demonstrated genetic changes could accelerate the progression of Wilms tumor (27,28). In 1990, Wilms tumor 1 (WT1) was the first gene to be implicated in tumorigenesis and as a strong candidate predisposition gene for the disease (29,30). Recent medical studies confirmed WT1, located at chromosome 11p13, is mutated or absent in 15% of cases (31,32), and in a recent study, WT1 mutations activated cell cycle genes and promoted proliferation of Wilms tumor (33). Some researchers have reported relevant studies that WT1 was high-expressed in many types of cancer, including breast cancer (34), myeloproliferative neoplasms (35), ovarian cancer (36), and endometrial cancer (37). In this study, we constructed a PPI network and identified the top 10 hub genes via cytoHubba as ALB, CDH1, EGF, AQP2, REN, SLC2A2, SPP1, UMOD, NPHS2, and FOXM1. ALB is involved in maintaining human homeostasis and is closely related to many physiological and pathological processes. In pediatric tumors, the ALB gene was detected to be positive in hepatoblastoma and fibrolamellar hepatocellular carcinoma using immunohistochemistry and in situ hybridization, and the sensitivity and specificity of the situ hybridization was 100% (38). Cadherin 1 (CDH1) is related to cell adhesion and is regarded as an invasion-suppressor gene (39), and adecrease in its expression has been shown to be solely responsible for pancreatic cancer metastasis (40). In cardiovascular progenitor cells, WT1 binds to promoters of CDH1 to inhibit its activity and promotes EMT (41). EGF, an epidermal growth factor, stimulates not only cell growth but also cell migration in cancer cells (42,43), and is a potential target gene of WT1 in kidney development (44). Moreover, significant expression of EGF may play a role in promoting proliferation of Wilms tumor, perhaps by an autocrine mechanism (45). Aquaporin 2 (AQP2) is the arginine vasopressin regulated water channel of the kidney collecting duct, and several studies have shown it plays a key role in nephrogenic diabetes insipidus (46-49). Niu et al. have shown that overexpression of AQP2 in pheochromocytoma was positively correlated with tumor size, indicated that AQP2 has the potential to be a diagnostic immunohistochemical biomarker for that disease (50). In 1990, a report have revealed that Renin gene (REN) was high-expressed in Wilms tumor (51), and using in situ hybridization, McKenzie et al. detected REN mutations in 9 of 12 tumor tissues (52). Knockdown of WT1 protein by siRNA significantly increased the expression of renin mRNA, while overexpression reduced REN expression (53). Previous GWAS have identified the rs8192675 in SLC2A2 as being significantly associated with metformin response (54,55). The C allele of rs8192675 in the intron of SLC2A2 was associated with a 0.17% greater metformin-induced reduction in hemoglobin A1c (HbA1c) in 10,577 patients (55), while the expression level of SLC2A2 in hepatocellular carcinoma was higher than part of the family members, suggesting it as a novel prognostic biomarker for that disease (56). The SPP1 gene codes for secreted phosphoprotein-1 (SPP1) protein, also named osteopontin (OPN). SPP1 is highly expressed in a variety of human cancers and involve in many physiological and pathological processes, including cell proliferation, invasion, survival, and tumor metastasis (57,58). In hepatocellular carcinomas, elevated SPP1 levels were able to effectively increase growth and metastasis (59), and were also higher in cervical cancer (60). Higher SPP1 expressed levels were associated with poor disease-free survival (DFS) and overall survival (OS) in cervical cancer patients (61), while in prostate cancer, SPP1 was involved in tumor recurrence and metastasis progression by mediating the Smad4/PTEN pathway (62). Uromodulin (UMOD), also known as Tamm-Horsfall protein, encodes the uromodulin glycoprotein. Under normal physiologic conditions, UMOD is the most abundant protein found in urine, and may be a factor involved in the pathogenic process of kidney disease (63). Remarkably, uromodulin urinary levels in patients after renal transplantation were associated with donor UMOD rs12917707 genotype (64). Mutations in NPHS2 were found worldwide in nephrotic syndrome patients from different countries and ethnic origins (65-67), and more than 100 mutations in NPHS2 have been identified in patients with the disease (68). By integration of cistromic and transcriptomic analyses, NPHS2 was identified as a target gene of WT1 in podocyte differentiation and maintenance (69). FOXM1 (forkhead box M1) is a member of the forkhead box family, and its dysregulation can significantly contribute to tumorigenesis and cancer progression. Overexpression of FOXM1 has been observed in many human cancers including breast cancer (70), pancreatic cancer (71), and hepatocellular carcinoma (72), and was shown to correlate with a poor prognosis in lung cancer patients (73). Knockdown of FOXM1 expression could inhibit cervical cancer cells migration, invasion, and promote sensitivity to cisplatin (74), while its knockdown alone could inhibit cell proliferation and invasion and induce apoptosis in esophageal squamous cell carcinoma (75). The hub genes in the present findings have been published in several articles and are of potential value in the study of Wilms tumor. From the perspective of epigenetics, cancer is the result of result of gene imbalance, methylation, acetylation and phosphorylation are the causes of gene dysregulation (76,77). We seeked potential targets via researched dysregulated genes for the prevention and treatment of nephroblastoma.

One interesting finding in this study was that the hub gene ALB was the most associated protein in the PPI network, compared with other hub genes. Therefore, we performed GSEA analysis in expression profiling data GSE66405 and GSE73209 for examples with high ALB expression and low ALB expression, and the results provide a theoretical basis for the study of ALB genes in Wilms tumor. Unfortunately, the expression profile data GSE11151 was in error due to insufficient examples of the dataset. In addition, the limitations of this study are lack of the validation of expression level and functional research on the main target genes.

Conclusions

Our study displayed robust gene signatures in Wilms tumor and dysregulation of the pathways was closely associated with the development and progression of the disease. The hub genes ALB, CDH1, EGF, AQP2, REN, SLC2A2, SPP1, UMOD, NPHS2, and FOXM1 may play an important role in diagnosing and treating Wilms tumor.

The article’s supplementary files as