As I mentioned in a prior post I am writing a book chapter on visualization for phylogenetic comparative biology. As I component of this, I discuss the basics of plotting a couple of different types of phylogenetic trees, including some instruction on programming these methods.
For one component of this I wrote a simplified tree plotting function in R. This is what it looks like. Excluding annotation, the function code is less than 20 lines:
simpleTreePlot<-function(tree){
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)
# open & size a new plot
plot.new(); par(mar=rep(0.1,4))
plot.window(xlim=c(0,1.1*max(X)),ylim=c(0,max(y)+1))
# 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)
# 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.1)
}
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)
# open & size a new plot
plot.new(); par(mar=rep(0.1,4))
plot.window(xlim=c(0,1.1*max(X)),ylim=c(0,max(y)+1))
# 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)
# 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.1)
}
Try it out:
> require(phytools)
Loading required package: phytools
> tree<-pbtree(b=1,d=0.2,n=30)
> simpleTreePlot(tree)
Loading required package: phytools
> tree<-pbtree(b=1,d=0.2,n=30)
> simpleTreePlot(tree)
Any suggestions?
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.