Changes in allele and genotype frequency

Setting up

For this tutorial, you will need to load the tidyverse package. Like so:

library(tidyverse)

Plotting the expected genotype frequencies

The first thing we need to do is generate a range of allele frequencies to visualize the expected genotype frequencies across. We’ll do this in the simplest way first, before going into how you might take a programming approach to the same problem. Either way, generating a range is simple - we just need to do it for a single allele, \(A_1\) and then we can very easily derive the frequency of \(A_2\). We will do this like so:

# generate a range for p
p <- seq(0, 1, 0.01)
# and also for q
q <- 1 - p

All we did here was use the seq function, which we first saw back in Chapter 1. Remember, we only have to do it once because \(q\) is the inverse of \(p\). Next we need to generate the expected frequencies. Because R is really great at handling vectors, the code for this is identical for if we were using a single value.

# generate the expected genotype frequencies
A1A1_e <- p^2
A1A2_e <- 2 * (p * q)
A2A2_e <- q^2

So now we basically have the expected frequencies at HWE of our different genotypes across a range of allele frequencies with just a little vector calculation. However, recalling from before that we often need to get our data into a data.frame in order to use ggplot2, we need to do a bit of work on this first. We’ll now get our data into a tibble ready for working on…

# arrange allele frequencies into a tibble/data.frame
geno_freq <- as.tibble(cbind(p, q, A1A1_e, A1A2_e, A2A2_e))

## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.

You also might remember from the previous session that we need to sometimes do some data manipulation to get the data into a form that we can easily plot. Luckily, this is quite straightforward in this case as we can use the gather function to do it quickly and easily.

# Use gather to reshape the data.frame for straightforward plotting
geno_freq <- gather(geno_freq, key = "genotype", value = "freq", -p, -q)

Now it is really straightforward to make a plot that shows the expected HW—frequencies across our allele frequencies.

# plot the expected genotype frequencies
a <- ggplot(geno_freq, aes(p, freq, colour = genotype)) + geom_line()
a <-a + ylab("Genotype frequency") + xlab("p frequency")
a + theme_light() + theme(legend.position = "bottom")

Some programming tricks for the HW model

The example above, where we derived the expected frequencies of the HW-model took advantage of R’s ability to easily handle vectors. In principle, this is the simplest way to achieve what we wanted and it was quick and clear to code. However, there are actually several different ways we could have approached this issue using some programming tricks. We will introduce these briefly here as a demonstration but these are topics we will return to throughout these tutorials. First up, we will delve into using a for loop.

for(x in 1:10){
  print(x)
}

for(x in 1:10){
  print(x + 10)
}

# set up p
p_vec <- seq(0, 1, 0.01)
# initiate null vectors for the expected frequencies
A1A1_e <-  NULL
A1A2_e <- NULL
A2A2_e <- NULL

# run the for loop
for(p in p_vec){
  # define q
  q <- 1 - p
  # calculate genotype frequencies
  A1A1_e <- append(A1A1_e, p^2)
  A1A2_e <- append(A1A2_e, 2 * (p * q))
  A2A2_e <- append(A2A2_e, q^2)
}

# combine them into a matrix (as an example here)
cbind(A1A1_e, A1A2_e, A2A2_e)

sapply(1:10, function(x) x)

sapply(1:10, function(x) x + 10)

# set up p
p_vec <- seq(0, 1, 0.01)

geno_freq <- sapply(p_vec, function(p){
  # define q
  q <- 1 - p
  # calculate probabilities
  A1A1_e <- p^2
  A1A2_e <- 2 * (p * q)
  A2A2_e <- q^2
  # write out
  c(A1A1_e, A1A2_e, A2A2_e)
})

# transpose the matrix to get it the right way up!
geno_freq <- t(geno_freq)

An easier way?

Since the chi-squared test is such a fundamental statistical tool, it probably won’t surprise you to find that there is actually already a chisq.test function for you to use. We didn’t explain it earlier because… well calculating things by hand is sometimes the best way to learn! The function is handy though since we can calculate the chi-squared value and also test the significance in one go. We can use it like so:

# perform a chi-squared test
mychi <- chisq.test(observed, p = expected_freq)

Here we provide the chisq.test function with our observed genotype counts. However, we also need to specify our expected frequencies, so there is still some calculation required. Nonetheless, the function returns the same chi-squared value we have previously calculated and shows us that we have a highly significant deviation from the expectation.

Note as well that the mychi object has a lot of information stored within it. We can take a closer look using the following:

str(mychi)

Using the dollar sign notation we learned about in Chapter 1, we can extract any of these at will. For example mychi$expected will produce the expected values whereas mychi$statistic will produce the chi-squared value.

A simple case of genetic drift over a single generation using the binomial distribution

To demonstrate a simple case of drift, we’ll use the example from the book. We have a population of \(N\) = 8 individuals and we will focus on a single locus \(A\) with two alleles, \(A_1\) and \(A_2\). We will assume both alleles are equally common in the population - i.e. \(p = q\) = 0.5. Our little population is in HWE such that there are two of each \(A_1A_1\) and \(A_2A_2\) homozygotes and four \(A_1A_2\) heterozygotes.

To simulate drift, we need to allow individuals to mate at random. In the book you will recall we used a coin-flipping experiment to sample our alleles to form genotypes, since both alleles have a frequency of 0.5. This is essentially a form of binomial sampling. We can easily recreate this in R using the rbinom function to sample the binomial distribution. We will recreate the coin-tossing experiment exactly here, using the following code:

# recreate coin flipping experiment using binomial sampling
rbinom(n = 16, size = 1, prob = 0.5)

What did we do here? Like with the pchisq function we used before, we made use of R’s suite of functions for statistical distributions (a topic we will cover in more detail later). The rbinom function just allows us to sample the binomial distribution. With the n = 16 flag, we sampled it 16 times. We set the size to 1 - meaning we sample either 0 or 1 (i.e. heads or tails) and finally we set the probability of sampling to the same for both outcomes with p = 0.5.

However, in the textbook we actually used pairs of coin flips to simulate genotypes. We could do this in R, but it is a bit clunky (it is much easier to just flip a coin!) so we could actually use binomial sampling in a different way to achieve this. We could set the size of the population as N and then just sample one value from the distribution. This value could then represent the number of \(A_1\) alleles in the next generation - so we could easily work out \(p´\) and \(q´\). We’ll also take this opportunity to set the value of \(p\) here.

# set population size
N <- 8
# set allele frequency
p <- 0.5 
# recreate coin flipping experiment using binomial sampling
nA1 <- rbinom(1, size = 2*N, prob = p)

Remember since we are sampling alleles for diploid individuals, we need to multiply N by 2. The value of nA1 will differ for everyone since this is a random sample of the binomial distribution. We can now easily work out the value of \(p´\) - i.e. the new frequency of the \(A_1\) in the next generation. Note that R doesn’t allow the prime notation so for the actual R code, we will use pa.

# calculate p´
pa <- nA1/(2*N)

OK so now we have a value for pa - again it will differ for everyone because we performed random sampling. You will see though that pa is likely different to p solely because of random sampling - i.e. we have demonstrated drift over a single generation. To properly understand drift over multiple generations, we will need to use some quite advanced R code - so for the next section, we will take a break from population genetics and instead focus on how to write a custom R function.

Writing your own R functions

Writing functions in R is an extremely useful thing to learn. It might seem a little daunting but it is actually quite straightforward once you have a handle on the basic format and structure necessary. Functions can be useful for cleaning up your code - instead of writing the same code block again and again, you can use your own function! You also might want to do it to then wrap it in a more complicated operation - like a for loop or a vectorised approach. We will actually build towards something like the latter, but for now we will just use a simple example. Let’s write a function that adds 10 to any numerical value we give it.

# write a custome function to add 10 to every numerical variable it is given
add_ten <- function(x){
  y <- x + 10
  return(y)
}
# test it out
add_ten(100)
add_ten(1000000)
add_ten(22)

Let’s breakdown what we did here in a bit more detail. Firstly, we need to give our function a name - in this case we called it add_ten. Crucially, this is also an object we assign. Once assigned, the function is callable in the R environment, just like any other function we load from a package or in the standard R distribution. We then use the statement function(x) to show it is a function. The x here is the name of the variable we pass to our function. The name here is completely arbitrary; we could just as easily use a human name or anything we choose. However, what we pass to the function will be referred to this throughout.

Inside the curly brackets - {} is the body of our function. Here we define a new object y that is simply the value of x plus 10. Finally, we use the return() function to tell R to produce the value of y when we use our add_ten function. We can really write a function for almost anything in R and it wouldn’t take us long to fill this tutorial with example after example. For now though, we will show you just one more example - a function that takes two arguments, instead of just one like our add_ten function did. This time, we will write a function to multiply two numerical values.

# write a custome function to multiply two numerical variables
multiply <- function(x, y){
  z <- x * y
  return(z)
}
# test it out
multiply(10, 10)
multiply(3, 5)
multiply(22, 44)

So this is exactly the same as before, but now we are passing two arguments - x and y. In the body of the function, we multiply them both and then return the answer. Hopefully this demonstration of how to write functions is a sufficient example. As with our previous programming topics, it is something we will return to throughout the course, building on it further and further. Again, we certainly don’t expect you to go from these examples to writing complex functions!

Using a custom function to simulate drift over multiple generations

Recalling our code for simulating genetic drift, let’s do the same thing again but this time wrap it in a simple function. To start with we do this only for a single generation.

# a custom function for simulating drift in a single generation
drift_sim <- function(N, p){
    pA <- rbinom(1, 2*N, p)
    p <- pA/(2*N)
  return(p)
}
# run drift sim on the example data above
drift_sim(N = 16, p = 0.5)

Now we have the drift_sim function, we can easily play around with the allele starting frequency(i.e. p) and also the population size. But we still haven’t managed to extend our simulations across multiple generations. Now we have the basic function, we can alter it to do just that:

# a custom function for simulating drift across multiple generations
drift_sim <- function(N, p, ngen){
  # initialise p
  p_init <- p
  # sample across all the generations
  pvec <- sapply(1:ngen, function(x){
    pA <- rbinom(1, 2*N, p)
    p <<- pA/(2*N)
  })
  # create a vector of p over time
  p <- c(p_init, pvec)
  # write out
  return(p)
}
# run drift sim on the example data above
drift_sim(N = 16, p = 0.5, ngen = 100)

OK so this will require a little explaining! Firstly you will notice we have added the ngen argument - which is the number of generations. Inside the function, we need to define the initial value of p - i.e. before the simulations begin - pinit. We then use an sapply call to use run our code over each generation. Each time we do this, we update the value of p - this explains the double arrow <<- we used. We have to do this updating because otherwise the value of p will also be reset to the initial value, which is not suitable when simulating drift over time. The final part of the function outputs our vector p and includes the original, initiating value first.

Now that we have a function to simulate data, across generations, we can easily do this for different population sizes. Let’s simulate 1000 generations of drift for three populations of 10, 100 and 1000 individuals. We’ll keep the initial allele frequency, p at 0.5 in all cases. We will then combine them all into a tibble alongside the number of generations

# simulate drift for three different population sizes
n10 <- drift_sim(N = 10, p = 0.5, ngen = 1000)
n100 <- drift_sim(N = 100, p = 0.5, ngen = 1000)
n1000 <- drift_sim(N = 1000, p = 0.5, ngen = 1000)
# get number of generations
g <- seq(0, 1000, 1)
# combine into a tibble
mydrift <- as.tibble(cbind(g, n10, n100, n1000))

Note that we could also have run drift_sim in a single line of code using sapply like so:

# simulate drift for three different population sizes
a <- sapply(c(10, 100, 1000), drift_sim, p = 0.5, ngen = 1000)
# rename matrix
colnames(a) <- c("n10", "n100", "n1000")
# get number of generations
g <- seq(0, 1000, 1)
# combine into a tibble
mydrift <- as.tibble(cbind(g, a))

Both ways work, but this is a nice demonstration of how you can use something like sapply to make your code more succinct. Anyway, what we ultimately want to do is visualise drift over time. So we need to get our data in a position to plot it using ggplot2 - once again, we turn to gather.

# use gather to get data ready for plotting
mydrift <- gather(mydrift, key = "pop_size", value = "p", -g)

All we did here was collapse the population sizes into a single factor for each generation. gather commands can take a while to get the hang of, the best thing is to continually practice them like this! Now we will plot our data using the following code:

# plot data
p <- ggplot(mydrift, aes(g, p, colour = pop_size))
p <- p + geom_line()
p <- p + xlab("No. generations") + ylab("Allele freq (p)")
p + theme_bw() + theme(legend.position = "bottom", legend.title = element_blank()) 

And there we have it! A plot of our simulated data showing allele frequency changes over time as a result of genetic drift. Although the plot you get will differ from the one you see here (because of the random sampling), what you should be able to see is that with a low population size alleles are either lost or go to fixation very rapidly. This is because the effect of drift is much stronger here, whereas in larger populations it may take a long time for drift to fix or remove an allele - indeed in some cases it will not happen at all.