Friday, January 27, 2012

New function to paint sub-trees with different states

I just posted a new function to "paint" subtrees on a phylogeny with different mapped states. I did this in response to a user request and the purpose is to allow users to specify arbitrary mapping for use in functions like brownie.lite, brownieREML, and evol.vcv. The name of this new function is paintSubTree and a direct link to the function code is here.

The way this function works is pretty straightforward. Basically, it takes the tree & the basal node number of a sub-tree on the phylogeny. It then paints all the tipward edges of that node (as well as the stem, if selected) in an arbitrarily specified state. Since the function also retains any existing mapping, it can be used multiple times to paint different clades, and even nested clades, in different states.

I accomplished this by first using a function called getDescendants (basically the same as phangorn:Descendants) to compute all the descendant node and tip numbers from the target node. Then I use this set of node numbers (and the handy base function which) to match the tipward nodes to the $edge matrix and then paint the branches accordingly.

Here is a quick illustration of use. First, let's simulate a small tree, set the tip & tip numbers equal, and then plot the tree with node numbers:

> require(phytools)
> tree<-pbtree(n=15); tree$tip.label<-1:15
> plot(tree,font=1); nodelabels(bg="white")




Now let's paint all the edges from node 20 in a derived state, including the stem leading to node 20:

> tree<-paintSubTree(tree,node=20,state="2",stem=T)
> cols<-c("blue","red","green"); names(cols)<-c(1,2,3)
> plotSimmap(tree,cols,lwd=3,pts=F)




Mission accomplished. Now, let's paint a separate part of the tree, say the clade descending from node 24 (this time excluding the stem), with state 3:

> tree<-paintSubTree(tree,node=24,state="3",stem=F)
> plotSimmap(tree,cols,lwd=3,pts=F)




Finally, we can also paint nested clades as well as only tip edges. If we paint a tip edge only, obviously stem must be set to TRUE. Let's do that:

> tree<-paintSubTree(tree,node=3,state="3",stem=T)
> plotSimmap(tree,cols,lwd=3,pts=F)




Cool.

3 comments:

  1. This is a very cool function

    ReplyDelete
  2. Hey Liam,
    Thanks for putting this function together, it will be very useful to me !!
    I just made a little wrapper for paintSubTree for those who want to paint specific clades on the tree (and then make rates or vcv computations on those clades). It is quick and dirty, so one may wanna add some check-up steps at the beginning of the function.
    It is using the bird orders tree, and comparing the output to the one of phyloch package.
    Seems to work fine.
    Thanks a lot
    Cheers
    Seb

    library(phytools)
    library(phyloch)

    #----------------------------------------------------------------
    #### painting subtrees with phyloch
    data(bird.orders)
    clade1 <- c("Struthioniformes", "Tinamiformes", "Craciformes", "Galliformes", "Anseriformes")
    clade2 <- c("Galbuliformes", "Bucerotiformes", "Upupiformes", "Trogoniformes", "Coraciiformes")
    clade3 <- c("Coliiformes", "Cuculiformes", "Psittaciformes", "Apodiformes", "Trochiliformes",
    "Musophagiformes", "Strigiformes", "Columbiformes", "Gruiformes", "Ciconiiformes", "Passeriformes")
    clade.list <- list(clade1, clade2, clade3)
    color <- c("blue", "red", "green")
    tcol <- tip.color(bird.orders, clade.list, color)
    ecol <- edge.color(bird.orders, clade.list, col = color)
    par(mfrow=c(1,2)) # a second plot will follow
    plot(bird.orders, tip.color = tcol, edge.color = ecol)


    #----------------------------------------------------------
    #### painting subtrees with phytools
    source('paintSubTree.r')
    paintClades <- function(myTree, clade.list){
    root.node <- as.character(length(myTree$tip.label)+1)
    myTreeSimmap <- paintSubTree(myTree, node=root.node, state="1", stem=F)
    for (i in 1:length(clade.list)){
    clade <- clade.list[[i]]
    mrcaClade <- mrca(myTree)[clade, clade]
    nodesClade <- unique(as.vector(mrcaClade))
    nodesClade <- nodesClade[nodesClade>length(myTree$tip.label)]
    allNodesAges <- branching.times(myTree)
    nodesCladeAges <- allNodesAges[as.character(nodesClade)]
    cladeAncestor <- names(which.max(nodesCladeAges))
    myTreeSimmap <-paintSubTree(myTreeSimmap, node=cladeAncestor, state=paste(i+1), stem=T)
    }
    return(myTreeSimmap)
    }

    bird.orders.Simmap <- paintClades(bird.orders, clade.list)
    nb.colors=length(clade.list)+1
    cols=c("black", "blue", "red", "green")
    names(cols)=seq(1:nb.colors)
    plotSimmap(bird.orders.Simmap, cols, lwd=2, pts=F)

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.