Clear cell renal cell carcinoma associated microRNA expression signatures identified by an integrated bioinformatics analysis

Dataset collection

Data processing

The quantified probing signals were background corrected using Normexp, with offset value set as 0. The background-subtracted data were normalized using the Quantile algorithm. Averages were derived from all quantile normalized data per miRNA for statistical analysis. Missing data were imputed with the k-nearest neighbors imputation approach (k = 5). We wrote the R scripts to run the processing procedures for all datasets. The mRNA expression data were analyzed using affy package from Bioconductor with RMA method. For the gene annotation, genes that correspond to multiple probes were removed and those with unique probe were retained for further analysis.

Comparison of the outlier detection methods

We compared the performance of five prevailing outlier detection algorithms, COPA (Cancer Outlier Profile Analysis) [21], MOST (Maximum Ordered Subset T-statistics) [22], ORT (Outlier Robust T-statistics) [23], OS (Outlier Sum) [24] and t-test.

All the algorithms were implemented in R scripts written by Lian [22] and Wang [25]. Outliers for each expression dataset were determined by five methods respectively. The threshold for outlier detection was set 0.05 (5%) for all the methods. The outliers detected by at least 3 methods were taken as the putative DE-miRNAs. The proportion of the putative DE-miRNAs in the outlier list found by each method was then calculated. Median value of the percentage among 5 datasets was then calculated as the accuracy for each method.

Determination of the differentially expressed microRNAs and mRNAs

Outlier microRNAs and mRNAs were detected with MOST, implemented in R scripts by Lian et al. [22]. Outliers that rank top 5% were extracted as differential candidates.

Target gene prediction for DE-microRNAs

The mRNAs targeted by the DE-microRNAs were obtained from three target prediction algorithms (TargetScan [26], PicTar [27] and miRanda [28]) as well as a database with experimental evidence (miRecord [29]). In order to obtain a more reliable result, we removed the targets found by only one prediction programs.

Functional enrichment of DE-miRNA targets

We then mapped the target genes of DE-miRNAs to different databases, such as GO (Gene Ontology), KEGG, GeneGo, for functional enrichment analysis. GO and KEGG pathway enrichment were performed using DAVID Bioinformatics Resources 6.7 [30]. GeneGo pathway analysis was performed by MetaCoreTM (GeneGo Inc). In MetaCoreTM, P-values were calculated by hypergeometric distribution to evaluate the statistical significance of the enriched pathways. MetaCoreTM used FDR (False Discovery Rate) adjustment for multiple test correction.

Detection of DE-miRNAs with a novel statistical method

T-statistics is most widely used in differential gene expression detection for microarray studies. However, there is a growing realization that the activation patterns of oncogenes are highly heterogeneous. Some oncogenes show altered expressions only in a minor fraction of samples. Tomlins et al. [31] showed that t-statistics has poor detection power for such unconventional circumstance. The problem with t-statistics has motivated a variety of novel analytical methods such as COPA (Cancer Outlier Profile Analysis) [21], MOST (Maximum Ordered Subset T-statistics) [22], ORT (Outlier Robust T-statistics) [23] and OS (Outlier Sum) [24]. In our previous study [32] we have applied these novel methods to prostate cancer microarray datasets. We have demonstrated that the choice of outlier detection method can greatly affect the genes that are identified.

We compared the performance of these methods in outlier detection based on microRNA expression data in ccRCC. Each method was applied to 5 public microRNA profiling datasets (details of these datasets are given in method section) to obtain the outliers. Outliers detected by at least 3 methods were designated as putative DE-miRNAs. The proportion of the putative outliers in the original gene list was then calculated and the median observation among 5 datasets was defined as the accuracy for the method, as illustrated in Figure 2.

Among the five algorithms for comparison, MOST gives best accuracy. ORT also displays favorable performance. COPA performs a little worse than ORT, but still offers some improvement over the t-statistics. The performance of the OS statistics is noticeably worst. The advantage over t-test is apparent for the majority of the novel methods. Thus we come to the conclusion that the newly developed statistics show superior performance in detecting outliers, and therefore provide promising alternatives to the traditional t-statistics. In this study we chose MOST to identify differentially expressed outliers from public datasets. MOST seems to have superior performance when the number of activated samples is unknown [22]. The threshold of outliers is set as 0.05 to select the top 5% of the miRNA outliers.

Refinement of DE-miRNA lists with the Pipeline of Outlier microRNA Analysis (POMA)

α: Number of outlier genes targeted exclusively by a specific microRNA; β: Number of all the outlier genes targeted by that microRNA; (α, β >1).

POMA improves inter-dataset consistency

POMA predicts a list of microRNA that might play regulatory role in ccRCC. These microRNAs serve as a filter that removes the DE-miRNAs without actual function. After filtration, the DE-miRNA list is greatly reduced, resulting in a robust set of functional microRNAs. In addition, a higher overlap between different miRNA datasets was observed. The inter-dataset overlapping percentage before and after POMA filtration is illustrated in Figure 3.

Furthermore, we performed a literature search of the 11 microRNA markers to validate their relevance in the regulation of ccRCC. It’s found that all of them have been reported for their roles in renal cell carcinoma. The numbers of supporting literatures for each microRNA are also given in Table 2. These “literature curated” microRNAs provide a focused and robust signature to separate normal from cancerous kidney tissues.

Targets prediction and functional enrichment

Identification of ccRCC related functions and pathways

DE-miRNA target genes are statistically enriched in GO terms of cell cycle and transcription. The top enriched GO terms include: regulation of cyclin-dependent protein kinase activity, transcription, DNA-dependent regulation of transcription, regulation of cell proliferation, sequence-specific DNA binding and transcription regulator activity.

KEGG pathways that are significantly enriched with the DE-miRNA targets were also identified, many of which are associated in cancer, e.g. colorectal cancer, Cell cycle, Neurotrophin signaling pathway, Renal cell carcinoma, Prostate cancer, MAPK signaling pathway, and p53 signaling pathway.

The top GeneGo pathway maps regulated by DE-miRNA converge on cell adhesion, cell cycle and cytoskeleton remodeling, most of which are known to play a part in tumor development. Among the significantly enriched GeneGo pathways (FDR < 0.001), 60 were shared by all of the 5 datasets (see Additional file 4). To evaluate the relevance of these pathway maps in ccRCC, we searched PubMed for the published papers describing their constituent network objects in ccRCC. The network objects with previous literature support were considered to be ccRCC-related.

After text mining, 36 out of the 60 pathways (60%) were found to be highly saturated with well-characterized ccRCC-related objects (enrichment ratio>0.15, p-value<0.0001). To visualize the most significantly enriched pathways we constructed a scatter plot (Figure 4) by plotting the -log10 of p-value versus gene enrichment ratio on the y- and x-axes, respectively. The most meaningful points that display both large enrichment ratio (>0.15, x-axis) as well as high statistical significance (P<0.0001, y-axis) were shown in red. These points could be the potential regulatory pathways in renal carcinogenesis. The top 10 significant GeneGo pathways enriched with ccRCC-related objects are listed in Table 4.

A further PubMed search highlighted 22 out of the 36 putative ccRCC-related pathways with literature support in ccRCC pathogenesis. The remaining 14 pathways without previous annotation are considered to be promising novel pathways contributing to ccRCC (See details in Additional file 5). Among the novel GeneGo pathways, TGF, WNT and cytoskeletal remodeling is most significantly enriched, the pathway map is illustrated in Additional file 6. The pathway is highly saturated with network objects previously found to associate with RCC, such as TCF, AKT, VEGF-A, WNT, Frizzled, TGF-beta, RhoA, Beta-catinin, c-Myc, Cyclin D1 and c-Jun. This pathway focuses on WNT protein family and its downstream effectors. The WNT signaling pathway plays a central role during tumorigenesis and inappropriate activation of this pathway has been observed in many human cancers. Wnt ligands first bind to Frizzled family of Wnt receptors to form Wnt-Frizzled complexes. Upon binding with the Axin-related protein, Wnt regulates the stability of catenin β, which is known to play essential roles at cell-cell adherence junctions. Catenin β then binds to TCF, a family of transcription factors, inducing the transcription of Wnt target genes, such as c-Myc, c-Jun and Cyclin D1. Both c-Myc and c-Jun are oncogenic transcription factors that function in cell cycle progression, apoptosis and cellular transformation. Cyclin D1 is a nuclear protein involved in cell cycle progression in G1/S transition. Activation of these genes contributes to early RCC development. These findings are consistent with the observation that WNT signaling pathway is deregulated during renal carcinoma development. Actin and a variety of actin-binding proteins are also central in the pathway. Remodeling of actins regulates the motility of cells and maintains the cytoskeleton. Cytoskeletal actin disruption is the key factor that triggers apoptosis.

TGF, WNT and cytoskeletal remodeling pathway also contains some objects without previous annotation in ccRCC carcinogenesis, such as ROCK, MEK1, p38, MAPK, axin. These objects could be putative therapeutic targets for novel treatment strategies against ccRCC.

Another pathway preferentially targeted by the DE-miRNAs is Brca1 as a transcription regulator, which belongs to the DNA damage category. The pathway map is illustrated in Figure 5 (drawn by MetaCore™, see Additional file 7 for the notation of each sign in this figure). It’s evident in Figure 5 that the pathway is enriched with RCC-related objects such as c-Myc, E2F1, Rb protein, IGF1-receptor, VEGF-A and cyclin D1, most of which are predicted targets of miRNAs differentially expressed in ccRCC. The pathway also includes several targets of DE-miRNAs whose association with ccRCC has not been reported before, e.g. Brca1, sp1, sp3, MSH2, p21, GADD45 alpha and stat1. Central in the pathway is the Breast cancer associated gene-1 (Brac1) which plays a central role in DNA repair. BRCA1 is a well-known tumor suppressor gene in breast cancer and ovarian cancers [36], but its role in other cancer types remains elusive [37]. This is the first study demonstrating the tumor suppressor activity of Brac1 in ccRCC. Brca1 regulates the transcription of proteins at DNA repair pathways via transcription factor p53, such as mismatch repair protein MSH2 [38]. Brca1 also participates in cell cycle regulation. In the absence of DNA damage, Brca1 is associated with ZBRK1 in a complex which inhibits transcription of GADD45 alpha. Upon DNA-damage, Brca1 is phosphorylated and dissociated from the Brca1-ZBRK1 repression complex [39]. The released Brca1 stimulates transcription of GADD45 alpha [40]. GADD45 alpha participates in DNA-damage-induced G1/S checkpoint arrest [41] and DNA-damage-induced G2/M checkpoint arrest [42]. In addition, Brca1 regulates the transcription of some other G1/S checkpoint arrest regulators, e.g. p21 and Cyclin D1. Transcription of p21 [43] may be activated through p53 whereas transcription of Cyclin D1 is regulated via c-Myc [44].

Signatures at pathway level are more consistent

We performed pair-wise comparison between five datasets at different observation levels, including DE-microRNA, target gene, GO-MF (Molecular Function), GO-BP (Biological Process), and GeneGo pathway, respectively. For 5 miRNA expression datasets, 10 pairs are available for comparison. P-values by paired t-test are well below 0.05, indicating that the overlapping percentages at functional level are significantly higher than that at individual DE-microRNA and target gene level.