I have been working on my long overdue 'online practical materials' for the Garamszegi book chapter. My intention is to post code so that readers could reproduce all the figures of the chapter. My chapter is available here. My OPM draft is below. The final OPM will be posted here.
Figures for “Chapter 4: Graphical Methods for Visualizing Comparative Data on Phylogenies”
In the following twelve sections I will describe and provide code that can be used to reproduce the twelve figures of my chapter in Modern Phylogenetic Comparative Methods and their Application in Evolutionary Biology: Concepts and Practice, L´szló Zsolt Garamszegi Ed. Please contact me if you have any questions about this material.
Figure 4.1:
Figure 4.1 is just a simple illustration of two different tree plotting styles that are employed repeatedly throughout the chapter: a right-facing square phylogram; and a circular (fan) phylogram. This exercise requires the file Centrarchidae.tre. Here's how the plots were produced:
library(phytools)
packageVersion("phytools")
## [1] '0.4.38'
text<-"((L:2.29,(K:1.33,(J:1.03,I:1.03):0.3):0.96):0.19,(((H:1.05,(G:0.75,(F:0.73,E:0.73):0.01):0.3):0.71,(D:1.07,(C:0.96,B:0.96):0.11):0.69):0.19,A:1.95):0.53);"
tree<-read.tree(text=text)
rm(text)
centrarchidae<-read.tree(file="Centrarchidae.tre")
par(mfcol=c(1,2))
plotTree(tree,mar=rep(1,4),fsize=1.5)
title(main="(a)",adj=0,cex.main=2)
plotTree(centrarchidae,type="fan",mar=rep(1,4))
## setEnv=TRUE for this type is experimental. please be patient with bugs
title(main="(b)",adj=0,cex.main=2)
Figure 4.2:
Figure 4.2 illustrates, step-by-step, the algorithm for drawing a rightward-facing, square phylogram. Here is the code illustrating how this figure was produced:
## preliminaries
n<-length(tree$tip.label)
# reorder cladewise to assign tip positions
cw<-reorder(tree,"cladewise")
y<-vector(length=n+cw$Nnode)
y[cw$edge[cw$edge[,2]<=n,2]]<-1:n
# reorder pruningwise for post-order traversal
pw<-reorder(tree,"pruningwise")
nn<-unique(pw$edge[,1])
# compute vertical position of each edge
for(i in 1:length(nn)){
yy<-y[pw$edge[which(pw$edge[,1]==nn[i]),2]]
y[nn[i]]<-mean(range(yy))
}
# compute start & end points of each edge
X<-nodeHeights(cw)
## end preliminaries
matrix(c(1,2,3,6,5,4),3,2,byrow=FALSE)->ii
matrix(1:6,3,2,byrow=FALSE)->ii
matrix(1:6,3,2,byrow=TRUE)->ii
layout(ii)
## plot (a)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(a)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")
## NULL
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)
## plot (b)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(b)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")
## NULL
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)
for(i in 2:tree$Nnode+length(tree$tip.label)){
lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),col="gray",
lty="dashed")
text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,i)<=length(cw$tip)]],
collapse=","),col="gray",pos=4,cex=2)
}
## plot (c)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(c)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",
lty="dashed")
## NULL
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)
for(i in 2:tree$Nnode+length(tree$tip.label)){
lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),col="gray",
lty="dashed")
text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,i)<=length(cw$tip)]],
collapse=","),col="gray",pos=4,cex=2)
}
for(i in 1:nrow(X))
points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)
## plot (d)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(d)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")
## NULL
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)
for(i in 2:tree$Nnode+length(tree$tip.label)){
lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
col="gray",lty="dashed")
text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}
# plot horizontal edges
for(i in 1:nrow(X))
lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)
for(i in 1:nrow(X))
points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)
## plot (e)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(e)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")
## NULL
## abline(h=1:max(y),col="gray",lty="dashed")
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray")
for(i in 2:tree$Nnode+length(tree$tip.label)){
lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
col="gray",lty="dashed")
text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}
# plot horizontal edges
for(i in 1:nrow(X))
lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)
# plot vertical relationships
for(i in 1:tree$Nnode+n)
lines(X[which(cw$edge[,1]==i),1],range(y[cw$edge[which(cw$edge[,1]==i),
2]]),lwd=2,lend=2)
for(i in 1:nrow(X))
points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)
## plot (f)
# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(f)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")
## NULL
## abline(h=1:max(y),col="gray",lty="dashed")
for(i in 1:n)
text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)
for(i in 2:tree$Nnode+length(tree$tip.label)){
lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
col="gray",lty="dashed")
text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}
# plot horizontal edges
for(i in 1:nrow(X))
lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)
# plot vertical relationships
for(i in 1:tree$Nnode+n)
lines(X[which(cw$edge[,1]==i),1],range(y[cw$edge[which(cw$edge[,1]==i),
2]]),lwd=2,lend=2)
for(i in 1:nrow(X))
points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)
# plot tip labels
for(i in 1:n)
text(X[which(cw$edge[,2]==i),2],y[i],tree$tip.label[i],pos=4,offset=0.3,
font=2,cex=2)
Figure 4.3:
Figure 4.3 shows a single stochastic character map with the aggregate posterior probabilities at all nodes from 100 stochastic character maps overlain. Here is how it was created:
# load anole 'ecomorph' tree
data(anoletree)
# this is a stochastic mapped tree - pull tip states
x<-getStates(anoletree,"tips")
# conduct stochastic mapping
trees<-make.simmap(anoletree,x,model="ER",nsim=100,message=FALSE)
# get the states at ancestral nodes
aa<-describe.simmap(trees,plot=FALSE)
# set colors to be used for mapping
# cols<-setNames(palette()[1:6],sort(unique(x)))
cols<-setNames(c("blue","gold1","chartreuse3","saddlebrown",
"orangered","purple"),sort(unique(x)))
# plotSimmap
plotSimmap(anoletree,cols,type="fan",lwd=3,fsize=0.8,ftype="i")
## setEnv=TRUE for this type is experimental. please be patient with bugs
# add node labels
par(fg="transparent")
nodelabels(pie=aa$ace,piecol=cols,cex=0.4)
# add legend
par(fg="black")
add.simmap.legend(colors=cols,fsize=0.8,x=-max(nodeHeights(anoletree)),
y=-max(nodeHeights(anoletree)),prompt=FALSE)
Figure 4.4:
Figure 4.4 has two parts. The first part shows a posterior density map of feeding mode on the tree of elopomorph eels, and the second part shows the variance among stochastic character maps for the same tree & trait. This exercise requires the files elopomorph.tre and elophomorph.csv. Here is how that figure was created:
tree<-read.tree("elopomorph.tre")
x<-as.matrix(read.csv("elopomorph.csv",row.names=1))[,1]
mtrees<-make.simmap(tree,x,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
## biting suction feeding
## biting -19.02 19.02
## suction feeding 19.02 -19.02
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## biting suction feeding
## 0.5 0.5
## Done.
dMap<-densityMap(mtrees,states=c("suction feeding","biting"),plot=FALSE)
## sorry - this might take a while; please be patient
## create 'variance' map
vMap<-dMap
aa<-colnames(vMap$tree$mapped.edge)
aa<-colnames(vMap$tree$mapped.edge)
tt<-vMap$tree
# collapse each prob with the same variance into each other
for(i in 0:499){
oo<-c(as.character(i),as.character(1000-i))
nn<-as.character(i)
chk<-oo%in%aa
if(sum(chk)==2) tt<-mergeMappedStates(tt,oo,nn)
else if(chk[2])
for(j in 1:length(tt$maps))
names(tt$maps[[j]])[which(names(tt$maps[[j]])==oo[2])]<-oo[1]
}
# compute the variances
p<-as.numeric(names(vMap$cols))/(length(vMap$cols)-1)
v<-(p*(1-p))[1:501]
vMap$tree<-tt
# this is a hack
vMap$cols<-setNames(as.vector(rbind(grey(v/0.25), grey(v/0.25))),
as.vector(rbind(0:500,0:500)))
vMap$lims<-c(0,0.25)
class(vMap)<-"contMap"
par(mfrow=c(1,2))
plot(dMap,fsize=c(0.7,0.9),outline=TRUE,legend=0.05,mar=c(0.1,0.5,2.1,0.1))
title(main="(a)",adj=0)
par(col="white")
plotTree(tree,fsize=0.7,lwd=3+2,offset=0.2+0.2/3,ftype="i",ylim=c(-6.32,62),
mar = c(0.1,0.5,2.1,0.1))
par(col="black")
plotSimmap(vMap$tree,vMap$cols,lwd=3,fsize=0.7,ftype="i",add=TRUE,
ylim=c(-6.32,62),mar=c(0.1,0.5,2.1,0.1))
add.color.bar(0.05,vMap$cols,title="variance",lims=as.numeric(vMap$lims),
digits=2,x=0,y=1-0.08*(62-1),lwd=3,fsize=0.9,prompt=FALSE)
title(main="(b)",adj=0)
Figure 4.5:
Figure 4.5 is a three-part illustration of one & two-dimensional traitgrams with simulated data. Since I didn't save the seed, results will vary. Here is how it was created:
tree<-pbtree(n=26)
tree$tip.label<-LETTERS
X<-sim.corrs(tree,matrix(c(1,0.7,0.7,1),2,2))
par(mfcol=c(1,3))
par(mar=c(4.1,4.1,4.1,2.1))
phenogram(tree,X[,1],spread.labels=TRUE)
title(main="(a)",adj=0,cex=3)
fancyTree(tree,"traitgram3d",X=X,method="static",angle=-30)
title(main="(b)",adj=0,cex=3)
par(mar=c(4.1,4.1,4.1,2.1))
fancyTree(tree,"phenogram95",x=X[,1],link=0.1*max(nodeHeights(tree)))
title(main="(c)",adj=0,cex=3)
Figure 4.6:
Figure 4.6 is a projection of the phylogeny into a two dimensional morphospace (aka. a phylomorphospace) with the mapping of a multistate discrete ecological trait overlain. The tree and data are for Caribbean Anolis lizards. This exercise requires the file anole.data.csv. Here is how it was created:
X<-read.csv("anole.data.csv",header=T,row.names=1)
X<-as.matrix(X)
rownames(X)<-paste("Anolis_",rownames(X),sep="")
data(anoletree)
X<-X[anoletree$tip.label,]
X<-cbind(X[,1],phyl.resid(anoletree,X[,1],X[,2:20])$resid,phyl.resid(anoletree,
X[,21],X[,22])$resid)
PC<-phyl.pca(anoletree,X)$S
PC<-cbind(-PC[,1],PC[,2])
colnames(PC)<-c("PC1 (limb length)","PC2 (size)")
par(mar=c(5.1,4.1,2.1,2.1))
colors<-setNames(c("blue","gold1","chartreuse3","saddlebrown","orangered","purple"),
sort(unique(getStates(anoletree,"tips"))))
phylomorphospace(anoletree,PC[,c(1,2)],label="off",node.by.map=TRUE,colors=colors)
add.simmap.legend(colors=colors,prompt=FALSE,x=-1.5,y=1.1)
Figure 4.7:
Figure 4.7 is a stochastic phylogeny and projection of that phylogeny into bivariate morphospace with a continuous color gradient showing time since the root overlain. The purpose of this visualization is to also convey temporal information in a phylomorphospace projection. This plot was created as follows (although I neglected to save the seed in the manuscript version, so this results will differ from the published one):
# simulate a tree
set.seed(1)
tree<-pbtree(n=26,scale=100)
tree$tip.label<-LETTERS[26:1]
# simulate data for our phylomorphospace
X<-fastBM(tree,nsim=2)
colnames(X)<-paste("V",1:2,sep="")
# now add the era map
tree<-make.era.map(tree,0:99)
# choose the color scale
colors<-rainbow(100,start=0,end=0.7)
colors<-colorRampPalette(c("blue", "red"))(100)
colors<-colorRampPalette(c("red", "orange", "blue"),space="rgb")(100)
names(colors)<-1:100
par(mfcol=c(1,2))
# check it
plotSimmap(tree,colors,pts=FALSE,lwd=6,mar=c(5.1,2.1,4.1,1.1),ftype="i")
title(main="(a)",adj=0)
# plot phylomorphospace
par(mar=c(5.1,4.1,4.1,1.1))
phylomorphospace(tree,X,colors=colors,node.by.map=TRUE,lwd=4,node.size=c(0.8,1.3),
label="horizontal")
add.color.bar(cols=colors,leg=10,title="time since the root",subtitle="",
lims=c(0,max(nodeHeights(tree))),prompt=FALSE,x=-7.5,y=16)
title(main="(b)",adj=0)
Figure 4.8:
Figure 4.8 shows observed and reconstructed body size (snout-to-vent length, SVL) mapped onto a phylogeny of Greater Antillean Anolis lizard species. This exercise also requires the file anole.data.csv.
X<-read.csv("anole.data.csv",header=TRUE,row.names=1)
data(anoletree)
svl<-as.matrix(X)[,1]
names(svl)<-paste("Anolis_",names(svl),sep="")
svl<-svl[anoletree$tip.label]
obj<-contMap(anoletree,svl,plot=FALSE)
plot(obj,type="fan",legend=5)
Figure 4.9:
Figure 4.9 shows a simlated four trait “phylogenetic scatterplot matrix” continuous character color mappings of each trait on the diagonal, and bivariate phylomorphospaces in off-diagonal positions. Once again, this result uses simulation so results will vary.
tree<-pbtree(n=26,tip.label=LETTERS[26:1])
V<-matrix(c(1,0.5,0.2,0.0,
0.5,1,0.7,0.2,
0.2,0.7,1,0.1,
0.0,0.2,0.1,1),4,4)
V
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0.5 0.2 0.0
## [2,] 0.5 1.0 0.7 0.2
## [3,] 0.2 0.7 1.0 0.1
## [4,] 0.0 0.2 0.1 1.0
X<-sim.corrs(tree,V)
fancyTree(tree,type="scattergram",X=X)
Figure 4.10:
Figure 4.10 shows in panel (a) a traitgram with a discrete character history overlain. Data were simulated with a high rate of evolution on the red branches. Panel (b) shows a posterior density map from stochastic character mapping projected onto a phylomorphospace plot. In this case, data were simulated with a higher evolutionary correlation on blue than red branches.
set.seed(4)
t1<-pbtree(n=26,scale=1)
t1$tip.label<-LETTERS
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-0:1
t1<-sim.history(t1,Q,anc="0")
## Done simulation(s).
colors<-setNames(c("blue","red"),0:1)
par(lend=2)
x<-sim.rates(t1,c(1,10))
## names absent from sig2: assuming same order as $mapped.edge
par(mfcol=c(1,2))
phenogram(t1,x,colors=colors,spread.labels=TRUE,lwd=3,spread.cost=c(1,0),ftype="i")
title(main="(a)",adj=0)
V<-setNames(
list(matrix(c(1,0.8,0.8,1),2,2),matrix(c(2,0,0,2),2,2)),
0:1)
X<-sim.corrs(t1,V)
colnames(X)<-paste("V",1:2,sep="")
trees<-make.simmap(t1,t1$states,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
## 0 1
## 0 -2.683 2.683
## 1 2.683 -2.683
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## 0 1
## 0.5 0.5
## Done.
xx<-densityMap(trees,plot=FALSE)
## sorry - this might take a while; please be patient
phylomorphospace(xx$tree,X,colors=xx$cols,lwd=3,node.by.map=TRUE,label="horizontal")
title(main="(b)",adj=0)
Figure 4.11:
Figure 4.11 shows two different projections of a tree onto a geographic map using simulated data. The following shows how this was created:
# 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: 216
## objective: 206
## objective: 206
## objective: 206
## objective: 206
## objective: 206
## objective: 202
## objective: 200
## objective: 200
## objective: 198
## objective: 198
## objective: 172
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
layout(matrix(c(1,2),2,1),heights=c(0.61,0.39))
plot(xx,type="phylogram",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
title(main="(a)",adj=0)
plot(xx,type="direct",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
title(main="(b)",adj=0)
Figure 4.12:
Figure 4.12 helps illustrate the structure of an object of class "phylo" in memory in R. Here is how it was
created:
tree<-list()
class(tree)<-"phylo"
tree$edge<-matrix(c(6,7,
7,1,
7,8,
8,9,
9,2,
9,3,
8,4,
6,5),8,2,byrow=TRUE)
tree$Nnode<-4
tree$tip.label<-LETTERS[1:5]
tree$edge.length<-c(4,8,5,1,2,2,3,12)
library(plotrix)
plotTree(tree,setEnv=TRUE,offset=1,fsize=1.5,ftype="i",ylim=c(0.5,5.5))
lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
boxed.labels(x=lastPP$xx,y=lastPP$yy,1:9,bg=grey(level=0.45),cex=1.6,xpad=2.0,
ypad=1.5)
Liam J. Revell, Nov. 26 2014.
Hello, it's very useful. However, I think I've found something in the script of the Figure 4.6. When trying the line "X< X[anoletree$tip.label,]" an error appears: "Error in X[tree$tip.label, ] : subÃndice fuera de los lÃmites".
ReplyDeleteThe data(anoletree) object has changed since this was written. Update to remove the line:
Deleterownames(X)<-paste("Anolis_",rownames(X),sep="")