Exercise 12: Plotting methods for phylogenies & comparative data in R

There are a range of different visualization methods in different packages for phylogenies in R; however, for comparative methods, phytools is probably the most powerful. This tutorial demonstrates a range of the functionality for plotting trees in the phytools package.

Due to the incredible flexibility of plotting in R, one of my goals has been to try to make figure production just as reproducible as the rest of our data analysis pipeline. Consequently, I will demonstrate some fairly complex ways of controlling our plots in R. In some cases it may be difficult to keep up with the exercise without copying & pasting some of the lines of code of this exercise.

An important word of caution - in many cases the code below will change the default plotting parameters. The easiest way to refresh these is by running dev.off() until all plotting windows have been closed.

1. Methods for continuous characters

First, we can start by loading our packages:

library(phytools)
packageVersion("phytools")
## [1] '0.6.59'

To demonstrate the simplest continuous character methods, I will use some simulated data in the following files:

  1. tree.tre
  2. X.csv

Let's read our data first:

tree<-read.tree("tree.tre")
x<-as.matrix(read.csv("X.csv",row.names=1))[,1]

The simplest class of visualization method simply plots the data at the tips of the tree. For instance:

dotTree(tree,x,length=10,ftype="i")

plot of chunk unnamed-chunk-3

plotTree.barplot(tree,x)

plot of chunk unnamed-chunk-3

plotTree.barplot(tree,x,args.barplot=list(col=sapply(x,
    function(x) if(x>=0) "blue" else "red"),
    xlim=c(-4,4)))

plot of chunk unnamed-chunk-3

dev.off()
## null device 
##           1

Similarly, we can also plot the data at the tips of the tree for multiple traits simultaneously. We will use the same data files we already downloaded.

X<-read.csv("X.csv",row.names=1)
dotTree(tree,X,standardize=TRUE,length=10)

plot of chunk unnamed-chunk-4

phylo.heatmap(tree,X,standardize=TRUE)

plot of chunk unnamed-chunk-4

dev.off()
## null device 
##           1

Aside from plotting the data at the tips, we can also visualize its reconstructed evolution. A simple method for doing this is called a 'traitgram' in which the phylogeny is projected into phenotype space.

phenogram(tree,x,spread.labels=TRUE,spread.cost=c(1,0))

plot of chunk unnamed-chunk-5

A variant of this function exists where the idea is to visualize the uncertainty around the reconstructed values at internal nodes, as follows:

fancyTree(tree,type="phenogram95",x=x,spread.cost=c(1,0))
## Computing density traitgram...

plot of chunk unnamed-chunk-6

Another method in phytools projects reconstructed states onto internal edges & nodes of the tree using a color gradient. I will also demonstrate a helper function called errorbar.contMap which computes & adds error bars to the tree using the same color gradient as was employed in the plot.

obj<-contMap(tree,x,plot=FALSE)
plot(obj,lwd=7,xlim=c(-0.2,3.6))
errorbar.contMap(obj)

plot of chunk unnamed-chunk-7

dev.off()
## null device 
##           1

This method creates a special object class that can then be replotted without needing to repeat all the computation. We can also replot this object using different parameters. For example, in the following script I will invert & then gray-scale the color gradient of our original plot & then replot these two trees in a facing manner.

obj
## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (-4.964627, 3.299418).
obj<-setMap(obj,invert=TRUE)
grey<-setMap(obj,c("white","black"))
par(mfrow=c(1,2))
plot(obj,lwd=7,ftype="off",xlim=c(-0.2,3.6),legend=2)
plot(grey,lwd=7,direction="leftwards",ftype="off",xlim=c(-0.2,3.6),
    legend=2)

plot of chunk unnamed-chunk-8

dev.off()
## null device 
##           1

Now let's try the same thing, but with some empirical data that we have already seen:

  1. Anolis.tre
  2. anole.data.csv
anole.tree<-read.tree("Anolis.tre")
anole.tree
## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
## 
## Rooted; includes branch lengths.
anole.data<-read.csv("anole.data.csv",header=TRUE,row.names=1)
svl<-setNames(anole.data[,1],rownames(anole.data))
obj<-contMap(anole.tree,svl,plot=FALSE)
obj<-setMap(obj,invert=TRUE)
plot(obj,fsize=c(0.4,1),outline=FALSE,lwd=c(3,7),leg.txt="log(SVL)")

plot of chunk unnamed-chunk-9

This plotting method is fairly flexible, as we've seen. We can also plot our tree in a circular or “fan” style. In the following demo, I will also show how to add raster images of some of the taxa to our plot, which is kind of cool.

It uses the following raster images:

  1. barbatus.png
  2. porcatus.png
  3. cristatellus.png
  4. equestris.png
  5. sagrei.png
plot(obj,type="fan",outline=FALSE,fsize=c(0.6,1),legend=1,
    lwd=c(4,8),xlim=c(-1.8,1.8),leg.txt="log(SVL)")
spp<-c("barbatus","porcatus","cristatellus","equestris","sagrei")
library(png)
w<-0.5
for(i in 1:length(spp)){
    xy<-add.arrow(obj,spp[i],col="transparent",arrl=0.45,lwd=3,hedl=0.1)
    arrl<-if(spp[i]%in%c("sagrei","porcatus")) 0.24 
        else if(spp[i]=="barbatus") 0.2 else 0.34
    add.arrow(obj,spp[i],col="blue",arrl=arrl,lwd=3,hedl=0.05)
    img<-readPNG(source=paste(spp[i],".png",sep=""))
    asp<-dim(img)[1]/dim(img)[2]
    rasterImage(img,xy$x0-w/2,xy$y0-w/2*asp,xy$x0+w/2,xy$y0+w/2*asp)
}

plot of chunk unnamed-chunk-10

dev.off()
## null device 
##           1

Another popular visualization method in phytools is something called a phylomorphospace, which is just a projection of the tree into (normally) a two-dimensional morphospace. For example:

phylomorphospace(tree,X[,1:2],xlab="trait 1",ylab="trait 2")

plot of chunk unnamed-chunk-11

For more than three continuous trait dimensions phytools also has a method that I really like, although I've literally never seen it used in publication. I call this method a 'phylogenetic scatterplot matrix' because it contains unidimensional contMap plots on the diagonal & bivariate phylomorphospaces for each pair of traits in off-diagonal positions. This is what that looks like for five traits:

obj<-fancyTree(tree,type="scattergram",X=X[,1:5])
## Computing multidimensional phylogenetic scatterplot matrix...

plot of chunk unnamed-chunk-12

dev.off()
## null device 
##           1

2. Methods for discrete characters

In addition to these continuous character methods, phytools also has a number of discrete character visualization tools for plotting discrete valued traits & trees.

For instance, phytools contains a range of different methods related to a statistical approach called stochastic character mapping, in which dicrete character histories are sampled from their posterior probability distribution.

In the following part of the exercise, I will map 'ecomorph' onto the tree of Anolis lizards, for which we will use the file ecomorph.csv.

ecomorph<-read.csv("ecomorph.csv",row.names=1)
ecomorph<-setNames(ecomorph[,1],rownames(ecomorph))
cols<-setNames(palette()[1:6],levels(ecomorph))
ecomorph.tree<-make.simmap(drop.tip(anole.tree,
    setdiff(anole.tree$tip.label,names(ecomorph))),ecomorph,
    model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            CG         GB         TC         TG         Tr         Tw
## CG -0.6942418  0.1388484  0.1388484  0.1388484  0.1388484  0.1388484
## GB  0.1388484 -0.6942418  0.1388484  0.1388484  0.1388484  0.1388484
## TC  0.1388484  0.1388484 -0.6942418  0.1388484  0.1388484  0.1388484
## TG  0.1388484  0.1388484  0.1388484 -0.6942418  0.1388484  0.1388484
## Tr  0.1388484  0.1388484  0.1388484  0.1388484 -0.6942418  0.1388484
## Tw  0.1388484  0.1388484  0.1388484  0.1388484  0.1388484 -0.6942418
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##        CG        GB        TC        TG        Tr        Tw 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
plot(ecomorph.tree,cols,type="fan",fsize=0.8,lwd=3,ftype="i")
add.simmap.legend(colors=cols,x=0.9*par()$usr[1],
    y=0.9*par()$usr[4],prompt=FALSE,fsize=0.9)