Diploid Selection Formula

In class we derived an equation for the change in alleles frequencies due to selection. If \(p\) and \(q\) are the frequencies of alleles \(A\) and \(a\) respectively and \(w_{11}\), \(w_{12}\) and \(w_{22}\) are the relative fitnesses of phenotypes AA, Aa and aa, we determine that the change in p is:

\[\Delta p = \frac{p \cdot q}{\bar{w}}[p \cdot(w_{11}-w_{12}) - q \cdot(w_{22}-w_{12})]\] where \[\bar{w} =w_{11} \cdot p^2 + w_{12} \cdot 2pq + w_{22} \cdot q^2 \]

We can use those formulas write an R function to compute the next generation’s frequency of A:

new.p <- function(p, w11, w12, w22){
      wbar <- w11*p^2 + 2*w12*p*(1-p) + w22*(1-p)^2
      p.change <- (p*(1-p)/wbar)*(p*(w11-w12) - (1-p)*(w22-w12))
      return(p + p.change)
}

Let’s try it out:

new.p(.5, 1, 1, 1)
new.p(.5, .8, .8, 1)
new.p(.5, 0, 0, 1)
new.p(.5, 1, 0, 0)

Simulating Many Generations

Of course, there’s no reason to stop after one generation, let’s look 100 generations in the future instead.

w11 <- 0.8; w12 <- 1.0; w22 <- 0.5; p <- 0.5

for (i in 1:100) p <- new.p(p, w11, w12, w22)

print(p)

Try out some different combinations of starting values and see what you see.

Phenotype Frequencies

We’ll also need a function that returns the frequencis of phenotypes AA, Aa and aa given p and, thanks to Hardy-Weinberg, we can do that readily as well:

pheno.freqs <- function(p){
return(c(p^2, 2*p*(1-p), (1-p)^2))
}

And here’s the function in action:

pheno.freqs(0.5)
pheno.freqs(0.9)
pheno.freqs(0.1)

Finally, let’s write a function that plots the phenotype and allele frequencies over time

library(ggplot2)


plot.allele.history <- function(p, w11, w12, w22, gens=100){

history <- data.frame(generation=numeric(), cat=character(), Type=character(), frequency=numeric())

for (i in 1:gens){
  phenos <- data.frame(generation=rep(i,3), cat=rep("phenotypes",3), Type=c("AA", "Aa", "aa"), freq=pheno.freqs(p))
  alleles <- data.frame(generation=rep(i,2), cat=rep("alleles",2), Type=c("p", "q"), freq=c(p, 1-p))
    history <- rbind(history, phenos,alleles)
  p <- new.p(p=p, w11=w11, w12=w12, w22=w22)
}

g <- ggplot(history, aes(generation, freq, color=Type))+geom_line()+ylim(c(0,1))+
  facet_grid(cat~.)

return(g)
}

Can you explain these three histories?

plot.allele.history(p=0.9, w11=0.8, w12=1, w22=0.8)

plot.allele.history(p=0.9, w11=0.9, w12=1, w22=0.8)

plot.allele.history(p=0.1, w11=0.9, w12=1, w22=0.8)

What about this pair of histories?

plot.allele.history(p=0.9, w11=0.8, w12=0.8, w22=1)

plot.allele.history(p=0.9, w11=0.9, w12=0.8, w22=1)

What about this trio, all with the same relative fitnesses and very similar initial p values?

plot.allele.history(p=0.34, w11=1, w12=0.8, w22=0.9)
plot.allele.history(p=0.33, w11=1, w12=0.8, w22=0.9)
plot.allele.history(p=1/3, w11=1, w12=0.8, w22=0.9)

Try out some different sets of starting values and try to make sense of the plots.

Challenges: