MicroRNAs (miRNAs) are post-transcriptional regulators that bind to their target mRNAs through base complementarity. Predicting miRNA targets is a challenging task and various studies showed that existing algorithms suffer from high number of false predictions and low to moderate overlap in their predictions. Until recently, very few algorithms considered the dynamic nature of the interactions, including the effect of less specific interactions, the miRNA expression level, and the effect of combinatorial miRNA binding. Addressing these issues can result in a more accurate miRNA: mRNA modeling with many applications, including efficient miRNA-related SNP evaluation. We present a novel thermodynamic model based on the Fermi-Dirac equation that incorporates miRNA expression in the prediction of target occupancy and we show that it improves the performance of two popular single miRNA target finders. Modeling combinatorial miRNA targeting is a natural extension of this model. Two other algorithms show improved prediction efficiency when combinatorial binding models were considered. ComiR (Combinatorial miRNA targeting), a novel algorithm we developed, incorporates the improved predictions of the four target finders into a single probabilistic score using ensemble learning. Combining target scores of multiple miRNAs using ComiR improves predictions over the naïve method for target combination. ComiR scoring scheme can be used for identification of SNPs affecting miRNA binding. As proof of principle, ComiR identified rs17737058 as disruptive to the miR-488-5p: NCOA1 interaction, which we confirmed in vitro. We also found rs17737058 to be significantly associated with decreased bone mineral density (BMD) in two independent cohorts indicating that the miR-488-5p/NCOA1 regulatory axis is likely critical in maintaining BMD in women. With increasing availability of comprehensive high-throughput datasets from patients ComiR is expected to become an essential tool for miRNA-related studies.
MicroRNAs (miRNAs) belong to a class of short (18–25 nucleotide) non-coding RNAs that regulate gene expression post-transcriptionally. Their regulatory activity depends heavily on the recognition of target sites located primarily on the 3′-untranslated regions (3′UTRs) of messenger RNAs [1] but also on ORFs or 5′-UTRs [2]. In general, a gene contains multiple miRNA binding sites. Computational miRNA target prediction depends on algorithms that typically use features like Watson–Crick base pair matching [3]–[5], thermostability of binding sites [3], [6]–[12], accessibility of target sites [3], [13], and phylogenetic conservation [3], [5]. Still, target prediction algorithms suffer from high number of false predictions and poor overlap in their predictions [14], [15]. It is worth noting that most existing algorithms only utilize site-specific features [16] ignoring factors like the relative expression of miRNAs that affects binding specificity and target combinatorial effects. The level of miRNA expression affects which targets will be occupied. When a miRNA is expressed at low levels it is expected to bind to only few, high affinity targets. As miRNA expression is increased, and all high affinity targets are occupied, the remaining miRNA molecules can bind to suboptimal targets of moderate affinity (due to target exclusivity). We believe that the 1∶1 stoichiometric binding model used by existing algorithms may be insufficient for the dynamic nature of real miRNA: mRNA interactions. This may be a major drawback of the current target prediction algorithms. Another drawback of most algorithms is that they only consider single miRNA: mRNA pairings, ignoring the fact that multiple miRNAs, each with a moderate effect, can collectively alter significantly the expression levels of a given mRNA. Extending single pairing predictions to the union or intersection of targets of multiple miRNAs does not solve this problem since these models do not consider multiple moderate binding effects and they also tend to be either very conservative or having large numbers of false positive predictions, respectively. Until recently, the only notable exceptions were PicTar [10] and GenMir++ [17]. PicTar uses a hidden Markov model to determine targets of multiple miRNAs, but it does not consider the miRNA expression in determining the relative binding. GenMir++ implements an elaborate Bayesian framework to model miRNA: mRNA dependencies using expression data and prior targeting information. The prior information is provided in binary form (putative targets/non-targets) and the effect of multiple miRNAs on a given mRNA is inferred using an expectation maximization (EM) algorithm. TaLasso [18], a recently published algorithm, combines miRNA expression with prior database information to infer regulation of mRNAs from a set of miRNAs using a lasso regression method to determine the activities of various miRNAs. Finally, Jayaswal et al. [19] proposed another method for finding multiple (many-to-many) miRNA: mRNA interactions by following a two-step approach. First, they identify miRNA and mRNA expression clusters and in the second step they find associations between them. In this case, there is no quantitative modeling of particular interactions. These algorithms have in common that they use miRNA expression to drive the target prediction for sets of miRNAs, but none of them has an underlying suitable thermodynamic model to account for target exclusivity. Furthermore, none of these algorithms has been used before to quantitatively rank single nucleotide polymorphisms (SNPs) affecting miRNA targeting. In general, so far, the evaluation of miRNA-related SNPs in disease has been done with more straightforward approaches [20], [21]. Until recently, a major obstacle for developing comprehensive miRNA binding models was the small number of available experimentally validated target pairs. This is now changing with the development of new experimental approaches that promise to generate validated miRNA: mRNA target data in a high-throughput fashion. One such technique is based on immunoprecipitation (IP) of miRISC proteins (RNA-induced silencing complex) [22], which probes the abundance of target mRNA bound to a mature miRNA. Although the actual miRNA: target pairs are not determined by this technique, the analysis of such datasets combined with the relative abundance of mRNAs and miRNAs [23] is critical for understanding the true cooperative interactions and model them in a quantitative way. Crosslinking immunoprecipitation (CLIP) is another recently developed and promising method. Both HITS-CLIP [24] and PAR-CLIP [25] allow for the miRNA target region to be determined in a narrow window on the mRNA. The high-throughput datasets these methods provide are ideal for developing methods to help a better understanding of the miRNA: mRNA targeting process. This paper addresses various issues related to miRNA targeting. First, we show that miRNA targets generally act additively in regulating mRNA expression. This was previously postulated but –to our knowledge– never put in test. Second, we show that using miRNA expression to appropriately weigh miRNA targets can result in efficient additive combinatorial models. To that extent we develop a novel thermodynamic model for miRNA binding, based on the Fermi-Dirac equation and we show that using this model improves prediction accuracy and target overlap of PITA [26] and miRanda [27]. Prediction efficiency of TargetScan [28] and mirSVR [29] is also improved by weighting multiple miRNA targets by miRNA expression and combining their target scores additively. Third, we use a Support Vector Machine (SVM) to combine the improved predictions of these four algorithms in ComiR (Combinatorial miRNA targeting), our novel algorithm that is designed to address the question of how likely is for a set of miRNAs with known expression levels to influence the expression of a given mRNA. The algorithm is tested on previously published miRISC protein IP independent datasets, ranging over three species, D. melanogaster, C. elegans and H. sapiens. Finally, we show that ComiR scoring scheme is suitable for ranking SNPs affecting miRNA binding. As a proof of principle, we used ComiR scoring for ranking the genes of the estrogen receptor (ER) pathway and we predicted that SNP rs17737058 would disrupt a miR-488-5p binding site in NCOA1. We subsequently confirmed the interaction in vitro and we found the SNP to be associated with decreased bone mineral density in two independent datasets. These results illustrate the predictive power of ComiR and that rs17737058 should be further studied as a risk factor for osteoporosis. This is an important ComiR application, since identifying functional DNA sequence variants will be critical in the era where plentiful information will become readily available for many diseases.
First, we examined whether cooperativity is an important factor in miRNA targeting. Let-7d, miR-30b and scrambled miRNA were transfected into human fetal lung fibroblast and microarray analysis was performed. Log2 of fold change (FC) in transcript abundance was measured with respect to scrambled miRNA transfection after let-7d transfection (FCl), miR-30b transfection (FCm) or let-7d/miR-30b co-transfection (FClm). Differentially expressed mRNAs were identified with SAM analysis (qval<0. 05) in each of the three experiments. In particular, 1,413 genes were significantly down-regulated post let-7d transfection (746 of them were only down-regulated in the let-7d transfection), 1,819 post miR-30b transfection (966 of them were only down-regulated in the miR-30b transfection) and 1,039 after both miRNAs were transfected. The 132 out of the 1,039 down-regulated genes were identified only in the co-transfection experiment, indicating that these were targets on which miRNAs act co-operatively. We used stepwise regression to test the potential miRNA additive effect on the fold change of the target genes. Our model is: where ℑ is the intercept. We performed a backward elimination by eliminating the less significant variables. The elimination of the variable from the regression model always causes a decrease of the resulting R2 value of the regressions. We found that the let-7d and miR-30b mRNA interactions are additive for most genes: small decrease of the R2 value was observed when the cross product term was eliminated (from 67. 1% to 67. 05%; both p-val<10−15). However, when the regression model was applied to the subset of the 132 genes that were found to be significantly down-regulated only in the co-transfection experiment, then elimination of the cross product term significantly impacts the R2 value (p-value changes from 10−3 to 0. 06). Notably, all weights obtained were <1, which probably reflects the saturation and competition effects on the miRISC machinery by the transfections [30]. The complete list of the regression coefficients is provided in Table S1. Figure S1 presents the distribution of binding sites for the genes that were down-regulated in the single transfection experiments only (red and yellow bars), in all three experiments (green bars) and in the co-transfection experiment only (black bars). We see that in all but the last category a higher percentage of genes has few (1 or 2) binding sites, whereas for the genes were co-operativity was observed had generally more sites (4 or 7). This indicates that co-operativity may be more important when a large number of sites is present in a gene. Based on the previous results, we developed ComiR (Combinatorial miRNA targeting), a novel algorithm that integrates the predictions of four top target prediction tools: PITA [26], miRanda [27], TargetScan [28] and mirSVR [29], into a single score. See Text S1 for more details about the implementation of the existing target prediction tools. First, each algorithm is run separately and for a given mRNA we identify all binding sites of each miRNA in its 3′UTR. Then, we incorporate miRNA expression and we additively combine the individual target scores using either Eq. 2 in the cases of PITA and miRanda (FD score), or Eq. 3 in the cases of TargetScan and mirSVR (WSUM score). By considering the miRNA expression in target score integration (Eq. 2 and Eq. 3) ComiR improves the efficiency of the corresponding algorithms as we showed above. The scores of the four tools for each mRNA are then combined through an SVM with linear kernel trained on the Drosophila RISC IP dataset [22] (see Materials and Methods). The target prediction score of ComiR is the class probability value computed by using the trained SVM model. Hence, ComiR scores range from 0 to 1 and higher scores correspond to higher probability of an mRNA being a functional target of the particular set of miRNAs. For the SVM implementation we used the ‘e1071’ R library. The general framework of ComiR is presented in Figure 2. The normalization step is a cross-species normalization of the score distributions (see Text S1 and Figure S2). ComiR is a multi-step algorithm that calculates the probability of an mRNA being targeted by a set of miRNA genes. First, depending on the prediction tool, it uses the FD score (Eq. 2) or the WSUM score (Eq. 3) to incorporate the expression of each miRNA and to combine the scores of individual targets of the miRNA set in a single score. The FD score combination is used when the primary target finding tool is either miRanda or PITA (where binding energies are used for the scoring); whereas the WSUM score combination is used for TargetScan or mirSVR. Subsequently, an SVM is used to incorporate the prediction scores of the four individual tools into a single probabilistic score characteristic of the probability that this set of miRNAs target a particular mRNA. We investigated how much the ComiR score combination of multiple targets contributes to the improved performance over a naïve combination of scores. Figure 4 shows that considering miRNA expression through the FD score or WSUM score always improves the SVM integration of any combination of target finders, although the degree of improvement depends on the particular tool combination and the dataset used. In all three datasets the performance of the SVM with the improved combined scores of PITA, miRanda and TargetScan (PMT in Figure 4) is almost as good as the SVM with all four algorithms (Figure S6). The smallest improvement the ComiR score offers is for the TargetScan/mirSVR combination. Also, in the Drosophila external dataset, the ComiR score combination substantially improves the performance of most dual tool combinations (except TargetScan/mirSVR) and all the 3- and 4-tool combinations. In the human PAR-CLIP dataset we see that in general when scores are combined with the naïve model, TargetScan has the best performance, followed by mirSVR, PITA and miRanda. However, with the ComiR model for incorporation of miRNA expression (FD score or WSUM score) PITA and miRanda become better than mirSVR. Finally, in all datasets we see that the improvement of prediction accuracy is higher when the FD score is used (i. e. , for of PITA and miRanda) than when the WSUM score is used (TargetScan and mirSVR). This indicates that the Fermi-Dirac model is indeed more accurate representation of the binding dynamics of miRNA: mRNA targeting, thus bringing the efficiency of PITA and miRanda closer to that of TargetScan and mirSVR. All the above indicate that incorporating miRNA expression in general and the Fermi-Dirac model in particular offer a very efficient way for combining individual target scores compared to the naïve model. The importance of the estrogen signaling, through the estrogen receptor α (ERα) pathway in bone maintenance is well established [35]–[38]. We asked the question whether ComiR could be applied to nine ERα pathway genes (Table S4) to identify SNPs in miRNA binding which could be associated with altered bone mineral density (BMD). Using dbSNP, we identified a total of 218 known SNPs in their 3′UTRs, 15 of which have minor allele frequency (MAF) greater than 10% (Table S4). There are about 400,000 SNP: miRNA pairs and nearly 29,000 of them correspond to the 15 SNPs with MAF>0. 1. We ranked the SNP: miRNA pairs based on the SNP induced change in ComiR binding probability score. We wanted to focus on high confidence targets (those with high binding probability), so we used Eq. 4 for the ranking. (4) Out of the ∼29,000 SNP: miRNA pairs we analyzed (those corresponding to the 15 SNPs of Table S4), 52 had >0. 01, with the miR-488-5p/rs17737058 (NCOA1) pair having the highest probability score among them. Importantly, rs17737058 is located in the center of the region matching the miR-488-5p seed sequence (Figure S7). We validated the effect of this SNP in the binding activity of miR-488-5p by first overexpressing a miR-488-5p mimic or mimic negative control (MNC) in U2OS-ERα cells and examining NCOA1 protein levels. Overexpression of miR-488-5p resulted in ∼50% relative reduction of NCOA1 levels (Figure 5A). As expected, no change was seen for NCOA3, a highly similar family member of NCOA1 that does not harbor a miR-488-5p target site (Figure 5A). Next, we examined if rs17737058 is sufficient to disrupt this regulation. Either the WT or rs17737058 3′UTR of NCOA1 was cloned downstream of the renilla luciferase CDS in the psiCHECK2 vector (Figure S8). The WT or SNP psiCHECK2 constructs were co-transfected with either an MNC or miR-488-5p. Consistent with the protein knockdown, WT renilla levels were reduced by ∼50% after overexpression of miR-488-5p (compared to MNC) (Figure 5B). However, SNP renilla showed a 30% attenuated miR-488-5p effect (Figure 5B). This indicates that the rs17737058 is sufficient to partially block the regulation of NCOA1 by miR-488-5p. Since NCOA1 is known to modulate the estrogenic effect in bone, we further investigated the role of this SNP in osteoporosis by examining existing GWAS data. While rs17737058 is not present in most SNP-chips, three other SNPs in the same linkage-disequilibrium (LD) block (rs719189, rs2083389, and rs9309308) are represented on the Affymetrix 100k SNP-chip used in the Framingham Heart Study bone mineral density (BMD) genome-wide association study. All three SNPs are significantly associated with decreased BMD specifically in women (Table S5), but not at genome-wide significance levels so they were not included in the original publication [39]. To assess further the association between these SNPs and BMD, we genotyped for three of these NCOA1 SNPs in an independent patient cohort. We utilized germline DNA from a prospective clinical trial (COBRA) in which BMD was measured as part of a comprehensive phenotype characterization [40]–[42]. In support of the GWAS studies, we detected significant association between the three NCOA1 SNP and decreased BMD in premenopausal women (Figure 5C). In postmenopausal women (excluding patients taking bisphosphonates for treatment of osteoporosis as potential cofounding factor), we observed the same trend although it did not reach significance. Interestingly, examination of the postmenopausal women revealed that SNP carriers were more likely than expected to have been on bisphosphonates at the time of study (Table S6) suggesting that women with this SNP may have an increased risk of loss of BMD. A likely explanation of these data is that SNP carriers may experience premature bone loss before menopause that is less evident after menopause onset. This may be due to the underlying function of NCOA1 as a coactivator of ERα and therefore the phenotypes may be more pronounced in the presence of estrogen. Together these data indicate that rs17737058 can be associated with decreased BMD likely through the disruption of miR-488-5p regulation of NCOA1. However, further experimentation with animal models is required to prove this association. Thus, we showed that applying ComiR to the ER pathway resulted in the identification of a SNP in a miRNA: mRNA pair with clinical significance in hormone response in bone.
In this study we presented a new method that advances the miRNA target prediction field in two key areas. One, it considers the quantitative effect of miRNA expression in target occupancy using a new thermodynamic model; and two, it quantitatively evaluates and combines the effect of target sites of multiple miRNAs on a given mRNA. We showed that our methodology improves the efficiency of popular target prediction algorithms as well as the overlap of their target datasets. Combining these improved predictions in a single probabilistic score (via SVM methodology) resulted in a new algorithm, ComiR, which when trained on Drosophila AGO1-IP data, it efficiently predicted targets of the differentially expressed set of miRNAs in Drosophila, C. elegans and human. By design, ComiR models the combinatorial effect of multiple miRNAs on a given mRNA. Thus, we expect that it will perform better in real-life examples, where multiple miRNAs are differentially expressed between two conditions. It is noteworthy that ComiR was proven to be more sensitive than any single tool. We attribute this to the nature of the SVM and the ComiR scoring system, which can elevate the targeting potential of a moderate affinity target predicted by any given tool if the miRNA is expressed in very high levels or if the mRNA contains multiple targets of this miRNA. Interestingly, TargetScan performed better than the other three algorithms and was competitive to ComiR on the single plotted point. However, ComiR remained the best of the algorithms tested. In addition, without evolutionary conservation TargetScan returns a binary outcome (target/no target) that does not allow for a threshold choice or for a scoring-based ranking of SNPs. The only exception is the human targets, where it provides a context score, which takes into consideration various features. Even in this case, ComiR was significantly more accurate in the human CLIP data and its better performance in terms of AUC was mostly in the region of high to medium false positive rate (Figure S3B). In any case, given the recent challenges of the assumption of evolutionary conservation of miRNA genes [43], [44] and their targets [45], methods that do not depend on evolutionary conservation may be proven a nice complement to the existing methods that do. Finally, in the high sensitivity area all algorithms seem to perform similarly, especially TargetScan and mirSVR. Notably, the SVM combination of TargetScan and mirSVR has the smallest improvement of ComiR vs. naïve combination of targets on this and the external Drosophila dataset (Figure 4). Silencing of miRNAs is usually considered a milder perturbation in the cell than the one caused by transfection, because the transfected miRNAs that are introduced en masse in the cell create a challenge to the capacity of the miRNA loading machinery [30]. There are currently no data to facilitate modeling of the miRNA affinity to the mRISC complex. So, although the transfection experiments were not the ideal test bed for ComiR, it was still proven to be more sensitive than the other four algorithms. We also showed that the ComiR score could be used to predict the effect of SNPs in single miRNA targets. Ideally, SNP and miRNA expression information should be obtained from the same individuals. We expect that such data will be routinely collected in the future. As a test case in this paper we analyzed the 3′UTR SNPs reported in dbSNP for the ERα pathway genes without having the benefit of knowing the miRNA expression levels. Based on ComiR top prediction, we postulated that miR-488-5p regulates NCOA1 and that rs17737058 reduces this regulatory effect. Since NCOA1 has a known role in maintaining BMD [35], [36], [38], we examined rs17737058 for an effect on BMD. Indeed, we found that this SNP was significantly associated with decreased BMD in two independent datasets. These results are strengthened by the observation that postmenopausal women carrying the SNP are more likely to be prescribed bisphosphonates than expected. This represents a rarely found example of a SNP disrupting a miRNA target site that results in a verified clinical phenotype. Interestingly, the clinical effect seems to be most evident in premenopausal women. We note that the FHS study was composed of two cohorts: the ‘Original Cohort’ (n = 159 women, mean age 77. 5) and the ‘Offspring Cohort’ (n = 487 women, mean age 58. 5) and the latter had three times more samples than the former. The other sub-studies within the osteoporosis GWAS meta-analysis [46] focused on older individuals, which might explain the failure to detect this association. Regardless, our data suggest that this SNP may identify premenopausal women at risk of osteoporosis and it should be considered a top candidate for further study and future development of personalized medicine therapeutic approaches. This is an important application of ComiR, because the identification of DNA sequence variants with a mechanistic functional role is becoming essential for the development of personalized medicine strategies. Notably, PITA and miRanda ranked the SNP very far down the list (their score didn' t practically change between wild-type and SNP sequence), mirSVR did not predict the pair, and TargetScan predicted the change but it offered no ranking. So, although in this case the lack of appropriate data did not allow us to take advantage of the full ComiR capabilities, it still provided a straightforward quantitative way to rank the SNPs affecting miRNA: mRNA interactions. We expect that in the future, when genotype and gene and miRNA expression data will be routinely collected from the same individuals, ComiR will be invaluable in identifying and ranking germline SNPs and somatic mutations that are associated to the disease. In summary, ComiR, the novel miRNA target prediction method we presented here, solves two important problems that hinder miRNA target prediction and offers a quantitative way to rank SNPs associated to miRNA binding. To our knowledge, this is the first algorithm that models the detailed thermodynamic interactions of miRNA binding dynamics in a cell and incorporates the quantitative effect of miRNA expression on multiple targets of multiple miRNA genes on the same mRNA. ComiR is by no means perfect. miRNA targeting is a complicated procedure and many characteristics still remain unknown. mRNA cellular localization or association with various RNA-binding proteins may influence miRNA binding. Interplay between miR-328 and RNA binding proteins has been previously reported [47], but the data are still scarce to allow for efficient modeling. RNA folding can also play an important role in miRNA targeting as it appears to do in other biological phenomena [48]. New high-throughput datasets will become available in the future and help elucidate these interactions. In that respect, ComiR is the first step towards a more complete modeling of miRNA: mRNA interactions, which is expected to be improved further as more types of high-throughput data become available.
The microarray data described in this work are deposited in the Gene Expression Omnibus database (GEO acc no. : GSE38530). A public web server for ComiR is available from the laboratory' s web page: http: //www. benoslab. pitt. edu/comir/.
MicroRNA genes (miRNAs) are small non-coding RNAs that regulate the expression levels of mRNAs post-transcriptionally. miRNAs are critical in many important biological processes, like development, and are important markers for many diseases. Identifying the targets of miRNAs is not an easy task. Recent developments of high-throughput data collection methods for identification of all miRNA targets in a cell are promising, but they still depend on computational algorithms to identify the exact miRNA: mRNA interactions. In this paper we present a novel algorithm, ComiR, which addresses a more general question, that is, whether a given mRNA is targeted by a set of miRNAs. ComiR uses miRNA expression to improve the targeting models of four target prediction algorithms. Then it combines their predicted targets using a support vector machine. By applying ComiR to single nucleotide polymorphism (SNP) data, we identified a SNP that is likely to be causally associated to osteoporosis in women.
genome analysis tools
genomics
genetic networks
biophysic al simulations
regulatory networks
biology
computational biology
genetics and genomics
transcriptomes
Novel Modeling of Combinatorial miRNA Targeting Identifies SNP with Potential Role in Bone Density