RNA-Sequencing, Connectivity Mapping, and Molecular Docking to Investigate Ligand-Protein Binding for Potential Drug Candidates for the Treatment of Wilms Tumor.

BACKGROUND Wilms tumor, or nephroblastoma, is a malignant pediatric embryonal renal tumor that has a poor prognosis. This study aimed to use bioinformatics data, RNA-sequencing, connectivity mapping, molecular docking, and ligand-protein binding to identify potential targets for drug therapy in Wilms tumor. MATERIAL AND METHODS Wilms tumor and non-tumor samples were obtained from high throughput gene expression databases, and differentially expressed genes (DEGs) were analyzed using the voom method in the limma package. The overlapping DEGs were obtained from the intersecting drug target genes using the Connectivity Map (CMap) database, and systemsDock was used for molecular docking. Gene databases were searched for gene expression profiles for complementary analysis, analysis of clinical significance, and prognosis analysis to refine the study. RESULTS From 177 cases of Wilms tumor, there were 648 upregulated genes and 342 down-regulated genes. Gene Ontology (GO) enrichment analysis showed that the identified DEGs that affected the cell cycle. After obtaining 21 candidate drugs, there were seven overlapping genes with 75 drug target genes and DEGs. Molecular docking results showed that relatively high scores were obtained when retinoic acid and the cyclin-dependent kinase inhibitor, alsterpaullone, were docked to the overlapping genes. There were significant standardized mean differences for three overlapping genes, CDK2, MAP4K4, and CRABP2. However, four upregulated overlapping genes, CDK2, MAP4K4, CRABP2, and SIRT1 had no prognostic significance. CONCLUSIONS RNA-sequencing, connectivity mapping, and molecular docking to investigate ligand-protein binding identified retinoic acid and alsterpaullone as potential drug candidates for the treatment of Wilms tumor.

Background

Wilms tumor, or nephroblastoma, is a malignant embryonal tumor of the kidney that most commonly occurs in the neonatal and pediatric population and has typical histological characteristics with undifferentiated blastema tissue, epithelial and mesenchymal elements [1-3]. Wilms tumor represents approximately 90% of pediatric renal tumors and 5% of all the pediatric tumors, and 98% of patients are less than 10 years old. In the UK, 8.1 out of 1,000,000 children were diagnosed with Wilms tumor as new cases annually. For example, in the UK, 90 new cases of Wilms tumor occurred annually, while in the US, the number was higher, with 650 children diagnosed annually [4,5]. In China, Xinjiang is an area with a high incidence of Wilms tumor in children, and most children present with late-stage malignancy. Several research centers in China have reported that the 5-year disease-free survival (DFS) rates of stage II and stage III Wilms tumor were 78% and 78.2%, respectively [6,7]. The overall survival (OS) rate of Wilms tumor is relatively greater, at >90% when compared with other tumors, but there is still a need for further studies [6,7]. Poor prognosis in approximately 25% of patients is associated with high-risk tumors, bilateral lesions, and tumor recurrence. Current treatments for Wilms tumor can have side effects, and up to 25% of survivors suffer from severe chronical health problems at 25 years following diagnosis [8,9]. Therefore, there is a need to continue to study the molecular mechanisms of Wilms tumor and to find molecular targets for future drug therapy.

Several recent studies have investigated the molecular pathogenesis of Wilms tumor to identify relevant genes, to discover new treatment targets and more effective therapy. For example, Wang et al. showed that EGF, CDK1, ENDRA, NGFR, OIP5, NUF2, and CDCA8 might have a role in the development of high-risk Wilms tumor by using weighted gene co-expression network analysis (WGCNA) to identify the hub genes involved in tumor progression [10]. Also, Wilms tumor was shown to be associated with the overexpression of miR-483-5p, mutation of TRIM28, and polymorphisms of miR-423 rs6505162 C>A [11-13]. However, despite these findings, molecular studies on the treatment for Wilms tumor are still limited.

Gene therapy for malignancy has gained increasing importance in clinical research [14]. Drug discovery and development is expensive and takes many years. However, a more economical strategy is to study existing drugs that potentially target cancer-related genes [15]. The Connectivity Map (CMap) database, developed by the Broad Institute in the USA, is an accessible and sophisticated resource that is available to help to establish the relationships between disease characteristics and differentially expressed genes (DEGs) in human disease, and following treatment [16]. Molecular docking is an approach to binding the predicted drugs to particularly selected proteins, which allows investigators to analyze the docking of the ligand to the target, and to understand the chemical and structural basis of the target specificity [14,17-19]. Molecular docking can detect the molecular initiating events (MIEs) that occur in the pathways that are associated with treatment side effects [14,17-19] An increasing number of studies have been conducted based on CMap and molecular docking to identify more effective drugs for the treatment of malignancy, including breast carcinoma, cervical carcinoma, prostate carcinoma, and carcinoma of the lung [20-22]. However, there have been few recent studies that mapped the potential drug targets in Wilms tumor.

Therefore, this study aimed to use bioinformatics data, the RNA-sequencing data resources Genotype-Tissue Expression (GTEx) and the TARGET targeted RNA-sequencing (RNA-Seq) database, connectivity mapping using CMap, molecular docking, and ligand-protein binding to identify potential targets for drug therapy in Wilms tumor.

Material and Methods

Acquisition of differentially expressed genes (DEGs) of Wilms tumor from different databases and enrichment analysis

The TARGET targeted RNA-sequencing (RNA-Seq) database was used that contained detailed data about from childhood tumors, including RNA, whole-genome sequencing (WGS), whole-exome sequencing (WES), copy number variants (CNVs), gene methylation, and clinical data. This study used the open controlled-access data in the TARGET targeted RNA-sequencing (RNA-Seq) database that was acquired from Data Matrix (). Non-tumor controls were obtained from Genotype-Tissue Expression (GTEx) (). The standards were |logFC| >1 and adjPVal <0.01.

The DEGs were identified using the voom method in the limma package. Also, this study included the WEB-based Gene SeT AnaLysis Toolkit (WebGestalt) () and the Metascape gene annotation and analysis resource () to perform Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of 980 DEGs with significant adjPVal for false discovery rate (FDR)-adjusted p-values, to investigate the biological process and molecular mechanism.

The prediction of candidate drugs based on DEGs using the Connectivity Map (CMap) database

DEGs were selected with a significant adjPVal value, and converted to probe number HG133A assigned by Affymetrix (). The CMap database (Broad Institute, Cambridge, MA, USA) () contained 6,100 cases and 1,309 small molecules, which were used to investigate the interactions between small molecule drugs, DEGs and Wilms tumor. Assisted by the annotations, it was possible to choose the candidate drugs according to the gene expression profile of each case. Because CMap could process up to 1,000 pieces of probe information, the probe numbers of the 980 DEGs with significant adjPVal were grouped into the upregulation and down-regulation categories. The data were then input into CMap, and the drugs with a score <−0.9 were deemed as potential drug candidates.

The prediction of target genes and the enrichment analysis

When the drugs with a score <−0.9 were input into PubChem (), and their corresponding simplified molecular-input line-entry system (SMILES) line notation structures were acquired. The SMILE structures were input into the search tool for interactions of chemicals (STITCH) database (), and the results were compared with 980 identified DEGs. The thermal maps were drawn with the drug target genes to illustrate their expression levels in both cancer and non-cancer cohorts.

The Benjamini–Hochberg (BH) procedure was used for the KEGG pathway analysis by the ORA method of WebGestalt and the pathways with an FDR <0.05 were compared with the DEGs-enriched pathways. Further analysis was performed if the DEGs and drug target genes overlapped, and the sequencing data were highly expressed in the Wilms tumor.

The molecular docking using systemsDock

Molecular docking assesses the behavior of the small molecule binding site of the target protein, and this study used the systemsDock website for drug prediction and analysis (). The overlapping DEGs and drug target genes were selected as the drug target genes for Wilms tumor. The upregulated genes in Wilms tumor were docked with the relevant drugs on systemsDock.

Clinical significance and prognosis analysis

Gene databases that were searched included the Gene Expression Omnibus (GEO), ArrayExpress, and Oncomine to identify the importance of overlapping genes and to obtain gene profiles associated with Wilms tumor. The criteria for screening included expression profiles from histological samples associated with the Wilms tumor and available clinical data. Box plots examined the clinical significance of the overlapping genes. The receiver operating characteristic (ROC) curves, and the standardized mean difference (SMD) were analyzed using GraphPad Prism version 8.0 (GraphPad Software, La Jolla, CA, USA) and STATA 12.0 (Stata Corporation, College Station, TX, USA). Of the 177 cases identified from TARGET and GTEx, 126 cases had corresponding prognostic information. A modification used in this study was the analysis of the prognostic information through survival package in R, including Kaplan-Meier curves and hazard ratios (HRs). The final expression value of each gene was obtained by the average, if the gene corresponded to different probes. According to the annotated information of each expression profile, all data not marked as standardized were processed using the log2 ratio. A P-value <0.05 was considered to be statistically significant.

Results

The identification of differentially expressed genes (DEGs) and the enrichment analysis

This study included 177 tumor and non-tumor gene sequencing data of Wilms tumors identified from the Genotype-Tissue Expression (GTEx) and TARGET databases. There were 980 DEGs chosen with significant adjPVal for false discovery rate (FDR)-adjusted p-values, which included 648 upregulated genes and 342 downregulated genes. The Gene Ontology (GO) enrichment analysis included biological process, cellular component, and molecular function. The top 10 enriched pathways are listed in Table 1. Also, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs and other relevant analysis are also shown in Figure 1.

Table 1

Gene Ontology enrichment analysis of differentially expressed genes in Wilms tumor.

Gene Set	Description	Size	Expect	Ratio	P Value	FDR
Biological process
GO: 0007049	Cell cycle	1739	81.91	2.7591	0	0
GO: 0022402	Cell cycle process	1274	60.008	3.0663	0	0
GO: 0051726	Regulation of cell cycle	1106	52.095	2.6106	0	0
GO: 0006259	DNA metabolic process	970	45.689	3.0204	0	0
GO: 0000278	Mitotic cell cycle	927	43.663	3.9163	0	0
GO: 0006974	Cellular response to DNA damage stimulus	806	37.964	2.7131	0	0
GO: 0007346	Regulation of mitotic cell cycle	578	27.225	3.5262	0	0
GO: 0045786	Negative regulation of cell cycle	559	26.33	2.9624	0	0
GO: 0006281	DNA repair	511	24.069	3.2822	0	0
GO: 0006260	DNA replication	268	12.623	5.2284	0	0
Cellular component
GO: 0044430	Cytoskeletal part	1621	60.433	2.1511	0	0
GO: 0015630	Microtubule cytoskeleton	1165	43.433	2.901	0	0
GO: 0005694	Chromosome	1014	37.803	4.8938	0	0
GO: 0044427	Chromosomal part	886	33.031	5.298	0	0
GO: 0005815	Microtubule organizing center	722	26.917	2.9349	0	0
GO: 0000228	Nuclear chromosome	573	21.362	5.1961	0	0
GO: 0044454	Nuclear chromosome part	535	19.945	5.1641	0	0
GO: 0000785	Chromatin	509	18.976	4.6374	0	0
GO: 0000790	Nuclear chromatin	341	12.713	4.4836	0	0
GO: 0000793	Condensed chromosome	223	8.3137	6.7359	0	0
Molecular Function
GO: 0003690	Double-stranded DNA binding	915	41.274	2.5682	0	0
GO: 0003682	Chromatin binding	520	23.456	3.1974	0	0
GO: 1990837	Sequence-specific double-stranded DNA binding	823	37.124	2.3974	9.10E-15	5.70E-12
GO: 0044212	Transcription regulatory region DNA binding	896	40.417	2.2515	1.87E-13	7.99E-11
GO: 0001067	Regulatory region nucleic acid binding	898	40.507	2.2465	2.13E-13	7.99E-11
GO: 0003697	Single-stranded DNA binding	107	4.8266	5.1796	7.64E-12	2.39E-09
GO: 0042393	Histone binding	192	8.6608	3.8103	2.98E-11	8.00E-09
GO: 0043565	Sequence-specific DNA binding	1097	49.484	1.9804	4.15E-11	9.73E-09
GO: 0001012	RNA polymerase II regulatory region DNA binding	735	33.155	2.232	6.25E-11	1.28E-08
GO: 0000976	Transcription regulatory region sequence-specific DNA binding	781	35.23	2.1857	6.80E-11	1.28E-08

GO – Gene Ontology; FDR – false discovery rate.

Figure 1

The bioinformatics analysis of differentially expressed genes (DEGs) in Wilms tumor. (A) Volcano plot. (B) Thermal map. (C) The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. FDR – false discovery rate.

The prediction of drug candidates for the DEGs

The Connectivity Map (CMap) database assisted in the comparison of the DEGs and identified a possible 21 therapeutic targets (score <−0.9) that were considered to be potential therapeutic targets for Wilms tumor (Table 2).

Table 2

Ten potential drugs for Wilms tumor acquired from the Connectivity Map database.

Rank	Batch	CMap name	Dose	Cell	Score	Up	Down	Instance_id
1	713	Menadione	23 μM	PC3	−1.000	−0.261	0.277	4662
2	702	Promazine	12 μM	PC3	−0.998	−0.306	0.230	4308
3	1066	Alsterpaullone	10 μM	MCF7	−0.981	−0.253	0.275	7051
4	502	Resveratrol	10 μM	MCF7	−0.981	−0.323	0.205	958
5	514	Tyrphostin AG-825	25 μM	MCF7	−0.97	−0.278	0.244	1114
6	701	Cloperastine	11 μM	PC3	−0.954	−0.272	0.241	4271
7	688	Fluvoxamine	9 μM	PC3	−0.944	−0.277	0.231	3995
8	701	Fenoprofen	7 μM	PC3	−0.937	−0.272	0.232	4274
9	764	1,4-chrysenequinone	15 μM	PC3	−0.927	−0.259	0.24	7139
10	505	Ionomycin	2 μM	MCF7	−0.916	−0.305	0.188	882
11	502	Quinostatin	10 μM	MCF7	−0.914	−0.297	0.195	973
12	714	Flupentixol	8 μM	PC3	−0.913	−0.296	0.195	6708
13	753	Zoxazolamine	24 μM	PC3	−0.908	−0.270	0.219	6290
14	665	Iopanoic acid	7 μM	HL60	−0.908	−0.26	0.229	2965
15	704	Sulpiride	12 μM	PC3	−0.906	−0.244	0.244	4566
16	772	8-azaguanine	26 μM	MCF7	−0.906	−0.283	0.204	7444
17	762	Nortriptyline	13 μM	PC3	−0.905	−0.294	0.192	7300
18	719	Gliclazide	12 μM	PC3	−0.903	−0.264	0.222	5089
19	1075	GW-8510	10 μM	PC3	−0.901	−0.271	0.214	7085
20	694	Benperidol	10 μM	MCF7	−0.901	−0.277	0.208	4781
21	746	Tretinoin	13 μM	MCF7	−0.901	−0.285	0.200	6243

CMap – the connectivity map.

The prediction of target genes for the drug candidates

The target genes predicted by the search tool for interactions of chemicals (STITCH) database as potential drug candidates were obtained. The overlapping 75 target genes were selected for further analysis (Figure 2A, 2B). Analysis of the overlap with the DEGs identified seven commonly overlapping genes that included CRAPB2, CDK2, SIRT1, CA12, MAP4K4, KCNJ1, and ADNP. The pathway analysis showed that the 75 target genes and 980 DEGs showed no significance in the KEGG pathways (Figure 2C).

Figure 2

The bioinformatics analysis of target genes for the drug candidates. (A) Venn diagram. (B) Thermal map. (C) The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. FDR, false discovery rate.

The docking of drug candidates to the overlapping genes

Of the seven overlapping genes, except CA12, six genes were highly expressed in Wilms tumors. The six upregulated genes in Wilms tumor were then docked to four drugs, retinoic acid, gliclazide, alsterpaullone, and trans-resveratrol by targeting these genes. The docking scores of ADNP and KCNJ1 could not be calculated as they did not have the corresponding numbers in the Protein Data Bank (PDB) database. The docking scores of the other four genes to the four drugs showed relatively high scores when retinoic acid and alsterpaullone were docked to the four genes (Figure 3). The results of docking are shown in three-dimension in Figure 4, and in two-dimensions in Figure 5.

Figure 3

The docking results shown as bar charts. (A) Retinoic acid. (B) Alsterpaullone.

Figure 4

The molecular interactions between drugs and genes are shown in three-dimension. (A) Retinoic acid and CDK2. (B) retinoic acid and MAP4K4. (C) retinoic acid and CRABP2. (D) Retinoic acid and SIRT1. (E) Alsterpaullone and CDK2. (F) Alsterpaullone and MAP4K4. (G) Alsterpaullone and CRABP2. (H) Alsterpaullone and SIRT1.

Figure 5

The molecular interactions between drugs and genes are shown in two-dimensions. (A) Retinoic acid and CDK2. (B) Retinoic acid and MAP4K4. (C) Retinoic acid and CRABP2. (D) Retinoic acid and SIRT1. (E) Alsterpaullone and CDK2. (F) Alsterpaullone and MAP4K4. (G) Alsterpaullone and CRABP2. (H) Alsterpaullone and SIRT1.

Clinical significance and prognostic significance of the four overlapping genes, CRAPB2, CDK2, SIRT1, and MAP4K4

The high docking scores between the two drugs and the four overlapping genes, CRAPB2, CDK2, SIRT1, and MAP4K4, showed the importance of the analysis. Four data series were screened that may be used for standardized mean difference (SMD) analysis, including GSE11024, GSE4530, GSE73209, and GSE2712. The results of the SMD are shown in Table 3. The clinical significance and the overlapping genes are shown in Figures 6 and 7. The Kaplan-Meier curves and the hazard ratio (HR) were used to investigate the mRNA expression of the four genes, as shown in Figure 8.

Table 3

The standardized mean difference of four overlapping genes.

Gene	SMD	95% CI	P-value	I2	P (I2)
CDK2	1.429	(0.376–2.482)	0.008	80.90%	0.001
MAP4K4	2.656	(0.717–4.595)	0.007	87.70%	<0.001
CRABP2	1.931	(0.910–2.952)	<0.001	77.50%	0.004
SIRT1	1.537	(−0.165–3.239)	0.077	87.70%	<0.001

SMD – standardized mean difference; CI – confidence interval.

Figure 6

The differential expression and the receiver operating characteristic (ROC) curves of the overlapping genes for molecular docking. (A) Differential expression of CDK2. (B) ROC of CDK2. (C) Differential expression of MAP4K4. (D) ROC of MAP4K4. (E) Differential expression of CRABP2. (F) ROC of CRABP2. (G) Differential expression of SIRT1. (H) ROC of SIRT1.

Figure 7

The forest plots of standardized mean difference for the overlapping genes. (A) CDK2 (P=0.008). (B) MAP4K4 (P=0.007). (C) CRABP2 (P<0.001). (D) SIRT1 (P=0.077).

Figure 8

The prognostic significance of overlapping genes shown by the Kaplan-Meier curves for the four genes, CDK2, MAP4K4, CRABP2, and SIRT1. (A) Kaplan-Meier curves of CDK2. (B) Kaplan-Meier curves of MAP4K4. (C) Kaplan-Meier curves of CRABP2. (D) Kaplan-Meier curves of SIRT1. (E) Hazard ratio (HR) of the four genes.

Discussion

The current treatment for children with Wilms tumor includes chemotherapy, which is given before or after surgery. Despite the increased cure rate for Wilms tumor following multimodal therapy, the recurrence rate remained relatively high (15%) for patients with low-grade Wilms tumor, and the five-year survival rate remains poor [23-26]. Therefore, an increasing number of studies have included the search for drugs to target the relevant genes in Wilms tumor, and have included salidroside and curcumin [27,28]. However, recent studies were based on in vitro or in vivo experiments on a specific drug, to identify or demonstrate the drug effects on Wilms tumor. However, the present study demonstrated the use of the Connectivity Map (CMap) database and molecular docking to conduct bioinformatics analysis of the sequencing data using several databases and showed that retinoic acid and alsterpaullone had potential treatment effects by affecting the differentially expressed genes (DEGs) of Wilms tumor, and identified the four genes, CRABP2, CDK2, MAP4K4, and SIRT1. Also, this study included enrichment analysis of DEGs, multi-factorial survival analysis, and weighted gene co-expression network analysis (WGCNA) analysis, which made this study the study more comprehensive.

Retinoic acid is the bioactive metabolite of Vitamin A, which is a compound that potentially prevents or treats malignancy and which participates in multiple biological processes [29,30]. Retinoic acid can suppress the growth of cancer cells in vitro and induces cell death in breast cancer, neuroblastoma, gastric carcinoma, prostate cancer, pancreatic cancer, and acute promyelocytic leukemia cells [31-33]. In 2005, molecular studies showed that the trans-retinoic acid gene (ATRA) had the potential for the treatment for Wilms tumor. In advanced-stage tumors treated with ATRA, upregulation of several genes was found, including RARRES1, RARRES3, CRABP2, IGFBP3, RARB, RARG, and TRIM22, and the transforming growth factor-β pathway was activated to induce the inhibition of cell growth [34]. The findings from the present study were consistent with those from previous studies and supported that CRABP2 functioned as a DEG in Wilms tumor and was a target gene of retinoic acid. Also, CDK2, MAP4K4, and SIRT1 might serve as targets for retinoic acid, which may be a future direction for studies on the molecular mechanisms of retinoic acid in the treatment of Wilms tumor.

Also, the findings from this study supported that alsterpaullone might have a role in the treatment of Wilms tumor. Alsterpaullone is an inhibitory small molecule that affects the cell cycle to suppress the cytoactivity of tumor cells and causes apoptosis of some tumor cells. Alsterpaullone may also act as a cyclin-dependent kinase (CDK) modulator that can compete with ATP for the binding site of CDK [35,36]. Lahusen et al. first showed that alsterpaullone could induce apoptosis and accelerate the loss of clonogenicity in the Jurkat cell line [37]. Also, previous studies showed that alsterpaullone was a potential treatment for hepatoblastoma, cervical carcinoma, and ovarian carcinoma [35,36,38]. Despite these findings, there have been no previous studies on the role of alsterpaullone in Wilms tumor. This study analyzed the sequencing data by using bioinformatics and showed that from all the upregulated genes in Wilms tumor, the docking CRABP2, CDK2, MAP4K4, and SIRT1 genes to alsterpaullone resulted in higher scores. These findings indicated that alsterpaullone could interfere with the biological processes of Wilms tumor, possibly by targeting these four genes. Therefore, alsterpaullone might have the potential as a new drug for the treatment of Wilms tumor.

The DEGs for Wilms tumor were enriched in three main pathways that included DNA replication, mismatch repair, and the cell cycle. Alsterpaullone had an inhibitory role by influencing the cell cycle. Based on the docking results, alsterpaullone might accelerate apoptosis by interfering with the cell cycle in Wilms tumor. However, further studies are needed to verify these molecular mechanisms. The standardized mean difference (SMD) of three overlapping genes, CDK2, MAP4K4, and CRABP2, were significant in this study. However, the heterogeneity of the SMD analysis was large, and no statistically significant prognostic results were obtained, possibly due to a comparatively small study sample. Therefore, further studies are required to validate these findings.

Conclusions

This study aimed to use bioinformatics data, RNA-sequencing, connectivity mapping, molecular docking, and ligand-protein binding to identify potential targets for drug therapy in Wilms tumor. RNA-sequencing, connectivity mapping, and molecular docking to investigate ligand-protein binding identified retinoic acid and alsterpaullone as potential drug candidates for the treatment of Wilms tumor.