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. Remember to close your plotting window between code chunks or you might inherit the plotting parameters of the previous exercise component (in which case your plot may look wrong).
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
library(phytools)
## Loading required package: ape
## Loading required package: maps
## simulate a tree
tree<-pbtree(n=26,tip.label=LETTERS)
plotTree(tree)
## simulate data under Brownian motion
x<-fastBM(tree)
x
## A B C D E F
## 1.1323803 1.3513710 2.2322002 2.8722665 3.2095505 2.8507352
## G H I J K L
## 3.3263623 2.3820773 2.1472923 -0.7320202 -1.0205397 -1.9494462
## M N O P Q R
## -2.3925898 -2.4977560 -1.8346364 -1.8935602 -2.2078588 -1.4278239
## S T U V W X
## 0.4921843 -1.9055746 -1.7026842 -2.9325151 -3.2121574 -2.6609568
## Y Z
## -2.8293651 -2.9252126
## plot a simple traitgram
phenogram(tree,x,spread.labels=TRUE)
## change the spread costs
phenogram(tree,x,spread.labels=TRUE,spread.cost=c(1,0))
We can visualize the states just at the tips of the tree using circles:
dotTree(tree,x)
## or
dotTree(tree,x,standardize=TRUE)
If we have multiple characters we have other options still:
X<-fastBM(tree,nsim=10)
dotTree(tree,X)
## or
dotTree(tree,X,standardize=TRUE)
## or
phylo.heatmap(tree,X)
## or
phylo.heatmap(tree,X,standardize=TRUE)
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))
## Computing density traitgram...
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)
The function contMap returns an object of class "contMap" which we can then
more easily replot using, for instance, different parameters:
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 (-3.212157, 3.326362).
## plot leftward
plot(obj,direction="leftwards")
## plot fan style tree
plot(obj,type="fan",lwd=5) ## for example
We can also change the color map, to invert it, or to set it to a gray scale:
obj<-setMap(obj,invert=TRUE)
plot(obj)
obj<-setMap(obj,colors=c("white","black"))
plot(obj)
Here are some real data using Anolis.tre and svl.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.
svl<-read.csv("svl.csv",header=TRUE,row.names=1)
svl<-setNames(svl[,1],rownames(svl))
svl
## ahli alayoni alfaroi aliniger
## 4.039125 3.815705 3.526655 4.036557
## allisoni allogus altitudinalis alumina
## 4.375390 4.040138 3.842994 3.588941
## alutaceus angusticeps argenteolus argillaceus
## 3.554891 3.788595 3.971307 3.757869
## armouri bahorucoensis baleatus baracoae
## 4.121684 3.827445 5.053056 5.042780
## barahonae barbatus barbouri bartschi
## 5.076958 5.003946 3.663932 4.280547
## bremeri breslini brevirostris caudalis
## 4.113371 4.051111 3.874155 3.911743
## centralis chamaeleonides chlorocyanus christophei
## 3.697941 5.042349 4.275448 3.884652
## clivicola coelestinus confusus cooki
## 3.758726 4.297965 3.938442 4.091535
## cristatellus cupeyalensis cuvieri cyanopleurus
## 4.189820 3.462014 4.875012 3.630161
## cybotes darlingtoni distichus dolichocephalus
## 4.210982 4.302036 3.928796 3.908550
## equestris etheridgei eugenegrahami evermanni
## 5.113994 3.657991 4.128504 4.165605
## fowleri garmani grahami guafe
## 4.288780 4.769473 4.154274 3.877457
## guamuhaya guazuma gundlachi haetianus
## 5.036953 3.763884 4.188105 4.316542
## hendersoni homolechis imias inexpectatus
## 3.859835 4.032806 4.099687 3.537439
## insolitus isolepis jubar krugi
## 3.800471 3.657088 3.952605 3.886500
## lineatopus longitibialis loysiana lucius
## 4.128612 4.242103 3.701240 4.198915
## luteogularis macilentus marcanoi marron
## 5.101085 3.715765 4.079485 3.831810
## mestrei monticola noblei occultus
## 3.987147 3.770613 5.083473 3.663049
## olssoni opalinus ophiolepis oporinus
## 3.793899 3.838376 3.637962 3.845670
## paternus placidus poncensis porcatus
## 3.802961 3.773967 3.820378 4.258991
## porcus pulchellus pumilis quadriocellifer
## 5.038034 3.799022 3.466860 3.901619
## reconditus ricordii rubribarbus sagrei
## 4.482607 5.013963 4.078469 4.067162
## semilineatus sheplani shrevei singularis
## 3.696631 3.682924 3.983003 4.057997
## smallwoodi strahmi stratulus valencienni
## 5.035096 4.274271 3.869881 4.321524
## vanidicus vermiculatus websteri whitemani
## 3.626206 4.802849 3.916546 4.097479
obj<-contMap(anole.tree,svl,fsize=c(0.5,1),outline=FALSE)
We can also plot this using a circular or fan tree:
plot(obj,type="fan",outline=FALSE,fsize=c(0.7,1),legend=0.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()
## png
## 2
The result can be viewed here.
Finally, another relatively simple plotting method in phytools is the function plotTree.wBars. That function pretty much does what it sounds like it does:
plotTree.wBars(anole.tree,exp(svl),type="fan",scale=0.002)
## or with tip labels:
plotTree.wBars(anole.tree,exp(svl),type="fan",scale=0.002,tip.labels=TRUE,
fsize=0.7)
It is not too difficult to combine this with a contMap plot. For example:
obj<-contMap(anole.tree,exp(svl),plot=FALSE)
plotTree.wBars(obj$tree,exp(svl),method="plotSimmap",
tip.labels=TRUE,fsize=0.7,colors=obj$cols,type="fan",scale=0.002)
add.color.bar(1.0,obj$cols,title="trait value",lims=obj$lims,prompt=FALSE,
x=0.9*par()$usr[1],y=0.9*par()$usr[3])
Or, conversely, to color just the bars!
## first, let's invert the color map
obj<-setMap(obj,invert=TRUE)
ii<-sapply(1:Ntip(obj$tree), function(x,tree) which(tree$edge[,2]==x),
tree=obj$tree)
cols<-setNames(obj$cols[sapply(ii,function(x,maps)
names(maps[[x]])[length(maps[[x]])],
maps=obj$tree$maps)],obj$tree$tip.label)
## I'm going to change the scale too
plotTree.wBars(obj$tree,exp(svl),method="plotTree",
tip.labels=TRUE,fsize=0.7,col=cols,type="fan",scale=0.004)
add.color.bar(1.0,obj$cols,title="trait value",lims=obj$lims,prompt=FALSE,
x=0.9*par()$usr[1],y=0.9*par()$usr[3])
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")
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")
Try the animated version on your own computer:
phylomorphospace3d(tree,X)
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)
## Computing multidimensional phylogenetic scatterplot matrix...
Stochastic character maps are easy to plot using phytools. For instance:
data(anoletree)
plot(anoletree,type="fan",fsize=0.8,lwd=3)
## no colors provided. using the following legend:
## CG GB TC TG Tr Tw
## "black" "red" "green3" "blue" "cyan" "magenta"
## 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,fsize=0.9)
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.11570723 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145
## GB 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145
## TC 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145
## TG 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145
## Tr 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145
## Tw 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723
## (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.
obj<-summary(trees,plot=FALSE)
obj
## 100 trees with a mapped discrete character with states:
## CG, GB, TC, TG, Tr, Tw
##
## trees have 24.05 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.22 0.33 0.14 0.1 0.29 0.58 0.84 1.28 0.37 1.26 1.27
## 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.77 0.41 0.83 1.18 1.05 3.13 1.82 0.88 1.48 0.1 0.18
## Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y 0.16 0.18 0.19 0.66 1.36 1.58 0.95 0.46
##
## mean total time spent in each state is:
## CG GB TC TG Tr Tw
## raw 13.0477375 45.811497 33.6124676 68.5566523 12.39801036 32.2405961
## prop 0.0634411 0.222746 0.1634315 0.3333382 0.06028197 0.1567612
## total
## raw 205.667
## prop 1.000
plot(obj,type="fan",fsize=0.8,cex=c(0.5,0.3))
add.simmap.legend(colors=cols,x=0.9*par()$usr[1],
y=0.9*par()$usr[4],prompt=FALSE,fsize=0.9)
We can also overlay the posterior probabilities on a random stochastic map in a similar way:
plot(trees[[1]],type="fan",fsize=0.8)
## no colors provided. using the following legend:
## CG GB TC TG Tr Tw
## "black" "red" "green3" "blue" "cyan" "magenta"
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,fsize=0.9)
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 using the elopomorph data from the ancestral state reconstruction part.
In case you no longer have them, these data are here:
tree<-read.tree("elopomorph.tre")
tree
##
## Phylogenetic tree with 61 tips and 60 internal nodes.
##
## Tip labels:
## Moringua_edwardsi, Kaupichthys_nuchalis, Gorgasia_taiwanensis, Heteroconger_hassi, Venefica_proboscidea, Anguilla_rostrata, ...
##
## Rooted; includes branch lengths.
X<-read.csv("elopomorph.csv",row.names=1)
fmode<-setNames(X[,1],rownames(X))
## stochastic maps
trees<-make.simmap(tree,fmode,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## bite suction
## bite -0.01582783 0.01582783
## suction 0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5 0.5
## Done.
## make densityMap
obj<-densityMap(trees,plot=FALSE)
## sorry - this might take a while; please be patient
plot(obj,lwd=4,outline=TRUE,states=c("suction","bite"),
fsize=c(0.7,0.9))
We can, for instance, overlay a discrete character history (or stochastic map) onto a traitgram or phylomorphospace.
## plot traitgram with single stochastic map
bsize<-setNames(X[,2],rownames(X))
phenogram(trees[[1]],bsize,colors=
setNames(c("blue","red"),c("suction","bite")),
spread.labels=TRUE,spread.cost=c(1,0),fsize=0.7,
ftype="i")
## Optimizing the positions of the tip labels...
## or with our posterior density
phenogram(obj$tree,bsize,colors=obj$cols,spread.labels=TRUE,
spread.cost=c(1,0),fsize=0.7,ftype="i",lwd=3,
ylab="body size")
## Optimizing the positions of the tip labels...
add.color.bar(50,obj$cols,title="PP(state=\"bite\")",subtitle="",
prompt=FALSE,x=0,y=240,fsize=0.8)
## or with our posterior probabilities
obj<-summary(trees)
phenogram(tree,bsize,fsize=0.7,ftype="i",ylab="body size")
## Optimizing the positions of the tip labels...
nodelabels(pie=obj$ace,piecol=setNames(c("red","blue"),colnames(obj$ace)),
cex=0.6)
tiplabels(pie=to.matrix(fmode[tree$tip.label],colnames(obj$ace)),
piecol=setNames(c("red","blue"),colnames(obj$ace)),cex=0.4)
It's also possible to do the same kind of thing with phylomorphospace.
How about projecting a tree onto a geographic map? Here's an example using
phytools function phylo.to.map:
# simulate tree & data
set.seed(1)
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: 126
## objective: 126
## objective: 126
## objective: 106
## objective: 106
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 98
## objective: 98
## objective: 96
## objective: 96
## objective: 94
## objective: 94
## objective: 94
## objective: 84
## objective: 76
## objective: 76
## objective: 72
## objective: 72
## objective: 70
## tree points to map coordinates
plot(xx,type="phylogram",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
## tree projected directly on map
plot(xx,type="direct",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
Co-phylogenetic plotting? Let's try it using data from Lopez-Vaamonde et al. (2001) from MPE:
Pleistodontes<-read.tree(text="(((1,(2,3)100)50,(((((((4,(5,(6,7)100)69)100,(8,(9,10))81)61,11),(12,13)100),(14,15)100),(16,17)96),18)),19);")
Pleistodontes$tip.label<-c("P._froggatti",
"P._blandus",
"P._near_blandus",
"P._proximus",
"P._athysanus",
"P._cuneatus",
"P._astrabocheilus",
"P._macrocainus",
"P._rigisamos",
"P._imperialis",
"P._nigriventris",
"P._xanthocephalus",
"P._greenwoodi",
"P._rieki",
"P._plebejus",
"P._schizodontes",
"P._nitens",
"P._spec._nov._1",
"P._regalis")
## rotate all nodes to match orientation in article
Pleistodontes<-rotateNodes(Pleistodontes,"all")
Sycoscapter<-read.tree(text="(1,((((2,3)100,((4,5),(6,7)78)100),((8,9)50,10)99),(((11,12),13)91,(14,15)77)));")
Sycoscapter$tip.label<-c("S._australis_{macrophylla}",
"S._3_{playpoda}",
"S._5_{cerasicarpa}",
"S._6_{lilliputiana}",
"S._7_{subpuberula}",
"S._1_{brachypoda}",
"S._15_{rubiginosa}",
"S._9_{triradiata}",
"S._10_{crassipes}",
"S._8_{pleurocarpa}",
"S._11_{glandifera}",
"S._12_{xylosycia}",
"S._14_{hesperidiiformis}",
"S._2_{aff._obliqua}",
"S._4_{obliqua}")
Sycoscapter<-rotateNodes(Sycoscapter,"all")
assoc<-cbind(
Pleistodontes$tip.label[c(7,8,5,6,18,1,3,4,10,15,13,16,12,14,19)],
Sycoscapter$tip.label)
assoc
## [,1] [,2]
## [1,] "P._greenwoodi" "S._4_{obliqua}"
## [2,] "P._xanthocephalus" "S._2_{aff._obliqua}"
## [3,] "P._plebejus" "S._14_{hesperidiiformis}"
## [4,] "P._rieki" "S._12_{xylosycia}"
## [5,] "P._blandus" "S._11_{glandifera}"
## [6,] "P._regalis" "S._8_{pleurocarpa}"
## [7,] "P._nitens" "S._10_{crassipes}"
## [8,] "P._schizodontes" "S._9_{triradiata}"
## [9,] "P._imperialis" "S._15_{rubiginosa}"
## [10,] "P._athysanus" "S._1_{brachypoda}"
## [11,] "P._astrabocheilus" "S._7_{subpuberula}"
## [12,] "P._proximus" "S._6_{lilliputiana}"
## [13,] "P._macrocainus" "S._5_{cerasicarpa}"
## [14,] "P._cuneatus" "S._3_{playpoda}"
## [15,] "P._froggatti" "S._australis_{macrophylla}"
obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
print=TRUE)
## Rotating nodes to optimize matching...
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 357.777777777778
## objective: 357.777777777778
## objective: 345.111111111111
## objective: 345.111111111111
## objective: 345.111111111111
## objective: 322.777777777778
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## Done.
plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
link.col=make.transparent("blue",0.25),fsize=0.8)
Interactive tree manipulation:
obj<-collapseTree(anoletree)
Collapse trees into triangles:
set.seed(10)
tree<-pbtree(n=40,scale=1)
obj<-phylo.toBackbone(tree)
Cool!
Finally, here is how I created the logo design for this course:
library(plotrix)
source("http://www.phytools.org/Bariloche2016/data/arctext.R")
set.seed(545)
tree<-pbtree(n=78,scale=0.8)
txt<-write.tree(tree)
txt<-strsplit(txt,";")[[1]]
txt<-paste("((Y:0.9,",txt,":0.1):0.1,Z:1);",sep="")
tree<-read.tree(text=txt)
tree<-phytools:::lambdaTree(tree,0.95)
tree<-make.era.map(tree,c(0,0.4,0.8))
par(bg="darkgrey")
par(ljoin=1,lmitre=30)
par(mar=rep(0,4))
col<-rgb(red=117/255,green=170/255,blue=219/255)
par(fg="transparent")
plot(tree,colors=setNames(rep("black",3),1:3),
lwd=10,type="fan",part=0.5,ftype="off",ylim=c(-0.2,1.2),
lend=1)
par(fg="black")
plot(tree,colors=setNames(c(col,"white",col),1:3),
lwd=8,type="fan",part=0.5,ftype="off",ylim=c(-0.2,1.2),
lend=1,add=TRUE)
text(mean(par()$usr[1:2]),-0.15,"Bariloche, Argentina 2016",
col="black",cex=2.7,font=2)
chars<-strsplit("Latin American Macroevolution Workshop","")[[1]]
stretch<-rep(1.1,length(chars))
arctext("Latin American Macroevolution Workshop",center=c(0,0),radius=1.1,
cex=3.5,font=2,stretch=stretch)
That's it!
Written by Liam J. Revell. Last updated 4 December 2016.