A refresher on plotting facing phylogenies with node labels

An R phylogenetics user today contacted me about plotting two identical facing trees, along with reconstructed ancestral states at internal nodes, but with only one set of tip labels between them.

This is easy enough to do, so let's see how.

First our tree & data are simulated - but with realistic looking tip labels.

library(phytools)
library(RColorBrewer)
cw<-untangle(tree,"read.tree") ## make sure tips are in order
layout(matrix(1:3,1,3),widths=c(0.45,0.1,0.45))
plotTree(cw,ftype="off",ylim=c(-1,Ntip(cw)))
obj<-rerootingMethod(cw,x,model="ER")
k<-length(levels(x))
cols<-setNames(brewer.pal(max(c(3,k)),"Set1")[1:k],levels(x))
nodelabels(pie=obj$marginal.anc,piecol=cols,cex=1)
tiplabels(pie=to.matrix(x,levels(x)),piecol=cols,cex=0.8)
legend(x="bottomleft",legend=levels(x),pch=21,pt.cex=2.5,pt.bg=cols,cex=1.5,
    bty="n",horiz=T)
ylim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim
plot.new(); plot.window(xlim=c(-0.1,0.1),ylim=ylim)
text(rep(0,Ntip(cw)),1:Ntip(cw),gsub("_"," ",cw$tip.label),cex=0.8,
    font=3)
plotTree(cw,ftype="off",ylim=c(-1,Ntip(cw)),direction="leftwards")
obj<-rerootingMethod(cw,y,model="ER")
k<-length(levels(y))
cols<-setNames(brewer.pal(max(c(3,k)),"Set2")[1:k],levels(y))
nodelabels(pie=obj$marginal.anc,piecol=cols,cex=1)
tiplabels(pie=to.matrix(y,levels(y)),piecol=cols,cex=0.8)
legend(x="bottomright",legend=levels(y),pch=21,pt.cex=2.5,pt.bg=cols,cex=1.5,
    bty="n",horiz=T)

Of course, some of the plotting parameters - such as the font size, line width, and the sizes of the plotted pies - will need to be adjusted depending on the size of our tree; but hopefully this is a good starting point!

FYI, here's how the tree & data were simulated:

N<-60
tips<-paste(sample(LETTERS,N,replace=T),"._",
    replicate(60,paste(sample(letters,sample(4:12,1)),
    collapse="")),sep="")
tree<-pbtree(n=N,tip.label=tips,scale=1)
Q1<-matrix(c(-1,1,1,-1),2,2,dimnames=list(0:1,0:1))
x<-sim.Mk(tree,Q1)
Q2<-matrix(c(-0.4,0.2,0.2,0.2,-0.4,0.2,0.2,0.2,-0.4),3,3,
    dimnames=list(LETTERS[1:3],LETTERS[1:3]))
y<-sim.Mk(tree,Q2)

Recursive bind.tree function

A phytools user contacted me today about coding a recursive version of ape::bind.tree - that is, a function that would take a list of trees & bind them successively together (based on his code, at the root) until only one tree remained.

A recursive function is nothing more than a function that calls itself internally. These are possible, though (in my experience) not that common, in R.

I'm not sure why we'd want to do this since a simple for or while loop could do the same thing, but I thought I'd help him out anyway.

Here's my recursive function:

bind.trees<-function(trees){
    if(length(trees)==2) return(bind.tree(trees[[1]],trees[[2]]))
    else { 
        trees<-c(bind.tree(trees[[1]],trees[[2]]),
            if(length(trees)>2) trees[3:length(trees)])
        trees<-bind.trees(trees) ## this is the recursive part
        return(trees)
    }
}

We can test it out as follows:

trees<-lapply(sample(2:20,9,replace=TRUE),rtree)
trees<-lapply(trees,function(x){
    x$root.edge<-runif(n=1)
    x})
class(trees)<-"multiPhylo"
print(trees,details=TRUE)

## 9 phylogenetic trees
## tree 1 : 7 tips
## tree 2 : 15 tips
## tree 3 : 12 tips
## tree 4 : 11 tips
## tree 5 : 17 tips
## tree 6 : 10 tips
## tree 7 : 12 tips
## tree 8 : 19 tips
## tree 9 : 3 tips

par(mfrow=c(3,3))
plotTree(trees,lwd=1)

tree<-bind.trees(trees)
plotTree(tree,fsize=0.6,lwd=1)

Or, to better see the structure of all our original trees:

tree$edge.length<-NULL
plotTree(tree,fsize=0.5,type="cladogram",nodes="centered",lwd=1)

That's it.