Wednesday, December 18, 2013

New di2multi function for stochastically mapped trees

I just posted some code the collapses internal branches of zero length (or, more specifically, branches with length shorter than some arbitrarily specified value tol) to polytomies for trees with mapped discrete characters created using (for instance) make.simmap or read.simmap. This is exactly the same as di2multi in the ape package (functionally - I programmed it totally differently, hopefully not at the peril of users); however it works for modified "phylo" objects with a mapped discrete trait.

The code to the function is here, but it uses internal phytools functions - so for it to work, you will probably have to install the latest version of phytools (here).

Here's a quick demo of the function at work:

> require(phytools)
> packageVersion("phytools")
[1] ‘0.3.81’
> ## first create a tree with some polytomies
> N<-26
> tree<-pbtree(n=N,tip.label=LETTERS[26:1],scale=1)
> tree$edge.length[tree$edge.length<0.1&tree$edge[,2]>N]<-0
> ## make the tree ultrametric
> d<-max(vcv(tree))-diag(vcv(tree))
> tree$edge.length[tree$edge[,2]<=N]<- tree$edge.length[tree$edge[,2]<=N]+d
> ## check is ultrametric
> is.ultrametric(tree)
[1] TRUE
> ## check is binary
> is.binary.tree(tree
[1] TRUE
> ## simulate discrete character history on tree
> Q<-2*matrix(c(-1,1,1,-1),2,2)
> tree<-sim.history(tree,Q)
> plotSimmap(tree,colors=setNames(c("blue","red"),1:2), lwd=3)
> ## collapse branches of zero length
> tree<-di2multi.simmap(tree)
> ## check is binary
> is.binary.tree(tree)
[1] FALSE
> plotSimmap(tree,colors=setNames(c("blue","red"),1:2), lwd=3)

The main reason this might be important is because densityMap runs into problems when your tree has internal edges that are very short. This can be addressed by first collapsing all zero/small branches in each of the stochastic mapped trees, and then computing the density-map on the collapses stochastic mapped trees. Code for that would look like the following:

trees<-lapply(trees,di2multi.simmap)
class(trees)<-"multiPhylo"
densityMap(trees)

1 comment:

  1. Hi Liam,

    Is it possible to collapse clades (in hard politomies) using support values as criterion? I could figure out how to do that using branch lenghts (di2multi function). How can I do that with support values? Any help is welcome.
    Cheers

    ReplyDelete

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