Recently, on Twitter a phytools user asked the following.
Hey #Phylogenetic #R peeps! Does anyone know how to fill in grids at nodes, similar to this figure?
— Luke Campillo (@luke_campillo) March 13, 2023
E.g., I have a df w/ tip/node names in col 1, and 9 discrete traits scored 0-10 in remaining cols. Can I use that to place colored grids on the tree?
Any ideas @phytools_liam? pic.twitter.com/00tOJHrPmW
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")
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")
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.