micrornas ( mirnas ) are nonprotein coding rnas of between 20 and 22 nucleotides that attenuate protein production by cleavage , translational inhibition , or sequestering of mrna in p bodies .
they are implicated in several different biological pathways , including plant and animal development , and cancer [ 24 ] . to better understand the role that mirnas play in these pathways , large datasets containing rna - seq , expressed sequence tags ( ests ) , and genomic sequences
are being investigated for new mirnas [ 5 , 6 ] . as these datasets grow in an ever increasing rate ,
the precursor folds back to base pair with itself to form a characteristic stem - loop structure .
the dicer protein cuts a short , double - stranded rna ( mirna : mirna * duplex ) from the precursor .
this double - stranded rna associates with the risc complex , where the mature mirna is retained while the mirna * is assumed to degrade .
the mirna - loaded risc complex is responsible for mrna cleavage , translational inhibition , or sequestering .  
many methods for predicting mirnas from sequence data have been reported [ 6 , 816 ] . some of these methods are directed specifically at ests   while others are specific for deep - sequencing data [ 10 , 11 , 17 ] .
several methods use machine learning approaches such as support vector machines ( svms ) [ 6 , 9 , 14 ] or decision trees [ 13 , 15 ] .
these methods typically incorporate information based on the predicted minimum free energy ( mfe ) of the precursor secondary structure .  
xue et al .   used a sliding triplet method to determine characteristic attributes within the triplets of mirna and non - mirna stem - loop structures for training an svm .
the authors reported prediction accuracy for mirna precursors of between 66.7% and 100% on a variety of plant and animal species , as well as viruses after training on known human mirna precursors .
it was trained using the same triplet element technique mentioned above , as well as other characteristics including dinucleotide shuffling .
the relationship between the mfe of the original sequence and its shuffled counterparts ( z statistics ) proved to be an important attribute for training this predictor .
the mipred method produced a sensitivity of 95.09% and specificity of 98.21% . both machine learning techniques described above agree on the high predictive value of specific triplet elements .
these are ( i ) three unpaired bases with a c in the middle position , ( ii ) three paired bases with a u in the middle , and ( iii ) three paired bases with an a in the middle .
approaches designed for deep sequencing data , such as mirtools , miranalyzer , and mirdeep2 [ 12 , 17 ] require the read count value to be known for predicting mirnas from rna - seq data .
although the read count may have good discriminative power , it prevents its application to genomic or est sequences .
the presence of the mirna * can be used by mirdeep2 . however , the mirna * is not always sequenced due to low expression level and degradation .   of the three de novo predictors of mirna precursors 
novomir   has the highest reported sensitivity and specificity of 83% and 99% , respectively .
hhmmir and triplet - svm report a sensitivity of 64% and 50% , respectively , for mirna precursors only .
novomir was trained using arabidopsis thaliana mirna for positive controls and uses trna , rrna , noncoding rna , mrna , and genomic sequences from a. thaliana as negative controls . the program reports the predicted precursor along with the start and end positions of any predicted mature mirna .  
one can ask , what data are necessary to correctly predict mirnas from sequences ? and ,
specifically , for deep sequencing data , can an accurate predictor be developed that does not require read count or rely on the presence of the mirna * ?
the predictive model described here is one of the most specific predictors ( 98.53% specificity ) , yet has the least requirements as it only requires that the candidate sequence occurs in a sequence region large enough to encompass the putative precursor .
fulfilling this single requirement allows the set of quantitative sequence properties ( attributes ) for the candidate sequence to be calculated and passed to the predictor . here ,
positive and negative controls from 18 plant species were used for training and model evaluation .
positive controls were taken from mirbase   while negative controls were taken from est sequences of each species .
the negative controls were collected based on the number and length of the known mirnas in the positive controls for corresponding species .
we chose randomly picked ests over specific subsets , such as mrna or various classes of noncoding rnas , to ensure that the predictor is not accidentally trained to detect the characteristics of the negative controls .
although this would increase sensitivity and specificity of the predictor , it would also introduce bias .
ests broadly represent protein - coding and nonprotein - coding genes including sequences resembling degraded mrna fragments that could exist in deep sequencing data [ 17 , 19 , 20 ] .
wet - lab validation is costly and time consuming ; therefore , accurately predicting that the candidate rna is truly a mirna is important .
our aim was to accurately predict when a sequence is not a mirna at the expense of missing a few true mirnas , which limits the number of false positives . to demonstrate the quality of any predictor
the positive controls are known mirnas and non - mirna sequences form the negative controls .
the negative controls were picked from ests of the respective species downloaded from tigr plant transcript assemblies .  
table 1 lists the 18 species by taxonomic group , and figure 1 is a flow chart of data collection and statistical validation . as stated above
, the formation of a hair - pin or stem - loop structure by the precursor is critical for mirna biogenesis . to obtain attributes for positive controls , the location of the known mirnas within the precursor stem loop
is required .   the known precursors from mirbase are only used to confirm two things about the mirnas used as positive controls : ( i ) the mirna location on the chromosome is inside a known precursor and not just a random match and ( ii ) the sequence is upstream or downstream of the mirna * ( in the 5 half or the 3 half of the precursor illustrated in figures 2(a ) and 2(b ) ) .
only attributes of mirnas for these correct chromosome locations and positions in the precursor are included as positive controls .
novoalign   was used to align mirna and precursor sequences to the respective genomes downloaded from plantgdb . at this point ,
the known precursors are no longer useful because all sequences will be treated equally from this point on regardless of whether they are positive controls , negative controls , or the assay sequences yet to be classified .   for any short sequence from the above three types ,
the precursor candidate was defined by first locating the mirna * candidate with the strongest duplex binding .
next , the region between and including the sequence and the mirna * candidate , along with 15  nt on both sides , was extracted .
this region defines the operative precursor region ( opr ) , which in real mirnas should form a stem - loop structure .
the search for the mirna : mirna * candidate duplex with the strongest binding was limited to a 300  nt window .
the 300  nt window was used because 95.80% of known plant mirna precursors are less than 300  nt ( table 2 ) . for positive controls
we have already defined that the correct window is upstream or downstream of the known mirna . for negative controls ,
both orientations are valid . after training , to use the model for prediction on sequences of unknown class , attributes from both upstream and downstream need to be collected and evaluated .
this is essentially the same as for negative controls but without prior knowledge of prediction outcome .
we use the opr instead of the precursor region as reported by mirbase to ensure equal treatment of all controls , and later for sequences of unknown class . simply , as unknowns and negative controls do not come with precursors
, we must define the opr equally for all , including the positive controls for training .  
figure 3 illustrates the relationship of the known precursor from mirbase to the opr defined using the method described above .
the computationally estimated mfe ( in kcal / mol ) for both the mirna : mirna * duplex and the opr structure are required attributes for training .
the rnaduplex function from the vienna rna package   was used to find the location of the mirna * and the binding energy of a mirna : mirna * duplex within the longer genomic sequences . with this information
rnafold from the vienna rna package was used to calculate the mfe for the opr , which is represented by deltag in table 3 .
rnafold also returns the secondary structure in dot - bracket format , representing the mismatch and base - pair matching , respectively , for the opr secondary structure .
the dot - bracket structure allows the longest run of uninterrupted matches and mismatches to be calculated , as well as the number of unpaired nucleotides that form the head of the loop .
all other attributes are based on characteristics of the mirna sequence , for example , the mirna base composition and the mirna : mirna * duplex measurements .
the rnaduplex function returns the number of matches and mismatches along with the duplex mfe .
it has previously been reported that the number of matching and mismatching base pairs within the precursor stem loop are important attributes for training the precursor predictors [ 13 , 14 ] .
the mirnas : mirnas * duplex information that includes base - pair matching and mismatching are also valuable attributes for predicting mature mirnas .
the attribute duplexenergynorm in table 3 is the mfe of the mirnas : mirnas * duplex normalized to the duplex length .
the minimum number of consecutive matches for each duplex is stored in the attribute minmatchpercent and is based on both sides of the duplex depending on the shortest side .  
attributes for negative controls were collected in a similar way as those for positive controls .
the randomly picked negative controls only qualified if they contained no ambiguity codes and had a length - normalized shannon entropy consistent with that of the positive controls .
the latter is used to avoid low complexity regions and to ensure that strong negative controls are collected .
the attribute shannonentropynorm in table 3 is the shannon entropy normalized by the sequence length .
not all sequences could be aligned ( e.g. , they may fall across an intron - exon boundary ) sufficient quantities were collected to make up for this expected loss . for negative controls , it does not matter whether the operational mirna is upstream or downstream of the mirna * as both orientations are equally valid and either the upstream or downstream attribute set is chosen randomly . in this way , the sequence is represented only once in that control set .
the final training set used for validation contains 2073 positive controls and 5306 negative controls .
more negative controls than positive were used so that the broadest set of non - mirnas is used to train against the known mirnas .
this increases the model 's ability to accurately recognize non - mirnas and to minimize false positive predictions .
once attributes have been collected for both positive and negative controls , validation sets can be produced .
the model was validated by calculating sensitivity and specificity based on leave - one - out cross - validation .
sensitivity is the ability of the classifier to identify positive results ( true mirnas ) , while specificity is the ability to distinguish negative sequences .
leave - one - out cross - validation can be described as putting all controls in a stack , taking the first one out and training with the rest .
only one of two possibilities will occur : the predictor is correct on the previously unseen control , or not .
the one previously removed is returned , and the next one is removed , and again training and testing are done until each control has been excluded and tested . when this is completed , the results are used to calculate sensitivity and specificity based on the following formulas : ( 1)sensitivity = number  of  true  positivesnumber  of  true  positives+number  of  false  negatives100,specificity = number  of  true  negativesnumber  of  true  negatives+number  of  false  positives100.only sequences that are not highly similar ( based on a cutoff value , see below ) to the test sequence will be allowed into the training set .  
it is important that sequences similar to the ones left out are also excluded to ensure the rigour of the leave - one - out approach .
a study on precursor prediction used the blastclust program to identify sequences of high similarity for exclusion from the training set .
however , we found that blastclust does not always ensure that sequences in two different clusters are sufficiently different , as any given sequence can only belong to one cluster .
we aligned each sequence to every other sequence and used their similarity to determine inclusion or exclusion from training .
in this way , each control sequence is excluded from training together with similar sequences .
all nonsimilar sequences are used to train , and this trained model is tested for its ability to accurately predict the class for the one excluded test case .
pools of positive and negative controls have each been pairwise aligned using the needleall program of the emboss package .
the results from needleall were used to determine the level of similarity between sequences in the positive controls set , then separately for negative controls .
needleall returns the similarity between two sequences , as well as other values . in this case ,
two sequences are considered similar if the similarity calculated by needleall is greater than 70% , based on the default setting .
the c5.0 program from rulequest incorporates a decision - tree machine learning method that we have used to train a mirna predictor using controls of known outcome . c5.0 and
the windows version see5 are improved versions of c4.5 that in turn descended from a program called id3 .
the training data used by c5.0 can be any combination of nominal attributes ( e.g. , the letter of the first nucleotide in the sequence ) and numeric attributes ( e.g. , mfe ) .
the output is a model in the form of if - then rules or decision trees for classifying cases of unknown outcome .
. producing models that are easy to understand can be useful when the goal is to discover the biological relationships between the attributes and class , rather than to classify unknown cases . in some investigations
the goal is to determine the attributes and cutoff values critical for separating the classes .
c5.0 can be applied to training data containing many thousands of cases with hundreds of attributes each .
the typical size of our training set is ~5294 cases , each with only 29 attributes .
if , after training , the model had a higher than acceptable false positive rate , a misclassification cost would be applied when retraining .
if the opposite was true ( i.e. , validation cheap and mirna rare ) , the misclassification cost could be adjusted accordingly .
for example , a misclassification cost would be applied to avoid missing true mirnas during training .
c5.0 supports a technique called boosting   that is used to improve classification results . boosting combines decision trees to produce and test new trees that may improve classification .
for example , the branches from one tree are swapped with those from another tree , and the two new classifiers are tested for increased predictive power .
different boosting values were tested to determine the level of improvement to the basic predictor .
figure s1 ( see supplementary material available online at doi : 10.1155/2012/652979 ) is a graph of roc space for 28 runs from 3 to 30 boosting trials .
this shows that , for this type of data , lower boosting values already produce good sensitivity and specificity , while higher boosting provides little improvement .
rulequest also provides an evaluator program that uses c5.0-trained classifiers to predict the outcome for data of unknown class .
we used this to record the predicted class for each case ( e.g. , mirna candidate ) along with the confidence value .
after training , the attribute usage information demonstrates the discriminative importance of each attribute .  
table 5 shows the attribute usage for one of many possible training runs of the classifier ; other training runs show similar usage .
several attributes , such as duplexenergy , minmatchpercent , and gc content , are required for all sequences to be classified .
genomic sequences and ests from cdna contain the base t , whereas mirnas from mirbase contain u. for this work , u was converted to t for processing attributes .
the attribute t% relates directly to the potential for g - u and a - u base pairing in the rna .
not surprisingly , the attribute named duplexenergy is important for training because it relates directly to the stability of the base pairing between the mirna and mirna * duplex .
this is evident by the 100% usage of duplexenergy and deltagnorm which represents the energy of the opr stem - loop structure normalized by its length .  
although this is an imperfect method , collecting data from all leave - one - out training runs can provide an overall view of an attribute 's discriminative power .
c5.0 has a built - in function called winnowing that , when applied during training , returns the percentage of increase in error if an attribute is removed from training .  
table 6 shows attributes with the largest average percentage of decline and therefore the level of importance for training the predictor . sensitivity and specificity
are critical values for assessing classifier accuracy ; values as high as 84.08% for sensitivity and 98.53% for specificity have been obtained .
the question remains whether this high specificity and sensitivity can also be achieved when the predictor is trained with different species than are used for prediction .
if all mirnas in each taxonomic category listed in table 7 are systematically excluded from training while including all others , how well does the predictor do when tested on the excluded category ?
results from these tests reveal how well the model predicts mirnas in plant species not included in the training set .
the ability of the classifier to correctly identify known mirnas ranged from 78% for the dicotyledon salicaceae to 100% for 7 of the 18 groups in table 7 .
all misclassified mirnas in the salicaceae are unique to the single species of poplar tree . in table 8 , four taxonomic groups are formed : monocotyledon , dicotyledon , lycopodiophyta , and embryophyta . when embryophyta is excluded from training ,
in contrast , when lycopodiophyta , embryophyta , and one of monocotyledon or dicotyledon are combined , the excluded set of monocotyledon or dicotyledon is recognized at 100% .
when lycopodiophyta and embryophyta are combined into one group named  primitive ,  training with the combination of any two groups produced a model that correctly classifies all mirnas in the excluded group , as shown in table 9 . for all exclusion experiments
this demonstrates that an accurate universal predictor of plant mirnas was produced when all groups , and therefore all mirnas , were combined for training .
identical mirna sequences are found in multiple plant species , often on multiple chromosomes where they are part of identical or different precursors .
although the mirna sequences can be identical , some of their attributes are not . at a minimum
leave - one - out sets are required for statistical validation by modelling the predictive accuracy of a control sequence that is unseen during training . for validation , excluding identical or near - identical mirnas is critical . when creating a predictor for unknown sequences , the classifier can be trained using all plant mirnas to ensure the best predictor possible .
the taxonomic exclusion tests demonstrate that training with the maximum number of known mirnas produces a classifier with the best cross - species recognition of mirnas .
in one training example , the set included 2278 positive controls , some of which are the same mirnas sequence but in different precursors and species .
the contrasting negative controls have 5680 short sequences randomly picked from ests . for this training set , from 3 to 300 boosting trials were run .
at boosting trials near 300 , only five negative controls were misclassified as mirnas during training .
although only wet - lab validation can confirm that a sequence is a mirna , two of the five are found in what look like very obvious mirna precursors .
the possible mirna precursors for these two est segments are shown in figures s3 and s4 of the supplementary data .
these two short sequences taken from the middle of randomly picked ests have yet to be validated as true mirnas .
one method of demonstrating that the predictor can be used for de novo prediction is to see how many known medicago truncatula ( mt ) mirnas can be accurately predicted from the precursor source .
this shows the quality of the predictor when applied to short segments taken from ests or chromosomes . for mt
, there are 342 unique mirnas from mirbase that occur in one or more of the 581 unique mt precursors .
therefore , short sequence segments were taken from the mirna precursors for only these lengths .
the segments were extracted by a sliding window method starting at the beginning of the sequence , taking 20  nt , moving one base over and again taking 20  nt .
this continues until the end where it becomes too short for a full 20  nt segment .
attributes are calculated for this set of short sequences , and then predictions are made by the trained classifier .
of the 342 unique mirnas , 309 are correctly classified as mirnas by the predictor for a score of 90.35% .
this demonstrates that the predictor can be applied to short segments taken from ests or genomic data as well as to short reads from rna - seq data .
we are currently working to achieve the speed increases required to analyze entire chromosomes . in some cases
of the 18 plant species used for training , only 17 known mirnas from seven species did not align to the respective genome . when these mirnas can be found in an est
, the est sequence can be used in place of a genomic sequence to collect attributes for prediction .
there is no genome listed for lja - mir167 in mirbase , and it is not found on any chromosomes tested for that species .
it is , however , found in an est [ genbank : bw598483 ] at position 43 .
a comparison between the stem loop from mirbase and our predicted opr for the mirna is shown in figure 3 .
this demonstrates that when a short test sequence is not found in a genome prediction can still occur if it is found in an est . during development
these negative controls were later pooled and the predictor applied to determine the number of potential mirnas in random samples of ests . out of 58,443 sequences , 633 ( 1.08% ) were classified as mirna .
as described in the methods section , only one mirna location ( either upstream or downstream of the mirna * ) was kept .
it is conceivable that if both were kept , the rate of mirna detection in randomly picked ests could double .
if some of these negative controls are true mirna , the specificity would be higher if they were removed from training .
while specificity is already high , this suggests that the predictor is more accurate than the value reported here
. the true sensitivity may also be higher than the reported value . despite using reported mirnas from mirbase for training
fifty mirnas across 10 species that were not classified as mirna were collected in a set .
twelve mirnas across four species out of the original 50 were classified as mirna when attributes from both upstream and downstream were tested .
closer inspection of the differences between attributes for the same mirna between the two sets revealed that in some cases the relative position of mirna and mirna * was different . in these cases
the precursor in mirbase was asymmetrical , and the automated location picking script returned the wrong location within the precursor .
when the correct location and attributes were tested by the predictor , it was classified as a mirna .
an example is gma - mir5034 mi0017906 , where the mirna is clearly located in the 3 end .
however , it has 69 bases on its downstream side and 60 bases on the upstream side , putting the mirna predominantly in the upstream half of the precursor . although some small set of erroneous attributes based on the wrong locations were included in the training dataset , the classifier correctly classified these as mirna when given both upstream and downstream attributes .
including attributes from incorrect locations for training produced a lower sensitivity value than if correct locations were used but did not diminish the predictor 's ability to correctly classify true mirna .
we have shown that a highly accurate universal plant mirna predictor can be produced by machine learning using c5.0 .
this predictor can be applied to any short sequence that aligns to a precursor candidate in a genome or transcriptome .
the source of the sequence for testing can be short reads from deep sequencing , or short segments taken from chromosomes or est sliding windows .
along with mirna prediction and prediction confidence level , the putative precursor is also produced ready for folding and visualization .
if used to scan a chromosome region , this predictor will reveal areas of high or low mirna density .
if applied to deep sequencing data the predictor reveals how often in a genome the short read exists , how many are predicted mirnas , and how many unique precursors contain that short - read mirna .
 micrornas ( mirnas ) are nonprotein coding rnas between 20 and 22 nucleotides long that attenuate protein production . different types of sequence data are being investigated for novel mirnas , including genomic and transcriptomic sequences . 
 a variety of machine learning methods have successfully predicted mirna precursors , mature mirnas , and other nonprotein coding sequences . 
 mirtools , mirdeep2 , and miranalyzer require  read count  to be included with the input sequences , which restricts their use to deep - sequencing data . our aim was to train a predictor using a cross - section of different species to accurately predict mirnas outside the training set . 
 we wanted a system that did not require read - count for prediction and could therefore be applied to short sequences extracted from genomic , est , or rna - seq sources . 
 a mirna - predictive decision - tree model has been developed by supervised machine learning . 
 it only requires that the corresponding genome or transcriptome is available within a sequence window that includes the precursor candidate so that the required sequence features can be collected . 
 some of the most critical features for training the predictor are the mirna : mirna duplex energy and the number of mismatches in the duplex . 
 we present a cross - species plant mirna predictor with 84.08% sensitivity and 98.53% specificity based on rigorous testing by leave - one - out validation . 
1. Introduction
2. Methods
3. Results and Discussion
4. Conclusion