src/inDelphi_Check_Data_Leakage.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from config import *\n",
    "from utils import *\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import copy\n",
    "import numpy as np\n",
    "import collections\n",
    "import multiprocessing\n",
    "import pickle\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "# Suppress pandas future warning, which messes tqdm\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "import pandas as pd\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from Bio import pairwise2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check for inDelphi training data leakage\n",
    "inDelphi (Shen et al. 2018) is trained on a dataset of 55-bp sequences (available from their [GitHub](https://github.com/maxwshen/indelphi-dataprocessinganalysis/blob/master/SupplementaryData.xlsx)), referred to as \"lib-A\" in its paper. We are evaluating the performance of inDelphi on our lib-SA library of 61-bp sequences (specifically, the dat-A subset). For the evaluation to be meaningful, we need to make sure inDelphi lib-A sequences do not overlap and are not homologous to sequences from our dat-A."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## inDelphi's Lib-A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "libA_df = pd.read_excel(os.path.join(DATA_DIR, 'indelphiLibA', 'SupplementaryData.xlsx'), header=1, sheet_name='Supplementary Table 2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Index number</th>\n",
       "      <th>Name</th>\n",
       "      <th>gRNA</th>\n",
       "      <th>Sequence Context</th>\n",
       "      <th>SpCas9 PAM</th>\n",
       "      <th>Total readcount, mESCs</th>\n",
       "      <th>Total number of unique high-confidence outcomes, mESCs</th>\n",
       "      <th>Total readcount, U2OS cells</th>\n",
       "      <th>Total number of unique high-confidence outcomes, U2OS cells</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0_0_0_0</td>\n",
       "      <td>CTTTCACTTTATAGATTTAT</td>\n",
       "      <td>GTGCTCTTAACTTTCACTTTATAGATTTATAGGGTTAATAAATGGG...</td>\n",
       "      <td>AGG</td>\n",
       "      <td>5498.984841</td>\n",
       "      <td>31</td>\n",
       "      <td>501.000000</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0_0_0_0</td>\n",
       "      <td>AATCCTCTAAATAGTTAAAT</td>\n",
       "      <td>CCTTTGTACCAATCCTCTAAATAGTTAAATAGGAAGGATAACCCAA...</td>\n",
       "      <td>AGG</td>\n",
       "      <td>30592.535156</td>\n",
       "      <td>62</td>\n",
       "      <td>5931.267377</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0_0_0_0</td>\n",
       "      <td>AAATATCTTTAACCTAAAAC</td>\n",
       "      <td>GTAAAAGTTAAAATATCTTTAACCTAAAACCGGGTTCAGCGACTTA...</td>\n",
       "      <td>CGG</td>\n",
       "      <td>13346.723123</td>\n",
       "      <td>42</td>\n",
       "      <td>524.000000</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>0_0_0_1</td>\n",
       "      <td>TGTCATCTTTATGCTTAAAG</td>\n",
       "      <td>ACTGTTAAATTGTCATCTTTATGCTTAAAGAGGCTTAAGTTACTTC...</td>\n",
       "      <td>AGG</td>\n",
       "      <td>197322.665139</td>\n",
       "      <td>206</td>\n",
       "      <td>14094.000000</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>0_0_0_1</td>\n",
       "      <td>TTTAAGCACAGCAAGGCTCT</td>\n",
       "      <td>GATATTATAATTTAAGCACAGCAAGGCTCTTGGCTAGAATTATACC...</td>\n",
       "      <td>TGG</td>\n",
       "      <td>129094.351837</td>\n",
       "      <td>121</td>\n",
       "      <td>15887.000000</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Index number     Name                  gRNA  \\\n",
       "0             0  0_0_0_0  CTTTCACTTTATAGATTTAT   \n",
       "1             1  0_0_0_0  AATCCTCTAAATAGTTAAAT   \n",
       "2             2  0_0_0_0  AAATATCTTTAACCTAAAAC   \n",
       "3             3  0_0_0_1  TGTCATCTTTATGCTTAAAG   \n",
       "4             4  0_0_0_1  TTTAAGCACAGCAAGGCTCT   \n",
       "\n",
       "                                    Sequence Context SpCas9 PAM  \\\n",
       "0  GTGCTCTTAACTTTCACTTTATAGATTTATAGGGTTAATAAATGGG...        AGG   \n",
       "1  CCTTTGTACCAATCCTCTAAATAGTTAAATAGGAAGGATAACCCAA...        AGG   \n",
       "2  GTAAAAGTTAAAATATCTTTAACCTAAAACCGGGTTCAGCGACTTA...        CGG   \n",
       "3  ACTGTTAAATTGTCATCTTTATGCTTAAAGAGGCTTAAGTTACTTC...        AGG   \n",
       "4  GATATTATAATTTAAGCACAGCAAGGCTCTTGGCTAGAATTATACC...        TGG   \n",
       "\n",
       "   Total readcount, mESCs  \\\n",
       "0             5498.984841   \n",
       "1            30592.535156   \n",
       "2            13346.723123   \n",
       "3           197322.665139   \n",
       "4           129094.351837   \n",
       "\n",
       "   Total number of unique high-confidence outcomes, mESCs  \\\n",
       "0                                                 31        \n",
       "1                                                 62        \n",
       "2                                                 42        \n",
       "3                                                206        \n",
       "4                                                121        \n",
       "\n",
       "   Total readcount, U2OS cells  \\\n",
       "0                   501.000000   \n",
       "1                  5931.267377   \n",
       "2                   524.000000   \n",
       "3                 14094.000000   \n",
       "4                 15887.000000   \n",
       "\n",
       "   Total number of unique high-confidence outcomes, U2OS cells  \n",
       "0                                                  1            \n",
       "1                                                 13            \n",
       "2                                                  3            \n",
       "3                                                 10            \n",
       "4                                                 24            "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "libA_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "libA_seqs = libA_df['Sequence Context'].unique().tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1998"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(libA_seqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our dat-A Target Sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Azimuth score</th>\n",
       "      <th>Cas type</th>\n",
       "      <th>Cutsite within sequence</th>\n",
       "      <th>Designed 61-bp target site (37i-24e, AG)</th>\n",
       "      <th>Designed 61-bp target site (37i-24e, CT reverse-complement)</th>\n",
       "      <th>Designed gRNA (NGG orientation, 19 and 20)</th>\n",
       "      <th>Has AG</th>\n",
       "      <th>MaxEntScan AG disruption frequency in MH dels and insertions genomic_context</th>\n",
       "      <th>MaxEntScan AG disruption frequency in MH dels and insertions library_context</th>\n",
       "      <th>Ontarget score</th>\n",
       "      <th>...</th>\n",
       "      <th>inDelphi Precision library_context</th>\n",
       "      <th>inDelphi expected indel length genomic_context</th>\n",
       "      <th>inDelphi expected indel length library_context</th>\n",
       "      <th>inDelphi highest outcome frequency genomic_context</th>\n",
       "      <th>inDelphi highest outcome frequency library_context</th>\n",
       "      <th>inDelphi logphi genomic_context</th>\n",
       "      <th>inDelphi logphi library_context</th>\n",
       "      <th>strand</th>\n",
       "      <th>total_vars</th>\n",
       "      <th>Designed gRNA (NGG orientation)</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Identifier number</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.685572</td>\n",
       "      <td>SpCas9</td>\n",
       "      <td>68</td>\n",
       "      <td>GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...</td>\n",
       "      <td>GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...</td>\n",
       "      <td>ATCCAGGTCAGGACTGCAGG</td>\n",
       "      <td>yes</td>\n",
       "      <td>0.796188</td>\n",
       "      <td>0.796141</td>\n",
       "      <td>0.858865</td>\n",
       "      <td>...</td>\n",
       "      <td>0.426247</td>\n",
       "      <td>7.552535</td>\n",
       "      <td>7.546523</td>\n",
       "      <td>16.796109</td>\n",
       "      <td>16.798237</td>\n",
       "      <td>-2.441151</td>\n",
       "      <td>-2.441337</td>\n",
       "      <td>+</td>\n",
       "      <td>12</td>\n",
       "      <td>GATCCAGGTCAGGACTGCAGG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.666555</td>\n",
       "      <td>SpCas9</td>\n",
       "      <td>69</td>\n",
       "      <td>GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...</td>\n",
       "      <td>GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...</td>\n",
       "      <td>ATCCAGGTCAGGACTGCAG</td>\n",
       "      <td>yes</td>\n",
       "      <td>0.760537</td>\n",
       "      <td>0.760490</td>\n",
       "      <td>0.844692</td>\n",
       "      <td>...</td>\n",
       "      <td>0.527270</td>\n",
       "      <td>8.262278</td>\n",
       "      <td>8.256890</td>\n",
       "      <td>37.627192</td>\n",
       "      <td>37.633988</td>\n",
       "      <td>-2.109904</td>\n",
       "      <td>-2.110079</td>\n",
       "      <td>+</td>\n",
       "      <td>12</td>\n",
       "      <td>TGATCCAGGTCAGGACTGCAG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.621866</td>\n",
       "      <td>SpCas9</td>\n",
       "      <td>70</td>\n",
       "      <td>GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...</td>\n",
       "      <td>GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...</td>\n",
       "      <td>TGATCCAGGTCAGGACTGCA</td>\n",
       "      <td>yes</td>\n",
       "      <td>0.793351</td>\n",
       "      <td>0.793299</td>\n",
       "      <td>0.854025</td>\n",
       "      <td>...</td>\n",
       "      <td>0.572962</td>\n",
       "      <td>6.545344</td>\n",
       "      <td>6.540914</td>\n",
       "      <td>28.429830</td>\n",
       "      <td>28.436377</td>\n",
       "      <td>-2.139149</td>\n",
       "      <td>-2.139360</td>\n",
       "      <td>+</td>\n",
       "      <td>12</td>\n",
       "      <td>ATGATCCAGGTCAGGACTGCA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.434994</td>\n",
       "      <td>SpCas9</td>\n",
       "      <td>71</td>\n",
       "      <td>GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...</td>\n",
       "      <td>GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...</td>\n",
       "      <td>ATGATCCAGGTCAGGACTGC</td>\n",
       "      <td>yes</td>\n",
       "      <td>0.699063</td>\n",
       "      <td>0.698923</td>\n",
       "      <td>0.781274</td>\n",
       "      <td>...</td>\n",
       "      <td>0.523115</td>\n",
       "      <td>8.410657</td>\n",
       "      <td>8.406734</td>\n",
       "      <td>37.406834</td>\n",
       "      <td>37.408835</td>\n",
       "      <td>-2.126282</td>\n",
       "      <td>-2.126381</td>\n",
       "      <td>+</td>\n",
       "      <td>12</td>\n",
       "      <td>GATGATCCAGGTCAGGACTGC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.535781</td>\n",
       "      <td>SpCas9</td>\n",
       "      <td>68</td>\n",
       "      <td>TGCCAGGCATGCCCACTGACTGGTGCCCTTGGCCCAGGGTGGAGAC...</td>\n",
       "      <td>AGTATTGGAATGAAGGTCTCCACCCTGGGCCAAGGGCACCAGTCAG...</td>\n",
       "      <td>TGACTGGTGCCCTTGGCCCA</td>\n",
       "      <td>yes</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.860908</td>\n",
       "      <td>...</td>\n",
       "      <td>0.443413</td>\n",
       "      <td>7.775675</td>\n",
       "      <td>7.788979</td>\n",
       "      <td>10.368532</td>\n",
       "      <td>10.361942</td>\n",
       "      <td>-2.355018</td>\n",
       "      <td>-2.354383</td>\n",
       "      <td>+</td>\n",
       "      <td>14</td>\n",
       "      <td>CTGACTGGTGCCCTTGGCCCA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 34 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                   Azimuth score Cas type  Cutsite within sequence  \\\n",
       "Identifier number                                                    \n",
       "2                       0.685572   SpCas9                       68   \n",
       "3                       0.666555   SpCas9                       69   \n",
       "4                       0.621866   SpCas9                       70   \n",
       "5                       0.434994   SpCas9                       71   \n",
       "6                       0.535781   SpCas9                       68   \n",
       "\n",
       "                            Designed 61-bp target site (37i-24e, AG)  \\\n",
       "Identifier number                                                      \n",
       "2                  GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...   \n",
       "3                  GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...   \n",
       "4                  GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...   \n",
       "5                  GGAAATGCCCTGCCTAAGCTTTGGGTGCCCCCTGCAGTCCTGACCT...   \n",
       "6                  TGCCAGGCATGCCCACTGACTGGTGCCCTTGGCCCAGGGTGGAGAC...   \n",
       "\n",
       "                  Designed 61-bp target site (37i-24e, CT reverse-complement)  \\\n",
       "Identifier number                                                               \n",
       "2                  GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...            \n",
       "3                  GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...            \n",
       "4                  GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...            \n",
       "5                  GCAGGATGATGATCCAGGTCAGGACTGCAGGGGGCACCCAAAGCTT...            \n",
       "6                  AGTATTGGAATGAAGGTCTCCACCCTGGGCCAAGGGCACCAGTCAG...            \n",
       "\n",
       "                  Designed gRNA (NGG orientation, 19 and 20) Has AG  \\\n",
       "Identifier number                                                     \n",
       "2                                       ATCCAGGTCAGGACTGCAGG    yes   \n",
       "3                                        ATCCAGGTCAGGACTGCAG    yes   \n",
       "4                                       TGATCCAGGTCAGGACTGCA    yes   \n",
       "5                                       ATGATCCAGGTCAGGACTGC    yes   \n",
       "6                                       TGACTGGTGCCCTTGGCCCA    yes   \n",
       "\n",
       "                   MaxEntScan AG disruption frequency in MH dels and insertions genomic_context  \\\n",
       "Identifier number                                                                                 \n",
       "2                                                           0.796188                              \n",
       "3                                                           0.760537                              \n",
       "4                                                           0.793351                              \n",
       "5                                                           0.699063                              \n",
       "6                                                           1.000000                              \n",
       "\n",
       "                   MaxEntScan AG disruption frequency in MH dels and insertions library_context  \\\n",
       "Identifier number                                                                                 \n",
       "2                                                           0.796141                              \n",
       "3                                                           0.760490                              \n",
       "4                                                           0.793299                              \n",
       "5                                                           0.698923                              \n",
       "6                                                           1.000000                              \n",
       "\n",
       "                   Ontarget score  ... inDelphi Precision library_context  \\\n",
       "Identifier number                  ...                                      \n",
       "2                        0.858865  ...                           0.426247   \n",
       "3                        0.844692  ...                           0.527270   \n",
       "4                        0.854025  ...                           0.572962   \n",
       "5                        0.781274  ...                           0.523115   \n",
       "6                        0.860908  ...                           0.443413   \n",
       "\n",
       "                  inDelphi expected indel length genomic_context  \\\n",
       "Identifier number                                                  \n",
       "2                                                       7.552535   \n",
       "3                                                       8.262278   \n",
       "4                                                       6.545344   \n",
       "5                                                       8.410657   \n",
       "6                                                       7.775675   \n",
       "\n",
       "                   inDelphi expected indel length library_context  \\\n",
       "Identifier number                                                   \n",
       "2                                                        7.546523   \n",
       "3                                                        8.256890   \n",
       "4                                                        6.540914   \n",
       "5                                                        8.406734   \n",
       "6                                                        7.788979   \n",
       "\n",
       "                  inDelphi highest outcome frequency genomic_context  \\\n",
       "Identifier number                                                      \n",
       "2                                                          16.796109   \n",
       "3                                                          37.627192   \n",
       "4                                                          28.429830   \n",
       "5                                                          37.406834   \n",
       "6                                                          10.368532   \n",
       "\n",
       "                   inDelphi highest outcome frequency library_context  \\\n",
       "Identifier number                                                       \n",
       "2                                                          16.798237    \n",
       "3                                                          37.633988    \n",
       "4                                                          28.436377    \n",
       "5                                                          37.408835    \n",
       "6                                                          10.361942    \n",
       "\n",
       "                  inDelphi logphi genomic_context  \\\n",
       "Identifier number                                   \n",
       "2                                       -2.441151   \n",
       "3                                       -2.109904   \n",
       "4                                       -2.139149   \n",
       "5                                       -2.126282   \n",
       "6                                       -2.355018   \n",
       "\n",
       "                   inDelphi logphi library_context strand total_vars  \\\n",
       "Identifier number                                                      \n",
       "2                                        -2.441337      +         12   \n",
       "3                                        -2.110079      +         12   \n",
       "4                                        -2.139360      +         12   \n",
       "5                                        -2.126381      +         12   \n",
       "6                                        -2.354383      +         14   \n",
       "\n",
       "                   Designed gRNA (NGG orientation)  \n",
       "Identifier number                                   \n",
       "2                            GATCCAGGTCAGGACTGCAGG  \n",
       "3                            TGATCCAGGTCAGGACTGCAG  \n",
       "4                            ATGATCCAGGTCAGGACTGCA  \n",
       "5                            GATGATCCAGGTCAGGACTGC  \n",
       "6                            CTGACTGGTGCCCTTGGCCCA  \n",
       "\n",
       "[5 rows x 34 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_design.head() # lib-SA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gRNA ID</th>\n",
       "      <th>Category</th>\n",
       "      <th>Genotype position</th>\n",
       "      <th>Inserted Bases</th>\n",
       "      <th>Length</th>\n",
       "      <th>Empirical frequency</th>\n",
       "      <th>Predicted frequency</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8495</td>\n",
       "      <td>ins</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "      <td>1</td>\n",
       "      <td>0.018119</td>\n",
       "      <td>0.009339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8495</td>\n",
       "      <td>ins</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.001777</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8495</td>\n",
       "      <td>ins</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000211</td>\n",
       "      <td>0.004559</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>8495</td>\n",
       "      <td>ins</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C</td>\n",
       "      <td>1</td>\n",
       "      <td>0.049616</td>\n",
       "      <td>0.070707</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8495</td>\n",
       "      <td>del</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>0.167597</td>\n",
       "      <td>0.032842</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3210325</th>\n",
       "      <td>3256</td>\n",
       "      <td>del</td>\n",
       "      <td>56.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>60</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3210326</th>\n",
       "      <td>3256</td>\n",
       "      <td>del</td>\n",
       "      <td>57.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>60</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3210327</th>\n",
       "      <td>3256</td>\n",
       "      <td>del</td>\n",
       "      <td>58.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>60</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3210328</th>\n",
       "      <td>3256</td>\n",
       "      <td>del</td>\n",
       "      <td>59.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>60</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3210329</th>\n",
       "      <td>3256</td>\n",
       "      <td>del</td>\n",
       "      <td>60.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>60</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3210330 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         gRNA ID Category  Genotype position Inserted Bases  Length  \\\n",
       "0           8495      ins                NaN              A       1   \n",
       "1           8495      ins                NaN              G       1   \n",
       "2           8495      ins                NaN              T       1   \n",
       "3           8495      ins                NaN              C       1   \n",
       "4           8495      del                0.0            NaN       1   \n",
       "...          ...      ...                ...            ...     ...   \n",
       "3210325     3256      del               56.0            NaN      60   \n",
       "3210326     3256      del               57.0            NaN      60   \n",
       "3210327     3256      del               58.0            NaN      60   \n",
       "3210328     3256      del               59.0            NaN      60   \n",
       "3210329     3256      del               60.0            NaN      60   \n",
       "\n",
       "         Empirical frequency  Predicted frequency  \n",
       "0                   0.018119             0.009339  \n",
       "1                   0.000000             0.001777  \n",
       "2                   0.000211             0.004559  \n",
       "3                   0.049616             0.070707  \n",
       "4                   0.167597             0.032842  \n",
       "...                      ...                  ...  \n",
       "3210325             0.000000             0.000000  \n",
       "3210326             0.000000             0.000000  \n",
       "3210327             0.000000             0.000000  \n",
       "3210328             0.000000             0.000000  \n",
       "3210329             0.000000             0.000000  \n",
       "\n",
       "[3210330 rows x 7 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datA_df = pd.read_csv(os.path.join(TABLES_DIR, 'datA_table.csv.gz'), compression='gzip')\n",
    "datA_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "datA_seqs = exp_design.loc[datA_df['gRNA ID'].unique()]['Designed 61-bp target site (37i-24e, AG)'].unique().tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1549"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(datA_seqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequence Identity Analysis\n",
    "For each target sequence in our dat-A, align it with every sequence in inDelphi's lib-A to determine the most similar sequence, and record the sequence identity. Plot distribution of such max sequence identities. If lib-SA sequences are dissimilar to inDelphi's lib-A sequences, then the distribution should be skewed towards lower max sequence identities.\n",
    "\n",
    "Local alignment (Smith Waterman) parameters: +1 match, -3 mismatch, -5 gap open, -2 gap extend. These are the same as the default parameters of BLAST's blastn-short program.\n",
    "\n",
    "Sequence identity is the definition used by BLAST: (# match positions in alignment(seq1, seq2))/(min(len(seq1), len(seq2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sequence_identity(seq1, seq2, alignment):\n",
    "    num_matches = pairwise2.format_alignment(*alignment).split('\\n')[1].count('|')\n",
    "    return num_matches / min(len(seq1), len(seq2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_seq_identity_libA(our_seq):\n",
    "    max_seq_identity = -1\n",
    "    for inDelphi_seq in libA_seqs:\n",
    "        # Using BLAST suite's blastn-short defaults:\n",
    "        # +1 match\n",
    "        # -3 mismatch\n",
    "        # -5 gap open\n",
    "        # -2 gap extend\n",
    "        alignment = pairwise2.align.localms(inDelphi_seq, our_seq, 1, -3, -5, -2)\n",
    "        identity = sequence_identity(inDelphi_seq, our_seq, alignment[0])\n",
    "        max_seq_identity = max(max_seq_identity, identity)\n",
    "    return max_seq_identity\n",
    "\n",
    "def compute_max_sequence_identities():\n",
    "    max_sequence_identities = []\n",
    "    \n",
    "    try:\n",
    "        p = multiprocessing.Pool(NUM_PROCESSES)\n",
    "        for max_seq_identity in tqdm(p.imap_unordered(max_seq_identity_libA, datA_seqs, chunksize=2), total=len(datA_seqs)):\n",
    "            max_sequence_identities.append(max_seq_identity)\n",
    "    finally:\n",
    "        p.close()\n",
    "        p.join()\n",
    "        \n",
    "    return max_sequence_identities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not pickle_exists(DAT_A_INDELPHI_SEQUENCE_IDENTITY):\n",
    "    max_sequence_identities = compute_max_sequence_identities()\n",
    "    save_var(max_sequence_identities, DAT_A_INDELPHI_SEQUENCE_IDENTITY)\n",
    "else:\n",
    "    max_sequence_identities = load_var(DAT_A_INDELPHI_SEQUENCE_IDENTITY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S2 FigA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_max_sequence_identities(max_sequence_identities):\n",
    "    plt.rcParams.update({'font.size': 12})\n",
    "    fig, ax = plt.subplots(figsize=(5,5))\n",
    "    sns.distplot(max_sequence_identities, kde=False, ax = ax)\n",
    "    ax.set(xlabel=\"Sequence Identity\",\n",
    "           ylabel='# of dat-A Target Sequences (' + str(len(datA_seqs)) + ' Total)',\n",
    "           title=\"Distribution of pairwise best aligned\\nsequence identity\\nbetween dat-A & inDelphi's Lib-A\")\n",
    "    median = np.median(max_sequence_identities)\n",
    "    plt.axvline(median, color='gray', linestyle='dotted')\n",
    "    plt.text(median + 0.01, 450, 'Median = ' + \"{:.2f}\".format(median))\n",
    "    plt.savefig(os.path.join(IMAGES_DIR, 'datA_indelphi_sequence_identity.png'), dpi=300, bbox_inches='tight')\n",
    "    plt.show()\n",
    "    print(\"Median sequence identity:\", np.median(max_sequence_identities))\n",
    "    print(\"Mean sequence identity:\", np.mean(max_sequence_identities))\n",
    "    print(\"Min sequence identity:\", np.min(max_sequence_identities))\n",
    "    print(\"Max sequence identity:\", np.max(max_sequence_identities))\n",
    "    print(\"Second largest sequence identity:\", np.sort(max_sequence_identities)[::-1][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVQAAAF0CAYAAAB8PsGDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5xU1f3/8debKl16CQIKFsRCFBV7b4ntq8bYIxprNNH80r62EM3X2BNTjBq7kagxajQmGo0ae0EFFQVEpXcpSwfh8/vj3IW7w+zOnd07e3dnP8/HYx47c+vnzsyeOffecz5HZoZzzrm6a5Z1AM45Vy68QHXOuZR4geqccynxAtU551LiBapzzqXEC1TnnEuJF6hFkDRS0qQSbXs/SSapb77XJdjfGZK+KsW2a0PSZpL+I2mZpJK25SvV5yjpJUl3pr3dnH3cK+n5Uu4jbdH3+NTY68mSLs8yprg0vw9NvkCNvqAWPb6StEDSG5J+LqlLzuI3AsOL2PYkSSMTLv460BuYmXT7CWPoGx3bfjmzHga+lua+6uhSoAcwlPA+lFJRn2MRjgV+WILt1qsavjNp2QX4dYm2nakWWQfQQLwCnED4gekM7Ab8BDhP0r5mNhHAzJYCS9PeuaRWZrYamJ32tqtjZiuAFfW1vwS2BN42s09LvaNCn6OkZoDMbG2R211Q19iaAjObl3UMJWNmTfoB3As8n2d6R+Az4IXYtJHApNjrvsDfgPmEwulz4MfRvJcAy3kMAPaLnn8TeBVYCVwYm943Wr/y9ZHA29Fy44CDY/uvsk5s+lfAGdHz3BgmR9PPAL7KWe8bwLvAKmAucCvQLve9As4BpgAVwN+B7gXe4w7A7cC86DhGA4fE5ufGeG812zkjOraDovdiZfTe7BRbpjPwZ2Bq9JlMAP4foYCs7nMcCUwCvg2Mj/axXRTLoNhyU4DpsdcDomW2jX3md8bm7wW8BiyJHmOBQ2Pze0bv6bxo/mvAPkm+r4Sa8AxgOeE72C1nuROBMdF7NBm4OeezrDa26r4z1cRzMvAWsJjwf/A0sFXOMgacGns9Gbg89ror8FdgGTAHuBq4j9j/ZeV7C1xBqHgsiN6LdkUed2vgj1G8C6Pnv4p/H+ryaPKn/NUxswrCm72fpO7VLHYr0InwDz4YOAuYHs07lvCB3kQ4he0NTIutexNwfbTeEzWEcjNwFfB14E3gSUnFnKrvFP09Lophl3wLSdoBeBJ4mXDa/R3gCOC2nEV3AfYn/CAcFi17Y4EY7gYOBU6NjuM14B+Stonm9wbeAEZFz39Qw7aaEd63C4BdCQX/05LaRvNbAx8CxwDbEv45f0EojGvSJ9rmGdF6UwiF8oEAkgYSLkl0krR1tM6BwGwz+zh3Y5KaE97PtwifwU6Egnt5NL8N8CLhx+bw6H35J/CcpMEFYt2V8GN6GOFHcAfCe1y57zMI392bomM5nfAdvS1JbCT8zkRaE97jnYCDgbWEz6NVgWOIuwfYkfB9O4BQUTkmz3LHA10Ix35ytMxPKmcWOu7ItdFxnQ7sTijEv1dErDWrj1pgQ35QTQ01mncY4dd11+j1SKrWbMYCI2vY9qTc+WyoVZ5WzfTcGupZsWVaEP7Rf5lvndhy8Rpq32iZ/XKWOYNYDRV4gHDKHV/maGAd0D/2Xs0DWseW+Rkwq4b3YFC0/2/kTH8PuDv2+iVitbtqtnVGtK0DY9M6E07fv1vDercAz8Ve536OI6Pj7Jfnu/FI9Pxs4D+EQu+CaNqDwKh8xxDFtdH7nnMs04EWOdNfAH5T4Pu6FOgUm3ZItK8to9eTgfNy1tsnWqZzgtjyfmcS/j91idbdMzat2hoq4VJP7mfaklD5yK2hfpCzr9uAN3K2W9NxtyPUXM/OWWY0XkOtF4r+VnfX+TfApZLeknSdpH2K2PbbCZd7o/KJmX0VrbdtEftJagihdhr3X8J7EN/fJ2a2KvZ6BuHUtTqV6+Zu++Von7URf08WAp9U7kdSM0k/kzRG0nxJS4HzgP4FtjnHzKbmTHsB2F+SCDWn/xBqlQdE8/ePltlIFNedwLOS/hXFtHVskV2AXsAiSUsrH8DehEKmJh+b2eLY69eiv4Ojs6n+wM052/1XtMygBLElJmmopMclfSFpCaFWD4Xf70qV3483KyeY2RpCIZdrTM7r9d+9JMcNDCTUqF/P2c6rCWMtyAvUmlVeR/s830wzu4fwId5GODX6l6Q/J9z2slrGpNjzdbnTotO52n6u1f1wxKevzjNPFE817K8226r0/4D/BX5HOAUdSig8Cp2C5vs8/gN0I5xSVxaelYXsEMJnnrdABTCzs4GdgeeAfYGPJJ0bzW5G+CEYmvMYTKgN11blZ/+DnO3uSCioP0wQWyLRZZZ/Ez7HMwmXInaJXhdzyg/Jvgv5vnuVx5vkuAtVkOrMC9RqSOoInA/8x8y+rG45M5tlZveY2emEa6inROtC+AI0r2Mo65v3SGpB+MJ+Ek2aG/3tE1t+KFULmMovYaE4xhH+seL2JXz5NrpGWIRx0d/c2vvesXnFir8nmwLbsOE92Qd4xszuMrP3zWwShWt8eZnZDOBT4CKgDaHW9D7h/b0EmGJmeX9sY9v4yMxuNrPDgbsIN/SItrUFUGFmk3IehZrODY59xwD2iP5+YmZzCKfLW+fZ7iQzW5kgtqTfmcFAd+AyM3vRzD4hnFoX8wNb+d3avXJC9D3fuYhtkPC4JxGObc+c1ffI3V5tebOpoJWkXoQvQmfCP+xPCKcH51e3kqTfE66pTQA2IdyImka4awrwBbCnpH6EC/61aVbzM0mzo239kHCK88do3iTCNdWRki4h1Kauoeov8HzCNbdDJI0DVkWnfLluAN6TdDNwB+EO9u+AB/OcCidmZp9J+itwa1QDmkJ4T7cj3FgoepPA9ZJ+SLhL+3+E2uWoaP4E4DRJ+xNOCU8nNIPLd8xJvED4oXwmuuSCpP8Sbto9UN1KkgYRappPEb4TfQg/Iu9FizxIKJSflnQZMJHw2R5AKBhrulFpwP1R4/guwB+Ap21Dk7PLgLskLSLc8FxDKPwON7NzE8SW9DszhdAi5CJJNxG+M9dSRA3QzD6V9BTwh+j7MY9wltGxmO1EajxuM1sm6Tbgl5LmEL4rZxF+kOdWt9FieA012BuYRfhyvUb4pR4FbBfVcKojwnXUjwjXBNsRPrzKL8LPCa0AJhC+KP1qEduPCHdRxxB+WY82s+mw/prqtwl3n98n/GNdxoZLAZjZOsJdzBOi43s/307M7APgKEKtdCyhsHiacP2xrr4LPEtozjQ2Oo4jzGx8Lba1jtAJ4HZCLa838E0zqzxlv5pw7ffvhGutnYHf1iH2/xAqHvFT+xfyTMu1jFAzfohQWP6NcO3uQoCoxrRvdAz3RMs8RjhtnlIgprcJ1/2eI7yv44ARlTPN7AHC5/3NaNl3CDfeZiSMLel3Zj6h5cbBUQw3Er6v6/ItX4MRhP+hfxFuPs2Ijm1lDevki6fQcUO4ifoE0U1YYFPC/00qtOF/37mGLWoWc6eZ+ZlVGYvuA4wHnjSz/5d1PMXwL6ZzLlNR65jKs6wOhEshAwhNxBoVL1Cdc1lrDlxOaNq0hnD6v7+ZfZhpVLXgp/zOOZcSvynlnHMp8QLVuSIkyZ1Z6ly2efbX6HKklisvUJ0rTqlyqRYk6dRqkm//APhWbLk7Jb1Ub4G59fymlHNFsBLlxK2LnH79LkNeQ3UASNpL0muSlkSPsZIOjc3vGZ1azovmv5abDEbS/pI+kLQy+ru/YsNfSBoQvd4rZ70qIxtIai/pFkkzJC2X9L6kY2PzK7dzgqSnomU+l3RaznbbS/qNpGmSVikMvXFpMceU533a6JRf0kWSpkdxPEueDhySdpb07yhpxzxJj0nqn7tdSUdLGq8wFMyLCmkDUcie/0D0vHKEiXuj1+tP+aP38Sxg39hyZ0i6T9K/88T1YuV2XN15gepSyd0pqQ/wD0KC6p0I3QdvqUUsInSJ3JHQC2w7QlfbhyQdmLP4tYRCZgfgEeAeSVvGtvMPQu+viwhdEE8n9Firaz7SeLxHE4bzuJmQR+ERQjfe+DLbEnpvvQEMI3QvXRvta5PYor0J3XJPIfQv35QNeU7X92RiQ37dfHljbyT08nsjttzDhAQ+B0naPBbXQEJvrT8lPV5XQBo5AP3RuB+kkLsT+CWhy2SL2PwjiOXCZEOG+71ytrM+bywhx+tKYvk+o+l3A0/kbOeHsfktCKfi50avD4yWGVbbY6pmvZFUzaX6KiHfQXyZG6ma2/Ze4KGcZVoTfrCOiW33K2KjHxCyz68DNolenxr+ZTeK6V6q5g69E3gpz3IfEOXSjV7/ChiX9fevnB5+DdVhZgsVRut8VtILhNrU42Y2IVoknrszvmprNoxLtS0hQXV8JNXa5JnchZD6bUbOvloRMj/Frc+PaWZfRQkvKnOz7gwsNLN8eTUr91PomJLYFvhLzrRXCTX0+L4GKeTnjNuEqpmwZlrV8ZZmEPJF9GBDntG6uJ2Qv/fn0XbPAK5LYbsu4gWqA0J+TEm3ELK/HwxcLelCM7udDbk7/yfPqpXDZuTLb5r7eqP8rZGWsefNCOP95Bt2IzcfZk35MfPtPy7JMSVVqHdMM8KliWvzzIunhsx3PJXrp+EBQgH6TTYMSHl/Stt2eIHqYszsI0K3v5ujNGfnsCGr0+mE3J3VpTkbR0ib19w2jBa6V84ylbWv9flbJfWg6nDWownXDjeJ4qmtd4EukoZVU0tNckxJfEzInnVrbFpuvs3RhOu8n1l0rl1LqyFc87aaR2TNm4fXzCokPURI3dcM+Jv5SK2p8ptSDkmDFIZw2UtSf0m7E1IaVib/fZCQj/VpSYdEd9l3k/S/kioHU/sjIdnwHZIGRzeQ/i++HwtDV78G/ETSjpJ2JtSQ4kOqvEAY1fMxSf8jaYvoDvlFkorJZP8CYXjwh6M755tL2lPSd4s4piRuAr4t6QeStpQ0AjgtZ5lrCDfF/ixp1yiW/aOWDFsUsa8vor9HSeouqX0Ny20jaYikbpJax+bdTrgJdygh761LkReoDlLI3Wkhu/2R0bQxhDv8P8yzrzMJN49ej/Z3ByEXLdF2jHBn/jHCnfPxhLys3yQM651ItJ1vEu7c30bISftnQhLuRMeUcD+PE66X/oRw0+cU4Kc5y3xCuGvfnpC/9GPCnfU2wKIi9vUO4X29jTDc8u+rWfQuQi7Q1wlnBSflbONDQm35v0n37ZLx5CiupBR69pxmZknH2nIlpDC8yBTgZjO7Ket4yo1fQ3WuCZDUjNBa4FxCTfnObCMqT16gOtc09CNcW50FjDDvrloSfsrvnHMp8ZtSzjmXEi9QnXMuJV6gFkHSS1EXTZdHub4/DeG48mW5SmMd1XMy7Gif6zOQ5XvdmHmBmgHlpKtrqiT1jf6Z9ityvd9LWivp+0Wss6OkFxTS9C2U9LykQQlXP5b8bWpr2t/IWPq8tZIWSXpX0vWSNitmWyX2OiEj1UyokhpxQG03qJAm8fIaFukNPFrb7efZXy+FlJGzJbUsvEbpeIHqGhVJbQlZl64hdI1N6hFCgpW9gN0I2atycwrkZWYLzKyiyFABJhMKj77RPm8gdCYYJ2mPWmwvdWa22sxmm9m6wkunts/ZUceKtJxJ6PzxJXB0itstmheoxWsm6VpJ8yVVKAw30Sa+QNRNcnz0q/mppMuiBtUoDE0xEPh5rAYzQNKrkn4Z28YvonkHxab9V9L1sdcHKyRFXqGQjPkeSV1zYjlR0pgolsmSbpbULjb/pegYroh+4RcoJCxuRw0Uuqg+E+17qqSL8ixzsqS3JC2O3q+nJW0VW2Ra9PfF6Fgn17TPyLeBzwnpAvsUUTB9BTxpZmPNbKKZjTKz3OxVeeWe8hfxnq2NCo9ZZjbBzB4iFOjjgPuitqGV2yz4WebEVJmQ+mSF5Noro1r35nmWzZu0OppX4ym/pJbRd2a6QpLuWQr5AGpN+U/xu0r6WxTjTEmJzgii9/Bs4L7oUcyPbOq8QC3e8UBXQl/3UwjdJNenQFM4lf8R8L+E/ts/IDSm/nm0yLGEmstNbEgAPI3Q9zyeQPkAQrfBA6PttiGMZfRC9PoA4O+E7ps7AMcQ8oQ+LoV8dJLOIPSxv4mQZu504CBC18XcY+pCyEV6crStn1T3BkTbfzx6H/aL3oOjCIml41oDV0fTDyYkVX5aUqtofuXyx0XvQ74MU7nOBe41s1XRsSf9B3oM+LGkHRIuX0hR71klM1tD+DwGER1/ks+yGr2BCwg/MnsTkmU/kbNOTUmrk7gIOIFwVrAl4XN+s4j1k/o58BIh0fd1wPWKjdJQg0OAdoQuxg8A+6m4/Ajpyjoha2N6ED7wyUDz2LRzCMk92gFtCanfDstZ73RgUez1+oTKsWn7EWpRHaPtrCL0EX87mn8wIYtQu1gs1+Zsox8h5dvQ6PVk4LycZfaJlukc284HOcvcBrxRw/twULSNrWLTuhPyiN5Zw3pdovX2jF73pYbE1nnW3zF6D7pHr3eN3u9NC6x3MrAEuITQsP2A2Lwtoxh2LPC535nzusb3jJxE1DnLbhPt84QiPssq24teGzAoNm2raNpBsWUKJa3ej1gy7Dyx3kL4EVcR/yeTgctrmL8+6Xjs9QM5y4wCXk2wr8eBX8de/xO4JmmsaT+8hlq8t61q6rTXCNfmBgJDCAkv/qYwdtBShaTCtwOdJHWvYbtvEAqLfQi1jSmETExfl9SJUGN928yWRcvvAlycs5/K7FBbRvvqT0jFF1/mX9Ey8RsyY6hqBhsSNeezLTDfzCZWTrCQGHlCfCFJQyU9LukLSUvYkCS5PzWIxyvpX7FZ5wL/jPaFmb1N6P1T7R3i6JTwRsI/2a+B7xBqfpXr7AgsIJyGF6PY96xKWNHfyl41NX6WNWxnnpmtv4sffR7zCZ9PpZqSVidxD7A9MEnSbZKOi51hpOmNnNevER2HpEtzvhN7R9N7E0aFuC+23r3ACEWX2Oqbdz2tu/jpVeUP1LcI2YtyVZt70sxWSXqdcIq/GnjBzOZJGk+oRRxAyFQU39d1RAO35ZhNqDFDuOTwYp5lpseeF0rUnCtfMumqC4SbR/8mZK8/M4oJQsFV6B9yaOz5imh77QinrR0kxUcFaEY4S6gu81IPwmnvuwBm9m9JJxLSA/YnvN93W9WRBpIo9j2L2y76W5k9q9BnWYzcSwR1SlptZmOi67IHA/sTaqxXSxputbtRl1T8OG4j3FSsNCP6exahDBudc2WkORsyltUrL1CLt4uqJvjdnfCl/YzwJVgJbGFm/6xhG3kTABNOrb4dzb8+Nu1/CEN6xNPCjQaGxGsoOZZKmgZsbWZpD8I2DuguaUuLbuxI6kY45axM5jyYcBngMgvp64huIMW/+ZX/7FXei2qO6UTCNdgdqVqYdwJejv7B813bW0C4LLAPoYDHzP4l6VuEf7gVhGu49UKhWc8PCT+4lbXcQp9ldbpLGmhmn0Xb3opwXfuTtOKF9UNnP06o2V9DuGyyL2EwxbQMp2qS7t2JjsNCEuwqlZHozOO7hNYeuUPQ/JTwI+sFaiPQFfiDwnAhWxBuuvyp8lQ8+sJdE/1iPkd4j7cHvm5mlQXiF8CekvoR/tkXWGi28gLh7vVaNtQqXyC02VtD1dOiK4F/S/o14ZRnCeH08FvAhRaSOV8G3CVpEfBEtI3BwOFmdm4d3oP/AGMJCZMvIhSM1xGu11WaQrgOfJGkmwg3Wa6lamE4n5Ab9RBJ44BVZrawmn2eSxjn6sPcGZJeI/wDbVSgmtlqSTcQxlKqLBg6AIcR3o/2hJpvdTXcumguqVf0vBPhhsslhGuoh9qGpkpJPst8lhNGer2E8EP1O0Ku0+fTOgBJPya0UR0T7e8kwvcz3xlYXC9JQ3OmzTez6XmXhiMkXUg4CzuMULE4sYbtH0a4zny7mVUZb0vSPYQRZQeY2eQCcaYrq4u3jfFBuHlwN6E94ZeEL/7dQNuc5c4ifAFXAgsJwzOfH5s/jHAKuoJQwAyIpjcnjKc0NrbspoSC6vk88exN+OdZQkgS/QnwG6qOPHoMoSBeDlREcV2Zc0x35mz3cmBygfdiAKHGt5Jw+eAHudsi3An/NFrmfUKt5ivgjNgypxN+YNZUt0/CJQAjFEL55n8vOv5O1cwXIYv+u9Fy8wm1msHAxYQC4rgCn3vuTaka3zM23DQywk2gxcB7hDOPjW4AFfosyX9TahLh+vFkwo/XC8DA3GVy9rNXznduP2q+KXVu9L5VEH783gGOLvDdmBw79vjjtmh+vptSFxN+9JcTasA/LrCPv1PNjVPC/9FsYiO81tfDs0051whFzfNONbOkvb1cPfC7/M45lxIvUJ1zLiV+yu+ccynxGqpzzqWk3ppNKSQFGc6GpjUzzGzraN7JwK8IQ/w+B5xpoe0ZkroQhsU9hHB39n/NbFSh/XXr1s0GDBiQ8lGkZ/HiMKRPp06dMo7EOVeMd999d76Z5e31WN/tUC80syqJeiUNIXTN/CahWckdhAa+lW3Q/kBo59iT0HzmaUljzazGroIDBgxg9OjRNS2SqeefD00FDzrooAJLOucaEklTqpvXEBr2nwI8ZWYvA0i6AvhEUgdC+73jgO0s9NZ4VdKThDaFP8sq4DR4Qepc+anva6i/UsiL+Zo2ZGkfQuh1A4CFbnSrCd0YtyLklIz3yhgbreOccw1KfdZQf0rIoLOacDr/VNQ1rT2hF0ncYkL3wLU1zNuIpHOI8mP269cvtcBL4e9//zsARx+daYJx51yK6q2GamZvmdkSM1tlZvcR0nN9g9CdrWPO4h0JXfBqmpdvH3eY2TAzG9a9e02Z8rLXsWNHOnbMPTTnXGOW5TVUI/SxHkfIIARAlG27NSH5wjqgRTyrUbRssbkrG5z9998/6xCccymrlxqqpE0lHSppE0ktJJ1CSKf2LPAgcKSkvaOcl1cBj0W12WWEFFxXSWonaU/CIFz58kY651ym6quG2pKQlm4bwnXR8cAxZjYBQNJ5hIK1KyHjzojYuhcQMjrNJWR4Or9Qk6nG4LHHQqrGY49NMmyOc64xqJcC1cIQDNUOwBY11M/bWD9q4H9MiULLTNeu1Q5o6ZxrpBpCO9Qmad999806BOdcyrwvv3POpaTGGmo00NqJhG6hOxKyxy8iNK7/F/AX2zAKpyvCo48+CsDxxx+fcSTOubRUW6BK+ilhWItXCGMI/Z7Q/rMDYeiIg4FfSPqdmV1bD7E2aKPemrrRtJN3q75zQa9evaqd55xrnGqqobYCtrH8Q8W+CNwajRd/UUkiK3N77bVX1iE451JWbYFqZlcXWtnMFhOaQznnXJNX0yn/Tkk2YGbvpRdO0/HII48AcMIJJ2QciXMuLTWd8v8nwfoGdEkplialb9++WYfgnEtZTaf8neszkKZmjz32yDoE51zKvB2qc86lJFFPKUltgB8D+xLGfVLlPDPboTShlbe//OUvAJx00kkZR+KcS0vSGurvgCOB+4GBwE2EwfYKDpbn8tt8883ZfPPN63Wfkpg0aRIA5513HldfXbAhh3OuCEkL1COAI6PE0F9Ff4+jDJOW1Jfhw4czfPjwvPMGDBhAq1atmD9/fpXpQ4cORRKTJ0+u8/5vu+02rrjiijpvJ02jRo2if//+tGvXjmOOOYYFCxbkXW7ixIkcffTRdO/enS5dunDooYcyYcKEvMsecMABSOKrr77KO9+5NCUtUFsR0ucBLIsG0JsMbFuKoFyowVZeFgD48MMPWbFiRYYRlda4ceM499xzeeCBB5gzZw5t27blggsuyLvsokWLOOqoo5gwYQJz5sxh1113zTuUzIMPPugFqatXSQvUj4DK29JvAjcA1xMKVVcLDz74IA8++GC180877TTuv//+9a/vu+8+Tj/99CrLrFq1ih/96Ef069ePnj17ct5551UpdG+44QZ69+5Nnz59uPvuu6use8YZZ3D55ZcDsHDhQo444gi6d+9O586dOeKII5g+ffr6Zffbbz+uuOIK9txzTzp06MAhhxyyUe25rh588EGOPPJI9tlnH9q3b8/VV1/NY489xpIlG492s+uuu3LWWWfRpUsXWrZsySWXXMKECRP48ssv1y+zePFifvGLX3D99denGqdzNUlaoF4ALIyeX0y4MbUDVRNBuyJstdVWbLXVVtXOHz58OBUVFXzyySesXbuWhx9+mFNPPbXKMj/96U+ZOHEiY8aMYdKkScyYMYOrrroKgGeeeYYbb7yR5557jk8//ZTnn3++2n2tW7eOESNGMGXKFKZOnUqbNm248MILqywzatQo7rnnHubOncvq1au58cYb825r6tSpbLrpptU+Ro3Kf9l93Lhx7Ljj+pFwGDhwIK1atWLixIl5l497+eWX6dWrV5Ucs5deeinnn3++50xw9SppPtQ1ldn1zWwacDyApK1LFVi522WXavNtr1dZS913333ZZptt+NrXvrZ+npnxpz/9iQ8++IAuXULfiksvvZSTTz6ZX/3qVzzyyCOMGDGC7bbbDoCRI0dWuYQQ17VrV4477rj1ry+77LKNxrwaMWLE+h+AE044gSeffDLvtvr168eiRYsKHluupUuX0qlTpyrTOnXqlLeGGjd9+nS+973vcfPNN6+fNnr0aF577TVuueWWKjVt50otaYH6DhuPPgrwBt5TqmROO+009tlnH7744ouNTvfnzZvH8uXL2XnnnddPMzPWrl0LwMyZM6vM69+/f7X7Wb58OZdccgnPPPMMCxeGE5ElS5awdu1amjdvDlTNjtW2bVuWLl1a9wOMad++PRUVVfPwVFRU0KFD3hHDgfAeHHLIIVxwwQXrm5+tW7eOCy64gFtuuYUWLTx/uqtfSU/5tdEEqQdhfChXC/fff3+Va6T59O/fn80335x//vOfG4091a1bN9q0acO4ceNYtGgRixYtYvHixesLut69ezNt2rT1y0+dunF6wUo33XQTEyZM4K233qKiooKXX34ZCAV0saZOnUr79u2rfVR33XjIkCGMHTt2/evPP/+cVatWVXtZZOHChRxyyCEcddRRXHbZZeunV1RUMHr0aB7Hr4UAACAASURBVL797W/Tq1ev9WcCffv25ZVXXin6eJwrRqEE00sI/fXbSspN49cW+G2pAit3Q4YMSbTcXXfdxcKFC2nXrl2VO9bNmjXj7LPP5pJLLuH3v/89PXr0YMaMGXz00UcceuihnHDCCYwYMYLTTz+dAQMG8Itf/KLafSxZsoQ2bdqw6aabsmDBghqXLaRfv361qr2ecsop7L777rzyyivstNNOXHnllRx77LF5a6gVFRUceuih7Lnnnlx7bdVUvJ06dWLmzJnrX0+bNo1dd92Vd999l+7duxd/QM4VoVANdS/CcM9Lgb1jj72Ar5nZD0sbXvnaeeedq5ySV2fgwIEMGzYs77zrrruOQYMGMXz4cDp27MhBBx20vj3m4YcfzsUXX8wBBxzAoEGDOOCAA6rdx8UXX8yKFSvo1q0bw4cP57DDDqvdQdXBkCFDuO222zjllFPo0aMHS5Ys4dZbb10///DDD+eaa64B4PHHH+edd97hnnvuqVL7nTp1KpLo1avX+kdlIdqzZ09atWpV78flmhYVe1onqY2ZNfgGkcOGDbPRo0fX2/6KzdjvnGucJL1rZnlrOYmuoUpqI+kWSQuApZIWRK/bphppE3Lvvfdy7733Zh2Gcy5FSW9K/RrYBjgM6BP93Qq4uaaVXPWGDh3K0KFDsw7DOZeipO1KjqLq+FJzJJ0IjAfOK0lkZc4LU+fKT9Ia6lqgXc60dnizqVpbu3bt+jajzrnykLSG+gDwjKRrgKlAf+Cn0XRXCw88EN66M844I9tAnHOpSVqgXgHMJySZ7gPMBP4M3FKiuMreTjslGgPROdeIFGrYf7GZ/cbM1hJuQPlNqJTssIMPdOBcuSl0DfWqeomiCVqzZg1r1qzJOgznXIoKFagb9eF36SiUD9U51/gUuobaQtJJ1FCwmpmPK1UL1XUndc41XoUK1JbA92uYb/hAfbVSmafUOVc+ChWoy81s93qJpIlZuXIlAJtssknGkTjn0pK0Yb9L2UMPPcRDDz2UdRjOuRQVqqHOLTDf1dJuu+2WdQjOuZTVWKCa2aD6CqSpGTx4cNYhOOdS5qf8GVm+fDnLly/POgznXIq8QM3II488wiOPPJJ1GM65FPmwkBnZfXdvPOFcuUlUoErqDWxHGJhvOjDezJaVMrByt/XWW2cdgnMuZYWSo/QA7gcOiSatA1YDqyXdYGb/V+L4ylblyKDt27fPOBLnXFoKXUO9k1AjHQhsSch/ehlhJNTDJf28tOGVr0cffZRHH3006zCccykqdMq/L9DNzNYASPoeMMnMfi3pNOAVoPaDuDdhe+21V9YhOOdSVqiGOp+Qnb9Sf2A5gJl9AXQqdoeStpS0UtKfY9NOljRF0jJJT0jqEpvXRdLj0bwpkk4udp8Nyai3pjLqram8/WUr3v6yVd7hp51zjVOhAvUm4CVJN0i6EfgP8DsASdsAX9Rin38A3ql8IWkIcDtwGtCTUGDfmrP86mjeKcAfo3UatTUrlrJmxdKsw3DOpahQT6lbJX0GHEFI4Xe2mT0dzZ4KHFDMzqKRUhcBrwOVvbBOAZ4ys5ejZa4APpHUgXAT7DhgOzNbCrwq6UlC4fuzYvbd0Mz54BUA+u52eMaROOfSUrDZlJk9CzybZ/pyotP/JCR1JIwAcCBwVmzWEEIBW7ndzyStBrYiFKhrzWxibPmxhGu7+fZxDnAOQL9+/ZKGlokuA3fMOgTnXMpq1VNK0lRJ3Ypc7WrgLjObljO9PbA4Z9pioEOBeRsxszvMbJiZDevevXuR4dWvtt360LZbn6zDcM6lqFA71H9XM6sH8Kik1WZ2SDXLxLczFDgI+Hqe2UuBjjnTOgJLCDXU6uY1amuWh0No2Tbvb4NzrhEqdMq/C+EUexQbhkERsCvwNFCRcD/7AQOAqZIg1DybS9oWeAZYf/4raQugNTCRUKC2kLSlmX0aLbIjMC7hfhusOR++Cvg1VOfKSaECdSvCnf6TgAvNbBxA1KD/PjNLmi/1DiCeTflHhAL2fEJt9w1JewPvEa6zPmZmS6J9PQZcJem7wFDgaGCPhPttsLpsma+y7pxrzArd5Z8HnC5pP2CUpOeAkYSxpBLLvYElaSmwMtr+PEnnAQ8CXYHngRGx1S8A7iYku/4SOL+yYG/M2nbplXUIzrmUJUqOYmYvSdqZULMcQy0a9Odsb2TO61FUM9ifmS0AjqnL/hqi1UvDvbZW7ev0VjrnGpDEd/nN7Cszu5bQ7OlMQntSV0tzx73O3HGvF17QOddo1CYf6izgCTNbnXYwTUnXrXbKOgTnXMpqrKFK+lHseWdJTwDLgKWS/iGpa6kDLFdtOvekTeeeWYfhnEtRoVP+K2PPryM0mRoUPVYD15corrK3aslCVi1ZmHUYzrkUFTrlV+z5ocBuZjYbILoz/16pAit38z5+E/B2qM6Vk0IFarx5VDNgzvoZZnMl+S3qWuq2zbCsQ3DOpaxQgdo+SlQC0BzYgdBzCkkDCflSXS1s0qlh5xpwzhWvUIE6OOf1rNjz7oAPgVJLqyq+BKB1R7+v51y5KNRTakIN894E3kw9oiZi3idvA34N1blyUqv0fQCSmkn6SZrBNCXdB+9K98G7Zh2Gcy5FtS5QgZbAr9IKpKlp3bGrn+47V2YK5UO9tYbZtell5SIrF88D/OaUc+WkUA31LKAVoXdU7sNHmKuD+eNHM3/86KzDcM6lqFAt8yNCv/1/5M6QtAlwcUmiagK6bzs86xCccykrVKA+QMien88aQndUVwutO3TOOgTnXMoKNZv6TQ3z1gL/m3pETcSKhaHTmSdIca581OUuv6uDLye+x5cTPRWCc+Wk2gJV0muSjpTUvJr5zSUdJemV0oVXvnoM2YMeQxr90FjOuZiaTvm/B/wSuEvSm8AEwvDNHQiD9+0GvAVcWOogy5EPfeJc+am2QDWzMcARkvoRUvdtD/QEFgJPAt8zs+n1EmUZWr5gNuCD9TlXTgo2zjezqcCf6iGWJmXBp+8D0Nb78jtXNry3U0Z6br9X1iE451LmBWpGWrbtkHUIzrmUebOpjCyfP5Pl82dmHYZzLkVeQ83Igs/GAtC2W5+MI3HOpSVRgSrpOOAjM5sQDX1yK7AWuMjMPitlgOWq5w57Zx2Ccy5lSU/5rwMqouc3ARMJY0vdVoqgmoKWbdrTsk37rMNwzqUo6Sl/DzObJak1sC/Qm5AcZV7JIitzy+aFJrztuvfNOBLnXFqSFqhfShpAaNz/rpmtlNQGv6lVaws//xDwAtW5cpK0QL0GeB8w4JRo2v7Ah6UIqinoNXTfrENwzqUsUYFqZn+S9Nfo+aJo8hjg5FIFVu5atG6bdQjOuZQVc8puwMGSfhC9/opwHdXVwtK5U1k6d2rWYTjnUpSoQJW0B/ApGzJQAWyH3+WvtUVfjGPRF+OyDsM5l6Kk11B/C5xuZs9IWhhNewMYVZqwyl/vr++fdQjOuZQlLVAHmtkz0XOL/q4ijIjqaqF5q02yDsE5l7Kk11AnSMqtUu0H+DlrLS2dPZmlsydnHYZzLkVJa6g/Bh6X9DdgE0m3AN8CjitZZGVu0ZRPAGjfa0C2gTjnUpO02dQrkoYBpwN/IWTt39PMvihlcOWs904HZh2Ccy5lSZOjNAOmmNlVsWmS1MzM1pUsujLWvKVffnau3CS9hvofIHeIzj2A59INp+lYMutzlsz6POswnHMpSnoNdSjwes60N4CvpxtO07F46gQAOvTeIuNInHNpSVqgLgG6UTW7VDdgReoRNRF9hh2cdQjOuZQlPeV/AnhA0iBJzSVtCdwLPFayyMpcs+YtaNbcB0xwrpwkLVB/BswAPgZWAx9Fr39aorjKXsXMz6iY6YMdOFdOEhWoZrbczM4C2gIDgHZmdraZLU+6I0l/ljRLUoWkiZK+G5t3oKTxkpZLelFS/9i81pLujtabLemHyQ+v4aqYNpGKaROzDsM5l6LE55yS2gKDgPbAZpIAMLPcm1XV+RVwlpmtkrQN8JKk94EphEsH3wWeAq4GHgaGR+uNBLYE+gO9gBclfRzrCtsofW2XQ7MOwTmXsqTtUE8mZJYyIF4rNSDRsJ1mFu+matFjILAzMM7M/hrtayQwX9I2Zjae0JlghJktBBZK+hNwBtCoC1Q188EOnCs3xQzSd7qZdTKz3rFHUWMgS7pV0nJgPDAL+CcwhDDgHwBmtgz4DBgiqTOhwB4b28zYaJ182z9H0mhJo+fNa9jDXVVM/5SK6Z9mHYZzLkVJC9TWhNPxOjGzC4AOwN6E0/xVhEsIi3MWXRwt1z72Ondevu3fYWbDzGxY9+7d6xpuSVXMmETFjElZh+GcS1HSAvXXhAQpdWZma83sVaAvcD6wFOiYs1hHQtvXpbHXufMatb67HU7f3Q7POgznXIqSFqhnAldJWhjdoV//qMO+WxCuoY4DdqycKKld5fTouums+PzouacNdM41OEnv8l9Yl51I6gEcAPyD0LvqIOAkwiB/rwM3SDoOeBq4EvgguiEFcD9wuaTRQE/gbGBEXeJpCBZPC11PO222dcaROOfSkjR937N13I8RTu9vI9SKpwAXm9nfAaLC9PfAn4G3gBNj6/4c+GO0zgrgusbeZApgyeyQ+dALVOfKR9JmUy0JvaVOArqZWQ9JBxGGRrm90PpmNg+odiB6M3se2KaaeasIlxzOTBJrY9F3l8OyDsE5l7Kk11BvIqTru4AN40hNoI6XApxzrpwkvYZ6PLC1mS2RtA7AzKZJ2qx0oZW3RVPDECib9huccSTOubQkraF+lTtBUlfCUCiuFpbNncayudOyDsM5l6KkNdTHgLsk/QAg6sH0G+CRUgVW7r427JCsQ3DOpayY9H0LgMnApsBsQqP7K0sTlnPONT5Jm02tBM6T9D1C3/pZZrbRZQCX3MLJoW9C5wF50xI45xqhpM2m4klQ1gI9Yun7ZpYgrrK34stZgBeozpWTpNdQpxMa5yt6bbF5zVONqInos/NBWYfgnEtZ0gK1Tc7r3sAVwIvphuOcc41X0muoq3ImTZZ0EWGMqT+nHlUTsPCLDwHovPn2GUfinEtLXYbd7E81eUldYSsWhQTYnTOOwzmXnqQ3pZ6j6nXTtsDXgetLEVRT0OfrB2QdgnMuZUlrqI/mvF4GjDGzj1KOxznnGq2k11ALZpRyxVnw+QcAdNlih4wjcc6lJekp/6VJljOza+oWTtOxqmJB1iE451KW9JR/Z+Ao4ANCoud+wA7Ak8CaaBnLv6rLp/fQ/bIOwTmXsqQF6nLgO2Y2qnKCpJOAw83s9JJE5pxzjUzS5ChHAQ/nTHsEODrdcJqOLyeN4ctJY7IOwzmXoqQF6hfAd3OmnUXIPuVqYc2yxaxZtjjrMJxzKUp6yn8O8LiknxD69fcFNgH+p1SBlbteO1Y7xJZzrpFK2mzqbUlbAPsQpe8DXo7S+jnnnCP5KT+Eu/mLgMVm9m8ASa1LElUT8OWn7/Hlp+9lHYZzLkVJ26EOBp4AWgNdCX34DyYMK31yyaIrY2tWLss6BOdcypJeQ70NuMHM7pRUOTDfi8AfSxNW+eu1/d5Zh+CcS1nSU/4dgLui5wZgZkuBdqUIyjnnGqOkBepUYMf4BEk7A5+lHlETMX/CaOZPGJ11GM65FCU95R8JPC3p90BLSZcAFwHfL1Vg5W7tmtyc3c65xi5ps6nHJc0itEd9m1BbPdXMXi9lcOWs53Z7Zh2Ccy5liTP2m9mbwJsljMU55xq1GgtUSQcAy8zsreh1f+BuYDvgDeBsM5tX8ijL0LzxbwPQfZtd884f9dbUjaadvFu/ksbknKubQjelfgW0jL2+i3CX/wLC8NE3liiusmfr1mLr1mYdhnMuRYVO+QcBlbXT7sC+wBZmNk3SK8D7JY6vbPXYdvesQ3DOpaxQDVXAuuj57sBUM5sGYGZz8VFPnXNuvUIF6nvAeVGf/TOAZytnSNoc8HE8amneJ28x75O3sg7DOZeiQqf8PwSeBn5LGPoknnPuVODlEsXlnHONTo0Fqpl9IKkfIWXfbDOL30W5HfD0fbXUffBuWYfgnEtZwXaoZmbAjDzT55YkIueca6QSN+xvqkrVHnTux28AfrffuXLiBWpG1Kx51iE451LmBWpGqush5ZxrvBKl75M0s5rpG58PO+dcE5U0H+pGDfgltQDapxtO0zHno9eY89FrWYfhnEtRoeQozxH67reW9O+c2ZsRdUt1xWve0sc3dK7cFLqG+mj0d2/gb7HpBswBnkmyk6in1a3AQUAXYBJwqZn9K5p/IPAHoB+hkD7DzKbE1v0jcDywHLjezG5Ost+GrNvWw7IOwTmXskIN+28HkPSWmY2p436mEXpaTQW+ATwiaXtgKfAY8F3gKeBq4GFgeLTuSGBLoD/QC3hR0sdmlqgwd865+pI0Y/8YSacRho3ubma7SNoT6GFmjydYfxmhYKz0D0lfADsThqUeZ2Z/BZA0EpgvaRszGw+cDowws4XAQkl/IuQVaNQF6uwPXwF89FPnyknSu/xXAhcDjxBqixBO+S+rzU4l9QS2AsYBQ4CxlfOiwvczYIikzoRur2Njq4+N1mnUWm7Sjpab+KCxzpWTpO1QzwJ2NbM5kn4dTfsM2KLYHUpqCTwI3Gdm4yW1B3Kz/i8mtCxoH3udOy/fts8hjHtFv34NO7t91y13yjoE51zKkjabagUsip5b9Lcd4SZRYpKaAQ8Aq4ELo8lLgY45i3YElkTzyJlfOW8jZnaHmQ0zs2Hdu3cvJjTnnKuzpAXqv4HrJMX7S14J/DPpjiSJMIRKT+A4M1sTzRpHGEW1crl2wEDCddWFwKz4/Oj5uKT7bahmj/0vs8f+N+swnHMpSlqgXky45rkI6ChpIbA98JMi9vVHYDBwpJmtiE1/HNhO0nGSNiEU1B9EN6QA7gcul9RZ0jbA2cC9Rey3QWrZrhMt23XKOgznXIqS3uVfCHwjGvW0HzDNzCYn3Um03rnAKmB2qKwCcK6ZPSjpOOD3wJ8J7VBPjK3+c0JhPAVYAVxXDk2mug4amnUIzrmUFZUcJWpsP6XYnUTrqYb5zwPbVDNvFXBm9HDOuQYrabOpNZJW53kskfSJpP+T1LbUwZaTWWNeYtaYl7IOwzmXoqQ11B8Sun5eT+jx1A/4f4SeTZ8DlxMa6J9XghjLUuuOXbIOwTmXsqQF6vcJ7VAXRq8/kPQG8JaZDZL0HuHapxeoCXXZYoesQ3DOpSzpXf7OQMucaS2i6QDTCe1SnXOuyUpaQx0FPCvpZsIpf1/gkmg6wIHAp+mHV75mvv8CAH2+fkDGkTjn0pK0QL0EuIhwSt+H0Nj+AUJTJ4A3gP1Tj66MtdnUe3I5V24KFqhR76ifAjeZ2W/yLWNmS/NNd9XrvPn2WYfgnEtZwWuoZrYW+BGh/71zzrlqJL0pNQoYUcpAmpqZ7z7PzHefzzoM51yKkl5DHQycI+knhJtSlRmnMLNDShFYuWvTtXfWITjnUpa0QH0keriUdB7Q6HNkO+dyJE2OcnupA3HOucYucXIUSV0IY0B1I5boxMxGVbuSq9aM0WFU7q8N8ysmzpWLRAWqpG8ADwEzCMmfPwMGAe+woXG/K0K7HptlHYJzLmVJa6jXAudHuUsXmtlgSecRhnZ2tbBpv8FZh+CcS1nSZlMDzOzB6HnlHf478Rylzjm3XtICdb6kHtHzaZJ2ATZj44QpLqHp7zzD9Hca/cADzrmYpKf89wD7An8Ffgu8BKwFbitNWOWvQ6/Nsw7BOZeypM2mro49v0vSy0B7M3u/ZJGVuU6bbZ11CM65lNVYoEqqMLOOudPNzFP1OedcjkLXUKsdWM/VzfS3/sX0t/6VdRjOuRQVOuW3AvNdLXX82qCsQ3DOpaxQgdpO0tSaFjCzfinG02R07Ltl1iE451JWqEBdBZxWH4E0NbZuHQBqlrTlmnOuoStUoH5lZv+tl0iamBnvPAtA390OzzgS51xaChWoflOqRDputlXWITjnUlaoQD23XqJogjr2GZh1CM65lNV4Ac9T85XOurVfsW7tV1mH4ZxLkd8RycjM0c8xc/RzWYfhnEtR4gTTLl2d+nnXU+fKTa1rqJI6S7ogzWCakg69t6BD7y2yDsM5l6KiClRJzSUdJelRYBZwfmnCKn9r16xm7ZrVWYfhnEtR0iFQdgJOB04C2gCtgePN7KkSxlbWZr33H8DboTpXTmqsoUr6kaQPgdeBLYAfAL2ABcBbpQ+vfG3afzCb9vdhUJwrJ4VqqNcDXxJqp381MwOQvL1/XbXvNSDrEJxzKSt0DfUA4CnC+FHTJd0kaWc8C1WdrV29krWrV2YdhnMuRYUa9r9kZmcCPYGfATsQTvV7AudK6lr6EMvTrPdfZNb7L2YdhnMuRUmHQFkBPAA8IGkzQgaq04D/BdqWLrzytenmQ7IOwTmXsqLboZrZNDO7xswGA/uXIKYmoX2PfrTv4alknSsnRReo0V1/AMzM7/TX0lerlvPVquVZh+GcS1Ftekr1Tz2KJmj2mP8ye4ynmnWunNSmL7+3mUpB5y22zzoE51zKalOgeteeFLTr3jfrEJxzKUt0yi/p/crnZvZqbProUgTVFKxZsZQ1K5ZmHYZzLkVJr6FuNOaxQnepxOmSJF0oabSkVZLuzZl3oKTxkpZLelFS/9i81pLullQhabakHybdZ0M254NXmPPBK1mH4ZxLUY2n/JLuj562ij2vNAAYV8S+ZgK/BA4lJFip3Ec34DHgu4ReWVcDDwPDo0VGAlsSbob1Al6U9LGZPVPEvhucLgN3zDoE51zKCl1D/aya5wa8Bvw16Y7M7DEAScOA+AXEY4FxZvbXaP5IYL6kbcxsPCGPwAgzWwgslPQn4AygUReobbv1SWU7o96autG0k3fz9q3OZaHGAtXMfgEg6U0ze7ZEMQwBxsb2uUzSZ8AQSXOAPvH50fNj8m1I0jnAOQD9+jXsQmXN8iUAtGzbIeNInHNpSdr19FlJrYCtgW7Emk6Z2Qt1jKE9MC9n2mKgQzSv8nXuvHxx3gHcATBs2LAGncBlzofh3p7nQ3WufCRNML0X4fS+NdARqCAUatMo4sZUNZZG24zrCCyJ5lW+Xpkzr1HrsuXXsw7BOZeypHf5fw1cb2ZdgCXR36uBW1OIYRyw/g6NpHbAQMJ11YWEoVbid3B2pLibYQ1S2y69aNulV9ZhOOdSlLRA3Qq4JWfatcAlSXckqYWkTYDmQHNJm0hqATwObCfpuGj+lcAH0Q0pgPuBy6NBAbcBzgbuTbrfhmr10sWsXrq48ILOuUYjaYG6mA2n5bMkbQt0ZsM1ziQuB1YQ8qqeGj2/3MzmAccB/wcsBHYDToyt93NCC4MpwH+BGxp7kymAueNeZ+6417MOwzmXoqRdTx8DvgGMAu4CXgTWUFyzqZGENqX55j0PbFPNvFXAmdGjbHTdaqesQ3DOpSzpXf6LY89vkvQ2oXZaqqZUZa9N555Zh+CcS1ltkqNgZt5nso5WLVkIQOsOnTOOxDmXlmoLVEmvkGAwPjPbJ9WImoh5H78JeDtU58pJTTXUO2PPBxKuYd5HuDnUD/gOcHfpQitv3bYZlnUIzrmUVVugmtl9lc8lvQkcambjYtNGEQrUn5c0wjK1SafuWYfgnEtZ0mZTg6maHAXgC6q5M+8KW1XxJasqvsw6DOdcipIWqP8F7pW0paQ2krYiNJ/ym1O1NO+Tt5n3ydtZh+GcS1HSAvWM6O84Qv/6DwkJUkaUIKYmofvgXek+eNesw3DOpShpO9QFwImSmgHdgXlmtq6kkZW51h27Zh2Ccy5lRbVDjQrROSWKpUlZuThkLPSbU86Vj6Sn/C5l88ePZv54H+PQuXJSq55Sru66bzu88ELOuUal2hqqpBtizw+on3CajtYdOnu3U+fKTE2n/OfEnj9R6kCamhUL57BioV+Odq6c1HTKP1bSo8DHQGtJV+VbyMyuLElkZe7Lie8B3pffuXJSU4F6PKGW2p/Q5nSzPMs06IHwGrIeQ/bIOgTnXMpq6ss/F/glhOFLzMwb8aeoVftOWYfgnEtZ0ob9IyR1Bo4EvgbMAP4RNfh3tbB8wWwAH6jPuTKSdBjp3YGngfGE9H1HAL+R9E0ze6OE8ZWtBZ++D0DberqGOuqtqRtNO3m3fvWyb+eaiqTtUH8DXGBmD1VOkPRt4LfALqUIrNz13H6vrENwzqWsmGGkH8mZ9igwKN1wmo6WbTvQsm2HrMNwzqUoaYH6KVWHdgb4FhvnSHUJLZ8/k+XzZ2YdhnMuRUlP+S8G/iHp+4RrqAOALQnXUl0tLPhsLABtu/XJOBLnXFqS3uV/XdJA4JtAH+Ap4J9+l7/2eu6wd9YhOOdSljg5ipktBP5cwlialJZt2mcdgnMuZZ6+LyPL5k1n2bzpWYfhnEuRp+/LyMLPPwSgXfe+GUfinEuLF6gZ6TV036xDcM6lLPEpv6T+pQykqWnRui0tWrfNOgznXIqKuYb6PkDUdMrV0dK5U1k6d+PuoM65xqvGU35J7wLvEgrT5tHkkYQup64OFn0xDoD2Pbw/vXPlotA11OOBnYFhQFtJ7xGSTe8PvGdmi0sdYLnq/fX9sw7BOZeyQqf8zczsUTP7GbAEOJqQbPoiYIykT0sdYLlq3moTmrfaJOswnHMpKlRDHSWpH2EYlE2AzsBKMzsWQFKXEsdXtpbOngxA+14DMo3DOZeeGgtUM9tNUgtge+BV4PdAB0l/BN6LHt79tBYWTfkE8ALVuXJSsB2qmX0FvC9ptZntI2kR8BLh2uq3gYNKG2J56r3TgVmH4JxLWTEN+y+J/pqZPQw8XIJ4mozmLVtlHYJzLmWJ26Ga2b3R0y1KE0rTsmTW5yyZ9XnWYTjnUlR019Mo61RZyjfuUqksnjoBgA69s/t98nGmnEuX9+XPSJ9hB2cdgnMuZV6gZqRZ88b11ntt1rnCChgpBQAACdhJREFUPB9qRipmfkbFTB+Sy7ly0riqSWWkYtpEADr2GZhxJM65tDSKGqqkLpIel7RM0hRJJ2cdU119bZdD+douh2YdhnMuRY2lhvoHYDXQExgKPC1prJmNyzas2lOzRvFb5pwrQoMvUCW1A44DtjOzpcCrkp4ETgN+lmlwdVAxPeSV6dh3y4wjSVfSm1d+k8uVowZfoAJbAWvNbGJs2ligUY8hUjFjElB+BapzDVV9/IjLzFLdYNok7Q381cx6xaadDZxiZvvlLHsOcE70cmtgQglC6gbML8F2s1JOx1NOxwJ+PA1VfzPrnm9GY6ihLgU65kzrSMjPWoWZ3QHcUcpgJI02s2Gl3Ed9KqfjKadjAT+exqgx3BmZCLSQFD833hFotDeknHPlqcEXqGa2DHgMuEpSO0l7EkYOeCDbyJxzrqoGX6BGLgDaAHOBvwDnZ9hkqqSXFDJQTsdTTscCfjyNToO/KeWcc41FY6mhOudcg+cFqnPOpcQL1BxJ8wZI+rGkjyQtkfSFpB/Xd6xJFHE8F0v6XFKFpJmSfh0N0NigFJvXQVIrSeMlTa+vGJMq4rMZKWmNpKWxR4MbOaOYz0bSTpJejo5ljqQf1GespeIF6sbieQNOAf4oaUie5QScThha+zDgQkkn1luUySU9nqeAncysI7AdoWna9+styuSSHk+lHxNuZjZExRzLw2bWPvZoiOPnJDoeSd2AZ4Dbga7AIODf9RhnyfhNqZgob8BCQt6AidG0B4AZZlZj3gBJvyW8nxeVPtJkans8kroSBmGcaGYX1EuwCRR7PJI2B/4J/BD4k5n1rc94a1LMsUgaCQwys1PrPdCEijyea4DNzOy0+o+0tLyGWlV1eQNqqgEhScDeNLzOBkUdj6STJVUQugfuSKhBNCTFfj6/Ay4FVpQ6sFoo9liOlLRA0jhJ55c+vKIVczzDgQWSXpc0V9JTksoiM44XqFW1BxbnTFsMdCiw3kjCe3lPCWKqi6KOx8xGRaf8WwG3AXNKG17REh+PpP8BWpjZ4/URWC0U89k8AgwGugNnA1dKOqm04RWtmOPpC3wH+AHQD/iC0L680fMCtarEeQMqSbqQcC31m2a2qoSx1UbRxwNgZp8Satu3liiu2kp0PNHp5/VAg7n8kkcxOSo+NrOZZrbWzF4HbgGOr4cYi1HMd20F8LiZvWNmK4FfAHtI6lTiGEvOC9SqisobIOlMQk7WA82swd1Fpm55EFoADW18lqTHsyUwAHhF0mxC1+XekmZLGlAPcSZRl8/GCDdFG5JijucDwjFUqnze0I6peGbmj9gDeIhw+tEO2JNw2jIkz3KnALOBwVnHnNLxfBfoET3flvCPcHPW8dfmeAg/Br1ij2OBmdHz5lkfQy0+m6MJrUkE7ArMAL6Tdfx1OJ4DCDewhgItgV8Dr2QdfyrvQdYBNLQH0AV4AlgGTAVOjqbvDSyNLfcFsIZwqlP5uC3r+OtwPPcQrpkuAyYDNwCbZB1/bY8nZ539gOlZx16Hz+YvwJfRd2w88P2sY6/rZwOcH/0wLCQ02dss6/jTeHizKeecS4lfQ3XOuZR4geqccynxAtU551LiBapzzqXEC1TnnEuJF6jOOZcSL1CdS4Gk/UqVc1VSvyhvaPNSbN+lxwtUVxRJe0VZghZH2Y9ek7RL1nH9//bONcSqKorjv7+PSg1NkTJzGOlFDyuIym/6RSl7QER+CHoMZdmHEMMoDAmjBgsSNKYXGWpOFJIgqUkPaKIHZgUikxqIY041RamTkRaYqw9r3Tpex+ZevXpHXT/YcPZZ+6y9z+HeddY++5y1jgeSTNKFx6Gf7ZImlepmtsM8BurfIW+TNO1YjyOpnj4XkT3pu0gaCqzGv3JZDpyGfwXT14LCJEldSA81qYaLAczsTfPIR/vM7H0z21hqIOleSZsl7Zb0nqTGgmxypCP5TVKLpI9Lnlak+WgttB0bHuGAqA+T9JqkLkk/SHq6NAWW1CTpU0nPRb8dkqYUdI2QtFie2mW3pJUF2c2SNkjqDs/7ykouhKRBkpaEvk3AtWXy0ZJWSPolxjOjIJsrabmk1+UpdL6RdE3IluEh7VbFNP/R4rWQ1IzfxFpC3iLpBUnzy/pfJWlmJeeS1JB6f/ua5cQpeDi2ncBSYAowvEx+K7AVj905AJgDfB6ykcAePOzcQOBhYD8wLeRzgdaCrrF4FKIBUV+JB7weApwNrAemh6wJj6twP9Af96B/5L+MFGvwDATDo++Jsf9qPD3K+DjuHjyOwemHOX/DI+cDPAN8gn+/3gC0E/ECcEfla+AJ3Is/H9gGXF841z+BG6PfecC6Qj/bgUn/cy3aStct6tfF+fYrXOu9wDn1/s2caqXuA8hyYpUwlkuA78MgvlP64wJrgfsKbfvFH7sRjxlbNBoKHb0aVDxH0V/AoIL8DuCj2G4CthZkg+PYUcC5wAHKjH+0ewl4qmzftyWD20P7okHdBtxQkD1QMKjjgR1lx84GFhfO9cOC7DJgX6FelUGNfZuBybH9EPBuvX8rp2LJKX9SFWa22cyazPMzjQNGAwtC3AgsjOlzN7ALN5znRbvOgh4r1nuhEfcsuwq6X8E91RI/FXTvjc0zce9xl5ntPozeWSWdobchxtobB50P8F2Z3tFleh/HbwyHjBe/6Zyho8syuxQo5Zy6E1h2FLqSIyQXpZIjxsy2SFoCTI9dnUCzmb1R3jYCDzcU6irW8ZBvgwv1UYXtTtxDHWlm+6scZicwQtJZZtbdg6zZzJqr1AnQhY+/FEC5mBOpE+gws4sOOaoyegsB15O8FWiXdBU+i1jZQ5vkGJMealIxki6RNEvSmKg34FPvddHkZWC2InVwLCRNDdka4HJJt4UnNoODjeYGYEK8czkMnyIDYGZdeJrh+ZKGSuon6QJJE3sbcxy7FnhR0nBJAyVNCPGrwIOSxssZIukmSb3lEAN/y2F26BzDwelW1gN7JD0Wi1f9JY2r4vWyn/HnrhXLzTNGfIl7pivMrC8mJjzpSYOaVMPv+PPBLyT9gRvSdmAWgHlCvGeBt+TZU9vxxSvM7FdgKr6YsxNPU/JZSbGZfYAvHG3EF3RWl/V9N77AswkPSvw2/ny0Eu7CF6224ItQM6PPr/CFrJbQuRV/HlsJT+LT/A7c2P87xTZ/X/QWPCJ9B55FdhFQac6kecCceFzwSA/yhcDt8YbB84X9S4EryOl+3cgA00ndkNSGL0QtqvdYTgbC824FxprZgXqP51QkPdQkOQmQNBBPy7wojWn9SIOaJCc4ki4FuvFHIAt6aZ4cQ3LKnyRJUiPSQ02SJKkRaVCTJElqRBrUJEmSGpEGNUmSpEakQU2SJKkRaVCTJElqxD/c9QSoKw9ycwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median sequence identity: 0.23636363636363636\n",
      "Mean sequence identity: 0.2343447385409942\n",
      "Min sequence identity: 0.18181818181818182\n",
      "Max sequence identity: 0.6727272727272727\n",
      "Second largest sequence identity: 0.36363636363636365\n"
     ]
    }
   ],
   "source": [
    "plot_max_sequence_identities(max_sequence_identities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:skipguide_data_processing] *",
   "language": "python",
   "name": "conda-env-skipguide_data_processing-py"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}