Plotting phylogenies & comparative data using phytools

In this tutorial I'm going to quickly overview a range of plotting methods for phylogenies & comparative data that are implemented in the phytools package.

Continuous character methods

The first plotting method I'll illustrate is called a 'traitgram' which is a projection of the tree into a space defined by trait values & time since the root. So, for example:

## first load phytools
require(phytools)
## simulate a tree
tree <- pbtree(n = 26, tip.label = LETTERS)
plotTree(tree)
plot of chunk unnamed-chunk-1
## simulate data under Brownian motion
x <- fastBM(tree)
x
##        A        B        C        D        E        F        G        H 
##  2.21950  0.25627  1.94721 -1.85028 -1.08259  0.64343  0.51627  1.29706 
##        I        J        K        L        M        N        O        P 
##  0.95156  1.76106 -1.51259 -1.43161  0.67923  1.23026  0.94695  1.02401 
##        Q        R        S        T        U        V        W        X 
##  0.52726  0.70021  0.51362  0.17891  0.28827  0.66855  0.75961  0.07452 
##        Y        Z 
## -0.19691  0.24588
## plot a simple traitgram
phenogram(tree, x)
plot of chunk unnamed-chunk-2
## tip labels are messy - spread them
phenogram(tree, x, spread.labels = TRUE)
plot of chunk unnamed-chunk-2
## change the spread costs
phenogram(tree, x, spread.labels = TRUE, spread.cost = c(1, 0))
plot of chunk unnamed-chunk-2

We can also plot a traitgram with the uncertainty about ancestral character states visualized using transparent probability density. This is implemented in the phytools function fancyTree. So, for instance:

## plot traitgram with 95% CI
fancyTree(tree, type = "phenogram95", x = x, spread.cost = c(1, 0))
plot of chunk unnamed-chunk-3

Next, we can explore the continuous character mapping function in phytools called contMap. Here is a quick demo using the same data:

## plot contMap
obj <- contMap(tree, x)
plot of chunk unnamed-chunk-4

contMap returns an object of class "contMap" which we can then more easily replot using, for instance, different parameters:

## plot leftward
plot(obj, direction = "leftwards")
plot of chunk unnamed-chunk-5
## plot fan style tree
plot(obj, type = "fan", lwd = 5)  ## for example
plot of chunk unnamed-chunk-5

Here are some real data using anole.tre and svl.csv:

anole.tree <- read.tree("anole.tre")
## or
## anole.tree<-read.tree('http://www.phytools.org/anthrotree/plot/anole.tre')
anole.tree
## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
## 	ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
## 
## Rooted; includes branch lengths.
svl <- read.csv("svl.csv", header = TRUE, row.names = 1)
## of
## svl<-read.csv('http://www.phytools.org/anthrotree/plot/svl.csv',header=TRUE,row.names=1)
svl <- as.matrix(svl)[, 1]
svl
##            ahli         alayoni         alfaroi        aliniger 
##           4.039           3.816           3.527           4.037 
##        allisoni         allogus   altitudinalis         alumina 
##           4.375           4.040           3.843           3.589 
##       alutaceus     angusticeps     argenteolus     argillaceus 
##           3.555           3.789           3.971           3.758 
##         armouri   bahorucoensis        baleatus        baracoae 
##           4.122           3.827           5.053           5.043 
##       barahonae        barbatus        barbouri        bartschi 
##           5.077           5.004           3.664           4.281 
##         bremeri        breslini    brevirostris        caudalis 
##           4.113           4.051           3.874           3.912 
##       centralis  chamaeleonides    chlorocyanus     christophei 
##           3.698           5.042           4.275           3.885 
##       clivicola     coelestinus        confusus           cooki 
##           3.759           4.298           3.938           4.092 
##    cristatellus    cupeyalensis         cuvieri    cyanopleurus 
##           4.190           3.462           4.875           3.630 
##         cybotes     darlingtoni       distichus dolichocephalus 
##           4.211           4.302           3.929           3.909 
##       equestris      etheridgei   eugenegrahami       evermanni 
##           5.114           3.658           4.129           4.166 
##         fowleri         garmani         grahami           guafe 
##           4.289           4.769           4.154           3.877 
##       guamuhaya         guazuma       gundlachi       haetianus 
##           5.037           3.764           4.188           4.317 
##      hendersoni      homolechis           imias    inexpectatus 
##           3.860           4.033           4.100           3.537 
##       insolitus        isolepis           jubar           krugi 
##           3.800           3.657           3.953           3.887 
##      lineatopus   longitibialis        loysiana          lucius 
##           4.129           4.242           3.701           4.199 
##    luteogularis      macilentus        marcanoi          marron 
##           5.101           3.716           4.079           3.832 
##         mestrei       monticola          noblei        occultus 
##           3.987           3.771           5.083           3.663 
##         olssoni        opalinus      ophiolepis        oporinus 
##           3.794           3.838           3.638           3.846 
##        paternus        placidus       poncensis        porcatus 
##           3.803           3.774           3.820           4.259 
##          porcus      pulchellus         pumilis quadriocellifer 
##           5.038           3.799           3.467           3.902 
##      reconditus        ricordii     rubribarbus          sagrei 
##           4.483           5.014           4.078           4.067 
##    semilineatus        sheplani         shrevei      singularis 
##           3.697           3.683           3.983           4.058 
##      smallwoodi         strahmi       stratulus     valencienni 
##           5.035           4.274           3.870           4.322 
##       vanidicus    vermiculatus        websteri       whitemani 
##           3.626           4.803           3.917           4.097
obj <- contMap(anole.tree, svl, fsize = c(0.6, 1), outline = FALSE)
plot of chunk unnamed-chunk-6

We can also plot this using a circular or fan tree:

plot(obj, type = "fan", outline = FALSE)
plot of chunk unnamed-chunk-7

Note that the intense aliasing is just a byproduct of plotting to .png format. To get a high quality image we can plot to .pdf or other formats. For instance:

pdf(file = "anole.contMap.pdf", width = 10, height = 10)
plot(obj, type = "fan", outline = FALSE)
dev.off()
## pdf 
##   2

The result can be viewed here.

Finally, a relatively simple new plotting method in phytools is the function plotTree.wBars. This function is not in the current CRAN release of phytools, so to run it you will need to first download & install phytools from source. That function pretty much does what it sounds like it does:

plotTree.wBars(anole.tree, exp(svl), type = "fan", scale = 0.002)
plot of chunk unnamed-chunk-9

This is also kind of cool:

obj <- contMap(anole.tree, exp(svl), plot = FALSE)
plotTree.wBars(obj$tree, exp(svl), method = "plotSimmap", colors = obj$cols,
    type = "fan", scale = 0.002)
add.color.bar(1, obj$cols, title = "trait value", lims = obj$lims, prompt = FALSE,
    x = 0.9 * par()$usr[1], y = 0.9 * par()$usr[3])
plot of chunk unnamed-chunk-10
pdf(file = "anole.plotTree.wBars.pdf", width = 10, height = 10)
plotTree.wBars(obj$tree, exp(svl), method = "plotSimmap", colors = obj$cols,
    type = "fan", scale = 0.002)
add.color.bar(1, obj$cols, title = "trait value", lims = obj$lims, prompt = FALSE,
    x = 0.9 * par()$usr[1], y = 0.9 * par()$usr[3])
dev.off()
## pdf 
##   2

(A non-aliased version here.)

We can also do two dimensional visualizations of phylogenies in morphospace. The is called a 'phylomorphospace'. E.g.:

X <- fastBM(tree, nsim = 3)  ## simulate 3 characters
phylomorphospace(tree, X[, c(1, 2)], xlab = "trait 1", ylab = "trait 2")
plot of chunk unnamed-chunk-12

With phytools we can do static or animated version of the same. I'll just do a static version here, but animated is the default:

phylomorphospace3d(tree, X, method = "static")
plot of chunk unnamed-chunk-13

Try the animated version on your own.

Finally, it's possible to combine the continuous character mapping and 2D phylomorphospaces using a type of phylogenetic scatterplot. Here's a demo:

fancyTree(tree, type = "scattergram", X = X)
plot of chunk unnamed-chunk-14

Discrete character methods

Stochastic character maps are easy to plot using phytools. For instance:

data(anoletree)
plotSimmap(anoletree, type = "fan")
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## setEnv=TRUE for this type is experimental. please be patient with bugs
## we can use any color scheme, but this is the default
cols <- setNames(palette()[1:length(unique(getStates(anoletree, "tips")))],
    sort(unique(getStates(anoletree, "tips"))))
add.simmap.legend(colors = cols, x = 0.9 * par()$usr[1], y = 0.9 * par()$usr[4],
    prompt = FALSE)
plot of chunk unnamed-chunk-15

This is just one stochastically mapped tree. Let's perform the method of stochastic mapping & plot the posterior probabilities at internal nodes.

x <- getStates(anoletree, "tips")
trees <- make.simmap(anoletree, x, model = "ER", nsim = 100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##          CG       GB       TC       TG       Tr       Tw
## CG -0.11571  0.02314  0.02314  0.02314  0.02314  0.02314
## GB  0.02314 -0.11571  0.02314  0.02314  0.02314  0.02314
## TC  0.02314  0.02314 -0.11571  0.02314  0.02314  0.02314
## TG  0.02314  0.02314  0.02314 -0.11571  0.02314  0.02314
## Tr  0.02314  0.02314  0.02314  0.02314 -0.11571  0.02314
## Tw  0.02314  0.02314  0.02314  0.02314  0.02314 -0.11571
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##     CG     GB     TC     TG     Tr     Tw 
## 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667
## Done.
obj <- describe.simmap(trees, plot = FALSE)
obj
## 100 trees with a mapped discrete character with states:
##  CG, GB, TC, TG, Tr, Tw 
## 
## trees have 25.59 changes between states on average
## 
## changes are of the following types:
##      CG,GB CG,TC CG,TG CG,Tr CG,Tw GB,CG GB,TC GB,TG GB,Tr GB,Tw TC,CG
## x->y  0.25  0.31  0.25   0.2  0.32  0.71  1.02   1.2  0.51  1.46  1.25
##      TC,GB TC,TG TC,Tr TC,Tw TG,CG TG,GB TG,TC TG,Tr TG,Tw Tr,CG Tr,GB
## x->y  0.67  0.36  0.69  1.01  1.13  3.25  1.98  1.11  1.66  0.18  0.34
##      Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y  0.23  0.22  0.45  0.65  1.26  1.49   0.8  0.63
## 
## mean total time spent in each state is:
##            CG      GB      TC      TG       Tr      Tw total
## raw  13.09759 46.1139 32.5794 69.7658 12.91992 31.1904 205.7
## prop  0.06368  0.2242  0.1584  0.3392  0.06282  0.1517   1.0
plot(obj, type = "fan")
## setEnv=TRUE for this type is experimental. please be patient with bugs
add.simmap.legend(colors = cols, x = 0.9 * par()$usr[1], y = 0.9 * par()$usr[4],
    prompt = FALSE)
plot of chunk unnamed-chunk-16

We can also overlay the posterior probabilities on a random stochastic map in a similar way:

plotSimmap(trees[[1]], type = "fan")
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## setEnv=TRUE for this type is experimental. please be patient with bugs
nodelabels(pie = obj$ace, piecol = cols, cex = 0.5)
add.simmap.legend(colors = cols, x = 0.9 * par()$usr[1], y = 0.9 * par()$usr[4],
    prompt = FALSE)
plot of chunk unnamed-chunk-17

For binary characters, we can also map the posterior density of stochastic maps onto the nodes & branches of a tree using the phytools function densityMap. Here's a demo of this using simulated data.

## first let's simulate data
Q <- matrix(c(-1, 1, 1, -1), 2, 2)
rownames(Q) <- colnames(Q) <- c(0, 1)
tree <- sim.history(tree, Q)
x <- tree$states
x
##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "0" "1" "0" "1" "0" "0" "1" "0" "0" "0" "0" "0" "1" "0" "1" "1" "1" "0" 
##   S   T   U   V   W   X   Y   Z 
## "0" "0" "0" "0" "0" "0" "0" "0"
## stochastic maps
trees <- make.simmap(tree, x, nsim = 100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##        0      1
## 0 -1.201  1.201
## 1  1.201 -1.201
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   0   1 
## 0.5 0.5
## Done.
## make densityMap
obj <- densityMap(trees, lwd = 4, outline = TRUE)
## sorry - this might take a while; please be patient
plot of chunk unnamed-chunk-18

Combining discrete & continuous characters

We can, for instance, overlay a discrete character history (or stochastic map) onto a traitgram or phylomorphospace.

## simulate data with a high rate on some branches
x <- sim.rates(tree, setNames(c(1, 10), c(0, 1)))
## plot traitgram
phenogram(tree, x, colors = setNames(c("blue", "red"), c(0, 1)), spread.labels = TRUE,
    spread.cost = c(1, 0))
plot of chunk unnamed-chunk-19
## we can do the same with phylomorphospaces
X <- cbind(sim.rates(tree, setNames(c(1, 10), c(0, 1))), sim.rates(tree, setNames(c(1,
    10), c(0, 1))))
phylomorphospace(tree, X, colors = setNames(c("blue", "red"), c(0, 1)))
plot of chunk unnamed-chunk-19
## it is even possible to overlay a posterior density from stochastic mapping
phylomorphospace(obj$tree, X, colors = obj$cols, lwd = 3, xlab = "x", ylab = "y")
add.color.bar(4, obj$cols, title = "PP(state=1)", prompt = FALSE, x = 0.9 *
    par()$usr[1], y = 0.9 * par()$usr[3])
plot of chunk unnamed-chunk-19

Other stuff

How about projecting a tree onto a geographic map? Here's an example using phytools function phylo.to.map:

# simulate tree & data
tree <- pbtree(n = 26, scale = 100)
tree$tip.label <- LETTERS[26:1]
lat <- fastBM(tree, sig2 = 10, bounds = c(-90, 90))
long <- fastBM(tree, sig2 = 80, bounds = c(-180, 180))
# now plot projection
xx <- phylo.to.map(tree, cbind(lat, long), plot = FALSE)
## objective: 132
## objective: 132
## objective: 132
## objective: 132
## objective: 132
## objective: 128
## objective: 128
## objective: 128
## objective: 128
## objective: 124
## objective: 120
## objective: 118
## objective: 116
## objective: 116
## objective: 108
## objective: 104
## objective: 102
## objective: 102
## objective: 100
## objective: 100
## objective: 94
## objective: 94
## objective: 94
## objective: 94
## objective: 94
## tree points to map coordinates
plot(xx, type = "phylogram", asp = 1.3, mar = c(0.1, 0.5, 3.1, 0.1))
plot of chunk unnamed-chunk-20
## tree projected directly on map
plot(xx, type = "direct", asp = 1.3, mar = c(0.1, 0.5, 3.1, 0.1))
plot of chunk unnamed-chunk-20

That's all for now!

Written by Liam J. Revell. Last updated May 29, 2014