Ok, my wrapper function, reroot(), also now can re-root the tree along terminal edges. The way I did this was by re-rooting the tree at the parent node of the target edge; dropping the target tip; then creating a new two taxon "phylo" object with the target edge split between an edge leading to the target tip and the new root edge for the remainder of the tree. Finally, I used bind.tree() to reattach the remaining tree to the two taxon "phylo" object.

The code is simple enough, so I have a reduced version (without error handling, etc.) below. I will post a completed version shortly to my R-phylogenetics page.

Code below:

reroot<-function(tree,node.number,position){

if(node.number>length(tree$tip)){

# first, re-root the tree at node.number

tr<-root(tree,node=node.number,resolve.root=T)

# now re-allocate branch length to the two edges descending

# from the new root node

b<-sum(tr$edge.length[tr$edge==(length(tree$tip)+1)])

tr$edge.length[tr$edge==(length(tree$tip)+1)]<-c(position, b-position)

} else {

# first, root the tree at the parent of node.number

tr1<-root(tree,node=tree$edge[match(node.number, tree$edge[,2]),1])

# now drop tip

tr1<-drop.tip(tr1,tree$tip.label[node.number])

# create phylo object

tr2<-list(edge=matrix(c(3L,1L,3L,2L),2,2,byrow=T), tip.label=c(tree$tip.label[node.number],"NA"), edge.length=c(tree$edge.length[match(node.number,tree$edge[,2])] -position,position),Nnode=1)

class(tr2)<-"phylo"

tr<-bind.tree(tr2,tr1,where=which(tr2$tip.label=="NA"))

}

return(tr)

}

Let's try it out. First, a random birth-death tree:

Here, the blue numbers are internal node numbers; the green numbers are edge lengths; and the numbers in [ square parentheses ] after each tip label are the node numbers for the tips.

First, let's try and re-root this tree 1.0 units of edge length along the branch leading to node #16 (marked with a red hashmark on the figure above):

> test1<-reroot(tree,node.number=16,position=1.0)

> plot(test1)

Cool - it seems to work.

Now, let's try re-rooting the tree along a terminal edge - say the edge leading to taxon "2". We'll do this 8.0 units along this edge (which is ~8.5 units long in total). We can actually do this without knowing the node number for taxon "2" as follows:

> test2<-reroot(tree,node.number=which(tree$tip.label=="2"), position=8.0)

> plot(test2)

And we're done.

I just wrote this simple function this morning, but it is for a project I am presently working on, so if you find any bugs - I would be happy to hear about it.

Nice! I tried it on a several trees and it appears to work well.

ReplyDelete