Friday, July 15, 2011

New function: map.overlap()

I just posted a new function that comes from a project I am currently working on. This function computes the fractional overlap between the mapped states of two different stochastic mapped trees. I will try to describe the algorithm that I used for this in a subsequent post (since I think it is pretty neat), but first, I thought I'd illustrate use of the function first.

First, load "phytools" and the function from source (direct link to source code for function here):

> require(phytools)
> source("map.overlap.R")

Now, let's read in a pair of simmap trees:

> tree1<-read.simmap(text="((((3:{1,0.2442},4:{1,0.2442}):{1,0.4139},2:{1,0.0207:2,0.6374}):{1,0.3572},5:{1,0.4678:2,0.3717:1,0.053:2,0.1227}):{1,0.4471},1:{1,0.4228:2,0.2277:1,0.0564:2,0.7554});")
> tree2<-read.simmap(text="((((3:{1,0.2442},4:{1,0.2442}):{1,0.4139},2:{1,0.0214:2,0.6367}):{1,0.3572},5:{1,0.6239:2,0.3914}):{1,0.0034:2,0.2051:1,0.0511:2,0.1639:1,0.0237},1:{1,0.963:2,0.4993});")

Below is a plot of these trees created using plotSimmap(), but in which I have overlayed the graphs so that they can be more easily compared. (Remember, the colors of the vertical branches are meaningless.)

Because the rate of evolution in the mapped character trait is low (and, thus, the character doesn't change very much), the overlap between the two different stochastic maps seems high. Let's compute its precise value using my function map.overlap():

> map.overlap(tree1,tree2,tol=1e-3)
[1] 0.7805543

We need to set the tolerance level quite high here because I have rounded the mapped edge lengths so that the Newick string for the tree doesn't fill up this blog window!

Now, for comparison, let's consider two simmap trees where the mapped state changes much more frequently (I won't give the Newick strings because they are far too long):

> map.overlap(tree3,tree4,tol=1e-3)
[1] 0.5369969

The overlap is much lower, just as we'd expect. Cool.

No comments:

Post a Comment