Tuesday, April 19, 2011

Problem with bind.tree()

Today I ran into a slight problem with the {ape} function bind.tree(). I posted this issue to the R-sig-phylo email listserve. Emmanuel Paradis, the author of {ape}, is normally excellent at responding to queries like this - but in the meantime I thought I would also describe the issue in this forum.

I said in my listserve post: "For some reason, when I bind a tree with a root.edge to a tip, the length of the edge leading to that tip is attached not only to the root.edge of the bound subtree, but also to another tip of the tree." Actually, this is imprecise. In fact, the edge length from the edge leading to the target tip is for some reason moved to another tip in the tree.

The easiest way to illustrate this problem is with an example.

First, take the following tree:

text="(1:3.785,(((((8:0.13,9:0.13):0.205,6:0.335):0.154,(10:0.278,7:0.278):0.211):0.1,(4:0.386,5:0.386):0.203):2.655,(2:1.08,3:1.08):2.164):0.54);"
tree<-read.tree(text=text)
plot(tree)
edgelabels(round(tree$edge.length,2),adj=c(0.5,-0.2),frame="none")


Which produces the plot below in which I have labeled edges with their (rounded) lengths:

Now to reproduce our issue, let's cut out the clade with tips "2" and "3", and then regraft it in the same place:

node<-19 # this is the node we will extract
position<-1.0 # this is the length of the root.edge
# extract clade and attach root edge
tr1<-extract.clade(tree,node=node); tr1$root.edge<-1.0
x11(); plot(tr1,root.edge=T) # plot the extracted clade
edgelabels(round(tr1$edge.length,2),adj=c(0.5,-0.2),frame="none")

The length of the root edge is not shown, but remember we have set it to be 1.0.

# now remove tips in tr1 from tree
tr2<-drop.tip(tree,tr1$tip.label,trim.internal=FALSE)
# subtract root.edge of tr1 from tip ending in "NA"
temp<-match(which(tr2$tip.label=="NA"),tr2$edge[,2])
tr2$edge.length[temp]<-tr2$edge.length[temp]-position
x11(); plot(tr2,root.edge=T) # plot the pruned tree
edgelabels(round(tr2$edge.length,2),adj=c(0.5,-0.2),frame="none")

Here is the tree with clade (2,3) extracted.
# now bind tr1 to tr2, where it was removed
tr.bound<-bind.tree(tr2,tr1,where=which(tr2$tip.label=="NA"))
# now plot, to visualize the result
x11(); plot(tr.bound)
edgelabels(round(tr.bound$edge.length,2),adj=c(0.5,-0.2),frame="none")


So, for some reason the edge leading to "NA" before we use bind.tree(), is actually moving to taxon "5" when we reattach the subtree with species "2" and "3". Weird.

Please let me know if you can repeat this error; or, even better, if you figure out that I have made some kind of mistake in my example, above. Thanks!

In a future post (hopefully soon, if I can get this figured out) I will describe what I have been using extract.clade() and bind.tree() for today as it may be of interest to some readers.

6 comments:

  1. Liam-
    Just letting you know: the code also produces an edge that is the wrong length on my machine also, so its not just you. It'll be interesting to know why this happens.
    -Dave

    ReplyDelete
  2. Thanks for giving it a whirl. Did you see any obvious (or, better yet, unobvious) errors in my code?

    If it is genuinely a bug in bind.tree() I'm sure that Emmanuel will respond fairly quickly to my R-sig-phylo query as he tends to be very good about that kind of thing.

    ReplyDelete
  3. Liam-
    Your code looks fine. Any response from EP yet?
    -Dave

    ReplyDelete
  4. Hi Liam,
    I know this post is very old and you probably already figured out what was wrong.. In any case, I tested your code and the tree seemed to be reassembled properly - looks like Emmanuel was able to fix the issue in the updated version of ape :)

    Regards,
    Ania

    ReplyDelete
  5. Hi Ania.

    Thanks for letting us know!

    All the best, Liam

    ReplyDelete