Tuesday, March 14, 2023

Adding grids to the nodes of a plotted tree

Recently, on Twitter a phytools user asked the following.

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
add_grid<-function(x,y,nr,nc,cols,box_size=0.75){
	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))

## read tree
tree<-read.tree(file="piranga_UCE.tre")

## 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))

## read data 1
X<-read.csv(file="patch_counts_jnd1.csv",row.names=1)

## add grids
for(i in 1:nrow(X)){
	node<-which(labs==rownames(X)[i])
	x<-pp$xx[node]
	y<-pp$yy[node]
	add_grid(x,y,2,3,cols[as.character(X[i,])])		
}
## add legend
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))
X<-read.csv(file="patch_counts_jnd4.csv",row.names=1)
for(i in 1:nrow(X)){
	node<-which(labs==rownames(X)[i])
	x<-pp$xx[node]
	y<-pp$yy[node]
	add_grid(x,y,2,3,cols[as.character(X[i,])])		
}
legend("bottomleft",names(cols),pch=22,pt.bg=cols,cex=1,
	pt.cex=2,bty="n")

plot of chunk unnamed-chunk-1

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)
par(fg="#F0EAD6")
for(i in 1:nrow(X)){
	node<-which(labs==rownames(X)[i])
	x<-pp$xx[node]
	y<-pp$yy[node]
	add_grid(x,y,4,4,cols[as.numeric(X[i,])],box_size=0.9)		
}
par(fg="black")
legend("topright",names(cols),pch=22,pt.bg=cols,cex=1,
	pt.cex=2,bty="n",col="#F0EAD6")

plot of chunk unnamed-chunk-3

No comments:

Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.