Wednesday, March 13, 2013

Applying countSimmap to a set of trees

I just received the following email request about the function countSimmap that counts the transitions on a mapped tree that's been created or read into memory:

I was wondering if this function (countSimmap) could work on SIMMAP files that contain multiple mappings (like the attached). I got it to work on single mappings like the example on the blog but get errors on my other files.

This is not too hard. Here's how we can do it using sapply:
library(phytools)
## first simulate a tree & some mapped histories
tree<-pbtree(n=100,1)
Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
trees<-sim.history(tree,Q,nsim=10)

## now let's run the analysis using sapply
foo<-function(tree){
  XX<-countSimmap(tree)
  setNames(c(XX$N,as.vector(t(XX$Tr))),c("N",
    sapply(rownames(XX$Tr),paste,colnames(XX$Tr),sep=",")))
}
XX<-t(sapply(trees,foo))

And here are what the results look like for this toy example:
> XX
        N c,c c,a c,b a,c a,a a,b b,c b,a b,b
 [1,] 212   0  38  28  36   0  41  33  36   0
 [2,] 179   0  25  36  28   0  23  31  36   0
 [3,] 188   0  29  31  33   0  31  33  31   0
 [4,] 184   0  36  33  30   0  30  24  31   0
 [5,] 189   0  36  30  27   0  32  28  36   0
 [6,] 196   0  30  31  36   0  36  28  35   0
 [7,] 177   0  24  25  37   0  33  26  32   0
 [8,] 175   0  30  19  42   0  38  19  27   0
 [9,] 184   0  24  27  31   0  31  36  35   0
[10,] 214   0  41  47  36   0  30  32  28   0

That's it.

1 comment:

  1. I just added this automatically within countSimmap if the input object is of class "multiPhylo" (described here).

    ReplyDelete