More on visualizing a Brownian motion (random walk) simulation on a phylogeny using phytools

Earlier today I tweeted about a plot I’m working on for an upcoming book chapter showing a phylogeny & Brownian motion simulation (technically, a discrete time random walk), as follows:

I'm working on making an illustrative figure for Brownian motion on a phylogeny using #Rstats #phytools. This is where I'm at on it so far.... 👇 pic.twitter.com/DVGMrjrLZ0

— Liam Revell (@phytools_liam) January 26, 2024

A follower asked to see the code, and (naturally) I’m happy to oblige. (He also asked about applying it to an OU process. I’m going to leave that for a future post.)

First, here is my code – with the only change being that I rearranged the two figure panels from horizontal to vertical so they’d better fit in my blog panel.

library(phytools)

set.seed(77)
tree<-read.tree(text="(C:0.7,(B:0.3,A:0.3):0.4);")
tt<-make.era.map(tree,seq(0,0.7,length.out=200))
tt<-map.to.singleton(tt)
x<-fastBM(tt,a=1,sig2=0.5,internal=TRUE)
par(mfrow=c(2,1))
mapped<-paintBranches(tree,edge=1,"1")
mapped<-paintBranches(mapped,edge=5,"2")
mapped<-paintBranches(mapped,edge=2,"3")
mapped<-paintBranches(mapped,edge=3,"4")
cols<-setNames(c("black",grey(1),grey(0.75),grey(0.25)),
  1:4)
plot(mapped,cols,ftype="i",offset=0.2,fsize=1.2,xlim=c(0,1),
  ylim=c(1,3.2),mar=c(5.1,3.1,2.1,1.1),lwd=2,split.vertical=TRUE,
  outline=TRUE)
mtext("a)",adj=0)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
points(pp$xx[4:5],pp$yy[4:5],pch=21,cex=2,bg="white")
text(1.1*pp$xx[1:3],pp$yy[1:3],round(x[1:3],2),pos=4)
text(pp$xx[4]+0.07,pp$yy[4]-0.05,"root",font=3)
text(pp$xx[5],pp$yy[5]-0.15,"internal",font=3,pos=2)
text(0.35,1.1,expression(paste('t'[A]," = 0.7")),cex=0.8)
text(0.2,2.6,expression(paste('t'[AB]," = 0.4")),cex=0.8)
text(0.55,2.1,expression(paste('t'[B]," = 0.3")),cex=0.8)
text(0.55,3.1,expression(paste('t'[A]," = 0.3")),cex=0.8)
phenogram(tt,x,lwd=4,ftype="i",fsize=1.2,
  spread.labels=FALSE,col="black",las=1,cex.axis=0.8,
  xlim=c(0,1),xlab="",ylab="")
mtext("time",side=1,at=0.35,line=2.5)
clip(-0.1,0.7,min(x),max(x))
grid()
pathA<-c(3,phangorn::Ancestors(tt,3,"all"))
pathB<-c(2,phangorn::Ancestors(tt,2,"all"))
pathAB<-c(findMRCA(tt,c("A","B")),
  phangorn::Ancestors(tt,findMRCA(tt,c("A","B")),"all"))
pathC<-c(1,phangorn::Ancestors(tt,1,"all"))
tt<-paintBranches(tt,intersect(tt$edge[,2],pathC),"1")
tt<-paintBranches(tt,intersect(tt$edge[,2],pathB),"3")
tt<-paintBranches(tt,intersect(tt$edge[,2],pathA),"4")
tt<-paintBranches(tt,intersect(tt$edge[,2],pathAB),"2")
phenogram(tt,x,colors=cols,add=TRUE,col="white",lwd=2)
mtext("b)",adj=0)

There’s a lot to unpack here – but the general idea is as follows. We first take an input tree (tree in this example). Then we “map” a set of regimes on it, evenly spaced from the roots to the tip using phytools::make.era.map, thus turning it into a "simmap" object. Next, we using phytools::map.to.singleton to convert it back to a "phylo" object, but with a whole series of unbranching nodes along each of it’s edges. Then we simulate a continuous character up the tree using phytools::fastBM, but saving the internal node conditions by setting internal=TRUE. Finally, we plot these tip and internal values using phytools::phenogram. Without getting into the details, I colored the different edges in both my tree plot and the phenogram graph so that they could be matched up between the two figure panels. (See if you can figure that out! It makes use of an interesting application of phangorn::Ancestors.) The outline effect is accomplished by overplotting – first in black, and then using my edge colors.

We can use the same principle for a tree containing any number of taxa. For example, here’s a simulation on a tree of salamanders that comes packaged with phytools but that is originally from Highton and Larson (1979).

set.seed(22)
data("salamanders")
h<-max(nodeHeights(salamanders))
mt<-make.era.map(salamanders,seq(0,h,length.out=100))
st<-map.to.singleton(mt)
x<-fastBM(st,internal=TRUE)
phenogram(st,x,ftype="i",fsize=0.7,
  mar=c(5.1,4.1,1.1,1.1),lwd=4)
clip(par()$usr[1],h,min(x),max(x))
grid()
phenogram(st,x,ftype="off",add=TRUE,lwd=2,
  color="white")

OK. This is all I’m going to post about that!