Simulating drift & selection with PopGen

I just added a new function to simulate the combined effects of natural selection & genetic drift for my GitHub package PopGen.

First, here's what the function actually looks like (though it can also be obtained by updating PopGen, as I'll show below):

drift.selection<-function(p0=0.5,Ne=100,w=c(1,1,1),ngen=400,nrep=10,
    colors=NULL,...){
    if(is.null(colors)) colors<-rainbow(nrep)
    w<-(w/max(w))[3:1]
    gametes<-rep(0,2*Ne)
    gametes[1:round(p0*2*Ne)]<-1
    gametes<-replicate(nrep,gametes,simplify=FALSE)
    p<-lapply(gametes,mean)
    for(i in 1:ngen){
        genotypes<-lapply(gametes,function(x) matrix(sample(x),
            length(x)/2,2))
        fitness<-lapply(genotypes,function(x,w) w[rowSums(x)+1],w=w)
        selected<-lapply(fitness,function(prob,x) 
            cbind(sample(x,prob=prob,replace=TRUE),
            sample(x,prob=prob,replace=TRUE)),x=Ne)
        copy<-replicate(nrep,matrix(sample(1:2,2*Ne,replace=TRUE),
            Ne,2),simplify=FALSE)
        gametes<-mapply(function(g,s,c) c(diag(g[s[,1],][,c[,1]]),
            diag(g[s[,2],][,c[,2]])),
            g=genotypes,s=selected,c=copy,SIMPLIFY=FALSE)
        for(j in 1:nrep) p[[j]][i+1]<-mean(gametes[[j]])
    }
    plot(0:ngen,p[[1]],type="l",col=colors[1],lwd=2,ylim=c(0,1),
        xlab="time (generations)",ylab="f(A)")
    nulo<-mapply(lines,x=replicate(nrep-1,0:ngen,simplify=FALSE),
        y=p[2:nrep],col=colors[2:nrep],lwd=2)
    invisible(p)
}

Now let's load update PopGen, load the package, & use it!

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

library(PopGen)

First, we can compare drift + selection (using the new function) to selection acting alone:

Ne<-200
w<-c(1,0.95,0.9)
drift.selection(p0=0.2,Ne=Ne,w=w,lwd=1,nrep=5,ngen=200)
selection(p0=0.2,w=w,time=200,lwd=6,color=phytools::make.transparent("blue",
    0.1),add=TRUE)

Clearly we can see the effect of drift, but it is not huge. As we increase N_e this effect should diminish. For example:

Ne<-1000
w<-c(1,0.95,0.9)
drift.selection(p0=0.2,Ne=Ne,w=w,lwd=1,nrep=5,ngen=200)
selection(p0=0.2,w=w,time=200,lwd=6,color=phytools::make.transparent("blue",
    0.1),add=TRUE)

A more realistic scenario, however, might be one in which we start the allele frequency of A at 1/(2N_e). Why? Because this will be the frequency of a new mutation should one arise.

Ne<-200
w<-c(1,0.95,0.9)
drift.selection(p0=1/(2*Ne),Ne=Ne,w=w,lwd=1,nrep=20,ngen=200)
selection(p0=1/(2*Ne),w=w,time=200,lwd=6,
    color=phytools::make.transparent("blue",0.1),add=TRUE)

Most likely, relatively few of our (now) 20 replicate simulations will have fixed for the allele A. Why? When a mutation first arises even a relatively large fitness advantage of about 6% in the heterozygote is normally insufficient to prevent it from being lost.

Neat.

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!