miRmap: Comprehensive prediction of microRNA target repression strength

Experimental data

Expression microarrays of miRNA-transfected HeLa cells from Grimson et al. (6) were downloaded from GEO (GSE8501) for miRNAs 122 a, 128 a, 132, 133 a, 142, 148 b, 181 a, 7 and 9. We used expression data at 24 h post-transfection and, similar to Grimson et al. (6), only selected probes with signal intensities above the median in the control transfection experiments to retain only the transcripts expressed enough to observe miRNA silencing.

Similarly, we downloaded the expression microarrays from Linsley et al. (24) from GEO (GSE6838). We used the experiments GSM156522, GSM156523, GSM156524, GSM156545, GSM156546, GSM156547, GSM156548, GSM156553, GSM156557, GSM156559, GSM156576, GSM156577, GSM156578, GSM156579 and GSM156581, measured at 24 h with the same experimental conditions. We applied the same selection filter as above (6).

We downloaded the Selbach et al. (7) proteomics over-expression data directly from the web site dedicated to the article. We included expression fold-changes measured at 32 h for miR-1, miR-155 and miR-16 but excluded let-7 b and miR-30 a as these miRNAs exert a negative feedback effect on the RNA silencing pathway (7,21).

HITS-CLIP data from the Chi et al. (9) study were also downloaded from the web site dedicated to the article. After cross-linking Argonaute (Ago) with its neighbouring RNAs, the authors immunopurified Ago and sequenced the pulled-down RNAs. We used the peak height as a measure of miRNA targeting for the 20 available most abundant miRNAs and filtered the relevance of the peaks using a biological complexity (BC, a measure of reproducibility between biological replicates) criterion strictly superior to 1 for medium stringency.

Hendrickson et al. (25) injected miR-124 into HEK293T cells and measured (i) the miR-RISC association with Ago IP, (ii) transcriptome expression with microarrays and (iii) translation activity with polysome fractionation. We used dataset number 5 from the Supplementary Information which includes all measurements for each transcript.

Sequence data

RefSeq 47 (26) mapped on the human (hg19) and mouse (mm9) genomes by the UCSC (27) were used to define mRNA annotations, restricted to ‘NM_’ transcripts. miRBase 16 (28) was used for miRNA annotations.

Target prediction features

Relative importance of features

We computed the relative importance of features in the multiple linear models with the CAR method (36) which decomposes the proportion of the variance explained by each variable of a model while taking the correlations among variables into account.

miRNA target prediction library

We developed a comprehensive prediction model implemented as the miRmap open-source Python library (Figure 1) with a total of 11 features covering a wide range of published and novel methods (Table 2). With our own implementation, we compared the different features without the biases inherent to comparison of pre-computed predictions. We evaluated the features' individual predictive power, measured their intercorrelations and examined different combinations of methods. Additionally, in order to facilitate the library usage, five features are implemented in pure Python.

Novel methods include (i) a more accurate way to compute the binding energy between the miRNA and the mRNA based on the ensemble free energy instead of the minimum free energy, (ii) an exact method to compute the probability that the seed match is an over-represented motif in the 3′-UTR and (iii) a non-empirical statistical test to assess the significance of target site evolutionary conservation.

Correlation among features

We identified potential miRNA target sites by searching for matches to canonical 7-mer seeds on all 3′-UTRs of the human transcripts and predicted their strengths using the 11 methods of our miRmap library (see above and ‘Materials and Methods’ section). We focused our analysis on 7-mer seeds rather than shorter 6-mer seeds as stronger mRNA repression is associated with longer seeds. While this choice results in greater confidence in our feature performance analysis, target prediction with increased sensitivity could be easily obtained by integrating shorter seeds (see below). To evaluate the target site rankings computed with each feature and ignore other differences, e.g. their variances, we computed the Spearman rank correlation between feature pairs (Supplementary Table S1). The absolute values are plotted in Figure 2A.

The three most highly correlated feature pairs are those that measure the same underlying parameters using slightly different approaches: ‘ΔG duplex’ and ‘ΔG binding’ with 0.962, ‘P.over exact’ and ‘P.over binomial’ with 0.806 and ‘ΔG open’ and ‘ΔG total’ with 0.725. ‘ΔG open’ and ‘AU content’ show a correlation of −0.635; as folding algorithms rely on pairing and stacking energies that are stronger for GC than AU pairs, AU-rich sequences form potentially less stable structures, which explains the inverse correlation between ‘ΔG open’ and ‘AU content’. Since these two features evaluate the accessibility of mRNA to miRNA repression, we grouped them in an ‘accessibility group’ together with ‘ΔG total’.

As only the top miRNA target predictions are often used in experimental studies, we measured the overlap among features for their best quartiles. On the first Venn diagram (Figure 2B), we present one feature per group (accessibility, conservation and probabilistic), revealing the low overlap among these methods. The second Venn diagram (Figure 2C) confirms that ‘ΔG open’ and ‘AU content’ features belong to the same accessibility group whereas ‘ΔG duplex’ is a distinct feature not related to the target accessibility. However, target prediction program comparisons (see ‘Introduction’ section) often include PITA (12) which combines both ‘ΔG open’ and ‘ΔG duplex’, making any conclusions made in these comparisons about individual feature performance inaccurate.

Individual feature performance

We evaluated the performance of each feature using data from seven experiments coming from five studies (Table 3) that cover different aspects of miRNA repression and use different assay techniques. (i) Chi et al. (9) performed an Ago-RNA cross-linking experiment followed by IP and sequencing from which miRNA binding sites were assayed. (ii) Hendrickson et al. (25) performed an Ago-IP without cross-linking that we included to underline the effect of the cross-linking step. To measure the effect on mRNA levels, we used studies based on miRNA transfections followed by microarray measurements from (iii) Grimson et al. (6), (iv) Linsley et al. (24) and (v) Hendrickson et al. (25). To assess the effect of miRNA on translation, we took advantage of polysome fractionation experiments from (vi) Hendrickson et al. (25), and of proteomics experiments from (vii) Selbach et al. (7) based on the pSILAC technology to obtain the final translation output.

We identified potential miRNA target sites by searching for matches to canonical 7-mer seeds on the transcripts involved in each experiment and predicted their strength with the 11 methods implemented in our miRmap library, and an additional feature derived from the PhastCons UCSC track (see ‘Materials and Methods’ section) to facilitate comparisons with Wen et al. (23) results. We then evaluated the correlations between the measured and predicted miRNA repression strengths.

We focused our first analysis on the transcriptomics data, as these experiments measure a predominant effect of miRNA repression (38,39) and have the largest scale (‘Trans.Grimson’, ‘Trans.Linsley’ and ‘Trans.Hendrickson’ involve a total of 24 miRNAs). Figure 3 shows the linear regressions and correlations between each feature and the observed reductions in mRNA levels for the ‘Trans.Grimson’ dataset (Supplementary Table S2). The correlation coefficients range from 0.000 for the worst performing feature, ‘ΔG duplex’, to −0.229 for the best feature, ‘AU content’. The next best features are ‘PhyloP’, ‘PhastCons’, ‘ΔG total’, ‘ΔG open’, followed by ‘P.over exact’ and ‘BLS’. Two of our novel features show better correlations than their related features: (i) ‘PhyloP’ is the best performing conservation method (−0.205) and (ii) ‘P.over exact’ performs better than ‘P.over binomial’, i.e. computing the exact probability distribution is better than using the binomial approximation (0.170 versus 0.147). In addition, (iii) considering the ensemble energy outperforms using only the MFE (‘ΔG binding’: 0.023 versus ‘ΔG duplex’: 0.000).

In our second analysis, we examined all the datasets in order to compare the performance of each feature across additional aspects of miRNA repression, assessed through IP, proteomics and polysome fractionation experiments. Correlations for each feature and each experimental dataset are plotted in Figure 4 (Supplementary Table S2). Remarkably, feature performances show high consistency between each of the experimental datasets: accessibility features (red) always perform well, while binding energies (light blue) are always poorly predictive. As target sites in our study contain a seed, the part of the binding energy discriminating target sites is due to the seed nucleotide composition and to the pairing outside the seed. This energy does not drive the miRNA repression strength, as confirmed by the low performance of ‘3′-pairing’. Moreover, the ranking of each feature performance is very similar between datasets using the same experimental techniques, e.g. the ‘Trans.Grimson’ and ‘Trans.Linsley’ datasets. While based on only a single miRNA, the ‘Trans.Hendrickson’ dataset shows better overall performance with only minor differences: ‘UTR position’ improved its ranking while ‘PhastCons’ is outperformed by ‘BLS’.

‘AU content’ consistently provides the best measure of target site accessibility. This is in agreement with findings from Wen et al. (23), but in contrast to results from Hausser et al. (21), which described better performance with ‘ΔG open’ for an IP experiment. However, for the ‘IP.Hendrickson’ dataset, which, like Hausser et al. (21) involved IP without cross-linking, ‘AU content’ and ‘ΔG open’ perform equally well. The ‘IP.Hendrickson’ experiment is also distinguished by the probabilistic (purple) and ‘UTR position’ (green) features that outperform the conservation features (grey), which may be explained by the lower precision of this method (i.e. IP without cross-linking), performed with a single miRNA.

The best conservation feature performance is generally slightly lower than the best accessibility feature, but it outperforms ‘AU content’ for the proteomics and HITS-CLIP datasets. ‘PhastCons’ performance on the HITS-CLIP dataset is consistent with findings from Wen et al. (23). Our novel conservation feature, ‘PhyloP’, shows the best or tied-best performance for five out of the seven datasets. When outperformed, it is only marginally outperformed implying that ‘PhyloP’ is the best overall conservation feature.

Hendrickson et al. (25) polysome fractionation measured the miRNA effects as ribosome occupancy (fraction of a given gene’s transcripts associated with ribosomes) and ribosome density (the average number of ribosomes bound per unit length of coding sequence). Effects caused by the miRNA on both parameters were detected by the authors, but were substantially higher on the ribosome density, in agreement with the absence of correlation with the ribosome occupancy that we observed, i.e. this measurement is not quantitative. However, the ribosome density is a quantitative measure of the miRNA effect, as the correlations were as high or higher than those of the large-scale transcriptomics experiments. We observed again, as for all Hendrickson et al. (25) datasets, a higher correlation for the ‘UTR position’ feature, probably caused by the experimental setup.

Combining prediction features

The features correlate linearly with experimentally measured miRNA repression levels. We combined 10 features of our miRmap library (we excluded ‘ΔG total’, as this feature is simply the sum of ‘ΔG duplex’ and ‘ΔG open’) with a multiple linear regression on the ‘Trans.Grimson’ dataset (P = 4.9 × 10⁻¹¹⁰; Supplementary Figure S7). This model explains 12.7% of the variance, close to a 2-fold increase over TargetScan context score (6): with the same type of regression, the three features of TargetScan context score (‘AU content’, ‘3′-pairing’ and ‘UTR position’) explain only 7.49% of the variance. This improved performance of our model is confirmed by the higher correlations with the experimental measurements, computed in the same manner as the individual feature correlations (Figure 5A). The contribution of each feature (on ‘Trans.Grimson’ dataset, Figure 5B) generally mirrors the rankings based on individual feature correlations: ‘AU content’ is the most explanatory feature, but ‘P.over exact’ contributes more in the regression model than its correlation rank suggests. Interestingly, the conservation features ‘PhyloP’ and ‘BLS’ contribute ∼14 and ∼11%, respectively, despite using the same input data (multiple genome sequence alignments) both contribute substantially to the explanation of the variance. Among the accessibility features, ‘ΔG open’ contributes only half as much as ‘AU content’ (15 and 30%, respectively). A model limited to the five features with the greatest contributions in the model with all the features (they represent 90.5% of the variance explanation of the full model) still explains 11.6% of the variance.

Instead of evaluating the model directly in terms of the explained variance, the quality of the ranking can be estimated by ordering the target sites by predicted strength, binning them and computing the mean expression fold-change of each bin. This approach, also used in (40) to evaluate the ranking of different tools for predicting miRNA repression strength on translation with proteomics data, was applied to 10 quantiles of the ordered predictions (Supplementary Figure S2). The overall distribution was shifted to lower fold-changes for miRmap compared with TargetScan context score, indicating a better ranking as a decrease in fold-change corresponds to greater repression. For the first quantile, the mean fold-change was reduced from −0.32 to −0.39 with miRmap.

Multiple linear regressions with the other datasets further support the conclusions from the analyses of individual feature performance (Supplementary Figures S1 and S3). They confirmed (i) the importance of ‘PhyloP’ for the ‘IPcross.Chi’ dataset (64% of R²) over 24% for ‘AU content’, (ii) the similar importance of ‘PhyloP’ and ‘AU content’ for proteomics (31% and 39% of R², respectively) and (iii) the relevance of polysome fractionation experiment (‘RibN.Hendrickson’ dataset) to measure miRNA repression strength compared with proteomics as 10.6% of the variance was explained by the model (5.75% for proteomics). We also observed that the model computed on the ‘Trans.Linsley’ dataset explains only 4.36% of the variance even though this dataset is larger and based on the same techniques as the ‘Trans.Grimson’ dataset (R² = 12.7%).

Shorter seeds may also promote miRNA repression, but usually with lower efficiencies (4). We therefore tested our approach on canonical 6-mer seeds by computing a model with these seed matches on the ‘Trans.Grimson’ dataset. While the global importance of each feature remained generally similar, with accessibility features being the most explanatory, R² dropped to 8.31% of the variance (Supplementary Figure S4A), which still outperforms TargetScan context score (R² = 4.70%). Interestingly, the importance of the ‘P.over exact’ probabilistic feature was reduced from 22 to 7%—falling from second position to fifth—as expected with shorter seeds where matches occur more frequently by chance and are therefore less statistically distinguishable from the background. We also evaluated the model by computing the distribution of fold-changes (Supplementary Figure S4B). As expected, the mean fold-changes were not as low as with the 7-mer seeds, nevertheless they confirmed the better ranking achieved with miRmap compared with TargetScan context score, e.g. the mean fold-change of the first quantile was reduced from −0.16 to −0.21. These results were further supported by the analysis of the other datasets (Supplementary Figures S5 and S6).

Combining multiple target sites

Each mRNA can contain many miRNA target sites. Although most experimental datasets focus on a single miRNA at a time (or all miRNAs for the ‘IPcross.Chi’ dataset), a framework that can capture the multiplicity of these interactions should improve the predictive power. We examined three simple functions to combine the individual scores of target sites into a global metric at the mRNA level: the best (minimum or maximum depending on the sign of the correlation), the sum and the log of the sum of the exponentials. For this analysis, we selected transcripts from the ‘Trans.Grimson’ dataset with exactly two target sites, resulting in a sample size of 370 mRNAs (only 53 mRNAs have exactly three target sites). For this study, only features predicting different strengths for each target site in a 3′-UTR are appropriate as they would show different correlations for each function, thereby allowing function comparison. As the probabilistic features compute the probability of a fixed number of seed matches in the 3′-UTR, and as the ‘BLS’ score is also computed for the entire 3′-UTR, they could not be used.

The log of the sum of the exponentials function is designed to approximate interaction kinetics on the principle that stronger sites would drive the observed repression at the mRNA level. However, this function performed poorly for every feature as opposed to the sum (Supplementary Figure S8), which means that every target site has the same importance, indicating that the quantity of miRNA molecules is not limiting the repression reaction in this experiment. Regarding the binding energy features, ‘ΔG duplex’ and ‘ΔG binding’, the minimum energies provided the best predictors, i.e. the best site drives the repression for these two features. In contrast to their relatively poor performance with single site predictions, their performance was substantially increased (with correlations from 0 to 0.094 (P = 0.072) and 0.023 to 0.119 (P = 0.022) for ‘ΔG duplex’ and ‘ΔG binding’, respectively) but they still did not outperform the other features. The performance ranking among the remaining features was not substantially different to the single site predictions and, as already observed before (7), summing was the best option for the majority of them.