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)

## 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)

## 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)

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!

Intensive short course on phylogenetic comparative methods in R

We are pleased to announce a new graduate-level intensive short course on the use of R for phylogenetic comparative analysis. The course will be four days in length and will take place at the Hotel Ilha Flata in Ilhabela, Sao Paulo State, Brazil, from the 2nd to the 5th of July, 2015. This course is funded by the National Science Foundation. The course is free of cost; and accommodation at the course venue, as well as breakfast & lunch on all course days, is included for all accepted students. There will be a small number of travel stipends available for qualified students and post-docs. Applicants are welcome from any country, but are especially encouraged from the Latin American region.

Topics covered will include: an introduction to the R environment and programming language, tree manipulation, independent contrasts and phylogenetic generalized least squares, ancestral state reconstruction, models of character evolution, diversification analysis, and community phylogenetic analysis. Course instructors will include Dr. Liam Revell (University of Massachusetts Boston), Dr. Luke Harmon (University of Idaho), and Dr. Mike Alfaro (University of California, Los Angeles). Instruction in the course will be primarily in English, thus all students must have a basic working knowledge of scientific English.

To apply for the course, please submit your CV along with a short (maximum 1 page) description of your research interests, background, and reasons for taking the course. Admission is competitive, and preference will go towards students with background in phylogenetics and a compelling motivation for taking the course. Applications should be submitted by email to ilhabela.phylogenetics.course@gmail.com by May 1st, 2015. Questions can be directed to liam.revell@umb.edu (or posted in the comments section, below).

Converting a phylogeny with node labels to a taxonomy

An R-sig-phylo query asked:

“I wondered if anyone had coded a method for converting a phylogenetic tree with polytomies and node labels into a taxonomy (in some form of data table).”

So far as I know, this has not been done - but nor is it very hard, at least in the highly hypothetical circumstances in which every taxonomic level (i.e., Order, Family, Genus, etc.) that is desired in our table is labeled using node labels and every path from the root to any tip has the same number of labels.

Here's a quick demo, using a balanced tree - although it could equally well be a polytomous tree, or a stochastic tree, so long as the aforementioned rule holds.

library(phytools)
## first here is our tree
tree<-stree(n=64,type="balanced")
tree$node.label<-rep("",tree$Nnode)
tree$node.label[1]<-"Order_1"
tree$node.label[66-Ntip(tree)]<-"Family_1"
tree$node.label[97-Ntip(tree)]<-"Family_2"
tree$node.label[68-Ntip(tree)]<-"Genus_1"
tree$node.label[75-Ntip(tree)]<-"Genus_2"
tree$node.label[83-Ntip(tree)]<-"Genus_3"
tree$node.label[90-Ntip(tree)]<-"Genus_4"
tree$node.label[99-Ntip(tree)]<-"Genus_5"
tree$node.label[106-Ntip(tree)]<-"Genus_6"
tree$node.label[114-Ntip(tree)]<-"Genus_7"
tree$node.label[121-Ntip(tree)]<-"Genus_8"
plotTree(tree,fsize=0.65,ftype="i",lwd=1,xlim=c(-0.06,1.1))
nodelabels(node=which(tree$node.label!="")+Ntip(tree),
    text=tree$node.label[which(tree$node.label!="")],cex=0.7)

Now here is the code to pull out the taxonomy:

getAncestors<-phytools:::getAncestors
foo<-function(tip,tree){
    node<-which(tree$tip.label==tip)
    a<-tree$node.label[getAncestors(tree,node,"all")-Ntip(tree)]
    a<-a[length(a):1]
    c(a[a!=""],tip)
}
T<-t(sapply(tree$tip.label,foo,tree=tree))
T

##     [,1]      [,2]       [,3]      [,4] 
## t1  "Order_1" "Family_1" "Genus_1" "t1" 
## t2  "Order_1" "Family_1" "Genus_1" "t2" 
## t3  "Order_1" "Family_1" "Genus_1" "t3" 
## t4  "Order_1" "Family_1" "Genus_1" "t4" 
## t5  "Order_1" "Family_1" "Genus_1" "t5" 
## t6  "Order_1" "Family_1" "Genus_1" "t6" 
## t7  "Order_1" "Family_1" "Genus_1" "t7" 
## t8  "Order_1" "Family_1" "Genus_1" "t8" 
## t9  "Order_1" "Family_1" "Genus_2" "t9" 
## t10 "Order_1" "Family_1" "Genus_2" "t10"
## t11 "Order_1" "Family_1" "Genus_2" "t11"
## t12 "Order_1" "Family_1" "Genus_2" "t12"
## t13 "Order_1" "Family_1" "Genus_2" "t13"
## t14 "Order_1" "Family_1" "Genus_2" "t14"
## t15 "Order_1" "Family_1" "Genus_2" "t15"
## t16 "Order_1" "Family_1" "Genus_2" "t16"
## t17 "Order_1" "Family_1" "Genus_3" "t17"
## t18 "Order_1" "Family_1" "Genus_3" "t18"
## t19 "Order_1" "Family_1" "Genus_3" "t19"
## t20 "Order_1" "Family_1" "Genus_3" "t20"
## t21 "Order_1" "Family_1" "Genus_3" "t21"
## t22 "Order_1" "Family_1" "Genus_3" "t22"
## t23 "Order_1" "Family_1" "Genus_3" "t23"
## t24 "Order_1" "Family_1" "Genus_3" "t24"
## t25 "Order_1" "Family_1" "Genus_4" "t25"
## t26 "Order_1" "Family_1" "Genus_4" "t26"
## t27 "Order_1" "Family_1" "Genus_4" "t27"
## t28 "Order_1" "Family_1" "Genus_4" "t28"
## t29 "Order_1" "Family_1" "Genus_4" "t29"
## t30 "Order_1" "Family_1" "Genus_4" "t30"
## t31 "Order_1" "Family_1" "Genus_4" "t31"
## t32 "Order_1" "Family_1" "Genus_4" "t32"
## t33 "Order_1" "Family_2" "Genus_5" "t33"
## t34 "Order_1" "Family_2" "Genus_5" "t34"
## t35 "Order_1" "Family_2" "Genus_5" "t35"
## t36 "Order_1" "Family_2" "Genus_5" "t36"
## t37 "Order_1" "Family_2" "Genus_5" "t37"
## t38 "Order_1" "Family_2" "Genus_5" "t38"
## t39 "Order_1" "Family_2" "Genus_5" "t39"
## t40 "Order_1" "Family_2" "Genus_5" "t40"
## t41 "Order_1" "Family_2" "Genus_6" "t41"
## t42 "Order_1" "Family_2" "Genus_6" "t42"
## t43 "Order_1" "Family_2" "Genus_6" "t43"
## t44 "Order_1" "Family_2" "Genus_6" "t44"
## t45 "Order_1" "Family_2" "Genus_6" "t45"
## t46 "Order_1" "Family_2" "Genus_6" "t46"
## t47 "Order_1" "Family_2" "Genus_6" "t47"
## t48 "Order_1" "Family_2" "Genus_6" "t48"
## t49 "Order_1" "Family_2" "Genus_7" "t49"
## t50 "Order_1" "Family_2" "Genus_7" "t50"
## t51 "Order_1" "Family_2" "Genus_7" "t51"
## t52 "Order_1" "Family_2" "Genus_7" "t52"
## t53 "Order_1" "Family_2" "Genus_7" "t53"
## t54 "Order_1" "Family_2" "Genus_7" "t54"
## t55 "Order_1" "Family_2" "Genus_7" "t55"
## t56 "Order_1" "Family_2" "Genus_7" "t56"
## t57 "Order_1" "Family_2" "Genus_8" "t57"
## t58 "Order_1" "Family_2" "Genus_8" "t58"
## t59 "Order_1" "Family_2" "Genus_8" "t59"
## t60 "Order_1" "Family_2" "Genus_8" "t60"
## t61 "Order_1" "Family_2" "Genus_8" "t61"
## t62 "Order_1" "Family_2" "Genus_8" "t62"
## t63 "Order_1" "Family_2" "Genus_8" "t63"
## t64 "Order_1" "Family_2" "Genus_8" "t64"

In this case, that's all there is to it. If the number of labels is different in different paths then we will end up with something much messier. For instance:

tree$node.label[122-Ntip(tree)]<-"Subgenus_1"
tree$node.label[125-Ntip(tree)]<-"Subgenus_2"
plotTree(tree,fsize=0.65,ftype="i",lwd=1,xlim=c(-0.06,1.1))
nodelabels(node=which(tree$node.label!="")+Ntip(tree),
    text=tree$node.label[which(tree$node.label!="")],cex=0.7)

T<-sapply(tree$tip.label,foo,tree=tree)
T

## $t1
## [1] "Order_1"  "Family_1" "Genus_1"  "t1"      
## 
## $t2
## [1] "Order_1"  "Family_1" "Genus_1"  "t2"      
## 
## $t3
## [1] "Order_1"  "Family_1" "Genus_1"  "t3"      
## 
## $t4
## [1] "Order_1"  "Family_1" "Genus_1"  "t4"      
## 
## $t5
## [1] "Order_1"  "Family_1" "Genus_1"  "t5"      
## 
## $t6
## [1] "Order_1"  "Family_1" "Genus_1"  "t6"      
## 
## $t7
## [1] "Order_1"  "Family_1" "Genus_1"  "t7"      
## 
## $t8
## [1] "Order_1"  "Family_1" "Genus_1"  "t8"      
## 
## $t9
## [1] "Order_1"  "Family_1" "Genus_2"  "t9"      
## 
## $t10
## [1] "Order_1"  "Family_1" "Genus_2"  "t10"     
## 
## $t11
## [1] "Order_1"  "Family_1" "Genus_2"  "t11"     
## 
## $t12
## [1] "Order_1"  "Family_1" "Genus_2"  "t12"     
## 
## $t13
## [1] "Order_1"  "Family_1" "Genus_2"  "t13"     
## 
## $t14
## [1] "Order_1"  "Family_1" "Genus_2"  "t14"     
## 
## $t15
## [1] "Order_1"  "Family_1" "Genus_2"  "t15"     
## 
## $t16
## [1] "Order_1"  "Family_1" "Genus_2"  "t16"     
## 
## $t17
## [1] "Order_1"  "Family_1" "Genus_3"  "t17"     
## 
## $t18
## [1] "Order_1"  "Family_1" "Genus_3"  "t18"     
## 
## $t19
## [1] "Order_1"  "Family_1" "Genus_3"  "t19"     
## 
## $t20
## [1] "Order_1"  "Family_1" "Genus_3"  "t20"     
## 
## $t21
## [1] "Order_1"  "Family_1" "Genus_3"  "t21"     
## 
## $t22
## [1] "Order_1"  "Family_1" "Genus_3"  "t22"     
## 
## $t23
## [1] "Order_1"  "Family_1" "Genus_3"  "t23"     
## 
## $t24
## [1] "Order_1"  "Family_1" "Genus_3"  "t24"     
## 
## $t25
## [1] "Order_1"  "Family_1" "Genus_4"  "t25"     
## 
## $t26
## [1] "Order_1"  "Family_1" "Genus_4"  "t26"     
## 
## $t27
## [1] "Order_1"  "Family_1" "Genus_4"  "t27"     
## 
## $t28
## [1] "Order_1"  "Family_1" "Genus_4"  "t28"     
## 
## $t29
## [1] "Order_1"  "Family_1" "Genus_4"  "t29"     
## 
## $t30
## [1] "Order_1"  "Family_1" "Genus_4"  "t30"     
## 
## $t31
## [1] "Order_1"  "Family_1" "Genus_4"  "t31"     
## 
## $t32
## [1] "Order_1"  "Family_1" "Genus_4"  "t32"     
## 
## $t33
## [1] "Order_1"  "Family_2" "Genus_5"  "t33"     
## 
## $t34
## [1] "Order_1"  "Family_2" "Genus_5"  "t34"     
## 
## $t35
## [1] "Order_1"  "Family_2" "Genus_5"  "t35"     
## 
## $t36
## [1] "Order_1"  "Family_2" "Genus_5"  "t36"     
## 
## $t37
## [1] "Order_1"  "Family_2" "Genus_5"  "t37"     
## 
## $t38
## [1] "Order_1"  "Family_2" "Genus_5"  "t38"     
## 
## $t39
## [1] "Order_1"  "Family_2" "Genus_5"  "t39"     
## 
## $t40
## [1] "Order_1"  "Family_2" "Genus_5"  "t40"     
## 
## $t41
## [1] "Order_1"  "Family_2" "Genus_6"  "t41"     
## 
## $t42
## [1] "Order_1"  "Family_2" "Genus_6"  "t42"     
## 
## $t43
## [1] "Order_1"  "Family_2" "Genus_6"  "t43"     
## 
## $t44
## [1] "Order_1"  "Family_2" "Genus_6"  "t44"     
## 
## $t45
## [1] "Order_1"  "Family_2" "Genus_6"  "t45"     
## 
## $t46
## [1] "Order_1"  "Family_2" "Genus_6"  "t46"     
## 
## $t47
## [1] "Order_1"  "Family_2" "Genus_6"  "t47"     
## 
## $t48
## [1] "Order_1"  "Family_2" "Genus_6"  "t48"     
## 
## $t49
## [1] "Order_1"  "Family_2" "Genus_7"  "t49"     
## 
## $t50
## [1] "Order_1"  "Family_2" "Genus_7"  "t50"     
## 
## $t51
## [1] "Order_1"  "Family_2" "Genus_7"  "t51"     
## 
## $t52
## [1] "Order_1"  "Family_2" "Genus_7"  "t52"     
## 
## $t53
## [1] "Order_1"  "Family_2" "Genus_7"  "t53"     
## 
## $t54
## [1] "Order_1"  "Family_2" "Genus_7"  "t54"     
## 
## $t55
## [1] "Order_1"  "Family_2" "Genus_7"  "t55"     
## 
## $t56
## [1] "Order_1"  "Family_2" "Genus_7"  "t56"     
## 
## $t57
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_1" "t57"       
## 
## $t58
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_1" "t58"       
## 
## $t59
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_1" "t59"       
## 
## $t60
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_1" "t60"       
## 
## $t61
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_2" "t61"       
## 
## $t62
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_2" "t62"       
## 
## $t63
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_2" "t63"       
## 
## $t64
## [1] "Order_1"    "Family_2"   "Genus_8"    "Subgenus_2" "t64"

Perhaps we could still post-process this - but it would obviously be more difficult. Note that there is no difficult if genera are at different temporal depths in the tree - so longer as there is the same number of named levels between the root & any tip.

That's all.

Phylogenetic PCA 'biplot' with choices argument

A phytools user recently asked the following:

“I am trying to plot the results of a phylogenetic PCA generated with your package Phytools. I would especially like to plot the third and fourth components using the "choices” option of the biplot() function in R. My feeling is that this option is not implemented while using your phyl.pca function since I get the message "choices" is not a graphical parameter. Am I right? Is there a way I can plot other components than the default ones?“

In fact, this is correct. Here is the function code for biplot.phyl.pca (originally suggested to me by Joan Maspons):

library(phytools)
biplot.phyl.pca

## function (x, ...) 
## {
##     biplot(x$S, x$L, ...)
## }
## <environment: namespace:phytools>

tree<-pbtree(n=26,tip.label=LETTERS)
X<-fastBM(tree,nsim=4)
obj<-phyl.pca(tree,X)
obj

## Phylogenetic pca
## Starndard deviations:
##       PC1       PC2       PC3       PC4 
## 1.0959134 0.9057384 0.8286637 0.6021671 
## Loads:
##              PC1        PC2        PC3          PC4
## [1,]  0.45746491  0.3875983  0.2376398 -0.764212515
## [2,] -0.05002709  0.3353159 -0.9366821 -0.087676734
## [3,] -0.68744007 -0.6538311 -0.1233330 -0.291066912
## [4,] -0.81927353  0.5421294  0.1867339  0.004131128

biplot(obj)

biplot(obj,choices=c(3,4)) ## doesn't work

## Warning in plot.window(...): "choices" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "choices" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "choices" is
## not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "choices" is
## not a graphical parameter

## Warning in box(...): "choices" is not a graphical parameter

## Warning in title(...): "choices" is not a graphical parameter

## Warning in text.default(x, xlabs, cex = cex[1L], col = col[1L], ...):
## "choices" is not a graphical parameter

## Warning in plot.window(...): "choices" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "choices" is not a graphical parameter

## Warning in title(...): "choices" is not a graphical parameter

## Warning in axis(3, col = col[2L], ...): "choices" is not a graphical
## parameter

## Warning in axis(4, col = col[2L], ...): "choices" is not a graphical
## parameter

## Warning in text.default(y, labels = ylabs, cex = cex[2L], col = col[2L], :
## "choices" is not a graphical parameter

It is non-trivial, but fairly straightforward, to modify this to permit use of the argument choices to select PC axes to be plotted. For this I made use of the do.call function:

biplot.phyl.pca<-function(x,...){
    to.do<-list(...)
    if(hasArg(choices)){ 
        choices<-list(...)$choices
        to.do$choices<-NULL
    } else choices<-c(1,2)
    to.do$x<-x$S[,choices]
    to.do$y<-x$L[,choices]
    do.call(biplot,to.do)
}

Let's try it:

biplot(obj) ## standard biplot

biplot(obj,choices=c(3,4))

That's it. Assuming I don't find (and no one points out) an error in this, I will add this to the next version of phytools.