Creating a phylomorphospace projection (using phytools) that pretty much looks like it was made with ggplot2

I recently stumbled on the following tweet:

Question for you R geniuses: I have a tree and associated x,y positions for the tips. I'd like to combine the two panels below so the tree is projected into the grid. Ideas? #phylogenetics #R pic.twitter.com/KjJHDzI3VK

— Zach B. Hancock (@zachbbionerd) November 4, 2022

My immediate interpretation was that the user wanted to create a phylomorphospace projection that looked like it had been made using ggplot2.

I'm not so sure that's correct, but we could do that – so why not?

Let's try.

First load phytools and some data.

library(phytools)
data(sunfish.tree)
sunfish.tree

## 
## Phylogenetic tree with 28 tips and 27 internal nodes.
## 
## Tip labels:
##  Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
##  non, pisc
## 
## Rooted; includes branch lengths.

data(sunfish.data)
head(sunfish.data)

##                      feeding.mode gape.width buccal.length
## Acantharchus_pomotis         pisc      0.114        -0.009
## Lepomis_gibbosus              non     -0.133        -0.009
## Lepomis_microlophus           non     -0.151         0.012
## Lepomis_punctatus             non     -0.103        -0.019
## Lepomis_miniatus              non     -0.134         0.001
## Lepomis_auritus               non     -0.222        -0.039

Because the tree is a "simmap" object, we'll next convert it to a simple "phylo" object before plotting. We'll also pull out the two quantitative characters from our input data frame.

tree<-as.phylo(sunfish.tree)
X<-sunfish.data[,2:3]

Now, let's create an ordinary phylomorphospace plot using the default settings.

phylomorphospace(tree,X,ftype="off")

OK, now we're ready to proceed.

## get the last plotted tree coordinates
## just for x & y limits
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## adjust axis label offsets
par(mgp = c(2, 0.4, 0))
## create empty plot
plot(NA,xlim=pp$x.lim,ylim=pp$y.lim,axes=FALSE,
    xlab="",ylab="")
## add grey polygon to simulate ggplot2 background
polygon(x=c(par()$usr[1:2],par()$usr[2:1]),
    y=c(rep(par()$usr[3],2),rep(par()$usr[4],2)),
    col="#dddddd",border=FALSE)
## set x & y tick positions
xticks<-pretty(par()$usr[1:2])
yticks<-pretty(par()$usr[3:4])
## add x-axis
axis(1,las=1,at=xticks,lwd.ticks=2,lwd=0,cex.axis=0.8,
    tcl=-0.25)
## simulate vertical white bars from ggplot2
xx<-seq(xticks[1],xticks[length(xticks)],
    length.out=2*length(xticks)-1)
abline(v=xx,col="white",lwd=rep(c(2,1),6))
## do the same for y-axis
axis(2,at=yticks,las=1,lwd.ticks=2,lwd=0,cex.axis=0.8,
    tcl=-0.25)
yy<-seq(yticks[1],yticks[length(yticks)],
    length.out=2*length(yticks)-1)
abline(h=yy,col="white",lwd=rep(c(2,1),6))
## add phylomorphospace without points
phylomorphospace(tree,X,ftype="off",add=TRUE,
    node.size=0)
## get ggplot2-like colors
hues<-seq(15,375,length=Ntip(tree)+1)
cols<-hcl(h=hues,l=65,c=100)[1:Ntip(tree)]
## add points
points(X[tree$tip.label,],col=cols,pch=16,cex=1.5)
## add x and y axis labels
title(xlab="gape width",ylab="buccal length",cex=0.9)

Finally, let's combine it with our phylogeny, as in the original tweet.

layout(matrix(c(1,2),1,2),width=c(0.45,0.55))
plotTree(tree,ftype="off",mar=c(5.1,1.1,4.1,1.1))
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
hues<-seq(15,375,length=Ntip(tree)+1)
cols<-hcl(h=hues,l=65,c=100)[1:Ntip(tree)]
points(pp$xx[1:Ntip(tree)],pp$yy[1:Ntip(tree)],
    pch=16,cex=1.5,col=cols)
par(mar=c(5.1,4.1,4.1,2.1),mgp=c(2,0.4,0))
phylomorphospace(tree,X,ftype="off",axes=FALSE,bty="n",
    xlab="",ylab="",node.size=0)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
polygon(x=c(par()$usr[1:2],par()$usr[2:1]),
    y=c(rep(par()$usr[3],2),rep(par()$usr[4],2)),
    col="#dddddd",border=FALSE)
xticks<-pretty(par()$usr[1:2])
yticks<-pretty(par()$usr[3:4])
axis(1,las=1,at=xticks,lwd.ticks=2,lwd=0,cex.axis=0.8,
    tcl=-0.25)
xx<-seq(xticks[1],xticks[length(xticks)],
    length.out=2*length(xticks)-1)
abline(v=xx,col="white",lwd=rep(c(2,1),6))
axis(2,at=yticks,las=1,lwd.ticks=2,lwd=0,cex.axis=0.8,
    tcl=-0.25)
yy<-seq(yticks[1],yticks[length(yticks)],
    length.out=2*length(yticks)-1)
abline(h=yy,col="white",lwd=rep(c(2,1),6))
phylomorphospace(tree,X,ftype="off",add=TRUE,
    node.size=0)
points(X[tree$tip.label,],col=cols,pch=16,cex=1.5)
title(xlab="gape width",ylab="buccal length",cex=0.9)