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.

Speed-up for collapseTree, and dev.hold to prevent plot animation 'blinking'

In the past 24 hrs or so I have done a few things to speed up the new phytools function collapseTree. This function gives an animated, interactive tree plotter through which users can collapse & expand subtrees of the phylogeny to & from their ancestral nodes. Code of this new version is here.

Since the animation requires that the tree be replotted many times in succession, most of the speed-ups I mentioned involves vectorizing various actions, as well as avoiding, where ever possible, the need to recalculate quantities. Specifically, I vectorized the drawing of edges & arcs on the tree (using segments & draw.arc, respectively), and I allow internally used plotting function, fan to accepted precalculated values for the "pruningwise" reordered "phylo" object and for the node heights of the tree, since these must be recomputed every time the same tree is plotted.

I also created a new, internally used function called circles to replace nodelabels from the ape package for the simple task of plotting circles at the nodes & tips of the phylogeny.

Finally, particularly when the tree is very large, I found that the animation would sometimes “blink” when running. This is because there is a discernible gap in time between different parts of the tree and node labels being plotted. To eliminate that phenomenon, I now use the function dev.hold to hold the elements being sent to the plotting device, and dev.flush to flush them.

This actually may turn out to be a handy feature to add to other plotting functions in the phytools package. For instance, compare (on your own machine, of course):

library(phytools)
tree<-pbtree(n=100)
x<-fastBM(tree)
contMap(tree,x,ftype="off") ## appears gradually

dev.hold()

## [1] 0

contMap(tree,x,ftype="off")

dev.flush() ## appears all at once

## [1] 0

Finally, here is another video showing the use of collapseTree:

packageVersion("phytools")

## [1] '0.4.51'

tree<-pbtree(n=150)
## don't run
# pruned<-collapseTree(tree)

Cool - no blinking!