BACKGROUND
ribonucleic acid  plays various roles in the cell.
many functions of rna are attributable to their structural particularities .
rna motifs have been extensively studied for noncoding rnas , such as transfer rna , ribosomal rna , small nuclear rna , small nucleolar rna , etc.
<cit> .
more recently, small interfering rna  and microrna  have been under intensive studies  <cit> .
less well characterized are the structures in the un-translated regions  of messenger rnas   <cit> .
however, biochemical and genetic studies have demonstrated a myriad of functions associated with the utrs in mrna metabolism, including rna translocation, translation, and rna stability  <cit> .
rna structure determination via biochemical experiments is laborious and costly.
predictive approaches are valuable in providing guide information for wet lab experiments.
rna structure prediction is usually based on thermodynamics of rna folding or phylogenetic conservation of base-paired regions.
the former uses thermodynamic properties of various rna local structures, such as base pair stacking, hairpin loop, and bulge, to derive thermodynamically favourable secondary structures.
a dynamic programming algorithm is used to find optimal or suboptimal structures.
the most well-known tools belonging to this group are mfold  <cit>  and rnafold in the vienna rna package  <cit> .
similar tools have been developed in recent years to predict higher order structures, such as pseudoknots  <cit> .
on the other hand, rna structure prediction using phylogenetic information infers rna structures based on covariation of base-paired nucleotides  <cit> .
it is generally believed that methods using phylogenetic information are more accurate.
however, their performance critically depends on the high quality alignment of a large number of structurally related sequences.
tools that align biosequences , such as fasta and blast, are valuable in identifying homologous regions, which can lead to the discovery of functional units, such as protein domains, dna cis elements, etc.
<cit> .
however, their success is more evident in the study of dnas and proteins than of rnas.
this is mainly because the sequence similarity among dnas and proteins can usually faithfully reflect their functional relationship, whereas additional structure information is needed to study the functional conservation among rnas.
therefore, it is necessary to take into account both structural and sequential information in comparing rna sequences.
several tools are available that carry out rna alignment and folding at the same time .
the pioneer work by sankoff  <cit>  involves simultaneous folding and aligning of two rna sequences, and has huge time and space complexity .
foldalign  <cit>  improves the sankoff's method by  scoring the structure solely based on the number of base pairs, instead of the stacking energies; and  disallowing branch structures .
dynalign  <cit>  reduces the time complexity by restricting the maximum distance allowed between aligned nucleotides in two structures.
by taking into account local similarity, stem energy and covariations, perriquet et al.
<cit>  proposed carnac for pairwise folding of rna sequences.
ji et al.
<cit>  applied a graph-theoretical approach, called comrna, to detect the common rna secondary structure motifs from a group of functionally or evolutionally related rna sequences.
one noticeable advantage of comrna is its capability to detect pseudoknot structures.
in addition, algorithms using derivative-free optimization techniques, such as genetic algorithms and simulated annealing, have been proposed to increase the accuracy in structure-based rna alignment  <cit> .
for example, notredame et al.
<cit>  presented raga to conduct alignment of two homologous rna sequences when the secondary structure of one of them was known.
as shown in table  <dig>  most of these methods suffer from high time complexities, making the structure-based rna alignment tools much less efficient than sequence-based alignment tools.
tools that search for optimal alignment for given structures include rnadistance  <cit> , rna_align  <cit> , and rnaforester  <cit> .
rnadistance uses a tree-based model to coarsely represent rna secondary structures, and compares rna structures based on edit distance.
in a similar vein, rna_align  <cit>  models rna secondary structures by nested and/or crossing arcs that connect bonded nucleotides.
with the crossing arcs, rna_align is able to align two rna secondary structures, one of which could contain pseudoknots.
rnaforester extends the tree model to forest model, which significantly improves both time and space complexities .
in addition, methods using stochastic context free grammars  have been developed to compare two rna structures.
original scfg models  <cit>  require a prior multiple sequence alignment  for the training purpose, thus their applicability is limited to rna types for which structures of a large number of sequences are available, such as snorna and trna  <cit> .
however, rsearch  <cit>  and stemloc  <cit> , both based on scfg, are capable of conducting pair-wise structure comparisons with no requirement for pre-alignment.
rsearch uses ribosum substitution matrices derived from ribosomal rnas to score the matches in single-stranded  and double-stranded  regions.
stemloc uses "fold envelope" to improve efficiency by confining the search space involved in calculations.
the time and space complexities of these two tools are also listed in table  <dig>  furthermore, pattern-based techniques such as rnamotif, rnamot and patsearch  <cit>  have been used in database searches to detect similar rna substructures.
these tools represent rna structures by a consensus pattern containing both sequence and structure information.
one important advantage of these pattern-based tools is the ability of dealing with pseudoknots.
we present here a computationally efficient tool, called rsmatch, capable of both globally and locally aligning two rna secondary structures.
rsmatch does not require any prior knowledge of structures of interest.
it can uncover structural similarities by means of direct aligning at the structure level.
we demonstrate its application to database search and multiple alignment.
we compared rsmatch with three widely used tools, patsearch  <cit> , stemloc  <cit>  and rsearch  <cit> , demonstrating that rsmatch is faster or achieves comparable or higher accuracy than the existing tools when applied to a number of known rna structures, including simple stem-loops and complex structures containing junctions.
implementation
secondary structure decomposition
rsmatch models rnas by a structure decomposition scheme similar to the loop model commonly used in the algorithms for rna structure prediction  <cit> .
with this model, pseudoknots are not allowed.
our method differs from the loop decomposition methods in that it completely decomposes an rna secondary structure into units called circles .
when the secondary structure is depicted on a plane, a circle is defined as a set of nucleotides that are reachable from one another without crossing any base pair.
as shown in figure 1a, all circles are closed or ended by a base pair except the first circle , which always contains the 5'-most and the 3'-most bases.
various types of rna structures, such as bulge, loop, and junction can be represented by circles, as shown in figure 1a.
circles of an rna structure can be organized as a hierarchical tree according to their relative positions in the secondary structure, where each tree node corresponds to a circle .
this tree organization is informative to deduce the structural relationship among circles and reflects the structure particularities of the given rna secondary structure.
if two circles reside on the same lineage  in the tree, the circle appearing higher in the tree is called an ancestor of the other, and the latter is a descendent of the former.
as a result, in the context of the hierarchical tree, two distinct circles fall into one of the following two categories, in the order of decreasing closeness:  the two circles maintain an ancestor/descendent relationship, or  they share a common ancestor in the tree.
for example, in figure 1b, circle  <dig> is an ancestor of circle  <dig>  whereas circle  <dig> does not have ancestor/descendent relationship with circle  <dig> since they are not on the same lineage.
the double-stranded region or stem of a structure is decomposed into a set of "degenerated" circles, each containing only two base pairs.
as such, a stem of n bases in length will result in n -  <dig> consecutive degenerated circles.
since a base pair may have two associated circles; we name one circle "the parent circle" and the other "the child circle" according to their positions in the hierarchical tree.
for example, for the boxed c-g base pair in figure 1a, circle  <dig> is its parent circle and circle  <dig> is its child circle.
structure alignment formalization
given an rna secondary structure, we consider two types of structure components, single bases and base pairs, in the secondary structure.
to integrate both sequence and structure information, we introduce two constraints among the structure components: precedence constraint and hierarchy constraint.
the precedence constraint is defined based on the precedence order among structure components and the hierarchy constraint specifies the inter-component relationship in the context of the hierarchical tree described above.
the precedence order is determined by the 3' bases of individual structure components: the one with its 3' base closer to the rna sequence's 5'-end precedes the other.
for example, in figure 1a, the single base component u  in circle  <dig> precedes the base pair component c-g  in circle  <dig>
to capture the inter-component relationship within the hierarchical tree context, we need to map each structure component to a circle in the tree.
it is obvious that each single base can be mapped to a unique circle.
however, a base pair could be mapped to two alternate circles: one parent circle and one child circle.
to resolve this ambiguity, we always require mapping to the parent circle.
the inter-component relationship is then reduced to the inter-circle relationship of three types:  ancestor/descendent,  common ancestor, and  identical circle.
given two rna secondary structures a and b, where a, referred to as the query structure, has m structure components {a <dig>  a <dig>  ..., am} and b, referred to as the subject structure, has n structure components {b <dig>  b <dig>  ..., bn}, the structure alignment between a and b is formalized as a conditioned optimization problem based on the above two constraints: given a scoring scheme consisting of two matrices, one for matching two single bases and the other for matching two base pairs, find an optimal alignment between the two sets of structure components such that the aforementioned precedence and hierarchy constraints are preserved for any two matched component pairs  and .
in other words, the two structure constraints between ai and aj must be respectively equivalent to that between bi and bj.
this formalization has an implicit biological significance in that a single stranded region in one structure, if not aligned to a gap as a whole, will always align with a single stranded region in the other structure.
this alignment requirement is important because single stranded regions are usually treated as functional units in binding to specific proteins.
algorithmic framework
a dynamic programming algorithm is employed in rsmatch.
as with sequence alignment, the structure alignment could be either global or local.
the difference lies only in the setup of initialization conditions; the algorithmic framework is the same since both global and local alignments must preserve the two constraints described above.
a scoring table is established with its rows/columns corresponding to the structure components of the two given rna secondary structures.
we organize the rows/columns in such a way that the precedence and hierarchy constraints are combined and easy to follow in the course of alignment computation.
specifically, we sort the structure components of each structure according to the precedence order defined above.
it is straightforward that this arrangement of rows/columns makes the precedence constraint automatically preserved.
however, preservation of the hierarchy constraint is much more complicated and can only be accomplished in the derivative analysis for each cell  in the scoring table.
we will discuss the derivation when filling in the scoring table.
each cell of the scoring table represents an intermediate comparison between two partial structures corresponding to the cell's row and column components  respectively.
the partial structure with respect to a structure component c  is a set of structure components sc such that for any component a ∈ sc, the following three structure constraints between c and a must be satisfied:  a precedes c;  by the hierarchy constraint, a is not an ancestor of c; and  c itself is included in sc.
furthermore, since a base pair could appear in two circles, its corresponding partial structure could be divided into two smaller substructures: parent structure and child structure.
formally, given a base pair component c, the parent structure of c is the set of structure components pc ⊆ sc  such that for any component a ∈ pc, a's 3'-base is always 5' upstream of c's 5'-base; the child structure of c is the set of structure components lc ⊆ sc  such that for any component a ∈ lc, a's 5'-base is always 3' downstream of c's 5'-base.
it can be shown that pc ∪ lc = sc and pc ∩ lc = φ. examples of partial structures are given in figure 1c–1e.
as shown in figure 1c, for a base pair, its child and parent structures together constitute the whole partial structure for the base pair.
as we will see in the following discussions, the concept of a partial structure and its byproducts  form the kernel of our algorithmic framework.
we can solve the rna structure alignment problem progressively by aligning small structures and expanding each of them one structure component at a time until all structure components are covered.
preliminaries
cells in the scoring table are processed row by row from top to bottom and from left to right within each row.
by considering the row/column components, we have three types of cells:  a cell corresponding to two single bases;  a cell corresponding to one single base and one base pair; and  a cell corresponding to two base pairs.
for , each cell stores the score of aligning the partial structures corresponding to the cell's row and column components respectively.
for  and , we need to consider alignments involving the partial and child structures induced by the base pair components.
notice that the parent structures of the base pair components are excluded.
it can be shown that each parent structure pc of component c can always be considered as the partial structure sx of some other component x, which means we only need to consider child and partial structures in the alignment computation.
consequently, the above three types of cells have one, two and four alignment scores respectively.
a scoring scheme is required to score the match of two structure components.
we define the scoring scheme as a function g where a and b represent two structure components that are matched with each other.
another important aspect of the alignment algorithm is to penalize the match involving gap.
in the course of computation, one structure component  could match with a gap or a whole small structure  could match with a large gap.
intuitively, the larger the gap is, the heavier the penalty will be.
in our implementation, we set an atomic penalty value, denoted as u, for the smallest gap equivalent to a single base.
the penalty value for a large gap is proportional to its size in terms of the number of bases matched with the gap.
let a* be a small structure in the query rna structure a and b* a small structure in the subject rna structure b. the score obtained by aligning the two structures a* and b*, denoted as f, is , where g represents the total number of gaps in aligning a* and b*.
initialization
we assume that the row components  are from the query rna structure a and the column components  from the subject rna structure b. we focus on global alignment here; initializations for local alignment can be derived similarly.
the initialization conditions deal with the cases where at least one of the structures under alignment is an empty structure φ. this is equivalent to setting up the 0th row/column in the scoring table.
as discussed above, each base pair component has two small structures to be considered: a child structure and a partial structure.
thus, the aforementioned three types of cells have one, two and four initialization scores respectively.
for a given structure component x , let sx represent its partial structure.
if x is a base pair, we also use lx to represent its child structure.
we have f =  <dig>  furthermore, for any structure components a and b, f = |sa|·u, f = |sb|·u, if a and b are base pairs, f = |la|·u and f = |lb|·u where |·| represents the cardinality of the respective set.
filling in the scoring table
the simplest cell type is the one whose row  component is a single base a .
let ap denote the structure component that precedes a by precedence order established before.
formally, in matching the partial structure sa with the partial structure sb there are only three possibilities:  a is aligned with b;  a is aligned with a gap; and  b is aligned with a gap.
thus the score of matching sa with sb can be calculated by equation .
the second cell type is the one formed by one single base and one base pair.
there are actually two symmetric subtypes where either a or b is a base pair.
since the analysis is identical, we only focus on the former case where a is a base pair.
as discussed before, besides the partial structure sa we have to consider the child structure la for the base pair a. thus, for this type of cells, we have to compute two alignment scores.
by the principle of dynamic programming, the smaller size problem needs to be solved before the larger size problem.
thus we first find the structure alignment between the child structure la and the partial structure sb.
there are only two possibilities:  the single base component b is aligned with a gap; and  the base pair a is aligned with a gap .
therefore we have
in aligning the partial structure sa with the partial structure sb, to preserve precedence and hierarchy constraints simultaneously, there are only three possibilities:  the single base b matches with a gap;  the partial structure sb matches with the child structure la;  the partial structure sb matches with the parent structure pa .
thus,
for the third cell type, a is a base pair and b is a base pair.
we need to compute four alignment scores because each base pair corresponds to two structures: one child structure and one partial structure.
while aligning the child structure la with the child structure lb, it is clear that
since both a and b are the last components in the respective child structures by precedence order.
equation  gives the alignment score between the partial structure sa and the child structure lb.
the first case corresponds to that b is aligned with a gap.
if b does not match with a gap, it can be shown that, to preserve both precedence and hierarchy constraints, the second and third cases in equation  cover all possible situations.
similarly, we can calculate the score of aligning the child structure la and the partial structure sb as shown in equation .
in aligning the partial structure sa with the partial structure sb, there are five possibilities:  the parent structure pa is matched with the parent structure pb and the child structure la is matched with the child structure lb;  the child structure la is matched with gaps;  the child structure lb is matched with gaps;  the parent structure pa is matched with gaps; and  the parent structure pb is matched with gaps.
therefore
data sets
all experiments  were carried out on a linux system with two  <dig>  ghz intel processors and  <dig> gb memory.
a human utr structure database was constructed as follows.
we downloaded  <dig>  human refseq mrna sequences  from national center for biotechnology information .
each refseq sequence containing utr regions, as indicated by refseq's genbank annotation, was processed to extract its 5'utr and 3'utr sequences.
for each utr sequence, we took a  <dig> nt subsequence at every 50th nucleotide position from 5' to 3', making consecutive subsequences overlap with one another on a  <dig> nt segment.
subsequences shorter than  <dig> nt, e.g.
at the 3' end, were also kept.
using the vienna rna package's rnasubopt function with setting "-e 0", we then folded all obtained sequences to form the structure database.
for any given rna sequence, the setting "-e 0" resulted in multiple rna structures all having the minimum free energy.
the final database contained ~ <dig>  rna secondary structures.
the structural patterns of a histone 3'utr stem-loop structure  and an iron responsive element  were used in this study, based on their specifications in the utrdb database  <cit> .
three tools, patsearch  <cit> , stemloc  <cit>  and rsearch  <cit> , were employed for comparison purposes.
the efficiency of these tools was measured by cpu running time.
the performance of each program was assessed by specificity and sensitivity.
specificity was calculated as tp/ and sensitivity as tp/, where tp was the number of true positives, fp the number of false positives, and fn the number of false negatives.
to test the applicability of rsmatch to complex structures, we downloaded rna families from rfam  <cit> .
we only chose those families that had more than  <dig> seed rnas and its consensus sequence length is no longer than  <dig> nucleotides.
we had  <dig> families in the final data set.
for each family, we randomly selected one member rna as the query rna and obtained its structure from rfam.
we then randomly chose  <dig> subject rnas in the same family.
here we intentionally introduced noise by extending each subject rna sequence with its adjacent sequences at both 3' and 5' ends to make the total length three times its original one.
RESULTS
studies of stem-loop structures in utrs
using our proposed algorithm, we first studied rna motifs in utr regions of human mrna sequences.
a well-known fact is that the accuracy and efficiency of rna folding programs will decrease significantly when the sequences to be folded become very long.
satisfactory performance is usually obtained when the sequences have moderate lengths, i.e.
one hundred nucleotides.
thus, we used a moving window scheme to get subsequences of  <dig> nt and folded them using the vienna rna package .
in the rsmatch package, this subsequence length is a user-defined parameter.
since the nucleotide conservation in the single-stranded region of an rna sequence may differ from that in the double-stranded region, we used two scoring matrices, one for substitutions among single bases and the other among base pairs.
this type of scoring scheme was also used in other studies  <cit> .
theoretically, the scoring matrix for single bases is a  <dig> ×  <dig> table for all types of substitutions of single nucleotides, and the one for base pairs is a  <dig> ×  <dig> table for all types of substitutions of base pairs.
however since we used the vienna rna package, only six types of base pairs were observed in our studies, i.e.
watson-crick base pairs a-u, u-a, g-c, c-g, and wobble base pairs g-u and u-g. values used in the two matrices were empirically chosen so as to conform to the general understanding of the sequence and structure conservation of rna motifs, as follows.
mutations in the double-stranded region may not be detrimental to rna's function if the mutated sequence still preserves the same secondary structure.
therefore base pair substitutions were rewarded with a positive score, instead of a penalty.
a sequence in the single-stranded region may be important for rna's function, such as binding to proteins, and thus mismatches were penalized.
to process gaps we used an arbitrary function u × l, where u was the atomic penalty value for a gap that is one single base long and l is the length of the gap in terms of the number of bases matched with the gap.
in our experiments otherwise stated explicitly, the u was empirically set to - <dig> and changing the u value did not change the qualitative conclusion made in the paper provided that the absolute value of u was greater than any positive score in the scoring matrices.
users can freely change the u value when applying rsmatch to their own data set.
we tested our program with a query sequence containing an iron response element .
the ire motif is a bipartite stem-loop structure containing ~ <dig> nucleotides.
two alternative types of ires have been found, which differ in the middle region  <cit> .
type i has a bulge, whereas type ii has a small internal loop.
ires have been found in both 5' and 3' utrs of genes that are involved in iron homeostasis in higher eukaryotic species.
they interact with iron regulatory proteins  and play a role in rna stability and translation.
using a subsequence in the 3'utr of transferrin receptor  that contains an ire motif, we searched the utr structure database described in implementation.
a list of top hits is shown in figure  <dig>  the best hit of the search is the query structure itself, as expected.
other regions of the same mrna and regions of other rnas are also found to have homologous structures with the query.
as clearly shown in the result, the region containing the ire motif, which is from about the 30th nucleotide to about the 60th nucleotide of the query structure, has been located by the rsmatch program, indicating that a local optimal alignment has been achieved.
among the top  <dig> hits, several sequences are known to have ires, such as several regions in the 3'utr of transferrin receptor  and the 5'utr of solute carrier family  <dig> protein .
other top hits have not been shown so far to have ires.
it is not known if some of them are novel ire-containing rnas and the definitive answer will await wet lab validation.
the output shows detailed alignment and related information, including the numbers of bases in the single-stranded and double-stranded regions, and the percentages of identity in single-stranded and double-stranded regions.
rsmatch can also accept pattern-based rna structures  to search a structure database.
since a pattern-based search method has an intrinsic primitive scoring scheme by using degenerate bases, we used simplified binary matrices as the equivalent to score an alignment.
in the matrices, the match of a pair of structure components , including those containing degenerated bases, was given a score of  <dig>  a mismatch was penalized by a score of - <dig>  and the atomic gap penalty u was set to - <dig>  to allow variability in single-stranded and/or double-stranded regions for a structure pattern, we introduced a wildcard "n " to represent optional single base component  and base pair component .
the meaning of "n" is identical to the iub code "n" except that the matching score for both structure components "n" and "n-n" is always zero regardless of whether they are aligned with a structure component or a gap.
two rna motifs were used to test our method, namely a histone 3'utr stem-loop structure  and ire.
hsl <dig>  which resides in the 3'utr region of histone mrnas, has a typical stem loop structure with two flanking tails .
both the stem and the flanking sequences are important to bind with a stem-loop binding protein , which controls the pre-mrna processing and stability of histone mrnas  <cit> .
in contrast to the hsl <dig> motif, ire is relatively flexible in length and in nucleotide composition in its stem region .
we compared our program with patsearch  <cit> , a widely used tool that searches a sequence database for sequence and structure patterns.
using the hsl <dig> motif and utr sequence database, patsearch found  <dig> hits whose locations were presented in table  <dig>  among them, one is a false positive .
therefore the specificity  of patsearch is very high.
this is attributable to the precise specification of the hsl <dig> pattern.
however, if a pattern description is too precise, it may lead to the "overfitting" problem.
this problem prevents the tool from finding slightly divergent structures, thus lowering the tool's sensitivity.
indeed, several histone genes were not detected by patsearch, including two histone genes  which were found by rsmatch among its top  <dig> hits.
several other histone genes appeared among the top  <dig> hits of rsmatch .
this indicates that by gaining specificity, patsearch loses sensitivity for hsl <dig>  since rsmatch gives a score to each alignment, different cutoffs can be used for selecting top hits .
it seems that newly detected true positives are heavily outnumbered by false positives as rsmatch relaxes its cutoff value.
however, with some properly chosen cutoff, i.e.
<dig>  rsmatch could still achieve a comparable specificity with patsearch.
one possible explanation of getting high false positives for rsmatch could be that, with respect to the particular case of the hsl <dig> motif, its secondary structure conformation might be too pervasive in rna sequences to be used as a discriminative feature.
this could point out a problem concerning rsmatch's current scoring matrices, which need to be fine tuned to improve the tool's specificity.
good tuning could be achieved by setting up the scoring matrices through learning from a training data set.
one interesting observation, however, was that rsmatch and patsearch agreed perfectly upon the hsl <dig> locations for almost all of the true positives they found.
using the ire motif, we performed further comparisons between rsmatch and three other tools: patsearch  <cit> , stemloc  <cit>  and rsearch  <cit> .
we used default parameters for rsearch; for stemloc, the fold envelope was set to  <dig>  instead of using the large utr structure database described in implementation we constructed a small test data set to expedite the comparison process.
first, we used patsearch to search human utr sequences for ire motifs.
then for each hit sequence we selected its corresponding mrna's 3' or 5' utr sequence.
following the same folding process as discussed in implementation, we folded these utr sequences to form the test data set.
totally, patsearch found  <dig> hits, among which  <dig> were known true positives.
therefore patsearch's specificity was ~33%.
these hits were from  <dig> distinct mrna sequences.
we assumed that patsearch had a 100% sensitivity.
we extracted the 5' and 3' utr sequences from the  <dig> distinct mrnas and obtained  <dig> utr sequences.
we then folded the  <dig> utr sequences to get a small test data set, which contained  <dig> structures.
using a known ire-containing structure , which was one of the  <dig> true positives found by patsearch, as the query, we searched the small test data set.
table  <dig> shows the results we obtained.
since rsearch accepts sequences only, it was tested using only the primary sequence information in the test data set.
except for the ire-containing structure nm_ <dig>  which was one of the  <dig> true positives found by patsearch, and the query itself , all tools agreed on the ire locations for the other seven true positives without salient discrepancy.
it was found that nm_ <dig> was not properly folded to exhibit the existence of the ire motif.
rsmatch has the best specificity by ranking all seven true positives within its top  <dig> hits with only one false positive .
rsearch is close to rsmatch by ranking all seven true positives within its top  <dig> hits with one false positive .
in contrast, stemloc gives five false positives within its top  <dig> hits.
setting different cutoff values yields different specificity and sensitivity for each tool.
the point of balanced specificity and sensitivity appears at the cutoff value of  <dig> for all three tools.
with this cutoff value, the specificity of rsmatch and rsearch tied at 7/ <dig> × 100% =  <dig> %.
this is better than the specificity of patsearch  and the specificity of stemloc .
the sensitivity of rsmatch, rsearch and stemloc is  <dig> %,  <dig> % and 50% respectively.
it is worth noting that rsmatch runs ~30% faster than rsearch; it took rsearch  <dig> seconds to search the whole data set of  <dig> structures while rsmatch used only  <dig> seconds.
consequently, rsmatch would be suitable for analyzing large data sets.
it should also be pointed out that rsmatch permits wildcards in database searching and structure matching, which are not supported by rsearch or stemloc.
performance on complex structures
we further tested how accurate rsmatch is for complex structures.
to this end, we downloaded rna structures and sequences from the rfam database .
we used  <dig> rna structure families, each of which has more than  <dig> seed sequences and has the consensus sequence length less than  <dig> nucleotides.
for each rna structure family, we randomly selected a structure and searched against  <dig> randomly selected sequences belonging to the same family.
to reflect real world scenarios, we extended rna sequences at both 5' and 3' ends so that the length of a subject sequence is three times that of the original one.
to ensure that the folded structures are long enough to fully contain the structure being investigated, we required the moving window size to be  <dig>  times the length of the query rna sequence.
furthermore, to include suboptimal structures, we used all structures with free energy within  <dig>  kcal/mol above the minimum one.
compared with hsl <dig> and ire, the  <dig> query structures we used were much more complex, with average length of ~ <dig> nt and more than 70% of them comprised of nested loops and conjunctions.
to assess the accuracy, we used a measure called structure coverage, denoted as p, which is calculated by the following formula: p = |qalign|/max, where |qalign| and |salign| are the lengths of aligned portion of query rna and subject rna, respectively, and |q| is the length of query rna sequence.
as shown in figure 5a, even though rsearch has slightly more points clustered around high coverage , the overall difference between rsmatch and rsearch is not significant.
in addition, the difference between rsmatch and rsearch do not seem to be related to structure size or complexity.
this result indicates that rsmatch has the ability to process complex structures.
we also selected 5s rrna for further detailed tests.
5s rrna has a length of ~ <dig> nt, which contains several types of rna structures, including hairpin, internal loop, bulge, and junction.
there are  <dig> sequences in the 5s rrna family, allowing us to carry out a thorough analysis.
we randomly chose one 5s rrna as query structure and ten others as subject sequences for alignment.
this process was repeated  <dig> times.
the performance comparison of rsearch and rsmatch is shown in figure 5b.
for 5s rrna, rsmatch outperforms rsearch in discovering the complete structure more frequently.
an exemplary alignment is shown in figure 5c–5e.
running efficiency
by dynamic programming, the running time of computing an alignment equals the number of writing operations needed to fill the scoring table.
thus the time complexity of rsmatch is o, where m  is the number of structure components in the query  rna structure.
to test the scalability, we downloaded the seed sequences for 5s rrna family from rfam and randomly selected one annotated structure as the query while folding the rest sequences to prepare the structure database as discussed above .
we plotted the rsmatch running time versus the database size.
the program was run  <dig> times and the result is shown in figure  <dig>  the nearly perfect linear growth of the running time gives an empirical proof that the algorithm's time complexity is bounded by o.
multiple structure alignment and iterative database search
we also extended rsmatch algorithm to conduct multiple structure alignment.
an example using ire is shown in figure  <dig>  while the alignment algorithm is the same, the multiple alignment function uses a position-specific scoring matrix .
for a given set of structures, the multiple alignment function first identifies the best alignment of two structures, and builds a pssm.
the pssm is then used to search for the closest structure in the rest of the set.
a flowchart of multiple structure alignment is shown in figure 7a.
if the alignment score of a structure to the pssm is above a cutoff , it is selected and its structure is used to update the pssm.
this step is iteratively conducted until no structures have alignment score above the cutoff.
in a sense, this method employs an implicit guided hierarchical tree using the average value for joining nodes.
as an example, from our human utr database we selected  <dig> ire-containing structures and randomly chose other  <dig> none-ire structures to form a small dataset and run rsmatch against it.
the output is shown in figure 7b.
the final result is in stockholm format for multiple structure alignment.
conceivably, when the given set of structures is a large database, the multiple structure alignment function of rsmatch in effect conducts iterative search for finding similar structures.
discussion and
CONCLUSIONS
the work presented here is intended to provide an efficient tool to directly perform structure alignment and search of rna secondary structure databases.
its capability to carry out multiple structure alignment and iterative database search can potentially be used to uncover rna motifs ab initio.
for example, one can use an rna structure of interest to search an rna structure database, and build pssm iteratively to build an rna motif, as demonstrated for ire in this study .
rsmatch bears similarities to rna_align and rnaforester in that the structural particularities are either explicitly captured using hierarchical tree/forest structures or implicitly represented using arc-annotated structures.
however, rsmatch differs from rna_align and rnaforester in two major aspects.
first, rsmatch keeps structural consistence by only allowing single bases matched with single bases and base pairs matched with base pairs whereas rna_align and rnaforest do not impose this restriction.
second, rsmatch keeps the integrity of single-stranded regions by matching one with another, instead of breaking a single-stranded region into pieces and aligning them with different single-stranded regions.
in addition, rsmatch has less time and space complexities than the other two tools.
the concept of circles introduced in this paper is reminiscent of the "k-loop" described in the classic rna structure prediction paper  <cit> .
the difference is that the circles can reflect the inter-base-pair relationship by focusing on two base-pairs at a time while the "k-loop" cannot.
by organizing all circles into a hierarchy tree, we can capture the overall structural particularity.
it should also be pointed that there is a major difference between the hierarchy tree introduced here and the parse trees of scfg  <cit> .
the hierarchy tree is constructed from circles and aims to obtain the panorama of the secondary structure of rna at a higher level than that of the scfg parse tree, while detailed information is still available within each circle in the tree.
with the introduction of partial structures, this two-level structure modeling  allows us to develop an efficient algorithm that runs at time o as we have shown in the paper.
our program takes full advantage of structure prediction techniques.
it separates rna folding from structure alignment.
simultaneous rna folding and alignment is believed to be the optimal solution for both finding the right structure and locating homologous sub-structures of rnas  <cit> .
unfortunately, it is computationally prohibitive for even a moderate number of rnas.
some improvements have been proposed, but extensive computing time still makes them infeasible for database searches  <cit> .
by separating the process into two steps, we greatly enhanced the computing efficiency, making it possible to process a large-scale pre-folded rna structure database for homologous motifs.
however, a drawback of using pre-folded rnas is that the prediction tools may not produce correct rna structures, as observed in our experiments.
it is estimated that the rna folding programs solely based on thermodynamic properties of rna can correctly predict rna structures with about 70% of chance  <cit> .
secondly, higher complex structures, such as pseudoknots, cannot be predicted in most commonly used programs, including the vienna rna package used in this study.
a solution to removing the first drawback is to choose suboptimal structures in addition to the optimal one to increase the chance of obtaining correct structures.
it has been reported that using suboptimal structures whose thermodynamic free energies are within 2% of that of the optimal one can greatly improve the structure prediction of rna  <cit> .
in our ire experiments, we found that the predicted structure for nm_001098/1– <dig> did not exhibit the existence of an ire motif.
by relaxing the free energy range, we finally detected the ire motif from one suboptimal structure whose free energy was  <dig>  kcal/mol higher than the optimal one.
because of the computing efficiency of our program, an increase of the number of rna structures does not impose big burden on database searching .
the cost will be at the database building stage, which is however done only once.
the moving window approach we used to extract and fold subsequences was aimed to make the folding process more accurate and efficient.
this is because rna folding programs are known to have pronounced difficulties in correctly predicting large rna structures.
furthermore, predicting the structure for a long sequence takes much longer time than predicting structures for its subsequences.
another advantage of using the moving window method is that small motifs falling in the overlapped subsequences could be folded twice, increasing the chance of their being detected.
pattern-based tools, such as patsearch and rnamotif, use descriptions of an rna structure as queries to search a sequence database for similar structures.
this type of search does not take into consideration the context of a hit sequence, which could influence the structure of the sequence.
for example, as shown in our experimental results, patsearch can achieve a satisfactorily high specificity when the structure of a pattern is not flexible and its description is relatively precise, such as the hsl <dig> motif.
however, the sensitivity of patsearch is low with rigid pattern descriptions.
for relatively flexible structures, such as ires, the specificity of patsearch drops because it does not take into account the context in which a motif is located.
on the other hand, using folded rna structures, the proposed rsmatch tool overcomes these shortcomings with a high specificity, thus complementing the pattern-based tool.
however, as also shown in our experimental results, the error existing in folding an rna sequence  can lower the sensitivity of rsmatch.
we suspect that the inaccuracy introduced by rna folding could be a bottleneck for our technique in achieving a very high sensitivity.
our scoring matrices for single-stranded and double-stranded regions and the gap penalty assignment are very primitive in the sense that they are not based on any probabilistic model or learned from any training data set.
one interesting observation in our hsl <dig> experiment was that rsmatch did find most hsl <dig> sites correctly.
however, the scoring scheme seemed not acute enough to filter out many false positives.
part of the problem is that there are not enough motifs that can be used to construct optimal scoring matrices.
in fact, we also tested the matrices  proposed by klein and eddy, which were built upon small subunit ribosomal rnas.
we did not find any discernible difference in our hsl <dig> experiment, in which both matrices were used .
another related question is whether different types of rna, such as trna, rrna, and utrs, need their own scoring matrices.
it is conceivable that large highly structured rnas, such as rrna, may be able to tolerate more mutations than short rna motifs that occur in utr regions.
if so, using different scoring matrices for different types of rnas will be warranted.
furthermore, it is possible that the mutation rate is different for nucleotides in different regions of an rna motif.
therefore, pssm might be more suitable in these cases.
to this end, the iterative search function of rsmatch, which searches a database using pssm, can be applied.
motivated by the statistical methods of assessing results in sequence alignment  <cit> , we tried to develop scores of our database search with known probabilistic distributions.
the score distribution seemed close to be normal .
however since our scoring scheme is still at its preliminary stage and much is to be learned about the rna structure database presented in the paper, we only presented search results in terms of ranking.
more elaborate statistical assessment of the search results will be developed in the future.
availability and requirements
the rsmatch package has been implemented in java and perl and is freely available for academic use at  or .
authors' contributions
jl, jtw, and bt participate in the design of the algorithm.
jl developed the software.
jl and bt did the study with various rna structures.
jh tested the software and participated in the study of hsl <dig>  jl, jtw and bt wrote the manuscript.

 ribonucleic acid plays various roles in the cell. many functions of rna are attributable to their structural particularities. rna motifs have been extensively studied for noncoding rnas, such as transfer rna, ribosomal rna, small nuclear rna, small nucleolar rna, etc.. more recently, small interfering rna and microrna have been under intensive studies. less well characterized are the structures in the un-translated regions of messenger rnas. however, biochemical and genetic studies have demonstrated a myriad of functions associated with the utrs in mrna metabolism, including rna translocation, translation, and rna stability. rna structure determination via biochemical experiments is laborious and costly. predictive approaches are valuable in providing guide information for wet lab experiments. rna structure prediction is usually based on thermodynamics of rna folding or phylogenetic conservation of base-paired regions. the former uses thermodynamic properties of various rna local structures, such as base pair stacking, hairpin loop, and bulge, to derive thermodynamically favourable secondary structures. a dynamic programming algorithm is used to find optimal or suboptimal structures. the most well-known tools belonging to this group are mfold and rnafold in the vienna rna package. similar tools have been developed in recent years to predict higher order structures, such as pseudoknots. on the other hand, rna structure prediction using phylogenetic information infers rna structures based on covariation of base-paired nucleotides. it is generally believed that methods using phylogenetic information are more accurate. however, their performance critically depends on the high quality alignment of a large number of structurally related sequences. tools that align biosequences, such as fasta and blast, are valuable in identifying homologous regions, which can lead to the discovery of functional units, such as protein domains, dna cis elements, etc.. however, their success is more evident in the study of dnas and proteins than of rnas. this is mainly because the sequence similarity among dnas and proteins can usually faithfully reflect their functional relationship, whereas additional structure information is needed to study the functional conservation among rnas. therefore, it is necessary to take into account both structural and sequential information in comparing rna sequences. several tools are available that carry out rna alignment and folding at the same time. the pioneer work by sankoff involves simultaneous folding and aligning of two rna sequences, and has huge time and space complexity. foldalign improves the sankoff's method by scoring the structure solely based on the number of base pairs, instead of the stacking energies; and disallowing branch structures. dynalign reduces the time complexity by restricting the maximum distance allowed between aligned nucleotides in two structures. by taking into account local similarity, stem energy and covariations, perriquet et al. proposed carnac for pairwise folding of rna sequences. ji et al. applied a graph-theoretical approach, called comrna, to detect the common rna secondary structure motifs from a group of functionally or evolutionally related rna sequences. one noticeable advantage of comrna is its capability to detect pseudoknot structures. in addition, algorithms using derivative-free optimization techniques, such as genetic algorithms and simulated annealing, have been proposed to increase the accuracy in structure-based rna alignment. for example, notredame et al. presented raga to conduct alignment of two homologous rna sequences when the secondary structure of one of them was known. as shown in table most of these methods suffer from high time complexities, making the structure-based rna alignment tools much less efficient than sequence-based alignment tools. tools that search for optimal alignment for given structures include rnadistance, rna_align, and rnaforester. rnadistance uses a tree-based model to coarsely represent rna secondary structures, and compares rna structures based on edit distance. in a similar vein, rna_align models rna secondary structures by nested and/or crossing arcs that connect bonded nucleotides. with the crossing arcs, rna_align is able to align two rna secondary structures, one of which could contain pseudoknots. rnaforester extends the tree model to forest model, which significantly improves both time and space complexities. in addition, methods using stochastic context free grammars have been developed to compare two rna structures. original scfg models require a prior multiple sequence alignment for the training purpose, thus their applicability is limited to rna types for which structures of a large number of sequences are available, such as snorna and trna. however, rsearch and stemloc, both based on scfg, are capable of conducting pair-wise structure comparisons with no requirement for pre-alignment. rsearch uses ribosum substitution matrices derived from ribosomal rnas to score the matches in single-stranded and double-stranded regions. stemloc uses "fold envelope" to improve efficiency by confining the search space involved in calculations. the time and space complexities of these two tools are also listed in table furthermore, pattern-based techniques such as rnamotif, rnamot and patsearch have been used in database searches to detect similar rna substructures. these tools represent rna structures by a consensus pattern containing both sequence and structure information. one important advantage of these pattern-based tools is the ability of dealing with pseudoknots. we present here a computationally efficient tool, called rsmatch, capable of both globally and locally aligning two rna secondary structures. rsmatch does not require any prior knowledge of structures of interest. it can uncover structural similarities by means of direct aligning at the structure level. we demonstrate its application to database search and multiple alignment. we compared rsmatch with three widely used tools, patsearch, stemloc and rsearch, demonstrating that rsmatch is faster or achieves comparable or higher accuracy than the existing tools when applied to a number of known rna structures, including simple stem-loops and complex structures containing junctions. implementation secondary structure decomposition rsmatch models rnas by a structure decomposition scheme similar to the loop model commonly used in the algorithms for rna structure prediction. with this model, pseudoknots are not allowed. our method differs from the loop decomposition methods in that it completely decomposes an rna secondary structure into units called circles. when the secondary structure is depicted on a plane, a circle is defined as a set of nucleotides that are reachable from one another without crossing any base pair. as shown in figure 1a, all circles are closed or ended by a base pair except the first circle, which always contains the 5'-most and the 3'-most bases. various types of rna structures, such as bulge, loop, and junction can be represented by circles, as shown in figure 1a. circles of an rna structure can be organized as a hierarchical tree according to their relative positions in the secondary structure, where each tree node corresponds to a circle. this tree organization is informative to deduce the structural relationship among circles and reflects the structure particularities of the given rna secondary structure. if two circles reside on the same lineage in the tree, the circle appearing higher in the tree is called an ancestor of the other, and the latter is a descendent of the former. as a result, in the context of the hierarchical tree, two distinct circles fall into one of the following two categories, in the order of decreasing closeness: the two circles maintain an ancestor/descendent relationship, or they share a common ancestor in the tree. for example, in figure 1b, circle is an ancestor of circle whereas circle does not have ancestor/descendent relationship with circle since they are not on the same lineage. the double-stranded region or stem of a structure is decomposed into a set of "degenerated" circles, each containing only two base pairs. as such, a stem of n bases in length will result in n - consecutive degenerated circles. since a base pair may have two associated circles; we name one circle "the parent circle" and the other "the child circle" according to their positions in the hierarchical tree. for example, for the boxed c-g base pair in figure 1a, circle is its parent circle and circle is its child circle. structure alignment formalization given an rna secondary structure, we consider two types of structure components, single bases and base pairs, in the secondary structure. to integrate both sequence and structure information, we introduce two constraints among the structure components: precedence constraint and hierarchy constraint. the precedence constraint is defined based on the precedence order among structure components and the hierarchy constraint specifies the inter-component relationship in the context of the hierarchical tree described above. the precedence order is determined by the 3' bases of individual structure components: the one with its 3' base closer to the rna sequence's 5'-end precedes the other. for example, in figure 1a, the single base component u in circle precedes the base pair component c-g in circle to capture the inter-component relationship within the hierarchical tree context, we need to map each structure component to a circle in the tree. it is obvious that each single base can be mapped to a unique circle. however, a base pair could be mapped to two alternate circles: one parent circle and one child circle. to resolve this ambiguity, we always require mapping to the parent circle. the inter-component relationship is then reduced to the inter-circle relationship of three types: ancestor/descendent, common ancestor, and identical circle. given two rna secondary structures a and b, where a, referred to as the query structure, has m structure components {a a..., am} and b, referred to as the subject structure, has n structure components {b b..., bn}, the structure alignment between a and b is formalized as a conditioned optimization problem based on the above two constraints: given a scoring scheme consisting of two matrices, one for matching two single bases and the other for matching two base pairs, find an optimal alignment between the two sets of structure components such that the aforementioned precedence and hierarchy constraints are preserved for any two matched component pairs and. in other words, the two structure constraints between ai and aj must be respectively equivalent to that between bi and bj. this formalization has an implicit biological significance in that a single stranded region in one structure, if not aligned to a gap as a whole, will always align with a single stranded region in the other structure. this alignment requirement is important because single stranded regions are usually treated as functional units in binding to specific proteins. algorithmic framework a dynamic programming algorithm is employed in rsmatch. as with sequence alignment, the structure alignment could be either global or local. the difference lies only in the setup of initialization conditions; the algorithmic framework is the same since both global and local alignments must preserve the two constraints described above. a scoring table is established with its rows/columns corresponding to the structure components of the two given rna secondary structures. we organize the rows/columns in such a way that the precedence and hierarchy constraints are combined and easy to follow in the course of alignment computation. specifically, we sort the structure components of each structure according to the precedence order defined above. it is straightforward that this arrangement of rows/columns makes the precedence constraint automatically preserved. however, preservation of the hierarchy constraint is much more complicated and can only be accomplished in the derivative analysis for each cell in the scoring table. we will discuss the derivation when filling in the scoring table. each cell of the scoring table represents an intermediate comparison between two partial structures corresponding to the cell's row and column components respectively. the partial structure with respect to a structure component c is a set of structure components sc such that for any component a ∈ sc, the following three structure constraints between c and a must be satisfied: a precedes c; by the hierarchy constraint, a is not an ancestor of c; and c itself is included in sc. furthermore, since a base pair could appear in two circles, its corresponding partial structure could be divided into two smaller substructures: parent structure and child structure. formally, given a base pair component c, the parent structure of c is the set of structure components pc ⊆ sc such that for any component a ∈ pc, a's 3'-base is always 5' upstream of c's 5'-base; the child structure of c is the set of structure components lc ⊆ sc such that for any component a ∈ lc, a's 5'-base is always 3' downstream of c's 5'-base. it can be shown that pc ∪ lc = sc and pc ∩ lc = φ. examples of partial structures are given in figure 1c–1e. as shown in figure 1c, for a base pair, its child and parent structures together constitute the whole partial structure for the base pair. as we will see in the following discussions, the concept of a partial structure and its byproducts form the kernel of our algorithmic framework. we can solve the rna structure alignment problem progressively by aligning small structures and expanding each of them one structure component at a time until all structure components are covered. preliminaries cells in the scoring table are processed row by row from top to bottom and from left to right within each row. by considering the row/column components, we have three types of cells: a cell corresponding to two single bases; a cell corresponding to one single base and one base pair; and a cell corresponding to two base pairs. for, each cell stores the score of aligning the partial structures corresponding to the cell's row and column components respectively. for and, we need to consider alignments involving the partial and child structures induced by the base pair components. notice that the parent structures of the base pair components are excluded. it can be shown that each parent structure pc of component c can always be considered as the partial structure sx of some other component x, which means we only need to consider child and partial structures in the alignment computation. consequently, the above three types of cells have one, two and four alignment scores respectively. a scoring scheme is required to score the match of two structure components. we define the scoring scheme as a function g where a and b represent two structure components that are matched with each other. another important aspect of the alignment algorithm is to penalize the match involving gap. in the course of computation, one structure component could match with a gap or a whole small structure could match with a large gap. intuitively, the larger the gap is, the heavier the penalty will be. in our implementation, we set an atomic penalty value, denoted as u, for the smallest gap equivalent to a single base. the penalty value for a large gap is proportional to its size in terms of the number of bases matched with the gap. let a* be a small structure in the query rna structure a and b* a small structure in the subject rna structure b. the score obtained by aligning the two structures a* and b*, denoted as f, is, where g represents the total number of gaps in aligning a* and b*. initialization we assume that the row components are from the query rna structure a and the column components from the subject rna structure b. we focus on global alignment here; initializations for local alignment can be derived similarly. the initialization conditions deal with the cases where at least one of the structures under alignment is an empty structure φ. this is equivalent to setting up the 0th row/column in the scoring table. as discussed above, each base pair component has two small structures to be considered: a child structure and a partial structure. thus, the aforementioned three types of cells have one, two and four initialization scores respectively. for a given structure component x, let sx represent its partial structure. if x is a base pair, we also use lx to represent its child structure. we have f = furthermore, for any structure components a and b, f = |sa|·u, f = |sb|·u, if a and b are base pairs, f = |la|·u and f = |lb|·u where |·| represents the cardinality of the respective set. filling in the scoring table the simplest cell type is the one whose row component is a single base a. let ap denote the structure component that precedes a by precedence order established before. formally, in matching the partial structure sa with the partial structure sb there are only three possibilities: a is aligned with b; a is aligned with a gap; and b is aligned with a gap. thus the score of matching sa with sb can be calculated by equation. the second cell type is the one formed by one single base and one base pair. there are actually two symmetric subtypes where either a or b is a base pair. since the analysis is identical, we only focus on the former case where a is a base pair. as discussed before, besides the partial structure sa we have to consider the child structure la for the base pair a. thus, for this type of cells, we have to compute two alignment scores. by the principle of dynamic programming, the smaller size problem needs to be solved before the larger size problem. thus we first find the structure alignment between the child structure la and the partial structure sb. there are only two possibilities: the single base component b is aligned with a gap; and the base pair a is aligned with a gap. therefore we have in aligning the partial structure sa with the partial structure sb, to preserve precedence and hierarchy constraints simultaneously, there are only three possibilities: the single base b matches with a gap; the partial structure sb matches with the child structure la; the partial structure sb matches with the parent structure pa. thus, for the third cell type, a is a base pair and b is a base pair. we need to compute four alignment scores because each base pair corresponds to two structures: one child structure and one partial structure. while aligning the child structure la with the child structure lb, it is clear that since both a and b are the last components in the respective child structures by precedence order. equation gives the alignment score between the partial structure sa and the child structure lb. the first case corresponds to that b is aligned with a gap. if b does not match with a gap, it can be shown that, to preserve both precedence and hierarchy constraints, the second and third cases in equation cover all possible situations. similarly, we can calculate the score of aligning the child structure la and the partial structure sb as shown in equation. in aligning the partial structure sa with the partial structure sb, there are five possibilities: the parent structure pa is matched with the parent structure pb and the child structure la is matched with the child structure lb; the child structure la is matched with gaps; the child structure lb is matched with gaps; the parent structure pa is matched with gaps; and the parent structure pb is matched with gaps. therefore data sets all experiments were carried out on a linux system with two ghz intel processors and gb memory. a human utr structure database was constructed as follows. we downloaded human refseq mrna sequences from national center for biotechnology information. each refseq sequence containing utr regions, as indicated by refseq's genbank annotation, was processed to extract its 5'utr and 3'utr sequences. for each utr sequence, we took a nt subsequence at every 50th nucleotide position from 5' to 3', making consecutive subsequences overlap with one another on a nt segment. subsequences shorter than nt, e.g. at the 3' end, were also kept. using the vienna rna package's rnasubopt function with setting "-e 0", we then folded all obtained sequences to form the structure database. for any given rna sequence, the setting "-e 0" resulted in multiple rna structures all having the minimum free energy. the final database contained ~ rna secondary structures. the structural patterns of a histone 3'utr stem-loop structure and an iron responsive element were used in this study, based on their specifications in the utrdb database. three tools, patsearch, stemloc and rsearch, were employed for comparison purposes. the efficiency of these tools was measured by cpu running time. the performance of each program was assessed by specificity and sensitivity. specificity was calculated as tp/ and sensitivity as tp/, where tp was the number of true positives, fp the number of false positives, and fn the number of false negatives. to test the applicability of rsmatch to complex structures, we downloaded rna families from rfam. we only chose those families that had more than seed rnas and its consensus sequence length is no longer than nucleotides. we had families in the final data set. for each family, we randomly selected one member rna as the query rna and obtained its structure from rfam. we then randomly chose subject rnas in the same family. here we intentionally introduced noise by extending each subject rna sequence with its adjacent sequences at both 3' and 5' ends to make the total length three times its original one.  studies of stem-loop structures in utrs using our proposed algorithm, we first studied rna motifs in utr regions of human mrna sequences. a well-known fact is that the accuracy and efficiency of rna folding programs will decrease significantly when the sequences to be folded become very long. satisfactory performance is usually obtained when the sequences have moderate lengths, i.e. one hundred nucleotides. thus, we used a moving window scheme to get subsequences of nt and folded them using the vienna rna package. in the rsmatch package, this subsequence length is a user-defined parameter. since the nucleotide conservation in the single-stranded region of an rna sequence may differ from that in the double-stranded region, we used two scoring matrices, one for substitutions among single bases and the other among base pairs. this type of scoring scheme was also used in other studies. theoretically, the scoring matrix for single bases is a × table for all types of substitutions of single nucleotides, and the one for base pairs is a × table for all types of substitutions of base pairs. however since we used the vienna rna package, only six types of base pairs were observed in our studies, i.e. watson-crick base pairs a-u, u-a, g-c, c-g, and wobble base pairs g-u and u-g. values used in the two matrices were empirically chosen so as to conform to the general understanding of the sequence and structure conservation of rna motifs, as follows. mutations in the double-stranded region may not be detrimental to rna's function if the mutated sequence still preserves the same secondary structure. therefore base pair substitutions were rewarded with a positive score, instead of a penalty. a sequence in the single-stranded region may be important for rna's function, such as binding to proteins, and thus mismatches were penalized. to process gaps we used an arbitrary function u × l, where u was the atomic penalty value for a gap that is one single base long and l is the length of the gap in terms of the number of bases matched with the gap. in our experiments otherwise stated explicitly, the u was empirically set to - and changing the u value did not change the qualitative conclusion made in the paper provided that the absolute value of u was greater than any positive score in the scoring matrices. users can freely change the u value when applying rsmatch to their own data set. we tested our program with a query sequence containing an iron response element. the ire motif is a bipartite stem-loop structure containing ~ nucleotides. two alternative types of ires have been found, which differ in the middle region. type i has a bulge, whereas type ii has a small internal loop. ires have been found in both 5' and 3' utrs of genes that are involved in iron homeostasis in higher eukaryotic species. they interact with iron regulatory proteins and play a role in rna stability and translation. using a subsequence in the 3'utr of transferrin receptor that contains an ire motif, we searched the utr structure database described in implementation. a list of top hits is shown in figure the best hit of the search is the query structure itself, as expected. other regions of the same mrna and regions of other rnas are also found to have homologous structures with the query. as clearly shown in the result, the region containing the ire motif, which is from about the 30th nucleotide to about the 60th nucleotide of the query structure, has been located by the rsmatch program, indicating that a local optimal alignment has been achieved. among the top hits, several sequences are known to have ires, such as several regions in the 3'utr of transferrin receptor and the 5'utr of solute carrier family protein. other top hits have not been shown so far to have ires. it is not known if some of them are novel ire-containing rnas and the definitive answer will await wet lab validation. the output shows detailed alignment and related information, including the numbers of bases in the single-stranded and double-stranded regions, and the percentages of identity in single-stranded and double-stranded regions. rsmatch can also accept pattern-based rna structures to search a structure database. since a pattern-based search method has an intrinsic primitive scoring scheme by using degenerate bases, we used simplified binary matrices as the equivalent to score an alignment. in the matrices, the match of a pair of structure components, including those containing degenerated bases, was given a score of a mismatch was penalized by a score of - and the atomic gap penalty u was set to - to allow variability in single-stranded and/or double-stranded regions for a structure pattern, we introduced a wildcard "n " to represent optional single base component and base pair component. the meaning of "n" is identical to the iub code "n" except that the matching score for both structure components "n" and "n-n" is always zero regardless of whether they are aligned with a structure component or a gap. two rna motifs were used to test our method, namely a histone 3'utr stem-loop structure and ire. hsl which resides in the 3'utr region of histone mrnas, has a typical stem loop structure with two flanking tails. both the stem and the flanking sequences are important to bind with a stem-loop binding protein, which controls the pre-mrna processing and stability of histone mrnas. in contrast to the hsl motif, ire is relatively flexible in length and in nucleotide composition in its stem region. we compared our program with patsearch, a widely used tool that searches a sequence database for sequence and structure patterns. using the hsl motif and utr sequence database, patsearch found hits whose locations were presented in table among them, one is a false positive. therefore the specificity of patsearch is very high. this is attributable to the precise specification of the hsl pattern. however, if a pattern description is too precise, it may lead to the "overfitting" problem. this problem prevents the tool from finding slightly divergent structures, thus lowering the tool's sensitivity. indeed, several histone genes were not detected by patsearch, including two histone genes which were found by rsmatch among its top hits. several other histone genes appeared among the top hits of rsmatch. this indicates that by gaining specificity, patsearch loses sensitivity for hsl since rsmatch gives a score to each alignment, different cutoffs can be used for selecting top hits. it seems that newly detected true positives are heavily outnumbered by false positives as rsmatch relaxes its cutoff value. however, with some properly chosen cutoff, i.e. rsmatch could still achieve a comparable specificity with patsearch. one possible explanation of getting high false positives for rsmatch could be that, with respect to the particular case of the hsl motif, its secondary structure conformation might be too pervasive in rna sequences to be used as a discriminative feature. this could point out a problem concerning rsmatch's current scoring matrices, which need to be fine tuned to improve the tool's specificity. good tuning could be achieved by setting up the scoring matrices through learning from a training data set. one interesting observation, however, was that rsmatch and patsearch agreed perfectly upon the hsl locations for almost all of the true positives they found. using the ire motif, we performed further comparisons between rsmatch and three other tools: patsearch, stemloc and rsearch. we used default parameters for rsearch; for stemloc, the fold envelope was set to instead of using the large utr structure database described in implementation we constructed a small test data set to expedite the comparison process. first, we used patsearch to search human utr sequences for ire motifs. then for each hit sequence we selected its corresponding mrna's 3' or 5' utr sequence. following the same folding process as discussed in implementation, we folded these utr sequences to form the test data set. totally, patsearch found hits, among which were known true positives. therefore patsearch's specificity was ~33%. these hits were from distinct mrna sequences. we assumed that patsearch had a 100% sensitivity. we extracted the 5' and 3' utr sequences from the distinct mrnas and obtained utr sequences. we then folded the utr sequences to get a small test data set, which contained structures. using a known ire-containing structure, which was one of the true positives found by patsearch, as the query, we searched the small test data set. table shows the results we obtained. since rsearch accepts sequences only, it was tested using only the primary sequence information in the test data set. except for the ire-containing structure nm_ which was one of the true positives found by patsearch, and the query itself, all tools agreed on the ire locations for the other seven true positives without salient discrepancy. it was found that nm_ was not properly folded to exhibit the existence of the ire motif. rsmatch has the best specificity by ranking all seven true positives within its top hits with only one false positive. rsearch is close to rsmatch by ranking all seven true positives within its top hits with one false positive. in contrast, stemloc gives five false positives within its top hits. setting different cutoff values yields different specificity and sensitivity for each tool. the point of balanced specificity and sensitivity appears at the cutoff value of for all three tools. with this cutoff value, the specificity of rsmatch and rsearch tied at 7/ × 100% = %. this is better than the specificity of patsearch and the specificity of stemloc. the sensitivity of rsmatch, rsearch and stemloc is %, % and 50% respectively. it is worth noting that rsmatch runs ~30% faster than rsearch; it took rsearch seconds to search the whole data set of structures while rsmatch used only seconds. consequently, rsmatch would be suitable for analyzing large data sets. it should also be pointed out that rsmatch permits wildcards in database searching and structure matching, which are not supported by rsearch or stemloc. performance on complex structures we further tested how accurate rsmatch is for complex structures. to this end, we downloaded rna structures and sequences from the rfam database. we used rna structure families, each of which has more than seed sequences and has the consensus sequence length less than nucleotides. for each rna structure family, we randomly selected a structure and searched against randomly selected sequences belonging to the same family. to reflect real world scenarios, we extended rna sequences at both 5' and 3' ends so that the length of a subject sequence is three times that of the original one. to ensure that the folded structures are long enough to fully contain the structure being investigated, we required the moving window size to be times the length of the query rna sequence. furthermore, to include suboptimal structures, we used all structures with free energy within kcal/mol above the minimum one. compared with hsl and ire, the query structures we used were much more complex, with average length of ~ nt and more than 70% of them comprised of nested loops and conjunctions. to assess the accuracy, we used a measure called structure coverage, denoted as p, which is calculated by the following formula: p = |qalign|/max, where |qalign| and |salign| are the lengths of aligned portion of query rna and subject rna, respectively, and |q| is the length of query rna sequence. as shown in figure 5a, even though rsearch has slightly more points clustered around high coverage, the overall difference between rsmatch and rsearch is not significant. in addition, the difference between rsmatch and rsearch do not seem to be related to structure size or complexity. this result indicates that rsmatch has the ability to process complex structures. we also selected 5s rrna for further detailed tests. 5s rrna has a length of ~ nt, which contains several types of rna structures, including hairpin, internal loop, bulge, and junction. there are sequences in the 5s rrna family, allowing us to carry out a thorough analysis. we randomly chose one 5s rrna as query structure and ten others as subject sequences for alignment. this process was repeated times. the performance comparison of rsearch and rsmatch is shown in figure 5b. for 5s rrna, rsmatch outperforms rsearch in discovering the complete structure more frequently. an exemplary alignment is shown in figure 5c–5e. running efficiency by dynamic programming, the running time of computing an alignment equals the number of writing operations needed to fill the scoring table. thus the time complexity of rsmatch is o, where m is the number of structure components in the query rna structure. to test the scalability, we downloaded the seed sequences for 5s rrna family from rfam and randomly selected one annotated structure as the query while folding the rest sequences to prepare the structure database as discussed above. we plotted the rsmatch running time versus the database size. the program was run times and the result is shown in figure the nearly perfect linear growth of the running time gives an empirical proof that the algorithm's time complexity is bounded by o. multiple structure alignment and iterative database search we also extended rsmatch algorithm to conduct multiple structure alignment. an example using ire is shown in figure while the alignment algorithm is the same, the multiple alignment function uses a position-specific scoring matrix. for a given set of structures, the multiple alignment function first identifies the best alignment of two structures, and builds a pssm. the pssm is then used to search for the closest structure in the rest of the set. a flowchart of multiple structure alignment is shown in figure 7a. if the alignment score of a structure to the pssm is above a cutoff, it is selected and its structure is used to update the pssm. this step is iteratively conducted until no structures have alignment score above the cutoff. in a sense, this method employs an implicit guided hierarchical tree using the average value for joining nodes. as an example, from our human utr database we selected ire-containing structures and randomly chose other none-ire structures to form a small dataset and run rsmatch against it. the output is shown in figure 7b. the final result is in stockholm format for multiple structure alignment. conceivably, when the given set of structures is a large database, the multiple structure alignment function of rsmatch in effect conducts iterative search for finding similar structures. discussion and  the work presented here is intended to provide an efficient tool to directly perform structure alignment and search of rna secondary structure databases. its capability to carry out multiple structure alignment and iterative database search can potentially be used to uncover rna motifs ab initio. for example, one can use an rna structure of interest to search an rna structure database, and build pssm iteratively to build an rna motif, as demonstrated for ire in this study. rsmatch bears similarities to rna_align and rnaforester in that the structural particularities are either explicitly captured using hierarchical tree/forest structures or implicitly represented using arc-annotated structures. however, rsmatch differs from rna_align and rnaforester in two major aspects. first, rsmatch keeps structural consistence by only allowing single bases matched with single bases and base pairs matched with base pairs whereas rna_align and rnaforest do not impose this restriction. second, rsmatch keeps the integrity of single-stranded regions by matching one with another, instead of breaking a single-stranded region into pieces and aligning them with different single-stranded regions. in addition, rsmatch has less time and space complexities than the other two tools. the concept of circles introduced in this paper is reminiscent of the "k-loop" described in the classic rna structure prediction paper. the difference is that the circles can reflect the inter-base-pair relationship by focusing on two base-pairs at a time while the "k-loop" cannot. by organizing all circles into a hierarchy tree, we can capture the overall structural particularity. it should also be pointed that there is a major difference between the hierarchy tree introduced here and the parse trees of scfg. the hierarchy tree is constructed from circles and aims to obtain the panorama of the secondary structure of rna at a higher level than that of the scfg parse tree, while detailed information is still available within each circle in the tree. with the introduction of partial structures, this two-level structure modeling allows us to develop an efficient algorithm that runs at time o as we have shown in the paper. our program takes full advantage of structure prediction techniques. it separates rna folding from structure alignment. simultaneous rna folding and alignment is believed to be the optimal solution for both finding the right structure and locating homologous sub-structures of rnas. unfortunately, it is computationally prohibitive for even a moderate number of rnas. some improvements have been proposed, but extensive computing time still makes them infeasible for database searches. by separating the process into two steps, we greatly enhanced the computing efficiency, making it possible to process a large-scale pre-folded rna structure database for homologous motifs. however, a drawback of using pre-folded rnas is that the prediction tools may not produce correct rna structures, as observed in our experiments. it is estimated that the rna folding programs solely based on thermodynamic properties of rna can correctly predict rna structures with about 70% of chance. secondly, higher complex structures, such as pseudoknots, cannot be predicted in most commonly used programs, including the vienna rna package used in this study. a solution to removing the first drawback is to choose suboptimal structures in addition to the optimal one to increase the chance of obtaining correct structures. it has been reported that using suboptimal structures whose thermodynamic free energies are within 2% of that of the optimal one can greatly improve the structure prediction of rna. in our ire experiments, we found that the predicted structure for nm_001098/1– did not exhibit the existence of an ire motif. by relaxing the free energy range, we finally detected the ire motif from one suboptimal structure whose free energy was kcal/mol higher than the optimal one. because of the computing efficiency of our program, an increase of the number of rna structures does not impose big burden on database searching. the cost will be at the database building stage, which is however done only once. the moving window approach we used to extract and fold subsequences was aimed to make the folding process more accurate and efficient. this is because rna folding programs are known to have pronounced difficulties in correctly predicting large rna structures. furthermore, predicting the structure for a long sequence takes much longer time than predicting structures for its subsequences. another advantage of using the moving window method is that small motifs falling in the overlapped subsequences could be folded twice, increasing the chance of their being detected. pattern-based tools, such as patsearch and rnamotif, use descriptions of an rna structure as queries to search a sequence database for similar structures. this type of search does not take into consideration the context of a hit sequence, which could influence the structure of the sequence. for example, as shown in our experimental results, patsearch can achieve a satisfactorily high specificity when the structure of a pattern is not flexible and its description is relatively precise, such as the hsl motif. however, the sensitivity of patsearch is low with rigid pattern descriptions. for relatively flexible structures, such as ires, the specificity of patsearch drops because it does not take into account the context in which a motif is located. on the other hand, using folded rna structures, the proposed rsmatch tool overcomes these shortcomings with a high specificity, thus complementing the pattern-based tool. however, as also shown in our experimental results, the error existing in folding an rna sequence can lower the sensitivity of rsmatch. we suspect that the inaccuracy introduced by rna folding could be a bottleneck for our technique in achieving a very high sensitivity. our scoring matrices for single-stranded and double-stranded regions and the gap penalty assignment are very primitive in the sense that they are not based on any probabilistic model or learned from any training data set. one interesting observation in our hsl experiment was that rsmatch did find most hsl sites correctly. however, the scoring scheme seemed not acute enough to filter out many false positives. part of the problem is that there are not enough motifs that can be used to construct optimal scoring matrices. in fact, we also tested the matrices proposed by klein and eddy, which were built upon small subunit ribosomal rnas. we did not find any discernible difference in our hsl experiment, in which both matrices were used. another related question is whether different types of rna, such as trna, rrna, and utrs, need their own scoring matrices. it is conceivable that large highly structured rnas, such as rrna, may be able to tolerate more mutations than short rna motifs that occur in utr regions. if so, using different scoring matrices for different types of rnas will be warranted. furthermore, it is possible that the mutation rate is different for nucleotides in different regions of an rna motif. therefore, pssm might be more suitable in these cases. to this end, the iterative search function of rsmatch, which searches a database using pssm, can be applied. motivated by the statistical methods of assessing results in sequence alignment, we tried to develop scores of our database search with known probabilistic distributions. the score distribution seemed close to be normal. however since our scoring scheme is still at its preliminary stage and much is to be learned about the rna structure database presented in the paper, we only presented search results in terms of ranking. more elaborate statistical assessment of the search results will be developed in the future. availability and requirements the rsmatch package has been implemented in java and perl and is freely available for academic use at or. authors' contributions jl, jtw, and bt participate in the design of the algorithm. jl developed the software. jl and bt did the study with various rna structures. jh tested the software and participated in the study of hsl jl, jtw and bt wrote the manuscript. 
rsmatch can find the optimal global or local alignment between two rna secondary structures using two scoring matrices, one for single-stranded regions and the other for double-stranded regions. the accuracy of rsmatch is tested by experiments using a number of known rna structures, including simple stem-loops and complex structures containing junctions. this tool shall be useful to researchers interested in comparing rna structures obtained from wet lab experiments or rna folding programs, particularly when the size of the structure dataset is large. however, existing tools either require a large number of prealigned structures or suffer from high time complexities. alignment of rna secondary structures is important in studying functional rna motifs. in recent years, much progress has been made in rna motif finding and structure alignment. therefore it can be used to identify functional rna motifs. the time complexity of rsmatch is o where m is the size of the query structure and n that of the subject structure. when applied to searching a structure database, rsmatch can find similar rna substructures, and is capable of conducting multiple structure alignment and iterative database search. motivated by widely used algorithms for rna folding, we decompose an rna secondary structure into a set of atomic structure components that are further organized by a tree model to capture the structural particularities. this makes it difficult for the tools to process rnas whose prealigned structures are unavailable or process very large rna structure databases. we present here an efficient tool called rsmatch for aligning rna secondary structures and for motif detection. with respect to computing efficiency and accuracy, rsmatch compares favorably with other tools for rna structure alignment and motif detection.
 alignment of rna secondary structures is important in studying functional rna motifs. in recent years, much progress has been made in rna motif finding and structure alignment. however, existing tools either require a large number of prealigned structures or suffer from high time complexities. this makes it difficult for the tools to process rnas whose prealigned structures are unavailable or process very large rna structure databases.  we present here an efficient tool called rsmatch for aligning rna secondary structures and for motif detection. motivated by widely used algorithms for rna folding, we decompose an rna secondary structure into a set of atomic structure components that are further organized by a tree model to capture the structural particularities. rsmatch can find the optimal global or local alignment between two rna secondary structures using two scoring matrices, one for single-stranded regions and the other for double-stranded regions. the time complexity of rsmatch is o where m is the size of the query structure and n that of the subject structure. when applied to searching a structure database, rsmatch can find similar rna substructures, and is capable of conducting multiple structure alignment and iterative database search. therefore it can be used to identify functional rna motifs. the accuracy of rsmatch is tested by experiments using a number of known rna structures, including simple stem-loops and complex structures containing junctions.  with respect to computing efficiency and accuracy, rsmatch compares favorably with other tools for rna structure alignment and motif detection. this tool shall be useful to researchers interested in comparing rna structures obtained from wet lab experiments or rna folding programs, particularly when the size of the structure dataset is large. 