Tuesday, March 14, 2023

Adding grids to the nodes of a plotted tree

Here is my solution

``````library(phytools)

## function to compute effective aspect ratio
asp<-function(){
dx<-diff(par()\$usr[1:2])
dy<-diff(par()\$usr[3:4])
asp<-(dy/dx)*(par()\$pin[1]/par()\$pin[2])
asp
}

## function to add grid to a node
box_size<-box_size/asp()
xleft<-x-box_size/2
ybottom<-y-box_size/2*asp()*(nr/nc)
xright<-x+box_size/2
ytop<-y+box_size/2*asp()*(nr/nc)
nn<-max(c(nr,nc))
m<-1
for(j in 0:(nr-1)) for(k in 0:(nc-1)){
xx<-c(xleft+k*box_size/nn,
xleft+(k+1)*box_size/nn)
yy<-c(ytop-j*box_size/nn*asp(),
ytop-(j+1)*box_size/nn*asp())
rect(xx[1],yy[1],xx[2],yy[2],col=cols[m])
m<-m+1
}
}

## slip plot area
par(mfrow=c(1,2))

## plot tree
plotTree(tree,ftype="i",xlim=c(0,0.009),offset=1.5,fsize=0.9)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)

## set colors
cols<-setNames(colorRampPalette(c("white","darkgreen"))(6),4:9)

## pull labels
labs<-c(tree\$tip.label,1:tree\$Nnode+Ntip(tree))

for(i in 1:nrow(X)){
node<-which(labs==rownames(X)[i])
x<-pp\$xx[node]
y<-pp\$yy[node]
}
legend("bottomleft",names(cols),pch=22,pt.bg=cols,cex=1,
pt.cex=2,bty="n")

## repeat for second panel
plotTree(tree,ftype="i",xlim=c(0,0.009),offset=1.5,fsize=0.9)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
cols<-setNames(colorRampPalette(c("white","darkred"))(9),0:8)
labs<-c(tree\$tip.label,1:tree\$Nnode+Ntip(tree))
for(i in 1:nrow(X)){
node<-which(labs==rownames(X)[i])
x<-pp\$xx[node]
y<-pp\$yy[node]
}
legend("bottomleft",names(cols),pch=22,pt.bg=cols,cex=1,
pt.cex=2,bty="n")
``````

In theory, this should scale to any arbitrary tree size or number of grid sells.

Let’s try with a simulated data & set of character traits.

``````tree<-rtree(n=40)
labs<-c(tree\$tip.label,1:tree\$Nnode+Ntip(tree))
Q<-matrix(10,7,7,dimnames=list(0:6,0:6))
diag(Q)<--rowSums(Q)+0.1
tt<-tree
n<-Ntip(tt)
for(i in 1:tree\$Nnode)
tt<-bind.tip(tt,as.character(i+n),0,where=Ntip(tt)+i)
X<-sim.Mk(tt,Q,nsim=16)
cols<-setNames(colorRampPalette(c("white","darkblue"))(7),0:6)
plotTree(tree,ftype="off",xlim=c(0,max(nodeHeights(tree))))
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:nrow(X)){
node<-which(labs==rownames(X)[i])
x<-pp\$xx[node]
y<-pp\$yy[node]