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)
## Loading required package: 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
## -2.5058532 -2.5353028 -1.5198028 -3.3444568 -2.6098775 -2.6830252
## G H I J K L
## -2.7673502 -0.5024255 -0.4028755 -2.1903261 -1.3840878 -0.7009695
## M N O P Q R
## -1.6265088 -1.8593068 -2.7533331 -3.1032355 -3.1114163 -2.2118431
## S T U V W X
## -3.7306722 -2.7624590 -4.5961496 -2.5471489 -2.4448663 -2.7218879
## Y Z
## 0.7788399 1.6670663
## 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 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:
## plot leftward
plot(obj,direction="leftwards")
## plot fan style tree
plot(obj,type="fan",lwd=5) ## for example
Here are some real data using anole.tre and svl.csv:
anole.tree<-read.tree("http://www.phytools.org/Ilhabela2015/data/anole.tre")
## or anole.tree<-read.tree("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("http://www.phytools.org/Ilhabela2015/data/svl.csv",
header=TRUE,row.names=1)
## or svl<-read.csv("svl.csv",header=TRUE,row.names=1)
svl<-as.matrix(svl)[,1]
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.6,1),outline=FALSE)
We can also plot this using a circular or fan tree:
plot(obj,type="fan",outline=FALSE)
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, a relatively simple new 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)
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])
pdf(file="anole.plotTree.wBars.pdf",width=10,height=10)
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])
dev.off()
## png
## 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")
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.
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...
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)
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<-describe.simmap(trees,plot=FALSE)
obj
## 100 trees with a mapped discrete character with states:
## CG, GB, TC, TG, Tr, Tw
##
## trees have 24.11 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.24 0.26 0.22 0.13 0.32 0.58 0.97 1.33 0.35 1.58 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.55 0.32 0.69 1 1.17 3.58 2.02 1 2.07 0.16 0.22
## Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y 0.18 0.15 0.3 0.55 0.83 1.24 0.5 0.43
##
## mean total time spent in each state is:
## CG GB TC TG Tr Tw
## raw 13.2553592 47.1436754 32.3207568 71.1973321 13.12972055 28.6201165
## prop 0.0644506 0.2292234 0.1571509 0.3461778 0.06383972 0.1391576
## total
## raw 205.667
## prop 1.000
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)
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)
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)
## Done simulation(s).
x<-tree$states
x
## A B C D E F G H I J K L M N O P Q R
## "1" "0" "0" "1" "1" "0" "0" "1" "0" "1" "0" "0" "0" "1" "1" "1" "1" "0"
## S T U V W X Y Z
## "0" "0" "0" "1" "1" "1" "0" "1"
## 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.223814 1.223814
## 1 1.223814 -1.223814
## (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
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))
## 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)))
## 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])
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: 190
## objective: 190
## objective: 190
## objective: 188
## objective: 188
## objective: 188
## objective: 188
## objective: 186
## objective: 184
## objective: 184
## objective: 184
## objective: 130
## objective: 130
## objective: 130
## objective: 128
## objective: 128
## objective: 128
## objective: 104
## objective: 98
## objective: 98
## objective: 94
## objective: 92
## objective: 90
## objective: 90
## objective: 88
## 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))
That's it.
Written by Liam J. Revell. Last updated July 5, 2015