Plotting trees under a Gaussian / bell curve

I just submitted my Evolution 2016 talk title, and so I'd thought it might be fun to post the code that I used to generate the title slide image. The method is on the effect of consensus/average tree methods on PCMs, so I decided to (metaphorically) plot some trees under a Gaussian/bell curve. Just in case you should ever care to do the same, here is how I did it:

library(phytools)

mt<-phytools:::make.transparent

trees<-rmtree(n=5,N=100)

plot.new()

h<-sapply(trees,function(x) max(nodeHeights(x)))

x0<-10
x<-seq(0,20,by=0.01)
y<-100*exp(-(x-x0)^2/25)
par(mar=rep(0.1,4))
plot(x,y,type="l",axes=FALSE,xlab="",ylab="",col="transparent")

for(i in 1:length(trees)){
    x<-runif(n=1)*20
    y<-Inf
    while(y>100*exp(-((x+h[i]/2)-x0)^2/20)){
        x<-runif(n=1)*20
        y<-sample(1:90,1)
    }
    plotTree(trees[[i]],xlim=c(0,20)-x,ylim=c(0,100),
        tips=setNames(1:Ntip(trees[[i]])+y,
        trees[[i]]$tip.label),ftype="off",add=TRUE,lwd=1,
        color=mt("grey",0.4))
}

x0<-mean(par()$usr[1:2])
x<-seq(par()$usr[1],par()$usr[2],by=0.01)
y<-100*exp(-(x-x0)^2/25)
lines(x,y,lwd=2,col=mt("navy",0.3))

text("Consensus trees & their\nsuitability for macroevolutionary\nanalysis",
    x=mean(par()$usr[1:2]),y=1.3*mean(par()$usr[3:4]),cex=3.4)

text("Liam J. Revell (UMass-Boston)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4]),cex=1.8)
text("Klaus Schliep (UMass-Boston)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4])-3.6*strheight("X"),cex=1.8)
text("Cristián Hernández Ulloa(Univ. de Concepción)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4])-7.2*strheight("X"),cex=1.8)

If we want to export to PDF, we just have to flank the code above with the following lines, and to comment out plot.new(). That even allows us to use alternative fonts (although to do so, we first have to install the extrafonts library and take some other preliminary steps that can easily be found online).

pdf(file="Evolution-title.pdf",font="Calibri",width=10,height=7.5)
## ... our code
dev.off()

That's all there is to it.

Update to cophylo to permit links in assoc for taxa not in the tree

I just added a small update (1, 2) to cophylo that permits links to be specified in assoc for taxa that are absent from the tree.

So, for instance:

library(phytools)
t1

## 
## Phylogenetic tree with 20 tips and 19 internal nodes.
## 
## Tip labels:
##  A, B, C, D, F, G, ...
## 
## Rooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 23 tips and 22 internal nodes.
## 
## Tip labels:
##  t16, t3, t15, t2, t13, t6, ...
## 
## Rooted; includes branch lengths.

assoc

##       LETTERS      
##  [1,] "A"     "t17"
##  [2,] "B"     "t7" 
##  [3,] "C"     "t6" 
##  [4,] "D"     "t26"
##  [5,] "E"     "t21"
##  [6,] "F"     "t9" 
##  [7,] "G"     "t20"
##  [8,] "H"     "t23"
##  [9,] "I"     "t8" 
## [10,] "J"     "t18"
## [11,] "K"     "t12"
## [12,] "L"     "t10"
## [13,] "M"     "t11"
## [14,] "N"     "t1" 
## [15,] "O"     "t5" 
## [16,] "P"     "t4" 
## [17,] "Q"     "t19"
## [18,] "R"     "t24"
## [19,] "S"     "t22"
## [20,] "T"     "t13"
## [21,] "U"     "t2" 
## [22,] "V"     "t16"
## [23,] "W"     "t15"
## [24,] "X"     "t3" 
## [25,] "Y"     "t14"
## [26,] "Z"     "t25"

Both assoc[,1] has species not present in t1, and assoc[,2] has species not present in t2:

setdiff(assoc[,1],t1$tip.label)

## [1] "E" "I" "O" "R" "W" "Z"

setdiff(assoc[,2],t2$tip.label)

## [1] "t21" "t12" "t25"

Nonetheless (and this wouldn't have worked before):

obj<-cophylo(t1,t2,assoc)

## Rotating nodes to optimize matching...
## Done.

plot(obj)

That's it.

The version of phytools with this update can be installed from GitHub as follows:

library(devtools)
install_github("liamrevell/phytools")

Control of font sizes in plotting an object of class "cophylo" with phytools

I received a user request recently about changing the font size in a plotted object of class "cophylo". This is straightforward using the optional argument fsize, which is passed to the phytools to internal function phylogram. This is different from the more standard cex (character expansion factor); however it is kind of a legacy attribute of phytools' plotting functions, which are many.

When I saw his plotted "cophylo" object, however, it was apparent to me that the problem was not merely that he didn't know how to control font size, but that the font sizes of the two plotted phylogenies needed to be able to be controlled independently. This is because it is allowed to plot trees in which the right tree & the left tree have different numbers of tips.

The solution to this is a little bit more complicated than it might seem because the optional arguments are passed internally in plot.cophylo using the three-dot ... argument. This makes it a little bit harder than it otherwise might be because we first have to put the ... arguments in a list, then we can check the length of the list vector fsize (if it exists), then we can supply the first or second of these numeric values to the first or second call of phylogram, respectively, by modifying our argument list values & then calling the function with a list of arguments using the R base function do.call.

Here is a code-chunk to sort of illustrate what that looks like:

obj<-list(...)
## ...
leftArgs<-rightArgs<-obj
if(!is.null(obj$fsize)){
    if(length(obj$fsize)>1){
        leftArgs$fsize<-obj$fsize[1]
        rightArgs$fsize<-obj$fsize[2]
        ## ...
    }
}
x1<-do.call("phylogram",c(list(tree=x$trees[[1]],part=0.4),leftArgs))
## ...  
x2<-do.call("phylogram",c(list(tree=x$trees[[2]],part=0.4,
    direction="leftwards"),rightArgs))

More details of the updates can be seen here and here.

The following is a quick demo to illustrate its application using simulated data & trees:

library(phytools)
t1

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 52 tips and 51 internal nodes.
## 
## Tip labels:
##  t16, t17, t9, t8, t7, t45, ...
## 
## Rooted; includes branch lengths.

head(assoc) ## top of association matrix

##    [,1] [,2] 
## I1 "I"  "t16"
## I2 "I"  "t17"
## A  "A"  "t9" 
## B  "B"  "t8" 
## F  "F"  "t7" 
## G1 "G"  "t45"

obj<-cophylo(t1,t2,assoc)

## Rotating nodes to optimize matching...
## Done.

obj

## Object of class "cophylo" containing:
## 
## (1) 2 (possibly rotated) phylogenetic trees in an object of class "multiPhylo".
## 
## (2) A table of associations between the tips of both trees.

Now if we plot using a single font size, then we either have text that is too small in one tree or too big in the other, for instance:

plot(obj) ## fsize=1

plot(obj,fsize=0.8)

Let's try with two different font sizes:

plot(obj,fsize=c(1.1,0.8))

We can even add a scale bar with the same or different font size:

h1<-max(nodeHeights(t1))
h2<-max(nodeHeights(t2))
plot(obj,fsize=c(1.2,0.7,1),scale.bar=round(0.5*c(h1,h2),2))

Neat.

This version of phytools can be installed from GitHub as follows:

library(devtools)
install_github("liamrevell/phytools")

The trees & associations were simulated as follows:

library(phangorn)
t1<-pbtree(n=26,tip.label=LETTERS)
t2<-pbtree(n=52)
tmp<-ltt(t2,plot=FALSE)
h<-tmp$times[tmp$ltt==26]
nn<-sapply(treeSlice(t2,h,trivial=TRUE),function(x) x$tip.label)
N<-sapply(nn,length)
mm<-mapply(function(x,n) rep(x,n),
    x=untangle(rNNI(t1,moves=10),"read.tree")$tip.label,n=N)
assoc<-cbind(unlist(mm),unlist(nn))