Tuesday, May 8, 2012

Getting the ancestral state estimate from the "same" node on many trees

An R-SIG-phylo user today asked how to extract the ancestral state estimate from the same node on a set of trees. First, we need to define what is meant by the "same" node. Here, I presume that the same node is defined by the set of taxa descended from the that node (so is equivalently the MRCA of a set of tip species). In a set of trees read into R using read.tree there is no guarantee that a node defined in this way will have the same node number in every tree.

I suggested that one could use the phytools function findMRCA, which returns the most recent common ancestor for a set of tip species provided in a vector. Say our set of tip species consists of "t1", "t3", and "t8", to extract the estimated ancestral value for their MRCA on a set of trees, we can just use the following code (trait vector is in x):

> sp<-c("t1","t3","t8")
> target.state<-sapply(trees,function(a,x,sp) ace(x,a)$ace[as.character(findMRCA(a,sp))],x=x,sp=sp)

Simple as that. Of course, we can also might want to check in each case as to whether our tip species represent a monophyletic group. To do that we could just do:

> m<-sapply(trees,is.monophyletic,tips=sp)

And then to pull out only the target node state estimates from trees in which our list of taxa are monophyletic, you can use:

> target.state<-target.state[m]

That's it.

No comments:

Post a Comment

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