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)

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.

Another function to visualize the evolution of a polygenic quantitative trait in R

Yesterday evening I posted about some new functions for computing multilocus and multiallelic Hardy-Weinberg frequencies, and for showing the phenotypic distribution expected for polygenic quantitative traits determined by many loci.

I also developed a function, phenotype.selection, which shows nicely how selection can result in a change in mean that is well beyond the original range of the phenotype, even in the complete absence of mutation.

The code update corresponding to this new functionality can be viewed here.

To get it, I recommend just installing PopGen from GitHub:

library(devtools)
install_github("liamrevell/PopGen")

Here's an example of how the function might be used:

phenotype.selection(nloci=8,p=rep(0.1,8),ngen=400)

FYI, since it is not possible to embed animations in knitted HTML, I first made the .gif using ImageMagick called from within R as follows:

png(file="ps-%03d.png",width=600,height=600)
phenotype.selection(nloci=8,p=rep(0.1,8),ngen=400)
dev.off()
system("ImageMagick convert -delay 10 -loop 0 *.png 6Mar18-post.gif")
file.remove(list.files(pattern=".png"))

then I embedded it manually.

Some new functions for teaching phenotypic evolution

For the Evolution class I'm teaching this semester I've been working on adding some basic population genetic functions to my GitHub R package PopGen.

Today, the first one is a function to compute Hardy-Weinberg frequencies, but for an arbitrary number of alleles.

First, here is what the function looks like:

hardy.weinberg<-function(p=c(0.5,0.5),alleles=c("A","a"),as.matrix=FALSE){
    if(sum(p)!=1) p<-p/sum(p)
    nalleles<-max(length(p),length(alleles))
    if(nalleles!=length(p)) p<-rep(1/nalleles,nalleles)
    if(nalleles!=length(alleles)) {
        if(length(p)<=26) alleles<-LETTERS[1:length(p)]
        else alleles<-paste("A",1:length(p),sep="")
    }
    p<-setNames(p,alleles)
    HW<-as.matrix(p)%*%t(as.matrix(p))
    if(as.matrix) return(HW)
    else {
        hw<-vector()
        k<-1
        for(i in 1:length(alleles)){
            for(j in i:length(alleles)){
                hw[k]<-if(i==j) HW[i,j] else 2*HW[i,j]
                names(hw)[k]<-paste(alleles[c(i,j)],collapse="")
                k<-k+1
            }
        }
        return(hw)
    }
}

Now let's use it. Obviously, it is straightforward in the biallelic case:

hardy.weinberg(p=c(0.7,0.3))

##   AA   Aa   aa 
## 0.49 0.42 0.09

## or
barplot(hardy.weinberg(p=c(0.7,0.3)))

but it also works nicely for multi-allelic situations in which Hardy-Weinberg frequencies are more laborious to compute:

hw<-hardy.weinberg(p=c(0.4,0.3,0.2,0.1))
hw

##   AA   AB   AC   AD   BB   BC   BD   CC   CD   DD 
## 0.16 0.24 0.16 0.08 0.09 0.12 0.06 0.04 0.04 0.01

barplot(hw)

The second function is designed merely to show that for multilocus genotypes in which allelic substitutions have additive effect, we quickly converge on a (approximately) normal phenotypic distribution for even a modest number of loci. Furthermore, that this is true regardless of our allele frequencies.

The way the function works is by computing all multi-locus genotypes and their HW frequencies, and then plotting the resultant phenotypic trait distribution.

Here's the function:

phenotype.freq<-function(nloci=10,p=NULL,effect=1){
    if(is.null(p)) p<-rep(0.5,nloci)
    genotypes<-t(apply(cbind(p,1-p),1,hardy.weinberg))
    COMBN<-permutations(n=3,r=nloci,set=T,repeats.allowed=T)
    PHEN<-rowSums(COMBN-2)*effect
    FREQ<-rep(1,nrow(COMBN))
    for(i in 1:nrow(COMBN)){
        for(j in 1:nloci)   FREQ[i]<-FREQ[i]*genotypes[j,COMBN[i,j]]
    }
    phen<-unique(PHEN)
    freq<-rep(0,length(phen))
    for(i in 1:length(phen)) freq[i]<-sum(FREQ[which(PHEN==phen[i])])
    plot(phen,freq,type="b",pch=21,bg="grey",
        xlab="phenotypic trait value",ylab="relative frequency",
        cex=1.5,xlim=range(phen)+c(-0.5,0.5)*effect)
    for(i in 1:length(phen))
        rect(phen[i]-0.4*effect,0,phen[i]+0.4*effect,
        freq[i],border="grey",
        col=make.transparent("blue",0.2))
}

and here's what I mean using random allele frequencies at each locus (note that this function uses gtools & phytools):

library(gtools)
library(phytools)
## 1 locus
nloci<-1
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)

## 2 loci
nloci<-2
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)

## 4 loci
nloci<-4
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)

## 8 loci
nloci<-8
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)

## 12 loci
nloci<-12
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)

Neat? I kind of think so.

Note that because the function finds all possible multilocus genotypes for nloci greater than 12 or 14 it is not going to work!