Friday, April 1, 2011

Adding a tip in all possible places on a tree

An R-sig-phylo list server, Matthew Vavrek, just posted to inquire as to whether anyone could suggest an algorithm to add a tip to the tree in all possible places (i.e., along each edge of the tree). I contributed my solution, which works by using the {ape} function bind.tree(); but I wonder if anyone else can suggest a more elegant version. The function bind.tree() is designed to merge two trees; the trick here is that we want to add a tip - not a tree.

This is how I did it.

Say, we start with a random, unrooted tree with 4 species. Of course, 4 is arbitrary here - we could start with a tree of any length.

tree<-rtree(n=4,rooted=FALSE,br=rep(1,5)) # random tree
# create a 5th species, here "t5", to add as a "phylo" object
# [I don't think this can be avoided with bind.tree()]
new.tip<-list(edge=matrix(c(2,1),1,2),tip.label="t5", edge.length=1,Nnode=1)
class(new.tip)<-"phylo"
# we have "tricked" R into thinking that new.tip is a tree
# add the new tip to all edges of the tree
trees<-list(); class(trees)<-"multiPhylo"
for(i in 1:nrow(tree$edge))
   trees[[i]]<-bind.tree(tree,new.tip,
    where=tree$edge[i,2],position=0.5)
# now plot them to see what we have done
plot(trees,type="unrooted",use.edge.length=F)


And that's it - we have created (and plotted) our four taxon starting tree with a fifth species added in all possible positions!

Of course, we could generalize this with the following function:

add.everywhere<-function(tree,tip.name){
   tree<-unroot(tree)
   tree$edge.length<-rep(1,nrow(tree$edge))
   new.tip<-list(edge=matrix(c(2,1),1,2),tip.label=tip.name,
    edge.length=1,Nnode=1)
   class(new.tip)<-"phylo"
   # add the new tip to all edges of the tree
   trees<-list(); class(trees)<-"multiPhylo"
   for(i in 1:nrow(tree$edge)){
      trees[[i]]<-bind.tree(tree,new.tip,
       where=tree$edge[i,2],position=0.5)
      trees[[i]]$edge.length<-NULL
   }
   return(trees)
}


After loading the function above, and {ape}, let's try it:

> tree<-read.tree(text="((George,Paul),(Ringo,John));")
> trees<-add.everywhere(tree,"Pete_Best")
> plot(trees,type="unrooted")


Cool.

1 comment:

  1. After making a small change suggested (here) by Emmanuel Paradis, I have now posted this function online. Direct link to code is here.

    ReplyDelete