Toward a combinatorial nature of microRNA regulation in human cells

Prediction databases

Database files that are included in miRror encompass most of the available miRNA-target prediction tools covering human miRNAs (and several other species). The resources used are collectively called MDBs: (i) TargetScan database (39); (ii) microCosm that is based on the miRanda algorithm (40); (iii) PicTar, according to species conservation scheme (35); (iv) DIANA–MicroT (41); (v) PITA (42); (vi) MirZ (43); (vii) microRNA.org that allows analysis on multiple miRNA acting on the same gene-target, based on miRanda (44); (viii) miRDB resource (45); (ix) TargetRank (46); (x) RNA22 (47); (xi) MAMI (48); and (xii) miRNAMap2 (49). We used the latest versions available (February 2012). miRNAMap2 is based on RNAhybrid, miRanda and TargetScan, while requiring a minimal score from miRanda. Reducing false positives is achieved by increasing the numbers of binding sites and their accessibility. We included a resource that applied statistical machine learning algorithms (i.e. miRDB, (22)) and a meta-server that integrates some of the major resources (i.e. MAMI (48)). The number of candidate genes and the number of miRNAs are indicated for each of the major MDBs (Supplementary Figure S1A, Supplementary Table S1). The union of the miRNAs and genes for all the reported MDBs is ∼2500 and ∼47 000, respectively. These numbers exceed the number of known miRNAs due to inconsistency in naming by the different MDBs. For this study, MDB refers to the entire collection of miRNA—gene matches that result from the underlying prediction algorithm.

Gene expression data

The complete list of experiments was extracted from the GEO (50) and ArrayExpress (51). The platforms that we considered are from Affymetrix (HG-U133 Plus 2.0, 22K probes, 14.5K genes) and Agilent spotted oligonucleotides (Human 3.0 A1, 23K genes). Data were collected from the SOFT files. See Supplementary Figure S1 and Supplementary Table S2 for a complete list.

Prediction performance

The MDBs are not designed to deal with hundreds and thousands of genes as input. In order to assess the performance of each MDB, we score each miRNA by counting the number of down-regulated genes that are predicted as targeted by it. The miRNA that is ranked first is the miRNA that is associated with the highest number of predictions from the down-regulated genes that were tested.

CLIP data were retrieved from StarBase v2.1 (52). The identification of miRNAs from the CLIP data is based on predictions from TargetScan, PicTar, RNA22, PITA and miRanda. We confirmed that 98% of the unified list of miRNAs derived from these 5 MDBs. We exported and analyzed all documented miRNA–target interactions that concerned human cells. For the CLIP data, the MDBs’ prediction accuracy was estimated by counting the number of cases the appropriate gene was predicted from the large number of its regulating miRNAs.

miRror application for miRNA cooperative regulation

The miRror platform (53) was used to match between the input gene list and the miRNAs that best explain such input. The statistical basis of miRror is implemented in the miRTegrate algorithm (53) that calculates the probability of matching a gene list while considering the coverage of the miRNA-targets for each of the prediction databases. Formally, the probability of the miRNA interaction with the input gene set as opposed to the rest of the genes in that MDB is calculated. We used a statistical threshold to ensure that the contribution of any MDB that covers only a small number of miRNAs with high specificity remains significant. On the other extreme, the impact of an MDB that over predicts (i.e. each miRNA matches almost any gene) to miRror results should be minimal. Calculating the P-value for genes as input was performed according to the hypergeometric distribution.

Where N is the total number of genes in the database, n is the number of genes predicted for the specific miRNA, m is the number of genes in the input set and k is the number of genes that appear in both the input set and the predicted miRNA targets, the probability of such k is:

P-value is obtained by summing over all X ≥ k (53). A correction for multiple testing was included using the false discovery rate (FDR) procedure. The P-value threshold that was used in this study is 0.05 (unless mentioned otherwise).

In this study, we used the miRror2.0 that includes number of new capabilities. Most notably, implementing an optional analysis on the basis of miRNA families, incorporating miRIS (miRror Internal Score) as a quality score, adding several MDBs (RNA22, MAMI and miRNAMap2) and updating the new scoring method presented by miRanda-microRNA.org.

Mapping genes to pairs of miRNAs

For each prediction results, we considered the 5 top predictions and analyzed by miRror the success of each pair (10 pairs). The ranking of miRror is according to miRIS. This score balances between (i) the proportions of predicting MDBs out of all tested MDBs and (ii) the fraction of the potentially regulated genes from the entire input genes. The 10 pairs are tested by their ability to increase the miRIS score beyond the miRIS obtained from the individual miRNA. Each pair is tested for the degree of complementing each of the individual miRNA and the degree of overlapping, redundant mode. These values range from 0 to 1. Formally:

An increase in miRIS to the level that was considered a standard deviation apart from the calculated value is considered a successful recovery of miRNA.

An inherent redundancy in miRNAs is evident in miRNA families (as defined by miRBase). The redundancy is based on target binding site and seed overlap. Still, several MDBs do not cover the same miRNA or the same targets, thus small differences are still found. These cases are often identified by high redundancy and low complementation. The consideration of pairs of miRNAs provides the notion of joint miRIS (JM), an intuitive generalization of miRIS:

The best five ranked candidates are chosen. We tested the joint miRIS values and the Jaccard Index (JI) to identify candidates of miRNA-Duos that explain best the gene expression profiles. Analyzing all combinations of the top 5 miRNAs provides 10 pairs per for each experiment.

miRNA family unification

A source for inconsistency and thus a reduced quantity of prediction is attributed to the use of slightly different names for the same miRNAs. A compilation of miRNAs to their families was applied according to the miRBase notations. Compressing the miRNAs to their families reduced the number of entries by 2.7-folds (Supplementary Table S3). Another conservative and systematic approach assigned miRNAs to their seed families. We cluster miRNAs into the same seed family if nucleotides 2–7 of the mature miRNA are identical. Seed families compress the number of miRNAs by 1.6-folds (on average).

Limited success in identifying the underlying miRNA from over-expression experiments

With the goal of revealing global trends in miRNA regulation, we considered the major miRNA resources. Data from MDBs that comprise the most stable resources were collected. We designed a systematic analysis of each MDB from an unbiased list of transcripts. Several criteria were used in order to select experimental data: (i) Experiments are conducted on human cell-lines that are transfected with a single miRNA. (ii) The entire cell transcriptome is compared with controlled cells (often introduced with a scrambled or mutated sequence). (iii) The platforms used for the gene expression profiling cover most known genes (See Materials and Methods). We tested the potential of recovering the appropriate miRNA by each MDB from the unbiased transcriptomic data. We selected 12 of the leading resources (PITA-Top, DIANA-microT3.0, MiRanda-microRNA.org, MicroCosm, TargetScan, TargetRank, PicTar-4ways, MirZ-ElMMo, RNA22, MAMI, miRDB and MAP2) for the rest of this study (Supplementary Figure S1).

The transcriptome from hsa-miR-124 overexpression experiment (54) (GSE6207, HepG2 cells, time series for 4–120 h after transfection) was used as a test case. We found that the success in predicting the subjected miRNA is restricted to data collected from ≥24 h after transfection. A maximal level of expression for the functional miRNA is reached between 24 and 48 h after transfection. Data from earlier time points (4–16 h after transfection) were uninformative and have not been further analyzed. We focused on the down-regulated genes at a permissive repression level of 1.2-fold relative to control cells (Supplementary Table S4). Different names for miR-124 (i.e. miR-124a and miR-124-3p) were unified to limit the inconsistency for these results. The 100 most significant down-regulated genes were used as the input for the 12 MDBs (see Materials and Methods). We found that 75% of the MDBs identified hsa-miR-124 as their top prediction and the rest of MDBs predicted hsa-miR-124 among the top 9 predictions. Notably, the actual degree (e.g. fold change) by which miRNA reduced the amount of its candidate targets is a poor indicator for a successful recovery (t-test, P-value = 0.8).

We collected raw transcriptomic data from additional 30 large-scale experiments that cover 26 different miRNAs from a variety of cell-lines (Supplementary Figure S1). Figure 1 illustrates the workflow for the miRNA-gene prediction task. All together we produced over 10 000 prediction results (30 experiments, 12 MDBs, 8 input sets, 4 selections of MDB internal scores). A common protocol for the 12 MDBs (Figure 1, right) and miRror (53) as a unified-predictor (Figure 1, Left) is shown. The miRror platform examines the combinatorial nature of miRNA regulation. miRror seeks the miRNAs that are the most likely explanation for the observed gene expression profile. For this study we applied miRror2.0, an improved version of the original platform (see Materials and Methods).

Figure 2 shows the performance for each MDB in view of the miRror2.0 results. We marked the success in predicting the appropriate miRNA at various levels of accuracy (ranked 1, 2 and 3–9). Figure 2A shows the results from 30 experiments and across 12 MDBs. Most MDBs fail to recover over 50% of the subjected miRNAs. miRror (Figure 1A, purple bar) outperformed all MDBs with 60% success (for the 1–9 ranked predictions, Supplementary Table S5).

Internal prediction scores fail to improve the quality of predictions

Most MDBs provide internal scores for ranking their prediction results. Such scoring is presented in order to reduce the inherent numbers of false predictions (discussed in (55)). We tested whether the miRNA predictions can be improved by applying the internal scoring as provided by each MDB (Figure 2B). The internal scores for the MDBs differ by their underlying principles, their distribution and range. We thus used percentiles (i.e. top 10, 25 and 50%) to select the suggested top predictions as defined by each MDB. Figure 2B shows prediction results, using the top 25% according to the internal scores (only 10 of the 12 MDBs support internal scorings). Surprisingly, only a few of the MDBs (microCosm, TargetScan and Miranda-microRNA.org) improved their success when the filter of the internal score was activated. For the other MDBs, the prediction success was reduced. Consequently, a slight reduction in the number of successful prediction is recorded for miRror (Figure 2B, dark purple, Supplementary Table S5).

An experiment-centric view shows that only some experiments are associated with a successful recovery of the relevant miRNA (Figure 2C, green line). The average number of miRNA predictions for the analyzed MDB is 720 (Supplementary Figure S1 and Supplementary Table S1). For about half of the experiments the average success by the MDBs is very poor (Figure 2C, orange and red lines).

The high performance of miRror in view of each MDB is intriguing (Figure 2). We thus assessed the consistency of the results by combining a large number of prediction tests (Supplementary Figure S2). The results confirmed the robustness of miRror predictions in view of the individual MDBs. The miRror2.0 protocol consistently outperformed the 12 other MDBs with 55–60% success. The performance of miRror is robust in view of thousands of predictions across a wide range of parameters (Supplementary Table S6).

miRror score unifies the combinatorial nature of the input and the coherence of predictions

Based on the reported results, we developed a scoring method that benefits from the individual MDBs. The score, called miRIS, balances between (i) the proportions of predicting MDBs out of all tested MDBs and (ii) the fraction of the potentially regulated genes from the entire input genes. A default of two successful MDBs is a minimal demand for consistency, as well as at least two genes from the input set—a minimal demand for cooperativity.

Figure 3A shows the contribution of the two components that comprise miRIS. The analysis is based on hsa-miR-124 (GSE6207, HepG2 cells, 24 h after transfection). We analyzed 240 genes as input (10% down-regulated genes, P-value = 0.05, Supplementary Table S4). Figure 3A shows the contribution of the fraction of MDBs and the fraction of hits from the input list. The two measures are not replaceable and only weakly correlate (r = 0.3–0.5). miRIS balances between these two complementary components (Figure 3A).

An additional statistical prerequisite is the P-value for the miRTegrate algorithm (P-value ≤ 0.05, see Materials and Methods). This threshold remains a useful parameter for tuning the stringency by changing the statistical significance of the analyses. The application of a stringent P-value tends to increase the prediction success (not shown).

Presenting the entire list of miRror predictions revealed that remarkably, only few miRNAs are associated with the highest miRIS values (Figure 3B and C). We illustrate the miRIS ranking for hsa-miR-7 (Figure 3B) and hsa-miR-124 (Figure 3C). Notably, all other 12 MDBs failed to recover the hsa-miR-7 (Supplementary Tables S5). For 75% of the successful predictions miRIS values drops after 1–3 predictions. This trend is not dependent on the size of the input or the length of the predicted list (compare Figure 3B with Figure 3C). Although the absolute miRIS values might be relatively low, an abrupt drop in miRIS often signifies a successful prediction.

Combinatorial nature of miRNA regulation is reflected by CLIP-Seq data

We further tested the combinatorial nature of miRNAs in human cells. The discussed miRNA experiments (Supplementary Figure S1) represent non-physiological settings that are prone to off-target noise (56). Therefore, we tested the combinatorial concept on CLIP technology experiments (24,27). Specifically, validated data for ∼6500 genes that are regulated by 653 miRNAs are compiled in StarBase (52). Note that based on the CLIP technology, the actual sequence from the mRNA that is bound to miRNA is identified by sequencing. Thus, it eliminates the largest source of error (e.g. the label of a site as true of false).

Figure 4A shows a cumulative representation of miRNA regulation according to the number of miRNAs per gene. Only 10% of the genes are regulated by a single miRNA. We noted that 50% of the targeted genes are associated with up to 8 miRNAs. About 90% of the genes are targeted by up to 35 miRNAs, with a few genes that are regulated by as many as >80 miRNAs (Supplementary Figure S3). The analysis includes the task of recovering the human genes by an input of the set of miRNAs that jointly regulate it (Figure 4B). We use the miRror procedure in a miR2Gene mode (Figure 4B). Figure 4C shows the analysis for the 12 MDBs and miRror (ranked by miRIS) for the genes that are regulated by 30–40 miRNAs (362 genes, Supplementary Table S7). The results show that for 33% of the genes, miRror correctly determined the subjected gene whereas the other MDBs failed to reach such success (success is considered a correct prediction among the top ranked 20 results). miRror also outperforms the other methods for recovering the genes at the top 1% of the predictions (Figure 4C, inset). Note that half of the MDBs failed to recover 90% of the relevant genes (Figure 4C).

Evidence for miRNA pairs governing cellular regulation

The tools and methodologies developed for the experimental data from miRNA overexpression (Figures 2 and 3) and CLIP datasets (Figure 4) are used to explore the concept of ‘miRNA-Duos’. A Duo (pair) of miRNAs is the simplest form of coordinated co-regulation by miRNAs and thus a starting point for the assessment. The number of human genes that are reported by the numerous CLIP data is ∼6500. A vast majority of them (90%) are regulated by multiple miRNAs (Figure 4A).

For a quantitative view of all possible pairs of miRNAs, we used the human CLIP data as an in vivo evidence for miRNA-Duos. There are 212 878 human miRNA-pairs in StarBase. Regulation by a miRNA-Duo can be generally described as one of two modes: (i) a miRNA-Duo that expands the set of gene-targets, thus allowing a higher coverage of the genes; and (ii) a miRNA-Duo with overlapping targets that can be considered as a backup mode (Figure 5A). The two extreme scenarios are formulated using the JI. Intuitively, JI is a straightforward measure for comparing the correspondence (intersection) and expansion (union) of two sets. Thus, a low JI value is indicative of the expansion mode whereas high JI indicates dominance by a backup mode (Figure 5A). We noted that the CLIP data have exceptionally low JI (average = 0.020, Supplementary Table S8).

The miRNA pairs (from the CLIP data) are tested by JI (Figure 5A). Each MDB is centered at a characteristic range of JI. This is strongly affected by the number of shared targets between the pair (Figure 5B). The average JI is 0.032, 0.027 and 0.096 for PITA-Top, TargetScan-conserved and Miranda-microRNA.org, respectively (Supplementary Table S8). The tendency of each MDB toward a restricted range of JI exposes the inherent bias of each MDB.

Inspecting transcriptome profiles for evidence of a combinatorial miRNA regulation mode

We tested the validity of miRNA-Duo regulation in a wide range of cellular contexts. As miRror provides a robust approach for a combinatorial regulation by miRNA, we tested miRror’s capacity toward miRNA-Duos. Figure 5C illustrates a simplified test case. Each miRNA is assigned with its calculated miRIS. An edge connecting two nodes is marked by JM value which is associated with a miRNA pair. Specifically, the JM considers the miRIS of each of the miRNAs under the assumption that the miRNA pair is a new entity. For example, a miRNA that covers 50% of the genes query and 50% of the MDBs will have a miRIS of 0.5. Similarly, the miRIS of its pair miRNA is also 0.5. For this example, the JM must range from 0.5 to 1.0. The JM is thus an indirect measure of the intersection and union (measured by JI) of the miRNA in the pair. We illustrate miRNA-Duos in which the overlapping, backup mode dominates (miR-a and miR-d; miR-b and miR-c). Other miRNAs adopt an expanding, complementation mode (miR-a and miR-c; miR-c and miR-d).

As miRNAs that belong to the same family regulate a similar list of genes, a bias toward a backup mode of regulation is expected (Figure 5A). Thus, we integrated a compression of miRNAs into their families (Supplementary Table S3) for eliminating miRNA-Duos that trivially overlap their gene lists. We selected 20 large-scale gene expression experiments (Supplementary Table S9), and set to evaluate the potential of miRNA-pairs to explain the observed transcriptomic profile. We focused on experiments that were conducted on human cell-lines that were exposed to a variety of stress conditions, molecular and pharmacological manipulations. Importantly, in none of these experiments the involvement of miRNAs regulation was proposed.

Figure 5D shows the relationship of two essential features for miRNA-Duos regulation: the joint miRIS (JM, Materials and Methods) and the JI. For each of the 20 experiments, the 5 top ranked predictions (by miRIS) were considered as miRNA-pairs. A correlation between the Joint miRIS and the JI is substantial (R² = 0.245). We focused on pairs that were specified by unusually high JI and low Joint miRIS. Such pairs are miR-25 and miR-32; miR-363 and miR-367. We found that these miRNAs share the same 6-mer seed and thus represent cases of overlapping, backup regulatory mode. Note that the family relations were not reported by miRBase or any other resource. The same is valid for the miRNA-Duo of miR-26 and miR-1297 whose overlap in the seed sequence is overlooked (Figure 5D).

The regulatory strategy of expanding, complementation mode dominates the CLIP data (Figure 5B, CLIP). We search pairs that are characterized by a low JI and a significant Joint miRIS. Such instances can be explained by cooperative activity by a small set of miRNAs. miRNAs that interact with hsa-miR-218 (Figure 5D, orange symbols) are promising candidates. The miRNA-Duos that pair with hsa-miR-218 maximized the likelihood of the observed profile of down-regulation genes. These pairs were predicted for the transcriptome of 293T cells treated with camptothecin (GSE2451). We found no evidence to support miRNA-Duos as a hidden layer of regulation for most (85%) of the tested experiments, as expected (Supplementary Table S10). However, for 3 experiments a significant coverage of the input genes can be explained by miRNA-pair acting together: (i) the RPTEC cells that were exposed to hypoxia (GSE12792); (ii) ultrasound treatment in human leukemia U937 cells (GSE10212); and (iii) 293T cells treated with camptothecin (GSE2451). Camptothecin is a natural cytotoxic quinoline alkaloid drug that inhibits DNA topoisomerase I and is used in cancer chemotherapy. Many of the down-regulated genes from this experiment are associated with cell-cycle regulation genes that are subjected to miRNA regulation. We anticipate that irrespectively to the primary target of the drug, the camptothecin treated cell reached homeostasis through the involvement of miRNA regulation. Indeed, several recent studies confirmed the participation of miRNAs in camptothecin treated cells (57).

Coping with noisy predictions by miRror

An estimation of prediction accuracy by established MDBs is reported in (21). The miRror2.0 system outperformed other MDBs, with the success attributable to: (i) the incorporation of statistical considerations for predictions above a predetermined threshold; and (ii) the concept that is applied to genes or miRNAs sets as input. We demonstrate that our protocol (Figure 1) successfully overcomes the noisy predictions of existing MDBs. Surprisingly, we found that the use of the internal scoring methods had a negative impact on miRror’s success rate (Figure 2B), dropping from 60% to 53% (success being defined as a correct entry in the top 9 predictions). Only Miranda-microRNA.org and microT3.0 improved their success (from 40% to 47% and from 40% to 50%, respectively). Both MDBs have recently refined their internal scoring methods (55,63). MAMI and TargetScan maintained their success rate. The success of the other MDBs had actually deteriorated. For PicTar-4 the success was reduced from 40% to 10% using its suggested internal score. Thus, it seems that the value of the scoring methods for some of the MDBs should be revisited.

A known pitfall in the field of miRNA prediction concerns the inconsistency in miRNA notation schemes. The use of 12 difference resources led to an artificial expansion of the miRNA input set (e.g. hsa-miR-19 is sometimes denoted as hsa-mir-19a, hsa-mir-19b or hsa-mir-19b-2). Furthermore, the hsa-mir-19a* and hsa-mir-19b* are disputably grouped as part of the has-miR-19 equivalence class. Due to the inconsistency in miRNA nomenclature, the union of miRNAs and genes submitted to miRror reached an unrealistic number of 2528 and 46 921, respectively (Supplementary Figure S1). Therefore, we applied the miRNA family compression scheme. The family representation according to miRBase definition compresses the MDBs by 2–4-folds. However, a family assignment that is based on seed sequence identity reduces the number of miRNAs by only 38% (Supplementary Table S3). Such representation improved our prediction success, with only a minimal reduction in sensitivity. Furthermore, the prevalence of the overlapping mode of miRNA regulation was reduced to a minimum following the compressing of miRNAs to their respective families.

Interestingly, many variants of miRNAs (called isomiRs) were found by deep-sequencing studies. Experimental analysis showed that isomiRs and their canonical counterparts act coordinately to target functionally related genes and pathways (64). Thus, in addition to the seed dependent family assignment (Supplementary Figure S3) isomiRs function as an additional mode for increasing the specificity of miRNA regulation.

Exposing a combinatorial mode of action from raw data

miRNAs constitute an additional layer of post-transcriptional regulation. Several studies aim to unveil the design principle underlying miRNA regulation (14). Systems biology has identified small miRNA modules that by simultaneous regulation alter a functional module (e.g. EGFR-driven cell-cycle pathway (65)). The hierarchy among miRNAs suggests that some miRNAs act as pathway hubs whereas others are interconnected with transcription factors to provide small functional modules (60). In this study, we assessed the capacity and limitations of the miRror platform in exposing miRNA cell regulation in healthy human cells. We applied the Joint miRIS (JM) as a measure for the validity of regulation by miRNA-Duos. The assumption being that a small number of coordinated miRNAs act to maintain robust homeostasis. Under this assumption, even a modest modulation executed in a combinatorial fashion can yield to a substantial change in cell phenotype. In our analysis we had not considered the spacing of the binding sites along the transcripts. It was shown that a combinatorial regulation at the level of one transcript is not necessarily cooperative but may be competing. Thus, the outcome of a combined action of miRNAs may be a lower efficacy than a single site (66).

An intriguing observation from the analyzed CLIP data is the dominant appearance of expanding, complementation mode. We confirm an identical distribution for simulated data (randomizing the pairs of miRNA and genes in a two-sided graph) and the observed CLIP data. Thus, we reject the possibility of a biased sampling of the in vivo situation (see Discussion in (24)). Instead, a strong preference in the CLIP data toward the expanded mode suggests that accessibility and competition with RNA binding proteins drives the preferential binding of miRNA on their targets (14). However, a test for each of the MDBs (e.g. TargetScan, Supplementary Figure S3) and the union of the MDBs differ from the properties of the CLIP data (24,27,67,68) (Supplementary Figure S3). Whether the expanding mode is the preferred mode of action in vivo under stress and extreme conditions in yet to be determined.

We anticipate that the miRror platform will be used successfully for screening for overlooked miRNA regulation (in either Gene2miR or miR2Gene modes). In this regards we intentionally excluded from our analysis any transcriptomic data from cancer tissues and transformed cells as many of these datasets are dominated overwhelmingly by cancer-related miRNAs (69). Notably, the majority of the transcriptomic studies (20 experiments, Supplementary Table S10) do not claim any miRNA regulation process. However, for 15% of cell manipulated experiments (3/20), a regulation by miRNAs remains a plausible hypothesis (Figure 5D).

The expression levels of miRNAs in cell-lines (70) and normal tissues are accumulating (2,71). In this view, the expression of hsa-miR-218 and hsa-miR-374 (Figure 5D) co-occur in most independent experiments. Actually, the miRNA-Duo of miR-218 and miR-374 doubled the number of successful hits from 27 to 54 genes (Supplementary Table S10). We expect to identify sets of miRNAs that have the potential to govern the observed transcriptomic profile in view of the co-occurrence in miRNAs expression.