Saturday, February 11, 2012

Quicker way to rescale the total length of a tree

In my function pbtree I provide the option to rescale the total length (i.e., height) of the tree to an arbitrary value. This is also available in the geiger function rescaleTree. The only problem with this latter option is that rescaleTree works by dividing tree$edge.length by the maximum value in the diagonal of vcv.phylo(tree). Unfortunately the calculation of vcv.phylo(tree) is not too fast for large trees.

In phytools there is a function nodeHeights which computes in a matrix of dimensions dim(tree$edge) the corresponding heights above the root for all the nodes in tree$edge. That means that to rescale the tree to have height scale we can just do:

tree$edge.length<-
   tree$edge.length/max(nodeHeights(tree)[,2])*scale


This is indeed much faster than rescaleTree.

For instance, on my VAIO laptop:

> require(phytools); require(geiger)
> tree<-pbtree(n=2000)
> # rescale tree to length 100
> system.time(tr1<-rescaleTree(tree,100))
   user  system elapsed
118.43    0.14  123.15
> # now let's submit an alternative version
> rescaleTree<-function(tree,scale){
tree$edge.length<-
   tree$edge.length/max(nodeHeights(tree)[,2])*scale
return(tree)
}
> system.time(tr2<-rescaleTree(tree,100))
   user  system elapsed
   2.49    0.06    2.56
> all.equal.phylo(tr1,tr2)
[1] TRUE


That's it.

No comments:

Post a Comment