Varscan Muliallelic check.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temperature experiment allele counts by filtering method\n",
    "---\n",
    "We are going to first count the number of multiallelic sites within unique plants, then within groups of plants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Dan Koboldt's 2012 comment on the `readcounts` utility](https://sourceforge.net/p/varscan/discussion/1073559/thread/c647b0ba/#2acb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load necessary package(s) and temperature experiment data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Erik\\anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py:3063: DtypeWarning: Columns (9,58) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "data = pd.read_csv(\"variant-calls_temperature-expt_2020-04-07_tables.csv\",\n",
    "                   header=0,\n",
    "                   index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>Note:</b> The above warning stems from incomplete data in rows from the 32$^{o}$ data, which we will remove.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove 32 degree data\n",
    "data = data.loc[data.temperature != 32]\n",
    "\n",
    "#fix the data types for the two problemmatic columns\n",
    "data.directVegParent = data.directVegParent.astype(str)\n",
    "data.DP = data.DP.astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to return a dictionary of counts of sites with $X$ alleles "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_allele_counts(data):\n",
    "    #initialize variables to hold working data\n",
    "    running_allele_counts = []\n",
    "    count_data = {}\n",
    "    \n",
    "    #count alleles per site by unique plant ID\n",
    "    for plant in data.plantID.unique():\n",
    "        plant_data =  data.loc[data.plantID == plant]\n",
    "        for position in plant_data.pos.unique():\n",
    "            running_allele_counts.append(plant_data.loc[plant_data.pos == position].shape[0])\n",
    "    \n",
    "    #count the number of sites with i _alternate_ alleles;\n",
    "    #add 1 to dictionary index to represent _total_ alleles per site\n",
    "    for i in range(1,max(running_allele_counts)+1):        \n",
    "        count_data[i+1] = running_allele_counts.count(i)\n",
    "        \n",
    "    #return the dictionary    \n",
    "    return count_data\n",
    "\n",
    "def print_counts(data):\n",
    "    \n",
    "    #print the allele counts nicely\n",
    "    print(\"Alleles/site\\t# Observed\")\n",
    "    for count in data.keys():\n",
    "        print(\"{}\\t\\t{}\".format(count,data[count]))\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_allele_counts = get_allele_counts(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alleles/site\t# Observed\n",
      "2\t\t14483\n",
      "3\t\t8799\n",
      "4\t\t5083\n",
      "5\t\t2336\n",
      "6\t\t750\n",
      "7\t\t165\n",
      "8\t\t17\n",
      "9\t\t1\n"
     ]
    }
   ],
   "source": [
    "print_counts(raw_allele_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove the alleles with no representation in one technical replicate: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_na_data = data.dropna(subset=['DPReplicateA','DPReplicateB'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_na_counts = get_allele_counts(no_na_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alleles/site\t# Observed\n",
      "2\t\t15008\n",
      "3\t\t5056\n",
      "4\t\t1258\n",
      "5\t\t120\n"
     ]
    }
   ],
   "source": [
    "print_counts(no_na_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter alleles by minimum 3% in at least one technical replicate, and present in the other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alleles/site\t# Observed\n",
      "2\t\t1158\n",
      "3\t\t12\n"
     ]
    }
   ],
   "source": [
    "min_3percent = pd.merge(no_na_data.loc[no_na_data.freqPropReplicateA >= 0.03],\n",
    "                        no_na_data.loc[no_na_data.freqPropReplicateB >= 0.03],\n",
    "                        how='outer')\n",
    "\n",
    "min_3percent.drop_duplicates(inplace=True)\n",
    "min_3percent_allele_counts = get_allele_counts(min_3percent)\n",
    "print_counts(min_3percent_allele_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at how minimum frequency cutoff affects max allele count per site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_snp_freq_data(data, minimum=0.01, freq1='freqPropReplicateA', freq2='freqPropReplicateB'):\n",
    "    \n",
    "    temp = pd.merge(data.loc[data[freq1] >= minimum],data.loc[data[freq2] >= minimum],how='outer')\n",
    "    return temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [i/1000 for i in range(0,31)]\n",
    "y = [max(get_allele_counts(min_snp_freq_data(no_na_data,j))) for j in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Alleles represented in both tech reps')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEnCAYAAAAq8Q2oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZhcVZ3/8fcnnaSbhHQIJOkAIQRFkU1BW1xQRERFwegoM6Ci4DKIK+7+cAHEFUVFhxkB2QVlUzSDCoIQlEGWIPu+Q4AsEAgJkP37++OcSi6V6u6q6upUdffn9Tz1dNe959577lL1rXvuWRQRmJmZDSUjmp0BMzOzRnNwMzOzIcfBzczMhhwHNzMzG3Ic3MzMbMhxcDMzsyGnz+Am6SBJkV8vrTB/98L8PQcik4Vt7D4Q629FlfZZ0ixJs+pYV+kcTm9cDluTpCMl7dHsfAwFkqbn6+agOpYdNtdcLSTtlK/RjfuxjtdJukbSs/kY75Snf13Sw5JWSrqxcbnuNS9HSmrJ9mS13LktBj5UYfqH87yB9C/gdfmvWW+OABzcGuNx0ufuT3Us+6e87OMNzdHgtxPpGq07uAEnAyOBd5GO8d2SdgG+B5wN7Ebl7+phpZbg9nvgAEkqTZC0AfA+4HeNzlhRRDwTEVdHxDMDuR2zRlEyutn56I+IWJY/dwvqWHZBXnbZQORtuJI0AtgG+FNEXJaP8XPAtjnJ8RFxVUTc0rxctoiI6PUFHAQE8BZgNfDGwrwPAEuAfXKaPQvzXg2cD8wBngfuAr4PbFBIs0Oe97OybX4fWAbsnN/vnte/eyHNLOBKYC/gxryeG4DXkH7VfJ/0q3EhcBowtrDsOusr29fphWkPAmeSfgndlbfzD+AlwFjgBOBJYB7wE2BkFcf026S70EXAE8BlwGvL0vS0z7PK0k0Efgk8mo/ZncDBfe1Xnv6fwE3A0pyPk4GN+8p/Ydl/5ePxFHAF8PrC/E2BM/J6lwE3AweUrePIdAmus+7TgAcL76fn/H8COCqf16eB/wWmFtJFhdeRvezD1sCvgQfyftyfj+WECmnfBFySz9mz+bh9rMJ18tF8DlYA/1bDsZgCnA48ltM8DlwITM7zRwLfAe4rnK8rgTf0cZ5mUf/npHTcDyo7N3OAnUmfg+eAe4BDBvKzVGl9PV1DOd13gS8BD+Xz9Sdgcn6dm8/jI8DXqrzeJwH/k5dZlv/+GmivdM1W+swW9qH8NT3P7wSOK1wDdwFfANTL8g/mbVR13QN/Bq6vMH1TYCXw+cL+ngDcnc/xI8BvgM17O/6Vrpk+vnPfC1ydt/E0cB4wrSzNB0jX7JJ83m4BPtHXORtJ9R4C/k66MP+Rp30YuCBvtNw00ofpNFKx5fbA4cCLgP1JR+RWSV8CjpP014j4i6Q3A18DvhIRN/SRp62BH5Nux5cAPwJm5tdI0sWwbU4zH/hqDftbtBvw4pyv0cCxpLvV+4F78/7sBnyT9OXzP32sb3PgZ6QvibHAAcDfJXVHxM3VZkpSJ/B/wAaki+wB4O3ALyW1R8R/9bLsD0kf/l8AX8l5+i6wg6TXR8SqXpY9Ji97MqmIZTXwWtI5v0rSWFKwmwB8nfTBOAD4taQxEXFitftY5jDgKlIAmUz6AjyLFHggFdH8k3TNnZCnzellfZvl+Z8nBegX5fz+Oa+rtL/vJp3v/yMF2CdI1/OWZet7M6nY6duk6+3BGo7Fr/P6vpLTdJF+UI7J879G+qL7Bulz1Ql0U13xVqM/J52kL7pjST82PkK65u6KiMv7WLbRn6WefAi4FfgU6VgeS/qBMQ74C3Ai8O/ADyXdEhF/7mlFkiaQrruNSZ+Rm0nX37vzPlR7d/qnvPw387ZL1+bj+Y7sT8ArSd+TtwB7Az8lBZqv5/lvIP1YORk4KW97GemaOowULB6n5+v+DOC3kraLiNsL0z+Q//42/92Y9CPqMGAB6bPyJeD/JL0sIpZWuc89knQI6cfkqaTraBzpe+wKSS+PiMWS3kD6QVT6nhoBvAzYqM8NVPGL5SBSxN2a9KXyFNDB2kj/VtZG5T17WIdIH6IDSF+Em5TN/wPp19r2pDuQi8i/VnqK+qRfKyuAFxWmzcjpLi1b/++BB6r4FVHa1/JfmwuB8YVpn8vpTipb/l/A5dX8Eiws05aPzV3Az6vY51mF998iXYAvKVvnr0hfwCMr7Rfp19Uq4PCy5XbN6d7TS363zsv+tJc0n+nh+F5K+vJsq/Srr5DuNCrfuV1Rlu7LefpmhWkBfLeWc1BYdiTpyyNYW2qgfA3MBkb0suyDpF+fU+o8FkuAz/Wy/guB39exT7Oo/3NSOu4HlZ2bAN5cmNaer7cTB+qzVGl9PV1DOd3dvPDO76d5+jfLzvd84NQ+juFR+ZrfuZc0L7hmy47/rAr7sXVZulLp10Fl00sBbGIhz+vcmQEfr3R8KuRnA9Ldzw/Kpt8I/LmX5dqALfI2/q2n41/pmsnTd6fwOQA2zPk4pSzddGA5a+8gvwwsrPW6j4iamwKcR7qQ3wV8EJgL/K1SQkmdko6WdF8+OStIv05FKoYo+liefz3p5B0Yec/6cHdE3F94f2f+e3FZujuBqcXnhTX6Z0QsqnI7W/S1Mkl7Srpc0pOkHwgrgJeSytJrsRdwDfCApJGlV87XJsB2PSz3VtIvoLPKlrsGeIb0y7kne+Zle7v72g14NCJmlU0/k/QrtKd89aW8YkPpucK0elYmaXSuYXanpOdJ56FUKrFN4e+WpC/f1X2s8uqImFs2rdpjcR3wFUmHStqxwrV6HfBOSd+T9IYan+c1+nPyXBTu0CI9V7uH6s5DQz9LvbgkIlb2tp08/94qtvM24LrouySpP3Yj/fD/bdn0M0l3h69bZ4k6RMTzpDvlD5bOs6QdgVeQ7urWkPRJSTdJWkL6nno4z6r1e6qS15FKAMq/g+aQzlXpO+g6YIKkMyXtI6nvO7aspuAWEYtJd1kfIhVJntXLB/5U4BDS7eRbSc/gPp3ndZSt90nSF1c78NuImFdllp4qe7+8l+kjSb8+6lHLdjrohaRXkoq9lpCC+mtJx+amvpatYDLpIlhR9jovz9+kl+UgfbDLl+3sZbniOnsr7tuYyrXk5hbm12Nh2ftScVCtx63kB6RfnmeSioB2IRXrFNdZzf6WVNrnao/FfqRiwq+Sir0elXR4Lq6C9GzsCNJd1z+AJyWdKmliFflq9OekfDlI56Ka89Cwz9J63M4mVHf++2Nj0h1KeRFnfz8zlZxBCui75/cfIj06+mMpgaTPkoqELyV9JnYhfVdB/85LSek76FLW/Q7akfy5i4grSEW4W5AegS2QdKmkl/e1gVqeuZWcQQpEI4D3V0ogqYNUHn1kRPy8MH3HHtLvSaqgMBv4lKQzI2J2HXmrVqm8uPzXb29f6o3yPtKvoPdGxIrSxFyu/3SN63qSVKxyaA/z7+plOUi/SCt9UT1ZYVrJE/nv5r2sfyGVf91NKVv/Ukh3UBGxvJBufZwHSM93zoiI75YmSNqwLE1xf/tSqbShqmMREfNJP/4+LWkb4EDSs7sFwC/ztXI0cLSkKaRirJ+SnsntV0XeBrtmfmafoO/zv5R18wYpf719nkoWAhtX+CyUf2Ya4QrSXdgBkq4gfY+fn+/qSvYH/hYRXypNkLRVFeuu9jyV9ucg4LYK61nTvCwizgfOz5/N3Umfg4skTe2tNKWeHkouIdU2Oj4iKmUK0h1YGykKFx1UnjD/8jyDdDfzelKtmN9U+JJppIfy3x3Kpr9zALdZMoZUfr/mizA3Oq6naO0i0sPVhyNidoVXT+0PLyEVgUzrYbkHetnmpXnZg3tJcwWpeGvXsukfIAXjO/L7dc5DLnZ4fS/r7sty0nOFaoxh3Wv0I2Xv7yY9K/p4ncXa1R6LNSLiroj4OumHR/k1SkTMjYiTSOdinflDVKVrZSTpB9pA+yuwi6RX9JLmIaCreCct6cWs+8OmdGdWfo1eQfo+/vey6R8kXdNX15rpnuRHPmcB+5K+86ZSViRJdZ+NSuaR9rH8uty77P1VpAC2dQ/fQev8cI6IJRFxIamy2Kb08cOm5ju3SLXoKt6xFdIsknQ18CVJj5N++XyUyr9+TiE9h/tIRKyQVKr2+V9UdzBrFhGP518sh0l6gvQlcwCpFtdAu4hUO+80SaeSnrV9i1SRplY/I/1q/4ekn5HupMaSAt4bI+LdlRaKiPskHU2qpboN6YO1lHTr/1bS86WKtd7ysj8DvihpHKkobRWp2OLOiDiH9HD9UOD3kr5BKtL5YF73J2JtTcy/kB4q/0rSEaQfRV+lcu3bat0O7C3pIlJweCwiHush7UXAgZJuIRXRvpeywBoRIenzpMoWl0k6nnQ3tS2pmv4RfeTnNPo4FpLGkwLVWaxtRvBuUg3LvwJI+iOp6Ppfeb92Jj1zPYHh4TpS7ckf56LaZaSakO3rYds/I/0YuVTSd0nPeieSztEh+UfkeaSmGmdJ+mmefxhr7/xLSjUUPy3pdNK5vpn0WbgSOF7SJNLdzDtJFUV+EBHl6+mvM3L+jifVzr2ibP5FwNckfR24ltQxwr59rTR/Xs4BPibpbtJ30t6sLQItpXtG0leA/877W/ou2JxU+3lWRPxG0lGk2q6Xk5pITCVVQrox+mp/2VeNE3qo3dNDTZhiO7fpOcOLScHjuLyTxRoznyHdBby1bH0H5HT7VappE2trIV1ZoaZNAB8vm35knl6sPTWV1E7qaVK59vepUOOI3Danr/2NtTWm5lRxTD/L2rZV15EqaczihbWqetrnWWXrmkD68D1A+oU3n/RM5vMVzuH0smU/RPpF+CwpoNyRz9PUKvbhENKHchmpSGUW8LrC/E1JFYh6bNuV070hH4PnSHdJB9Bzbcny81rpGO1Kqpi0lL7buU0k9ejwVH6dRXr+Wam21x6kD9iS/LqJ9IOsx+uk2mNB+oI+gfSFtoRUqec64AOFNF/K5+pJ1rYbPRIY1cd5mkWdnxN6aefWw3aK1+8611ylY0QNnyVSbepZ+Rg9DHyRXtq5lU0r5ae8luI6x6eH4ziZVInqcdLn7BFSu8T2Qpr3kJofPJ+vj7eVH5ec7gjSj9lSCc70PL3Uzq20jbsptHPLafpVW7JsmevyMt+vMG8DUjX9BaTv8AuBrcq33cPx34i11/tCUgB9wXd/Ie07SZ+rZ/Jxu5d0w7Ndnr83qRLQ46xtX3gyhRrSPb1KjQPNzMyGDI8KMExIOi0XqZQ6Za6q9pdSx6hnDmzuDNZ0NnxlD/OmSVoiqeYav7Wc7yrXt+ZaWp98LVotHNyGGKWRA56StD6eRVhBPvYfH4h1R8TDEbFh9NJzzEDoLeBWuXxpZIF6amab1c3BbQhRGl7kjaSy7RlNzUwNGv3F5y9SK+drYvhxcBtaPkyqdHAaqZ1UVSRtJul3khZIekDS53pJ+1pJV0l6OvdesHth3kGS7pe0OK/ngz2s40hJ5+deB54BDpI0QtL/k3SfpCclnas85lXh1//Bkh6T9LhSn6T1rq8jp30y78d1krryvPGSTs7beFTSd0tFgaW7GEnH5LvjByS9I8/7HumHxXG5+PC4PP1lki6RtFDSXZL+o5DvTSTNlPSMpGvppbZu+R1Qvkv8jqT/y8f7r+qjQbdSbyxPSHqweG7yPp+Rz/9Dkr6Zj9+2pMoAr8v7VGyHOUHSn/K2r1Gq9l7J3/Pfp/M6Xpe3+VFJd+TjeLGkLQv52b5wzOYp1dgrGZ3zuljSbZK6e9nfkPRpSfeQelDp63ycJun4PH+xpCtK+VLyM0nzJS2SdLOk4dIMY3CqtmaNX63/ItU0+hTwKlIV467CvNPINchINdTm5P9HkGoXHk5qePkiUie2b8/zjyTXcCNV032SVMNpBKk6+5OkbqTGkmo8bZPTbgps30M+j8z5e09ezwak5hFXk2qxlmoP/jann066G/1t3s6OpFpce9a5vk+QasqOIbXHfBXQmef9IacdS6ohdy25B3JSjbsVpA4H2oBPkqonlypmzaJQAzGv4xFSk5aRpE5xnygdF1JNzXNzuh1INegq1twrHIORhW3dR2pKskF+/8Melt2d1HHAT/OxeBOphmzpXJ1B6p1iXN7O3eQRD/I+l9e2PI1UC26XvF9nAWdXk+887T2ka3XbvPw3gavyvHGkmnFfIvWEMQ54TeE8LyVdf22kHmau7uXzEKQ2nRvnY9TX+TiNVDNwt3ycfl7ad1KH5NeTagIq533TZn/m/erl+7DZGfCrQScyValfwdoOVu8EvlCYfxqVg9trSI3Ai+s6jNyZLC8Mbl8Dfl2W9mLSXeJYUrOK91EY1qiHvB4J/L1s2h3AWwrvN837M7LwBfmywvwfASfXub6PkhqRvrxsmS5SdePisEzvJ3fgS/qiv7cwb0zO15T8fhYvDG77Af8o28YJpKrgpU4Oivv0fWoLbsVOgD8FXNTDsruTgltxOJtzSe0r2/I+b1eY9wleOExLpeB2UuH9O0ltHPvMd572F144XNAIUlOQLfPxvqGX6+bSwvvtgOd7uc4C2KOa81HYr7ML8zYkVdffgtQU5G5SF1Q9dqDtV+u8XCw5dBwI/DXWNvb8DdUVTW4JbJaL557ORU9fJ33RV0r772Vp30D6Bfss6cvjENIQHn+S9LJetvtIhXVfUFjvHaQvlq4elnmINAxHPev7NSkon52LOX8kaVReZlTOf2m5E1jbDx6s7euPSINEQvoSrGRL4DVlx+uDpC6VJpECbfk+1aLYSfNzveQD4Kl8jorb2ozU1m902bYfou/upmrZdrktgZ8XjslC0t3Q5qRAcl8N2+1Q78/Tise3t/OxTvqIWJLztllEXEZqg/bfwDxJJyoNOWUtyg9ZhwClEdH/A2iTVPrwtwMbSXpFRNzUy+KPkIY5KR+poae0v46I/6w0MyIuBi7O+fkuaeidN/awrvIGlo8AH42I/ytPqFRRBtIXX6l392mkIsGa15d9G/h2XvefSY2i/8za4UVW9rBcbyrl4YqIeGt5wvwcbyXr7tNAmSBpbCHATSM1OH6CdAe5JWt7z5jG2h5z+tsQttLyjwDfi4izymfkZ1y99oDUj+33eD4K1owQoNQF4Mbk6ywifgH8QlJpwNOvkO5+rQX5zm1oeA/prmQ70mCZO5GeCfyDVMmkN9cCz0j6mqQNJLVJ2kHSqyukPRN4l6S353QdSm2opkrqkjRDaXDOZaReJGqptn488L3CA/xJSoOEFn1L0hhJ25Oem5xTz/okvVlpWJk20nPCFcCqiHic1N3VT5SGbBoh6cWS3lTlPswjPbMsuRB4qaQPSRqVX6+WtG2kKv2/B47M+7QdNVQCqtO3lYb5eSOp4+Xzcj7OJR2rcfl4fZF0rkv7NFW1DbFTtIDUC1HxuBxP6vpue1hToaXUp+KFwBRJn5fUnvP0mjq3Xa7H81FI806tHVLoO8A1EfFITveafIf/LOnZ33ptlmG1cXAbGg4kPSN7OFKnunMjjSt2HGncph7v0POX27tIAfEB0i/5k4DxFdI+QupP7+ukL61HWDs67ghSJYDHSEU5byI9B6rWz0n9VP5V0mJSZZDyL7UrSBUR/gYcExF/rXN9U4DzSYHtjrze0pf5h0nFdLeTuuQ6n/S8rtp92FepBuAvIvU5+DZSD+uPkYrUjmZtf4ifIRXnzSU97zm1yu3UYy65r01SBZBDIqJ0x/hZ0hf2/aT+DX9D6gIJ4DJSt2BzlfphrUkuuv0eaQTnpyW9NiIuIB2Hs5Vqt94KvCOnX0yqqPSunOd7SCOc91sV5wPSvh9BuoZfRSq2hNQ11q9Ix/AhUkWqYxqRLxsY7n7LWl4uOnyA1I9iPcWFZn2SdBqpotU3m50X6z/fuZmZ2ZDj4GZmZkOOiyXNzGzI8Z2bmZkNOcOmndvEiRNj+vTpzc6Gmdmgcf311z8REZOanY96DJvgNn36dGbPnt3sbJiZDRqSau01p2W4WNLMzIYcBzczMxtyHNzMzGzIcXAzM7Mhx8HNzMyGnKYFN6Wh7m+RdKOkdaox5mHdfyHp3jyk+ysL8w6UdE9+DXRP6mZmNsg0uynAmwuDa5Z7B/CS/HoN8EvSQIMbk3rt7iaN1XS9pJkR8dT6yLCZmbW+Vi6WfDdwRiRXkwbe3BR4O3BJRCzMAe0SYK+BysSJf7+Pv942t++EZmbWMpoZ3II01tb1kg6uMH9zXjhE/Jw8rafp65B0sKTZkmYvWLCgrkyeftVDXHSrg5uZ2WDSzOC2a0S8klT8+GlJu5XNV4Vlopfp606MODEiuiOie9Kk+nqQ6epsZ97ipXUta2ZmzdG04BYRj+W/84ELgF3KkswBtii8n0oaPben6QOiq7ODuYsc3MzMBpOmBDdJYyWNK/1PGvr91rJkM4EP51qTrwUWRcTjwMXA2yRNkDQhL3vxQOW1q7OD+c8sG6jVm5nZAGhWbcku4AJJpTz8JiIuknQIQEQcD/wZeCdwL/Ac8JE8b6Gk7wDX5XUdFRELByyjnR0sXraSZ5etZGx7syuXmplZNZrybR0R9wOvqDD9+ML/AXy6h+VPAU4ZsAwWdHW2AzB/8TK2cnAzMxsUWrkpQEvo6uwA8HM3M7NBxMGtD6XgNt81Js3MBg0Htz6UiiXnPePgZmY2WDi49WHD9pGMGd3G3EWuMWlmNlg4uPVBEl2dHW7IbWY2iDi4VaGrs535LpY0Mxs0HNyq0NXZwTw35DYzGzQc3KrQ1dnB3GeWkpremZlZq3Nwq0JXZwfLV65m0fMrmp0VMzOrgoNbFdY2B3DRpJnZYODgVoU1vZS4UomZ2aDg4FaFrnEpuLkht5nZ4ODgVoXJpc6THdzMzAYFB7cqdIxqY6Mxo/zMzcxskHBwq1LXuA4/czMzGyQc3KrUNb7DxZJmZoNEU0fflNQGzAYejYh9yub9DHhzfjsGmBwRG+V5q4Bb8ryHI2LGQOe1a1w7d89dPNCbMTOzBmj20NKHAncAneUzIuILpf8lfRbYuTD7+YjYaeCzt1ZXZwcLlixj1eqgbYTW56bNzKxGTSuWlDQV2Bs4qYrk7wd+O7A56l1XZzurVgdPLnGlEjOzVtfMZ27HAl8FVveWSNKWwFbAZYXJHZJmS7pa0nt6WfbgnG72ggUL+pXZUkNu15g0M2t9TQlukvYB5kfE9VUk3x84PyJWFaZNi4hu4APAsZJeXGnBiDgxIrojonvSpEn9yvPa4OZKJWZmra5Zd267AjMkPQicDewh6cwe0u5PWZFkRDyW/94PzOKFz+MGhLvgMjMbPJoS3CLisIiYGhHTScHrsog4oDydpG2ACcA/C9MmSGrP/08kBcrbBzrPEzcczQi5lxIzs8Gg2bUlX0DSUcDsiJiZJ70fODteOJDatsAJklaTgvMPI2LAg9vIthFM3LDdz9zMzAaBpge3iJhFKlokIg4vm3dkhfRXATuuh6ytozRoqZmZtTb3UFKDrs52VygxMxsEHNxq0NXZwfzFLpY0M2t1Dm416OrsYOGzy1m2clXfic3MrGkc3GrQtWZcN9+9mZm1Mge3GkzObd3mL/ZzNzOzVubgVoMp7oLLzGxQcHCrwZpeShb5zs3MrJU5uNVgwphRjG4bwTwXS5qZtTQHtxpIYnJnuyuUmJm1OAe3GnV1drght5lZi3Nwq1FXZ7u74DIza3EObjWaPK7DxZJmZi2u5uAmaZWkyRWmbyJpyHfdMWV8B0uWrWTJspXNzoqZmfWgnjs39TC9HVjej7wMCqVeSvzczcysdVU95I2kL+Z/AzhE0pLC7DbgjcCdDcxbS+oaV2rIvZQXT9qwybkxM7NKahnP7bP5r4CPA8UiyOXAg8AhtWxcUhswG3g0IvYpm3cQ8GPg0TzpuIg4Kc87EPhmnv7diDi9lu32R9f43AWXn7uZmbWsqoNbRGwFIOly4L0R8VQDtn8ocAfQ2cP8cyLiM8UJkjYGjgC6SXeR10ua2aD89Kmrc+2dm5mZtaaan7lFxJsbEUgkTQX2Bk6qcdG3A5dExMKcj0uAvfqbn2pt2D6SsaPb3BzAzKyF1VIsuYak/YC3AJMpC5ARMaPK1RwLfBUY10ua90naDbgb+EJEPAJsDjxSSDMnT1tvujrdHMDMrJXV0xTgx8CZwHTgaeDJslc169gHmB8R1/eS7H+B6RHxcuBSoPRcrVJtzehhOwdLmi1p9oIFC6rJWlXcS4mZWWur587tw8D7I+L8fmx3V2CGpHcCHUCnpDMj4oBSgogoBspfAUfn/+cAuxfmTQVmVdpIRJwInAjQ3d1dMQDWo6uzndkPrZdHfGZmVod62rmNAG7sz0Yj4rCImBoR04H9gcuKgQ1A0qaFtzNIFU8ALgbeJmmCpAnA2/K09aZULBnRsHhpZmYNVE9wOxE4oM9UdZB0lKTSM7vPSbpN0k3A54CDACJiIfAd4Lr8OipPW2+6OjtYvmo1Tz+3Yn1u1szMqlRPseRGwAckvRW4GXjBN3xEfK6WlUXELHKxYkQcXph+GHBYD8ucApxSy3Yaac2gpc8sZcLY0c3KhpmZ9aCe4LYda4slX1Y2b1iU0xW74Np2056a6JmZWbPUHNwi4s0DkZHBpHTn5uYAZmatyUPe1GGyO082M2tpNd+5SZrZ2/waGnEPWu0j25gwZpR7KTEza1H1PHMrb6g9CngFsAXw+37naJBIDbldLGlm1orqeeb2kUrTJf0EWNzvHA0Skzs7mL/Yd25mZq2okc/cTgA+1cD1tbQpne3MXeTgZmbWihoZ3LZp4LpaXldnB08sWcbKVaubnRUzMytTT4WSX5RPAjYF3kETG1avb5M7O1gd8OSzy9c0DTAzs9ZQT4WSHcverwYWAF9gGAW3KYVBSx3czMxaixtx16nUS8ncRUt5+dQmZ8bMzF6grsFKASR1AFuTuty6LyKGVe2K0t3avMVuDmBm1mrqGax0VB6w9CngJuAW4ClJP5I0qtEZbFWbjB3NCMF8N+Q2M2s59dy5HQ28HzgEuDJPeyPwA1Kw/HJjstbaRraNYNI4NwcwM2tF9QS3DwAfjYg/F6bdJ2kBcBLDJLhB7qXExZJmZi2nnnZu44H7Kky/jzTW27AxecoLXNcAAB8YSURBVFyHiyXNzFpQPcGtNDJ2uUNZO85bVSS1SbpB0oUV5n1R0u2Sbpb0N0lbFuatknRjfvXakfNAmjK+3SMDmJm1oHqKJb8K/DmPxP1PUm3J1wGbkRpy1+JQ4A6g0oifNwDdEfGcpE8CPwL2y/Oej4id6sh7Q3WN6+Cp51awdMUqOka1NTs7ZmaW1XznFhF/J3W1dR6wISkwnQdsExFX9rZskaSpwN6k53SVtnN5RDyX314NtFxrslJzgAV+7mZm1lLqaucWEY8C3+jnto8l3QWOqyLtx4C/FN53SJoNrAR+GBF/qLSQpIOBgwGmTZvWv9xWUBy0dIuNxzR8/WZmVp962rl9RtIBFaYfIKmqUQEk7QPMj4jrq0h7ANAN/LgweVpEdJNqbh4r6cWVlo2IEyOiOyK6J02aVE3WajJlfLpz86ClZmatpZ4KJZ8HHqkw/UFS/5LV2BWYIelB4GxgD0lnlieStCfpDnFGRKwp+4uIx/Lf+4FZwM7VZ79xusaV+pd0saSZWSupJ7hNBR6qMH0OVT4Xi4jDImJqREwH9gcui4gX3A1K2pk0RtyMiJhfmD5BUnv+fyIpUN5ex37020ZjRjF65Ag3BzAzazH1BLe5QKWaiq8EnuhPZiQdJWlGfvtjUoWV88qq/G8LzJZ0E3A56ZlbU4KbJLo63RzAzKzV1FOh5DfALyQ9SyoSBHgzqYLIWbWuLCJmldYTEYcXpu/ZQ/qrWHfYnabpGtfhZ25mZi2mnuB2BLAVcDGwKk8bQWoO8K0G5WvQ6Ors4I7Hn2l2NszMrKCe8dxWAO+X9C1SRQ4B/4qIexuducFgcmc7s+7ynZuZWSupezy3HMyGZUArmtLZwbPLV7F46QrGdQybEX/MzFpaPRVKrGDNoKVuDmBm1jIc3Pqp1EuJmwOYmbUOB7d+WnPnttjBzcysVTi49VMpuM1d5GJJM7NWUVdwk9Ql6cuSfpl7CUHSrpK2amz2Wt+G7SPZsH2kG3KbmbWQejpOfhVwF/BBUm/9pbHY3gp8r3FZGzwmd7Yz38WSZmYto547t2OAn0fEzkCxLO5iUj+Pw86Uzg7mLnJwMzNrFfUEt1cBp1eY/jjQ1b/sDE5dnR1uCmBm1kLqCW7PAxMqTH8ZML/C9CGvVCwZEc3OipmZUV9w+yNwRGnYGSAkTQeOBn7XoHwNKl3jOlixKnjquRXNzoqZmVFfcPsysDGwABgDXEnqhutp4JuNy9rgsWZEbj93MzNrCfV0nPwM8AZJe5DGcBtB6jj50kZnbrDoyr2UzFu8lO3WVB41M7Nm6U/HyZcBlzUwL4PW5HHpzs1dcJmZtYaqgpukL1a7woj4abVpJbUBs4FHI2KfsnntwBmk2plPAvtFxIN53mGkNnargM9FxMXVbnMglPqXdC8lZmatodo7t89WmS6AqoMbcChwB1Qsy/sY8FREbC1pf1KFlf0kbQfsD2wPbAZcKumlEbGqwjrWi/aRbWw8drT7lzQzaxFVBbeIaHi3WpKmAnuTejWpdGf4buDI/P/5wHGSlKefHRHLgAck3QvsAvyz0XmsxeRx7Vx2x3w++ez1zcyGDUJv2baLfV81tdnZMBtS6n7m1gDHAl8FxvUwf3PgEYCIWClpEbBJnn51Id2cPG0dkg4GDgaYNm1aY3Ldgxk7bcYfbniU+xYsGdDt2NDy+KKl3DN/iYObWYPVHNzy3dMngU8DWwE7RMT9kv4fcH9EnFvFOvYB5kfE9ZJ27ylZhWnRy/R1J0acCJwI0N3dPaAtrD+1+9Z8avetB3ITNgQd8cdb+f0NjzY7G2ZDTj3t3A4ltWc7kRcGmkeBz1S5jl2BGZIeBM4G9pB0ZlmaOcAWAJJGAuOBhcXp2VTgsdp2waw1TO7sYPHSlTy3fGWzs2I2pNQT3A4B/jMifg4UP5H/IlXy6FNEHBYRUyNiOqlyyGURcUBZspnAgfn/fXOayNP3l9Seh9h5CXBtHfth1nRTOkvNSFzT1qyR6nnmtiVwa4XpK4AN+pMZSUcBsyNiJnAy8OtcYWQhKQgSEbdJOhe4nRRcP93MmpJm/bFmsNtnljJ94tgm58Zs6KgnuN1P6pnkobLp7yQFnJpExCxgVv7/8ML0pcC/97DM9ximY8fZ0LKmdxt3AGDWUPUEt2NI1fLHkJ65vU7Sh0g1Hz/ayMyZDXWTXSxpNiDq6Vvy1FzB4/ukjpN/TapM8rmIOKfB+TMb0jo7RrLBqDbm+s7NrKHqaucWEb8CfiVpIjAiIoblOG5m/SWJrs52F0uaNVi/GnFHxBONyojZcDW5s8PFkmYNVm3HybfQQ0PpchHx8n7lyGyYmdLZwU1znm52NsyGlGrv3M4f0FyYDWNdne3MXbSUiCB1AGRm/VVtx8nfHuiMmA1XXZ0dLFu5mmeeX8n4MaOanR2zIaGeHkrMrIFKzQE8ZJJZ4/iZm1mTlbrgmvfMUl7a1dMgGWZWCz9zM2uyrjUjufvOzaxR/MzNrMlK/UvOX+zmAGaN4mduZk3WMaqN8RuMckNuswaqqxG3pI8A7wemAaOL8yLiRQ3Il9mwUmoOYGaNUfOdm6SvAD8BrgemA38gDYGzMXBKIzNnNlx0dXYwz8WSZg1TT7HkfwIHR8RhpDHcjouIGaSAt2UjM2c2XEwe18F8F0uaNUw9wW0qa0e+fh7ozP//FnhfNSuQ1CHpWkk3SbpN0joVViT9TNKN+XW3pKcL81YV5s2sYx/MWsqU8e3MX7yM1auranFjZn2o55nbXGAi8DBpwNLXATcCW1NlWzhgGbBHRCyRNAq4UtJfIuLqUoKI+ELpf0mfBXYuLP98ROxUR97NWlJXZwerVgdPPLuMyeM6mp0ds0Gvnju3y4AZ+f+TgZ9Kuhw4B/h9NSuIZEl+Oyq/eguM7yfdGZoNSaWA5tEBzBqjnju3g8lBMSKOl/QUsCvwO+CEalciqY1UKWVr4L8j4poe0m0JbEUKqiUdkmYDK4EfRsQf6tgPs5YxZfzaXkp22Hx8k3NjNvjVMxL3amB14f05pLu2WtezCthJ0kbABZJ2iIhbKyTdHzg/py+ZFhGPSXoRcJmkWyLivvIFJR1MCsZMmzat1iyarTelXkrm+c7NrCGa3og7Ip4GZgF79ZBkf8qKJCPisfz3/rzszusuBhFxYkR0R0T3pEmTGpVls4abuGE7Esx1jUmzhmhKcJM0Kd+xIWkDYE/gzgrptgEmAP8sTJsgqT3/P5FUJHr7+si32UAZ1TaCTca2uzmAWYPU1UNJA2wKnJ6fu40Azo2ICyUdBcyOiFL1/vcDZ0dEsbLJtsAJklbnZX8YEQ5uNuhNGd/uLrjMGqQpwS0ibqZCUWJEHF72/sgKaa4CdhywzJk1Sde4Dh5zF1xmDVFP91vTe5n3+v5kxmw4m9zpXkrMGqWeZ243STqgOEHSiFykeHljsmU2/HR1tvPks8tZvnJ134nNrFf1BLevAcdL+o2kTkkvJlX4+Ciwd0NzZzaMlEbkXrDEzQHM+qvm4BYRxwPdwDak0QBuAOYAL4+ISxubPbPhozRoqYe+Meu/epsCzAUeBLqADYCLImJhozJlNhxNzg25/dzNrP/qqVCyG3AzaXSA7UnFkT+WdIGkTRqcP7Nho1Qs6eYAZv1Xz53bpcAZwK4RcW9E/BrYCZgE3NLIzJkNJxPGjGZUmzxoqVkD1NPObc+I+HtxQkQ8mO/ovt6YbJkNPyNGiMnjOpjnZ25m/VZPx8l/72H6auC7/c6R2TA2ubOdeYsd3Mz6q64eSiRtTOroeBowujgvIo5qQL7MhqUpnR3cM39J3wnNrFc1BzdJrwX+RBpNexLwKKmvyGWkGpQObmZ16urs4Mp7nmh2NswGvXoqlPwYOAvYHFgK7EG6g5sNHN24rJkNP5M721m8bCXPLlvZ7KyYDWr1BLeXA8flnvpXAe0RMY/Uc8mRDcyb2bBTag4w3zUmzfqlnuC2vPD/PGDL/P8SYLN+58hsGOtyWzezhqinQsm/gFcDd5NGwf6upC7gAFLjbjOrU1fupcTBzax/6rlz+wbwWP7/m8AC4L9II2Yf3KB8mQ1Lk33nZtYQ9XScPDsiLs//L4iId0REZ0R0R0RVPZRI6pB0raSbJN0m6dsV0hwkaYGkG/Pr44V5B0q6J78OrHUfzFrVuPaRjBndxrxn/MzNrD+aMhI3qdnAHhGxRNIo4EpJf4mIq8vSnRMRnylOyG3sjiCNTBDA9ZJmRsRT6yXnZgNIEl2dHb5zM+unqoObpJnVpIuIGVWkCVIFFIBR+RVVZuXtwCWlUQgkXUJqUP7bKpc3a2mTx7U7uJn1Uy3FkvsAOwJP9vGqiqQ2STcC80nB6poKyd4n6WZJ50vaIk/bHHikkGZOnlZpGwdLmi1p9oIFC6rNmllTTRnf4WJJs36qpVjyGFKNyN2AU4HTImJOvRuOiFXATpI2Ai6QtENE3FpI8r/AbyNimaRDgNNJDcZVaXU9bONE4ESA7u7uau8MzZqqVCwZEUiVLncz60vVd24R8VVgC+ALpOdd90j6i6R983OzukTE06QmBXuVTX8yIko/X38FvCr/Pyfno2Qqa2tvmg16k8e1s2zlahY9v6LZWTEbtGqqLRkRqyJiZkS8B9gKuJw0EsCjkjasdj2SJuU7NiRtAOwJ3FmWZtPC2xnAHfn/i4G3SZogaQLwtjzNbEhY25DbRZNm9epPbcmxwEbAhqTKIbUU+20KnC6pjRRgz42ICyUdBcyOiJnA5yTNAFYCC4GDACJioaTvANfldR1VqlxiNhRMGb+2rds2U8Y1OTdmg1NNwS3fZf0H8DFS0eQFwIER8bda1hMRNwM7V5h+eOH/w4DDelj+FOCUWrZpNlh0jXNDbrP+qqUpwInAfsA9wMnAjPy8zMwaaLK74DLrt1ru3D4OPAw8DrwDeEelmlzVtHMzs551jGpjozGj/MzNrB9qCW5nUNtzNTOrU9c491Ji1h9VB7eIOGgA82FmBZM725nnMd3M6lbPqABmNsC6OjuYt8h3bmb1cnAza0FTOjtYsGQZq1b7SYBZPRzczFpQV2c7q1YHTz7rokmzeji4mbWgNYOWLnJwM6uHg5tZC5riEbnN+sXBzawFrelfcrGDm1k9HNzMWtDEDUcjufNks3o5uJm1oJFtI5i4YbubA5jVycHNrEVN6exwsaRZnRzczFpUV2e7iyXN6uTgZtaiJnd2MN+1Jc3q4uBm1qKmdHbw5LPLWbZyVbOzYjboNCW4SeqQdK2kmyTdJunbFdJ8UdLtkm6W9DdJWxbmrZJ0Y37NXL+5N1s/uvK4bgvcgbJZzWoaibuBlgF7RMQSSaOAKyX9JSKuLqS5AeiOiOckfRL4EWmwVIDnI2Kn9Zxns/VqTS8lzyxj6oQxTc6N2eDSlDu3SJbkt6PyK8rSXB4Rz+W3VwNT12MWzZqua5x7KTGrV9OeuUlqk3QjMB+4JCKu6SX5x4C/FN53SJot6WpJ7+llGwfndLMXLFjQoJybrR9Txju4mdWracEtIlblosWpwC6SdqiUTtIBQDfw48LkaRHRDXwAOFbSi3vYxokR0R0R3ZMmTWrwHpgNrAljRjGqTW4OYFaHpteWjIingVnAXuXzJO0JfAOYERHLCss8lv/en5fdeX3k1Wx9ksTkcW4OYFaPZtWWnCRpo/z/BsCewJ1laXYGTiAFtvmF6RMktef/JwK7Arevr7ybrU9dne3MdXAzq1mzaktuCpwuqY0UYM+NiAslHQXMjoiZpGLIDYHzJAE8HBEzgG2BEyStzsv+MCIc3GxImjK+g7vmLm52NswGnaYEt4i4mQpFiRFxeOH/PXtY9ipgx4HLnVnrmDyug3/c/USzs2E26DT9mZuZ9ayrs4PFy1by7LKVzc6K2aDi4GbWwqaMT72UuDmAWW0c3Mxa2NqG3G4OYFYLBzezFlbqgmu+x3Uzq4mDm1kLK3WePNcjcpvVxMHNrIWN6xjF2NFtLpY0q5GDm1mL6+rsYJ6LJc1q4uBm1uImd7a7Cy6zGjm4mbW4KZ0d7oLLrEYObmYtrquzg3nPLCMi+k5sZoCDm1nLm9zZwfKVq1n0/IpmZ8Vs0HBwM2txpeYArjFpVj0HN7MWNyU35PZzN7PqObiZtbiuzlIXXA5uZtVycDNrcZPGpWJJNwcwq56Dm1mL6xjVxoQxo1wsaVaDpgQ3SR2SrpV0k6TbJH27Qpp2SedIulfSNZKmF+YdlqffJent6zPvZs1Qag5gZtVp1p3bMmCPiHgFsBOwl6TXlqX5GPBURGwN/Aw4GkDSdsD+wPbAXsD/SGpbbzk3a4LJnR0uljSrwchmbDRSa9Ql+e2o/Cpvofpu4Mj8//nAcZKUp58dEcuAByTdC+wC/HOg823WLF3j2vnnfU/w1p9e0eys2CAzqm0Efz70jc3OxnrXlOAGkO+2rge2Bv47Iq4pS7I58AhARKyUtAjYJE+/upBuTp5WaRsHAwcDTJs2raH5N1uf9nv1Fjy3fBWxzm9As96NHDE8q1Y0LbhFxCpgJ0kbARdI2iEibi0kUaXFepleaRsnAicCdHd3+1vBBq3u6RvTPX3jZmfDbNBoekiPiKeBWaTnZ0VzgC0AJI0ExgMLi9OzqcBjA55RMzMbNJpVW3JSvmND0gbAnsCdZclmAgfm//cFLsvP6mYC++falFsBLwGuXT85NzOzwaBZxZKbAqfn524jgHMj4kJJRwGzI2ImcDLw61xhZCGphiQRcZukc4HbgZXAp3MRp5mZGQAaLsNodHd3x+zZs5udDTOzQUPS9RHR3ex81KPpz9zMzMwazcHNzMyGHAc3MzMbchzczMxsyBk2FUokLQAeqnPxicATDcxOMw2VfRkq+wHel1Y0VPYD+rcvW0bEpEZmZn0ZNsGtPyTNHqw1hsoNlX0ZKvsB3pdWNFT2A4bWvtTCxZJmZjbkOLiZmdmQ4+BWnRObnYEGGir7MlT2A7wvrWio7AcMrX2pmp+5mZnZkOM7NzMzG3KGXXCTtJekuyTdK+n/VZjfLumcPP8aSdML8w7L0++S9PZq1zlQBmhfHpR0i6QbJa23zjjr3RdJm0i6XNISSceVLfOqvC/3SvpFHsl9MO7HrLzOG/Nr8kDvRz/35a2Srs/H/npJexSWWe/nZAD3ZbCdl10Keb1J0r9Vu85BKSKGzQtoA+4DXgSMBm4CtitL8yng+Pz//sA5+f/tcvp2YKu8nrZq1jlY9iXPexCYOIjOy1jgDcAhwHFly1wLvI40wO1fgHcM0v2YBXQPonOyM7BZ/n8H4NFmnZMB3pfBdl7GACPz/5sC80kjwzTlO2ygX8Ptzm0X4N6IuD8ilgNnA+8uS/Nu4PT8//nAW/Kvy3cDZ0fEsoh4ALg3r6+adQ6WfWmWuvclIp6NiCuBpcXEkjYFOiPin5E+zWcA7xnQvRiA/Wii/uzLDRFRGkD4NqAj300045zAAOzLeshzT/qzL89FxMo8vQMoVbho1nfYgBpuwW1z4JHC+zl5WsU0+UJYBGzSy7LVrHMgDMS+QLrg/5qLYA4egHxX0p996W2dc/pYZ6MNxH6UnJqLk761noryGrUv7wNuiIhlNOecvCCfvWy31n0pGVTnRdJrJN0G3AIckuc36ztsQA234Fbp4iuvLtpTmlqnD7SB2BeAXSPilcA7gE9L2q3+LFatP/vSn3U22kDsB8AHI2JH4I359aE68larfu+LpO2Bo4FP1LDOgTAQ+wKD8LxExDURsT3wauAwSR1VrnPQGW7BbQ6wReH9VOCxntJIGgmMJ40E3tOy1axzIAzEvlAqgomI+cAFrJ/iyv7sS2/rnNrHOhttIPaDiHg0/10M/IZBcE4kTSVdPx+OiPsK6df3OXlBPnvZbq37MijPS0lE3AE8S3qO2KzvsAE13ILbdcBLJG0laTTpYevMsjQzgQPz//sCl+XnAzOB/fOzg62Al5AejlezzkGxL5LGShoHIGks8Dbg1hbfl4oi4nFgsaTX5uKiDwN/bHzWX6Dh+yFppKSJ+f9RwD60+DmRtBHwJ+CwiPi/UuImnRMYgH0ZpOdlqxzskLQlsA2pAlmzvsMGVrNrtKzvF/BO4G5S7aBv5GlHATPy/x3AeaRKFtcCLyos+4283F0UanlVWudg3BdSbamb8uu2QbQvD5J+mS4h/QrdLk/vJn3h3AccR+60YDDtB6kW5fXAzfmc/Jxcs7VV9wX4Jumu4MbCa3KzzslA7MsgPS8fynm9EfgX8J7e1jnYX+6hxMzMhpzhVixpZmbDgIObmZkNOQ5uZmY25Di4mZnZkOPgZmZmQ46Dmw0oSUdKqqn9j6TTJF04UHkaqvKxnicpJB3U7PyYNZODm9UkB56QdFKFeT/K84qB6RjgTTVu5lDggP7kc32T9ApJf5Q0V9JSSQ9L+l1uLFtKE5KWS3pR2bIvCOY5SEV+rZL0iKSTJE3qZfs7AEeQRhXYFDin8XtpNng4uFk9HgH2y72YAGu6+fkQ8HAxYUQsiYgna1l5RCyKiKcbktP1IAedv5EaX+8NvIx0LO4DOsuSrwK+V8Vq7yIFqWnAJ4F3kXrR78nW+e8fImJuRDxfIZ+jq9iu2ZDg4Gb1uBm4B/iPwrS9ScO1zComLC+WLN2lSDpU0qOSnpJ0qqQx5WkK72dJ+qWkn0haKGlBXr5d0n9LejrfKX2osMz0fOfTXZafkLRvWZr9JV0h6XlJN0h6uaQdJF0l6VlJV+ZuynqyKzAB+EhEXB8RD0bEFRHx1Yi4pSztfwH/IelVvR5hWJmD1KMRcSHwC+BtkjYoTyjpSFLfhwCrJUXZsf6apDnkHvkljZZ0tKQ5ef+uU2HA2pxmL0l35rvQf0j6QD5W0/P8gyQtKVtm95xmYmHa6/OxfS6f719K6izMnyXpfyR9X9ITkuZLOkbSiEKa0Xn+Q5KWSbpf0ueU3Cvpy2X5eEnOxyv7OMY2hDm4Wb1OBj5aeP9R4FSq6038jaQOW/cE9gP+jVQU2ZsPAouB1wA/BI4F/kDqMqibNH7VSZI2q34X1vg2qcf3nYGnSZ3g/hepi7JdSN0Z/aKX5eeSPkv7Sn0Oe3It8DvgRzXm8fm8jZEV5h0D/Gf+f9P8KnkT8HJgL+AtedqpefoHgB1Jx+5/Jb0CQNIWpGN7CbAT6VjUml8k7Qj8ldRP4SuA9+b1nVKW9IPASuD1wGeAz5Oui5LTSf1QfhHYFvgY8HSk7pXKr0Py+xsj4l+15tmGkGb3/+XX4HoBpwEXku5Unid1ujwFWEYqQjsNuLCQ/kjg1rLlHyGPCJyn/Qq4tHwbhfezgH8W3gtYAMwsTBsFLAf2ze+nkwJtd1n+o0KaTxTm75Onvbcw7SBgSR/H5XvACuAp0hf614EtK22bVIS4HNirh/0tP2YvI90pX9PL9vdNH+d1ztUCoL0w7cXAamBaWdo/AP+T//8+6UeDCvO/mfM/vadjAuye00zM788ATi5Ls1NOU+pr8gXnNk+7BDgp//+SnH6vHvZ7Sj7ur83v24BHgc80+7PiV3NfvnOzukTEU6SisI+SeiCfFREP977UGrfH2hGBIQ2vMbmPZW4ubDuA+aQBF0vTSoGlr/X0um5gXv57S9m0scWi03IR8Q3SF+3BedmPAbdLekuFtPeSAvoPi8VvZbaVtETS88DtpB8EH6xyf4pujRcOrvlK0o+D2/P6l+Tixb1JgQ/S3dHV+TiX/LOObb8KOKBsO6We9V9cSHdz2XLF62FnUjC+vNIGImIu6cdW6e5tL9LAnGfVkV8bQioVcZhV6xRSkdES4PAalltR9j7ou4i80jK9rWd1/rummFBpaJK+1h29TOs1j5EqzpwHnCfpMOAG4Fukyiblvk2qcNJTwLqP1FP7KuCxsgBVi2fL3o8g7c+rWff4lSqhVDOi9OoK6cqP7wjgJOBnFZZ/tPB/b+exmrycBPxG0udJQe73+ceXDWMObtYffyMVr00kFWu1kgX5b/H5007ra+MRsVzSfUDFZ4ARMV/SMcB3qHxXtDzf4TXaDaSAMSUiKt4Nke4U3ydJhbu315alWQCMkdQZEc/kaeXH91/A9v3cj3+RAt2bgYt6SHMR8AypGcS7SD8KbJhzsaTVLX/xvRzYqh93FgMiUlX4q4GvSdpe0utJFS8aTtI+ks7Mf18qaZtcg++drK3FWMlPSJVV3jMQ+aokIu4mFdmdJmlfSS+S1C3py5Lem5MdT3oeeWzel31JgaPoGtJd4Q8kbS3pfcCnytIcDewi6XhJO+d0+0g6oYb83gOcS6os9D6lATXfWKwZGxGrSKUIPyDdEVa6U7ZhxsHN+iUiFhd+ubea0nOY64ATSJUiBsLtpKLZY0h3RteSGqF/mVQ5o6KIWEIqnuwYoHz15COkGpM/Au4kPbPaDXgo5+thUs3GvUgD134B+H/FFUTEQlKR6ltJzxgPJhXBFtPcnNc7Hbgir+sHrH2uWa0Pk2qw/iLn9zRgfFmaU4DRwKllzwptmPJgpWbWp9xe8DrSXfqDTc7OOiS9hlRZ5UU1VGyyIczP3Mxs0JLUDmwBfBe4wIHNSlwsaWaD2ftJXZVtQmrkbQa4WNLMzIYg37mZmdmQ4+BmZmZDjoObmZkNOQ5uZmY25Di4mZnZkOPgZmZmQ87/B8H01NV57gOuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y)\n",
    "plt.xlabel(\"Minimum SNP frequency\",fontsize=14)\n",
    "plt.ylabel(\"Max allele count\",fontsize=14)\n",
    "plt.suptitle(\"Maximum allele count across minimum cutoff values\",fontsize=16)\n",
    "plt.title(\"Alleles represented in both tech reps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_r = [i/1000 for i in range(0,31)]\n",
    "y_r = [max(get_allele_counts(min_snp_freq_data(data,j))) for j in x_r]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Alleles represented in either tech rep')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEnCAYAAAAq8Q2oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5xcZb3H8c83ZRN2JSQLSQglhIgCARQwUkQQkHhREAtcC0UQFbk27IoN9IKKYkcFBEGKiGLjAiI1CIJAQHqPhA4JJJRQQgi/+8fzTDIMs7szu7NzZna/79drXrtzynOeU+b8TnnO7ygiMDMzG0pGFF0BMzOzRnNwMzOzIcfBzczMhhwHNzMzG3Ic3MzMbMhxcDMzsyGnz+AmaT9JkT+vrtJ/+7L+Ow1GJcumsf1glN+Kqs2zpNmSZvejrNI6nNa4GrYmSYdK2rHoegwFkqbl7Wa/fow7bLa5ekjaNG+j3QMoY2tJV0p6Oi/jTXP3r0i6V9ILkq5rXK17rcuhklryebJ6ztyeAvap0v0Dud9guhbYOv81680hgINbYzxE+t2d3Y9xz87jPtTQGrW/TUnbaL+DG3A8MAp4O2kZ3yFpC+Bw4HfAdlTfVw8r9QS3PwF7S1Kpg6SVgN2BPza6YuUi4smI+FdEPDmY0zFrFCUdRddjICJiSf7dLejHuAvyuEsGo27DlaQRwPrA2RFxUV7GzwAb5kGOjojLI+LG4mrZIiKi1w+wHxDAm4EXgW3L+u0JLAZ2zcPsVNbv9cAZwP3As8DtwLeBlcqG2Tj3+1HFNL8NLAE2y9+3z+VvXzbMbOAyYGfgulzOv4EtSUc13yYdNS4ETgS6ysZ9WXkV8zqtrNs84BTSkdDteTqXAq8CuoBjgMeAR4AfAKNqWKbfJJ2FPgE8ClwEbFUxTE/zPLtiuNWAXwIP5GV2G3BAX/OVu38EuB54LtfjeKC7r/qXjXttXh6LgEuAN5T1nwKclMtdAtwA7F1RxqFpE3xZ2ScC88q+T8v1/yjwrbxeHwf+D1irbLio8jm0l3lYDzgZuDvPx3/yspxQZdg3AefndfZ0Xm4fqrKd7J/XwVLgXXUsi9WB3wAP5mEeAs4CJuX+o4D/BeaWra/LgDf2sZ5m0//fSWm571exbu4HNiP9Dp4B7gQOHMzfUrXyetqG8nCHAZ8D7snr62xgUv78Pq/H+4Av1bi9TwR+kcdZkv+eDIypts1W+82WzUPlZ1ruPw44qmwbuB34DKBexp+Xp1HTdg+cA1xTpfsU4AXg02XzewxwR17H9wG/BdbsbflX22b62Oe+G/hXnsbjwB+AqRXD7EnaZhfn9XYj8NG+1tkoancP8A/Shnlp7vYB4M95opWmkn5MJ5IuW24EfAOYDryPtERukvQ54ChJ50XE3yTtAHwJ+EJE/LuPOq0HfJ90Or4Y+B5wZv6MIm0MG+Zh5gNfrGN+y20HvDLXqwP4Mels9T/AXXl+tgO+Rtr5/KKP8tYEfkTaSXQBewP/kDQzIm6otVKSxgH/BFYibWR3A/8F/FLSmIj4WS/jfpf04/8p8IVcp8OAjSW9ISKW9TLukXnc40mXWF4EtiKt88sldZGC3QTgK6Qfxt7AyZI6I+LYWuexwsHA5aQAMom0AzyVFHggXaK5grTNHZO73d9LeWvk/p8mBejpub7n5LJK8/sO0vr+JynAPkrantepKG8H0mWnb5K2t3l1LIuTc3lfyMNMJh1Qdub+XyLt6L5K+l2NA2ZS2+WtRv9OxpF2dD8mHWx8kLTN3R4RF/cxbqN/Sz3ZB7gJ+BhpWf6YdICxMvA34Fjgv4HvSroxIs7pqSBJE0jbXTfpN3IDaft7R56HWs9Oz87jfy1Pu7RtPpTPyM4GNiftJ28EdgF+SAo0X8n930g6WDkeOC5PewlpmzqYFCweouft/iTgNEkzIuKWsu575r+n5b/dpIOog4EFpN/K54B/StogIp6rcZ57JOlA0sHkCaTtaGXSfuwSSa+JiKckvZF0QFTaT40ANgDG9zmBGo5Y9iNF3PVIO5VFwFhWRPpZrIjKO/VQhkg/or1JO8JVK/r/hXS0thHpDORc8tFKT1GfdLSyFJhe1m23PNwFFeX/Cbi7hqOI0rxWHm0uBFYp6/apPNxxFeNfC1xcy5Fg2Tgj87K5HfhJDfM8u+z710kb4KsqyvwVaQc8qtp8kY6ulgHfqBhvmzzcO3up73p53B/2Mswneli+F5B2niOrHfWVDXci1c/cLqkY7vO5+xpl3QI4rJ51UDbuKNLOI1hx1UB5G5gDjOhl3Hmko8/V+7ksFgOf6qX8s4A/9WOeZtP/30lpue9XsW4C2KGs25i8vR07WL+lauX1tA3l4e7gpWd+P8zdv1axvucDJ/SxDL+Vt/nNehnmJdtsxfKfXWU+1qsYrnT1a7+K7qUAtlpZnV92ZgZ8uNryqVKflUhnP9+p6H4dcE4v440E1s7TeFdPy7/aNpO7b0/Z7wB4Ra7HryuGmwY8z4ozyM8DC+vd7iOi7kcB/kDakN8O7AU8DFxYbUBJ4yQdIWluXjlLSUenIl2GKPeh3P8a0srbN/Kc9eGOiPhP2ffb8t+/Vwx3G7BW+f3COl0REU/UOJ21+ypM0k6SLpb0GOkAYSnwatK19HrsDFwJ3C1pVOmT67UqMKOH8WaRjoBOrRjvSuBJ0pFzT3bK4/Z29rUd8EBEzK7ofgrpKLSnevWlsmFD6b7C1P4UJqkjtzC7TdKzpPVQuiqxftnfdUg73xf7KPJfEfFwRbdal8XVwBckHSRpkyrb6tXA2yQdLumNdd7Pa/Tv5JkoO0OLdF/tTmpbDw39LfXi/Ih4obfp5P531TCdtwBXR99XkgZiO9KB/2kV3U8hnR1u/bIx+iEiniWdKe9VWs+SNgFeSzqrW07S/0i6XtJi0n7q3tyr3v1UNVuTrgBU7oPuJ62r0j7oamCCpFMk7Sqp7zO2rK7gFhFPkc6y9iFdkjy1lx/8CcCBpNPJWaR7cB/P/cZWlPsYacc1BjgtIh6psUqLKr4/30v3UaSjj/6oZzpj6YWkzUmXvRaTgvpWpGVzfV/jVjGJtBEsrfj8IfdftZfxIP2wK8cd18t45WX2drmvm+qt5B4u698fCyu+ly4H1bvcSr5DOvI8hXQJaAvSZZ3yMmuZ35Jq81zrsngv6TLhF0mXvR6Q9I18uQrSvbFDSGddlwKPSTpB0mo11KvRv5PK8SCti1rWQ8N+S02czqrUtv4Hopt0hlJ5iXOgv5lqTiIF9O3z931It47+WhpA0idJl4QvIP0mtiDtq2Bg66WktA+6gJfvgzYh/+4i4hLSJdy1SbfAFki6QNJr+ppAPffcSk4iBaIRwPurDSBpLOl69KER8ZOy7pv0MPxOpAYKc4CPSTolIub0o261Kl0vrjz67W2n3ii7k46C3h0RS0sd83X9x+ss6zHSZZWDeuh/ey/jQToirbajeqxKt5JH8981eyl/IdWP7lavKP85SGdQEfF82XDNWA+Q7u+cFBGHlTpIekXFMOXz25dqVxtqWhYRMZ908PdxSesD+5Lu3S0Afpm3lSOAIyStTrqM9UPSPbn31lC3dlfkb/ZR+l7/z/HyukGqX2+/p5KFQHeV30Llb6YRLiGdhe0t6RLSfvyMfFZX8j7gwoj4XKmDpHVrKLvW9VSan/2Am6uUs/zxsog4Azgj/za3J/0OzpW0Vm9XU/qToeR8UmujoyOiWqUgnYGNJEXhcvtVDpiPPE8inc28gdQq5rdVdjKNdE/+u3FF97cN4jRLOknX75fvCPNDx/25tHYu6ebqvRExp8qnp+cPzyddApnaw3h39zLNC/K4B/QyzCWky1vbVHTfkxSMb83fX7Ye8mWHN/RSdl+eJ91XqEUnL99GP1jx/Q7SvaIP9/Oydq3LYrmIuD0ivkI68KjcRomIhyPiONK6eFn/IaratjKKdIA22M4DtpD02l6GuQeYXH4mLemVvPzApnRmVrmNXkLaH/93Rfe9SNv0v+qtdE/yLZ9TgT1I+7y1qLgkSW2/jWoeIc1j5Xa5S8X3y0kBbL0e9kEvO3COiMURcRapsdgU+jiwqfvMLVIruqpnbGXDPCHpX8DnJD1EOvLZn+pHP78m3Yf7YEQslVRq9vkzaluYdYuIh/IRy8GSHiXtZPYmteIabOeSWuedKOkE0r22r5Ma0tTrR6Sj9ksl/Yh0JtVFCnjbRsQ7qo0UEXMlHUFqpbo+6Yf1HOnUfxbp/lLVVm953B8Bn5W0MulS2jLSZYvbIuJ00s31g4A/Sfoq6ZLOXrnsj8aKlph/I91U/pWkQ0gHRV+keuvbWt0C7CLpXFJweDAiHuxh2HOBfSXdSLpE+24qAmtEhKRPkxpbXCTpaNLZ1IakZvqH9FGfE+ljWUhahRSoTmXFYwTvILWwPA9A0l9Jl66vzfO1Geme6zEMD1eTWk9+P1+qXUJqCTmmCdP+Eelg5AJJh5Hu9a5GWkcH5oPIP5Ae1ThV0g9z/4NZceZfUmqh+HFJvyGt6xtIv4XLgKMlTSSdzbyN1FDkOxFRWc5AnZTrdzSpde4lFf3PBb4k6SvAVaTECHv0VWj+vZwOfEjSHaR90i6suARaGu5JSV8Afp7nt7QvWJPU+nl2RPxW0rdIrV0vJj0isRapEdJ10dfzl321OKGH1j09tIQpf85tWq7wU6TgcVSeyfIWM58gnQXMqihv7zzce6u1tIkVrZAuq9LSJoAPV3Q/NHcvbz21Fuk5qcdJ17W/TZUWR+Rnc/qa31jRYur+GpbpJ1nxbNXVpEYas3lpq6qe5nl2RVkTSD++u0lHePNJ92Q+XWUdTqsYdx/SEeHTpIBya15Pa9UwDweSfpRLSJdUZgNbl/WfQmpA1OOzXXm4N+Zl8AzpLGlvem4tWbleqy2jbUgNk56j7+fcViNldFiUP6eS7n9Wa+21I+kHtjh/ricdkPW4ndS6LEg76GNIO7TFpEY9VwN7lg3zubyuHmPFc6OHAqP7WE+z6efvhF6ec+thOuXb78u2uWrLiDp+S6TW1LPzMroX+Cy9POdW0a1Un8pWii9bPj0sx0mkRlQPkX5n95GeSxxTNsw7SY8fPJu3j7dULpc83CGkg9nSFZxpuXvpObfSNO6g7Dm3PMyAWktWjHN1HufbVfqtRGqmv4C0Dz8LWLdy2j0s//Gs2N4XkgLoS/b9ZcO+jfS7ejIvt7tIJzwzcv9dSI2AHmLF84XHU9ZCuqdP6eFAMzOzIcNvBRgGJJ2YL6eUEjLX1PJLKSnqKYNbO4PliYYv66HfVEmLJfW3tW9v071ZvSQkV0rW/eFGT7fR6tmubXhwcBtC8o5okaRm3IewMoMZBCLi3oh4RfSSNWYAZW8U+Rm8Ig5mJM3TIL1NxIY3B7chQunVItuSrmvvVmhl6pBbvLVsedYYSlpyf9PKdbP+8wodOj5AanBwIukZqZpIWkPSHyUtkHS3pE/1MuxWki6X9HjOXLB9Wb/9JP1H0lO5nL16KONQSWfkjANPAvtJGiHpy5LmSnpM0u+V33elFe8UO0DSg5IeUspH2t/yxuZhH8vzcbWkybnfKpKOz9N4QNJhpUuBpcuGko7MZ8d3S3pr7nc46cDiqHz58KjcfQNJ50taKOl2Se8pq/eqks6U9KSkq+ilpW7ZMhiVv8+W9L+S/pmX93nq5WFupcwO1+X5vVxlD8CWzpwk7UzKX/jePA/XlxWxTk/T6mObmK2UUeWfpAZD0yvqdTLpEZj/y9P8Yg1ldis9vP5gXg9/qSjzc5Lm53XYY2vranWrYf3/U9LPJD2hlNXmzT2Vby2g1pY1/rT2h9TK6GPA60jNiyeX9TuR3HqM1Drt/vz/CFLLwm+QHrqcTkpg+1+5/6Hk1m2kJrqPkVo3jSA1ZX+MlEKqi9Taaf087BRgox7qeWiu3ztzOSuRHo34F6kFa6nl4Gl5+Gmks9HT8nQ2IbXg2qmf5X2U1Eq2k/Qs5uuAcbnfX/KwXaTWcVeRs4+TWtstJSUbGAn8D6lpcqlR1mzKWh/mMu4jPc4yipQQ99HSciG10vx9Hm5jUuu5qq32ypbBqLJpzSU9RrJS/v7dHsbdnNSCdstc731JrRZL2eznVSzLytaMPU6LXraJsnHvJbVyHEWVlp3l06+xzLOB00mthEcDbyrbrl8g5YEcncd/hipveOipbjWs/xdIrRdHkx7BeYIa36LhTwH7xKIr4E8DVmJqTr+UFclVbwM+U9b/RKoHty1JD4CXl3UwOZEsLw1uXwJOrhj273ln2UV6pGJ3yl5p1ENdDwX+UdHtVuDNZd+n5PkZxYod+wZl/b8HHN/P8vYnPUD6mopxJpOaGpe/kun95OS9eed2V1m/zlyv1fP32bw0uL0XuLRiGseQmoGXEhyUz9O3qS+4lScA/hhwbg/j/hL434put7MiKMyj7+BWdVq9bRNl436rj+1h+fRr2M6mkB4dqvZKou1JTcnLH/eZT8WrpCrm61tl32tZ/8sPZnK3q4B9Bvr79WdwPr4/MTTsC5wXKx70/G3u9qM+xlsHWENSedqvkaxIHlw57H9LentZt9GkH//Tkt5LyuB9fL7U87mIuK1KOZDOaCrL/rOk8lQ6y0g7nGrj3EM6g+tPeSeTHlb/nVI2lFNIr5FZJ8/PQ1qRiGRERdnLkyJHxDN5uJ4y6awDbFmxbEfl6U/M/1fOUz3KEzQ/00c99lXKFVjSQXqFyUCn1eM2Ufa9ct30pbcy1yblX6yWMg7gsXhpsuTelktl3WpZ/w9EjmrZPdS3HK2JHNzanNLb0N8DjJRU2gmNAcZLem1EXN/z2NxHesVJ5Vsaehr25Ij4SLWeEfF34O+5PoeRXruzbQ9lVT5ceR+wf0T8s3JApYYykHZspWA5lXQUXXd52TeBb+ayzyGdyZzDileLvNDDeL2pVodLImJW5YD5Ps4LvHyeBsN9wOERcXgNw9b70Guv20SNZVZbblXLlDSFlH9xfETUm4e1r2mXXkLa2/pfU5LKAtxUUoYea0FuUNL+3kk6K5lBelHmpqTUUJeSGpn05irgSUlfkrSSpJGSNpb0+irDngK8XdJ/5eHGKj1btJakyZJ2U3ox5xJSBol6mq0fDRwuaR0ASROVXhBa7uuSOiVtRLqPdXp/ypO0g9IrZUaS7hMuBZZFxEOkVFc/UHpd0whJr5T0phrn4RFe2mDiLODVkvaRNDp/Xi9pw0hN+v8EHJrnaQZ1NAKq06+AAyVtqaRL0i5KqdOqzcM01d5ysMdtoo76VS63HsvM6+hvwC8kTcjLtLfXM9WsxvU/CfhUnu5/k35nPb7k1Irl4Nb+9iXdI7s3UkLdhyO9U+wo0jubejw7zzvZt5MC4t2kBg/HAatUGfY+Ui69r5AadNzHijfjjiClh3qQlG7nTaR7M7X6CekI+DxJT5Eag2xZMcwlpEYzFwJHRsR5/SxvdeAMUmC7NZdberbrA6RLdreQ0nGdQbrPU+s87JFb8P00Ur7Bt5Cyqz9IurR3BCtyIX6CdMnsYdI90RNqnE5dIr1d4yOk7WERaRnu18PgpVclPSbp2hrK7m2bqNV3gK/llpGfr6HMfUgHJLeR7ql9uo5p9aWv9X8l6V2Uj5Lear5HpNd1WQty+i1rafnS4d2klnb9uVxoNmCS9iM1GHpj0XWx2vjMzczMhhwHNzMzG3J8WdLMzIYcn7mZmdmQM2yec1tttdVi2rRpRVfDzKxtXHPNNY9GxMSi69Efwya4TZs2jTlz5hRdDTOztiGp3sw5LcOXJc3MbMhxcDMzsyHHwc3MzIYcBzczMxtyHNzMzGzIadvgJukgSTdJullSI5OnmplZm2vL4CZpY1Km8y2A1wK7SqrlnWRmZjYMtGVwI71H6V8R8UzOFH8J8K7BmNCx/5jL329+uO8BzcysZbRrcLsJ2E7SqpI6gbeR3mr8EpIOkDRH0pwFCxb0a0K/ufweBzczszbTlsEtIm4lvfjxfOBc4HrgZe/6iohjI2JmRMycOLF/GWQmdI3m8WeWDqS6ZmbWZG0Z3AAi4viI2DwitiO9/fnOwZjOhM4OFj79/GAUbWZmg6Rtg5ukSfnvVODdwGmDMZ0JnR0sesbBzcysnbRz4uQ/SloVWAp8PCIWDcZEurs6WOQzNzOzttK2wS0itm3GdCZ0dvDkcy+wdNmLjB7Ztie6ZmbDivfWfZjQNRrAjUrMzNqIg1sfJnR2APi+m5lZG3Fw60N3Vw5uvu9mZtY2HNz6ML4zXZb0mZuZWftwcOtD6cxt4dO+52Zm1i4c3Prge25mZu3Hwa0PY0ePZKXRI33PzcysjTi41aC7q4OFPnMzM2sbDm41cPJkM7P24uBWAydPNjNrLw5uNXDyZDOz9uLgVoPuLp+5mZm1Ewe3Gkzo7OCpnDzZzMxan4NbDZw82cysvTi41cAPcpuZtRcHtxo4ebKZWXtxcKuBz9zMzNqLg1sNSvfcnDzZzKw9OLjVwGduZmbtxcGtBmNHj6Szw8mTzczahYNbjSZ0OnmymVm7cHCr0YSu0T5zMzNrEw5uNUr5Jd2gxMysHTi41ai7y8mTzczaRdsGN0mfkXSzpJsknSZp7GBOz6+9MTNrH20Z3CStCXwKmBkRGwMjgfcN5jSdPNnMrH20ZXDLRgErSRoFdAIPDubEup082cysbbRlcIuIB4AjgXuBh4AnIuK8yuEkHSBpjqQ5CxYsGNA0x/tBbjOzttGWwU3SBOAdwLrAGkCXpL0rh4uIYyNiZkTMnDhx4oCmWUqe7PtuZmatry2DG7ATcHdELIiIpcCfgDcM5gRLKbge95mbmVnLa9fgdi+wlaROSQLeDNw6mBNccebme25mZq2uLYNbRFwJnAFcC9xImo9jB3Oa4ztTgxLfczMza32jiq5Af0XEIcAhzZqekyebmbWPtjxzK4qTJ5uZtQcHtzo4ebKZWXtwcKtDOnNzgxIzs1bn4FaH7q4OPwpgZtYGmh7cJC2TNKlK91UlLWt2ferh5MlmZu2hiDM39dB9DNDSkcPJk83M2kPTHgWQ9Nn8bwAHSlpc1nsksC1wW7Pq0x/lyZMnrjym4NqYmVlPmvmc2yfzXwEfBsovQT4PzAMObGJ96jaha0XyZAc3M7PW1bTgFhHrAki6GHh3RCxq1rQbpZRf0vfdzMxaW9MzlETEDs2eZqM4ebKZWXsoJP2WpPeSkh1PoqJRS0TsVkSdauHkyWZm7aHpwU3S94FPAxeT3p4dza5Dfzl5splZeyjizO0DwPsj4owCpj0gpeTJvudmZtbainjObQRwXQHTbYgJnR0+czMza3FFBLdjgb0LmG5DdHd1OHmymVmLK+Ky5HhgT0mzgBuAl7TOiIhPFVCnmo3vHO3kyWZmLa6I4DaDFZclN6jo1/KNS7q7Orh34TNFV8PMzHrh59zq5OTJZmatz6+8qZOTJ5uZtb4innM7s7f+rfwQN6xInrzomeeZtPLYgmtjZmbVFHHm9ljF50lgXWA74NEC6lOXUvLkx92oxMysZRVxz+2D1bpL+gHwVJOrU7duJ082M2t5rXTP7RjgY0VXoi/jc3Dzs25mZq2rlYLb+kVXoBbdy9/p5suSZmatqogGJT+t7ARMAd4K/LrGMtYHTi/rNB34RkT8uCGV7IWTJ5uZtb4iHuLepOL7i8AC4DPUGNwi4nZgUwBJI4EHgD83sI49cvJkM7PWNxQe4n4zMDci7mlwuT1y8mQzs9ZWyMtKASSNBdYjpdyaGxHP9bOo9wGn9TCNA4ADAKZOndrP4l/OyZPNzFpb0xuUSBqdX1i6CLgeuBFYJOl7kkbXWVYHsBvwh2r9I+LYiJgZETMnTpw40Kov5+TJZmatrYgztyOA9wMHApflbtsC3yEF28/XUdZbgWsj4pGG1rAP3V0d3POYkyebmbWqIoLbnsD+EXFOWbe5khYAx1FfcHs/PVySHEy+52Zm1tqKeM5tFWBule5zSe96q4mkTmAW8KcG1atm3V1Onmxm1sqKCG7XA9VeSHoQK97z1qeIeCYiVo2IJxpWsxpN8LNuZmYtrYjLkl8Ezslv4r6C1Fpya2AN0j20lleePNlvBjAzaz1NP3OLiH+QUm39AXgFMC7/v35EXNbbuK3CyZPNzFpbIc+5RcQDwFeLmHYjOHmymVlrK+I5t09I2rtK970ltfxbAWBF8uSFvudmZtaSimhQ8mngvird55HyS7a8UvJkv7DUzKw1FRHc1gKq5YG8P/dreWNHj6TLyZPNzFpWEcHtYXJG/wqbA482uS79Nr7T+SXNzFpVEQ1Kfgv8VNLTwOzcbQfgx8CpBdSnX7q7nKXEzKxVFRHcDgHWBf4OLMvdRpAeB/h6AfXplwldHU6ebGbWoop4n9tS4P2Svg5sRnoT97URcVez6zIQEzpHM+/Rp4uuhpmZVVHY+9xyMGurgFbOyZPNzFpXEQ1KhgQnTzYza10Obv1Uyi/pszczs9bj4NZPy98M8LQblZiZtRoHt34qJU/2mZuZWespJLhJmizp85J+KWm13G0bSesWUZ/+WH5Z0g9ym5m1nCISJ78OuB3YC/gQ6ZU3kN6qfXiz69NfEzqdPNnMrFUVceZ2JPCTiNgMWFLW/e/ANgXUp1+cPNnMrHUVEdxeB/ymSveHgMlNrku/OXmymVnrKiK4PQtMqNJ9A2B+k+syIBO6nDzZzKwVFRHc/gocImlM/h6SpgFHAH8soD79NqGzw/fczMxaUBHB7fNAN7AA6AQuI6Xhehz4WgH16bcJXR0s8j03M7OWU0Ti5CeBN0rakfQOtxGkxMkXNLsuA9Xt5MlmZi2pyMTJFwEXFTX9RvALS83MWlNTgpukz9Y6bET8sMYyxwPHARsDAewfEVf0r4b9093VwVNLUvLk0SOd7MXMrFU068ztkzUOF0BNwQ34CXBuROwhqYN0/66pypMnT1p5bLMnb2ZmPWhKcIuIhqbVkjQO2A7YL5f/PND064PlyZMd3MzMWke7XkubTmpteYKkf0s6TlJX5UCSDpA0R9KcBQsWNLwSTp5sZtaaisgtKUkfk3SzpGckTc/dvyzpPTUWM4rU0vKXOY3X08CXKweKiGMjYmZEzJw4cWLD5qHEyZPNzFpTEWduB7S182IAABaWSURBVJGeZzsWUFn3B4BP1FjG/cD9EXFl/n4GKdg1VXeXkyebmbWiIoLbgcBHIuInwAtl3a8FNqqlgIh4GLhP0vq505uBWxpayxqMX37PzcHNzKyVFPGc2zrATVW6LwVWqqOcTwKn5paS/wE+2IC61WXMqJQ82VlKzMxaSxHB7T+kS4j3VHR/G3WcfUXEdcDMBtarX5w82cys9RQR3I4EjpLUSbrntrWkfYAvAvsXUJ8BcfJkM7PWU0RuyRMkjQK+TXrw+mRSY5JPRcTpza7PQDl5splZ6ykkt2RE/Ar4laTVgBER0VbvcSvn5MlmZq2nsMTJABHxaJHTbwTfczMzaz3NSpx8IylvZJ8i4jWDXJ2GmtCZkic//8KLdIxq14QvZmZDS7PO3M5o0nSarpSl5PFnnTzZzKxVNCtx8jebMZ0iLM8v6eTJZmYtw9fRBqj0ZoCFvu9mZtYyfM9tgJZflvSzbmZmLcP33AbIyZPNzFqP77kNkJMnm5m1Ht9zGyAnTzYzaz2FPMQt6YPA+4GpQEd5v4iYXkSdBsIPcpuZtZYi3sT9BeAHwDXANOAvpFfgdAO/bnZ9GqG7y8mTzcxaSRGXJT8CHBARB5Pe4XZUROxGCnjrFFCfARvf6TM3M7NWUkRwWwu4Kv//LDAu/38asHsB9Rmw7s7RvudmZtZCighuDwOr5f/vAbbO/69Hjc/CtRrfczMzay1FBLeLgN3y/8cDP5R0MXA68KcC6jNg3WXJk83MrHhFtJY8gBxUI+JoSYuAbYA/AscUUJ8BG+/kyWZmLaWIN3G/CLxY9v100llb23LyZDOz1uKHuBtgQpeTJ5uZtRIHtwaYUDpz87NuZmYtwcGtAUrJkx3czMxag4NbAzh5splZayki/da0Xvq9oY5y5km6UdJ1kuY0om79VUqevPBpP8htZtYKijhzu17S3uUdJI2Q9C3g4jrL2iEiNo2ImY2rXv9M6OrwC0vNzFpEEc+5fQk4WtLbgAOBicBvgTWBXQqoT0N0d3VwyR0L2Of4K5s+7XFjR/Pd3Tdh5bGjmz5tM7NWVMRzbkdLmg2cSnobwHjgfOCtEbGwnqKA8yQFcExEHFs5gKQDSA+NM3Xq1IFWvVe7b74Wf7nuARYveWFQp1PphWXBpXc+yo4bTGL3163V1GmbmbUqRTQ/naOk8aTUW7uSLo1+LCJ+VWcZa0TEg5ImkYLjJyPiHz0NP3PmzJgzp9Bbc4PixReDrb97IZutPYGj93ld0dUxsyFE0jWtcNunP4poULIdcAPp7QAbAfsD35f0Z0mr1lpORDyY/84H/gxsMQjVbXkjRoidNpzMP+5cwHNLlxVdHTOzllBEg5ILgJOAbSLirog4GdiUdO/txloKkNQlaeXS/8BbSJc4h6VZMybzzPPLuGLuY0VXxcysJRQR3HaKiK9FxPKbUxExD9gO+EWNZUwGLpN0PendcGdHxLkNr2mb2PqVq9LVMZLzbnmk6KqYmbWEIhqUVL0vlhMqH1ZjGf8BXtvIerWzMaNG8qb1J3LBrY9w+IsbM2KEiq6SmVmhingUAEndwM7AVKCjvF9EfKuIOrW7WTMmc86ND3P9/Y+z2dQJRVfHzKxQTQ9ukrYCzgaWkO6zPQBMyd/nAQ5u/bDD+pMYOUKcf8sjDm5mNuwVcc/t+6Rn3NYEngN2JJ3BzQGOKKA+Q8L4zg62mNbN+b7vZmZWSHB7DXBUpAfslgFjIuIRUuaSQwuoz5Axa8Zk7py/mHmPPl10VczMClVEcCtPwPgIsE7+fzGwRvOrM3TMmjEZgAtu9dmbmQ1vRQS3a4HX5/9nA4dJ2hf4Kenhbuuntbs72WD1lf1IgJkNe0UEt68CD+b/vwYsAH4GTCDngbT+mzVjMnPmLWSh3y1nZsNY04NbRMyJiIvz/wsi4q0RMS4iZkZETRlKrGezZkzmxYCLbptfdFXMzArjN3EPMZusuQqrjxvL+bc8XHRVzMwK07Tn3CSdWctwEbHbYNdlKJPETjMm8cdrHuC5pcsYO3pk0VUyM2u6Zp657QpsAjzWx8cGaNaM1Xl26TIun/to0VUxMytEMzOUHAnsTUqQfAJwYkTc38TpDxtbTe/mFWNGcf4tj7DjBpOLro6ZWdM17cwtIr4IrA18BpgJ3Cnpb5L2kDS6WfUYDlYkUp7Piy82/2W0ZmZFa2qDkohYFhFnRsQ7gXWBi0lvAnhA0iuaWZehbtaGk1nw1BKuu//xoqtiZtZ0RbaW7ALGA68gZSfxKUYDlSdSNjMbbpoa3CStJGlfSf8gvXV7HWDfiJgeEU6I2ECrdI5my3W7ucDBzcyGoaYFN0nHAg8DnwROA9aIiL0i4sJm1WG4cSJlMxuumnnm9mFgEfAQ8FbgJElnVn6aWJ8hr5RI2ZcmzWy4aeajACfh+2pNtdaETjacMo7zb3mEj2w3vejqmJk1TdOCW0Ts16xp2QqzZkzmqIvuZOHTz9Pd1VF0dczMmsK5JYe4WRumRMoX+h1vZjaMOLgNcRuvOY4pq4z1fTczG1Yc3IY4Sey04WQuvfNRnlu6rOjqmJk1hYPbMDBrxmSeXbqMf97lRMpmNjy0dXCTNFLSvyWdVXRdWtlW01dl5ZxI2cxsOGjr4AYcBNxadCVaXceoEU6kbGbDSjOfc2soSWsBuwCHA58tuDotb9aMyZx1w0Ps9vPLGD2y3Y9phpadNpzMx3dYr+hqmA0pbRvcgB8DXwRW7mkASQcABwBMnTq1SdVqTTttOJldXzOFJ55dWnRVrMx9C5/h5xffxYfeuK7fmm7WQG0Z3CTtCsyPiGskbd/TcBFxLHAswMyZM4f19biuMaM4as/Ni66GVZh9+3z2O+FqLp/7qF8sa9ZA7Xp9ahtgN0nzgN8BO0o6pdgqmdVv61euuvyt6WbWOG0Z3CLi4IhYKyKmAe8DLoqIvQuullndxowayZte7cY+Zo3WlsHNbCiZNSO9Nf16vzXdrGHaPrhFxOyI2LXoepj1l9+abtZ4bR/czNrdKp2j2WJat4ObWQM5uJm1AL813ayxHNzMWoDfmm7WWA5uZi1g7e5ONlh9ZQc3swZxcDNrEW+ZMZk59yxk4dPPF10Vs7bn4GbWImbNWJ0XAy66bX7RVTFrew5uZi1i4zXHsfq4sZx/y8NFV8Ws7Tm4mbUISew0YxL/uMNvTTcbKAc3sxYya8bqfmu6WQM4uJm1kK2mdzuRslkDOLiZtZAxo0b6relmDeDgZtZi3jJjMo8uXsJ1TqRs1m8ObmYtZvtXO5Gy2UA5uJm1mFU6R7Pluk6kbDYQDm5mLWjWjMncNX8xdzuRslm/OLiZtaAViZT9QLdZfzi4mbWgtSZ0suGUcb40adZPDm5mLWrWhpO45p5FPLZ4SdFVMWs7Dm5mLcqJlM36z8HNrEVtvOY4pqwy1pcmzfrBwc2sRUlipw0nc+mdTqRsVi8HN7MWNmvGZJ5duozL7nQiZbN6OLiZtbCtpq/qRMpm/eDgZtbCOkaN4E3rT+TC2x5hmRMpm9WsLYObpLGSrpJ0vaSbJX2z6DqZDZaUSPl5rrvPiZTNatWWwQ1YAuwYEa8FNgV2lrRVwXUyGxTbrz+JUU6kbFaXUUVXoD8iIoDF+evo/PE1GxuSVllpNFtO7+akK+Zx4a0OcFaf0SNHcM5B2xZdjaZry+AGIGkkcA2wHvDziLiyyjAHAAcATJ06tbkVNGugT+34Kk664h7Cx3BWp1Ej2vUC3cAonQS1L0njgT8Dn4yIm3oabubMmTFnzpzmVczMrM1JuiYiZhZdj/5o+5AeEY8Ds4GdC66KmZm1iLYMbpIm5jM2JK0E7ATcVmytzMysVbTrPbcpwG/yfbcRwO8j4qyC62RmZi2iLYNbRNwAbFZ0PczMrDW15WVJMzOz3ji4mZnZkOPgZmZmQ46Dm5mZDTlt/xB3rSQtAO7p5+irAUPlhVpDZV6GynyA56UVDZX5gIHNyzoRMbGRlWmWYRPcBkLSnHZ9Sr/SUJmXoTIf4HlpRUNlPmBozUs9fFnSzMyGHAc3MzMbchzcanNs0RVooKEyL0NlPsDz0oqGynzA0JqXmvmem5mZDTk+czMzsyFn2AU3STtLul3SXZK+XKX/GEmn5/5XSppW1u/g3P12Sf9Va5mDZZDmZZ6kGyVdJ6lpL8Dr77xIWlXSxZIWSzqqYpzX5Xm5S9JPJalN52N2LvO6/Jk02PMxwHmZJemavOyvkbRj2ThNXyeDOC/ttl62KKvr9ZLeVWuZbSkihs0HGAnMBaYDHcD1wIyKYT4GHJ3/fx9wev5/Rh5+DLBuLmdkLWW2y7zkfvOA1dpovXQBbwQOBI6qGOcqYGtAwN+At7bpfMwGZrbROtkMWCP/vzHwQFHrZJDnpd3WSycwKv8/BZhPSp5fyD5ssD/D7cxtC+CuiPhPRDwP/A54R8Uw7wB+k/8/A3hzPrp8B/C7iFgSEXcDd+XyaimzXealKP2el4h4OiIuA54rH1jSFGBcRFwR6dd8EvDOQZ2LQZiPAg1kXv4dEQ/m7jcDY/PZRBHrBAZhXppQ554MZF6eiYgXcvexQKnBRVH7sEE13ILbmsB9Zd/vz92qDpM3hCeAVXsZt5YyB8NgzAukDf68fAnmgEGodzUDmZfeyry/jzIbbTDmo+SEfDnp6026lNeoedkd+HdELKGYdfKSevYy3XrnpaSt1oukLSXdDNwIHJj7F7UPG1TDLbhV2/gqm4v2NEy93QfbYMwLwDYRsTnwVuDjkrbrfxVrNpB5GUiZjTYY8wGwV0RsAmybP/v0o271GvC8SNoIOAL4aB1lDobBmBdow/USEVdGxEbA64GDJY2tscy2M9yC2/3A2mXf1wIe7GkYSaOAVYCFvYxbS5mDYTDmhdIlmIiYD/yZ5lyuHMi89FbmWn2U2WiDMR9ExAP571PAb2mDdSJpLdL284GImFs2fLPXyUvq2ct0652XtlwvJRFxK/A06T5iUfuwQTXcgtvVwKskrSupg3Sz9cyKYc4E9s3/7wFclO8PnAm8L987WBd4FenmeC1ltsW8SOqStDKApC7gLcBNLT4vVUXEQ8BTkrbKl4s+APy18VV/iYbPh6RRklbL/48GdqXF14mk8cDZwMER8c/SwAWtExiEeWnT9bJuDnZIWgdYn9SArKh92OAqukVLsz/A24A7SK2Dvpq7fQvYLf8/FvgDqZHFVcD0snG/mse7nbJWXtXKbMd5IbWWuj5/bm6jeZlHOjJdTDoKnZG7zyTtcOYCR5GTFrTTfJBaUV4D3JDXyU/ILVtbdV6Ar5HOCq4r+0wqap0Mxry06XrZJ9f1OuBa4J29ldnuH2coMTOzIWe4XZY0M7NhwMHNzMyGHAc3MzMbchzczMxsyHFwMzOzIcfBzQaVpEMl1fX8j6QTJZ01WHUaqvKyfkRSSNqv6PqYFcnBzeqSA09IOq5Kv+/lfuWB6UjgTXVO5iBg74HUs9kkvVbSXyU9LOk5SfdK+mN+WLY0TEh6XtL0inFfEsxzkIr8WSbpPknHSZrYy/Q3Bg4hvVVgCnB64+fSrH04uFl/3Ae8N2cxAZan+dkHuLd8wIhYHBGP1VN4RDwREY83pKZNkIPOhaSHr3cBNiAti7nAuIrBlwGH11Ds7aQgNRX4H+DtpCz6PVkv//1LRDwcEc9WqWdHDdM1GxIc3Kw/bgDuBN5T1m0X0utaZpcPWHlZsnSWIukgSQ9IWiTpBEmdlcOUfZ8t6ZeSfiBpoaQFefwxkn4u6fF8prRP2TjT8pnPzIr6hKQ9KoZ5n6RLJD0r6d+SXiNpY0mXS3pa0mU5TVlPtgEmAB+MiGsiYl5EXBIRX4yIGyuG/RnwHkmv63UJwws5SD0QEWcBPwXeImmlygElHUrKfQjwoqSoWNZfknQ/OSO/pA5JR0i6P8/f1Sp7YW0eZmdJt+Wz0Esl7ZmX1bTcfz9JiyvG2T4Ps1pZtzfkZftMXt+/lDSurP9sSb+Q9G1Jj0qaL+lISSPKhunI/e+RtETSfyR9Ssldkj5fUY9X5Xps3scytiHMwc3663hg/7Lv+wMnUFs28W1JCVt3At4LvIt0KbI3ewFPAVsC3wV+DPyFlDJoJun9VcdJWqP2WVjum6SM75sBj5OS4P6MlKJsC1I6o5/2Mv7DpN/SHlKfrz25Cvgj8L066/hsnsaoKv2OBD6S/5+SPyVvAl4D7Ay8OXc7IXffE9iEtOz+T9JrASStTVq25wObkpZFvfVF0ibAeaQ8ha8F3p3L+3XFoHsBLwBvAD4BfJq0XZT8hpSH8rPAhsCHgMcjpVeq3A7J36+LiGvrrbMNIUXn//KnvT7AicBZpDOVZ0lJl1cHlpAuoZ0InFU2/KHATRXj30d+I3Du9ivggspplH2fDVxR9l3AAuDMsm6jgeeBPfL3aaRAO7Oi/lFlmI+W9d81d3t3Wbf9gMV9LJfDgaXAItIO/SvAOtWmTbqE+Dywcw/zW7nMNiCdKV/Zy/T3SD/nl62rBcCYsm6vBF4EplYM+xfgF/n/b5MOGlTW/2u5/tN6WibA9nmY1fL3k4DjK4bZNA9TyjX5knWbu50PHJf/f1Uefuce5nv1vNy3yt9HAg8Anyj6t+JPsR+fuVm/RMQi0qWw/UkZyGdHxL29j7XcLbHijcCQXq8xqY9xbiibdgDzSS9cLHUrBZa+yum1bOCR/PfGim5d5ZdOK0XEV0k72gPyuB8CbpH05irD3kUK6N8tv/xWYUNJiyU9C9xCOiDYq8b5KXdTvPTlmpuTDg5uyeUvzpcXdyEFPkhnR//Ky7nkin5M+3XA3hXTKWXWf2XZcDdUjFe+PWxGCsYXV5tARDxMOtgqnb3tTHox56n9qK8NIdUucZjV6tekS0aLgW/UMd7Siu9B35fIq43TWzkv5r/LLxMqvZqkr7Kjl2691jFSw5k/AH+QdDDwb+DrpMYmlb5JanDSU8CaS8rUvgx4sCJA1ePpiu8jSPPzel6+/EqNUGp5o/SLVYarXL4jgOOAH1UZ/4Gy/3tbj7XU5Tjgt5I+TQpyf8oHXzaMObjZQFxIury2GumyVitZkP+W33/atFkTj4jnJc0Fqt4DjIj5ko4E/pfqZ0XP5zO8Rvs3KWCsHhFVz4ZIZ4q7S1LZ2dtWFcMsADoljYuIJ3O3yuV7LbDRAOfjWlKg2wE4t4dhzgWeJD0G8XbSQYENc74saf2Wd3yvAdYdwJnFoIjUFP5fwJckbSTpDaSGFw0naVdJp+S/r5a0fm7B9zZWtGKs5gekxirvHIx6VRMRd5Au2Z0oaQ9J0yXNlPR5Se/Ogx1Nuh/54zwve5ACR7krSWeF35G0nqTdgY9VDHMEsIWkoyVtlofbVdIxddT3TuD3pMZCuyu9UHPb8paxEbGMdBXhO6QzwmpnyjbMOLjZgETEU2VH7q2mdB/mauAYUqOIwXAL6dLskaQzo6tID6F/ntQ4o6qIWEy6PDl2kOrVkw+SWkx+D7iNdM9qO+CeXK97SS0bdya9uPYzwJfLC4iIhaRLqrNI9xgPIF2CLR/mhlzuNOCSXNZ3WHFfs1YfILVg/Wmu74nAKhXD/BroAE6ouFdow5RfVmpmfcrPC15NOkufV3B1XkbSlqTGKtPraNhkQ5jvuZlZ25I0BlgbOAz4swOblfiypJm1s/eTUpWtSnrI2wzwZUkzMxuCfOZmZmZDjoObmZkNOQ5uZmY25Di4mZnZkOPgZmZmQ46Dm5mZDTn/Dxk3zJn/yZM1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_r,y_r)\n",
    "plt.xlabel(\"Minimum SNP frequency\",fontsize=14)\n",
    "plt.ylabel(\"Max allele count\",fontsize=14)\n",
    "plt.suptitle(\"Maximum allele count across minimum cutoff values\",fontsize=16)\n",
    "plt.title(\"Alleles represented in either tech rep\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate a mean frequency between technical replicates and filter alleles by minimum 3% of the mean\n",
    "First, use raw data and fill NA's with zeros then filter by minimum 3% average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alleles/site\t# Observed\n",
      "2\t\t1110\n",
      "3\t\t12\n"
     ]
    }
   ],
   "source": [
    "#fill the NAs of the frequency columns with zeros\n",
    "data.freqPropReplicateA.fillna(0)\n",
    "data.freqPropReplicateB.fillna(0)\n",
    "\n",
    "#generate an average frequency proportion column (simple)\n",
    "data['meanFreqProp'] = data[['freqPropReplicateA','freqPropReplicateB']].mean(axis=1)\n",
    "\n",
    "#count the sites with X alleles\n",
    "mean_3percent = data.loc[data.meanFreqProp >= 0.03]\n",
    "mean_3percent_allele_counts = get_allele_counts(mean_3percent)\n",
    "print_counts(mean_3percent_allele_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate a mean frequency between technical replicates after having dropped alleles with no representation in one technical replicate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Erik\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: UserWarning: Pandas doesn't allow columns to be created via a new attribute name - see https://pandas.pydata.org/pandas-docs/stable/indexing.html#attribute-access\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alleles/site\t# Observed\n",
      "2\t\t1090\n",
      "3\t\t12\n"
     ]
    }
   ],
   "source": [
    "#generate an average frequency proportion column (must be present in both tech reps)\n",
    "no_na_data.meanFreqProp = no_na_data[['freqPropReplicateA','freqPropReplicateB']].mean(axis=1)\n",
    "\n",
    "#count the sites with X alleles\n",
    "mean_3percent_no_na = no_na_data.loc[no_na_data.meanFreqProp >= 0.03]\n",
    "mean_3percent_no_na_allele_counts = get_allele_counts(mean_3percent_no_na)\n",
    "print_counts(mean_3percent_no_na_allele_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(*to do*)\n",
    "### Now let's count alleles by groups of plants\n",
    "Maintain separation of species and segment, but group by round, temperature, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#raw data, no 32 degree samples\n",
    "acmv_a = data.loc[data.chrom == 'ACMV DNA-A']\n",
    "acmv_b = data.loc[data.chrom == 'ACMV DNA-B']\n",
    "eacmcv_a = data.loc[data.chrom == 'EACMCV DNA-A']\n",
    "eacmcv_b = data.loc[data.chrom == 'EACMCV DNA-B']\n",
    "\n",
    "#NAs dropped\n",
    "acmv_a_no_na = no_na_data.loc[no_na_data.chrom == 'ACMV DNA-A']\n",
    "acmv_b_no_na = no_na_data.loc[no_na_data.chrom == 'ACMV DNA-B']\n",
    "eacmcv_a_no_na_data = no_na_data.loc[no_na_data.chrom == 'EACMCV DNA-A']\n",
    "eacmcv_b_no_na_data = no_na_data.loc[no_na_data.chrom == 'EACMCV DNA-B']\n"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}