- Open Access
RNAstrand: reading direction of structured RNAs in multiple sequence alignments
© Reiche and Stadler; licensee BioMed Central Ltd. 2007
- Received: 12 April 2007
- Accepted: 31 May 2007
- Published: 31 May 2007
Genome-wide screens for structured ncRNA genes in mammals, urochordates, and nematodes have predicted thousands of putative ncRNA genes and other structured RNA motifs. A prerequisite for their functional annotation is to determine the reading direction with high precision.
While folding energies of an RNA and its reverse complement are similar, the differences are sufficient at least in conjunction with substitution patterns to discriminate between structured RNAs and their complements. We present here a support vector machine that reliably classifies the reading direction of a structured RNA from a multiple sequence alignment and provides a considerable improvement in classification accuracy over previous approaches.
RNAstrand is freely available as a stand-alone tool from http://www.bioinf.uni-leipzig.de/Software/RNAstrand and is also included in the latest release of RNAz, a part of the Vienna RNA Package.
- Support Vector Machine Model
- Reading Direction
- Reverse Complement
- Pairwise Identity
- Consensus Structure
Genome wide computational screens for structured ncRNA genes in mammals [1–3], urochordates , nematodes , and drosophilids  resulted in tens of thousands putative structured ncRNAs. Functional and structural annotation of these predictions thus becomes a pressing problem. Evidence for evolutionary conservation of RNA structure alone usually does not distinguish very well between the two possible reading directions. This information, however, is crucial already for the most basic annotation information. Direction information is needed e.g. to determine whether a conserved RNA motif is intronic, located within a coding sequence or an untranslated exon, an independent ncRNA, or one of the many classes of small RNAs associated with other transcripts .
The RNAstrand tool is designed specifically to predict the reading direction of a multiple sequence alignment under the assumption that the alignment contains an evolutionary conserved RNA secondary structure. Our task at hand is a conceptually simple two class prediction problem for which we employ a support vector machine (SVM) . The basic idea is to devise descriptors that utilize both the small asymmetry in the energy rules  and the asymmetric effect of GU base pairs.
1.1 Selection of descriptors
Small differences in the measured folding energies between an RNA molecule and its reverse complement are captured by corresponding small asymmetries in the standard energy model used by thermodynamic folding algorithms [9, 10]. These differences distinguish the two reading directions even in the absence of GU pairs. In addition, GU pairs have an asymmetric effect in multiple sequence alignments: Suppose a particular pair of alignment columns exhibits a GC → GU substitution in one reading direction; this preserves base pairing and hence is consistent with a conserved structure. The reverse complement of the same alignment, however, displays a GC → AC substitution which is inconsistent with a conserved base pair. The patterns of structure conservation, and hence the consensus structure and its associated average folding energy, as computed by the RNAalifold algorithm , thus differ between the reading directions. In contrast, compensatory mutations, such as GC → AU do not provide strand-specific information.
The effects of both the asymmetries of the energy rules and of the GU base pairs are conveniently captured in terms of thermodynamic quantities, more precisely, in terms of the folding energies of the consensus structure and the individual folding energies of a set of aligned RNAs. These parameters can be computed much more reliably than quantities that have to be derived from predicted base pairs due to the limited accuracy of the structure prediction algorithms on individual sequences . We avoid the use of sequence motifs (e.g. ), since this bears the danger that the SVM is biased to the ncRNA families in the training set and fails to distinguish plus and minus strands of other structured ncRNAs.
Here we use:
1) Average of the folding energies of the individual sequences contained in the alignment, computed by the minimum energy folding program RNAfold of the Vienna RNA Package, version 1.6  (meanmfe).
2) Mean of the energy z-scores of the individual sequences contained in the alignment (meanz). The z-score is defined as z = (E - )/σ, where and σ are mean and standard deviation of the folding energy distribution of shuffled (permuted) sequences. We use here the same SVM-regression procedure as RNAz  to estimate the z-scores from the sequence composition to avoid the time consuming sampling of shuffled alignments.
3) Folding energy of the consensus secondary structure of the alignment computed by RNAalifold (consmfe). The parameter is defined as the optimal average of the folding energies that can be achieved when all aligned sequences simultaneously fold into the same structure.
4) Structure conservation index (sci), which is defined as the ratio of the consensus folding energy and the average of the folding energies of the individual sequences, i.e. sci = consmfe/meanmfe, . An sci close to 1 indicates perfect structure conservation, while alignments without structural conservation yield values close to 0. A more detailed discussion of the sci can be found in  in the context of RNA alignment.
The first two descriptors assess the thermodynamic stability of the folds, while the last two evaluate structural conservation.
The reading direction of a structured ncRNA can be identified by evaluating the differences of the above descriptors between both strands. To be precise, the difference Δx of descriptor x is defined as Δx = x+ - x-, where x+ denotes the value of x in reading direction of the input alignment and x- the value of x in the reverse complementary alignment. Hence, Δmeanmfe and Δmeanz capture the energetic differences between both strands, while Δconsmfe and Δsci describe the differences in structure conservation.
Note, that although sci = consmfe/meanmfe, i.e., these three quantities are not independent, this is not the case for their differences. Δsci cannot be computed from Δconsmfe and Δmeanmfe. Furthermore, for alignments where the structural conservation is very high in both reading directions the strand of the ncRNA cannot be inferred by Δsci alone. But the difference of consensus structure stability, which is measured by Δconsmfe may still predict the strand correctly.
Same holds for Δmeanz and Δmeanmfe. Both measure the folding energy differences of the individual sequences, but do not capture identical features of the input alignment nor can be transformed into each other. The mean z-score compares the average stability of individual sequences to a random control set. Whereas the mean of minimum free energies of individual sequences specifies the actual observed minimum free energies. The difference in z-scores describes the relative loss of stability compared to a random control set. It quantifies that the input alignment swaps from very stable to unstable between both strands. The difference in minimum free energy, on the other hand, is able to specify small changes in energies, which is needed to find the correct reading direction of the ncRNA in case both reading directions result in very stable structures. An example are miRNAs, which are very stable on both strands but are nevertheless successfully classified by RNAstrand. Hence, all four descriptors carry different information.
The significance of differences in folding energies depends on the number of sequences in the input alignment, denoted by n, and on sequence variation. The latter is conveniently quantified as the average pairwise sequence identity H of both reading directions.
We regard GU base pair fraction rather of the consensus structure than of the predicted structures of the single sequences, as the structure prediction of RNAalifold is based on evolutionary information of a set of sequences and hence produces a fold more similar to the real structure than RNAfold is able to predict from one single sequence. We did not introduce the difference of GU base pairs as a descriptor, because the error rate of such an descriptor depends largely on the correctness of the predicted secondary structure. Small errors in structure prediction have a large impact on the difference of GU base pairs. In contrast, the difference in structure stability and conservation regards all base pairs and hence depends only very weakly on the correctness of individual base pairs.
In summary, the SVM classification is based on seven descriptors, of which four, Δmeanmfe, Δmeanz, Δconsmfe and Δsci directly measure differences between the reading directions, while the remaining three, n, H, and λ GU provide information on the structure of the input alignment that allow the SVM to interpret the significance of strand differences.
1.2 Training of Support Vector Machine
Alignments for training were taken from the same sources as in  including representatives for rRNAs, spliceosomal RNAs, tRNAs, miRNAs, small nucleolar RNAs, nuclear RNaseP and SRP RNA. Sequence similarity in this data set ranges from 47% to 99% mean pairwise identity in alignments of 40 nt to 400 nt length and of 2 to 6 sequences. The detailed distributions of mean pairwise identity, length, number of sequences and GU base pair content are given in the supplementary material (see Additional file 1). A total of 5886 ClustalW alignments, approximately equally representing these ncRNA families, were used for training after removing alignments that were not recognized as structured RNA by RNAz in both reading directions. This data set was splitted into two subsets of equal size, namely the positive and negative training set. Alignments in the negative training set were transformed to the reverse complement and realigned with ClustalW as opposed to take just the reverse complementary alignment of the structured RNA.
The number of sequences a training alignment contains is limited to 6 as the SVM regression procedure to estimate the z-scores is trained with alignments of maximal 6 sequences . In case an alignment has more than 6 sequences a subalignment with optimal mean pairwise identity may be chosen with the perl script rnazWindow.pl  of the RNAz package.
We use libsvm 2.8  with SVM type C_SVC, a radial basis function (RBF) kernel, probability estimates and descriptor vectors scaled linearly to the interval [-1, +1]. The scaling avoids that descriptors which have a large variance dominate the classification. The values for the RBF kernel parameters γ and C were identified by a grid search in the parameter space applying grid.py of the libsvm 2.8 package with a 5-fold cross-validation on the training data. Maximal prediction accuracy is achieved with parameters C = 128 and γ = 0.5.
The SVM returns an estimated class probability p, that the ncRNA is found in the reading direction of the input alignment. We convert p into a score D = 2p - 1, so that D ≈ +1 means "RNA in reading direction of input alignment" while D ≈ -1 means "RNA is reverse complement of input alignment".
2.1 Testing the classifier
Classification performance is evaluated using 30920 automatically generated ClustalW alignments of 313 of the 503 ncRNA families from RFAM (version 7.0). All sequences attending at the training alignments were excluded from the test set. For each family at most 500 ClustalW alignments were randomly constructed each for 2 to 6 sequences, resulting in maximal 2500 alignments for a family. Since the alignments which were taken to train the SVM are no longer than 400 nt, have a minimal pairwise sequence identity of 60% and contain maximal six sequences, test alignments were created which meet the same criteria. For alignments which do not fall into those ranges probability estimates of the SVM need to be regarded with certainty. 8 families had no alignments between 40 and 400 nt and were hence discarded from the test set. 67 families are not included because they consist of only one or two sequences. 2 families had no sampled alignments with a mean pairwise sequence identity larger than 60%. Lastly, the sampled alignments of 113 families were not recognized as ncRNA by RNAz on at least one reading direction and were also discarded from the test data set. A list of families excluded from the test data can be found in the supplementary material (see Additional file 1). All alignments in the test set were used as positive test cases and their realigned reverse complements as negative test cases.
Evaluation of RNAstrand.
c = 0
c = 0.5
c = 0.9
Alignments classified as structured RNA by RNAz
euk. SRP RNA
Alignments not classified as structured RNA by RNAz
euk. SRP RNA
To evaluate the performance of RNAstrand on alignments which have not been identified as structured RNA by RNAz, we constructed a second test set which only consists of alignments not classified as structured RNA by RNAz in both reading directions. This resulted in 207 families meeting the criteria described in the first paragraph of this section. The corresponding distributions are shown in the supplementary material (see Additional file 1). For those alignments a dramatic decrease of structure stability and conservation is observed which leads to smaller descriptor values (see Additional file 1). Hence, the classification performance is worse compared to RNAz-positive alignments (Table 1). However, for the majority of alignments the correct reading direction was inferred.
Accuracies depending on different alignment features.
c = 0
N S = 2
N S = 3
N S = 4
N S = 5
N S = 6
50 ≤ sequence identity < 70
70 ≤ sequence identity < 80
80 ≤ sequence identity < 90
90 ≤ sequence identity < 100
40 ≤ length ≤ 100
101 ≤ length ≤ 200
201 ≤ length ≤ 300
301 ≤ length ≤ 400
The results highlight that our classification task has an intrinsic symmetry: the fraction of correctly classified alignments for the "plus strand" of a ncRNA should be similar to the accuracy of the "minus strand". However, we observe a small but noticeable bias to predict that the ncRNA lies in same reading direction as the input alignment (Table 1). The SVM model was trained with different alignments in the positive and negative training sets, which results in an asymmetric model. If the same alignments, but in different directions, were taken for training, the SVM model would be exactly symmetric. But training data should be independent in the different classes, hence we refrained from enforcing this exact symmetry to avoid potential overtraining artifacts. Another possibility to avoid asymmetry would be to take the averaged SVM decision values of both reading directions as the final decision. But this has an unknown effect on the probability estimates.
The distribution of decision values of the SVM is shown in Fig. 3. The majority of alignments were classified correctly. Most of them have large absolute decision values stating that they belong to the corresponding class with high probability. If RNAstrand is applied to shuffled alignments the decision values are more concentrated around 0, but most of them are still classified correctly. To explain this observation we checked which combination of descriptors performs best on shuffled alignments. We trained a SVM model for each possible descriptor combination and calculated the true and false positive rates at different decision levels by using plotroc.py of the libsvm 2.8 package . The corresponding ROC curves are given in Fig. 4 and indicate that except of Δmeanmfe all descriptors classify shuffled alignments randomly. Individual shuffled sequences, presumably by virtue of their base composition (see Additional file 1), still contain information on the reading direction of the structured RNA which is captured by Δmeanmfe. This observation implies that RNAstrand must not be used for alignments that do not contain structured RNAs. In other words, RNAstrand cannot be used to infer an ncRNA on the grounds that it returned a preferred reading direction for a non-structured input alignment. We could have also removed Δmeanmfe from the set of descriptors, because of this bias. However, due to its high sensitivity (Fig. 1) it seems preferable to keep it as descriptor, in particular since RNAstrand is designed to operate on structured RNAs only.
2.2 Comparison to naïve approaches
A naïve way to determine the likely reading direction is to score an alignment and its reverse complement using RNAz, EvoFold, or another tool for recognizing structured RNAs. This approach was taken e.g. in [1, 2, 4, 5]. A manual inspection of the data, however, showed that this approach is problematic in particular in those cases where RNAz scores are high for both reading directions. This is the case in particular for microRNA precursors, but also for many other small house-keeping ncRNAs.
Table 1 gives the accuracy of RNAstrand compared to this simple approach, i.e., taking the strand with the larger RNAz probability. RNAstrand yields for all ncRNA types an improvement. The largest increase of classification accuracy is observed for miRNAs, RNase MRP, tRNAs, nuclear RNaseP and IRES. Table 3 shows that the reading direction is classified correctly in the majority of test alignments by RNAstrand. The misclassification rate of the naïve approach is two times higher than that of RNAstrand.
Comparison of classification accuracies versus RNAz.
Naïve RNAz-based classification
Hence, we use instead the subset of known ncRNAs among the 48479 EvoFold predictions in human assembly hg17 .
A blast search with E < 1e - 10 against NonCode , Rfam , mirBase  and snoRNA-LBME-db  identified only 248 unique known ncRNA loci in human. (Note, that tRNAs and most snRNAs are multi-copy genes and hence were deliberately excluded from the data in ). To compare strand predictions of EvoFold with RNAstrand the multiz8way alignments of 202 loci, which are completely covered by a blast hit, were reconstructed. The majority (177) were identified to be miRNA precursors as most of the EvoFold predictions in ref.  are short conserved hairpins. The direction of the blast hit indirectly determines the strand of the known ncRNA when it is compared to the strand prediction of EvoFold. For 14 (13 miRNAs and 1 U6atac) loci the multiple alignments could not be reconstructed. The remaining 188 alignments were realigned and all which did not meet the prerequisites of RNAstrand were discarded: 15 alignments were shorter than the minimum length for which RNAstrand was trained with, 5 alignments had a mean pairwise identity smaller than 50%, and one alignment contained of too many gaps. This leaves 167 alignments for which the strand prediction of RNAstrand is compared to the strand prediction of EvoFold. Alignments containing more than 6 sequences were reduced to 6 sequences by rnazWindow.pl which optimizes the final alignment for a mean pairwise identity.
Comparison of classification accuracies versus EvoFold.
Naïve EvoFold-based classification
17 [17; 0]
19 [19; 0]
36 [36; 0]
RNA molecules and their reverse complements in general form fairly similar secondary structures . For individual sequences, small differences between plus and minus strand arise from small asymmetries in the energy model . In a multiple sequence alignment, GU pairs in an evolutionary conserved stem provide information on the correct reading direction since their reverse complement, AC, is not a canonical base pair. Nevertheless, it is a surprisingly hard problem to recognize the correct reading direction of a structured RNA from a multiple sequence alignment in practise. This is an important task in genome annotation, however, since without reliable strand information it is not even possible to determine whether an evolutionarily conserved secondary structure is located in an UTR or intron, or in an antisense transcript. The reading direction is also of obvious importance in context of recognizing class membership by means of short sequence motifs such as SMN-binding sites  or a Cajal body localization signal .
The RNAstrand tool presented in this contribution uses a SVM to predict strand information from a set of four thermodynamic features that can readily be computed for any multiple sequence alignment based on well-established energy parameters and dynamic programming algorithms. We show here that, together with basic information on the size, sequence and GU base pair variation in the input alignment, these features are sufficient to determine the reading direction of an RNA motif with an evolutionary conserved secondary structure. The tool RNAstrand achieves classification accuracies of 90% and above for most ncRNA families. On microRNAs, its performance is comparable to that of EvoFold. In applications to data from organisms for which not much genomic DNA has been sequenced, RNAstrand has the advantage that it does not require fairly accurate estimates of evolutionary distances as input.
The main area of application for a tool like RNAstrand is of course in large scale surveys for evolutionary conserved ncRNAs. RNAstrand achieves a 2-fold reduction of misclassifications on known ncRNAs compared to the naïve approach of determining the likely reading direction by comparing the scores of ncRNA detectors in both directions in the case of RNAz. It has therefore been integrated into the current release 1.0 of the RNAz package .
Project name: RNAstrand
Project homepage: http://www.bioinf.uni-leipzig.de/Software/RNAstrand/
Operating system(s): platform independent
Programming language: C
License: GNU GPL.
Restrictions to use by non-academics: Note that a license is needed to include source code from the Vienna RNA Package in commercial software projects.
We thank Ivo L. Hofacker and Stefan Washietl for valuable comments, and an anonymous referee for his ideas concerning the "mysterious" residual bias. This work was supported in part by the DFG Bioinformatics Initiative (BIZ-6/1-2) and the State of Saxony.
- Washietl S, Hofacker IL, Stadler PF: Mapping of conserved RNA secondary structures predicts thousands of functional noncoding RNAs in the human genome. Nat Biotechnol. 2005, 23 (11): 1383-1390.PubMedView ArticleGoogle Scholar
- Pedersen JS, Bejerano G, Siepel A, Rosenbloom K, Lindblad-Toh K, Lander ES, Kent J, Miller W, Haussler D: Identification and Classification of Conserved RNA Secondary Structures in the Human Genome. PLoS Comput Biol. 2006, 2 (4): e33-PubMedPubMed CentralView ArticleGoogle Scholar
- Washietl S, Pedersen JS, Korbel JO, Gruber A, Hackermüller J, Hertel J, Lindemeyer M, Reiche K, Stocsits C, Tanzer A, Ucla C, Wyss C, Antonarakis SE, Denoeud F, Lagarde J, Drenkow J, Kapranov P, Gingeras TR, Guigó R, Snyder M, Gerstein MB, Reymond A, Hofacker IL, Stadler PF: Structured RNAs in the ENCODE Selected Regions of the Human Genome. Gen Res. 2007Google Scholar
- Missal K, Rose D, Stadler PF: Non-coding RNAs in Ciona intestinalis. Bioinformatics. 2005, 21 (Suppl 2): ii77-ii78.PubMedView ArticleGoogle Scholar
- Missal K, Zhu X, Rose D, Deng W, Skogerbø G, Chen R, Stadler PF: Prediction of Structured Non-Coding RNAs in the Genome of the Nematode Caenorhabitis elegans. J Exp Zoolog B Mol Dev Evol. 2006, 306 (4): 379-392.View ArticleGoogle Scholar
- Rose DR, Hackermüller J, Washietl S, Findeiß S, Reiche K, Hertel J, Stadler PF, Prohaska SJ: Computational RNomics of Drosophilids. BMC Genomics. 2007Google Scholar
- Kapranov P, Cheng J, Dike S, Nix D, Duttagupta R, Willingham AT, Stadler PF, Hertel J, Hackermüller J, Hofacker IL, Bell I, Cheung E, Drenkow J, Dumais E, Patel S, Helt G, Madhavan G, Piccolboni A, Sementchenko V, Tammana H, Gingeras TR: RNA Maps Reveal New RNA Classes and a Possible Function for Pervasive Transcription. Science. 2007Google Scholar
- Cristianini N, Shawe-Taylor J: An Introduction to Support Vector Machines. 2000, Cambridge UK: Cambridge University PressGoogle Scholar
- Mathews DH, Sabina J, Zuker M, Turner H: Expanded Sequence Dependence of Thermodynamic Parameters Provides Robust Prediction of RNA Secondary Structure. J Mol Biol. 1999, 288: 911-940.PubMedView ArticleGoogle Scholar
- Lu ZJ, Turner DH, Mathews DH: A set of nearest neighbor parameters for predicting the enthalpy change of RNA secondary structure formation. Nucleic Acids Res. 2006, 34: 4912-4924.PubMedPubMed CentralView ArticleGoogle Scholar
- Hofacker IL, Fekete M, Stadler PF: Secondary Structure Prediction for Aligned RNA Sequences. J Mol Biol. 2002, 319: 1059-1066.PubMedView ArticleGoogle Scholar
- Doshi KJ, Cannone JJ, Cobaugh CW, Gutell RR: Evaluation of the suitability of free-energy minimization using nearest-neighbor energy parameters for RNA secondary structure prediction. BMC Bioinformatics. 2004, 5: 105-PubMedPubMed CentralView ArticleGoogle Scholar
- Carter RJ, Dubchak I, Holbrook SR: A computational approach to identify genes for functional RNAs in genomic sequences. Nucl Acids Res. 2001, 29 (19): 3928-3938.PubMedPubMed CentralGoogle Scholar
- Hofacker IL, Fontana W, Stadler PF, Bonhoeffer LS, Tacker M, Schuster P: Fast Folding and Comparison of RNA Secondary Structures. Monatsh Chem. 1994, 125: 167-188. 10.1007/BF00818163.View ArticleGoogle Scholar
- Washietl S, Hofacker IL, Stadler PF: Fast and reliable prediction of noncoding RNAs. Proc Natl Acad Sci USA. 2005, 102: 2454-2459.PubMedPubMed CentralView ArticleGoogle Scholar
- Gardner PP, Wilm A, Washietl S: A benchmark of multiple sequence alignment programs upon structural RNAs. Nucl Acids Res. 2005, 33: 2433-2439.PubMedPubMed CentralView ArticleGoogle Scholar
- Washietl S: rnazWindow.pl. http://www.tbi.univie.ac.at/~wash/RNAz/man/rnazWindow.html
- Chang CC, Lin CJ: LIBSVM: a library for support vector machines. 2001, http://www.csie.ntu.edu.tw/~cjlin/libsvmGoogle Scholar
- Youden WJ: Index for rating diagnostic tests. Cancer. 1950, 3: 32-35.PubMedView ArticleGoogle Scholar
- Liu C, Bai B, Skogerbœ G, Cai L, Deng W, Zhang Y, Bu D, Zhao Y, Chen R: NONCODE: an integrated knowledge database of non-coding RNAs. Nucleic Acids Res. 2005, 33: D112-D115.PubMedPubMed CentralView ArticleGoogle Scholar
- Griffiths-Jones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A: Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res. 2005, 33: D121-D124.PubMedPubMed CentralView ArticleGoogle Scholar
- Griffiths-Jones S: The microRNA Registry. Nucl Acids Res. 2004, 32: D109-D111.PubMedPubMed CentralView ArticleGoogle Scholar
- Lestrade L, Weber MJ: snoRNA-LBME-db, a comprehensive database of human H/ACA and C/D box snoRNAs. Nucleic Acids Research. 2006, 34: D158-D162.PubMedPubMed CentralView ArticleGoogle Scholar
- The ENCODE Project Consortium: Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature. 2007Google Scholar
- Fontana W, Stadler PF, Bornberg-Bauer EG, Griesmacher T, Hofacker IL, Tacker M, Tarazona P, Weinberger ED, Schuster P: RNA Folding Landscapes and Combinatory Landscapes. Phys Rev E. 1993, 47: 2083-2099. 10.1103/PhysRevE.47.2083.View ArticleGoogle Scholar
- Yong J, Golembe TJ, Pellizzoni DJBL, Dreyfuss G: snRNAs Contain Specific SMN-Binding Domains That Are Essential for snRNP Assembly. Mol Cell Biol. 2004, 24: 2747-2756.PubMedPubMed CentralView ArticleGoogle Scholar
- Richard P, Darzacq X, Bertrand E, Jády BE, Verheggen C, Kiss T: A common sequence motif determines the Cajal body-specific localization of box H/ACA scaRNAs. EMBO J. 2003, 22: 4283-4293.PubMedPubMed CentralView ArticleGoogle Scholar
- Washietl S: RNAz. http://www.tbi.univie.ac.at/~wash/RNAz/
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.