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)
