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.
This is a very cool function
ReplyDeleteThanks. Glad it is helpful.
ReplyDeleteHey Liam,
ReplyDeleteThanks 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)