lecture-01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 01: Introduction to population genetics\n",
    "\n",
    "### Random mating in dioecious populations\n",
    "\n",
    "Previously, we considered a population of hermaphroditic individuals that mate randomly. If the population contains two alleles, A with frequency $p$ and B with frequency $1-p$, the frequency of the genotypes AA, AB, and BB in the next generation are $p^2$, $2p(1-p)$, and $(1-p)^2$, respectively. This is [Hardy-Weinberg equilibrium](https://en.wikipedia.org/wiki/Hardy%E2%80%93Weinberg_principle).\n",
    "\n",
    "Does the same law hold if a species is [dioecious](https://en.wikipedia.org/wiki/Dioecy)? In other words, the population consists of two sexes, F and M, that reproduce with each other. To explore this question, we'll numerically solve problem 1.3 in Gillespie's _Population Genetics_ textbook. We'll assume that half the population is F and half is M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's write down the genotype frequencies for F and M in the first generation\n",
    "# Choose your own values here\n",
    "\n",
    "F_AA = 0.1  # frequency of the AA genotype in F\n",
    "F_AB = 0.3  # frequency of the AB genotype in F\n",
    "F_BB = 1 - F_AA - F_AB  # frequency of BB genotype in F\n",
    "\n",
    "M_AA = 0.4  # frequency of the AA genotype in M\n",
    "M_AB = 0.2  # frequency of the AB genotype in M\n",
    "M_BB = 1 - M_AA - M_AB  # frequency of BB genotype in M\n",
    "\n",
    "\n",
    "# Let's summarize our choices\n",
    "\n",
    "print('sex\\tAA\\tAB\\tBB')\n",
    "print('F\\t%.2f\\t%.2f\\t%.2f' % (F_AA, F_AB, F_BB))\n",
    "print('M\\t%.2f\\t%.2f\\t%.2f' % (M_AA, M_AB, M_BB))\n",
    "\n",
    "F_A = F_AA + F_AB/2  # frequency of the A allele in F\n",
    "F_B = F_BB + F_AB/2  # frequency of the B allele in F\n",
    "\n",
    "M_A = M_AA + M_AB/2  # frequency of the A allele in M\n",
    "M_B = M_BB + M_AB/2  # frequency of the B allele in M\n",
    "\n",
    "print('')\n",
    "print('sex\\tA\\tB')\n",
    "print('F\\t%.2f\\t%.2f' % (F_A, F_B))\n",
    "print('M\\t%.2f\\t%.2f' % (M_A, M_B))\n",
    "print('avg\\t%.2f\\t%.2f' % ((F_A + M_A)/2, (F_B + M_B)/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the Hardy-Weinberg law, we should find that the genotype frequencies are given by products of the allele frequencies. Does this happen in dioecious populations? How many generations does it take? Let's calculate below to find out. First we'll write out H-W expectations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Averaged across both sexes, the A and B allele frequencies are\n",
    "\n",
    "avg_A = (F_A + M_A)/2\n",
    "avg_B = (F_B + M_B)/2\n",
    "\n",
    "# Then, the expected genotype frequencies at H-W equilibrium would be\n",
    "\n",
    "HW_AA = avg_A**2\n",
    "HW_AB = 2 * avg_A * avg_B\n",
    "HW_BB = avg_B**2\n",
    "\n",
    "print('\\tAA\\tAB\\tBB')\n",
    "print('HW\\t%.2f\\t%.2f\\t%.2f' % (HW_AA, HW_AB, HW_BB))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's compute genotype frequencies over the next (few) generation(s). We'll write a function to make this easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start by loading the genotype frequencies into lists\n",
    "\n",
    "F_gts_start = [ F_AA,  F_AB,  F_BB]  # F genotype frequencies\n",
    "M_gts_start = [ M_AA,  M_AB,  M_BB]  # M genotype frequencies\n",
    "HW_gts      = [HW_AA, HW_AB, HW_BB]  # H-W equilibrium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, let's write a function to take those frequencies and get new genotype frequencies after random mating\n",
    "\n",
    "def rm(F_gts, M_gts): \n",
    "    # Gillespie suggests to enumerate all possible pairings and then to compute how they contribute to\n",
    "    # genotype frequencies in the next generation\n",
    "    \n",
    "    F_gts_new = [0, 0, 0]  # F genotype frequencies in the next generation\n",
    "    M_gts_new = [0, 0, 0]  # M genotype frequencies in the next generation\n",
    "    \n",
    "    # F AA, M AA\n",
    "    # only produces AA offspring\n",
    "    f = F_gts[0] * M_gts[0]\n",
    "    F_gts_new[0] += f\n",
    "    M_gts_new[0] += f\n",
    "    \n",
    "    # F AA, M AB -- same as F AB, M AA\n",
    "    # 50% chance AA, 50% AB\n",
    "    f = (F_gts[0] * M_gts[1]) + (F_gts[1] * M_gts[0])\n",
    "    F_gts_new[0] += f/2\n",
    "    M_gts_new[0] += f/2\n",
    "    F_gts_new[1] += f/2\n",
    "    M_gts_new[1] += f/2\n",
    "    \n",
    "    # F AA, M BB -- same as F BB, M AA\n",
    "    # only produces AB offspring\n",
    "    f = (F_gts[0] * M_gts[2]) + (F_gts[2] * M_gts[0])\n",
    "    F_gts_new[1] += f\n",
    "    M_gts_new[1] += f\n",
    "    \n",
    "    # F AB, M AB\n",
    "    # 25% AA, 50% AB, 25% BB\n",
    "    f = F_gts[1] * M_gts[1]\n",
    "    F_gts_new[0] += f/4\n",
    "    M_gts_new[0] += f/4\n",
    "    F_gts_new[1] += f/2\n",
    "    M_gts_new[1] += f/2\n",
    "    F_gts_new[2] += f/4\n",
    "    M_gts_new[2] += f/4\n",
    "    \n",
    "    # F AB, M BB -- same as F BB, M AB\n",
    "    # 50% AB, 50% BB\n",
    "    f = (F_gts[1] * M_gts[2]) + (F_gts[2] * M_gts[1])\n",
    "    F_gts_new[1] += f/2\n",
    "    M_gts_new[1] += f/2\n",
    "    F_gts_new[2] += f/2\n",
    "    M_gts_new[2] += f/2\n",
    "    \n",
    "    # F BB, M BB\n",
    "    # only produces BB offspring\n",
    "    f = F_gts[2] * M_gts[2]\n",
    "    F_gts_new[2] += f\n",
    "    M_gts_new[2] += f\n",
    "    \n",
    "    # Now that we've computed all the contributions, return the vectors\n",
    "    return F_gts_new, M_gts_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's call the function to get new genotype frequencies, then print them out\n",
    "\n",
    "F_gts_2, M_gts_2 = rm(F_gts_start, M_gts_start)\n",
    "\n",
    "print('sex\\tAA\\tAB\\tBB')\n",
    "print('F\\t%.2f\\t%.2f\\t%.2f' % (F_gts_2[0], F_gts_2[1], F_gts_2[2]))\n",
    "print('M\\t%.2f\\t%.2f\\t%.2f' % (F_gts_2[0], F_gts_2[1], F_gts_2[2]))\n",
    "print('H-W\\t%.2f\\t%.2f\\t%.2f' % (HW_gts[0], HW_gts[1], HW_gts[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's check the results after one more generation\n",
    "\n",
    "F_gts_3, M_gts_3 = rm(F_gts_2, M_gts_2)\n",
    "\n",
    "print('sex\\tAA\\tAB\\tBB')\n",
    "print('F\\t%.2f\\t%.2f\\t%.2f' % (F_gts_3[0], F_gts_3[1], F_gts_3[2]))\n",
    "print('M\\t%.2f\\t%.2f\\t%.2f' % (F_gts_3[0], F_gts_3[1], F_gts_3[2]))\n",
    "print('H-W\\t%.2f\\t%.2f\\t%.2f' % (HW_gts[0], HW_gts[1], HW_gts[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hardy-Weinberg law also suggests that overall allele frequencies should stay constant under random mating. Is that true? Let's check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = F_gts_3[0] + F_gts_3[1]/2  # frequency of the A allele (now same in F and M)\n",
    "B = F_gts_3[2] + F_gts_3[1]/2  # frequency of the B allele (now same in F and M)\n",
    "\n",
    "print('')\n",
    "print('\\tA\\tB')\n",
    "print('avg\\t%.2f\\t%.2f' % (A, B))\n",
    "print('old\\t%.2f\\t%.2f' % (avg_A, avg_B))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Genetic drift in finite populations\n",
    "\n",
    "We saw that in finite populations, chance in reproduction leads to random changes in allele frequencies, termed [genetic drift](https://en.wikipedia.org/wiki/Genetic_drift). A simple model for this is the Wright-Fisher model, in which the alleles in the next generation of a population are chosen at random with replacement from the total pool of $2N$ alleles that exist in the current generation. \n",
    "\n",
    "Let's call two types of alleles A and B, with frequencies $p$ and $1-p$. Under the Wright-Fisher model, the probability of getting $n$ alleles of type A in the next generation is binomial,\n",
    "\n",
    "$$ P(n) = \\binom{2N}{n} p^n \\left(1-p\\right)^{2N-n}\\,. $$\n",
    "\n",
    "Below, let's simulate evolution of the A allele frequency over time. To do this easily, we can use a `NumPy` library that generates random numbers from a binomial distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np          # here we import numpy\n",
    "import numpy.random as rng  # and here we import the random number generation (sub-)library\n",
    "\n",
    "\n",
    "# As a simple test, we can see how many heads we would get if we flip a fair coin 10 times\n",
    "\n",
    "p = 0.5\n",
    "N = 10\n",
    "heads = rng.binomial(N, p)\n",
    "\n",
    "print('We got %d heads out of %d coin flips' % (heads, N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, let's use this to simulate genetic drift\n",
    "# First, let's set the parameters\n",
    "\n",
    "p   = 0.5  # starting frequency of the A allele\n",
    "N   = 100  # population size\n",
    "n_g = 100  # number of generations of evolution\n",
    "p_t = np.zeros(n_g)  # a blank vector of allele frequencies over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next, let's fill in the vector of allele frequencies over time by simulating the WF model\n",
    "\n",
    "p_t[0] = p  # set the starting frequency to p\n",
    "\n",
    "for i in range(n_g-1):\n",
    "    n_next   = rng.binomial(2*N, p_t[i])  # get As in next generation\n",
    "    p_t[i+1] = n_next/(2*N)  # divide by number of alleles to get frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, let's make a plot\n",
    "\n",
    "import seaborn as sns            # import seaborn\n",
    "import matplotlib.pyplot as plt  # and matplotlib\n",
    "\n",
    "sns.lineplot(x=np.arange(n_g), y=p_t)\n",
    "plt.xlabel('Time (generations)')\n",
    "plt.ylabel('Allele frequency')\n",
    "plt.ylim(0, 1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}