Wednesday, April 6, 2011

All furcating trees

Brian O'Meara posted a while ago to R-sig-phylo that he was looking for code to find all bi- and multifurcating rooted trees for a set of taxa. I realized that a function to generate all unrooted trees would be straightforward to program using my add.everywhere() function described previously. [Adding rootedness should then merely involve going through each tree and placing the root along each edge.]

I ran into a few difficulties doing this, one being that the ordering of $edge was somehow incompatible with plot.phylo() for some trees (but I couldn't figure out how). To get around this, I added the option to.plot, which if set to TRUE uses:

trees<-read.tree(text=write.tree(trees))

Unfortunately, my code for this (called allFurcTrees() and available on my R-phylogenetics page, direct link here) runs very slowly. For instance, even with to.plot=FALSE it takes 7 minutes to run for 8 taxa. I have passed the code along to Brian, and hopefully if he intends to use it, he will be able to make it more efficient. I did use the code to create the neat movie (showing all possible bi- and multifurcating trees for 7 taxa), below.



To make your own, just use the following code (or some variant). This code will create a large number (2,752 in this case) .png files, which can then be combined as frames in your movie using a standard movie-making software (for instance, the Windows Live Movie Maker that is free with Windows 7). I'm not sure if there is anyway to generate the animation directly from R.

source("allFurcTrees.R")
trees<-allFurcTrees(n=7,c("A","B","C","D","E","F","G"), to.plot=TRUE)
for(i in 1:length(trees)){
plot(trees[[i]],type="unrooted")
savePlot(file=paste("tree",i,".png",collapse="",sep=""), type="png")
}

6 comments:

  1. Thanks Liam, the little video is quite a nice illustration of the size of the sample space! If you haven't seen it already, you might also be interested in the animation package in R, see, for instance: http://animation.yihui.name/animation:start#gif_animation

    ReplyDelete
  2. Hi Carl.
    Thanks for the hint. I will check it out.
    By the way, the quality of the video created by Movie Maker is actually quite high - the low quality of the posted version is due to blogger.com.

    ReplyDelete
  3. Liam-
    Neat!

    I've run into a number of problems lately with some phylo objects not wanting to run in plot.phylo or prop.part. I've been using the same read.tree(write.tree()) trick, which seems to fix it, but it isn't clear to me what is going wrong.
    -Dave

    ReplyDelete
  4. Hi Dave. I don't know either, but I think it has something to do with the order of the rows in $edge. I believe, for example, that write.tree() will work for any order; but some other functions require a particular ordering - for instance, "ancestor" edges before their "descendants" (even though ancestor & descendant are not meaningful on an unrooted tree). I have tried to use the "ape" function reorder.phylo(), but it doesn't seem to work in this case. The work-around read.tree(write.tree(...)) is quite slow.

    ReplyDelete
  5. This post was featured on Bioephemera, a popular blog about "the intersection of art and biology." Link to blog post here

    ReplyDelete
  6. I also posted a higher quality version of this to youtube (link here).

    ReplyDelete