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.
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)
## 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
## -0.009685797 0.509191927 0.619909980 -0.382375084 -0.759424172
## F G H I J
## -1.017865681 -1.423530384 -1.390218699 0.071477307 0.016275700
## K L M N O
## 0.093906234 -0.333070445 -1.584253105 0.197771724 0.809517227
## P Q R S T
## 0.948086563 -0.238809121 0.205425189 -0.165778994 -0.225704617
## U V W X Y
## -0.084803920 1.624216958 1.905642115 1.253396087 1.421246681
## Z
## 1.527713040
## 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 (-1.584253, 1.905642).
## 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 23.99 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.31 0.23 0.26 0.13 0.36 0.63 0.65 1.32 0.4 1.41 1.17
## 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.41 0.29 0.61 0.85 0.98 3.28 1.96 1.08 1.76 0.18 0.17
## Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y 0.29 0.12 0.26 0.61 1.37 1.66 0.72 0.52
##
## mean total time spent in each state is:
## CG GB TC TG Tr Tw
## raw 13.46162668 46.0141794 31.1704542 69.9454383 13.03485493 32.0404070
## prop 0.06545352 0.2237315 0.1515579 0.3400908 0.06337846 0.1557878
## 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:
t1<-rtree(n=26,tip.label=LETTERS)
t2<-rtree(n=26,tip.label=LETTERS)
obj<-cophylo(t1,t2)
## Rotating nodes to optimize matching...
## Done.
plot(obj)
## or non-matching:
t2<-rtree(n=40)
assoc<-cbind(t1$tip.label,sample(t2$tip.label,26,replace=TRUE))
assoc
## [,1] [,2]
## [1,] "C" "t29"
## [2,] "W" "t8"
## [3,] "J" "t12"
## [4,] "K" "t36"
## [5,] "P" "t3"
## [6,] "F" "t29"
## [7,] "U" "t39"
## [8,] "N" "t20"
## [9,] "S" "t26"
## [10,] "A" "t14"
## [11,] "I" "t10"
## [12,] "X" "t9"
## [13,] "R" "t28"
## [14,] "V" "t34"
## [15,] "D" "t14"
## [16,] "E" "t19"
## [17,] "B" "t13"
## [18,] "O" "t4"
## [19,] "T" "t35"
## [20,] "Z" "t29"
## [21,] "G" "t38"
## [22,] "Y" "t29"
## [23,] "M" "t6"
## [24,] "Q" "t7"
## [25,] "L" "t35"
## [26,] "H" "t2"
obj<-cophylo(t1,t2,assoc=assoc)
## Rotating nodes to optimize matching...
## Done.
plot(obj)
Interactive tree manipulation:
obj<-collapseTree(anoletree)
Collapse trees into triangles:
set.seed(10)
tree<-pbtree(n=40,scale=1)
obj<-phylo.toBackbone(tree)
Cool!
That's all I can think of right now!
For this one you will have to install the packages mapdata and png.
library(mapdata)
library(png)
## set seed
set.seed(289)
## shuffle function
shuffle<-function(ind,n){
for(i in 1:n){
## pick 2
ii<-sample(1:length(ind),2)
## swap 'em
x<-ind[ii[1]]
y<-ind[ii[2]]
ind[ii[1]]<-y
ind[ii[2]]<-x
}
ind
}
## shadow text function
shadowtext<-function(x,y,labels,...,
offset=0.05,shadow.col="grey",col="black"){
if(length(offset)==1) offset<-rep(offset,2)
if(hasArg(cex)) cex<-list(...)$cex else cex<-1
asp<-par()$din[2]/par()$din[1]
offset.x<-offset[1]*strwidth("W")*cex
offset.y<-offset[2]*strheight("W")*cex
text(x-offset.x,y+offset.y,labels,...,col=shadow.col)
text(x,y,labels,...,col=col)
}
## read in city coordinates
coords<-read.csv(
file="http://www.phytools.org/SanJuan2016/data/Latin-American-cities.csv",
row.names=1)
coords<-coords[1:12,] ## first twelve
coords<-as.matrix(coords)
rownames(coords)<-gsub(" ","_",rownames(coords))
coords<-coords[order(coords[,1]),]
## stretch the tree a little bit to make it look better
tree<-phytools:::lambdaTree(pbtree(n=nrow(coords),
tip.label=rownames(coords),scale=1),0.8)
tree$tip.label<-shuffle(tree$tip.label,1)
## create an object (first pass)
obj<-phylo.to.map(tree,coords,database="worldHires",plot=FALSE,
regions=c("Argentina","Bolivia","Brazil","Chile","Colombia",
"Costa Rica","Cuba","Dominican Republic","Ecuador","El Salvador",
"French Guiana","Guatemala","Haiti","Honduras",
"Mexico","Nicaragua","Panama","Paraguay","Peru",
"Puerto Rico","Saint Bartelemy","Saint Martin","Uruguay",
"Venezuela",
"Belize","Guyana","Suriname","Jamaica"),rotate=FALSE)
## get elements that correspond with the main island of PR
ii<-grep("Puerto Rico",obj$map$names)-1
jj<-which(is.na(obj$map$x))
pr.x<-obj$map$x[(jj[ii[1]]+1):(jj[ii[1]+1]-1)]
pr.y<-obj$map$y[(jj[ii[1]]+1):(jj[ii[1]+1]-1)]
## move San Juan (we're going to make the island big
coords["San_Juan",]<-20*(coords["San_Juan",]-c(mean(pr.y),mean(pr.x)))+
c(28,-50)
## create phylo.to.map object to plot
obj<-phylo.to.map(tree,coords,database="worldHires",plot=FALSE,
regions=c("Argentina","Bolivia","Brazil","Chile","Colombia",
"Costa Rica","Cuba","Dominican Republic","Ecuador","El Salvador",
"French Guiana","Guatemala","Haiti","Honduras",
"Mexico","Nicaragua","Panama","Paraguay","Peru",
"Puerto Rico","Saint Bartelemy","Saint Martin","Uruguay",
"Venezuela",
"Belize","Guyana","Suriname","Jamaica"),rotate=FALSE)
## layout for our heading
layout(matrix(c(1,2),2,1),heights=c(0.16,0.84))
par(mar=rep(0,4))
## create header
plot.new()
plot.window(xlim=c(-1,1),ylim=c(-1,1))
shadowtext(x=0,y=0.2,
"Latin American Macroevolution Workshop\nSan Juan, Puerto Rico",
cex=1.8)
shadowtext(x=0,y=-0.9,"June 28 - July 1, 2016",cex=1)
## plot phylo.to.map object
par(lwd=2)
plot(obj,direction="rightwards",ftype="off",lwd=2,split=c(0.17,0.83),
xlim=c(-120,-36),colors="blue",pch=21)
ii<-grep("Puerto Rico",obj$map$names)-1
jj<-which(is.na(obj$map$x))
pr.x<-obj$map$x[(jj[ii[1]]+1):(jj[ii[1]+1]-1)]
pr.y<-obj$map$y[(jj[ii[1]]+1):(jj[ii[1]+1]-1)]
## this is for the blowout box
box1.x<-c(-67.20739,-67.17,-65.62316,-65.92952)
box1.y<-c(17.95173,18.47613,18.38505,17.96829)
box2.x<-20*(box1.x-mean(pr.x))-50
box2.y<-20*(box1.y-mean(pr.y))+28
## lines for PR blowout
for(i in 1:4) lines(c(box1.x[i],box2.x[i]),c(box1.y[i],box2.y[i]),lwd=1,
lty="dashed")
## blowout PR (although we we actually plot our PNG
mm.x<-mean(pr.x)
mm.y<-mean(pr.y)
pr.x<-20*(pr.x-mm.x)-50
pr.y<-20*(pr.y-mm.y)+28
## read PNG
download.file("http://www.phytools.org/SanJuan2016/data/PR_flag_island.png",
"PR_flag_island.png",mode="wb")
pr.flag<-readPNG(source="PR_flag_island.png")
## add image to map
rasterImage(pr.flag,min(pr.x),min(pr.y),max(pr.x),max(pr.y))
X<-read.csv("http://www.phytools.org/SanJuan2016/data/PR_MAP_Coordinates.csv",
row.names=1)
polygon(19.9*(X[7605:24846,1]-mm.x)-50,
19.9*(X[7605:24846,2]-mm.y)+28,
lwd=2)
## re-plot the cities and the line to San Juan to it overlies PNG
lastPP<-get("last_plot.phylo", envir = .PlotPhyloEnv)
lines(c(lastPP$xx[which(obj$tree$tip.label=="San_Juan")],
coords["San_Juan",2]),
c(lastPP$yy[which(obj$tree$tip.label=="San_Juan")],
coords["San_Juan",1]),lty="dashed",col="blue")
points(coords[,2],coords[,1],pch=21,fg="black",bg="blue",cex=1.1,lwd=1)
points(lastPP$xx[1:Ntip(tree)],lastPP$yy[1:Ntip(tree)],pch=21,fg="black",
bg="blue",cex=0.8,lwd=1)
Written by Liam J. Revell. Last updated 1 Jul. 2016