Introduction

Since you are now familiar with R, we can proceed to our first lab in evolutionary processes. Today, we will focus on the fundamental entity of evolution: genes. We will first review how the laws of probability can be used to build very simple (but powerful!) models of genotypic evolutionary change. Then we will calculate some properties of and test hypotheses about such change. We will finish the lab by introducing methods to measure the genetic diversity generated by such processes.

It is important to note that we will change conventions quickly in this tutorial, meaning that the same letter might have different meanings in different equations, so be mindful of which section you are in.

Simple math notation

The first notation we will introduce is the summation, denoted by the uppercase Greek letter sigma (∑). Usually, the subscript of this letter tells you from which element to start the operation, and the upperscript tells you at which element the summation ends. For instance the notation:

indicates a summation of elements starting at the first (i = 1) and going until the n-th element (the upperscript). It is left implicit in this notation that the summation is made element by element, but sometimes this information might appear explicitly. If you are familiar with R’s for loops, it may be useful to think of the summation as representing a for loop that iterates from the sub- to the upperscript, each time applying the operation and adding its result to the running sum.

A familiar example is the formula to calculate the arithmetic mean:

As another example, take the following equation, which describes the expected number of heterozygotes in a diploid population with k alleles in Hardy-Weinberg equilibrium:

where k is the number of alleles, and p_i and p_j are the frequencies of the ith and jth alleles. Note that on the right side of the equality, we sum the product of all allele frequencies. The inequality i ≠ j tells you to not apply the summation when i = j. This kind of explicit information may appear in a summation.

So, in the case of k = 3, we would have the following equation:

$H_{k=3} = \sum_{i=1 ; i\neq{j}}^{3}p_ip_j = (p_1 \times p_2) + (p_1 \times p_3) + (p_2 \times p_1) + (p_2 \times p_3) + (p_3 \times p_1) + (p_3 \times p_2)$

Just for the sake of completion, one could also say that:

$H = 1-G = 1- (\sum_{i=1}^{k}p_i^2)$

There is an equivalent of the summation for multiplication, and it is symbolized by the uppercase Greek letter pi (∏). So, the probability P of n independent events v, each having the same probability P_v is:

With v representing each event. Note that the upper- and subscript are read in the same way as in the summation.

Probability of independent events

The joint probability of independent events is the multiplication of each respective probability (a quick video recap on that here).

So, the probability of seeing “snake eyes” (i.e. two “ones”) when you throw two dice is:

The notation in the middle equality above might seem overly complicated when compared to the left side of the equality, but the compactness of this notation becomes clear when we want to multiply many elements. R has functions that allow you to “throw dice” and study the probability of events. These are the functions that generate random numbers from probability functions.

Random number generators and density/mass probability functions

Using functions that generate random numbers, we can simulate much more than just coins and dice. R has functions that generate random draws from famous probability distributions.

To explore this concept, we will focus on one of these famous distributions: the binomial distribution. The binomial distribution applies for samples of events that have a binomial outcome, such as coin flips which result in either a “heads” or “tails” outcome. In the jargon of the field, we call one of these outcomes a “success” and the other one a “failure”; this makes more sense in the coin flip analogy if we imagine that we are using the coin to determine the result of some outcome, and one of the sides has the more favorable outcome (e.g. eat a piece of chocolate only when the coin lands on heads, no chocolate on tails).

The binomial distribution has two major properties:

1. it counts the number of “successes” resulting from a sample of binomial events
2. the probability of success is the same between events

The binomial distribution is therefore a distribution of integers that each represent how many successes there were within a sample of binomial events. The distribution is made by taking many repeated such samples. Each time we generate a random number from a binomial distribution, we are generating an integer x between 0 (no successes occurred) and the total number of coin flips n (only successes occurred). In the coin tossing example, x is equivalent to counting heads out of n coin tosses.

Like all probability functions, the binomial distribution has certain parameters. These are: p, the probability of success for a single binary event, and n, the total number of binary events. In the coin flip example, p would be 0.5, representing a 50% chance of getting heads. n would be however many times we flip the coin.

Let’s simulate a coin flip. Recall that the binomial distribution outputs the number of successes from a sample of events, so this is what each generated random number represents.

# Define how many samples we want:
rr <- 10000

We use the rbinom function to randomly generate numbers from the binomial distribution. rbinom is in the r-family of functions, the collection of functions that randomly generate numbers from probability distributions.

binom_nbrs <- rbinom(n = rr, size = 10, prob = 0.5)  # taking 10000 samples of 10

A potential point of confusion is that, in rbinom, the argument size represents the parameter n discussed above, while the argument n actually represents the total number of times we repeated sampling. The argument prob represents the parameter p discussed above, the probability of success – much more straightforward. The line of code above represents flipping a coin (with a 50% chance of heads) 10 times and counting the number of heads flipped a total of rr times.

If we check what the object binom_nbrs contains, we will see…

head(binom_nbrs, 10)

##  [1] 2 6 6 5 4 6 2 4 4 8

… a bunch of integers between 0 and 10. Each of these integers represents the number of heads flipped in a sample of 10 coin tosses.

As the binomial is a discrete probability function, it is meaningful to tabulate the frequencies of each number we got from the random draw we just made:

table(binom_nbrs)

## binom_nbrs
##    0    1    2    3    4    5    6    7    8    9   10 
##    7   84  430 1142 2055 2534 1972 1221  444   96   15

We can also plot what we did very easily.

plot(
  table(binom_nbrs),
  ylab="Frequency", xlab="Number of successes")

Above are the random numbers we generated – i.e., the stochastic realizations of the process that follows the distribution we chose. Each number represents the summation of the number of successes per sample. But what if instead of having a stochastic realization of the binomial process, we wanted to know what is the fixed probability associated with each specific number of successes?

For this, we would use the d-family of functions! These calculate the probability density associated with an event. An event in this case is a number of successes. For instance, to directly calculate the probability of observing, say, 3 heads out of a sample of 10 coin tosses, we would write:

P_x3_s10_p0.5 <- dbinom(x = 3, size = 10, prob = 0.5) 
# x is the integer value we want to know the probability density of

P_x3_s10_p0.5

## [1] 0.1171875

This function basically just applies the probability mass function of the binomial distribution, which is

with: $\binom{n}{x} = \frac{n!}{x!(n-x)!}$

Knowing that dinom is simply using the function written above, you could manually calculate the expected probability density for x = 3, and check to make sure it is the same as what dbinom gives you. To do that, use equation (6) and substitute n with 10 and x with 3 – you should get exactly the same number as the output of dbinom, which is 0.1171875.

There is a difference here between the stochastic realization of a process, acquired with rbinom, and the theoretically expected outcome of a process, acquired with dbinom. If we draw many random numbers, we should be able to approximate the expected probability by calculating the proportion of each outcome.

# comparing the stochastic realization of x = 3 with the expectation for x = 3
table(binom_nbrs)[4] / rr #why the division?

##      3 
## 0.1142

P_x3_s10_p0.5

## [1] 0.1171875

What do you think? Is it a good approximation? How could we improve it?

We can also plot the random draws we made:

plot(
  table(binom_nbrs)/rr, 
  ylab="Probability", xlab="Number of heads")

# Finally, we can mark the probability we are interested in
# by using a red dot in our plot:
points(x = 3, y = P_x3_s10_p0.5, col = "red", cex = 2)

Malthusian growth

The idea of geometric/exponential/Malthusian growth is of major importance for evolutionary biology. You have probably learnt of this concept before, but we will reiterate it here. Although often called “Malthusian growth”, the basic analytic approach was first developed by Euler and published in 1748. He modeled population growth with equation (7):

where x is the growth rate of the population or the fractional increase per year, N is the population size (in number of individuals), and t is time.

popsize <- 10

time <- 0

R <- 1.05

tmax <- 100 # this is the number of generations

all_generations <- seq(from = 1, to = tmax, by = 1)

# For each generation, beginning in the generation = 1 and:
for(generation in all_generations){ 
  
  N_t <- popsize[length(popsize)] # by indexing 
  # the vector in this way (i.e. using "[length()]"),
  # we guarantee we are always taking the last value of 
  # this vector.
  
  N_t_plus_one <- N_t * R # this is the application of Euler's formula
  
  #Now, let's store the population size at that generation:
  popsize <- c(popsize, N_t_plus_one)
  
  #let's not forget to record the time that we are on:
  time <- c(time, generation) 
}

#
plot(NA, xlim=c(0,tmax), ylim=c(0,1500),
     xlab="time", ylab="Population size",
     main="Malthusian growth")
lines(x = time, y = popsize, col = "blue", lwd = 3)

Mendelian genetics and Hardy-Weinberg Equilibrium (HWE) at a single locus

We will now use simulations to test if the observed set of genotype counts is consistent with what we would see if we sampled genotypes at random. We will first create a null pattern using R functions, then we will contrast this pattern with the observed genotype frequencies and compare them to expectations under Hardy-Weinberg equilibrium (hereafter, HWE, Hardy, 1908; Weinberg, 1908). For a historical perspective on HWE, we recommend reading Mayo (2008).

For the purposes of this lab, we will look to a bear population in Canada.

Kermode bear from British Columbia

Hedrick & Ritland (2012) studied a population of Ursus americanus where a “single nucleotide change from G to A, resulting in the replacement of Tyr to Cys at codon 298 in the melanocortin 1 receptor gene (mc1r)” created a whole new phenotype related to color: the “Spirit bear” variety (see figure 1, from Hedrick & Ritland (2012)).

A white mother Spirit bear and a black cub offspring; the father must have been black and the cub is a heterozygote for the coat color polymorphism (photo mod. from Hedrick and Ritland (2012)

In a previous study (Ritland et al, 2001), the following frequencies of genotypes were observed: 42 individuals with the dominant homozygote phenotype, 24 heterozygotes, and 21 recessive genotypes. Hedrick & Ritland (2012) asks: are the Spirit bear frequencies a result of neutral processes? This is equivalent of asking is this genotype frequency a simple reflex of Hardy-Weinberg Equilibrium (HWE)?

sim_pops <- OneGenHWSim(n.ind = 100, n.sim = 5, p = 0.467)

#now, checking the result of our simulations:
sim_pops

##       A1A1 A1A2 A2A2
## sim_1   19   56   25
## sim_2   21   52   27
## sim_3   19   56   25
## sim_4   22   49   29
## sim_5   19   55   26

# To remember how this works, let's imagine you want to 
# arbitrarily calculate (f(A1) + 3) / 4.5 for all your simulations. 
#You should run:

freqs_A1 <- sim_pops$A1A1 / (sim_pops$A1A1 + sim_pops$A1A2 + sim_pops$A2A2)

result <- (freqs_A1 + 3) / 4.5
result

## [1] 0.7088889 0.7133333 0.7088889 0.7155556 0.7088889

# the above has no biological "meaning". 
# it is just to remind you how vectorized calculation works

Heterozygosity

As defined in the lecture, “heterozygosity” might mean two things:

1. the probability that 2 alleles sampled from a population are different.
2. The actual frequency of heterozygotes in the population.

HWE, deleterious alleles, and mutation

Mendelian genetics at multiple loci

Describing genetic variation in DNA segments

As you saw in the lecture, some of the most remarkable first studies of genetic variation are Hubby & Lewontin (1966) and Lewontin & Hubby (1966). To deepen your perspective on these studies and their historical and current applications, we recommend reading Charlesworth et al (2016).

Genetic diversity has had an important role in the study of biological diversity. For the purposes of this class, we will take an extremely simplifying approach and calculate different, but simple statistics of “genetic diversity”.

Consider the following aligned DNA sequences from 3 genes (genes separated by space):

Individual 1    ACCGTA   AAAAAT   CTTATA     
Individual 2    AGCGGA   CATAAT   CTTATA     
Individual 3    ACCGTA   AAAAAT   CTACTA     
Individual 4    ACCGGA   AAAAAT   CTACTA