Thursday, January 30, 2014

Function for midpoint rooting

Today Todd Oakley asked:

"Does anyone know of an existing midpoint rooting routine? I am displaying trees and would like to show them as midpoint rooted. I've been using the phangorn package, which does the midpoint rooting perfectly, but it has some dependencies that make it unstable on the linux machines I've been using. When I looked last, phangorn was the only package I could find with midpoint rooting."

Well, midpoint rooting is not theoretically very difficult. All one needs to do is find the longest path between any pair of tips & then locate the root midway along that path. I just posted code for this here (as well as a new phytools build, which can be downloaded & installed from source).

Here is a brief description of some of the tricks that I used:

(1) I used cophenetic.phylo from ape to get all the distances between tips, choose the longest one, and then identify the two species on either end of that path.

(2) I used reroot in phytools to re-root the tree immediately below one of these two tips. I did this so that I could use a new custom function getAncestors (which does the same thing as Ancestors in phangorn, but using no phangorn code) to find the set of all internal nodes ancestral to the other tip subtending the longest path. The new position of the root will be between two of these nodes.

(3) I used dist.nodes in ape to compute the distances between the tip of interest and all the internal nodes I had found in (2).

(4) Finally, I found the two nodes subtending the new root position & re-rooted the tree (using reroot from phytools) in the correct position between those nodes.

That's it. The function has not yet been thoroughly tested, so please give me feedback if it doesn't work as intended. Here's a quick demo:

> library(phytools)
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.89’
> tree<-rtree(n=12)
> mpt1<-midpoint.root(tree)
> plotTree(mpt1)
> require(phangorn)
Loading required package: phangorn
> mpt2<-midpoint(tree)
> plotTree(mpt2)

These trees may not look exactly the same, but they are (just with different rotations of internal nodes):

> all.equal.phylo(mpt1,mpt2)
[1] TRUE

I have not tested this function against phangorn's midpoint, but if history is any indication, midpoint is probably faster & more elegantly programmed. Hopefully for Todd's purposes, this will work.

No comments:

Post a Comment