(G) PacBio Illumina Comparison - Empirical Base Pair Recall Score Positions to Drop.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "%matplotlib inline\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.ticker as ticker\n",
    "from pylab import plot, show, savefig, xlim, figure, hold, ylim, legend, boxplot, setp, axes\n",
    "\n",
    "#for exporting to Adobe Illustrator\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Empirical Base Recall scores from Max Marin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "emp_basepair_recall_df = pd.read_csv('/n/data1/hms/dbmi/farhat/mm774/Projects/PacBio_Evaluation_Project/200328_PBvsI_VCeval_AnalysisDir/200328_EBR_H37rv_31CI_MM2vsPilon_V1/200328_EBR_H37rv_31CI_MM2vsPilon_V1.tsv' , sep = '\\t')\n",
    "emp_basepair_recall_df = emp_basepair_recall_df.loc[: , ['H37rv_RefPos_1based' , 'EmpiricalBasepairRecall_31CI']]\n",
    "emp_basepair_recall_df.rename(columns = {'H37rv_RefPos_1based':'pos'}, inplace = True) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pos</th>\n",
       "      <th>EmpiricalBasepairRecall_31CI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.83871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>0.83871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>0.83871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>0.83871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>0.83871</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pos  EmpiricalBasepairRecall_31CI\n",
       "0    1                       0.83871\n",
       "1    2                       0.83871\n",
       "2    3                       0.83871\n",
       "3    4                       0.83871\n",
       "4    5                       0.83871"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emp_basepair_recall_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = np.arange(0.0 , 1.01 , 0.01)\n",
    "num_positions_dropped = []\n",
    "\n",
    "for threshold_i in thresholds:\n",
    "    \n",
    "    num_positions_dropped.append( np.shape(emp_basepair_recall_df[emp_basepair_recall_df.EmpiricalBasepairRecall_31CI < threshold_i])[0] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')\n",
    "plt.rcParams['lines.linewidth']=3\n",
    "plt.rcParams['axes.facecolor']='1.0'\n",
    "plt.rcParams['xtick.color']='black'\n",
    "plt.rcParams['axes.grid']=True\n",
    "plt.rcParams['axes.edgecolor']='black'\n",
    "plt.rcParams['grid.color']= '1.0'\n",
    "plt.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnwAAAFuCAYAAAAbG8VbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XtcVHX+P/DXAbkPN0FBLiK6IoobfVHUSk1ZZM3yWobmhdKMzNxvZaWZF1qtTY3ais37FS+wWrll6SaiaWtmLQW/yAt5QYVEGdThJjbw+f3hg/N1lMEzwBkOw+v5ePBY5nPOmXkfzvaYt5/L+yMJIQSIiIiIyGbZNXcARERERKQuJnxERERENo4JHxEREZGNY8JHREREZOOY8BERERHZOCZ8RERERDaOCR8RERGRjWPCR0RERGTjmPARERER2TgmfBaSJKm5QyAiIiKyCBM+IiIiIhvHhI+IiIjIxjHhIyIiIrJxTPiIiIiIbBwTPiIiIiIbx4SPiIiIyMYx4SMiIiKycUz4iIiIiGwcEz4iIiIiG8eEj4iIiAiAXq9HTk4O9Hq9RW11naM1bZo7ACIiIqLmVFlZieTkZGRmZkIIAUmSMGDAAADAoUOHzLbV1NTAwcEBRqMRkiRBkiTExMRg1qxZcHFxac5bugMTPiIiImrVapO9+Ph4REVFISsrC0uWLIG3tzeeeeYZs22LFy9GXl4eJk6ciCFDhiArKwvp6ekAgHnz5jXzXZliwkdEREStll6vl5O9xMREAEDnzp2RkpKCoUOHYsqUKXBwcLijzWAwoLS0FAkJCYiLi0PPnj0RHR0NAEhPT8f06dPh4+PTnLdmgnP4iIiIWjBJau4IWraCggIIIRAVFWXS5ubmhu7du6OqqqrOttrr+vTpAwDyeVFRURBCoLCw0Po3Uw8mfERERNRqBQYGQpIkZGVlmbSVl5fj2LFjcHJyqrOt9rqjR48CgHxeVlYWJElCQECA9W+mHhzSJSIiolbLx8cHMTEx8ty72vl65eXl2LNnDzp27Gi2zd3dHRs3bkR1dTWMRqM8hy8mJkZTw7kAIAkhRHMH0ZJIkgT+yYiISCskCeDXUuO0hlW6TPgsxISPiIi0hAlf09Hr9SgsLERAQIDcQ6ekra5ztIYJn4WY8BERkZYw4SMluGiDiIiIyMYx4SMiIiKycUz4iIiIiGwcEz4iIiIiG2e1hG/QoEFwdnaGTqeDTqdDt27d5GNbt25FSEgI3NzcMGrUKJSUlMjHSkpKMHr0aLi5uSEkJARbt241eV+1riUiIiKyFVbt4UtJSUFZWRnKyspw4sQJAEBubi4SExORmpqKoqIiuLq64rnnnpOvmTFjBhwdHVFUVIQtW7Zg+vTpyM3NVf1aIiIiIlthtbIsgwYNwsSJE/H000+btM+dOxdnz56Ve99OnTqF7t27Q6/Xw87ODt7e3vj5558RFhYGAJg0aRICAwPx9ttvq3atu7u72ftgWRYiItISlmUhJazaw/faa6/B19cXDzzwAA4cOADgZk9bZGSkfE6XLl3g6OiIkydP4uTJk7C3t5cTNgCIjIw06aVT41oiIiIiW2K1vXSXLFmCHj16wNHREWlpaRg+fDh++uknlJWVwdPT0+RcT09PlJaWwt7e3uwxAKpde7uMjAxkZGQ07MaJiIiImpnVEr6+ffvKvyckJGDbtm348ssvodPpYDAYTM41GAxwd3eHnZ2d2WMAVLv2drGxsYiNjQVwM3ElIiIiakmarSxL7Vy4iIgIZGdny+2nT59GVVUVwsLCEBYWBqPRiLy8PPl4dnY2IiIiAEC1a4mIiIhsiVUWbVy9ehXfffcdHnzwQbRp0wbp6el45plnkJWVBaPRiPvuuw9ffPEFoqKikJiYCKPRiLS0NADAuHHjIEkS1qxZg59++gnDhg3D4cOHERERgdzcXNWuNfsH46INIiLSEC7aIEWEFVy6dEn07t1b6HQ64enpKfr27Su++uor+fiWLVtEcHCwcHV1FSNGjBB6vV4+ptfrxciRI4Wrq6sIDg4WW7ZsMXlvta41x0p/MiIiIkX4tURKWK0si61gDx8REWkJe/hICW6tRkRERGTjmPARERER2TgmfEREREQ2jgkfERERkY1jwkdERERk45jwEREREdk4JnxERERENo4JHxEREZGNY8JHREREZOOY8BERERHZOCZ8RERERDaOCR8RERGRjWPCR0RERGTjmPARERER2bg25g7U1NQoegM7O+aMRERERFpmNuFr06YNJEm66xtUV1c3aUBERERE1LTMJnxnzpyRf//iiy+wY8cOvPbaawgJCUF+fj6WLFmCRx991CpBEhEREVHDSUIIcbeT/vCHP+CHH36Al5eX3HblyhX07t0bp06dUjVArZEkCQr+ZERERFYhSQC/luhuFE3Au3btGioqKkzaKioqcO3aNVWCIiIiIqKmY3ZI91YJCQmIjY3FCy+8gODgYJw/fx4ffPABEhIS1I6PiIiIiBpJ0ZBuTU0NVq1ahe3bt6OwsBAdOnTA448/jmnTpsHe3t4acWoGh3SJiEhLOKRLSihK+Oj/MOEjIiItYcJHSiiawyeEwOrVq/GnP/0J99xzDwDg4MGD+Oc//6lqcERERETUeIoSvgULFmDt2rWYNm0azp07BwAICgrCkiVLVA2OiIiIiBpP0ZBucHAwfvzxR/j6+sLb2xtXrlyBEAJt27bFlStXrBGnZnBIl4iItIRDuqSEoh6+6upq6HQ6AJB33ygrK5PbiIiIiEi7FCV8w4YNw0svvYSqqioAN+f0zZ8/H8OHD1c1OCIiIiJqPEUJ37vvvovCwkJ4enri2rVr0Ol08vZqRERERKRtFpVluXTpEvLz8xEcHAx/f38149IszuEjIiIt4Rw+UkLRThsAcPXqVezduxeFhYUICAjAsGHD4O3trWZsRERERNQEFPXwZWZmYsyYMejWrRtCQkJw7tw5HD9+HB9//DH+9Kc/WSNOzWAPHxERaQl7+EgJRQlfjx49kJSUhMcff1xu2759O+bPn4/jx4+rGqDWMOEjIiItYcJHSihK+Ly8vKDX6032zTUajfD19cXVq1dVDVBrmPAREZGWMOEjJRSt0p08eTL+8Y9/mLQtX74ckydPViUoIiIiImo6inr4+vfvj++++w5+fn4IDAxEQUEBLl26hL59+8qFmIGb++vaOvbwERGRlrCHj5RQlPBt3LhR0ZslJCQ0OiCtY8JHRERawoSPlLCoDh8x4SMiIm1hwkdKKJrDBwDr1q3DkCFDEBERgSFDhmDt2rVMfIiIiIhaAEWFl1999VX861//wgsvvCDX4XvnnXdw4sQJLF26VO0YiYiIiKgRFA3ptm/fHllZWQgKCpLbzp8/j6ioKFy+fFnVALWGQ7pERKQlHNIlJRQN6bq7u8Pd3f2ONg8PD1WCIiIiIqKmo2hI94UXXsCYMWMwZ84cBAUF4fz581i2bBlefPFFnD59Wj6vc+fOqgVKRERERA2jaEjXzu7uHYGSJKG6urpJgtIyDukSEZGWcEiXlFDUw1dTU6N2HERERESkEsVlWQDg3Llz+Pbbb3H+/PkGf2BeXh6cnZ0xceJEuW3r1q0ICQmBm5sbRo0ahZKSEvlYSUkJRo8eDTc3N4SEhGDr1q0m76fWtURERES2QlHC99tvv+HBBx/EH/7wB4wZMwZdunTBwIEDUVhYaPEHzpgxA9HR0fLr3NxcJCYmIjU1FUVFRXB1dcVzzz1ncr6joyOKioqwZcsWTJ8+Hbm5uapfS0RERGQrFM3hGzVqFDp27Ii//e1vcHNzQ3l5OebOnYszZ87gs88+U/xhaWlp+OSTT9CjRw/8+uuv2Lx5M+bOnYuzZ8/KvW+nTp1C9+7dodfrYWdnB29vb/z8888ICwsDAEyaNAmBgYF4++23Vbv29hXJJn8wzuEjIiIN4Rw+UkLRHL5vvvkGv/32GxwcHAAAbm5uWLp0KQIDAxV/kMFgwIIFC7Bv3z6sXbtWbs/NzcX9998vv+7SpQscHR1x8uRJ2NnZwd7eXk7YACAyMhJff/21qtf26tVL8X0RERERaZ2ihM/b2xu//PILIiMj5bYTJ07Ay8tL8QfNnz8fU6dORXBwsEl7WVkZPD09Tdo8PT1RWloKe3t7s8fUvPZ2GRkZyMjIUHyvRERERFqieGu12NhYTJ06FSEhIcjPz8f69euxaNEiRR/y008/ISMjAz/++OMdx3Q6HQwGg0mbwWCAu7s77OzszB5T89rbxcbGIjY2FgCwZMkSJbdMREREpBmKEr5p06ahS5cu2Lp1K3JychAQEIBt27YhJiZG0YccOHAAZ8+eRceOHQHc7F2rrq7GL7/8gqFDhyI7O1s+9/Tp06iqqkJYWBjs7OxgNBqRl5eHrl27AgCys7MREREBAIiIiFDlWiIiIiJbctdFG9XV1ZgyZQpWrVoFJyenBn1IRUWFSW/aO++8g7Nnz2L58uW4dOkS7rvvPnzxxReIiopCYmIijEYj0tLSAADjxo2DJElYs2YNfvrpJwwbNgyHDx9GREQEcnNzVbvW7B+MizaIiEhDuGiDlLhrWRZ7e3t89dVXinbbMMfV1RX+/v7yj06ng7OzM9q1a4eIiAisWLECEyZMQPv27VFaWoqPPvpIvvajjz5CZWUl2rdvj/Hjx2P58uUmvXRqXUtERERkKxSVZVm6dCmuXr2KN954Q16p21qxh4+IiLSEPXykhKKELzg4GBcvXoS9vT3atWsHSZLkY+fOnVM1QK1hwkdERFrChI+UULRoY/PmzWrHQUREREQqUdTDR/+HPXxERKQl7OEjJe7aw1dUVITk5GQcOnQIJSUlaNu2LQYOHIgXX3wR/v7+1oiRiIiIiBqh3h6+ixcvolevXmjXrh1GjhyJgIAAFBQU4PPPP0dRURH++9//okOHDtaMt9mxh4+IiLSEPXykRL0J38yZM3Hx4kWkp6eblGURQmDcuHFo164dUlJSrBKoVjDhIyIiLWHCR0rUm/CFh4fj008/Rffu3e84duzYMYwaNQonTpxQNUCtYcJHRERawoSPlKg34fP09ERJSQns7e3vOGY0GuHt7Y3S0lJVA9QaJnxERKQlTPhICUU7bdSlTZs2JvX4iIiIiEib6l2lW1lZicmTJ9d5TAiBqqoqVYIiIiIioqZTb8L3+uuv13vx3LlzmzQYIiIiImp6LLxsIc7hIyIiLeEcPlLirnP4iIiIiKhlY8JHREREZOOY8BERERHZOCZ8RERERDZOUcK3bds2HDt2DABw4sQJDBw4EDExMTh+/LiqwRERERFR4ylapdulSxccPnwYfn5+GD58OLp16wadToeDBw8iMzPTGnFqBlfpEhGRlnCVLimhKOHz8PCAwWDA9evX0aFDB1y8eBEODg7w9fVFSUmJNeLUDCZ8RESkJUz4SIl6Cy/XateuHX799Vf8v//3/xAdHQ0nJydUVFQw8SEiIiJqARQlfPPnz0evXr1gb2+P9PR0AMC+ffsQGRmpanBERERE1HiKd9qoqKgAALi6ugIALl26hJqaGvj7+6sXnQZxSJeIiLSEQ7qkhKIePuBmonfp0iVcvHhRzXiIiIiIqIkpSvj27NmDqVOn4rfffjNplyQJ1dXVqgRGRERERE1DUR2+GTNmYP78+SgvL0dNTY38w2SPiIiISPsUzeFr27Yt9Ho9JEmyRkyaxjl8RESkJZzDR0oo6uGbOnUq1q9fr3YsRERERKQCRT18AwYMwNGjRxESEnLHqtyDBw+qFpwWsYePiIi0hD18pISihG/jxo1mjyUkJDRpQFrHhI+IiLSECR8pobgOH93EhI+IiLSECR8poWgOHwCsX78eMTEx6NatG2JiYjinj4iIiKiFUFSH780338SmTZswa9YshISEID8/H0uXLkVhYSFef/11tWMkIiIiokZQNKQbGhqKAwcOICQkRG7Lz8/HwIEDkZ+fr2qAWsMhXSIi0hIO6ZISioZ0y8vL0a5dO5M2Hx8fVFZWqhIUERERETUdRQnf0KFDMWHCBJw4cQKVlZU4fvw4EhIS8Oc//1nt+IiIiIiokRQlfCkpKXB3d0dkZCR0Oh3uvfdeuLm54cMPP1Q7PiIiIiJqJIvKstTU1KC4uBi+vr6ws1O8wNemcA4fERFpCefwkRJmE76zZ8+iU6dOAIDTp0+bfYPOnTurEphWMeEjIiItYcJHSphN+Nzd3VFaWgoAsLOzqzPRkSQJ1dXV6kepIUz4iIhIS5jwkRLcacNCTPiIiEhLmPCREoom4v3lL3+ps/2FF15o0mCIiIiIqOkp6uHz8PCAwWC4o93Hxwd6vV6VwLSKPXxERKQl7OEjJerdWm3dunUAAKPRKP9e6/Tp0/D19VUvMiIiIiJqEvUmfKmpqQCAGzduyL8DN3u5/Pz8sHHjRnWjIyIiIqJGq3cO3/79+7F//37MmTNH/n3//v3IzMzEtm3b0K9fP8UfNHHiRHTo0AEeHh4ICwvDmjVr5GP79u1DeHg4XF1dMXjwYJP9eauqqjBlyhR4eHjA398f7777rsn7qnUtERERkc0QZtTU1Mi/V1dXm/1R6ueffxbXr18XQghx7Ngx4efnJ3744Qdx+fJl4eHhIf75z3+KyspK8fLLL4u+ffvK182ZM0f0799flJSUiF9++UX4+fmJ3bt3CyGEqteaU8+fjIiIyOr4tURKmF20cetCjdo6fLclig2uw3fixAkMGjQI77//Pq5evYoNGzbg8OHDAIDy8nL4+vrixx9/RHh4OAIDA7F+/XrExcUBAObPn4+8vDykpaVh1apVql1rDhdtEBGRlnDRBilhdkg3NzdX/v3MmTM4ffq0yU9tmyWee+45uLq6Ijw8HB06dMCwYcOQm5uLyMhI+Rw3Nzd06dIFubm5uHLlCgoLC02OR0ZGyrGpdS0RERGRLTG7aCM4OFj+PSQkxORYZWUl7O3t4ejoaNGHffTRR/jwww/x7bff4sCBA3ByckJZWRnatWtncp6npydKS0tRVlYmv779GADVrr1dRkYGMjIyLLpXIiIisj69Xo+CggIEBgbCx8fHojZbVu8q3Vovv/wyHn/8cfTp0wdffPEFHnvsMUiShPT0dAwfPtyiD7S3t0f//v2xefNmLF++HDqd7o4afwaDAe7u7tDpdPJrZ2dnk2MAVLv2drGxsYiNjQUALFmyxKL7JSIiIvVVVlYiOTkZmZmZ8rSzAQMGAAAOHTpUb1tMTAxmzZoFFxeX5rwFVSnaaWPLli3o2bMnAOCvf/0rNm/ejM8++wxz585t8AcbjUacOnUKERERyM7OltvLy8vldm9vb3To0MHkeHZ2NiIiIgBAtWuJiIioZalN9uLj47F06VLEx8cjNTUVu3btumtbZmYmkpOTm/sW1KVkZYeHh4cQQoji4mLh6+srt7u7uytaGVJUVCS2bdsmSktLhdFoFHv27BGurq5i586d4tKlS8LDw0Ps2LFDVFZWildffdVktezs2bPFwIEDRUlJiTh27Jjw9/eXV9qqea05Cv9kREREVsGvpZv5yeDBg8WKFStM2nr27Cnmzp0rbty4YbZNCCFWrFghBg8eLIqLi60eu7Uo6uELCwvDli1bkJKSgiFDhgAAiouLFXd9SpKE5cuXIygoCN7e3nj55Zfx97//HSNHjkS7du3w8ccf4/XXX4e3tze+++47pKWlyde+8cYb6NKlC0JCQvDggw/ilVdewdChQwFA1WuJiIioZSgoKIAQAlFRUSZtbm5u6N69O6qqqsy2AUBUVBSEECgsLLR67NaiaC/d77//Hv/7v/8LBwcHrFu3Dl26dMGWLVuwZ88ekx04WgOWZSEiIi1hWZabCzDGjh2L+Ph4JCYmym2DBg3CiBEjkJSUBAcHhzrbAGDlypVIT0/H9u3bbXYBh6JFG9HR0XK9uloTJkzAhAkTVAmKiIiISCkfHx/ExMQgPT0dwM0eu6ysLJSXl2PPnj3o2LFjvW3p6emIiYmx2WQPUNjDB9zcZi01NVVewjxx4kTExMSoHZ/msIePiIi0hD18N3GVbv0UJXxr1qzB3Llz8fTTTyMkJATnzp3D2rVrsWjRIkybNs0acWoGEz4iItISJnym9Ho9CgsLERAQYFJzT0mbLVOU8IWFhWH79u0mO1Pk5OTg0UcfRV5enqoBag0TPiIi0hImfKSEooTPx8cHFy9elCc3AkBVVRUCAgKg1+tVDVBrmPAREZGWMOEjJRSVZenfvz9eeuklVFRUALhZpPiVV17B/fffr2pwRERERNR4ihK+FStWICcnB56envDz84OXlxeys7OxcuVKteMjIiIiokZSvEoXAC5cuCBPcAwKClIzLs3ikC4REWkJh3RJiXoTvoqKCixevBg///wzoqKi8Nprr8HJycma8WkOEz4iItISJnykRL1Dus8//zw+//xzhIeHY8eOHXj55ZetFRcRERERNZF6e/g6dOiArKwsdOjQAefPn8fAgQNx5swZa8anOezhIyIiLWEPHylRb8Ln4eEBg8Egv27bti1KSkqsEphWMeEjIiItYcJHStS7l67RaMT+/fvlBOf21wBa5fZqRERERC1JvT18nTp1giRJ5i+WJJw+fVqVwLSKPXxERKQl7OEjJSwqy0JM+IiISFuY8JESigovExEREVHLxYSPiIiIyMYx4SMiIiKycUz4iIiIiGycooTvhRdewPfff692LERERESkAkUJnxACI0eORNeuXbFw4UKcOHFC7biIiIiIqIkoLstSU1ODffv2Ydu2bfj000/RuXNnTJgwAS+99JLaMWoKy7IQEZGWsCwLKdGgOnwFBQV46qmnsG/fPlRXV6sRl2Yx4SMiIi1hwkdKKF60UVZWhs2bN+Phhx9GWFgY2rRpg40bN6oZGxERERE1AUU9fGPHjsXu3bsRFRWF8ePHY+zYsfD19bVGfJrDHj4iItIS9vCREm2UnNS7d28kJyejY8eOasdDRERERE2swXvp3rhxA6tXr8aMGTOaOiZNYw8fERFpCXv4SIm7zuHbt28fkpOT8a9//QsAYDQa8cEHHyA0NBQrVqxQPUAiIiIiapx6h3SXLFmCRYsWISIiArm5uXjuuedw4MABODk5YdWqVXj44YetFScRERERNVC9Q7qdO3fG9u3b0atXLxw5cgQPPPAA3nnnHbz44ovWjFFTOKRLRERawiFdUqLehM/DwwMGg0F+7erqivLyckiSZJXgtIgJHxERaQkTPlLirqt0hRDyj7OzM4Cbu27UsrNTXMqPiIiIiJpBvT18dnZ2Jr15Qgj5de3v3GmDiIio+bCHj5Sot4fvzJkz1oqDiIiIiFTS4Dp8rRV7+IiISEvYw0dK1DsBb+TIkdiyZQsqKyutFQ8RERERNbF6e/gcHBzQsWNHXL58GY899hieeuopDBgwwJrxaQ57+IiISEvYw0dK1NvD5+LiglOnTuHzzz+HJEl45JFH0LlzZyQlJXF+HxEREVELYVEdvsrKSnz88cfYtGkT9u/fjwceeAAHDhywRpyawR4+IiLSEvbwkRIWJXy3unDhAjZv3ow5c+aoFpwWMeEjIiItYcJHStSb8PXs2RM///yzNePRPCZ8RESkJUz4SAmWZbEQEz4iItISW0z49Ho9CgoKEBgYCB8fH4vaqG533VrtdjU1NTAYDPDy8lIjHiIiImqlKisrkZycjMzMTHlHr9rqIIcOHaq3LSYmBrNmzYKLi0tz3oJm1ZvwFRYWYvr06Th16hRmz56NsLAwjBw5EpcvX0ZUVBQ++eQTBAcHWytWIiIismG1yV58fDyioqKQlZWFJUuWwNvbG88880y9benp6QCAefPmNfNdaFO9ZVlmzpyJ8PBwDB8+HE8//TQOHjyI7Oxs5OfnIzQ0FLNnz1b0IVVVVZg6dSpCQkLg7u6O//mf/8Hu3bvl4/v27UN4eDhcXV0xePBg5Ofnm1w7ZcoUeHh4wN/fH++++67Je6t1LREREVmPXq+Xk73ExERER0fjscceg5ubG4YOHYopU6aYbUtMTER8fDwyMzOh1+ub+1Y0qd6E7+DBg1i8eDGSkpJgNBoxbdo0+Pn5ISgoCP/4xz8Ul2QxGo0IDg7G119/jWvXrmHRokV4/PHHcfbsWRQXF2PMmDFYtGgRSkpK0Lt3b8THx8vXJiUlIS8vD/n5+di/fz+WLl2KPXv2AICq1xIREZH1FBQUQAiBqKgokzY3Nzd0794dVVVVZtsAICoqCkIIFBYWWj32lqDehM9oNMLBwQFOTk7w8PAwmbfn6+uL0tJSRR/i5uaGpKQkdOrUCXZ2dnjkkUcQGhqK//73v/jkk08QERGBsWPHwtnZGUlJScjOzsbx48cBAJs2bcL8+fPh7e2N7t27Y9q0adiwYQMAqHotERERWU9gYCAkSUJWVpZJW3l5OY4dOwYnJyezbQCQlZUFSZIQEBBg9dhbgnrn8AUHB+PixYvw9/dHTk6OybG8vDy0a9euQR9aVFSEkydPIiIiAsuXL0dkZKR8zM3NDV26dEFubi78/PxQWFhocjwyMhI7d+4EAOTm5qpybXh4eIPui4iIiBrGx8cHMTEx8ly82rl55eXl2LNnDzp27FhvW3p6OmJiYrha14x6E74PPvgA9vb2AHDH4ozc3Fy89NJLFn/g77//jgkTJiAhIQHh4eEoKyu7I3H09PREaWkpysrK5Ne3HwOg2rW3y8jIQEZGhsX3SkRERMrNmjULAJCeno60tDRIkoRJkyYpaqtdpUt1qzfhGzRokNljo0ePtvjDampqMGnSJDg6OiIlJQUAoNPp7tjNw2AwwN3dHTqdTn7t7OxsckzNa28XGxuL2NhYAMCSJUssvm8iIiK6OxcXF8ybNw/Tp09HYWEhAgICTGruKWmjutWb8GVmZt71DWJiYhR9kBACU6dORVFREb788ks4ODgAACIiIrBx40b5vPLycpw6dQoRERHw9vZGhw4dkJ2djSFDhgAAsrOzERERoeq1RERE1Hx8fHzuSOCUtpEZoh6dOnUy+bG3tzd5HRoaWt/lJhITE0Xfvn1FaWmpSfulS5eEh4eH2LFjh6isrBSvvvqq6Nu3r3x89uzZYuDAgaKkpEQcO3b2+QClAAAgAElEQVRM+Pv7i927d6t+rTl3+ZMRERFZFb+WSAmL/m/i5eXVoA85e/asACCcnJyEm5ub/LN582YhhBB79+4V3bp1E87OzuLBBx8UZ86cka+9fv26eOqpp4S7u7to3769SE5ONnlvta41hwkfERFpCb+WSAmL9tJt27YtSkpKVOttbAm4ly4REWmJLe6lS02v3jp8RERERNTyMeEjIiIisnH1rtIdMGAAJEmSX5eWlmLgwIEm5xw8eFCdyIiIiIioSdQ7h+/WsiXmJCQkNGlAWsc5fEREpCWcw0dKWLRog5jwERGRtjDhIyU4h4+IiIjIxjHhIyIiIlXo9Xrk5ORAr9fX+dpcGzW9ehdtEBEREVmqsrISycnJyMzMhBACNTU1cHBwgNFohCRJkCQJAwYMAAAcOnQIQghIkoSYmBjMmjULLi4uzXwHtsdsD1+/fv3k39944w2rBENEREQtX22yFx8fj6VLl8LLywuFhYWIi4vD0qVLER8fj9TUVOzatUs+Jz4+HpmZmUhOTm7u8G2S2UUbbdu2RWFhIZydneHh4QGDwWDt2DSJizaIiEhLtLZoQ6/XY+zYsYiPj0diYqL8+s9//jPi4uLQs2dPGAwGDBo0CCNGjEBSUhIcHBwAACtXrkR6ejq2b98OHx+fZr4T22J2SHfkyJEICwtDp06dUFlZeUf9vVqsw0dERES1CgoKIIRAVFSUyes+ffoAAKqqqlBQUAA3Nzd0794dVVVVcsIXFRWFtLQ0FBYWMuFrYmYTvvXr1+Obb77B2bNn8f3332Pq1KnWjIuIiIhaoMDAQEiShKysLERHR8uvjx49iri4ODg5OSEwMBDl5eU4duwY4uPj5WuzsrIgSRICAgKa8Q5sU72LNvr374/+/fvjxo0bra7AMhEREVnOx8cHMTExSE9PB3Cz187d3R0bN25EdXU1jEYjsrKyUF5ejj179qBjx46IiopCVlYW0tPTERMTw949FSguvLx//36kpqaioKAAgYGBmDhxImJiYtSOT3M4h4+IiLREa3P4AK7S1SJFCd+aNWswd+5cPP300wgJCcG5c+ewdu1aLFq0CNOmTbNGnJrBhI+IiLSkKRM+vV4vd+zU9rI1tq2wsBABAQHw8fG543Vd55A6FCV8YWFh2L59OyIjI+W2nJwcPProo8jLy1M1QK1hwkdERFrSFAnf7T1y5nrglLaxp057FBVe1uv16NGjh0lbt27dUFJSokpQREREZD231s2rnU+3ZMkSeHt745lnnrG4rXb+3rx585r5zqiWoq3V+vfvj5deegkVFRUAgPLycrzyyiu4//77VQ2OiIiI1KXX6+VkLzExEdHR0Xjsscfg5uaGoUOHYsqUKRa1JSYmykWUuV2adihK+FasWIGcnBx4enrCz88PXl5eyM7OxsqVK9WOj4iIiFR0e9282rZb6+RZ0gbcXJkrhEBhYaF1b4bMUpTwdejQAV9//TXOnDmDzz//HGfOnMHXX3/NOjlEREQt3K11825tq62T5+TkZFEbwHp6WqRoDl+toKAgBAUFqRULERERWVlddfPM1clT2sZ6etqjuA4f3cRVukREpCW3r9JtSBkVV1dXrtK1cUz4LMSEj4iItKQ24WtMaZXaBK2iokJRnTzW02t57prw1dTU4MCBA+jfvz8cHR2tFZdmMeEjIiItqU34Fi9e3ODSKrVDsCyjYrsU9fC5u7ujtLTUGvFoHhM+IiLSEkkCiov1GDt2rFxaBbjZ4zZo0CCMGDECSUlJcHBwqLMNAFauXIn09HRs376dvXM2StEq3YEDB+LIkSNqx0JEREQN0JjSKgDLqLQGilbphoSE4KGHHsLIkSMRHBwMSZLkY3/9619VC46IiIju7tbSKtHR0XJbbcmU+Ph4s20Ay6i0BooSvsrKSowaNQoAcOHCBVUDIiIiIss0trQKy6jYPq7StRDn8BERkdosKa3i6+uD4mJ9o0qrsIyK7VOc8B07dgw7duxAUVERUlJScOLECVRVVeGee+5RO0ZNYcJHRERqaUhplf37MzF4cEyTlFYh26VoSHf79u147rnn8Oijj2Lr1q1ISUlBaWkp5syZg4yMDLVjJCIiahVqkz1LSqvs3w/Ex8fLw7nz5s27I4Hz8fFR1Ea2S1EPX/fu3bFt2zbce++98Pb2xpUrV/D7778jICAAly9ftkacmsEePiIiUoNe37DSKo6ODhCCpVWoforKsly6dAmRkZEAIK/QlSTJZLUuERERNRxLq5CaFCV8vXr1QmpqqklbWloa+vTpo0pQREREtkCv1yMnJwd6vf6ubVevXkV1dTWysrLk9lvLqDg5OZltA1haheqnaA7fBx98gLi4OKxduxbl5eX485//jJMnT+Krr75SOz4iIqIWp6H72l64cAGrVq1CTU0Nevfurai0CpAoD+eytAqZo3iVbkVFBXbt2oX8/HwEBwfjkUcegU6nUzs+zeEcPiIiupuG7mt75MgRvPfee3BycoK/v3+DVumytArVxaI6fAUFBfIS7sDAQDXj0iwmfEREtsuS+nfm2lxcXJCYmNiofW03b96Mt956Cz169LhraZXaOnzs2aP6KBrSPXfuHCZMmIAjR47Iq3T79u2LLVu2ICQkRO0YiYiIVNXQIdi62srLy1FZWSkvdgTuXGjh4OBQZxtwc/FFWloavLy8TJI4c6VVbv1fInMULdpISEhAr169cPXqVVy6dAlXrlxBdHQ0EhIS1I6PiIhIdbfWv1u6dCni4+ORmpqKXbt2NaitpKQEmZmZ8vtz8QU1N0U9fP/973/x1Vdfyf/60Ol0WLJkCf9FQUREmtOQIdjaZK92CLZz585ISUnB0KFDMWXKFDg4OChui46OxsGDB7F582Z4eXkhOjqa+9pSs1OU8PXr1w9Hjx7FAw88ILf98MMPuO+++1QLjIiIyBINHZZtzBCsuWHZl19+GZMnT0Zqaiq2b98OSZIwadIkAEB6ejrS0tLMttUuviBqSmYTvgULFsi/d+nSBcOGDcPDDz+M4OBgnD9/Hl9++SWeeOIJqwRJRER0Nw3ZliwqKgoHDx7Eu+++i8zMTPTr1w+A6XBrfHy8RW0A8MsvvyA0NBQrV67E9evXua8tNTuzCd/58+dNXo8ZMwbAzV03nJycMHr0aFy/fl3d6IiIqNVSMjTbFMOyjR2CrW9YtmvXrnfcF/e1peZgUVkWYlkWIiK1KRmarampgYODA4xGIyRJkodlV69eLffS5eTk4JlnnsHzzz+PUaNGQafT1dlW+76TJ0+Gv78/nJ2dG7VK99ZhWWvUxJMkgF9LdFdCofLycpGdnS3+85//mPwo9eGHH4pevXoJR0dHkZCQYHIsIyNDdOvWTbi4uIhBgwaJs2fPyseuX78unnrqKeHu7i78/PxEcnKyVa41x4I/GRFRq1dcXCyys7NFcXGx4rZFixaJwYMHixUrVoijR4+KFStWiNDQUBEVFSW3jRgxQnTv3l28+eab4ujRo+Kdd94RAQEB4s033zR5z549e4q5c+eKGzdumG0TQogVK1aIwYMHi5MnT4qcnJw7YmtomzXwa4mUULRoY9OmTXj++efh6Oho8q8VSZJw7tw5RYllQEAA5s2bh3//+9+orKyU24uLizFmzBisWbMGw4cPx/z58xEfH48jR44AAJKSkpCXl4f8/HxcvHgRgwcPRo8ePTB06FBVryUiooZr6AIKo9GI4uJi/OUvfzE7NGswGFBaWoqEhATExcWhZ8+eTbYytjFDsByWJS1TNKTr7++P1NRUDBkypNEfOG/ePFy4cAEbNmwAAKxatQobNmzA4cOHAQDl5eXw9fXFjz/+iPDwcAQGBmL9+vWIi4sDAMyfPx95eXlIS0tT9VpzOKRLRHR3Dd1a7LPPPsOWLVvw3nvvYeTIkQDuHJo9ffo0/vd//xcLFiyAl5cXunbtCp1O16hh2Za8LRmHdEkJRT18jo6OGDRokCoB5ObmmiyFd3NzQ5cuXZCbmws/Pz8UFhaaHI+MjMTOnTtVvba+hI+IqLWxdl27nTt34ujRoxg2bBgcHBzuWAkbGBgISZJw9OhRxMXFyUWLuTKWyDxFCd+iRYvw0ksvYeHChfD19W3SAMrKytCuXTuTNk9PT5SWlqKsrEx+ffsxNa+9XUZGBjIyMhp6i0RELVJz1LXz8fHBAw88gE8//RR+fn6477776hyGdXd3x8aNG1FdXQ2j0ciVsUR3oSjhCwsLw4IFC/DRRx/JbbX/oVdXVzcqAJ1OB4PBYNJmMBjg7u4ur54yGAxwdnY2OabmtbeLjY1FbGwsAGDJkiWNul8iImswV75EaS9dYGAgli9fbvW6dgAQERGBr776Cjt27MCnn35aZ4HimpoaBAQE4KuvvsLevXtZsJjoLhQlfJMmTcLkyZMRHx/f5PMbIiIisHHjRvl1eXk5Tp06hYiICHh7e6NDhw7Izs6W5w9mZ2cjIiJC1WuJiFqq23vlbi9f0pSLJ9Sqa/fpp59iypQpmD59+l2HYTksS6SMokUb3t7eKCkpgSRJDf4go9EIo9GIN954AxcuXMDq1avRpk0bXLlyBX/4wx+wbt06PPzww1i4cCG+/vprebXsnDlz8O2332Lnzp0oKirC4MGDsX79egwdOhSXL19W7VqzfzAu2iAiK2mKXrnFixcjLy8PEydOxJAhQ5p08YRade1a8gKK5sBFG6SIktotL774oti4cWOj6r8sXLhQADD5WbhwoRBCiL1794pu3boJZ2dn8eCDD4ozZ87I191aS699+/Z31NJT61pzFP7JiIgarKKiQq5FN2jQIDF48GCxYMECsWDBgnrb+vfvL8LDw8VHH30khLhZF27w4MHi7bffFllZWeLGjRuKa9M1pk2IpqlrR8rwa4mUUDSke/ToUaSkpODNN9+En5+fybGDBw8qSiyTkpKQlJRU57HY2FgcP368zmNOTk5Yt24d1q1bZ9VriYiagjXnztX2ygUEBAC4uVBCCIE+ffoAAKqqqpp88QS3FiNqGRQlfNOmTcO0adPUjoWISJMakrS5urqqVnhYaUmTusqXNPXiCXNtXEBBpC2KEr6EhAS14yAi0pyGliWRJAl2dnaorq7GuHHjGtVLBzSupEld5UvUWDxhro2ItEHRog1zQ6IAMGXKlCYNSOu4aIOoZWhIr9ztbWlpaThy5EiDy5LMmDEDc+fOld9z0KBBGDFiBJKSkuDg4NDkbQDw4Ycf4r333kNQUBDs7e0bvEqXiydaDi7aICUU9fClpqaavL548SJOnTqFBx54oNUlfERkPdYcSr297caNGzhz5gxmzJiBp59+Gvb29oqHVp2cnJCamoqOHTvi999/N9sjp8bcOXO9ckp75NhLR2SbFCV8+/fvv6Nt3bp1OHbsWJMHRES2xdpJW0OHUm9vc3d3x9KlS9GxY0ecP38enTp1Upy0BQYGwsnJCceOHTNps+bcORcXF5OETelCCS6eILJNihK+ujz55JPw9fXFsmXLmjIeItIoJYlbcydttw6lNmTBw61tBoMBa9asQUFBAa5cuYLAwEDFSZuPjw/8/f3x6aefIjAwUNXCw+baiIhupSjhq6mpMXldUVEhV1AnItumZOFCXfPEmiNpa8xQ6u1tPj4+iImJwWeffQaDwQCDwYDjx48rTtquXr2KgIAA7NixA9u3b2+SXjqAvXJE1DCKEr42bdrcsctGYGAgVq9erUpQRKQOtWrC3b6bQ3MlbY0ZSq2rbdasWfj++++xYcMGHDp0CG3atFGctA0ZMgSzZs1CRUUFe+mIqNkpSvjOnDlj8trNzQ2+vr6qBERE/6cpVpqqXRPOYDCgtLQUCQkJiIuLQ8+ePZstaWvsUGpdbbX3Nn78+AYlbeylIyItUJTwhYSEqB0HkU2rb66bmitNrVETrq7dHJozaWvMUOrdhldvxaSNiFqSehO+wYMH3zGUeytJkrBv374mD4pIi5qiZ01pTbSmWmnaFAsZbt25oa7Era7dHHQ6XbMlbU0xlMrhVSKyNfUmfBMnTqyzvaCgAB988AEqKipUCYpIS5pyt4Xb57qpvdLUWjXh6trNobmTtsYMpbKnjohsTb0J39SpU01e6/V6/O1vf8Pq1asRHx+PBQsWqBocUVOz5kb2tyduer3+jrluaq80tVZNuJqaGgQEBOCrr77C3r17NZG0ERHR/1E0h89gMGDZsmVISUnBI488gqysLHTp0kXt2IiaTEN76Rqzkf3tiVtdc92ssdLUmjXhmLQREWlTvQlfZWUl/v73vyM5ORmDBg3CN998g4iICGvFRtRkq1Qb2kvXmI3sb0/c6prrZq2VptaqCcekjYhIm+pN+EJDQ1FdXY1XX30VvXv3RlFREYqKikzOiYmJUTVAap0aM2+uKXvplCxaAJT3rNU1180aK01ZE46IqHWrN+FzdnaGJElYvnx5ncclScLp06dVCYxsi6W9cmlpaThy5IjFPXJN3UvX2I3sb0/c6prrZs2VphxeJSJqnSQhhGjuIFoSSZLAP5lyDempu3HjBs6cOYMZM2Zg9uzZsLe3h16vx6BBgzBixAgkJSXBwcHBKm0A8OGHH+K9995DUFAQ7O3tLeptrB0OvT1xY68aETUVSQL4tUR3o2jRBlFD1SZ7lvTUubu7Y+nSpejYsSPOnz+PTp06NWqVamN76ZpiI/vbe9bYq0ZERNbEhI9Uo9fr5WTPkrlzBoMBa9asQUFBAa5cuYLAwMBGrVIFlJUW4Ub2RERkq5jwkWKWzsO7evUqhBCIioqS30NJr5yPjw9iYmLw2WefwWAwwGAw4Pjx441apdoUvXREREQtFefwWag1zuFrihp206dPBwDFc+cqKysxbtw4/PDDD+jSpQvatGnTqL1kze2HSkTU0nEOHynBHj66q4bMw7u1bc2aNbCzs7O4V652V4rx48dzP1QiIqJGYA+fhVpbD59er8fYsWNN5uFZe4Ure+WIiMxjDx8pwR4+qlftdmCWzsOrTfj69euHkJAQ/PWvf4WXlxd75YiIiJqBXXMHQNal1+uRk5MDvV6vqM3FxQWSJCErK0s+dutKWCcnJ7NtAJCVlQVJktCjRw/88Y9/vKM0iZI2IiIiahz28LUSjdmqzM7ODmlpaQBg8Ty89PR0xMTEMIEjIiJqRpzDZ6GWOodv8eLFDV54sXXrVjg4OKCmpoarY4mINIZz+EgJ9vC1UEpq4tW+dnFxaVABZAcHB0RHRwO4WYh45cqVuH79OufhERERtTBM+FoYJUOzNTU1cHBwgNFohCRJKC8vR2VlJSIjI+X3sWThRVRUFNLS0nD9+nX88Y9/NImHO00QERFpHxdttDC31sRbunQp4uPjkZqail27dsltXl5eKCwsRFxcnHxOSUkJMjMz5fdpyMKLgIAAq98vERERNR57+FoQJXvTGgwGuWBxXFwcevbsiejoaBw8eBCbN2+Gl5cXoqOjufCCiIioFWHC14IoqYlXe06fPn0AQB6affnllzF58mSkpqZi+/btkCQJkyZNAnBzfl5aWprZttqFF0RERNQyMeFrQQIDA+WaeLWLKW4dho2Pj5fPOXr0KOLi4uSh2V9++QWhoaFceEFERNQKMeFrQXx8fBATE4P09HQA5mviubu7Y+PGjaiurobRaDQZlu3atWud78uFF0RERLaLdfgs1Nx1+BqySpf18IiIbBfr8JESTPgs1NwJXy0lw7AcliUisn1M+EgJJnwW0krCR0REBDDhI2VYh4+IiIjIxjHhIyIiIrJxTPiIiIiIbBwTPiIiIiIbx4SPiIiIyMYx4SMiIiKyca0+4SspKcHo0aPh5uaGkJAQbN26tblDuoNer0dOTg70en1zh0JEREQtUKvfWm3GjBlwdHREUVERfvrpJzz88MOIjIxEREREc4dW564a3DGDiIiILNWqCy+Xl5fD29sbP//8M8LCwgAAkyZNQmBgIN5+++06r7Fm4eXFixcjMzMT8fHx8r65tXvizps3zyoxEBGRtrHwMinRqnv4Tp48CXt7eznZA4DIyEh8/fXXzRjVTXq9Xk72EhMTAQDR0dEAgPT0dEyfPp3bpREREZEirTrhKysrg6enp0mbp6cnSktLTdoyMjKQkZFhzdBQUFAAIQSioqJM2qOiopCWlobCwkImfERERKRIq160odPpYDAYTNoMBgPc3d1N2mJjY/H222+bHeZVQ2BgICRJQlZWlkl7VlYWJElCQECA1WIhIiLt4nAuKdGqe/jCwsJgNBqRl5eHrl27AgCys7M1sWDDx8cHMTExSE9PB4A75vCxd4+IiIiUatWLNgBg3LhxkCQJa9aswU8//YRhw4bh8OHDZpM+ay7a4CpdIiIiagqtPuErKSnBlClTsHfvXvj4+ODtt9/GE088YfZ8ayZ8tfR6PQoLCxEQEMCePSIiIrJYq0/4LNUcCR8RERFRY7TqRRtERERErQETPiIiIiIbx4SPiIiIyMYx4SMiIiKycUz4iIiIiGwcEz4iIiIiG8eEj4iIiMjGMeEjIiIisnGtei/dhpIkqblDICIiIgIARRtCMOGzUHPtstG7d2/88MMPzfLZ9H/4HLSBz0Eb+By0gc9BG7T+HDikS0RERGTjmPARERER2Tj7pKSkpOYOgpTp1atXc4dA4HPQCj4HbeBz0AY+B23Q8nOQRHNNSiMiIiIiq+CQLhEREZGNY8JHREREZOOY8GlESUkJRo8eDTc3N4SEhGDr1q11nieEwOzZs+Hj4wMfHx+8+uqrzVYqxhYpfQ7Lli1Dz5494e7ujtDQUCxbtszKkdo2pc+h1o0bNxAeHo6goCArRdg6WPIcsrKyMHDgQOh0Ovj5+eH999+3YqS2T+mzqKqqwrPPPgs/Pz+0bdsWw4cPR0FBgZWjtU0pKSno3bs3nJyc8OSTT9Z77nvvvQd/f394enpiypQpqKqqsk6Q9WDCpxEzZsyAo6MjioqKsGXLFkyfPh25ubl3nLdq1Srs3LkT2dnZyMnJwa5du7By5cpmiNg2KX0OQghs2rQJV65cwZ49e5CSkoK0tLRmiNg2KX0OtZYtW4b27dtbMcLWQelzKC4uxtChQ5GYmAi9Xo9ff/0VcXFxzRCx7VL6LN5//318++23yMnJQWFhIby8vDBz5sxmiNj2BAQEYN68eZgyZUq95/373//G22+/jX379uHs2bM4ffo0Fi5caKUo6yGo2ZWVlQkHBwdx4sQJuW3ixIli9uzZd5x73333iZUrV8qv16xZI/r27WuVOG2dJc/hdjNnzhTPP/+8muG1GpY+h9OnT4vw8HDx5ZdfisDAQGuFafMseQ6vvfaamDhxojXDa1UseRbPPvuseOWVV+TXu3btEmFhYVaJs7V4/fXXRUJCgtnj48ePF6+99pr8OiMjQ/j5+Vkhsvqxh08DTp48CXt7e4SFhcltkZGRdf7rLTc3F5GRkXc9jyxnyXO4lRAChw4dQkREhNohtgqWPoeZM2firbfegouLi7VCbBUseQ5HjhxB27Ztcf/996N9+/YYPnw4zp07Z81wbZolz2Lq1Kn4z3/+g8LCQlRUVGDLli146KGHrBluq1fX93RRURH0en0zRsUhXU0oKyuDp6enSZunpydKS0vveq6npyfKyso4j68JWPIcbpWUlISamho89dRTaobXaljyHD799FMYjUaMHj3aWuG1GpY8hwsXLmDjxo14//33ce7cOYSGhmL8+PHWCtXmWfIswsLC0LFjRwQGBsLDwwPHjh3DggULrBUqoe7vaQB3/S5RGxM+DdDpdDAYDCZtBoMB7u7udz3XYDBAp9NBkiTV47R1ljyHWikpKdi0aRO++OILODk5qR1iq6D0OZSXl+PVV1/Fhx9+aM3wWg1L/ntwcXHB6NGjER0dDWdnZyxcuBCHDx/GtWvXrBWuTbPkWUyfPh3Xr1+HXq9HeXk5xowZwx4+K6vrexpAvd8l1sCETwPCwsJgNBqRl5cnt2VnZ9c5RBgREYHs7Oy7nkeWs+Q5AMC6devkiblcHdp0lD6HvLw8nD17FgMGDIC/vz/GjBmD3377Df7+/jh79qyVo7Y9lvz3cM8995j8o7P2d448NA1LnkV2djaefPJJtG3bFk5OTpg5cyaOHj2K4uJia4bcqtX1Pe3n5wcfH59mjApctKEV8fHxYty4caKsrEx88803wsPDQ/z88893nLd8+XIRHh4uLly4IAoKCkSPHj3E8uXLmyFi26T0OWzevFn4+fmJX375pRmitH1KnsPvv/8ufvvtN/nn448/Fh06dBC//fabMBqNzRS5bVH638O+ffuEl5eX+PHHH8WNGzfECy+8IPr3798MEdsupc/iySefFGPGjBFXr14VN27cEG+++aYICAhohohtz++//y4qKyvFnDlzxMSJE0VlZaX4/fff7zhv9+7dws/PT+Tm5oqSkhIxePBgRYv/1MaETyP0er0YOXKkcHV1FcHBwWLLli1CCCEOHjwo3Nzc5PNqamrEK6+8Iry9vYW3t7d45ZVXRE1NTXOFbXOUPodOnTqJNm3aCDc3N/knMTGxucK2OUqfw63279/PVbpNzJLn8NFHH4mAgADh5eUlHnnkEXHu3LnmCNlmKX0WxcXF4oknnhDt2rUTnp6e4oEHHhDfffddc4VtUxYuXCgAmPwsXLhQ5OfnCzc3N5Gfny+fm5ycLNq3by/c3d3Fk08+Ka5fv96Mkd/EvXSJiIiIbBzn8BERERHZOCZ8RERERDaOCR8RERGRjWPCR0RERGTjmPARERER2TgmfEREREQ2jgkfEWmGTqfD6dOnzR5/9tlnsWjRokZ9xoEDB7gzSgtXWVmJ4cOHw9PTE2PHjgUAzJs3D76+vvD392/m6Ii0iQkfEdWrU6dOcHFxgU6nk3+ef/55VT6rrKwMnTt3Nnt8xYoVmD9/viqfXUuSJLi5uUGn08HX1xfjx4/H1aFC2pkAAAp9SURBVKtXVf3M2124cAGPPvoofH194enpiT/+8Y/YsGGDVWOoy40bN5CUlISuXbvCzc0NnTp1wpQpUxRtZXf27FlIkgSj0VjveUIILFu2DF27doWLiws6duyIOXPmoKqqSj5nx44dKCoqgl6vx/bt23H+/HkkJyfjl19+wcWLFxt7m0Q2iQkfEd3V559/jrKyMvknJSXF6jFUV1db7bOys7NRVlaG06dP48qVK0hKSrLaZwPApEmTEBwcjPz8fOj1emzatAl+fn5N+hl3S7zq8thjj+Gzzz7D1q1bce3aNWRnZ6NXr17Yt29fk8X1l7/8BatWrcKmTZtQWlqK3bt3IzMzE48//rh8Tn5+PsLCwtCmTRv5tY+PD9q3b99kcRDZnGbe6YOINC4kJETs3bu3zmPr168X999/v3jhhReEp6enCA0NFf/5z3/E+vXrRVBQkGjXrp3YsGGDfH5CQoJITEwUsbGxQqfTiYEDB4qzZ8/KxwGIvLw8+dxnn31WPPTQQ8LV1VXs3btXJCQkiNdff10+f+fOnSIyMlK4u7uLzp07i927dwshhFi3bp0IDw8XOp1OhIaGihUrVsjX3G0LtltjEEKIf/zjH2LIkCHy6/re+/Lly+Lhhx8Wnp6ewtvbW/Tv319UV1cLIYQoKCgQY8aMEb6+vqJTp07i/fffNxuDm5ub+PHHH80eP3TokLjvvvuEp6enCAoKEuvXrxdCCHH16lUxadIk4evrKzp27CgWLVokf/6tz8rb21v+O65du1aEh4cLLy8vERcXZ/I8brV3717h7Oxc75Zpt/9/ZeHChWLChAlCCCGCg4MFAHkrwsOHD99x/cmTJ4Wdnd0dW4GdO3dOODo6in379okFCxYIBwcHeWvDFStWCGdnZyFJknBzcxMJCQlm4yNqzdjDR0SN8t133+Gee+6BXq/HE088gXHjxuH777/Hr7/+is2bN+P5559HWVmZfP6WLVswf/58FBcX495778WECRPMvvfWrVvx+uuvo7S0FP379zc5dvToUUyePBnLli3D1atXcfDgQXTq1AkA0L59e+zatQsGgwHr16/Hiy++iKysLIvv7cqVK9i5cyf69esnt9X33snJyQgKCsLly5dRVFSEt956C5IkoaamBsOHD0dkZCQKCgqwb98+/P3vf8e///3vOj+3X79+mDFjBtLS0nDu3DmTY+fOncNDDz2EmTNn4vLly/jpp59w7733AgBmzpyJa9eu4fTp0/j666+xadMmrF+/Xr72u+++Q+fOnXHp0iW8/vrr2LlzJ9566y188sknuHz5MgYMGIDx48fXGVNGRgb69OmD4OBgi/+OAHDw4EEAwNWrV1FWVob77rvvjnP27duHoKAg9OnTx6Q9ODgY/fr1w969e/HGG29g7ty5iI+PR1lZGRITE7F7924EBASgrKxME0PfRFrEhI+I7mrUqFHw8vKSf1avXi0fCw0NxVNPPQV7e3vEx8fj/PnzWLBgAZycnBAXFwdHR0f8+uuv8vkPP/wwBg4cCCen/9/O3YU02b8BHP/Ouf6Zm+7NA7VlTokiodGBvdhQakoIvQzBSj0IU7IgOsiQioqw0g70oILyoCSisDSKjAmRRwaJL2HgICiLuZxIKhPDt1z7Hzw8w9WWJc8L7Lk+sIP77bqu372Ti9/9u+//cenSJV69eoXb7Q6Zd8+ePWRlZREVFcXy5cuDjt26dYvS0lJyc3OJiooiOTmZtWvXBnKkpaWhUCjIzs4mLy+Pjo6OXx7vxo0b0Wq1GI1GBgcHOXz4cFD94WKrVCqGh4dxuVyoVCqsVisKhYLu7m4+f/7MuXPnWLZsGWazmfLycpqamkLmb25uxmq1Ul1dTWpqKhaLhe7ubuCPhtlms3HgwAFUKhUGgwGLxYLP5+PBgwfU1NSg0WhYvXo1J06c4O7du4G4SUlJHDt2jOjoaGJiYmhoaODUqVOsW7eO6OhoTp8+TV9fHy6X64eaxsbGSExM/OV7uBSjo6NhcyQmJjI6Ovq35hcikknDJ4RY1JMnT/B6vYFfeXl54NjCtWUxMTEh9y2c4Vs4Q6RWq9Hr9Xg8npB5fzab5Ha7SUtLC3msra2NzZs3o9fr0Wq1OByO32oWXr9+jdfrZWZmhiNHjmC1WpmZmVk09smTJ0lPTycvLw+z2UxtbS3wxxozj8cT1DRfvnyZkZGRkPl1Oh21tbU4nU5GRkawWCzs3bsXv98fdtyjo6PMzc2RkpIS2JeSksLQ0FBg+/v76XK5OH78eKAmvV6P3+8PuuZPBoOB4eHhX76Hv2L9+vWBF4E6OjowGo1hcwwPD2M0Gv/S/EL8l0jDJ4T4Ry2czfvy5Qvj4+MkJSWFPFehUISNYzKZGBgY+GH/7OwsBQUFVFZWMjIygtfrJT8/H7/f/9u1qlQqysrK+PjxI/39/YvG1mg01NXV8eHDB1pbW6mvr6e9vR2TyURqampQ0zw5OYnD4Vi0BqPRSGVlJR6Ph/Hx8bDjNhqNqFSqoNm5wcFBkpOTA9vf30+TyURDQ0NQXdPT02zduvWH+Dabja6uLj59+hS21tjYWKampgLbC9+YDfVfOp3OwItAVquV7du343a76erqCjrP7XbT2dnJjh07wuYWQvycNHxCiH+Uw+Hg5cuXzM3NcfbsWTZt2rSkdWGHDh2isbGR9vZ2vn37xtDQEG/fvmVubo7Z2VkSEhKIjo6mra2N58+fL6lWn89HY2MjMTExmM3mRWM/e/aM9+/f4/f7iYuLQ6lUolQqyczMJC4ujitXrjA9PY3P56O/vz/wmPZ7VVVV9Pf3Mz8/z+TkJDdu3CA9PR2DwUBxcTEvXrzg4cOHzM/PMzY2Rl9fH0qlksLCwsCaR5fLRX19PSUlJWHHV1FRQU1NDU6nE4CJiQmam5tDnmuz2cjNzcVut9Pb2xuo7ebNm9y+fRsAi8VCU1MTX79+paenh5aWlsD1CQkJREVF/fQ7i2vWrKGiooLi4mI6Ozvx+Xw4nU4KCgqw2WzYbLbwf5YQ4qek4RNCLGrXrl1B3+Gz2+1LjlVUVMSFCxfQ6/X09vZy7969JcXJzMwMvDQRHx9PdnY2LpcLjUbD1atXKSwsRKfTcf/+fXbv3v1bsTds2IBarUan03Hnzh0eP36MXq9fNPa7d++w2Wyo1Wq2bNnC0aNHycnJQalU0traSl9fH6mpqRiNRsrKypiYmAiZf2pqCrvdjlarxWw243K5ePr0KQCrVq3C4XBQV1eHXq/HYrHw5s0bAK5du0ZsbCxms5lt27ZRVFREaWlp2HHa7XaqqqrYv38/cXFxZGRk0NbWFvb8lpYW8vPz2bdvH/Hx8WRkZNDT0xNoxKqrqxkYGECn03H+/HmKiooC165YsYIzZ86QlZWFVquls7MzZI7r169TVlZGSUkJarWanTt3kpOTw6NHj8LWJYRYnMK/lOccQgixBAcPHmTlypVcvHjx3y5FCCH+U2SGTwghhBAiwknDJ4QQQggR4eSRrhBCCCFEhJMZPiGEEEKICCcNnxBCCCFEhJOGTwghhBAiwknDJ4QQQggR4aThE0IIIYSIcNLwCSGEEEJEuP8DV8/I46W1eOAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x378 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig , ax = plt.subplots()\n",
    "\n",
    "#SNP found in isolate denoted with a '1', sum across columns (isolates)\n",
    "plt.scatter(thresholds , num_positions_dropped , color = '1.0' , edgecolor = 'k' , linewidth = 1.5 , alpha = 0.75)\n",
    "plt.plot([0.8 , 0.8] , [0 , np.max(num_positions_dropped)] , color = 'blue' , linewidth = 1.0)\n",
    "\n",
    "#plt.title('Distribution of pairwise SNP distances between 10,024 isolates', fontsize = 12, color = 'k')\n",
    "plt.xlabel('Empirical Base Score Cut-Off', fontsize = 12, color = 'k')\n",
    "plt.ylabel('Number of H37Rv Positions Dropped' , fontsize = 12, color = 'k')\n",
    "\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.grid(False)\n",
    "\n",
    "ax.tick_params(labelcolor = 'k')\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(9.0, 5.25)\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.tick_params(axis='y', which='major', labelsize=12 , labelcolor = 'k')\n",
    "plt.tick_params(axis='x', which='major', labelsize=12 , labelcolor = 'k')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bad Regions for IL mapping/variant calling: sites with empirical score < 0.80"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(129215, 2)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape( emp_basepair_recall_df[emp_basepair_recall_df.EmpiricalBasepairRecall_31CI < 0.80] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "H37Rv_positions_to_drop = list( emp_basepair_recall_df[emp_basepair_recall_df.EmpiricalBasepairRecall_31CI < 0.80].pos )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "129215"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(H37Rv_positions_to_drop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Export pickled list of H37Rv positions to drop when calling variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/H37Rv_sites_to_drop.pkl', 'wb') as f:\n",
    "    pickle.dump(H37Rv_positions_to_drop, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use to load the list for use in variant calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/H37Rv_sites_to_drop.pkl', 'rb') as f:\n",
    "    H37Rv_positions_to_drop = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}