Reversing Gene Erosion—Reconstructing Ancestral Bacterial Genomes from Gene-Content and ...

13 Pages·2012·0.22 MB·English

Reversing Gene Erosion(cid:151)Reconstructing Ancestral
Bacterial Genomes from Gene-Content and Order Data
1 2 1
JoelV.Earnest-DeYoung ,EmmanuelleLerat ,andBernardM.E.Moret
1
Dept.ofComputerScience,Univ.ofNewMexico,Albuquerque,NM87131,USA,
joeled,[email protected]
2
Dept.ofEcologyandEvolutionaryBiology,Univ.ofArizona,Tucson,AZ85721,USA,
[email protected]
Abstract. Inthelastfewyears,ithasbecomeroutinetousegene-orderdatato
reconstructphylogenies,bothintermsofedgedistances(parsimonioussequences
of operations that transform one end point of the edge into the other) and in
termsofgenomesatinternalnodes,onsmall,duplication-freegenomes.Current
gene-order methodsbreakdown,though,whenthegenomescontainmorethan
afewhundredgenes,possesshighcopynumbersofduplicatedgenes,orcreate
edge lengths in the tree of over one hundred operations. We have constructed
aseriesofheuristicsthatallowustoovercome theseobstacles and reconstruct
edgesdistancesandgenomesatinternalnodesforgroupsoflarger,morecomplex
genomes.Wepresentresultsfromtheanalysisofagroupofthirteenmodern(cid:13)-
proteobacteria,aswellasfromsimulateddatasets.
1 Introduction
Althoughphylogeny,theevolutionaryrelationshipsbetweenrelatedspeciesortaxa,is
a fundamentalbuildingblockin muchof biology,it has beensurprisinglydif(cid:2)cultto
automate the process of inferring these evolutionary relationships from modern data
(usually molecular sequence data). These relationships include both the evolutionary
distanceswithinagroupofspeciesandthegeneticformoftheircommonancestors.
Inthelastdecade,anewformofmoleculardatahasbecomeavailable:gene-content
andgene-orderdata;thesenewdatahaveprovedusefulinsheddinglightonthesere-
lationships [1(cid:150)4].The orderandthe orientationin which genes lie ona chromosome
changesveryslowly,inevolutionaryterms,andthustogetherprovidearichsourceof
informationforreconstructingphylogenies.Until recently,however,algorithmsusing
suchdatarequiredthatallgenomeshaveidenticalgenecontentwithnoduplications,re-
strictingapplicationstoverysimplegenomes(suchasorganelles)orforcingresearchers
toreducetheirdatabyequalizingthegenecontent(deletingallgenesnotpresentinev-
erygenomeandall(cid:147)copies(cid:148)ofeachgene,e.g.,usingtheexemplar strategy[5]).The
formerwas frustratingto biologistswantingtostudymorecomplexorganisms,while
thelatterresultedindatalossandconsequentlossofaccuracyinreconstruction[6].
Ourgrouprecentlydevelopedamethodtocomputethedistancebetweentwonearly
arbitrarygenomes[7,?]andanothertoreconstructphylogeniesbasedongene-content
and gene-orderin the presence of mildly unequal gene content [6]. In this paper, we
bringtogetherthesemethodsinaframeworkthatenablesustoreconstructthegenomes
2
Pasteurella multocida
Haemophilus influenzae
Yersinia pestis−CO92
Yersinia pestis−KIM
Salmonella typhimurium
Escherichia coli
Wigglesworthia brevipalpis
Buchnera aphidicola
Vibrio cholerae
Pseudomonas aeruginosa
Xylella fastidiosa
Xanthomonas axonopodis
Xanthomonas campestris
Fig.1.The13gamma-proteobacteriaandtheirreferencephylogeny[8].
ofthecommonancestorsofthe13modernbacteriashowninFig.1(from[8]).Gamma-
proteobacteriaareanancient groupofbacteria,at least 500millionyearsold[9];the
groupincludesendosymbiotic,commensal,andpathogenicspecies,withmanyspecies
playinganimportantmedicaloreconomicrole.Theevolutionaryhistoryofthegroupis
quitecomplex,includinghighlevelsofhorizontalgenetransfer[10(cid:150)12]and,inthecase
ofB.aphidicolaandW.brevipalpis,massivelevelsofgeneloss.Thesefactorsmakea
phylogeneticanalysisofthisgroupbothinterestingandchallenging.
Therestofthispaperisorganizedasfollows.Section2presentstheproblem.Sec-
tion 3 summarizes prior work on phylogeneticreconstructionfrom gene-content and
gene-orderdata.Section4presentsourframeworkfortacklingtheproblemofancestral
genomereconstructiongivenareferencephylogeny;itisitselfdividedintothreesub-
sections,oneon eachofourthreemaintools: median-(cid:2)nding,contentdetermination,
andgeneclustering.Section5discussesourapproachtothetestingofourframework:
giventhatwehaveonlyonedatasetandthatancestralgenomesforthatdatasetareen-
tirelyunknown,ourtestingwas ofnecessitybasedonsimulations.Section6 presents
theresultsofthistesting.
2 The Problem
Wephrasethereconstructionproblemintermsofaparsimonycriterion:
Giventhegeneordersofagroupofgenomesandgivenarootedtreewiththese
genomesat theleaves,(cid:2)ndgeneordersforthe internalnodesofthetreethat
minimizethesumofalledgelengthsinthetree.
Thelengthofanedgeisde(cid:2)nedintermsofthenumberofevolutionaryevents(permis-
sibleoperations)neededtotransformthegenomeatoneendoftheedgeintothegenome
attheotherend.Thepermissibleoperationsinourcaseareinversions,insertions(and
duplications),anddeletions;alloperationsaregiventhesamecostincomputingedge
lengths.Restrictingrearrangementstoinversionsfollowsfrom(cid:2)ndingsthattheinver-
sion phylogenyis robust evenwhenotherrearrangements,such as transpositions,are
3
usedincreatingthedata[13].Ourassignmentofunitcoststoalloperationssimplyre-
(cid:3)ectsinsuf(cid:2)cientbiologicalknowledgeabouttherelativefrequencyoftheseoperations.
Inoursetting, one insertionmay addanarbitrarynumberof genesto a single lo-
cationandonedeletionmayremoveacontiguousrunofgenesfromasinglelocation,
a convention consistent with biological reality. Gene duplications are treated as spe-
cialized insertionsthat onlyinsert repeats.Finally,oneach edgea genecan either be
insertedordeleted,butnotboth;thesameholdsformultiplecopiesofthesamegene.
Allowingdeletionandinsertionofthesamegenesonthesameedgewouldleadtobio-
logicallyridiculousresultssuchasdeletingtheentiresourcegenomeandtheninserting
theentiretargetgenomeinjusttwooperations.
Finding internal labels that minimize edge distances over the tree has been ad-
dressed byourgroup(cid:151)thisis the mainoptimizationperformedbyour softwaresuite
GRAPPA [14]. However, even the most recent version of GRAPPA [6] is limited to
relativelysmallgenomes(typicallyoforganellarsize,withfewerthan200genes),with
modestly unequal content and just a few duplications. In stark contrast, the bacterial
genomes in our dataset contain 3,430 different genes and range in size from 540 to
2,987genes,withsevencontainingover2,300genes;moreover,thesegenomescontain
alargenumberofduplications,rangingfrom3%to30%ofthegenome.Thus,inour
model,mostpairwisegenomicdistancesareverylarge:asimplepairwisecomparison
alongthetreeofFig.1indicatesthatsomeedgesofthetreemustrepresentatleast300
events.SuchlengthsareatleastanorderofmagnitudelargerthanGRAPPAcanhandle.
Thelargegenomesize,vastlyunequalgenecontent,largenumberofduplications,and
largeedgelengthsallcombinetomakethisdatasetordersofmagnitudemoredif(cid:2)cult
toanalyzethanpreviouslyanalyzedgenomesets.
3 PriorWork
Athoroughrecentreviewofthecurrentworkinphylogeneticreconstructionbasedon
genecontentandgeneorderappearsin[15];wereviewonlytherelevantpointshere.
TheGRAPPAsoftwarepackage[16]computesinternallabelsintwophases.First,
itinitializesinternallabelsofthetreebysomemethod.Thenititerativelyre(cid:2)neslabels
until convergence:each newly labeled (or relabeled) node is pushed on a queue and,
while the queue is not empty, the node at the head of the queue is removed, a new
labelcomputedforit(bycomputingthemedianofitsthreeneighbors),and,ifthenew
labelreducesthetotaldistancetothethreeneighbors,theexistinglabelisreplacedwith
the improved label and the three neighbors are placed on the queue. Thus GRAPPA
relies on the computation of the median of three genomes, that is, a fourth genome
whichminimizesthesumofthenumberofoperationsneededtoconvertitintoeachof
thethreegivengenomes.GRAPPA(cid:2)ndsoptimalinversionmedianswithanalgorithm
that runs in worst-case exponential time, but (cid:2)nishes quickly when the edge lengths
aresmall (10to40operationsperedge)[6,17].GRAPPA treats groupsofgenesthat
occur in the same order and orientation in all genomes as a single genetic unit; this
condensationstepreducescomputationalcostsanddoesnotaffectthe(cid:2)nalresult[18].
Ourgroupdevelopedamethodto(cid:2)ndthedistancebetweentwogenomeswithar-
bitrarygenecontent[7,?];thismethodreliesonaduplication-renamingheuristicthat
4
matches multiple copies of genes between genomes and renames each pair and each
unmatchedcopytoanew,unusedgenenumber.Thusarbitrarygenomesareconverted
intoduplication-freegenomes.Weprovedthat,giventwogenomeswithunequalgene
contentandnoduplications,anyoptimalsortingsequencecanberearrangedtocontain
(cid:2)rstallinsertions,thenallinversions,and(cid:2)nallyalldeletions(cid:151)atypeofnormalform
for edit sequences [7].(Deletions here are genes uniqueto the source genome,while
insertionsaregenesfoundonlyinthetargetgenome.) Usingthegenomesproducedby
the duplication-renamingmethod,an optimalinversionsequencecanbe calculatedin
timequadraticinthesizeoftheconsensusgenomes[19,20].Thenumberofdeletions
iscalculatedbycountingthenumberofHannenhalli-Pevznercyclesthatcontaindele-
tions,asdescribedin[21].Finally,thenumberofinsertionsisestimatedbycalculating
allpossiblepositionsinthesourcegenometowhichtheinversionsequencecouldmove
insertions,thenchoosingthe(cid:2)nalpositionforeachinsertionthatminimizesthenumber
ofgroupsofinsertedgenes.
Insome genomes,especiallybacterial ones,geneswith similar functionareoften
locatedtogetherononestrandofachromosome;thesefunctionalunitsarecalledoper-
ons.Inbacteria,atleast, whiletheorderofgenesin anoperonmaychange,thegene
contentoftheoperonis muchless likelytodoso[22].Ingene-orderdata,anoperon
appearsas a cluster of genenumberswith thesame sign, with content,but notorder,
preservedacrossgenomes.HeberandStoyedevelopedalinear-timecluster-(cid:2)ndingal-
gorithmtoidentifytheseoperon-likeclusterswithinequal-contentgenomes[23].
McLysaght et al. [4] reconstructed ancestral genomes for a group of poxviruses;
shedeterminedgenecontentbyassumingthatthephylogenetictreecontainedasingle
pointoforiginforeachgenefamilyinthemoderngenomes.Eachpointoforiginwas
assignedtothatinternalnodewhichminimizedthenumberoflosseventsnecessaryto
achievethegenecontentoftheleafgenomes.
4 DesigninganAlgorithmicFramework
Toaddresstheproblemofreconstructingancestralgenomesatthelevelofcomplexity
ofgamma-proteobacteria,weusecondensationofgeneclustersinordertoreducethe
sizeofthegenomes,describeaproceduresimilartothatofMcLysaghtetal.todeter-
minethegenecontentofeveryinternalnode,andpresenta newheuristic to compute
themedianofthreeverydifferentgenomes.
4.1 Medians
Weusethequeue-basedtree-labelingheuristicdescribedinSection3.Sinceleavescon-
taintheonlylabelsknowntobecorrect,weupdatethenodesinorderoftheirdistance
fromtheleaves,as showninFig.2.Theheartofthetop-levelheuristicis themedian
computation.Exactmedian-(cid:2)ndingalgorithmsarelimitedtosmallgenomes,smalledge
lengthsinthetree,andfewchangesincontent(cid:151)andnoneofthesepropertiesholdsin
ourproblem.Wethereforepursueasimpleheuristicinspiredbygeometry.Themedian
ofatriangleintheplanecanbefoundbydrawingalinefromonevertextothemiddle
oftheoppositesegment,thenmovingtwothirdsofthewayalongthisline.Byanalogy,
5
R
4
3
2
1 1 1 1
L L L L L L L L L
Fig.2.Internalnodesorderedbytheirdistancefromtheleaves.Nodeswithlowerindiceswillbe
labeled(cid:2)rst;nolabelisgeneratedfortheroot.
wegenerateasortingsequencefromonegenometoanother(anedgeofthetriangle),
thenchooseagenomehalfwayalongthissortingsequenceandgenerateanewsorting
sequencefromittothethirdgenome,stoppingone-thirdalongtheway.
We extend the method of Marron et al. [7] to enumerate all possible positions,
orientations,andorderingsofgenesaftereachoperation.Deletedgenesattheendpoint
ofaninversionaremovedtotheotherendpointifdoingsoavoids(cid:147)trapping(cid:148)thedeleted
genes between two consensus genes that are adjacent in the target genome. Inserted
genesaremovedsoastoremainadjacenttooneofthetwoconsensusgenesbetween
which they lie in the target genome.We can thus generate the genomesproducedby
(cid:147)running(cid:148)a portionofthe sortingsequence,thenuse these intermediategenomesfor
themedianheuristicjustdescribed,allinpolynomialtime.
Thishandlingofinsertedgenesleadstoanoverestimateoftheeditdistance,which
Marron et al. showed at most doubles the number of operations [7]. Their original
method calculates all possible positions in the source genome to which the inversion
sequence could moveinsertions and chooses the (cid:2)nal position (for each insertion)to
minimizethenumberofgroupsofinsertedgenes;itmayunderestimatetheeditdistance
becausethegroupingofinsertedgenesmayrequireaninversiontojoininsertedgenes
andsimultaneouslysplit deleted genes,which is not possible. We comparedpairwise
distancesproducedbyourmethodandbytheirstogetanupperboundontheoveresti-
mation:averageandmaximumdifferencesbetweentheoverestimateandunderestimate
were11.3%and24.1%,respectively.
4.2 GeneContent
Wepredeterminethegenecontentofeveryinternalnodebeforecomputinganymedian:
oncethegenecontentofaninternalnodeisassigned,itremainsunchanged.Sincethe
tree is rooted,we knowthe directionoftime (cid:3)owoneach treeedge;we also assume
thatdeletionsarefarmorelikelythaninsertions,Thenumberofcopiesofeachgeneg
isdecidedindependentlyofallothergenes;atinternalnodei,itissettothemaximum
numberofcopiesofgfoundinanyoftheleavesini’ssubtreeif:(i)thereareleavesboth
insideandoutsidei’ssubtreethatcontaing;or(ii)thereareleavescontaininggineach
halfofi’ssubtree.Otherwisethenumberofcopiesofgeneginnodeiissettozero.
ThisvaluecanbecalculatedinO(NG) time,whereN is thenumberofnodesin
thetreeandGisthenumberofdistinctgenesinalltheleaves,asfollows.Foreachnode
inthetree,wedeterminethemaximumnumberofcopiesofeachgenefromamongthe
6
0
0
2
0 1
0 0 0 2 1 0
Fig.3.Thenumberofcopiesofageneininternalnodes.
leavesofthatnode’ssubtree,usingasingledepth-(cid:2)rsttraversal.Weuseaseconddepth-
(cid:2)rsttraversaltosettheactualnumberofcopiesofeachgeneateachinternalnode.If
eitheroftheroot’schildrenhasasubtreemaximumofzero,thenwesettheroot’sactual
numbertozeroaswell.Foreachinternalnodeotherthantheroot,ifitsparent’sactual
numberofcopiesiszeroandatleastoneofitstwochildren’ssubtreemaximumiszero,
thenwesetthenumberofcopiesforthegenetozero;otherwisewesetthenumberof
copiestothenode’ssubtreemaximumforthegene.
Internalnodes will thus possess at least as manycopies ofa geneas the majority
consensus of their neighbors’ gene contents. An internal node will always possess a
copy of a gene if two or more of its neighbors do. (We consider the two children of
the rootto be neighbors.) Moreover,if a nodeis the nearest commonancestorof all
genomespossessingthegene,itmayhavemorecopiesofthegenethanitsparentand
oneofitschildren,asinthecaseoftheblacknodeinFig.3.Thegenecontentofinter-
mediategenomesalongsortingsequencesisaunionofthegenecontentsofthestarting
genomes,becausethesortingsequenceofoperationsthatweusealwaysinvolves(cid:2)rst
insertions,theninversions,and(cid:2)nallydeletions.Therefore,whencalculatingmedians
fromsortingsequences,weface threecases in whichthenumberofcopiesofa gene
differbetweentheintermediategenome,themediangenome,andthemedian’sparent(cid:151)
see Fig.4. InFig.4a,the intermediategenomehas thesame numberofcopies as the
median,butfewerthantheparent,aswiththeblacknode’srightchildinFig.3.Each
copy in the parentthat is not matchedby the duplication-renamingalgorithmwill be
excludedfromthe mediangenome.The case of Fig. 4b onlyarises whenthe median
genome is the nearest common ancestor of all genomes containingthe gene in ques-
tion,aswiththeblacknodeinFig.3.Genomesalongtheintermediatesequencehave
thesamenumberofcopiesasthemedian,whiletheparentofthemediancontainsno
>m 0 0.
PSfragreplacements
PSfragreplacements PSfragreplacements
m m 0
m m (cid:20)m m m (cid:20)m 0 >0 >0
0
(a) (b) (c)
Fig.4.Caseswherethemediananditsneighborshavedifferentnumbersofcopiesofagene.Solid
1 1
linesaretreeedges;dashedanddottedlinesarefractions( and ,resp.)ofsortingsequences.
2 3
7
copyatall.Finally,thecaseofFig.4ccanonlyarisewhentherightchildofthemedian
isthenearestcommonancestorofallgenomescontainingthegene,aswiththeparent
oftheblacknodeinFig.3.
Biologically,thisprocessof(cid:2)ndingwhichduplicatestoincludeinthemediancor-
responds to matching orthologous duplicates of each gene between genomes and to
discarding unmatchedparalogousduplicates. Since the original nucleotidesequences
are abstracted away beforethe analysis begins,this orthologmatchingis decided en-
tirely on the basis of which other genes are located next to the different homologs.
Fortunately, orthologs and paralogs that can be distinguished by a nucleotide-based
analysisareassigneddifferentgenenumbersbeforeouranalysisbegins.Therefore,our
methodrepresentsareasonablewaytointegratebothnucleotideandgene-orderdatain
differentiatingorthologousandparalogoushomologsofgenes.
4.3 ClusterCondensation
Toextractinformationfromlargerandmorecomplexbiologicaldatasets,weneedfast
algorithmswithfastimplementations;fasterprocessingalsoenablesa morethorough
analysisandthusproducesresultsofhigherquality.Thekeyfactorhereisthesizeof
thegenomes(cid:151)theirnumberisamuchsmallerissue.Wethusdevelopedatechniqueto
identifyandcondensegeneclustersinordertoreducethesizeofthegenomes.Ourap-
proachgeneralizesthatusedingenomeswithequalcontent[23];incontrast,GRAPPA
onlycondensesidenticalsubsequencesofgenes,becauseitaimstopreservetheidentity
ofeditsequences.Ourmethodallowsthecondensationofclustersbasedonlyoncontent
(notorder,atleastaslongasgenesstayonthesamestrand)andalsohandlesthedif(cid:2)cult
casesthatariseoutofunequalgenecontent(suchasaninsertionwithinacluster).
Toidentifyclusters,we(cid:2)rstusetheduplication-renamingtechniqueofMarronetal.
tocreateduplication-freegenomes.Afterrenaming,weremoveanygenesnotpresent
in all of the genomes under examination.This step creates a group of genomes with
equalgenecontent.Wethenusethecluster-(cid:2)ndingalgorithmofHeberandStoye[23]
to (cid:2)nd equivalent clusters of genes within the equal-content genomes. Once clusters
are identi(cid:2)ed, each one is condensed out of the original genomes and replaced with
a single marker (as if it were a single gene). In a set of genomes with unequal gene
content, there can be genes inside a cluster that are not present in the corresponding
equal-contentgenomes.Wedealwiththesegenesinoneoftwoways.Ifeveryoccur-
rence of that gene is located inside the cluster in each of the genomes that possesses
thegene,thenthegeneiscondensedalongwiththerestofthecluster.Otherwise,the
extrageneismovedtoonesideoftheclusterandtheclustercondensed.Whena me-
diangenomeiscomputed,amedianforeachclusterisalsocomputedandeachcluster’s
marker in the median genome is replaced with the cluster’s median. At this point, if
anyextragenesmovedtothesideoftheclusterarestillbesideit,theyaremovedback
insidetheclustertoapositionsimilartotheiroriginalone.
4.4 PuttingItAllTogether
Ancestral genomereconstructionsare performedusing these three main components.
Initialization of the internal nodes of the tree is done from the leaves up by taking
8
eitherthemidpointoroneofthetwoendpoints(alongtheinversionportionofanedit
sequence)ofaninternalnode’stwochildrenanddiscardinganygenesnotallowedby
themediangenecontent.ThismethodaccountsforallthreeofthecasesinFig.4and
produceslabelswiththedesiredgenecontent.Newmediansarecomputedlocallynode
by node in a postorder traversal of the tree, so as to propagate information from the
leaves towards the root. Whenever a median is found that reduces the local score at
a node, it immediately replaces the previous label at that node; that node and all its
neighborsarethenmarkedforfurtherupdate.
5 Testing
Weusedourlabelreconstructionmethodonthebacterialdatasetaswellasonsimulated
datasets. With simulated datasets, we know the true labels for the internal nodes as
wellastheexactevolutionaryeventsalongeachedge,sothatwecantesttheaccuracy
of the reconstruction. The goal of our experiments was to generate datasets roughly
comparabletoourbiologicaldatasetsothatourexperimentalresultswouldenableus
topredictarangeofaccuracyfortheresultsonthebiologicaldataset.
ThesimulateddatawerecreatedusingthetreeofFig.1;edgelengthswereassigned
tothetreebasedonourbestestimateoftheedgelengthsforthebacterialgenomes.To
keepthedataconsistent,edgelengthswereinterpretedasthenumberofoperationsper
generatherthanasanabsolutenumber,allowingustousethesamerelativevaluefor
genomesof differentsizes. The tree was labeled by (cid:2)rst constructinga root genome.
Thenumberofgenesgandthetotalsizenoftherootgenomeweresetasvariableuser
parameters.Oneofeachgenefrom1tog wasaddedtotherootgenome,afterwhich
n(cid:0)gadditionalgeneswerechosenuniformlyatrandomintherange[1;g]andadded
totherootgenome.Therootgenomewasthenrandomlypermutedandeachgeneas-
signedarandomsign.Theothernodeswerethenlabeledfromtherootbyevolvingthe
genomesdownthetreeaccordingtotheprescribednumberofoperations.Theallowed
operationswereinsertions,deletions,andinversions.Althoughthetotalnumberofop-
erationswas(cid:2)xed,theproportionofeachofthethreetypesofoperationswasleftasa
variableparameterbysettingtheratioofinversionstoinsertionstodeletions.Thismix
ofoperationswasusedoveralledgesofthetree.
Thecharacteristicsofeachtypeofoperationweredeterminedseparately.Thelength
ofeachinversionwas chosenuniformlyat randombetween1andhalfthesize ofthe
genome, with a start point chosen uniformly at random from the beginning of the
genometo the size of the genomeminus the lengthof the inversion.The average in-
sertionlengthwassetviaauserparameterasaportionofthesizeoftherootgenome
andwasusedunchangedovertheentiretree,whiletheactuallengthofeachinsertion
wasdrawnfromaPoissondistributionwiththisexpectationanditslocationwaschosen
uniformlyatrandomfromthebeginningtotheendofthegenome.Inmovingfromthe
roottotheleaves,itwasassumedthataparticulargenecouldonlybeinsertedalongone
edgeofthetree(cid:151)multipleinsertionsofthesamegene,evenalongseparatepaths,were
notallowed.Theaveragedeletionlengthwaschosenasauser-speci(cid:2)edportionofthe
genomefromwhichgeneswerebeingdeleted,thusvaryingfromedgetoedgeaswell
asalongeachedgewitheachsuccessivedeletion,whiletheactualsizeofeachdeletion
9
was drawnfroma Poisson distributionwith this expectationandwith a start location
chosenuniformlyatrandomfromthebeginningofthegenometothesizeofthegenome
minusthelengthofthedeletion.Withtheconstantexpectedinsertionlength,genomes
growlinearlyintheabsenceofdeletions,while,withaproportionalexpecteddeletion
length,genomesshrinkexponentiallyintheabsenceofinsertions.Whenbothinsertions
anddeletionsareused,genomesfartherfromtheroottendtowardsastablesize.Along
eachedge,theprescribednumberofinsertionsareperformed(cid:2)rst,theninversions,and
(cid:2)nallydeletions.Onceall nodes havebeenbeenassignedgenomes,theresultingleaf
genomesare fed into our reconstructionprocedure.The results of the reconstruction,
intermsofgenecontentandgeneorderateachinternalnode,arecomparedwiththe
(cid:147)true(cid:148)tree,i.e.,thatgeneratedinthesimulation.
Weconstructedtreesusing(cid:2)vedifferentmodels:an(cid:147)inversion-only(cid:148)model,a(cid:147)no-
deletions(cid:148) modelwith a 6:1inversion-to-insertionratio,a (cid:147)no-insertions(cid:148)modelwith
a 6:1 inversion-to-deletionratio, a (cid:147)low-insertion/deletion(cid:148)model with a 40:4:1 ratio
of inversions to deletions to insertions, and a (cid:147)high-insertion/deletion(cid:148) model with a
30:10:3ratioofinversionstodeletionstoinsertions.Theaverageinsertionlengthwas
setto2%oftherootgenomeandtheaveragedeletionlengthto3%ofthelocalgenome.
In order to test the ef(cid:2)cacy of cluster condensation, we tested the technique on
triples amongthe bacterial genomesthat lie close to each otheron the tree in Fig. 1.
Tripleswerechosenbyselectinginternalnodes,then,foreachofthethreeedgesleading
outfromtheinternalnode,bychoosinganearbyleafreachablebyfollowingtheedge.
Foreachsetofthreegenomes,wemeasuredthesumofthelengthsofallclustersthat
werefound.
6 Results
Our discussion and summaries of results refer to Fig. 5. Reconstruction of ancestral
genomesforthebacterialgenomestakesaround24hoursonatypicaldesktopcomputer.
Themidpoint-initializationprovedquitestrong:theonlygenomestobeupdatedinthe
14
6 Pasteurella multocida
7 15
Haemophilus influenzae
16
10 Yersinia pestis−CO92
4 8 17
3 8 Yersinia pestis−KIM
5 18
11 Salmonella typhimurium
7 9 19
2 4 Escherichia coli
2 20
9 Wigglesworthia brevipalpis
10 21
1 Buchnera aphidicola
5
Vibrio cholerae
3
1 Pseudomonas aeruginosa
12
Xylella fastidiosa
6 13 22 Xanthomonas axonopodis
11 23
Xanthomonas campestris
Fig.5.Thebacterialtreewithnumberededgesandinternalnodes.
10
Table1.ErrorPercentageinTreeScores
AvgerrorMinerrorMaxerror
Inversiononly 63.2% 57.3% 67.4%
Nodeletions 62.6% 54.8% 70.7%
Noinsertions 45.2% 37.6% 54.3%
Lowinsertion/deletion 56.4% 46.7% 64.8%
Highinsertion/deletion 34.9% 25.1% 46.4%
subsequentlocalimprovementprocedurewerethetwochildrenoftheroot(nodes1and
6inFig.5),theonlyneighboringgenomesinwhichoneneighborwasnotusedtocre-
atetheother.Whenweusedendpoint-initialization,threeinternalnodeswereupdated
(nodes3,4,and6inFig.5)andthescoreoftheentiretreewasloweredby2.8%.This
(cid:2)ndingmayindicatethattheinitializationisverygood,butitmayalsore(cid:3)ectthelarge
numbersoflocaloptimainthesearchspace(cid:151)asimilar(cid:2)ndingwasreportedforthesim-
plerGRAPPA[18].Itshouldbenotedthat,whencalculatingmedians,onlyfourdiffer-
entmidpointsinthechild-to-childsortingsequenceareused;fromeachofthesemid-
points,onlythreemidpointsinthesortingsequencefromtheintermediategenometo
theparentaretested.Thusweonlyperformaveryshallowsearchandcouldeasilymiss
a better solution. Interestingly,though,when we did a slightly more thoroughsearch
withtenmidpointsfromchildtochildandfourmidpointsfromintermediatetoparent,
using endpoint-initialization, the tree score was slightly worse than in the shallower
analysis,althoughthesearch,whichtookabout3.5timeslonger,updatedthesamethree
internalnodes.Ofcourse,thislargersearchremainsveryshallow;goingbeyonditwill
requireamuchmoreef(cid:2)cientimplementationoftheduplication-renamingheuristicof
Marronetal.[7](cid:151)inourcurrentversion,itusesupover90%ofthecomputingtime.
We simulated 100 labelings of the tree with a root genomesize of 200 genes for
eachofthe(cid:2)vepreviouslydescribedscenarios.Endpoint-initializationwasusedinall
scenarios.Theleafgenomesproducedinoursimulationsrangedinsizefrom70genes
to 400 genes.We comparedthepredictedgenecontentof the internalnodes with the
actual gene content.As expected(due to our restriction on generation),the predicted
genecontentalwaysmatched,exceptwhenagenecopythatwaspresentataninternal
nodewaslostinallleaves.Failuretodetectthiskindofmissinggeneisunavoidablein
ananalysissincethedeletionfromallleavesmeansthatnohistoricalrecordis leftto
attestthepresenceofthatgeneinancestralgenomes.Whenwecomparedthenumber
of operations over all edges in reconstructed trees versus the original simulated tree,
thescoreforthetreewasfairlyinaccurate,consistentlyoverestimatingthetruescore,
asillustratedinTable1.Therathertightdistributionforthetreeindicatesthattheerror
isnotarandomprocess,butaresultofsomeaspectofourreconstructionmethod,one
thatmaylenditselftoreversemapping.
We compared edge lengths in the reconstructed trees with those in the true trees
bycalculatingtheratioofthelengthsforeachedge(Fig. 6).Aperfectreconstruction
wouldgivearatioof1.0;asthe(cid:2)gureshows,mostratiosarehigher,withedgesfurther
fromleaveshavinglargerratios(andalsolargervariances).Abouthalfofthe23edges
arewithinafactoroftwoandanotherquarterarewithinafactoroffour.