IEEE Transactions on Computational Biology and Bioinformatics (January-March)

78 Pages·2005·4.953 MB·English

IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005 1
Guest Editorial: WABI Special Section Part ll
Junhyong Kim and Inge Jonassen
(cid:1)
THE Fourth International Workshop on Algorithms in The second paper is “Topological Rearrangements and
BIoinformatics (WABI) 2004 was held in Bergen, Nor- LocalSearchMethodforTandemDuplicationTrees”andis
way, September2004.The programcommittee consisted of authoredbyDenisBertrandandOlivierGascuel.Thepaper
33membersandselected,among117submissions,39tobe approaches the problem of estimating the evolutionary
presentedattheworkshopandincludedintheproceedings history of tandem repeats. A tandem repeat is a stretch of
from the workshop (volume 3240 of Lecture Notes in DNA sequence that contains an element that is repeated
Bioinformatics, series edited by Sorin Istrail, Pavel Pevzner, multipletimesandwheretherepeatoccurrencesarenextto
and Michael Waterman). each other in the sequence. Since the repeats are subject to
The WABI 2004 program committee selected a small mutations,theyarenotidentical.Therefore,tandemrepeats
number of papers among the 39 to be invited to submit occur through evolution by “copying” (duplication) of
extendedversionsoftheirpaperstoaspecialsectionofthe repeat elements in blocks of varying size. Bertrand and
IEEE/ACM Transactions on Computational Biology and Bioin- Gascuel address the problem of finding the most likely
formatics. Four papers were published in the October- sequenceofeventsgivingrisetotheobservedsetofrepeats.
Each sequence of events can be described by a duplication
December 2004 issue of the journal and this issue contains
tree and one searches for the tree that is the most
anadditionalthreepapers.Wewouldliketothankboththe
parsimonious, i.e., one that explains how the sequence has
entire program committee for WABI and the reviewers of
evolved from an ancestral single copy with a minimum
the papers in this issue for their valuable contributions.
number of mutations along the branches of the tree. The
Thefirstofthepapersis“ANewDistanceforHighLevel
main difference with the standard phylogeny problem is
RNASecondaryStructureComparison”authoredbyJulien
that linear ordering of the tandem duplications impose
Allali and Marie-France Sagot. This paper describes algo-
constraints the possible binary tree form. This paper
rithmsforcomparingsecondarystructuresofRNAmolecules
describes a local search method that allows exploration of
wherethestructuresarerepresentedbytrees.Theproblemof
thecompletespaceofpossibleduplicationtreesandshows
classifyingRNAsecondarystructureisbecomingcriticalas that the method is superior to other existing methods for
biologists are discovering more and more noncoding func- reconstructing the tree and recovering its duplication
tionalelementsinthegenome(e.g.,miRNA).Mostlikely,the events.
major functional determinants of the elements are their The third paper is “Optimizing Multiple Seeds for
secondary structure and, therefore, a metric between such Homology Search” authored by Daniel G. Brown. The
secondary structures will also help delineate clusters of paper presents an approach to selecting starting points for
functional groups. In Allali and Sagot’s paper, two tree pairwise local alignments of protein sequences. The
representations of secondary structure are compared by problem of pairwise local alignment is to find a segment
analysing how one tree can be transformed into the other from each so that the two local segments can be aligned to
obtain a high score. For commonly used scoring schemes,
using an allowed set of operations. Each operation can be
this can be solved exactly using dynamic programming.
associatedwithacostandthedistancebetweentwotreescan
However,pairwisealignmentisfrequentlyappliedtolarge
then be defined as the minimum cost associated with a
datasetsandheuristicmethodsforrestrictingalignmentsto
transformofonetreetotheother.AllaliandSagotintroduce
be considered are frequently used, for instance, in the
two new operations that they name edge fusion and node
BLAST programs. The key is to restrict the number of
fusion and show that these alleviate limitations associated
alignments as much as possible, by choosing a few good
with the classical tree edit operations used for RNA
seeds,without missing highscoring alignments. The paper
comparison. Importantly, they also present algorithms for
shows that this can be formulated as an integer program-
calculating the distance between trees allowing the new
mingproblemandpresentsalgorithmforchoosingoptimal
operationsinadditiontotheclassicalones,andanalyzethe seeds. Analysis is presented showing that the approach
performanceofthealgorithms. givesfourtimesfewerfalsepositives(unnecessaryseeds)in
comparison with BLASTP without losing more good hits.
. J. Kim is with the Department of Biology, University of Pennsylvania,
3451WalnutStreet,Philadelphia,PA19104. Junhyong Kim
E-mail:[email protected]. Inge Jonassen
. I. Jonassen is with the Department of Informatics and Computational
Guest Editors
BiologyUnit,UniversityofBergen,HIBN5020Bergen,Norway.
E-mail:[email protected].
For information on obtaining reprints of this article, please send e-mail to:
[email protected].
1545-5963/05/$20.00(cid:1)2005IEEE PublishedbytheIEEECS,CI,andEMBSocieties&theACM
2 IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005
Junhyong Kim is the Edmund J. and Louise Inge Jonassen is a professor of computer
Kahn Term Endowed Professor in the Depart- scienceintheDepartmentofInformaticsatthe
mentofBiologyattheUniversityofPennsylvania. University of Bergen in Norway, where he is
HeholdsjointappointmentsintheDepartmentof memberofthebioinformaticsgroup.Heisalso
ComputerandInformationScience,PennCenter affiliated with the Bergen Center for Computa-
for Bioinformatics, and the Penn Genomics tional Science atthe sameuniversitywhere he
Institute. He serves on the editorial board of headstheComputationalBiologyUnit.Heisalso
Molecular Development and Evolution and the vicepresidentoftheSocietyforBioinformaticsin
IEEE/ACMTransactionsonComputationalBiol- theNordicCountries(SocBiN)andamemberof
ogyandBioinformatics,thecounciloftheSociety theboardoftheNordicBioinformaticsNetwork.
for Systematic Biology, and the executive committee of the Cyber Hecoordinatesthetechnologyplatformforbioinformaticsfundedbythe
Infrastructure for Phylogenetics Research. His research focuses on Norwegian Research Council functional genomics programme FUGE.
computational and experimental approaches to comparative develop- Hehasworkedinthefieldofbioinformaticssincetheearly1990s,where
ment. The current focus of his lab is in three areas: computational he has primarily focused on methods for discovery of patterns with
phylogenetics, in silico gene discovery, and comparative development applicationstobiologicalsequencesandstructuresandonmethodsfor
usinggenome-widegeneexpressiondata. theanalysisofmicroarraygeneexpressiondata.
. For more information on this or any other computing topic,
pleasevisitourDigitalLibraryatwww.computer.org/publications/dlib.
IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005 3
A New Distance for High Level RNA
Secondary Structure Comparison
Julien Allali and Marie-France Sagot
Abstract—WedescribeanalgorithmforcomparingtwoRNAsecondarystructurescodedintheformoftreesthatintroducestwonew
operations,callednodefusionandedgefusion,besidesthetreeeditoperationsofdeletion,insertion,andrelabelingclassicallyusedin
theliterature.Thisallowsustoaddresssomeseriouslimitationsofthemoretraditionaltreeeditoperationswhenthetreesrepresent
RNAsandwhatissearchedforisacommonstructuralcoreoftwoRNAs.Althoughthealgorithmcomplexityhasanexponentialterm,this
termdependsonlyonthenumberofsuccessivefusionsthatmaybeappliedtoasamenode,notonthetotalnumberoffusions.The
algorithmremainsthereforeefficientinpracticeandisusedforillustrativepurposesonribosomalaswellasonothertypesofRNAs.
IndexTerms—Treecomparison,editoperation,distance,RNA,secondarystructure.
(cid:1)
1 INTRODUCTION
RNASareoneofthefundamentalelementsofacell.Their 1. hairpinloopswhicharesequencesofunpairedbases
role in regulation has been recently shown to be far closing a helix;
more prominent than initially believed (20 December 2002 2. internal loops which are sequences of unpaired
issue of Science, which designated small RNAs with bases linking two different helices;
regulatory function as the scientific breakthrough of the 3. bulgeswhichareinternalloopswithunpairedbases
year). It is now known, for instance, that there is massive on one side only of a helix;
4. multiloopswhichareunpairedbaseslinkingatleast
transcriptionofnoncodingRNAs.Yetcurrentmathematical
three helices.
and computer tools remain mostly inadequate to identify,
analyze, and compare RNAs. Stems are successions of one or more among helices,
An RNA may be seen as a string over the alphabet of internal loops, and/or bulges.
nucleotides (also called bases), {A, C, G, T}. Inside a cell, The comparison of RNA secondary structures is one of
RNAsdonotretainalinearform,butinsteadfoldinspace. themainbasiccomputationalproblemsraisedbythestudy
Thefoldisgivenbythesetofnucleotidebasesthatpair.The of RNAs. It is the problem we address in this paper. The
maintypeofpairing,calledcanonical,correspondstobonds motivationsaremany.RNAstructurecomparisonhasbeen
of the type A(cid:1)U and G(cid:1)C. Other rarer types of bonds
used in at least one approach to RNA structure prediction
may be observed, the most frequent among them is G(cid:1)U,
that takes as initial data a set of unaligned sequences
also called the wobble pair. Fig. 1 shows the sequence of a
supposed to have a common structural core [1]. For each
foldedRNA.Eachboxrepresentsaconsecutivesequenceof
sequence, a set of structural predictions are made (for
bonded pairs, corresponding to a helix in 3D space. The
instance, all suboptimal structures predicted by an algo-
secondary structure of an RNA is the set of helices (or the
rithm like Zucker’s MFOLD [15], or all suboptimal sets of
list of paired bases) making up the RNA. Pseudoknots,
whichmaybedescribedasapairofinterleavedhelices,are compatiblehelicesorstems).Thecommonstructureisthen
in general excluded from the secondary structure of an found by comparing all the structures obtained from the
RNA.RNAsecondarystructurescanthusberepresentedas initial set of sequences, and identifying a substructure
planargraphs.AnRNAprimarystructureisitssequenceof commontoall,ortosomeofthesequences.RNAstructure
nucleotides while its tertiary structure corresponds to the comparison is also an essential element in the discovery of
geometric form the RNA adopts in space. RNA structural motifs, or profiles, or of more general
Apartfromhelices,theothermainstructuralelementsin
models that may then be used to search for other RNAs of
an RNA are:
the same type in newly sequenced genomes. For instance,
generalmodelsfortRNAsandintronsofgroupIhavebeen
derived by hand [3], [10]. It is an open question whether
. J. Allali is with the Institut Gaspard-Monge, Universite´ de Marne-la-
Valle´e,Cite´Descartes,Champs-sur-Marne,77454,Marne-la-Valle´eCedex models at least as accurate as these, or perhaps even more
2,France.E-mail:[email protected]. accurate,couldhavebeenderivedinanautomaticway.The
. M.-F.SagotiswithInriaRhoˆne-Alpes,Universite´ClaudeBernard,LyonI,
43BdduNovembre1918,69622Villeurbannecedex,France. identification of smaller structural motifs is an equally
E-mail:[email protected]. important topic that requires comparing structures.
Manuscriptreceived11Oct.2004;accepted20Dec.2004;publishedonline As we saw, the comparison of RNA structures may
30Mar.2005.
concernknownRNAstructures(thatis,structuresthatwere
For information on obtaining reprints of this article, please send e-mail to:
[email protected],andreferenceIEEECSLogNumberTCBB-0164-1004. experimentally determined) or predicted structures. The
1545-5963/05/$20.00(cid:1)2005IEEE PublishedbytheIEEECS,CI,andEMBSocieties&theACM
4 IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005
Fig.1.PrimaryandsecondarystructuresofatransferRNA.
objectiveinbothcasesisthesame:tofindthecommonparts
of such structures.
In[11],ShapirosuggestedtomathematicallymodelRNA Fig.2.Exampleofdifferenttreerepresentations((b),(c),(d),and(e))of
secondary structures without pseudoknots by means of thesameRNA(a).
trees. The trees are rooted and ordered, which means that
theorderamongthechildrenofanodematters.Thisorder nodeislabeledwith,respectively,abaseorapairofbases.A
corresponds to the 5’-3’ orientation of an RNA sequence. node of the tree in Fig. 2c represents a set of successive
GiventwotreesrepresentingeachanRNA,therearetwo unpairedbasesorofstackedpairedones.Thelabelofanode
main ways for comparing them. One is based on the isanintegerindicating,respectively,thenumberofunpaired
computation of the edit distance between the two trees basesortheheightofthestackofpairedones.Thenodesofthe
while the other consists in aligning the trees and using the tree in Fig. 2d represent elements of secondary structure:
scoreofthealignmentasameasureofthedistancebetween hairpinloop(H),bulge(B),internalloop(I),ormultiloop(M).
the trees. Contrary to what happens with sequences, the The edges correspond to helices. Finally, the tree in Fig. 2e
two, alignment and edit distance, are not equivalent. The contains only the information concerning the skeleton of
alignment distance is arestrained form of the edit distance multiloopsofanRNA.Thelastrepresentation,thoughgiving
between two trees, where all insertions must be performed ahighlysimplifiedviewofanRNA,isimportantnevertheless
before any deletions. The alignment distance for general as it is generally accepted that it is this skeleton which is
treeswasdefinedin1994byJiangetal.in[9]andextended usually the most constrained part of an RNA. The last two
to an alignment distance between forests in [6]. More models may be enriched with information concerning, for
recently, Ho¨chsmann et al. [7] applied the tree alignment instance,thenumberof(unpaired)basesinaloop(hairpin,
distance to the comparison of two RNA secondary
internal,multi)orbulge,andthenumberofpairedbasesina
structures. Because of the restriction on the way edit
helix.Thefirstlabelthenodesofthetree,theseconditsedges.
operations can be applied in an alignment, we are not
Other types of information may be added (such as overall
concerned in this paper with tree alignment distance and
compositionoftheelementsofsecondarystructure).Infact,
wethereforeaddressexclusivelyfromnowontheproblem
one could consider working with various representations
of tree edit distance.
simultaneouslyorinaninterlocked,multilevelfashion.This
OurwayforcomparingtwoRNAsecondarystructuresis
goesbeyondthescopeofthispaperwhichisconcernedwith
thentoapplyanumberoftreeeditoperationsinoneorbothof
comparingRNAsecondarystructuresusinganyoneamong
the trees representing the RNAs until isomorphic trees are
the many tree representations possible. We shall, however,
obtained. The currently most popular program using this
commentfurtheronthismultilevelapproachlateron.
approachisprobablytheViennapackage[5],[4].Thetreeedit
Concerningtheobjectivesofthispaper,theyaretwofold.
operations considered are derived from the operations
Thefirstistogivesomeindicationsonwhytheclassicaledit
classically applied to sequences [13]: substitution, deletion,
operationsthathavebeenconsideredsofarintheliterature
andinsertion.In1989,ZhangandShasha[14]gaveadynamic
forcomparingtreespresentsomelimitationswhenthetrees
programming algorithm for comparing two trees. Shapiro
standforRNAstructures.Threecasesofsuchlimitationswill
and Zhang then showed [12] how to use tree editing to
beillustratedthroughexamplesinSection3.InSection4,we
compareRNAs.Thelatteralsoproposedvarioustreemodels
that could be used for representing RNA secondary struc- then introduce two novel operations, so-called node-fusion
tures.Eachsuggestedtreeoffersamoreorlessdetailedview and edge-fusion, that enable us to address some of these
of an RNA structure. Figs. 2b, 2c, 2d, and 2e present a few limitationsandthengiveadynamicprogrammingalgorithm
examplesofsuchpossibleviewsfortheRNAgiveninFig.2a. forcomparingtwoRNAstructureswiththesetwoadditional
In Fig. 2, the nodes of the tree in Fig. 2b represent either operations. Implementation issues and initial results are
unpairedbases(leaves)orpairedbases(internalnodes).Each presentedinSection4.InSection5,wegiveafirstapplication
ALLALIANDSAGOT: ANEWDISTANCEFORHIGHLEVELRNASECONDARYSTRUCTURECOMPARISON 5
Letaninsertionoradeletioncostoneandtherelabelingof
anodecostzeroifthelabelisthesameandoneotherwise.For
thetwotreesofthefigureontheleft,theseriesrelabelðA!
FÞ:deleteðBÞ:insertðGÞrealizestheeditingofthelefttreeinto
the right one and costs 3. Another possibility is the series
deleteðBÞ:relabelðA!GÞ:insertðFÞ which also costs 3. The
distancebetweenthesetwotreesis3.
Fig.3.Editoperations:(a)theoriginaltreeT,(b)deletionofthenode
labelled D, (c) insertion of the node labeled I, and (d) relabeling of a
nodeinT (thelabelAoftherootischangedintoK).
of our algorithm to the comparison of two RNA secondary
structures. Finally, in Section 6, we sketch the main ideas
behindthemultilevelRNAcomparisonapproachmentioned Given a series of operations S, let us consider the nodes
above.Beforethat,westartbyintroducingsomenotationand of T that are not deleted (in the initial tree or after some
byrecallinginthenextsectionthebasicsaboutclassicaltree relabeling).SuchnodesareassociatedwithnodesofT0.The
editoperationsandtreemapping. mapping MS relative to S is the set of couples ðu;u0Þ with
Thispaperisanextendedversionofapaperpresentedat u2T and u0 2T0 such that u is associated with u0 by S.
the Workshop on Algorithms in BioInformatics (WABI) in Theoperationsdescribedabovearethe“classicaltreeedit
2004,inBergen,Norway.Afewmoreexamplesaregivento operations”thathavebeencommonlyusedintheliterature
illustrate some of the points made in the WABI paper, forRNAsecondarystructurecomparison.Wenowpresenta
complexity and implementation issues are discussed in fewresultsobtainedusingsuchclassicaloperationsthatwill
more depth as are the cost functions and a multilevel allowustoillustrateafewlimitationstheymaypresentwhen
usedforcomparingRNAstructures.
approach to comparing RNAs.
2 TREE EDITING AND MAPPING 3 LIMITATIONS OF CLASSICAL TREE EDIT
OPERATIONS FOR RNA COMPARISON
Let T be an ordered rooted tree, that is, a tree where the
Assuggestedin[12],thetreeeditoperationsrecalledinthe
order among the children of a node matters. We define
three kinds of operations on T: deletion, insertion, and previous section can be used on any type of tree coding of
relabeling (corresponding to a substitution in sequence an RNA secondary structure.
Fig.4showstwoRNAsePsextractedfromthedatabase[2]
comparison). The operations are shown in Fig. 3. The
(they are found, respectively, in Streptococcus gordonii and
deletion (Fig. 3b) of a node u removes u from the tree. The
Thermotoga maritima). For the example we discuss now, we
childrenofubecomethechildrenofu’sfather.Aninsertion
code the RNAs using the tree representation indicated in
(Fig.3c)isthesymmetricofadeletion.Givenanodeu,we
Fig. 2b where a node represents a base pair and a leaf an
remove a consecutive (in relation to the order among the
unpaired base. After applying a few edit operations to the
children) set u1;...;up of its children, create a new node v,
trees,weobtaintheresultindicatedinFig.4,withdeleted/
makevachildofubyattachingitattheplacewheretheset
insertedbasesingray.Wehavesurroundedafewregionsthat
was,and,finally,makethesetu1;...;up (inthesameorder) match in the two trees. Bases in the rectangular box at the
thechildrenofv.Therelabelingofanode(Fig.3d)consists
bottomoftheRNAontheleftarethusassociatedwithbasesin
simply in changing its label.
thebottomrightmostrectangularboxoftheRNAontheright.
GiventwotreesT andT0,wedefineS ¼fs1...segtobe Thesameisobservedforthebasesintheovalboxesforboth
a series of edit operations such that, if we apply succes-
RNAs.Suchmatchesillustrateoneofthemainproblemswith
sivelytheoperationsinS tothetreeT,weobtainT0 (i.e.,T
theclassicaltreeeditoperations:BasesinoneRNAmaybe
and T0 become isomorphic). A series of operations like S
mapped to identically labeled bases in the other RNA to
realizes the editing of T into T0 and is denoted by T !S T0.
minimise the total cost, while such bases should not be
We define a function cost from the set of possible edit
associatedintermsoftheelementsofsecondarystructureto
operations(deletion,insertion,relabeling)totheintegers(or which they belong. In fact, such elements are often distant
thereals)suchthatcosts isthescoreoftheeditoperations. fromoneanotheralongthecommonRNAstructure.Wecall
IfSisaseriesofeditoperations,wedefinebyextensionthat this problem the “scattering effect.” It is related to the
P
costS is s2Scosts.Wecandefinetheeditdistancebetween definitionoftreeeditoperations.Inthecaseofthisexample
two trees as the series of operations that performs the andoftherepresentationadopted,theproblemmighthave
editing of T into T0 and such that its cost is minimal: been avoided if structural information had been used.
distanceðT;T0Þ¼fminðcost ÞjT !S T0g. Indeed, the problem appears also because the structural
S
6 IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005
Fig.4.Illustrationofthescatteringeffectproblem.ResultofthematchingoftwoRNAsePs,ofStreptococcusgorgoniiandofThermotogamaritima,
usingthemodelgiveninFig.2b.
location of an unpaired base is not taken into account. It is followedbyaninternalloopandanotherhelixofsize5.By
therefore possible to match, for instance, an unpaired base definition (see Section 2), the algorithm can only associate
fromahairpinloopwithanunpairedbasefromamultiloop. one element in the first tree to one element in the second
Usinganothertypeofrepresentation,asweshalldo,would, tree.Inthiscase,wewouldliketoassociatethehelixofthe
however,notbeenoughtosolveallproblemsasweseenext. left tree to the two helices of the second tree since it seems
Indeed,tocomparethesametwoRNAs,wecanalsousea clear that the internal loop represents either an inserted
more abstract tree representation such as the one given in element in the second RNA, or the unbonding of one base
Fig.2d.Inthiscase,theinternalnodesrepresentamultiloop, pair. This, however, is not possible with classical edit
internal-loop,orbulge,theleavescodeforhairpinloopsand operations.
edgesforhelices.TheresultoftheeditionofTintoT0forsome A third type of problem one can meet when using only
costfunctionispresentedinFig.5(weshallcomebacklaterto thethreeclassicaleditoperationstocomparetreesstanding
thecostfunctionsusedinthecaseofsuchmoreabstractRNA for RNAs is similar to the previous one, but concerns this
representations;forthesakeofthisexample,wemayassume time a node instead of edges in the same tree representa-
anarbitraryoneisused). tion.Often,anRNAmaypresentaverysmallhelixbetween
The problem we wish to illustrate in this case is shown two elements (multiloop, internal-loop, bulge, or hairpin-
bytheboxesinthefigure.Considertheboxesatthebottom. loop) while such helix is absent in the other RNA. In this
IntheleftRNA,wehaveahelixmadeupof13basepairs.In case, we would therefore have liked to be able to associate
the right RNA, the helix is formed by seven base pairs one node in a tree representing an RNA with two or more
Fig.5.Illustrationoftheone-to-oneassociationproblemwithedges.ResultofthematchingofthetwoRNAsePs,ofSaccharomycesuvarumandof
Saccharomyceskluveri,usingthemodelgiveninFig.2d.
ALLALIANDSAGOT: ANEWDISTANCEFORHIGHLEVELRNASECONDARYSTRUCTURECOMPARISON 7
Fig. 6. Illustration of the one-to-one association problem with nodes. The two RNAs used here are RNAsePs from Pyrococcus furiosus and
Metallosphaerasedula.Trianglesstandforbulges,diamondstandforinternalloops,andsquaresforhairpinloops.
nodesinthetreefortheotherRNA.Onceagain,thisisnot replacingeciandeuwithanewsingleedgee.Theedgeelinks
possible with any of the classical tree edit operations. An thefatherofutoci.Itslabelthenbecomesafunctionofthe
illustration of this problem is shown in Fig. 6. (numerical)labelsofeu,uandeci.Forinstance,ifsuchlabels
WeshalluseRNArepresentationsthattaketheelements indicatedthesizeofeachelement(e.g.,forahelix,thenumber
ofthestructureofanRNAintoaccounttoavoidsomeofthe ofitsstackedpairs,andforaloop,themin,maxortheaverage
scattering effect. Furthermore, in addition to considering ofitsunpairedbasesoneachsideoftheloop),thelabelofe
information of a structural nature, labels are attached, in could be the sum of the sizes of eu, u and eci. Observe that
general,tobothnodesandedgesofthetreerepresentingan mergingtwoedgesimpliesdeletingallsubtreesrootedatthe
RNA. Such labels are numerical values (integers or reals). childrencjofuforjdifferentfromi.Thecostofsuchdeletions
They represent in most cases the size of the corresponding isaddedtothecostoftheedgefusion.
element, but may also further indicate its composition, etc. AnexampleofnodefusionisgiveninFig.7b.Letubea
Such additional information is then incorporated into the nodeandc oneofitschildren.Performinganodefusionof
i
costfunctionsforallthreeeditoperations.Itisimportantto uandci consistsinmakinguthefatherofallchildrenofci
observe that when dealing with trees labeled at both the and in relabeling u with a value that is a function of the
nodesandedges,anynodeandtheedgethatleadstoit(or, values of the labels of u, ci and of the edge between them.
in an alternative perspective, departs from it) represent a Observe that a node fusion may be simulated using the
single object from the point of view of computing an edit classical edit operations by a deletion followed by a
distance between the trees. relabeling. However, the difference between a node fusion
It remains now to deal with the last two problems that andadeletion/relabelingisinthecostassociatedwithboth
are a consequence of the one-to-one associations between operations. We shall come back to this point later.
nodes and edges enforced by the classical tree edit Obviously, likeinsertions ordeletions, edge fusionsand
operations. To that purpose, we introduce two novel tree nodefusionshaveofcoursesymmetriccounterparts,which
edit operations, called the edge fusion and the node fusion. are the edge split and the node split.
Given two rooted, ordered, and labeled trees T and T0,
wedefinethe“editdistancewithfusion”betweenT andT0
4 INTRODUCING NOVEL TREE EDIT OPERATIONS
4.1 Edge Fusion and Node Fusion
Inordertoaddresssomeofthelimitationsoftheclassicaltree
editoperationsthatwereillustratedintheprevioussection,
weneedtointroducetwonoveloperations.Thesearetheedge
fusionandthenodefusion.Theymaybeappliedtoanyofthe
treerepresentationsgiveninFigs.2c,2d,and2e.
AnexampleofedgefusionisshowninFig.7a.Leteubean
edge leading to a node u, ci a child of u and eci the edge
between u and ci. The edge fusion of eu and eci consists in Fig.7.(a)Anexampleofedgefusion.(b)Anexampleofnodefusion.
8 IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005
Fig.8.ZhangandSasha’sdynamicprogrammingalgorithm:thetreedistancepart.Therightboxcorrespondstotheadditionaloperationsaddedto
takefusionintoaccount.
asdistancefusionðT;T0Þ¼fminðcostSÞjT !S T0gwithcoststhe index of the leftmost child of the subtree rooted at ti. Let
cost associated to each of the seven edit operations now Tði...jÞ denote the forest composed by the nodes ti...tj
considered (relabeling, insertion, deletion, node fusion and
(T (cid:3)Tð0...jTjÞÞ.Tosimplifynotation,fromnowon,when
split, edge fusion and split).
thereisnoambiguity,iwillrefertothenodeti.Inthiscase,
Proposition 1. If the following is verified:
distanceði1...i2;j1...j2ÞwillbeequivalenttodistanceðTði1
. costmatchða;bÞ is a distance, ...i2Þ;T0ðj1...j2ÞÞ.
. cost ðaÞ¼cost ðaÞ(cid:2)0, The algorithm of Zhang and Sasha is fully described by
ins del
. costnodefusionða;b;cÞ¼costnodesplitða;b;cÞ(cid:2)0, and the following recurrence formula:
. cost ða;b;cÞ¼cost ða;b;cÞ(cid:2)0,
edgefusion edgesplit
then distance is indeed a distance.
fusion
Proof.Thepositivenessofdistancefusion isgivenbythefact
if ðði ¼¼lði ÞÞandðj ¼¼lðj ÞÞÞ
that all elementary cost functions are positive. Its 1 2 1 2
symmetry is guaranteed by the symmetry in the costs MIN
of the insertion/deletion and (node/edge) fusion/split 8
distanceð i ...i (cid:1)1 ;j ...j Þþcost ði Þ
operations. Finally, it is straighforward to see that >< 1 2 1 2 del 2
distanceð i ...i ;j ...j (cid:1)1 Þþcost ðj Þ
distancefusion satisfies triangular inequality. tu >:distanceð i1...i2(cid:1)1 ;j1...j2(cid:1)1 Þþcostins 2ði ;j Þ
Besidestheabovepropertiesthatmustbesatisfiedbythe 1 2 1 2 match 2 2
cost functions in order to obtain a distance, others may be ð1Þ
introducedforspecificpurposes.Somewillbediscussedin
Section 5. else
We now present an algorithm to compute the tree edit
MIN
distance between two trees using the classical tree edit
4op.2eratAiolgnosrpitluhsmthe two operations just introduced. 8>>>>>>>><þddiisscttoaasnntdcceeelððði2Þ ii1......ii2Þ(cid:1)1 ;;jj1......jj2Þ(cid:1)1 ÞÞ ð2Þ
1 2 1 2
The method we introduce is a dynamic programming
aTlhgeoirritahlgmorbiathsemdiosndtihveidoendeipnrtowpoospeadrbtsy:TZhheayngfirasntdcoSmhapsuhtae. >>>>>>>>:þdisctoasnticnesððj2Þ i1...lði2Þ(cid:1)1 ;j1...lðj2Þ(cid:1)1 Þ
theeditdistancebetweentwotrees(thispartisdenotedby þdistanceð lði2Þ...i2 ;lðj2Þ...j2 Þ
TDist)andthenthedistancebetweentwoforests(thispart
Part(1)oftheformulacorrespondstoFig.8,whilepart(2)
is denotedbyFDist).Fig.8illustrates inpictorial form the
part TDist and Fig. 9 the FDist part of the computation. corresponds to Fig. 9. In practice, the algorithm stores in a
Inordertotakeourtwonewoperationsintoaccount,we matrixthescorebetweeneachsubtreeofT andT0.Thespace
need to compute a few more things in the TDist part. complexityisthereforeOðjTj(cid:4)jT0jÞ.Toreachthiscomplexity,
Indeed, we mustaddthe possibilityfor eachtree to havea
the computation must be done in a certain order (see
nodefusion(inversely,nodesplit)betweentherootandone
Section4.3).Thetimecomplexityofthealgorithmis
of its children, or to have an edge fusion (inversely edge
split) between the root and one of its children. These OðjTj(cid:4)minðleafðTÞ;heightðTÞÞ
additionaloperationsareindicatedintherightboxofFig.8. (cid:4)jT0j(cid:4)minðleafðT0Þ;heightðT0ÞÞÞ;
Wepresentnowaformaldescriptionofthealgorithm.Let
T beanorderedrootedtreewithjTjnodes.Wedenotebyti where leafðTÞ and heightðTÞ represent, respectively, the
the ith node in a postfix order. For each node ti, lðiÞ is the number of leaves and the height of a tree T.
ALLALIANDSAGOT: ANEWDISTANCEFORHIGHLEVELRNASECONDARYSTRUCTURECOMPARISON 9
Fig.9.ZhangandSasha’sdynamicprogrammingalgorithm:theforestdistancepart.
Theformulatocomputetheeditscoreallowingforboth leading to, respectively, nodes u and v. The symmetric
node and edge fusions follows. operationsaredenotedby,respectively,node splitðu;vÞand
edge splitðu;vÞ.
The distance computation takes two new parameters
if ðði1 (cid:2)lðikÞÞ and ðj1 (cid:2)lðjk0ÞÞÞ path and path0. These are sets of pairs ðeoru;vÞ which
indicate, for node ik (respectively, jk), the series of fusions
MIN that were done. Thus, a pair ðe;vÞ indicates that an edge
8
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><dfddiiiosssrtttaaaedpdnnnaaiiccccsstheeetthaaððð0þþþcfffnnÞhiiiccdcc;111iooeei;l...sdððss;...ffttt;...aineii;iiicd1lonkkkgÞd.o.(cid:1)(cid:1)egce.f.e11;f.f.ggðpuiuiif;;sacskc(cid:1)ii;;i(cid:1)tooci1h;;n1nþn;ffð;;ð1iijjiiffkc.cc11jg;þi;.1..1i;i.1k..;k.pi.Þ...Þk.a.jjð..þ(cid:1)t.okk.jh00i;(cid:1)bkgdk:10ið1;(cid:1)gsik;egps;:1g;;;tagp;:;ia;i;tac;ksnhÞ;Þt;e;h0Þcdþ;Þtfe:aþÞðþðjictuf1laoc;i¼c.so1iao.tsc..srmÞlt.jðet;iadni.kftcec0scigjlÞhðlhð(cid:1)1;jðia1pik.knk0g.aÞÞ;;g.tjjhekkd000ÞÞÞg; cfi4Tcncthanouuao.oahnf3slmrdseoscTioroouerehpdmmnblnpelIeaemeatrahetlnxeenviytanpivonihdtoshgtiyllioaniegabntpteuosmhehikdiefose-esbeelrgnnteedfhenestvoephdpeetwneraacueifmlattturtismihltifidosgemhoo:iceindðomosnoerpirtngk;nialeg.atvieecndhnÞamdnaacmvtbindeenewered.nnbdeCtTitwaletitaahhcohsttteewatpiemiisedntoloeenionprntesrte.denliodqtekcehWtuetuxaotawitihirnkatertr.oeeyidetnlshsnntosvrepcoetd,eehtwaoeweaypsphtfieaa,onidirktllrigailhenmsotacaiwftuoouotnhtnflirsedænesnðpðgucetatiat;;ffoihhoosvvrneeesrrÞÞ
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:foreddpaaiicssthtthaa0ððþþþþðcnn:ooohðdccdccubbbiooeeiil;sssssdððssj:::ffttttcaajnei0i:::Þdc11onn0gÞjjid..ekcckkeo..00eefsds..ððddppaii;;llaajikkitt;;ktttggðað0;;aaj;;j;;icappcnaa0fj0r;aa;crrjje0jfttþ1eekhhkjc01.0ccÞ1;;Þh..hh;ff.a..aajjj.n.1lnnl0j0.g(cid:1)..kgg;e.0.1ee(cid:1)jd.g.kdd1jÞj;0;ÞÞgcc;;00;(cid:1);ÞÞj1s;ke0j;tcp0jþal10t¼h..0:.lððjejk;0c;j0Þc0ÞÞ oiSciattTFtinnshohonufuðefdeendlcrlLtðiutthchtcichehehsle2oeaetdeeprÞnemtanrueeat.cvmrt.dccpdseeae.oetuioeiscihismol2rntoesteÞaneeafntpwtdf,o(nasiiu2otoicwcnieht)tdrinerdtaedeehdaeetneornbieTifhtotrsfeid.ht0tnrahtðatreswvilh(oentonð1eeojoecor)l2feteeabtdÞwttnossoi.ipttbcehs.hcraetarlo.ieaoaivcclnjbnnrhba2eytercÞucdldayoe.eisiottmmnbsWoewfobtya,puea(oewn2pZtlnmtwe)mcehsh.xtearauheaeWuifylteobnengsybnrtgeothrecmreawotttoamihwwevftnmehesaetdoetrmhypheeoacSefunfeoofonhottbmrhaaroaTetisetasmpwlpsieghtotemruaosopanvrintsr.aoeereidTtfmretttirrtthifhooahetmeoTmhennemrary0eesss-tt..
ð3Þ distance between two subtrees U and U0, the distance
between any subtrees of U and U0 must already have been
else setil ¼lðikÞandjl0 ¼lðjk0Þ measured. This explains the space complexity which is in
OðjTj(cid:4)jT0jÞandcorrespondstothesizeofthetableusedfor
MIN
storing such distances in memory.
8distanceðfi ...i g;;;fj ...j g;path0Þþdelði Þ
>>><distanceðfi11...ikk(cid:1)g1;path;f1j1...kj0k0(cid:1)1g;;Þþinsðjkk0Þ ð4Þ calcIfuwlateelosoekpaarta(t1e)lynotwhe, wdeisstaenectehabteittwiseennottnheecessusabrtryeetos
>>>:distaþncdeiðsftia1n.c.e.ðifl(cid:1)i1g.;.;.i;fgj;1p.a.t.hj;l0(cid:1)fj1g;.;.Þ.j g;path0Þ risoootendtahtei0paantdh jf0riofmi0 ilsðjoÞntothej,pafothr firoamndlðiÞj tnooidaensdojf,0
l k l0 k0 respectively, T and T0.
Given two nodes u and v such that v is a child of u, We define a set LRðTÞ of the left roots of T as follows:
node fusionðu;vÞ is the fusion of node v with u, and
LRðTÞ¼fkj1(cid:5)k(cid:5)jTjand 69k0 >ksuchthatlðk0Þ¼lðkÞg
edge fusionðu;vÞ is the edge fusion between the edges
10 IEEE/ACMTRANSACTIONSONCOMPUTATIONALBIOLOGYANDBIOINFORMATICS, VOL.2, NO.1, JANUARY-MARCH2005
Thealgorithmforcomputingtheeditdistancebetweent our algorithm is thus in Oðð2dÞ‘(cid:4)ð2d0Þ‘(cid:4)jTj(cid:4)jT0jÞ, with d
and T0 consists then in computing the distance between and d0 the maximum degrees of, respectively, T and T0.
each subtree rooted at a node in LRðTÞ and each subtree Thecomputationofthetimecomplexityofouralgorithm
rooted at a node in LRðT0Þ. Such subtrees are considered is done in a similar way as for the algorithm of Zhang and
fromtheleavestotherootofT andT0,thatis,intheorder Shasha. For each node of T and T0, one must compute the
of their indexes. number of subtree distance computations the node will be
Zhang and Shasha proved that this algorithm has a involved in by considering all subtrees rooted in, respec-
time complexity in OðjTj(cid:4)minðleafðTÞ; heightðTÞÞ(cid:4)jT0j(cid:4) tively, a node of LRðTÞ and a node of LRðT0Þ. In our case,
minðleafðT0Þ; heightðT0ÞÞÞ, leafðTÞ designating the num-
onemustalsotakeintoaccountforeachnodethepossibility
ber of leaves of T and heightðTÞ its height. In the worst
of applying a fusion. This leads to a time complexity in
case (fan tree), the complexity is in OðjTj2(cid:4)jT0j2Þ.
Taking fusion and split operations into account does Oðð2dÞ‘(cid:4)jTj(cid:4)minðleafðTÞ;heightðTÞÞ(cid:4)ð2d0Þ‘(cid:4)jT0j(cid:4)
not change the above reasoning. However, we must now minðleafðT0Þ;heightðT0ÞÞÞ:
store in memory the distance between all subtrees
Tðlði2Þ...i2Þ and T0ðlðj2Þ...j2Þ, and all the possible values Thiscomplexitysuggeststhatthefusionoperationsmay
of path and path0. be used only for reasonable trees (typically, less than
We must therefore determine the number of values that 100nodes)andsmallvaluesofl(typically,lessthan4).Itis
path can take. This amounts to determine the total number however important to observe that the overall number of
ofsuccessivefusionsthatcouldbeappliedtoagivennode. fusions one may perform can be much greater than l
We recall that path is a list of pairs ðeoru;vÞ. Let path¼ without affecting the worst-case complexity of the algo-
fðeoru;v1Þ;ðeoru;v2Þ;...;ðeoru;v‘Þgbethelistfornodei rithm. Indeed, any number of fusions can be made while
ofT.Thefirstfusioncanbeperformedonlywithachildv1 still retaining the bound of
of i. If d is the maximum degree of T, there are d possible
Oðð2dÞl(cid:4)jTj(cid:4)minðleafðTÞ;heightðTÞÞ(cid:4)jT0j(cid:4)minðleafðT0Þ;
choices for v1. The second fusion can be done with one of
the children of i or with one of its grandchildren. Let v2 be heightðT0ÞÞÞ
the node chosen. There are d + d2 possible choices for v2. so long as one does not realize more than l consecutive
Following the same reasoning, there are Pk¼‘dk possible
k¼1 fusions for each node.
choices for the ‘th node v‘ to be fusioned with i. In general, also, most interesting tree representations of
an RNA are of small enough size as will be shown next,
together with some initial results obtained in practice.
5 APPLICATIONTORNASECONDARYSTRUCTURES
COMPARISON
The algorithm presented in the previous section has been
coded using C++. An online version is available at http://
www-igm.univ-mlv.fr/~allali/migal/.
We recall thatRNAs arerelatively small molecules with
sizes limited to a few kilobases. For instance, the small
ribosomal subunit of Sulfolobus acidocaldarius (D14876) is
Furthermore, we must take into account the fact that a
madeupof1,147bases.Usingtherepresentationshownin
fusion can concern a node or an edge. The total number of
Fig. 2b, the tree obtained contains 440 internal nodes and
values possible for the variable path is therefore:
567leaves,thatis1,007nodesoverall.Usingtherepresenta-
Yk¼‘Xj¼k Yk¼‘dkþ1(cid:1)1 tioninFig.2d,thetreeiscomposedof78nodes.Finally,the
2‘(cid:4) dj ¼2l ; tree obtained using the representation given in Fig. 2e
d(cid:1)1
k¼1j¼1 k¼1 containsonly48nodes.Wethereforeseethatevenforlarge
RNAs,any of the knownabstracttree-representations (that
that is:
is,representationswhichtaketheelementsofthesecondary
2‘(cid:4)(cid:1) 1 (cid:2)‘Yk¼‘ðdkþ1(cid:1)1Þ<2l(cid:4)(cid:1) 1 (cid:2)l(cid:4)dð‘þ1Þ2ð‘þ2Þ: structureofanRNAintoaccount)thatwecanuseleadstoa
d(cid:1)1 d(cid:1)1 treeofmanageablesizeforouralgorithm.Infact,forsmall
k¼1
values of l (2 or 3), the tree comparison takes reasonable
A node i may then be involved in Oðð2dÞlÞ possible time (a few minutes) and memory (less than 1Gb).
successive (node/edge) fusions. Aswealreadymentioned,afusion(respctively,split)can
As indicated, we must store in memory the distance be viewed as an alternative to a deletion (respectively,
betweeneachsubtreeTðlði2Þ...i2ÞandT0ðlðj2Þ...j2Þforall insertion) followed by a relabeling. Therefore, the cost
possible values of path and path0. The space complexity of function for a fusion must be chosen carefully.