genetic_drift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EEB214: Understanding genetic drift through computer-based simulations in [`R`](https://www.r-project.org/)\n",
    "This tutorial (i.e. Jupyter Notebook) will guide you through some hands-on example on how *Genetic Drift* operates within populations. Read and follow the instructions below. Enter the information requested and run each block of code separately. Some blocks of code will give you an output, such as text or a plot, others will not.\n",
    "\n",
    "## Background\n",
    "Imagine you have a population of a species of bacteria. Individuals (i.e. cells) in the population will have either the `A` or the `a` allele at a frequency set initially by you. Then you will let the population reproduce and every generation you will select randomly which individuals pass on their genes. Now, all population genetic models have a set of (normally simplifying) assumptions. In our case, we will assume that the population size (total number of individuals at any given generation) will remain constant and that `A/a` is a **neutral polymorphism**. That is, it is not affected by selection. Also, remember that bacteria are [haploid](https://simple.wikipedia.org/wiki/Ploidy) (only one or a single set of chromosomes) as opposed to humans being [diploid](https://simple.wikipedia.org/wiki/Ploidy) (two sets of chromosomes), meaning that every individual will only have one  of the alleles.\n",
    "\n",
    "The way the simulation works is that every generation we will randomly sample \"cells\" with either `A` or `a`. We'll keep on doing that until one allele gets \"fixed\" or goes \"extinct\" in the population. For tractability, we will only focus on the fate of allele `A`.\n",
    "\n",
    "## Goals\n",
    "* Run through the exersizes\n",
    "* Try to understand what's happening in each step\n",
    "* Use your results (text and plots) to answer the quiz on [Quercus](https://q.utoronto.ca/). How does changing the initial parameters of the simulation help you understand how *Genetic Drift* operates.\n",
    "\n",
    "### Let's get started\n",
    "First load the libraries and functions needed for the exersize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load libraries and functions\n",
    "source('genetic_drift.R')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's set the initial parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_size = 10\n",
    "initial_frequency = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are [variables](https://en.wikipedia.org/wiki/Variable_(computer_science)) or placeholders. They are basically names that occupy a space in the memory of the computer and store a value.\n",
    "\n",
    "We can recall these at anytime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the print function prints to screen\n",
    "print(population_size)\n",
    "print(initial_frequency)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run our first simulation with these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation = genetic_drift(population_size, initial_frequency)\n",
    "# first 6 rows of the table\n",
    "head(simulation)\n",
    "# last 6 rows of the table\n",
    "tail(simulation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `x` and `y` are randomly generated coordinates (**don't worry about them! They are for plotting**). The last set of rows under `generation` are also telling us at which generation the simulation stopped, as well as the frequency of allele `A`.\n",
    "\n",
    "### Let's look at how our population of bacteria might look\n",
    "We will plot the population at generation `1`. Imagine the dark circle is our [Petri dish](https://en.wikipedia.org/wiki/Petri_dish) and the colored dots are our colonies with either the `A` or the `a` allele. Since our `initial_frequency` is set to `0.5`, this means that bacteria with `A` or with `a` will be at equal numbers at generation 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = plot_population(simulation, 1, population_size)\n",
    "p1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the last generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_generation = simulation[nrow(simulation),\"generation\"]\n",
    "# print(last_generation)\n",
    "p_last = plot_population(simulation, last_generation, population_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_last"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What happend here? Where did the other allele go?**\n",
    "\n",
    "### Plot the frequecy of allele `A` through time\n",
    "Now let's look at the frequency of allele `A`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_freq = plot_frequency(simulation, population_size)\n",
    "p_freq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the frequency of allele `A` at the last generation? Did it fix or did go extinct? How many generations did it take?\n",
    "\n",
    "### Play with the parameters: population size\n",
    "Now let's see what happens if we change the parameters a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_size = 100\n",
    "simulation = genetic_drift(population_size, initial_frequency)\n",
    "plot_frequency(simulation, population_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What happened now? What changed?**\n",
    "\n",
    "Let's try repeating these two population sizes multiple times, then average the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for storing plots\n",
    "plots = vector(\"list\",12)\n",
    "# for the last generation in each simulation\n",
    "last_gens = vector(\"numeric\",12)\n",
    "last_freq = vector(\"numeric\",12)\n",
    "# vector of population sizes\n",
    "popsize = 10\n",
    "# loop through and save plot\n",
    "for (i in 1:12){\n",
    "    simulation = genetic_drift(popsize, initial_frequency)\n",
    "    last_gens[i] = simulation[nrow(simulation),\"generation\"]\n",
    "    last_freq[i] = simulation[nrow(simulation),\"frequency\"]\n",
    "    plots[[i]] = plot_frequency(simulation, popsize)\n",
    "}\n",
    "# this might take a while... Be patient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first 4\n",
    "do.call(plot_grid, c(plots[1:4], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[5:8], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[9:12], ncol=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What do the** <span style=\"color:red\">red</span> **and** <span style=\"color:green\">green</span> **colors stand for? Can you note any patterns?**\n",
    "\n",
    "For a population size of 10 let's look at the averge time it takes an allele to get fixed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(median(last_gens))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the number of times the allele got fixed or got lost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = table(last_freq)\n",
    "names(tab) = c(\"lost\",\"fixed\")\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on these numbers we can calculate the probability an allele will get fixed as a function of it's initial frequency and a assuming a random sampling process (*Genetic Drift*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab/sum(tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try increasing the population size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for storing plots\n",
    "plots = vector(\"list\",12)\n",
    "# for the last generation in each simulation\n",
    "last_gens = vector(\"numeric\",12)\n",
    "last_freq = vector(\"numeric\",12)\n",
    "# vector of population sizes\n",
    "popsize = 100\n",
    "# loop through and save plot\n",
    "for (i in 1:12){\n",
    "    simulation = genetic_drift(popsize, initial_frequency)\n",
    "    last_gens[i] = simulation[nrow(simulation),\"generation\"]\n",
    "    last_freq[i] = simulation[nrow(simulation),\"frequency\"]\n",
    "    plots[[i]] = plot_frequency(simulation, popsize)\n",
    "}\n",
    "# this might take a while... Be patient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first 4\n",
    "do.call(plot_grid, c(plots[1:4], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[5:8], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[9:12], ncol=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the average number of generations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(median(last_gens))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compare this result with the previous one (i.e. population size = 10). What's different?**\n",
    "\n",
    "Let's now look at how many times it fixed/went extinct and the probability of fixation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = table(last_freq)\n",
    "names(tab) = c(\"lost\",\"fixed\")\n",
    "tab\n",
    "# probability\n",
    "tab/sum(tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compare the probability of fixation with the initial frequency. Is there a relationship?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Play with the parameters: initial frequency\n",
    "Now let's see what happens if we changes the initial frequency a bit.\n",
    "\n",
    "Let's run the same simulation, but this time let's start with a frequency of `0.1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for storing plots\n",
    "plots = vector(\"list\",12)\n",
    "# for the last generation in each simulation\n",
    "last_gens = vector(\"numeric\",12)\n",
    "last_freq = vector(\"numeric\",12)\n",
    "# vector of population sizes\n",
    "initfreq = 0.1\n",
    "popsize = 100\n",
    "# loop through and save plot\n",
    "for (i in 1:12){\n",
    "    simulation = genetic_drift(popsize, initfreq)\n",
    "    last_gens[i] = simulation[nrow(simulation),\"generation\"]\n",
    "    last_freq[i] = simulation[nrow(simulation),\"frequency\"]\n",
    "    plots[[i]] = plot_frequency(simulation, popsize)\n",
    "}\n",
    "# this might take a while... Be patient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first 4\n",
    "do.call(plot_grid, c(plots[1:4], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[5:8], ncol=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next 4\n",
    "do.call(plot_grid, c(plots[9:12], ncol=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Is there a pattern?**\n",
    "\n",
    "Let's look at the average time to fixation or loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(median(last_gens))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**If it's the same population size (population_size = 100), why is the time lower now?**\n",
    "\n",
    "Let's now look at the times it fixed/went extinct and the probability of fixation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = table(last_freq)\n",
    "names(tab) = c(\"lost\",\"fixed\")\n",
    "tab\n",
    "# probability\n",
    "tab/sum(tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's round those numbers to their nearest decimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "round(tab/sum(tab),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__What is the relationship between the probability of fixation of an allele and it's initial frequency?__"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}