Wednesday, April 8, 2015

Finding the closest set of node rotations to a given tip ordering

An R-sig-phylo user asked the following:

“Is there an easy way to get R to automatically rotate the nodes of a phylogeny to match an arbitrary ordering of the tips?…. Say I have a particular taxonomic order, such as: SpeciesA, SpeciesC, SpeciesB…. And I want to rotate the nodes of ((C,B),A) to match it - ie to automatically rotate the nodes to give (A(C,B))”

Well, there are some functions in ape that do something related to this (or perhaps this precisely) - I'm thinking rotateConstr and perhaps cophyloplot, but phytools also has a function, minRotate, used primarily internally in phylo.to.map, which attempts to do this. It does it via a simple, greedy algorithm of performing a pre-order traversal of the tree, rotating each node, and accepting the rotation if it improves the objective function which is the match between the desired order and the realized order.

Remarkably, this seems to be surprisingly effective at finding the original order if a set of random rotations are applied to the nodes of a tree.

So, for example:

library(phytools)
tree<-pbtree(n=26,tip.label=LETTERS)
plotTree(tree)

plot of chunk unnamed-chunk-1

## random set of 100 rotations
nn<-sample(1:tree$Nnode+Ntip(tree),100,replace=TRUE)
for(i in 1:length(nn)) tree<-read.tree(text=write.tree(rotate(tree,nn[i])))
## tree all scrambled up
plotTree(tree)

plot of chunk unnamed-chunk-1

## original order
x<-setNames(1:Ntip(tree),LETTERS)
unscrambled<-minRotate(tree,x)
## objective: 48
## objective: 48
## objective: 48
## objective: 44
## objective: 40
## objective: 40
## objective: 40
## objective: 36
## objective: 36
## objective: 34
## objective: 22
## objective: 20
## objective: 18
## objective: 14
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 8
## objective: 6
## objective: 6
## objective: 2
## objective: 2
## objective: 0
## objective: 0
plotTree(unscrambled)

plot of chunk unnamed-chunk-1

This even seems to work for larger trees:

tree<-pbtree(n=200)
x<-setNames(1:200,tree$tip.label)
## random rotations
nn<-sample(1:tree$Nnode+Ntip(tree),100,replace=TRUE)
for(i in 1:length(nn)) tree<-read.tree(text=write.tree(rotate(tree,nn[i])))
## the objective function going to zero indicated fully
## unscrambled
unscrambled<-minRotate(tree,x)
## objective: 6842
## objective: 6842
## objective: 6034
## objective: 6034
## objective: 6032
## objective: 6032
## objective: 6006
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5948
## objective: 5946
## objective: 5942
## objective: 5942
## objective: 5942
## objective: 5942
## objective: 5872
## objective: 5872
## objective: 5872
## objective: 5872
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5864
## objective: 5860
## objective: 5860
## objective: 5848
## objective: 5848
## objective: 5846
## objective: 5840
## objective: 5840
## objective: 5840
## objective: 5838
## objective: 5836
## objective: 5836
## objective: 5836
## objective: 5836
## objective: 5836
## objective: 5836
## objective: 5836
## objective: 5834
## objective: 5822
## objective: 5820
## objective: 5820
## objective: 5820
## objective: 5818
## objective: 4994
## objective: 4994
## objective: 4994
## objective: 4856
## objective: 4856
## objective: 4854
## objective: 4854
## objective: 4844
## objective: 4844
## objective: 4842
## objective: 4842
## objective: 4842
## objective: 4842
## objective: 4840
## objective: 4836
## objective: 4836
## objective: 4834
## objective: 4834
## objective: 4834
## objective: 4834
## objective: 4834
## objective: 4832
## objective: 4828
## objective: 4828
## objective: 4828
## objective: 4828
## objective: 4826
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4824
## objective: 4820
## objective: 4820
## objective: 4820
## objective: 4820
## objective: 4818
## objective: 4810
## objective: 4810
## objective: 4810
## objective: 4810
## objective: 2498
## objective: 2498
## objective: 2498
## objective: 2498
## objective: 2498
## objective: 2492
## objective: 2492
## objective: 2490
## objective: 2488
## objective: 2486
## objective: 2324
## objective: 2324
## objective: 2324
## objective: 2324
## objective: 2324
## objective: 2320
## objective: 2320
## objective: 2320
## objective: 2320
## objective: 2320
## objective: 2316
## objective: 2316
## objective: 2314
## objective: 2314
## objective: 2314
## objective: 2308
## objective: 2308
## objective: 2306
## objective: 2304
## objective: 2300
## objective: 2300
## objective: 2294
## objective: 2294
## objective: 2294
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 172
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 170
## objective: 168
## objective: 168
## objective: 168
## objective: 164
## objective: 164
## objective: 160
## objective: 160
## objective: 158
## objective: 158
## objective: 158
## objective: 28
## objective: 28
## objective: 28
## objective: 28
## objective: 26
## objective: 24
## objective: 24
## objective: 24
## objective: 24
## objective: 24
## objective: 24
## objective: 24
## objective: 24
## objective: 22
## objective: 20
## objective: 18
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 10
## objective: 6
## objective: 6
## objective: 0
## objective: 0
## objective: 0
## objective: 0

I have no idea whether this will work in general - nor if this strategy will minimize the objective function if a perfect match does not exist. Nonetheless….

That's all!

No comments:

Post a Comment