Tuesday, March 6, 2018

Coalescent genealogy plot in R

Nowadays every evolutionary genetics seems to feature the same plot illustrating genealogical coalescence. You know the one - it usually consists of different colored circles or squares connected by lines meant to represent parent-daughter relationships? I thought it would be handy to reproduce this in R. Here's my function which will shortly join my PopGen package on GitHub:

coalescent.plot<-function(n=10,ngen=20,colors=NULL,...){
    if(hasArg(sleep)) sleep<-list(...)$sleep
    else sleep<-0.2
    if(is.null(colors)) colors<-rainbow(n=n)
    popn<-matrix(NA,ngen+1,n)
    parent<-matrix(NA,ngen,n)
    popn[1,]<-1:n
    for(i in 1:ngen){
        parent[i,]<-sort(sample(1:n,replace=TRUE))
        popn[i+1,]<-popn[i,parent[i,]]
    }
    plot.new()
    par(mar=c(2.1,4.1,2.1,1.1))
    plot.window(xlim=c(0.5,n+0.5),ylim=c(ngen,0))
    axis(2)
    title(ylab="time (generations)")
    cx.pt<-2*25/max(n,ngen)
    points(1:n,rep(0,n),bg=colors,pch=21,cex=cx.pt)
    for(i in 1:ngen){
        dev.hold()
        for(j in 1:n){
            lines(c(parent[i,j],j),c(i-1,i),lwd=2,
                col=colors[popn[i+1,j]])
        }
        points(1:n,rep(i-1,n),bg=colors[popn[i,]],pch=21,
            cex=cx.pt)
        points(1:n,rep(i,n),bg=colors[popn[i+1,]],pch=21,
            cex=cx.pt)
        dev.flush()
        Sys.sleep(sleep)
    }
}

OK, let's try it!

library(RColorBrewer)
n<-20
cols<-colorRampPalette(brewer.pal("Accent",n=8))(n)
par(bg="grey",fg="white")
coalescent.plot(n=n,colors=cols,ngen=15)

plot of chunk unnamed-chunk-2

The color scheme is arbitrary & the grey background / white foreground is just for effect - but I like it!

Note that the plot shown is static, but in R it is animated. Here's what I mean:

coalescent.plot(n=n,colors=cols,ngen=20)

(Note that I also made the above animation using ImageMagick, but I also had to modify the code of coalescent.plot to include a call to dev.print for every generation of the simulation.)

That's it.

No comments:

Post a Comment