Data Collection and Preprocessing
The workflow of the study was represented in Fig. 1. A large amount of RNA-seq data and clinical information were obtained from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. We downloaded the GSE11151 and GSE11024 datasets from the NCBI GEO database (http://www.ncbi.nih.gov/geo). The functions of the datasets in our manuscript were shown in Table 1. In each GEO dataset, we only extracted the samples of the Wilms tumor samples and fetal normal kidney samples for subsequent analysis.
Table 1. Characteristics of the included datasets.
Dataset
|
Country
|
Platforms
|
No. of samples
|
Usage here
|
GSE11151
|
Netherlands
|
GPL570
|
Wilms tumor(n=4), fetal normal kidney (n=2)
|
Identification of hub genes
|
GSE11024
|
USA
|
GPL6671
|
wilms tumor(n=27), fetal normal kidney(n=4)
|
Identification of hub genes
|
Target-WT
|
USA
|
RNA seq
|
DAWT(n=22), FHWT(n=114)
|
Verification of hub genes, Construction and identification of the prognostic model
|
Differential Gene Expression Analysis
The differentially expressed genes between Wilms tumor(WT), and fetal normal kidney tissue from the GEO database were downloaded. The raw data were downloaded as MINiML files. Using the limma package in the R software to study the differentially expressed genes. The adjusted P-value was analyzed to correct the false positive results in GEO datasets. Adjusted P < 0.05 and Log (Fold Change) >1 or Log (Fold Change)< −1 were defined as differentially expressed genes (DEGs). Then, we compared the DEGs identified from GSE11151 with those DEGs from the GSE11024 using a Venn diagram and calculated the overlap coefficient between the two gene sets. The key hub genes were validated in the Target-WT dataset.
The differential gene expression analysis was performed to identify genes that are significantly upregulated or downregulated in DAWT compared to FAWT. We used the R package limma package (version: 3.40.2) to fit a negative binomial model for each gene and estimate the log2 fold change and p-value of differential expression. We considered genes with P < 0.05 were defined as the threshold for the differential expression of mRNAs.
GO and KEGG pathway enrichment analysis
The GO enrichment analysis and KEGG pathway enrichment analysis were used to identify the biological functions and pathways associated with the intersection of differentially expressed genes (DEGs) derived from TARGET and GEO data. The cluster profile package was used for GO enrichment analysis and KEGG pathway enrichment analysis, adjusted p-values for multiple testing using the Benjamini-Hochberg method, and less than 0.05 as significantly enriched. The BP, MF, and CC categories separately and applied a filter of a minimum count of 10 genes per GO term. GO biological processes gene sets, GO cellular components gene sets, and GO molecular functions gene sets were obtained from the Molecular Signatures Database (MSigDB) as reference gene sets[1]. The R package enrich plot was conducted to visualize the enriched GO terms and KEGG pathways using dot plots and bar plots [2]. In addition, we downloaded Gene Set Enrichment Analysis (GSEA) software from the Broad Institute (http://software.broadinstitute.org/gsea/msigdb).
PPI network analysis
A PPI network between DEGs was constructed based on the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db.org) with a confidence level ≥0:400[3]. Then, the PPI file was imported into Cytoscape 3.9.1 (http://cytoscape.org/) to visualize and analyze the PPI network[4]. The hub genes were screened with the MCODE algorithm using default settings.
Survival Analysis
Survival analysis was used to evaluate the prognostic value of the key hub genes in WT. The median expression value of each gene or module eigengene was used as the cutoff to divide the patients into high and low-expression groups. The relationship between the mRNA expression levels of hub genes and the prognosis (OS/PFS) was analyzed through the Kaplan-Meier analysis. The genes or modules with log-rank p-values less than 0.05 were considered as significantly associated with survival. The hazard ratio (HR) and 95% confidence interval (CI) for each gene were calculated by using the log-rank tests and univariate Cox proportional hazards regression.
Construction and identification of the prognostic risk model
To confirm the potential prognosis-related hub genes. The expression matrix integrating the initial genes of the model with patient survival status and survival time was constructed. The LASSO regression algorithm was used for feature selection, and 10-fold cross-validation was used to determine the parameters among which the key genes associated with the patient survival cycle were screened. We calculated the risk score of each patient based on the regression coefficient of the hub genes in the signature and the corresponding expression value of the hub genes. The risk score was calculated using the following formula:
Risk score=expression of Gene1∗α1+expression of Gene2∗α2+…expression of Genen∗αn,
where α represents the regression coefficient of the hub genes in the signature. Based on the median risk score, the patients were divided into high-risk and low-risk groups. Kaplan–Meier survival curve analysis was carried out to compare the OS and PFS between the high-risk group and the low-risk group. A p-value < 0.05 was selected as the significant cutoff value. The time ROC (v 0.4) analysis was used to compare the predictive accuracy of the risk score.
Cell culture
The HEK 293T cells and WiT49 cell line were donated by Dr. Tanpeng Chen. It is a Wilms’ tumor (WT) cell line that is derived from the first-generation xenograft of a human WT lung metastasis. Some differentiation potential is retained by WiT49 cells, displaying the so-called “triphasic” histology when grown in tissue culture plates[5-7]. The WiT49 and 293 cells were cultured as described previously[8]. All cell lines were proven to be mycoplasma negative.
Transfection
Transfection was conducted using Lipofectamine 2000 (Invitrogen, USA) following the manufacturer’s instructions. Si RNA-CCNA1, si-NC, EMCN mimic, and mimic-NC were obtained from Rubibio Company (Guangzhou, China). The CCNA1 targeting siRNA, negative control (NC) siRNA, EMCN mimic, and OE-NC sequences were showed in Table 2.
Table 2 Oligo sequences.
Gene
|
Target sequence
|
Negative Control (si-NC)
|
TTCTCCGAACGTGTCACGT
|
Si RNA-CCNA1-1
|
GTGTTATTCTGGATCAGAAAATG
|
Si RNA-CCNA1-2
|
GACATCTACATGGATGAACTAGA
|
EMCN mimic
|
EMCN-F:ATGGAACTGCTTCAAGTGACCATT
|
EMCN-R:TCAGTTCTTGGTTTTTCCTTGTGCAG
|
Western blotting
Total proteins were extracted from WiT49 and HEK 293T. The protein concentration was detected with the BCA protein assay. Then, 25 μg total protein samples were subjected to SDS-PAGE, and the separated bands were transferred to 0.22 μm PVDF membranes. Protein was blocked for 1 hour with a blocking solution. The membrane was incubated with the primary antibody overnight at 4 °C and with the secondary antibody at room temperature for 1 hour. Finally, the gel was imaged. Anti-EMCN (1:1000 dilution, PA5-21395, Thermo Fisher Scientific), Anti-CCNA1(1:1000 dilution, 13295-1-AP, Proteintech) and anti-α-Tubulin (1:5000 dilution, 11224-1-AP, Proteintech) . α-Tubulin was used as a loading control.
RT-PCR
Cellular RNA was extracted from WT cells (WiT-49) and normal renal epithelial cells (293T). We then synthesized cDNA from the total RNA samples using an M‐MLV reverse transcription kit (M1705, Promega). Quantitative PCR was performed on the resulting cDNA samples using the SYBR Master Mix (DRR041B, TAKARA). The expression of CCNB1 and GAPDH was based on the formula 2^-ΔΔCt. Table 3 lists the primers that were utilized.
Table 3 Primers used for quantitative real time PCR.
Oligonucleotides
|
Sequence (5′–3′)
|
GAPDH-F
|
GTGGGCAAGGTATCCTG
|
GAPDH-R
|
GATTCAGTGTGGTGGGGGAC
|
EMCN-F
|
TGCAGGACTTTCTCCTTTTC
|
EMCN-R
|
ATTTGTTCTGGTGGGTTTGT
|
CCNA1-F
|
GCACACTCAAGTCAGACCTGCA
|
CCNA1-R
|
ATCACATCTGTGCCAAGACTGGA
|
Cell viability and plate clone formation assay
Transfected WiT49 was seeded in 96-well plates at a density of 3,500 cells/well for 24 h. Ten microliters of Cell Counting Kit-8 (CCK-8, Dojindo, Japan) reagent was added into each well at 24, 48, and 72 h posttransfection. After incubating at 37°C and 5% CO2 for 2 h, the absorbance (450 nm) was recorded using a plate reader (Pulang New Technology, Beijing, China).
For the plate clone formation assay, 600 cells per well were seeded in a 6-well plate and cultured for 12 days. The culture medium was changed every 4 days. Then, the cells were fixed in 4% formaldehyde and stained with crystal violet. Cell clones were counted and analyzed.
Transwell assay
Transfected WiT49 was resuspended in a serum-free medium and seeded to the top chamber (Corning, USA) with or without precoating of Matrigel (BD Bioscience, USA). A complete medium with 15% FBS was added to the bottom chamber. After culturing for 24 h at 37°C and 5% CO2, the invaded or migrated cells at the bottom side of the transwell membrane were fixed with 4% paraformaldehyde for 10 min and stained with 0.1% crystal violet (Solarbio, China) at room temperature for 30 min. After washing 3 times with phosphate-buffered saline (PBS), the number of migrated cells was counted under a phrase contrast microscope (Nikon, Japan).
Statistical Analysis
R software (version 4.0.3) for data analysis and visualization. The Kruskal–Wallis test was used for continuous variable data, and the wilcoxon test was used to compare the differences between the two groups. The risk ratio (HR) and 95% confidence interval (CI) were estimated using the survival package of the Cox regression model. Kaplan-Meier method and log-rank test to compare the survival curves between different groups. Two-tailed P values to determine the statistical significance of the differences, and considered them significant when P < 0.05 (*P < 0.05, **P < 0.01, ***P < 0.001).