MTar: a computational microRNA target prediction architecture for human transcriptome

Survey of existing solutions

A number of computational tools are available for animal and plant miRNA target identification. Of these, only MiRanda and RNAhybrid provide source code. Most of these approaches are based on evolutionary conservation and the presence of miRNA target sites in 3' UTRs of target mRNAs and their relatively better complementarities to 5' end of miRNAs. At the initial stages of microRNA target identification, researchers used near-perfect complementarities to predict miRNA targets for model species from plants. Tools like miRCheck [20], findmiRNA [21], PatScan [9] and mirU [22] are used for rapid prediction of miRNA targets in plants where perfect complementarities of miRNA and mRNA make the task easier. Though the targets for plant miRNAs can be identified on a genome-wide scale by searching for the ones that require a high degree of sequence complementarities, this cannot be used as such to find targets for animal miRNAs. The animal miRNAs tested till date pair imperfectly with their targets and act to control translation. Also the systematic analysis of the complete miRNA complement has confirmed the absence of targets with perfect or near-perfect sequence complementarities. So, target prediction in animal transcriptomes calls for more complex algorithms due to the imperfect complementarities of miRNA: mRNA pairs.

PicTar [23–25] predict miRNA targets in Drosophila and other species based on complementarities between miRNA and 3' UTR of mRNA sequence. PicTar used techniques like seed match, free energy calculation and species conservation. Its false positive rate has been estimated to be 30.0%.

TargetScan [26] is a tool used to predict miRNAs which bind to 3' UTRs of vertebrate transcriptomes. TargetScan could predict more than 451 human microRNA targets. TargetSanS [10], a modified version of TargetScan, omits multiple sites in each target and further filters the targets using thermodynamic stability criterion. Using this modified method more than 5300 human genes were predicted as possible targets of miRNAs [10]. The false positive rate varies between 22% to 31%.

MiRanda [12, 27, 28], a target prediction tool, relies on the evolutionary relationships between miRNAs and their targets. This tool focused to sequence matching of miRNA: mRNA pairs, by estimating energy of physical interaction. MiRanda was initially developed for predicting miRNA targets in Drosophila [27] and was later extended to find miRNA targets in mammals (human, mouse and rat) and Zebrafish [12]. The miRanda algorithm works by scanning for miRNA complementary pairs in the 3' UTR of an mRNA. Using this software, a large number of targets were identified including protein-coding genes in Homo sapiens. The false positive rate was estimated to be 24%.

The DIANA-microT [29] is a method based on the rules of single miRNA: mRNA pairing. It predicts targets which contain a single complementary site based on binding energies. MiTarget algorithm [30] combines thermodynamics based processing of RNA: RNA duplex interactions with the sequence analysis to predict miRNA targets. RNAhybrid is another computer program for predicting miRNA targets based on complementarities between miRNA and 3' UTR of coding sequence. This program was used to predict targets in Drosophila [11]. MovingTarget [31] is a program used to detect miRNA targets satisfying a set of biological constraints. Using this program more than 83 potential targets was predicted in Drosophila. MicroTar [32] is a program used to detect target sites in C.elegans, Drosophila and mouse by target complementarities and thermodynamic data. This algorithm uses predicted free energies of unbounded mRNA and putative mRNA:miRNA hetero dimmers, implicitly addressing the accessibility of the mRNA 3' UTR. This software is able to predict both conserved and non-conserved targets.

Anyway most of these existing tools make use of the complementarities in the 5' end of the miRNA alone. But MTar, the proposed computational method, can trap all the three types of targets (5' seed-only, 5' dominant and 3' canonical) and hence found to be more accurate. Multiplicity and cooperation which are common in animal miRNA: mRNA interactions are also handled effectively by MTar.

Training set

Experimentally verified microRNAs and their targets are required for training dataset preparation. In the data collection step, we excluded the sequences that were not verified by wet lab experiments in order to ensure the quality of the training data. We also excluded all the targets whose exact binding site could not be verified accurately. The dataset was selected based on two criteria. 1) The binding site of miRNA:target duplex should be known. 2) The target site sequences should match its corresponding references mRNA sequence in NCBI gene database.

The miRNAs were downloaded from miRBase database [33]. There are 706 reported human microRNA entries in the miRBase registry. The experimentally verified human microRNA targets were downloaded from Tarbase and miRecords registries. There are 609 human target sites for 443 genes by 107 miRNAs in Tarbase database [17]. The miRecords [34], contains 778 human records for 651 genes by 125 miRNAs. After filtered the target sites from these sources, the combined dataset consists of 882 human records for 741 genes by 138 miRNAs. The positive dataset needs three types of target classes (5' seed-only, 5' dominant, and 3' canonical). The classification was done based on the complementarity in the seed region. Seed region is the region of the nucleotides from 2-8 or 2-9 positions of the miRNA from the 5' end.

Randomly generated negative examples were not included in the training set, as such sequences are often found to interact with miRNAs due to their low signal to noise ratios as it is evident from previous studies [15, 26, 29]. Deletion of target positions on the target miRNA sequence can give a large number of negative examples. We generate a common negative dataset contains non miRNA:mRNA target sites with different parameter score. We collected examples with more than 4 mer matched at their seed part. Alignments of sequences in the training datasets were thoroughly checked in order to avoid ambiguities. The training set consisted of 350 examples with 150 positives and 200 negatives. The selection of positive dataset was based on the availability of experimentally verified target sites for the three class of targets in Tarbase registry. From these sequences three separate training datasets were created for three target classes (5' seed-only, 5' dominant, and 3' canonical). The training dataset contains 40 positives and 56 negatives for 3' canonical target class, 58 positives and 74 negatives for 5' dominant target class and 52 positives and 70 negatives for 5' seed-only target class.

Parameter selection

The method

Figure 2 summarizes the computational structure of MTar. To search for all possible alignments in each miRNA:mRNA pairs, segment of mRNA with length equal to the length of miRNA plus 10 nucleotide, starting from the first position, are selected. To locate miRNA targets, the miRNA sequence input is first aligned with the given mRNA target sequence using modified Smith-Waterman local alignment algorithm [41]. In the algorithm, gaps are allowed between the miRNA: mRNA pairs, but mismatches are preferred to gaps by giving a higher penalty for gaps. A scoring scheme in which each Watson-Crick pair (G:C and A:U) enjoys a score of 5, each G:U pair , a score of 1 and all others a score of -3, is employed. The gap opening amounts to -8 and a gap extension is less penalized with a score of -2. Based on this rule, a score S is computed for an alignment.

If the selected region does not belong to any of these classes, it is not considered for further processing. The predictions are strictly based on the class of the targets. Each of these segments is aligned with the miRNA to locate potential target candidates belonging to all the three different categories (5' seed-only, 5' dominant and 3' canonical) based on the nature of complementarities. After this, all overlaps are removed by filtering. Then, using the parameters listed in Table 1, a feature vector is formed by giving appropriate weights to the parameters for each of the candidates. Then they are submitted to an ANN classifier. The threshold levels for the parameters are different for targets of different categories. The validated potential targets are displayed along with their class to which they belong.

ANN classifier

In this study, artificial neural networks (ANNs) were chosen as the tool for miRNA target verification as they are powerful classifiers whose ability to cope with complex data and their potential for modeling data of high non-linearity [42, 43]. We used a feed forward three layer multi layer perceptron (MLP) for the classification of target sites. Calculated value of the parameters (Table 1) of miRNA:mRNA pair was selected as a 16 dimensional vector and fed into 16 input nodes. Output is set as either '1', if the output pattern is true and -1', if the output pattern is false. Hence there is one unit in the output layer of the ANN. The number of units in the hidden layer was chosen as nine, by trail and error. We used a sigmoid transfer function as the activation function. Back propagation algorithm was used to train the network. The training can be performed with use of several optimization schemes and there is access to exact partial derivatives of network outputs versus its inputs. The over fitting was avoided and this package makes automatic normalization of input data [34]. The learning rate and momentum were initially set at 0.2 and 0.8 respectively.

The training dataset is divided into two subsets. First subset (70% of the total training data) were used to train the neural network. Second subset was used to stop the training process once the model had reached the performance conditions like optimal error value thus preventing over training. Once the training is stopped, the efficiency of the model was further assessed by presenting another data subset, to determine the performance for unseen cases which were not involved in the training process. Optimization was done by repeating the process with different data subsets. The optimization needs nearly 500 epochs for this network. Three separate ANNs for each target class (5' seed-only, 5' dominant and 3' canonical) were trained to validate the target candidates of three different classes. Each ANN was trained with their training set and optimized the network as discussed above.

Performance evaluation

Extensive evaluation of MTar architecture was carried out using human genome. We could computationally predict 2663 target sites including 819 experimentally verified targets of 129 miRNAs (MFE below 17.0 Kcal/mol). For evaluation, the miRNA test data was downloaded from miRBase registry and the mRNA sequences from NCBI RefSeq and Biomart [44] cites.

We analyzed the performance of MTar using Receiver Operator Characteristic (ROC) curve which is shown in Figure 3. ROC is a plot of the true positive rate (sensitivity) against the false positive rate (1-specificity) for the different possible cutoffs of a diagnostic test. The test dataset contains experimentally verified target sites from Tarbase and miRecords databases. We confirmed that, any of the test data was not included in the training dataset. There are 190 positives and 200 negative examples selected in the test dataset. First, we tested MTar with all three features (positional, thermodynamic and structural features) combined. The area under the ROC curve was found to be 96%. The performance of MTar in terms of sensitivity (94.5%) and specificity (90.5%) is obvious. The threshold cutoff of MTar is 0.98 at this point. Then we investigated the effect of combining two features. For that, we tested MTar by taking only positional and thermodynamic properties at a time. The ROC area was decreased by 9.6%. Next we combined structural and thermodynamic properties together and then the ROC area decreased by 12.6%. Then we combined structural and positional properties together and tested MTar by the same dataset. The ROC area was decreased by 11.2%. From these experiments we could establish that all the three features (positional, thermodynamic and structural features) should be taken into account for improved miRNA target identification.

Role of parameter in target prediction

We also investigated the role of various parameters in miRNA target prediction. For ranking the features, we used Weka software [45]. This software is a machine learning algorithm based data mining software used for classification and visualization of dataset. Using a dataset, the features are ranked and are shown in Table 2. Positional features got high ranks compared to structural and thermodynamic parameters. We investigated the performance of MTar with different combinations of various parameters. The same training and testing dataset (used for ROC curve analysis) was used for analyzing top 3, 7, 12 and 16 features respectively (Table 1). The ROC curve for each test is shown in Figure 4. The prediction tool gave its best performance when all the features were included in dataset. Specificity and sensitivity were also significantly increased when the training was done with all the features given in Table 1.

Rank	Score	Parameter
1	84.3	Seed score
2	83.4	Positional pair score
3	81.0	Matrix score
4	80.6	WC pairs
5	78.8	Deviation positional score
6	77.2	Deviation matrix score
7	75.0	Out seed score
8	72.5	Wobble pairs
9	70.4	Normalized free energy
10	68.6	Mismatches
11	68.3	Hybridization energy
12	68.0	Free energy
13	65.5	Difference in hybridization energy
14	57.0	Length - bulge
15	52.4	Number-bulges
16	50.1	Proportion

Comparison with existing approaches

The MiTarget is one of the latest web based miRNA target prediction tool, which predict most of the 5' seed-only and 5' dominant miRNA target regions but it fails in identifying 3' canonical targets. This may be the due to picking the features from 5' part of the miRNA only while it is binding to the 3' part of the mRNA and ignoring the 3' part of the miRNA target.

After testing the MTar with a set of data, the test was repeated for MiRanda, TargetScan, RNA22, PicTar and MiTarget, with the same dataset. Figure 3 depicts the receiver operator characteristics curve for MTar along with a comparison of other three tools (MiTarget, MiRanda and TargetScan). The other three methods not provide cutoffs, so ROC generation was difficult. MiRanda shows a specificity of 82% by the test dataset. The specificity of PicTar nearly 70% and that of TargetScan is comes around 80% by the same dataset. We are not ploted the position of these two tools in Figure 3, due to their low specificity. The sensitivity of these tools are seen in Table 3. MTar gives an average accuracy of 92.8%, sensitivity as high as 94.5% and a specificity of 90.5% for the miRNA targets in the testing dataset.

The relative merit and comprehensiveness of MTar can be attributed to the following facts. First, the approach identifies all the three types of targets (5' seed-only, 5' dominant, and 3' canonical) in a single framework unlike other approaches. The second reason is the use of all the three properties (positional, thermodynamic and structural features) of the miRNA:mRNA duplex for the target identification. Third, multiple target sites are treated differently by set of user-defined biological constraints. A user can fix the free energy cutoff criteria of a miRNA:mRNA pair. Finally, our method is a desktop application for human transcriptome and also extensible to the other species.

A sample output of MTar with folding energy is shown in Figure 5. The user selected folding energy cutoff is below -17.0 Kcal/mol. The average MFE from the experimentally verified human miRNA target was calculated as -17.4 Kcal/mol.