Tuesday, June 28, 2011

Slicing a tree to create all subtrees

Someone recently emailed me about slicing a tree at a particular height from the root and then keeping all the subtrees. This seemed pretty easy, so I thought I would implement it using the {ape} function extract.clade(). This is how I did it. First, I reordered the tree "cladewise", using reorder.phylo(); next, I computed the heights of all the internal and terminal nodes in the tree and put these values into a matrix to match the $edge matrix from my "phylo". Then, I found all the edges intersected by my slice and identified the tipward nodes on those edges, excluding all the nodes that were also tips. Finally, I extracted each subtree descended from the nodes using extract.clade() and added a $root.edge to each tree depending on where its root branch was sliced. I then returned all of these trees as a "multiPhylo" object.

The function is posted to my website, here. Please be warned that I'm not sure what will happen if the slice point is greater than the tree length or if it only intersects terminal edges.

First, let's create a random tree:

> require(geiger)
> tree<-birthdeath.tree(b=1,d=0,taxa.stop=20)




Now, let's slice it at slice=0.6, shown as the dotted line above (having verified, of course, that this will indeed produce some non-trivial subtrees):

> source("treeSlice.R")
> trees<-treeSlice(tree,0.6)
> plot(trees,root.edge=T)


which, in this case, produces the following set of three trees:



We're done!

10 comments:

  1. Note that this is *not* in phytools 0.0-3.

    ReplyDelete
  2. Cool!

    If it comes up, I've written code for doing the opposite: taking a slice and obtaining the tree up to that slice.

    ReplyDelete
  3. Hi Liam,

    Thanks a lot. This is really helpful function and I am going to use it very often from now on.

    best
    Rama

    ReplyDelete
  4. @Rama (and others)

    Can you share how you will be using this function? I have some ideas but wanted to hear other ideas.

    Laurence

    ReplyDelete
  5. A neat function, will be using it now to cut a tree of sequence motifs. I'm wondering though whether it's also possible to make sure that single-tip trees are also returned?
    Thanks!
    Mikhail

    ReplyDelete
  6. Done. I will post about this shortly.

    ReplyDelete
  7. This comment has been removed by the author.

    ReplyDelete
  8. This comment has been removed by the author.

    ReplyDelete
  9. Thank you very much for this post. It works.

    ReplyDelete
  10. Hi I tried your tree slicing program and so far it is giving me different answers. My tree has 255 tips and when I do this:
    trees<-treeSlice(nj_tree_root,0.07)
    I get only one subtree:
    str(trees)
    Class "multiPhylo"
    List of 1
    $ :List of 5
    ..$ edge : int [1:4, 1:2] 4 5 5 4 5 1 2 3
    ..$ edge.length: num [1:4] 0.0358 0.0569 0.0485 0.0163
    ..$ tip.label : chr [1:3] "gene1" "gene2" "gene3"
    ..$ Nnode : int 2
    ..$ root.edge : num 0.0393
    ..- attr(*, "class")= chr "phylo"
    ..- attr(*, "order")= chr "cladewise"
    I am not sure what happens to the rest of the tree. I thought this program will slice the tree at a specific height as shown in the example above and then return different parts of the sliced trees.
    Also, when I tried this:trees<-treeSlice(nj_tree_root,0.007)
    I got 31 trees but the sum of all of these tree tip.labels is way more than total tip.labels in the original tree. Please comment. I can send you all the info you need of my tree.

    ReplyDelete