Monday, October 29, 2012

Finding the MRCA for two or multiple species across a set of trees

A user recently* contacted me with a question about finding the MRCA for two or multiple species across a set of trees in an R "multiPhylo" object. (*Unfortunately, "recently" is running about two weeks behind, so this email was actually sent on 10/15. To that user, I say "sorry!") She nearly had it correct, but there was a small error in how she was using sapply. Here is the correct way, by way of demonstration:

> require(phytools)
Loading required package: phytools
Loading required package: ape
...
> # simulation some random trees
> trees<-rmtree(N=100,n=10)
> sp<-c("t1","t2") # for example
> # now find the MRCA of sp across the set of trees
> ca<-sapply(trees,function(x,sp) findMRCA(x,sp),sp=sp)
> ca
  [1] 11 19 12 12 12 14 15 14 16 13 12 12 11 16 19 14 17 11
[19] 11 11 11 11 11 11 13 18 19 11 15 12 11 19 12 12 11 12
[37] 11 12 13 11 11 12 12 11 18 12 11 11 15 16 11 16 12 12
[55] 11 14 11 11 15 17 11 14 16 19 11 13 18 11 16 13 12 11
[73] 12 14 13 11 13 11 11 14 17 11 12 11 11 11 11 12 11 13
[91] 11 19 19 11 13 12 12 11 11 11
> # now verify (for the first 3 trees)
> par(mfrow=c(3,1),mar=c(1.5,1.5,1.5,1.5))
> plot.phylo(trees[[1]],cex=1.5); nodelabels(cex=1.5)
> title("tree 1",cex.main=2)
> plot.phylo(trees[[2]],cex=1.5); nodelabels(cex=1.5)
> title("tree 2",cex.main=2)
> plot.phylo(trees[[3]],cex=1.5); nodelabels(cex=1.5)
> title("tree 3",cex.main=2)


Cool. It works.

1 comment:

  1. Aha now I get it, thankds for the information...

    ReplyDelete