3.Recombining_genomes_and_treesequences_STUDENT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical 3. Simulating recombining genomes - a simulation journey to the happy land of treesequences \n",
    "\n",
    "At the end of this practical you should:\n",
    "\n",
    "- Know what ARGs and treesequences are, how they are constructed and how they differ\n",
    "- Appreciate that not all recombination events are detectable\n",
    "- Understand the difference between map and physical length of a sequence\n",
    "- Know that the span of trees along the genome is a random variable and that nodes are shared between many trees\n",
    "- Understand that the duality between branch lengths and popgen measures extends to correlated trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import msprime\n",
    "import tskit\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from IPython.display import Image, SVG, set_matplotlib_formats\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "#jupyter magic\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Three theoretical flies! A simple coalescent simulation with recombination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen:</b>\n",
    "So far we have thought about the coalescent as a model of ancestry that applies to non-recombining \"loci\". While <i>a locus</i> is a useful concept in population genetics precisely because it is rather vague, it is not particularly helpful when thinking about genome sequence data. Some types of sequence (e.g. organelle genomes such as mitochdria) are non-recombining. However, most DNA sequences undergo recombination. In this case, the most basic/indivisible unit of inheritance, i.e. the only meaningful locus, is a single nucleotide.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have heard (lecture 3) how we can include recombination into the standard coalescent. In this practical, we will use simulations to explore basic properties of the coalescent with recombination.\n",
    "\n",
    "To simulate the ancestry of a sample of recombining sequences with `msprime`, we have to specify a `sequence_length` and a `recombination_rate`. It is alsways useful to make simulations both as simple as possible but also concrete, i.e. to choose parameters that apply at least roughly to some real population in nature. We will focus on two cases in this practical: *Drosophila melanogaster* and *Homo sapiens*. Population genomic variation has been studied excessively in both models and, together, they span a broad range of population genetic parameter space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 1:</b> The length of the genetic map of <i>D. melanogaster</i> is <a href=\"https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002905\">282 cM</a>. Given a genome size of $\\approx$ 180 Mb, what is the per base, per generation rate of recombination? How does this compare to the mutation rate estimated by <a href=\"https://genome.cshlp.org/content/19/7/1195\">Keightley et al (2009)</a>?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 2:</b> Set the per base recombination rate to $1/2$ of your estimate (given that recombination in <i>Drosophila</i> is restricted to females) to simulate the ancestry of a neutrally evolving sequence of length 50kb sampled in 3 individuals. How many marginal trees contribute to the ancestry of this sequence? Assume $N_e = 10^6$ which corresponds roughly to estimates for <i>D. melanogaster</i> (see practical 2).\n",
    "\n",
    "<b>Hint:</b> You can use <a href=\"https://tskit.dev/tskit/docs/stable/python-api.html?highlight=num_trees#tskit.TreeSequence.num_trees\"><code>num_trees</code></a> to count the trees in a ts (use random seed = 87).\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 3:</b> Use <code>draw_svg()</code> to look at the tree sequence. Going from left to right, can you work out which genealogical nodes each marginal tree does not share with its neighbour? On which branches could the recombination events be located that define the first 5 trees?\n",
    "\n",
    "<b>Hint:</b> to keep the figure manageable, use the <code>x_lim = [0,100]</code> argument to plot the first 100 bases. Note down the number of the node above each reccombination event. \n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen:</b>\n",
    "    \n",
    "Since each successive tree in the tree sequence is separated from its neighbours to the left and right by a single recombination breakpoint, we might expect that the average number of trees in a tree sequence is equal to the expected number of recombination breakpoints $E[R]+1$.\n",
    "\n",
    "The expected number of recombination events can be derived using the total tree length under the standard coalescent using the same logic we've used to derive $E[S]$, the expected number of mutations (see mini lectures 2 and 3)!\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 4:</b> Work out the population scaled recombination rate $\\rho_L = 4 N_e r \\times L$ for the region you have simulated. Analogous to $\\theta$, this can be interpreted as the expected number of recombination events in a sample of $n=2$.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 5: </b> What is the expected number of recombination events for the <i>D. melanogaster</i> like simulation of $n=6$ you have simulated above (see lecture 3)? \n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Use the formula for E[R] given in lecture 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 6: </b> How does this expectation compare to the number of trees in a 50 kb stretch of sequence? Is it likely that your <code>ts</code> has fewer trees than expected just by chance? Use the code below to simulate 100 replicate treesequences for the same  <i>D. melanogaster </i> parameters. Check the minimum and maximum number of trees and plot the distribution of the number of trees across sim replicates using <a href=\"https://seaborn.pydata.org/tutorial/distributions.html\"><code>displot</code></a>.\n",
    "\n",
    "<b>Hint:</b> The functions <code>min()</code> and <code>max()</code> are handy. \n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is using the numpy generator recipe: we are first defining an empty array: \n",
    "num_replicates = 100\n",
    "tree_count = np.zeros(num_replicates)#\n",
    "#Then setting up an ancestry simulation:\n",
    "repl_ts = msprime.sim_ancestry(3, population_size = 1_000_000,\n",
    "                              recombination_rate = r_dmel_sex_avg,\n",
    "                              sequence_length = 50_000, num_replicates = num_replicates, random_seed = 87)\n",
    "#For each replicate we record the number of trees:\n",
    "for repl_index, ts in enumerate(repl_ts):\n",
    "    tree_count[repl_index] = ts.num_trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "#this syntax is useful for labelling plot axes:\n",
    "#ax = sns.displot(x = data, stat=\"density\");\n",
    "#ax.set(xlabel='number of trees',\n",
    "#       ylabel='freq')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 7:</b> What is the distribution of sequence lengths trees span in your 50kb stretch of neutrally evolving <i>D. melanogaster </i> sequence? What is the mean tree span? Record the span of each tree as a <code>np.array</code> and plot the distribution of spans using <a href=\"https://seaborn.pydata.org/tutorial/distributions.html\"><code>displot</code></a>\n",
    "\n",
    "**Hint:** You need to set up a <code>for</code> loop to iterate over all trees in the ts. To do this, use the <a href=\"https://tskit.dev/tskit/docs/stable/python-api.html#tskit.TreeSequence.trees\"><code>trees()</code></a> method, which returns an iterator over the trees in a tree sequence.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "#Another case for the numpy generator plus for loop recipe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 8:</b> Record the $t_{mrca}$ of each tree (see practical 1!). What is the mean $t_{mrca}$ across trees? What are the units? How well does this agree with the analytic expectation $E[t_{mrca}]$ for this sample?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 9:</b> What is the relationship between tree span and tree depth, i.e. the $t_{mrca}$?\n",
    "\n",
    "**Hint:** A useful seaborn plot for this purpose is [regplot](https://seaborn.pydata.org/tutorial/regression.html#functions-to-draw-linear-regression-models)\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "#this plot sytax is useful:\n",
    "#fig, ax = plt.subplots(figsize=(10,10))\n",
    "#sns.regplot(x=span_list, y=tmrca_list, ax =ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 10:</b> Can you use the spans of trees to find the correct mean $t_{mrca}$ for your treesequence?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.2 From treesequences to haplotypes\n",
    "While a tree sequence records the exact correlation in ancestry along a genome, when analysing real data we generally do NOT know the underlying tree sequence. Measures of LD and tests for recombination events are generally based on pairs of variants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 11:</b> Sprinkle mutations on your <i>D. melanogaster</i> treesequence simulation according to the rate estimate of Keightley et al (2009). How many mutations are there on average per tree in this ts?\n",
    "\n",
    "**Hint:** set <code>discrete_genome = False</code> and <code>random_seed = 54</code>.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 12:</b> Check out the region 300-500 bases. What mutations can you see? Would mutations 18 and 19 be compatible with the same tree? What about 13 and 19? Does either pair of loci fail the 4-gamete test (see mini lecture 3)?\n",
    "\n",
    "<b>Hint:</b> use <code>.draw_svg(x_lim = [300,500])</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "#Since we have used the ancestry simulated in the previous section here, the ts must remains the same but now includes some mutations..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Once you have assigned the <i>D. melanogaster</i> ts (with mutations!) to a new variable, say <code>dmel_ts_muts</code>, you can store the haplotype alignment as an <code>np.array</code>:\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmel_array = np.zeros((dmel_ts_muts.num_sites, dmel_ts_muts.num_samples))\n",
    "for var_index, var in enumerate(dmel_ts_muts.variants()):\n",
    "    dmel_array[var_index] = var.genotypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "e.g. mutation 13 is shared by samples 3 and 5; mutation 19 is shared by all samples apart from 5: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This pulls out the list of haplotypes for variant 13. Samples 3 and 5 share the derived allele, 1:\n",
    "dmel_array[13]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmel_array[19]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 13:</b> Tweak the above code to simulate a sample of 50kb of sequence sampled for three modern humans. Assume a recombination rate of 1.2 cM/Mb and a mutation rate of $2.5 \\times 10^{-8}$. Compare the <i>D. melanogaster</i> and <i>Homo</i> simulations using the summary table of the ts. How do toy treesequences for the two taxa differ?\n",
    "    <div/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 14:</b> Define the haplotype array for the <i>Homo</i> simulation as a new variable. Use the function below (I borrowed this code from <a href=\"http://alimanfoo.github.io/2017/02/14/mendelian-transmission.html\">A Miles</a>) to plot the <i>Homo</i> and the <i>D. melanogaster</i> alignment! Is it possible to see any treesequence structure in these haplotypes?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAABECAYAAACF4e8fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAACkklEQVR4nO3ZsUuUYQDH8Ve9xKKrtIJCwggco73aijb/iJrbIhragoYIt7agwNXRqYamgtZodMgjpCDD7CLF1Lcp6IZy6bknfnw+yzO8z/C74YXv3Y20bds2AADBRmsPAAAoTfAAAPEEDwAQT/AAAPEEDwAQT/AAAPE6+124/Ohls7qxNYwtQ/fu7tXm/ZPZ2jOKOXNjuTm3eL/2jCKevX7czM6vNMu3ztaeUsTs/Erz9sGr2jOKOX/nUvPm6e3aM4q5cP1hs/P8Y+0ZxXSunWoWFhZqzyhi7spiMzm91KyvztWeUsTk9FJz88Va7RlFTE2MNvcuTv3x+b7Bs7qx1fTWN//pqP/Jztde7QlF9b6t155QxI+13sCZaPtL5heNX7b7n2tPKOv7bu0FRfX7/doTitjb/TBwJvq0uVd7QhX+0gIA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4gkeACCe4AEA4nX2uzB9dGIYO6rpHJmpPaGomcOTtScUceDEzMCZaPxY9rs33j1ee0JZh8ZqLyiq2+3WnlDE6NjpgTPRyYOZv3VMTfz9c420bdsOaQsAQBWZmQcA8BvBAwDEEzwAQDzBAwDEEzwAQDzBAwDE+wkXsFJol2wNmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x72 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#This creates a colour palette\n",
    "palette = sns.color_palette('colorblind');\n",
    "sns.palplot(palette);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sequence(painting, title=None):\n",
    "    # set figure height depending on number of haplotypes\n",
    "    fig, ax = plt.subplots(figsize=(20, .4 * painting.shape[1]))\n",
    "    \n",
    "    # map painting codes onto colours\n",
    "    cmap = mpl.colors.ListedColormap([\n",
    "        palette[0],  # 0 = an\n",
    "        palette[3],  # 1 = allele inherited from first parental haplotype\n",
    "    ])\n",
    "    # plot painting\n",
    "    ax.pcolormesh(painting.T, cmap=cmap, vmin=0, vmax=1)\n",
    "    # tidy up axes\n",
    "    ax.set_yticks(np.arange(painting.shape[1]) + .5)\n",
    "    ax.set_yticklabels(range(painting.shape[1]))\n",
    "    ax.set_ylabel('Sample haplotypes')\n",
    "    ax.set_xlabel('Variants')\n",
    "    if title:\n",
    "        ax.set_title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIgAAACyCAYAAAApgvQbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfCElEQVR4nO3dfXQU1f3H8c+GJJIIlMcEsMiDxqhYgSpIgj8RkCDE8JCmCEEQKMdSFSRVKYUQ0YoCBVMFLVKBCioSlADhQcHyIDFRClIoT4poAmhOSCUgS0JINvv7w+O2tIZZk5nZdff9OsdzMrMz937v3Lmzw9c7sw632+0WAAAAAAAAglaIrwMAAAAAAACAb5EgAgAAAAAACHIkiAAAAAAAAIIcCSIAAAAAAIAgR4IIAAAAAAAgyJEgAgAAAAAACHIkiAAAAAAAAIJcqK8DqMkXT9+uqtKTl90moekrNkVjnc2nxxluEwjt9BVvjq9k7jH2tk4z+eocMbOt3rbBF33qLX+OzVtmtyEQxqCZbTW77/35+Przdcns89dM/hybt/y57yXv4vPVee7v48sb/nxvGwjXTF/ce6F2fHU/Egj8+Vr4Y+/Xq35SX7kP317j536bIKoqPanKfxVedptCR7lN0VjHqI1SYLTTV7w5vpK5x9jbOs3kq3PEzLZ62wZf9Km3/Dk2b5ndhkAYg2a21ey+9+fj68/XJbPPXzP5c2ze8ue+l7yLz1fnub+PL2/4871tIFwzfXHvhdrx1f1IIPDna2Gg9yuPmAEAAAAAAAQ5EkQAAAAAAABBjgQRAAAAAABAkCNBBAAAAAAAEOQsfUn1yJEjdfr0aYWGflvNU089pU6dOllZJQAAAAAAAH4gyxJEbrdbBQUF2rZtmydBBAAAAAAAAP9j2SNmn3/+uSRp7NixGjhwoF577TWrqgIAAAAAAEAdWDa155tvvlFcXJymT5+uyspKjRo1Su3bt1ePHj2sqhIAAAAAAAC1YFmCqEuXLurSpYtnOSUlRTt27CBBBAAAAAAA4Gcse8Rs9+7dys/P9yy73W7eRQQAAAAAAOCHLEsQnTt3TnPmzFFFRYWcTqeys7PVt29fq6oDAAAAAABALVk2padXr17at2+fBg8erOrqaqWmpl7yyBkAAAAAAAD8g1cJIqfTqQYNGujQoUP69NNPlZiYqLCwMMP9Jk2apEmTJtU1RgAAAAAAAFjIMEH0/PPP6/jx43r00Uc1btw4XXvttfr73/+umTNn2hEfAAAAAAAALGb4DqIdO3bo6aef1ubNm5WYmKhly5bpyJEjdsQGAAAAAAAAG3j1kuqIiAjl5eWpe/fukqSLFy9aGhQAAAAAAADsY5ggatKkiWbMmKEDBw4oPj5ec+fOVVRUlB2xAQAAAAAAwAaGCaLZs2crKipKL7/8siIiIuRwODR79mw7YgMAAAAAAIANDBNEzZs3V3Jysk6fPi2Xy6Xhw4erefPmdsQGAAAAAAAAGzjcbrf7chts375dM2bMUEhIiN58803169dPf/zjH3XXXXdZGlj7me+psLS8zuUcKBloQjT/dlOLdaaW5w1ftcHber0pz+w2eMvstppZp5nM7CsryvMFf25DIFyXzObNMfFVO31xLvnq/DWzH/z5uu/P14dA4KvvJG8FQr8Gy32LtwKhT33B38eqmQLhHjgQ7kcC4Zzz5+uNN8cjrHlbxcwrqPFzwxlEL774orKystSoUSNFRUXpjTfe0AsvvPCDAgUAAAAAAID/MkwQuVyuS15KfcMNN8jhcFgaFAAAAAAAAOxjmCCKiIjQV1995UkK7d69W1dccYXlgQEAAAAAAMAeoUYbPPbYYxo7dqxKSkp07733qqCgQPPnz7cjNgAAAAAAANjAMEHUpUsXZWVlae/evaqurlanTp3UtGlTO2IDAAAAAACADQwTRJK0d+9e7dy5U/Xq1VNkZKRuu+02q+MCAAAAAACATQzfQTR//nzNmjVLDRs2VP369ZWRkaFly5Z5XcHs2bM1ZcqUOgUJAAAAAAAA6xjOIFq3bp1Wr16thg0bSpLGjh2rYcOGadSoUYaF5+fnKzs7W3feeWedAwUAAAAAAIA1DGcQNW7cWFdeeaVnuVGjRoqMjDQs+MyZM8rMzNT48ePrFiEAAAAAAAAsZTiD6JZbbtGDDz6oe++9V/Xq1dO6devUunVrbd68WZKUkJDwvftlZGQoLS1NRUVF5kYMAAAAAAAAUxkmiA4ePChJWrJkySXrly9fLofD8b0JolWrVqlVq1aKi4vT6tWrTQoVAAAAAAAAVjBMEKWnpys2NvYHFbpx40aVlJRo0KBBOnv2rMrKyvTMM89o6tSptQ4UAAAAAAAA1jBMEE2cOFFNmjRRamqq7r77boWHhxsWunTpUs/fq1ev1q5du0gOAQAAAAAA+CnDl1S/++67euSRR7Rt2zb16dNHc+bMUWFhoR2xAQAAAAAAwAaGM4gkKS4uTnFxccrPz1d6erqWLVum+Ph4TZ8+XW3atLnsvsnJyUpOTjYlWAAAAAAAAJjPMEHkdDqVk5OjVatWqby8XPfff7+GDBmi999/Xw8++KBycnLsiBMAAAAAAAAWMUwQ9ezZU3FxcXrssccUHx/vWZ+YmKg333zT0uAAAAAAAABgPcME0cKFC9W1a9dL1uXl5Sk+Pl7Lly+3LDAAAAAAAADYo8YE0aFDh+R2u5WRkaF58+bJ7XZLkqqqqjRjxgxt3rzZtiABAAAAAABgnRoTRCtWrNAHH3ygU6dO6eGHH/73DqGh6tu3ry3BAQAAAAAAwHoO93dTg2qQmZmptLQ0u+LxaD/zPRWWll92mwMlAw3LuanFOrNC+kG8iU3yXXxmoh9qV6e3/PkcMbOt3rbTV2PLF231lj/HZiZ/73tv6g2ENphdp9l80Q++OEe8xZiuXb1m9723fNFfgdAGswXL96q3fHXNNLNOs/HvH9/X+UPqNZM/fyd5o22TCH0x7a4aPzd8B1FaWpo2btyonTt3qrKyUrfffrsGDx5sWoAAAAAAAADwrRCjDZYsWaKXX35ZsbGx6tixo5YuXaqXXnrJjtgAAAAAAABgA8MZRNnZ2VqxYoUaNGggSUpJSdHQoUP14IMPWh4cAAAAAAAArGc4g0iSJzkkSQ0bNlRoqGFeCQAAAAAAAD8Shgmiq666Sq+++qoqKytVWVmpv/71r2rdurUdsQEAAAAAAMAGhgmiJ598Uu+99546d+6szp07a/PmzcrIyLAjNgAAAAAAANjA8Fmx6OhoLV++XOXl5aqurtaVV15pR1wAAAAAAACwSY0JovHjx192x4ULFxoW/vzzz+vdd9+Vw+FQSkqKxowZ88MjBAAAAAAAgKVqTBD169evTgXv2rVLH374odatW6eqqioNGDBAPXv2VIcOHepULgAAAAAAAMxVY4JoyJAhnr9LS0u1e/duhYSEqFu3bmrYsKFhwd26ddOyZcsUGhqq4uJiuVwuRUZGmhM1AAAAAAAATGP4kuotW7YoISFBr776ql555RX17dtXH374oVeFh4WF6YUXXlBiYqLi4uIUHR1d54ABAAAAAABgLsMEUWZmpl577TW99tprWrFihRYvXqzZs2d7XcHEiROVn5+voqIiZWVl1SlYAAAAAAAAmM8wQVS/fn3FxsZ6ljt27CiHw2FY8LFjx3T48GFJUkREhBISEvTJJ5/UIVQAAAAAAABYwTBBdMcdd2jRokUqKytTRUWFVq5cqZiYGJ09e1Znzpypcb+TJ08qPT1dFy9e1MWLF/W3v/1Nt9xyi5mxAwAAAAAAwAQ1vqT6O3/5y1/kcrn03HPPXbJ+7dq1cjgcnllC/61nz57av3+/Bg8erHr16ikhIUGJiYnmRA0AAAAAAADTGCaIDh48WOvCJ0yYoAkTJtR6fwAAAAAAAFjPMEF08eJF7dixQ+fPn5ckuVwuHT9+XGlpaZYHBwAAAAAAAOsZJojS0tJ04sQJlZSU6MYbb9S+ffvUrVs3O2IDAAAAAACADQxfUn348GGtXr1affr00dSpU7VixQqdPXvWjtgAAAAAAABgA8MEUVRUlEJDQ9WuXTt9+umniomJ0blz5+yIDQAAAAAAADYwTBBFRkYqJydH119/vTZt2qRPPvlEZWVldsQGAAAAAAAAGzjcbrf7chsUFBQoKytLjz/+uNLS0vTBBx8oLS1NqamplgbWfuZ7Kiwtr3M5B0oGmhDNv93UYp2p5XnDV23wtl5vyjO7Dd4yu61m1mkmM/vKivJ8wZ/bEAjXJbN5c0x81U5fnEu+On/N7Ad/vu778/UhEPjqO8lbgdCvwXLf4q1A6FNf8PexaqZAuAcOhPuRQDjn/Pl6483xCGveVjHzCmr83PAl1e3atdPkyZMlSX/605+8Dg4AAAAAAAA/DoYJoj179mjBggX6+uuv9Z+TjXJyciwNDAAAAAAAAPYwTBBNnz5dQ4cO1Q033CCHw2FHTAAAAAAAALCRYYIoPDxco0ePtiEUAAAAAAAA+ILhr5h16NBB//znP+2IBQAAAAAAAD5Q4wyipKQkSdL58+c1fPhwtWnTRqGh/96cdxABAAAAAAAEhhoTRNOnT69z4QsWLNCmTZskST179vT8GhoAAAAAAAD8R40Jom7dutWp4Ly8POXm5io7O1sOh0Pjxo3Tli1b1Ldv3zqVCwAAAAAAAHMZvqS6tlq0aKEpU6YoPDxcknTNNdfoq6++sqo6AAAAAAAA1JJlCaKYmBjP3wUFBdq0aZNWrFhhVXUAAAAAAACoJcNfMfvON998U6sKjh49qrFjx2ry5Mlq165drcoAAAAAAACAdQwTRJ9//rkGDBigxMREFRcXq3///jp27JhXhe/Zs0ejR4/Wo48+qiFDhtQ5WAAAAAAAAJjPMEH09NNPa9q0aWrWrJmio6N13333KSMjw7DgoqIiPfTQQ5o7d64SExNNCRYAAAAAAADmM0wQnTlzRj169PAsjxgxQk6n07DgxYsXq6KiQrNmzdKgQYM0aNAg3kEEAAAAAADgh7x6SXVFRYUcDockqaSkRNXV1Yb7pKenKz09vW7RAQAAAAAAwHKGCaLU1FT96le/0tdff6158+Zpw4YNGjdunB2xAQAAAAAAwAaGCaKUlBS1bdtW27dvV1VVlf7whz9c8sgZAAAAAAAAftxqTBCdOXPG83dMTIxiYmIu+axx48ZWxgUAAAAAAACb1Jgg6t69uxwOh9xut2fdd8sOh0OHDx+2JUAAAAAAAABYq8YE0ZEjR+yMAwAAAAAAAD5i+A4il8ulN998U7m5uapXr5569+6t5ORkO2IDAAAAAACADRzu/3yG7HvMmDFDx44d06BBg+R2u5Wdna2uXbsqLS3N0sDaz3xPhaXll93mQMlAw3JuarHOrJC8rvOHMDs+b3jbBl/EFgj8+RzxVd+beUx8FRvj4VL+fNx8NQb9+TspEK4j3jKzH8w+l7zhz7F5y5/7XuJ6HogCoe8DoQ3e8sV9YTAdX2/54r7FbP78PW0mM8/ftk0i9MW0u2r83HAGUV5enjZs2KCwsDBJ0sCBAzVw4EDLE0QAAAAAAACwR4jRBk2bNpXL5fIsOxwONWrUyNKgAAAAAAAAYB/DGUTXX3+9UlNTlZycrHr16mnjxo1q0qSJli5dKkkaM2aM5UECAAAAAADAOoYJooqKCsXGxurgwYOSpJ/+9KeSpE8//dTayAAAAAAAAGALwwTRs88+a0ccAAAAAAAA8BHDBNFHH32kRYsW6ezZs5esf+uttywLCgAAAAAAAPYxTBClp6dr5MiRuvrqq39w4U6nU8OGDdPChQs9j6YBAAAAAADAvxgmiJo1a6ZRo0b94IL37dun9PR0FRQU1CYuAAAAAAAA2MTwZ+579+6t119/XcePH9dXX33l+c9IVlaWnnjiCUVFRZkSKAAAAAAAAKxhOIOotLRUzz33nCIiIjzrHA6HPv7448vuN3PmzLpHBwAAAAAAAMsZJoi2bdum3NxcNW/e3I54AAAAAAAAYDPDR8yaNWumpk2b2hELAAAAAAAAfMBwBtF1112n1NRU9erVS+Hh4Z71Y8aMsTQwAAAAAAAA2MMwQXThwgW1b9+eXyMDAAAAAAAIUIYJomeffbZOFWzdurVO+wMAAAAAAMBahgmivXv3atGiRSorK5Pb7VZ1dbVOnjyp7du32xAeAAAAAAAArGb4kur09HR16dJFTqdTSUlJatCggRISEuyIDQAAAAAAADYwnEHkcDj0wAMPqLS0VB06dFBSUpJ+8Ytf2BEbAAAAAAAAbGA4g+jKK6+UJF199dU6evSo6tevr5AQw90AAAAAAADwI2E4g+jmm2/WpEmT9Mgjj+jXv/61CgoKFBpquBsAAAAAAAB+JAynAk2dOlWjR49W+/btNW3aNFVXV2vevHl2xAYAAAAAAAAbONxut7umD91ut1wul0JDQ+V0OpWXl6frrrtO7dq1szyw9jPfU2FpeZ3LOVAy0IRo/u2mFutMLc8bvmqDt/V6U57ZbfCW2W01s04zmdlXVpTnC/7chkC4LpnNm2Piq3b64lzy1flrZj/483Xfn68PgcBX30neCoR+DZb7Fm8FQp/6gr+PVTMFwj1wINyPBMI558/XG2+OR1jztoqZV1Dj5zXOIPrss8/Up08f7dy5UxcuXNAvf/lLZWZmauTIkfrggw9qFTAAAAAAAAD8T40Jojlz5mjSpEnq1auXNmzYIEnasGGDsrKyNH/+fNsCBAAAAAAAgLVqTBAVFRVp4MBvpyh99NFH6tOnj0JCQtSqVSs5nU7bAgQAAAAAAIC1akwQ/edP2e/du1ddu3b1LFdUVFgbFQAAAAAAAGxT4+/V/+QnP9GRI0fkdDpVUlLiSRB9/PHHio6Oti1AAAAAAAAAWKvGBNFvf/tbjR49Wk6nU4899pgiIyO1ePFiLVy4UC+++KKdMQIAAAAAAMBCNSaIOnfurPfff18XLlxQo0aNJEldunTRqlWrvP6Z+5ycHP35z39WVVWV7r//fo0YMcKUoAEAAAAAAGCeGhNEkhQeHq7w8HDP8s9//nOvCy4uLlZmZqZWr16t8PBwDRs2TLfddpuuvfba2kcLAAAAAAAA09X4kuq6ysvLU/fu3dW4cWNFRkaqX79+euedd6yqDgAAAAAAALVkWYLo1KlTatGihWc5KipKxcXFVlUHAAAAAACAWrIsQVRdXS2Hw+FZdrvdlywDAAAAAADAP1iWIGrZsqVKSko8yyUlJYqKirKqOgAAAAAAANSSZQmi+Ph45efn6/Tp0yovL9fmzZt1xx13WFUdAAAAAAAAaumyv2JWF9HR0UpLS9OoUaNUWVmplJQU3XzzzVZVBwAAAAAAgFqyLEEkSUlJSUpKSrKyCgAAAAAAANSRZY+YAQAAAAAA4MeBBBEAAAAAAECQI0EEAAAAAAAQ5EgQAQAAAAAABDlLX1JdF1f9pL4p5YS525pSznfaNokwtTxv+KoN3tbrTXlmt8FbZrfVzDrNZGZfWVGeL/hzGwLhumQ2b46Jr9rpi3PJV+evmf3gz9d9f74+BAJffSd5KxD6NVjuW7wVCH3qC/4+Vs0UCPfAgXA/EgjnnD9fb7w5HqFNfnrZzx1ut9ttVkAAAAAAAAD48eERMwAAAAAAgCBHgggAAAAAACDIkSACAAAAAAAIciSIAAAAAAAAghwJIgAAAAAAgCBHgggAAAAAACDIkSACAAAAAAAIciSIAAAAAAAAghwJIgAAAAAAgCDnVwminJwcDRgwQAkJCXr99dd9HQ5guwULFigxMVGJiYmaM2eOJCkvL09JSUlKSEhQZmamjyME7Dd79mxNmTJFEuMBwWnr1q1KTk5W//799fTTT0tiLCB4rV271nOvNHv2bEmMBwQfp9Ope+65RydPnpRU8xg4fPiwkpOT1a9fP02bNk1VVVW+Chk/En6TICouLlZmZqbeeOMNrVmzRitXrtRnn33m67AA2+Tl5Sk3N1fZ2dlas2aNDh48qPXr12vq1Kl66aWXtHHjRh04cEA7duzwdaiAbfLz85WdnS1JunDhAuMBQefEiRN64okn9NJLL2ndunU6dOiQduzYwVhAUCovL9fMmTO1fPlyrV27Vrt379bWrVsZDwgq+/bt0/Dhw1VQUCDp8vdHjz/+uDIyMvTuu+/K7XYrKyvLh5Hjx8BvEkR5eXnq3r27GjdurMjISPXr10/vvPOOr8MCbNOiRQtNmTJF4eHhCgsL0zXXXKOCggK1bdtWbdq0UWhoqJKSkhgXCBpnzpxRZmamxo8fL0nav38/4wFBZ8uWLRowYIBatmypsLAwZWZmKiIigrGAoORyuVRdXa3y8nJVVVWpqqpKDRo0YDwgqGRlZemJJ55QVFSUpJrvj7788ktduHBBnTt3liQlJyczNmAo1NcBfOfUqVNq0aKFZzkqKkr79+/3YUSAvWJiYjx/FxQUaNOmTbrvvvv+Z1wUFxf7IjzAdhkZGUpLS1NRUZGk7/+eYDwg0BUWFiosLEzjx49XUVGR7rzzTsXExDAWEJQaNGigRx55RP3791dERIS6du3KdwOCzsyZMy9ZrmkM/Pf6Fi1aMDZgyG9mEFVXV8vhcHiW3W73JctAsDh69KjGjh2ryZMnq02bNowLBKVVq1apVatWiouL86zjewLByOVyKT8/X88884xWrlyp/fv368SJE4wFBKUjR47o7bff1rZt27Rz506FhISooKCA8YCgVtP9EfdNqA2/mUHUsmVL7d6927NcUlLimTYHBIs9e/Zo4sSJmjp1qhITE7Vr1y6VlJR4PmdcIFhs3LhRJSUlGjRokM6ePauysjJ9+eWXqlevnmcbxgOCQfPmzRUXF6emTZtKku666y698847jAUEpdzcXMXFxalZs2aSvn1kZvHixYwHBLWWLVt+778X/nv9v/71L8YGDPnNDKL4+Hjl5+fr9OnTKi8v1+bNm3XHHXf4OizANkVFRXrooYc0d+5cJSYmSpI6deqkL774QoWFhXK5XFq/fj3jAkFh6dKlWr9+vdauXauJEyeqd+/eeuWVVxgPCDq9evVSbm6uvvnmG7lcLu3cuVN33303YwFB6frrr1deXp7Kysrkdru1detW7pUQ9GoaA1dddZWuuOIK7dmzR9K3vwDI2IARv5lBFB0drbS0NI0aNUqVlZVKSUnRzTff7OuwANssXrxYFRUVmjVrlmfdsGHDNGvWLE2YMEEVFRXq2bOn7r77bh9GCfjOFVdcwXhA0OnUqZPGjRun1NRUVVZWqkePHho+fLg6dOjAWEDQuf3223Xo0CElJycrLCxMP/vZzzRhwgT16NGD8YCgdbn7o7lz5yo9PV1Op1MdO3bUqFGjfBwt/J3D7Xa7fR0EAAAAAAAAfMdvHjEDAAAAAACAb5AgAgAAAAAACHIkiAAAAAAAAIIcCSIAAAAAAIAgR4IIAAAAAAAgyJEgAgAAQWHkyJFatGjR/6xfsmSJfvOb33hVxooVK763DG/t379fGRkZtd4fAADAKiSIAABAUEhNTdXbb7/9P+uzsrI0YsQIr8oYPny4HnjggVrH8Nlnn6m4uLjW+wMAAFjF4Xa73b4OAgAAwGpVVVXq1auXMjMzdeutt0qSdu3apenTp2vw4MHaunWrLly4oPLycv3ud79T3759NX/+fP3jH//QqVOnFBsbq7Zt26q0tFQZGRnatm2bXn75ZV28eFGnT5/W4MGDNWnSJH300UfKzMxUmzZtdPToUVVVVenJJ59U69atNXz4cJ07d04JCQlKT0/X73//exUWFiokJEQdO3bUU089pZAQ/v8dAACwH3cgAAAgKISGhmro0KF66623POtWrlyppKQk5efna/ny5crJyVFaWppeeOEFzzZffvmlsrOzNXfuXM86t9utJUuWaNasWVq9erVWrlypRYsW6fTp05K+fZRs7NixWrNmjZKTk5WZmalWrVpp4sSJuvXWW/Xss89qy5YtOn/+vNauXeuJ6cSJEzYdDQAAgEuF+joAAAAAuwwdOlSJiYlyOp2qqqpSbm6uZsyYoZSUFOXk5KiwsFD79u3T+fPnPft07txZoaGX3jI5HA4tXLhQ27dv1/r163Xs2DG53W6Vl5dLklq3bq0bbrhBknTjjTcqOzv7f2K55ZZblJmZqZEjRyo+Pl7333+/2rZta2HrAQAAasYMIgAAEDSio6MVHx+vjRs3as2aNerXr5+OHz+ue++9V06nUz169NC4ceMu2ScyMvJ/yikrK9OQIUN08OBB3XjjjZo8ebJCQ0P13ZP79evX92zrcDj0fU/0t2nTRlu2bNEDDzwgp9OpMWPGaOvWrSa3GAAAwDvMIAIAAEFlxIgRWrBggc6dO6fZs2crPz9fN910k8aMGSOXy6Unn3xSLpfrsmUUFhbK6XRq0qRJCg8P15o1a3Tx4kVVV1dfdr969eqpqqpKkvTGG29oz549mjt3rv7v//5PX3/9tQ4dOqTevXub1lYAAABvMYMIAAAEldtuu01nzpxRgwYNFBsbq3vuuUelpaXq37+/BgwYoMjISJ09e1ZOp7PGMmJjY3XnnXeqf//+6t+/v7Zt26Zrr71WhYWFl627c+fOOnHihB5++GENHjxYLpdLAwYMUHJyss6dO6eRI0ea3VwAAACv8CtmAAAAAAAAQY4ZRAAAAAAAAEGOBBEAAAAAAECQI0EEAAAAAAAQ5EgQAQAAAAAABDkSRAAAAAAAAEGOBBEAAAAAAECQI0EEAAAAAAAQ5EgQAQAAAAAABLn/B3aDy9Y8dCjxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x172.8 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "homo_array = np.zeros((homo_ts_muts.num_sites, homo_ts_muts.num_samples))\n",
    "for var_index, var in enumerate(homo_ts_muts.variants()):\n",
    "    homo_array[var_index] = var.genotypes\n",
    "plot_sequence(homo_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIkAAACyCAYAAADGQJ8lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArT0lEQVR4nO3dfXwV1Z3H8e+FJEpEqwhBVApaERVWoD4gYEWwBCEGMLIoICiUl0ttBWItooaIT+Vh0bRKt0hFWsUiSEEIgoUtYHmqVkplBVSKJaCwGOVBAiGPs3/0lbuTYc6cuUnuDZbP+/Xal3dmzsPvnDnnzPTszSXiOI4jAAAAAAAAnNYa1HcAAAAAAAAAqH9sEgEAAAAAAIBNIgAAAAAAALBJBAAAAAAAALFJBAAAAAAAALFJBAAAAAAAALFJBAAAAAAAAElJ9R2AyT+evlE9I5MkSSsPjlJ6k5eqXfc7V3VektKbvHTSZ7/03jzussPUERRTUH53XVW8dYbNH7Zdthi97fbGZkrj919TO233xe9a0P0wxW8q31uPt3/C9qFfPlt/x9oH7rKCzoUZ37b7aUvvzefXXlOf2sarrRzT/PSLOWzbw44lU5zevgk7H8PM7bDl+eVz5/WLubZ1hpnf7mtB/RgmfdD8CTsGbOV4ywxag0xMcyXMmLKV6zcHbG0MitFvnrnj9cYey9wJOw9jWfODhFlLgtbWoOdW2DaFyWNqj23++Knp/Az7nmCaH2Fj8Cs37Ph1p7etB6bYw74zBN2jWNbcmq6bYfs97Jpra5u3bG/5prkQJo6gOWSK0dQ/pv7z6yu//jK1OZZ5EUuasGu4LV8sY76mMcdyzTQuwjzbw8QR5npYYesNmj9h3wWCynKXFya+oOdCLO8GYfswKM4w60xQ3EHPO9t6Z4rRLZY12S8eU3lhxor3uK7mmh/b2muqvyb9ttqZpEty1huvS6fwJlH5oc9UECmWJJV9WRD9XMXvXNV5SSqIFJ/02S+9N4+77DB1BMUUlN9dVxVvnWHzh22XLUZvu72xmdL4/dfUTtt98bsWdD9M8ZvK99bj7Z+wfeiXz9bfsfaBu6ygc2HGt+1+2tJ78/m119SntvFqK8c0P/1iDtv2sGPJFKe3b8LOxzBzO2x5fvncef1irm2dYea3+1pQP4ZJHzR/wo4BWzneMoPWIBPTXAkzpmzl+s0BWxuDYvSbZ+54vbHHMnfCzsNY1vwgYdaSoLU16LkVtk1h8pjaY5s/fmo6P8O+J5jmR9gY/MoNO37d6W3rgSn2sO8MQfcoljW3putm2H4Pu+ba2uYt21u+aS6EiSNoDpliNPWPqf/8+sqvv0xtjmVexJIm7BpuyxfLmK9pzLFcM42LMM/2MHGEuR5W2HqD5k/Yd4GgstzlhYkv6LkQy7tB2D4MijPMOhMUd9DzzrbemWJ0i2VN9ovHVF6YseI9rqu55se29prqr0m/lTufGa9V4c/NAAAAAAAAwCYRAAAAAAAA2CQCAAAAAACA2CQCAAAAAACA4vzD1cOGDdPBgweVlPTPap588kl16NAhnlUCAAAAAACgBuK2SeQ4jnbv3q01a9ZEN4kAAAAAAABwaorbn5t9+umnkqSRI0eqX79+mjt3bryqAgAAAAAAQC3F7Ss+X3/9tbp06aKJEyeqrKxMw4cP1yWXXKJu3brFq0oAAAAAAADUUNw2iTp16qROnTpFjwcOHKh33nmHTSIAAAAAAIBTUNz+3Oz999/Xpk2boseO4/DbRAAAAAAAAKeouG0SHT16VNOmTVNJSYmKioq0ePFi9erVK17VAQAAAAAAoBbi9tWeHj166IMPPtCAAQNUWVmpIUOGVPvzMwAAAAAAAJw6Qm0SFRUVqXHjxtq+fbs++eQTZWRkKDk52Zpv3LhxGjduXG1jBAAAAAAAQJxZN4l+8YtfaM+ePfrJT36iUaNG6bLLLtNf/vIXPfPMM4mIDwAAAAAAAAlg/U2id955R08//bRWrlypjIwMvfLKK/roo48SERsAAAAAAAASJNQPVzdq1EgbN27UDTfcIEkqLS2Na1AAAAAAAABILOsm0XnnnadJkybpww8/VNeuXTV9+nSlpaUlIjYAAAAAAAAkiHWTaOrUqUpLS9OLL76oRo0aKRKJaOrUqYmIDQAAAAAAAAli3SRq2rSpsrKydPDgQVVUVGjw4MFq2rRpImIDAAAAAABAglj/dbO1a9dq0qRJatCggV5//XVlZGToP//zP/X9738/roGlN3lJH37SK3r8YWE/tW+2NPrf9s2WVkv/YWE/Sap23u+zLZ372HveHYe7rKq07vPeY790QXVXfXbX5y3LncZ0zpvXm96vDlO73ddM5Xr7wO9+edP4le9tk18sXt7YvOn82u9Xpt89MNXljtuv3/zyePPb8nnLCJoD7utB98zLr/1B48AdS1hVaf3GplvQmPCbZ+5Ygv7rrcM0R4JiMJVnitUkzD0PKsO2HvrVY5v77rLc7THNDXcsJqb1KSh9UPzez94x4I7Jtv4ExWKbM2HHsCmfO53fGu+t329d9Z735jWlD3o+mWIL6uew7TCN5zBj3JvWtu6Eef6Z2uQXn2l99Gu/O63t2RrL+ulNH+bdwJ3OdH9ief54y3PXU5P8fvcz6DntTWPid49ijdHGNNb98vv1f5jnWqyfbWuhOxa/doRtr+162LUiqO1+ZZvmbNjnrikWv/psafzK9Jvf3rJMa45prQx65prW1bD31PTO7hdPmDK8dZv6w1ufae4ErTOm93jTmmZqQ9A7gqkt3vXKmybMXPe7x2HHr03Y9wz3db/56Bdb0Hu1V9j1zhSjt01h3jP8zvvFYxLmuVjTe2R7jwpTX9A8NdVpKif98lX6hyXmBpbr+uUvf6kFCxbonHPOUVpamn73u9/p+eeft2UDAAAAAADAN4h1k6iioqLaD1VfeeWVikQicQ0KAAAAAAAAiWXdJGrUqJH27dsX3Rh6//33dcYZZ8Q9MAAAAAAAACSO9TeJHnroIY0cOVKFhYW68847tXv3br3wwguJiA0AAAAAAAAJYt0k6tSpkxYsWKAtW7aosrJSHTp0UJMmTRIRGwAAAAAAABLEukkkSVu2bNG6devUsGFDpaamqnPnzvGOCwAAAAAAAAlk/U2iF154QVOmTNHZZ5+tM888U7m5uXrllVdCVzB16lRNmDChVkECAAAAAAAgvqzfJFq6dKkWLVqks88+W5I0cuRI3XXXXRo+fLi18E2bNmnx4sW6+eabax0oAAAAAAAA4sf6TaJzzz1XZ511VvT4nHPOUWpqqrXgw4cPKy8vT6NHj65dhAAAAAAAAIg76zeJrrnmGt1///2688471bBhQy1dulQXXnihVq5cKUlKT0/3zZebm6vs7Gzt37+/biMGAAAAAABAnbNuEm3btk2S9PLLL1c7/+qrryoSifhuEr3xxhtq0aKFunTpokWLFtVRqAAAAAAAAIgX6yZRTk6O2rZtG1Ohy5cvV2Fhofr3768jR47o+PHj+tnPfqZHH320xoECAAAAAAAgfqybRGPGjNF5552nIUOG6NZbb1VKSoq10Dlz5kQ/L1q0SO+99x4bRAAAAAAAAKcw6w9X/+EPf9DYsWO1Zs0a3XLLLZo2bZoKCgoSERsAAAAAAAASxPpNIknq0qWLunTpok2bNiknJ0evvPKKunbtqokTJ6ply5aBebOyspSVlVUnwQIAAAAAACA+rJtERUVFys/P1xtvvKHi4mLdc889uv322/WnP/1J999/v/Lz8xMRJwAAAAAAAOLIuknUvXt3denSRQ899JC6du0aPZ+RkaHXX389rsEBAAAAAAAgMaybRDNnztR1111X7dzGjRvVtWtXvfrqq3ELDAAAAAAAAIlj3CTavn27HMdRbm6unn32WTmOI0kqLy/XpEmTtHLlyoQFCQAAAAAAgPgybhLNmzdPGzZs0BdffKEf//jH/58hKUm9evVKSHAAAAAAAABIDOMm0VNPPSVJysvLU3Z2dsICcmvfbGm14w8L+0XPuT+b0lels5Xlp+q6O7+t/qpzQXH5leGNserYL253+qB2VKVzX/e2wdYOvzbY+sAdv7d93vZ4P/ulc5dh6wtv/N7+NJVj6hNv2X7l2vIFjaGaCjMvgs4F3Q/vvTPV7x2DJn59Ebb9tnljuu433sPE4S3Pe85bflC8tbnX3vvgPVfVdm+fmNoR1M6azoWg8t3x+q2H3nL95lTQnAlaE0ztiqegddD2jApa12J51tnGtN/zIGjeuuv1W8+9c8F2zR2DX9x+bTK11XQu6LnpPme6B3WxNleV79emsO8s7jL8/ustK+xzM2zdpnj9xktN1zfbddu6EzZ20/M+1viCxq07T5j3zaA1Lag9fjEGPSeqjsPc/6Cx5heLbX33xuxNFzYmvzR1xbSGxprftBaa6jK991ddM6X3K6+2a5bfffSr27R2e9Pb3h3CjgFTXd6Y/JjWSb9noLvsmrybmuax33ptymd7nnrPh0lnev6FfUYG/W8GP371B70bmmILc0+97Y1l/Aetu2HvXW3Ynqdh4ja9C8T6rGzfbKk+PjhK0u7Auq2/SZSdna3ly5dr3bp1Kisr04033qgBAwbYsgEAAAAAAOAbpIEtwcsvv6wXX3xRbdu2Vbt27TRnzhz913/9VyJiAwAAAAAAQIJYv0m0ePFizZs3T40bN5YkDRw4UIMGDdL9998f9+AAAAAAAACQGNZvEkmKbhBJ0tlnn62kJOveEgAAAAAAAL5BrJtEF110kX7729+qrKxMZWVl+s1vfqMLL7wwEbEBAAAAAAAgQaybRE888YT++7//Wx07dlTHjh21cuVK5ebmJiI2AAAAAAAAJIj178aaN2+uV199VcXFxaqsrNRZZ52ViLgAAAAAAACQQMZNotGjRwdmnDlzprXwX/ziF/rDH/6gSCSigQMHasSIEbFHCAAAAAAAgLgzbhL17t27VgW/9957+vOf/6ylS5eqvLxcffv2Vffu3XXppZfWqlwAAAAAAADUPeMm0e233x79fOjQIb3//vtq0KCBrr/+ep199tnWgq+//nq98sorSkpK0oEDB1RRUaHU1NS6iRoAAAAAAAB1yvrD1atWrVJ6erp++9vf6qWXXlKvXr305z//OVThycnJev7555WRkaEuXbqoefPmtQ4YAAAAAAAAdc+6SZSXl6e5c+dq7ty5mjdvnmbPnq2pU6eGrmDMmDHatGmT9u/frwULFtQqWAAAAAAAAMSHdZPozDPPVNu2baPH7dq1UyQSsRa8a9cu7dixQ5LUqFEjpaen6+OPP65FqAAAAAAAAIgX6ybRTTfdpFmzZun48eMqKSnR/Pnz1aZNGx05ckSHDx825vvss8+Uk5Oj0tJSlZaW6o9//KOuueaauowdAAAAAAAAdcT4w9VVfv3rX6uiokLPPfdctfNLlixRJBKJflvIq3v37tq6dasGDBighg0bKj09XRkZGXUTNQAAAAAAAOqUdZNo27ZtNS78gQce0AMPPFDj/AAAAAAAAEgM6yZRaWmp3nnnHR07dkySVFFRoT179ig7OzvuwQEAAAAAACAxrJtE2dnZ2rt3rwoLC3XVVVfpgw8+0PXXX5+I2AAAAAAAAJAg1h+u3rFjhxYtWqRbbrlFjz76qObNm6cjR44kIjYAAAAAAAAkiHWTKC0tTUlJSWrdurU++eQTtWnTRkePHk1EbAAAAAAAAEgQ6yZRamqq8vPzdcUVV2jFihX6+OOPdfz48UTEBgAAAAAAgASx/iZRbm6uFixYoJ/+9KdauHCh7r777nr70er2zZb6fq7yYWG/k877pas675fem89bZ1Ueb76q8x8W9gus253GVEbVeXdZtRXU1poytcXbt95+qSumMr395z72i8V039xtcuf1nq9KG5TPFLcpHr/y/eIIuq/uvFWfbX3mLTsonTdWb91VfWKaE35tMx17y3G3O8zYCjMG/eL01mPqc9M9N811bzrTee+98FsDTWPS3W6/+xp2rQzi186gNteUdz1x12UaU7Z2e/OGiddv/TfNYdNYMd0fbx2mttnuo+lZZHvW+Y07b/m2uv3aZ+pvU51BcztovIVlii9oPgStW6b0QeMuTJmm8oPGrXfdDepvv+Mgpr42zasw48rWl2HXk7DjKOwcD5PG1O9+5die2ba4TM8/v/KC1gPTvbfF5xeLaS76jX0b0xwxpQ1aB+LxrusV1E+mdcAdW9DcNM2xMM8yv3sRti1+8YcRdP/9yjOtf2GfhbYxYnumh4nHe82UL+wzw7Q2muoNitd7Leh9MaisMGti0Fi2xViTOei3ZrrLMj1bTM9P29po6j/Ts9b0nhvLmmPr9zBtrMkY9mtr+uWr9A9LvNZNotatW2v8+PGSpJ///Oe25AAAAAAAAPgGsm4Sbd68WTNmzNBXX30lx3Gi5/Pz8+MaGAAAAAAAABLHukk0ceJEDRo0SFdeeaUikUgiYgIAAAAAAECCWTeJUlJSdO+99yYgFAAAAAAAANSXBrYEl156qf7nf/4nEbEAAAAAAACgnhi/SZSZmSlJOnbsmAYPHqyWLVsqKen/k/ObRAAAAAAAAP86jJtEEydOrHXhM2bM0IoVKyRJ3bt3j/4raQAAAAAAADi1GDeJrr/++loVvHHjRq1fv16LFy9WJBLRqFGjtGrVKvXq1atW5QIAAAAAAKDuWX+4uqaaNWumCRMmKCUlRZL0ne98R/v27YtXdQAAAAAAAKiFuG0StWnTJvp59+7dWrFihebNmxev6gAAAAAAAFAL1n/drMrXX39dowp27typkSNHavz48WrdunWNygAAAAAAAEB8WTeJPv30U/Xt21cZGRk6cOCA+vTpo127doUqfPPmzbr33nv1k5/8RLfffnutgwUAAAAAAEB8WDeJnn76aT322GM6//zz1bx5c919993Kzc21Frx//3796Ec/0vTp05WRkVEnwQIAAAAAACA+rJtEhw8fVrdu3aLHQ4cOVVFRkbXg2bNnq6SkRFOmTFH//v3Vv39/fpMIAAAAAADgFBXqh6tLSkoUiUQkSYWFhaqsrLTmycnJUU5OTu2iAwAAAAAAQEJYN4mGDBmiH/zgB/rqq6/07LPP6q233tKoUaMSERsAAAAAAAASxLpJNHDgQLVq1Upr165VeXm5nnrqqWp/fgYAAAAAAIBvPuMm0eHDh6Of27RpozZt2lS7du6558YzLgAAAAAAACSQcZPohhtuUCQSkeM40XNVx5FIRDt27EhIgAAAAAAAAIg/4ybRRx99lMg4AAAAAAAAUI+sv0lUUVGh119/XevXr1fDhg3Vs2dPZWVlJSI2AAAAAAAAJIh1k+ipp57Srl271L9/fzmOo4ULF6qgoEDZ2dmJiE+S9GFhP7VvtjTw2C9Plaq03nTtmy09qSxbDO6y3PliKcuUJihvVezuGNx5/Npmis+vb2yxutvvjsWvbHc6b33uPgrbX95+MN17b32mtnrPB/WDqa1BeYJ48/j1QdDY9Y6zoGO/eeF3Psx9MJXnvZ9VaU1jwVaHt0y//va7z+54vOX5xWw6DpvHlCbs/Ky6ZhqHpnHtN3+C5nTVefd/g9rovea+D964gtYxb/x+ZbrTuvP4ra0m3nEQtCb7PUNMzxa/+eaXJ+xaYLr/Ya4HrbPumIPSmeaQtw229vj1pa3dNVknwwpTvqlNsawntW1DmGe797NpzfXmNT3z/Ma3aczYYg3q51jWozDl2co3rVNBbXXXG0udpnymNtrWG7+ywzLNY2/53rhM9Zri8zsXtI74xej3DDCVb1qDTcfedwu//vEKmt+mMvzuqbdMU1tM4z/M2mx6tsQy/muy7vrVa6o76Floi8PvPgWNUe87jKkPvPz6N+g+BrXNXaY3X9i6vddqy6/PbM98W7yx1h/03hnmurcdQe9+fvfFNGeqzgXNizDvcba1wu+8je1++MVkq6MqTdDYdaer+vzxwVGSdgfGa90k2rhxo9566y0lJydLkvr166d+/foldJMIAAAAAAAA8dXAlqBJkyaqqKiIHkciEZ1zzjlxDQoAAAAAAACJZf0m0RVXXKEhQ4YoKytLDRs21PLly3Xeeedpzpw5kqQRI0bEPUgAAAAAAADEl3WTqKSkRG3bttW2bdskSRdffLEk6ZNPPolvZAAAAAAAAEgY6ybR5MmTExEHAAAAAAAA6pF1k+jdd9/VrFmzdOTIkWrnFy5cGLegAAAAAAAAkFjWTaKcnBwNGzZM3/72t2MuvKioSHfddZdmzpwZ/TM1AAAAAAAAnHqsm0Tnn3++hg8fHnPBH3zwgXJycrR79+6axAUAAAAAAIAEamBL0LNnT7322mvas2eP9u3bF/0/mwULFujxxx9XWlpanQQKAAAAAACA+LF+k+jQoUN67rnn1KhRo+i5SCSiv/71r4H5nnnmmdpHBwAAAAAAgISwbhKtWbNG69evV9OmTRMRDwAAAAAAAOqB9c/Nzj//fDVp0iQRsQAAAAAAAKCeWL9JdPnll2vIkCHq0aOHUlJSoudHjBgR18AAAAAAAACQONZNohMnTuiSSy7hXykDAAAAAAD4F2bdJJo8eXKtKli9enWt8gMAAAAAACD+rJtEW7Zs0axZs3T8+HE5jqPKykp99tlnWrt2bQLCAwAAAAAAQCJYf7g6JydHnTp1UlFRkTIzM9W4cWOlp6cnIjYAAAAAAAAkiPWbRJFIRPfdd58OHTqkSy+9VJmZmbrjjjsSERsAAAAAAAASxPpNorPOOkuS9O1vf1s7d+7UmWeeqQYNrNkAAAAAAADwDWL9JtHVV1+tcePGaezYsfqP//gP7d69W0lJ1mwAAAAAAAD4BrF+JejRRx/Vvffeq0suuUSPPfaYKisr9eyzzyYiNgAAAAAAACRI4FeCHMdRRUWFOnbsqKKiIpWUlGjIkCFq3bp1QoL7sLCf2jdbetJn97mq46rr3nNupnPu+rznvOmq6vGWVXXOL25TfX7X3XlN57z1BsXtjdmUzt1/7v8GleNtl18ev3Tez6Z8frGYjsP2hV8cfnlsbfXGHabtiD+/OZ7ouoPGlmnNCEofdh4Gtd1vjXTXY+K3rnqv+cXrri+oDtPc9lvb/WKylW8S5lniThtmzTHVE3ZtsJVlum+m8eGXvjZ95s7vjcnv2ZKINdH2PDU9l93Hfvfam9ddtu056/3sN36C2uAt78PCftb3grDPJNu8MsUTVI87jen9yNYnptiC8vh9DvsOYpojpvxBZZrmYtD8q4kw98s0xr3HNY3H1Gfe8oOeG+68ftfCzo+q+vzShPnfDkHvd2H7xa+coPd/U36/9tRmzQ6q3/Z8tfWJaS3zm/Om++I3NoLKNbXPXW6YssKsE37P+arjoHHmjcUUe02Y3iOD5rZpffaLza9c0xrvxxSXX/qa/u8m0zMjaJ0Juh60bgfNAW8sYeeaqU6/fKZzsYwn2/tKWA1MF/7+97/rlltu0bp163TixAn9+7//u/Ly8jRs2DBt2LChVpUCAAAAAADg1GLcJJo2bZrGjRunHj166K233pIkvfXWW1qwYIFeeOGFhAUIAAAAAACA+DNuEu3fv1/9+v3z60rvvvuubrnlFjVo0EAtWrRQUVFRwgIEAAAAAABA/Bk3idz/zP2WLVt03XXXRY9LSkriGxUAAAAAAAASyvjD1d/61rf00UcfqaioSIWFhdFNor/+9a9q3rx5wgIEAAAAAABA/Bk3iR588EHde++9Kioq0kMPPaTU1FTNnj1bM2fO1C9/+ctExggAAAAAAIA4M24SdezYUX/605904sQJnXPOOZKkTp066Y033lDr1q1DFZ6fn69f/epXKi8v1z333KOhQ4fWSdAAAAAAAACoW8ZNIklKSUlRSkpK9Pi73/1u6IIPHDigvLw8LVq0SCkpKbrrrrvUuXNnXXbZZTWPFgAAAAAAAHFh/OHq2tq4caNuuOEGnXvuuUpNTVXv3r319ttvx6s6AAAAAAAA1ELcNom++OILNWvWLHqclpamAwcOxKs6AAAAAAAA1ELcNokqKysViUSix47jVDsGAAAAAADAqSNum0QXXHCBCgsLo8eFhYVKS0uLV3UAAAAAAACohbhtEnXt2lWbNm3SwYMHVVxcrJUrV+qmm26KV3UAAAAAAACohcB/3aw2mjdvruzsbA0fPlxlZWUaOHCgrr766nhVBwAAAAAAgFqI2yaRJGVmZiozMzOeVQAAAAAAAKAOxO3PzQAAAAAAAPDNwSYRAAAAAAAA2CQCAAAAAAAAm0QAAAAAAABQnH+4ujYu+taZSm7aSq3Oa6Rkp5UkqdV5jaqlcZ9Pdk5O65fPVJbtmjtN1fWq9H51+KUJaoO3He5y3OdMeYLi9uYPSmOr0y9+b7x+ZZrab2tn0Hm/2Lxxhm27N7ZY+iwova0cv+t+/eU31k1l2OZC0Jjzy29K59cH3piD+sc2pvzaYWqfqdwwczFofJr6MsxcC7qPfnV5Yw9Tnl+cQXMjaD4G1RV2joSZE7Z4bf1tyhtUj43t/obp17Ax+N1jd1nua6YYbW2J9Z6FqT9MW2xpw4yXWOZs0NrpvR52/IdZV/zOhXkWxiLM+hfr/LPNK792B8UWdN72LA67VgWtUWH7xC1oLYzlPoWZK+6y/drsd87WxqD6w74beI/DzMMwcdie3bUV5t2gStj2BNVjSxtmbQ/zTmHKb6ojzHrrLs8tlvcZU9y2MR+mv8L0VVhhywq7dnrLdZcdNk9dPcODxnmsz5gw7am6FmYOhInL1Da/cmvyLmV7Zod5rwobQ5jnp1vQ+4Ypn+m5ZIvfG5upTne6JOdi33a6RRzHcaypAAAAAAAA8C+NPzcDAAAAAAAAm0QAAAAAAABgkwgAAAAAAABikwgAAAAAAABikwgAAAAAAABikwgAAAAAAABikwgAAAAAAABikwgAAAAAAABikwgAAAAAAAA6xTaJ8vPz1bdvX6Wnp+u1116r73BQz2bMmKGMjAxlZGRo2rRpkqSNGzcqMzNT6enpysvLi6bdsWOHsrKy1Lt3bz322GMqLy+vr7BRT6ZOnaoJEyZIYpzgZKtXr1ZWVpb69Omjp59+WhLjBP6WLFkSffZMnTpVEmMF/1RUVKTbbrtNn332maTYx8W+ffs0dOhQ3XrrrfrhD3+oY8eO1Us7EH/esTJ//nzddtttyszM1COPPKLS0lJJjJXTnXecVJk7d66GDRsWPWacnN6842TLli0aNGiQMjIy9OCDD8ZnPXFOEf/7v//r9OjRwzl06JBz7NgxJzMz09m5c2d9h4V6smHDBufOO+90SkpKnNLSUmf48OFOfn6+0717d2fPnj1OWVmZM3LkSGft2rWO4zhORkaGs2XLFsdxHOeRRx5xXnvttXqMHom2ceNGp3Pnzs7DDz/sFBcXM05QzZ49e5wbb7zR2b9/v1NaWuoMHjzYWbt2LeMEJzl+/Lhz3XXXOV999ZVTVlbmDBw40PnjH//IWIHzt7/9zbntttucdu3aOXv37q3Rs+a+++5zli1b5jiO48yYMcOZNm1avbQF8eUdK59++qnTq1cv5+jRo05lZaUzfvx4Z86cOY7jMFZOZ95xUmXnzp3O9773Pefuu++OnmOcnL684+To0aNOt27dnB07djiO4zjZ2dnR8VCX4+SU+SbRxo0bdcMNN+jcc89Vamqqevfurbfffru+w0I9adasmSZMmKCUlBQlJyfrO9/5jnbv3q1WrVqpZcuWSkpKUmZmpt5++219/vnnOnHihDp27ChJysrKYuycRg4fPqy8vDyNHj1akrR161bGCapZtWqV+vbtqwsuuEDJycnKy8tTo0aNGCc4SUVFhSorK1VcXKzy8nKVl5ercePGjBVowYIFevzxx5WWliYp9mdNWVmZ/vKXv6h3797VzuNfj3espKSk6PHHH1fjxo0ViUR0+eWXa9++fYyV05x3nEhSaWmpcnNzNWbMmOg5xsnpzTtONmzYoI4dO+qKK66QJOXk5KhXr151Pk6S4tOc2H3xxRdq1qxZ9DgtLU1bt26tx4hQn9q0aRP9vHv3bq1YsUJ33333SWPkwIEDJ42dZs2a6cCBAwmNF/UnNzdX2dnZ2r9/vyT/tYRxcnorKChQcnKyRo8erf379+vmm29WmzZtGCc4SePGjTV27Fj16dNHjRo10nXXXceaAknSM888U+041nFx6NAhNW7cWElJSdXO41+Pd6xcdNFFuuiiiyRJBw8e1GuvvabJkyczVk5z3nEiSc8++6zuuOMOXXzxxdFzjJPTm3ecFBQUKDU1VdnZ2fr000/13e9+VxMmTND27dvrdJycMt8kqqysVCQSiR47jlPtGKennTt3auTIkRo/frxatmzpO0YYO6evN954Qy1atFCXLl2i50zjgXFy+qqoqNCmTZv0s5/9TPPnz9fWrVu1d+9exglO8tFHH+n3v/+91qxZo3Xr1qlBgwbavXs3YwUnifVZ4zc+GC+nlwMHDuiee+7RHXfcoc6dOzNWUM2GDRu0f/9+3XHHHdXOM07gVlFRofXr1+vBBx/UokWLVFxcrFmzZtX5ODllvkl0wQUX6P33348eFxYWVvv6HU4/mzdv1pgxY/Too48qIyND7733ngoLC6PXq8bIBRdcUO38l19+ydg5TSxfvlyFhYXq37+/jhw5ouPHj+vzzz9Xw4YNo2kYJ2jatKm6dOmiJk2aSJK+//3v6+2332ac4CTr169Xly5ddP7550v659eyZ8+ezVjBSbz33zYumjRpoqNHj6qiokINGzbkPfc0s2vXLo0aNUrDhg3TyJEjJZ08hhgrp7dly5Zp586d6t+/v44fP64vv/xS48aN009/+lPGCaKaNm2qDh06qGXLlpKkPn36aO7cucrKyqrTcXLKfJOoa9eu2rRpkw4ePKji4mKtXLlSN910U32HhXqyf/9+/ehHP9L06dOVkZEhSerQoYP+8Y9/qKCgQBUVFVq2bJluuukmXXTRRTrjjDO0efNmSf/8l2kYO6eHOXPmaNmyZVqyZInGjBmjnj176qWXXmKcoJoePXpo/fr1+vrrr1VRUaF169bp1ltvZZzgJFdccYU2btyo48ePy3EcrV69mmcPfMU6LpKTk3Xttddq+fLlkqQ333yT8XKaKCoq0g9+8AONHTs2ukEkibGCaiZPnqwVK1ZoyZIlevrpp9W+fXv9/Oc/Z5ygmhtvvFHbtm2L/szGmjVr1K5duzofJ6fMN4maN2+u7OxsDR8+XGVlZRo4cKCuvvrq+g4L9WT27NkqKSnRlClToufuuusuTZkyRQ888IBKSkrUvXt33XrrrZKk6dOnKycnR0VFRWrXrp2GDx9eX6Gjnp1xxhmME1TToUMHjRo1SkOGDFFZWZm6deumwYMH69JLL2WcoJobb7xR27dvV1ZWlpKTk/Vv//ZveuCBB9StWzfGCqqpybPm8ccf14QJE/SrX/1KLVq00HPPPVefTUCCLFy4UF9++aXmzJmjOXPmSJJ69uypsWPHMlYQCuMEVVq0aKEnn3xSo0ePVklJia688ko9/PDDkup2nEQcx3Hi2hIAAAAAAACc8k6ZPzcDAAAAAABA/WGTCAAAAAAAAGwSAQAAAAAAgE0iAAAAAAAAiE0iAAAAAAAAiE0iAABwmhg2bJhmzZp10vmXX35ZP/zhD0OVMW/ePN8ywtq6datyc3NrnB8AACCe2CQCAACnhSFDhuj3v//9SecXLFigoUOHhipj8ODBuu+++2ocw9///ncdOHCgxvkBAADiKeI4jlPfQQAAAMRbeXm5evTooby8PF177bWSpPfee08TJ07UgAEDtHr1ap04cULFxcV6+OGH1atXL73wwgv629/+pi+++EJt27ZVq1atdOjQIeXm5mrNmjV68cUXVVpaqoMHD2rAgAEaN26c3n33XeXl5ally5bauXOnysvL9cQTT+jCCy/U4MGDdfToUaWnpysnJ0ePPPKICgoK1KBBA7Vr105PPvmkGjTg/4cHAADqB28hAADgtJCUlKRBgwZp4cKF0XPz589XZmamNm3apFdffVX5+fnKzs7W888/H03z+eefa/HixZo+fXr0nOM4evnllzVlyhQtWrRI8+fP16xZs3Tw4EFJ//yzspEjR+rNN99UVlaW8vLy1KJFC40ZM0bXXnutJk+erFWrVunYsWNasmRJNKa9e/cmqDcAAABOllTfAQAAACTKoEGDlJGRoaKiIpWXl2v9+vWaNGmSBg4cqPz8fBUUFOiDDz7QsWPHonk6duyopKTqr0yRSEQzZ87U2rVrtWzZMu3atUuO46i4uFiSdOGFF+rKK6+UJF111VVavHjxSbFcc801ysvL07Bhw9S1a1fdc889atWqVRxbDwAAEIxvEgEAgNNG8+bN1bVrVy1fvlxvvvmmevfurT179ujOO+9UUVGRunXrplGjRlXLk5qaelI5x48f1+23365t27bpqquu0vjx45WUlKSqv+I/88wzo2kjkYj8/rq/ZcuWWrVqle677z4VFRVpxIgRWr16dR23GAAAIDy+SQQAAE4rQ4cO1YwZM3T06FFNnTpVmzZtUvv27TVixAhVVFToiSeeUEVFRWAZBQUFKioq0rhx45SSkqI333xTpaWlqqysDMzXsGFDlZeXS5J+97vfafPmzZo+fbq+973v6auvvtL27dvVs2fPOmsrAABALPgmEQAAOK107txZhw8fVuPGjdW2bVvddtttOnTokPr06aO+ffsqNTVVR44cUVFRkbGMtm3b6uabb1afPn3Up08frVmzRpdddpkKCgoC6+7YsaP27t2rH//4xxowYIAqKirUt29fZWVl6ejRoxo2bFhdNxcAACA0/nUzAAAAAAAA8E0iAAAAAAAAsEkEAAAAAAAAsUkEAAAAAAAAsUkEAAAAAAAAsUkEAAAAAAAAsUkEAAAAAAAAsUkEAAAAAAAAsUkEAAAAAAAASf8HDzIXW0ingJIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x172.8 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#SOL:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. 3 Pairwise measures of LD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen:</b>\n",
    "Since we generally cannot directly observe the tree sequence underlying a real samples of genomes, measures of LD are generally based on mutations. In the simplest case, we can measure the correlation between a pair of loci using the correlation coefficient:\n",
    "\n",
    "$r^2$ <a href= \"https://en.wikipedia.org/wiki/Linkage_disequilibrium\">(Hill \\& Robertson 1968)</a>\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 15:</b> Calculate $r^2$ between the two pairs of mutations you've considered in exercise 12, i.e. 13 and  19 in the <i>D. melanogaster</i> ts?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dmel_array[13])\n",
    "print(dmel_array[19])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We can use <code>LdCalculator()</code> to compute pairwise measures of LD such as $r^2$:\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.24999999999999992"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ld_calc = tskit.LdCalculator(dmel_ts_muts)\n",
    "ld_calc.r2(0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 16:</b> Does <code>r2()</code> give the same answer as your calculation above?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL ld_calc needs two arguments, specifying the pair of mutations\n",
    "#Check whether computing D^2 (and from this r^2 by hand) gives the same result...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We can use <code>LdCalculator()</code> to efficiently compute the entire matrix of pairwise $r^2$ values:\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmel_ld_calc = tskit.LdCalculator(dmel_ts_muts)\n",
    "dmel_A = ld_calc.r2_matrix()\n",
    "\n",
    "# Now plot this matrix (don't worry about the plotting code, it is simply there to visualize the matrix):\n",
    "x = dmel_A.shape[0] / plt.rcParams[\"figure.dpi\"]\n",
    "x = max(x, plt.rcParams[\"figure.figsize\"][0])\n",
    "fig, ax = plt.subplots(figsize=(x, x))\n",
    "fig.tight_layout(pad=0)\n",
    "im = ax.imshow(dmel_A, interpolation=\"none\", vmin=0, vmax=1, cmap=\"Blues\")\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for s in \"top\", \"bottom\", \"left\", \"right\":\n",
    "    ax.spines[s].set_visible(False)\n",
    "plt.gcf().colorbar(im, shrink=0.5, pad=0)\n",
    "plt.savefig(\"ld.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 17:</b> Use this code to plot the matrix of $r^2$ values for the <i>Homo</i> simulation. How does this differ?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL Recomputing the matrix of r2 for the homo_ts_muts\n",
    "homo_ld_calc = \n",
    "homo_A = \n",
    "\n",
    "# I am using the same plotting code as above\n",
    "x = homo_A.shape[0] / plt.rcParams[\"figure.dpi\"]\n",
    "x = max(x, plt.rcParams[\"figure.figsize\"][0])\n",
    "fig, ax = plt.subplots(figsize=(x, x))\n",
    "fig.tight_layout(pad=0)\n",
    "im = ax.imshow(homo_A, interpolation=\"none\", vmin=0, vmax=1, cmap=\"Blues\")\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for s in \"top\", \"bottom\", \"left\", \"right\":\n",
    "    ax.spines[s].set_visible(False)\n",
    "plt.gcf().colorbar(im, shrink=0.5, pad=0)\n",
    "plt.savefig(\"ld.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen: Visualizing the decay of LD</b> \n",
    "\n",
    "When analyzing WGS data, we are often interested in how quickly LD decays with physical distance. We can think about this decay both in terms of the correlations between trees and pairwise summaries.\n",
    "\n",
    "Note that arguments about LD decay are often used to justify treating SNPs (or loci) that are more than some threshold distance apart as unlinked. While this is standard practice in the field (for a lack of better alternative and may be fine for a particular analysis), it is **WRONG!**\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code chunk below bins the $r^2$ values in the matrix <code>dmel_A</code> we've defined above by the physical distance between loci:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the list of SNP positions from the ts sites table\n",
    "pos = dmel_ts_muts.tables.sites.position\n",
    "# Sets up a matrix of distances for all pairs of SNPs\n",
    "pos_i, pos_j = np.meshgrid(pos, pos)\n",
    "distance_ij = abs(pos_i-pos_j)\n",
    "# Defines a new array which has distance and r^2 as separate columns\n",
    "dis_r2_arr = np.column_stack((distance_ij.flatten(), dmel_A.flatten()))\n",
    "# Sorts based on the distance column and remove the first len(pos) entries. These are the diagnoal elements\n",
    "# which by definition have distance 0 and r^2=1\n",
    "dis_r2_arr_srtd = dis_r2_arr[np.argsort(dis_r2_arr[:, 0])][len(pos):]\n",
    "# Defining bins of 5 bases\n",
    "dis_bins = np.arange(0,1_000,5)\n",
    "# Binning by distance\n",
    "dis_binned = dis_bins.searchsorted(dis_r2_arr_srtd[:,0])\n",
    "bin_means_r2 = (np.bincount(dis_binned, weights=dis_r2_arr_srtd[:,1]) /\n",
    "             np.bincount(dis_binned, minlength=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 18:</b> Plot the mean LD (as measured by $r^2$) for each distance bin (<code>bin_means_r2</code>) against distance <code>dis_binned</code>!\n",
    "\n",
    "How does LD decay in <i>D. melanogaster</i> and differ from that in <i>Homo</i>? What is the half-distance of LD? What is the background level of LD?\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#SOL\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "sns.lineplot(dis_bins, bin_means_r2[:-1], label='The decay of LD')\n",
    "plt.xlabel('bases')\n",
    "plt.ylabel('r^2')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}