## Friday, January 26, 2024

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

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)
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)
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,
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")
``````

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