Tuesday, July 23, 2024

An animated GIF of Brownian evolution on a tree using phytools

Yesterday, I tweeted an animated GIF of Brownian motion evolution on a phylogeny.

Here’s how this GIF was created (modified a bit so that the graphic stands vertically for this blog format, rather than horizontally).

## load phytools
library(phytools)
## set simulation parameters
N<-10 ## number of taxa
nlev<-200 ## number of steps
## simulate tree
tree<-pbtree(n=N,scale=1)
## simulate evolution on the tree 
mapped<-make.era.map(tree,
  limits=seq(0,1,length.out=nlev+1))
x<-fastBM(tt<-map.to.singleton(mapped),
  internal=TRUE)
## compute plot without plotting
plotTree(tt,plot=FALSE,ftype="off")
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## set tolerance
tol<-1e-6
## function to mix colors
mix_colors<-function(cols) colorRampPalette(cols)(3)[2]
## all colors
cols<-setNames(hcl.colors(n=Ntip(tree)),tree$tip.label)
## plotting loop
for(i in c(1:nlev,rep(nlev,20))){
  dev.hold()
  ## subdivide plotting area
  layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
  ## plot tree from 1 to i
  par(mar=c(1.1,3.1,2.1,1.1))
  plot(NA,xlim=c(0,1),ylim=c(range(pp$yy)),bty="n",
    axes=FALSE)
  ii<-intersect(which(pp$xx[pp$edge[,2]]>((i+tol-1)/nlev)),
    which(pp$xx[pp$edge[,2]]<=((i+tol)/nlev)))
  for(j in 1:nrow(pp$edge)){
    if(pp$xx[pp$edge[j,2]]<=(i+tol)/nlev){
      lines(pp$xx[pp$edge[j,]],pp$yy[pp$edge[j,]],type="s")
    }
  }
  ## add points at tips
  for(j in 1:length(ii)){
    daughters<-getDescendants(tt,pp$edge[ii[j],2])
    COL<-mix_colors(cols[daughters[daughters<=Ntip(tree)]])
    points(pp$xx[pp$edge[ii[j],2]],pp$yy[pp$edge[ii[j],2]],
      pch=16,col=COL,cex=1.5)
  }
  ## plot random walk from 1 to i
  par(mar=c(3.1,3.1,1.1,1.1))
  plot(NA,xlim=c(0,1),ylim=range(x),bty="n",las=1)
  for(j in 1:nrow(pp$edge)){
    if(pp$xx[pp$edge[j,2]]<=(i+tol)/nlev){
      lines(pp$xx[pp$edge[j,]],x[pp$edge[j,]])
    }
  }
  ## add points at tip
  for(j in 1:length(ii)){
    daughters<-getDescendants(tt,pp$edge[ii[j],2])
    COL<-mix_colors(cols[daughters[daughters<=Ntip(tree)]])
    points(pp$xx[pp$edge[ii[j],2]],x[pp$edge[ii[j],2]],
      pch=16,col=COL,cex=1.5)
  }
  dev.flush()
  Sys.sleep(0.05) ## sleep R
}