Identification of the miRNA-mRNA regulatory pathways and a miR-21-5p based nomogram model in clear cell renal cell carcinoma

Gathering of relevant microarray data

The microarray expression profiling data of miRNAs (GSE16441–GPL8659, GSE71302 and GSE116251) and genes (GSE15641, GSE16441–GPL6480, GSE16449, GSE53757 and GSE68417) were obtained from Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). The details of the microarray datasets are provided in Table 1.

Data preprocessing and DEG/DEM analysis

Within each GEO dataset, using the “Limma” package for R v3.5.1, we selected miRNAs and genes with |log fold change (FC) |>1 and adjusted p-value of <0.05 as DEMs and DEGs. Three Venn diagrams were constructed with the online Venn diagram drawing tool Draw Venn Diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). Genes found in interactions in at least three gene datasets were considered the primary key DEGs, whereas miRNAs appearing in more than two miRNA datasets were considered the key DEMs.

Enrichment analysis of DEGs

Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed in Metascape (Zhou et al., 2019) (http://metascape.org/gp/index.html#/main/step1). Key DEGs were uploaded to Metascape, and related GO and KEGG pathway enrichment analyses were performed. The cut-off criterion was set at a p-value of <0.05, a minimum count of 3 and an enrichment factor (the ratio between the observed counts and the counts expected by chance) >1.5.

Construction of the miRNA-target gene regulatory network

To determine target genes of the key DEMs and the interactions among them, we downloaded data from miRWalk (Dweep et al., 2011) (version 3.0; http://mirwalk.umm.uni-heidelberg.de/). For each miRNA, key DEM–target gene intersection data were downloaded with standards of “min-p-value” = 1; “position” = 3UTR, 5UTR and CDS, and three databases (TargetScan, miRDB and miRTarBase) were used. Those target genes were regarded as candidate genes. The interactions between the candidate genes and the primary key DEGs were denoted key DEGs within the list; genes targeted by at least three miRNAs were screened as target nodes, whereas the corresponding miRNAs were considered source nodes. Cytoscape 3.6.1 was used for visualization of the miRNA–mRNA regulatory network.

Further validation by TCGA data

On the basis of the expression data and clinical data from the “TCGA-KIRC” dataset, we constructed Kaplan–Meier plots with the “survival” package in R v3.5.1. We further selected key DEMs that could be considered independent prognostic factors via univariate and multivariate Cox regression analyses. We then combined the independent prognostic miRNAs with stage and grade to generate a model for better prediction of patient prognosis. Moreover, two-tailed t-tests were performed on the expression levels of these prognostic miRNAs to validate whether our results were reasonable. Box-plots were constructed in GraphPad Prism 7.0.4 for visualization.

Validation by qRT-PCR from clinical specimens of ccRCC

Further validation was performed via qRT-PCR of the four prognostic key DEMs in 15 pairs of tumor and matched normal tissues (n = 30; the clinical and pathologic characteristics of these 15 patients are listed in Table S1.

Detailed procedures of this step are described in our previous article (Chen et al., 2016). The primers for PCR are shown in Table S2. U6 was used as an internal control to normalize the results. Threshold cycle (Ct) values were calculated after each PCR reaction. Each sample was tested in triplicate, and the relative quantification equation (RQ = 2^−ΔΔCt) was used for evaluating relative miRNA expression. The expression data of all 15 pairs of ccRCC samples and matched normal samples are listed in Table S3. This study was approved by Shengjing Hospital Ethics Committee (2017PS012J).

Screening of key DEMs and primary key DEGs according to interactions among GEO datasets

A total of 7,022 DEGs and 131 DEMs (Fig. 2A) were identified. Within these DEMs, miRNAs that appeared in more than two datasets, such as miR-21, miR-142-3p, miR-155-5p and miR-342-5p, were filtered and considered 26 key DEMs (Supplementary file 1). After the miRWalk 3.0 procedure, 3903 candidate genes were obtained. These candidate genes were found to interact with 3506 upregulated DEGs (Fig. 2B) and 3516 downregulated DEGs (Fig. 2C). A total of 167 key upregulated and 140 key downregulated DEGs were obtained.

GO and KEGG enrichment analyses of primary key DEGs

GO analysis results showed that the up-regulated DEGs were significantly enriched in biological processes including response to wounding, cytokine production, epithelial cell proliferation, animal organ regeneration and blood vessel development (Fig. 3A). Downregulated DEGs were significantly enriched in biological processes including kidney development, regulation of neurogenesis, positive regulation of sodium ion transport, nephron development and regulation of phosphatidylinositol 3-kinase signaling (Fig. 3E). For cell components, upregulated DEGs were significantly enriched in extracellular matrix, platelet alpha granule membrane and endoplasmic reticulum lumen (Fig. 3B). Downregulated DEGs were significantly enriched in bicellular tight junction, axon part, apical part of cell and receptor complex (Fig. 3F). Moreover, 16 GO molecular function were over-represented in the upregulated DEGs, including extracellular matrix structural constituent, growth factor activity, organic acid transmembrane transporter activity, cysteine-type endopeptidase inhibitor activity involved in apoptotic process, and RNA polymerase II regulatory region sequence-specific DNA binding (Fig. 3C), whereas the downregulated DEGs were significantly enriched in 18 GO molecular functions, including small GTPase binding, steroid binding, semaphorin receptor binding cell adhesion molecule binding and structural constituent of cytoskeleton (Fig. 3G). KEGG pathway analysis showed that the upregulated DEGs were significantly enriched in the HIF-1 signaling pathway, transcriptional misregulation in cancer and focal adhesion (Fig. 3D). Downregulated DEGs were significantly enriched in tight junction, Rap1 signaling pathway and hematopoietic cell lineage (Fig. 3H).

Construction of an miRNA-mRNA network of primary key DEGs and key DEMs

Target genes were screened with miRTarBase, miRDB and TargetScan. A total of 20 mRNAs, which were targeted by at least three key DEMs, were selected for the construction of the network. For example, BCL6 had connections with miR-124, miR-181c and miR-34a. The miRNA-mRNA network (Fig. 4) was mapped in Cytoscape 3.6.1.

Overall survival analysis of key DEMs and construction of a nomogram model on the basis of TCGA data along with qRT-PCR validation

Kaplan–Meier plots were drawn on the basis of the survival data downloaded from TCGA. Four of 26 key DEMs—miR-21-5p (Fig. 5A), miR-142-3p (Fig. 5B), miR-155-5p (Fig. 5C) and miR-342-5p (Fig. 5D)—were believed to be associated with overall survival (OS) in patients with ccRCC. According to the results, low expression of miR-21-5p and elevated expression of the other three miRNAs in patients with ccRCC might lead to poorer OS.

Further validation using TCGA clinical data was then performed with t-tests. The expression levels of all four miRNAs in tumor tissues were notably different from those in normal tissues (p < 0.001). The expression levels of miR-21-5p (Fig. 6A), miR-142-3p (Fig. 6B), miR-155-5p (Fig. 6C) and miR-342-5p (Fig. 6D) were dramatically higher in tumor tissues. Box-plots were constructed, and the primary data were log2 standardized to facilitate the visualization of differential expression. Only hsa-miR-21-5p had p < 0.05 in both univariate and multivariate analyses (Table S4). We then built a nomogram model (Fig. S1) for better prediction. The areas under the model’s three- and five-year OS receiver operating characteristic curve were 0.759 and 0.714, respectively (Figs. S2A, S2B); the concordance-index (C-index) was 0.717 (95% confidential interval: 0.713–0.721); and the calibration curves (Figs. S3A, S3B) showed high agreement between the predicted and observed OS. Together, these results indicated the accuracy of our nomogram model. The outcomes of qRT-PCR revealed higher expression levels of miR-21-5p, miR-142-3p, miR-155-5p and miR-342-5p in ccRCC specimens (Figs. 7A–7D). These findings were consistent with results based on bioinformatics methods.