03-MHC-HLA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install the PGR-TK package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If neccessary, install the package into your python environment \n",
    "```\n",
    "pip install --force /code/pgrtk-0.3.2-cp38-cp38-linux_x86_64.whl\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pgrtk\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import SVG\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'pgrtk 0.3.2 (main:00d8eb3, release build, linux [x86_64] [rustc 1.62.0 (a8314ef7d 2022-06-27)])'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pgrtk.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"/results\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lazy load the the AGCfile for fetching sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_db =pgrtk.AGCFile(\"/data/pgr-tk-HGRP-y1-evaluation-set-v0.agc\") # lazy load an agc file of the reference without any SHIMMER index "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the index into memory\n",
    "The loads the index into the memory. It takes about one minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.7 s, sys: 2min 20s, total: 2min 40s\n",
      "Wall time: 2min 2s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sdb = pgrtk.SeqIndexDB()\n",
    "sdb.load_from_agc_index(\"/data/pgr-tk-HGRP-y1-evaluation-set-v0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the specification about how SHIMMERs are computed using the `get_shmmr_spec()` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(80, 56, 4, 64, False)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sdb.get_shmmr_spec() # window_size, k_mer_size, hieraical_reduction_factor, min_span, use_k_mer_sketch_other_than_minimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Query the database with HLA-Class II region sequences "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a copy of the annotation for IGV hg19 to look up the chromosome and locations of the gene of interest.\n",
    "We can fetch the region of interest with `ref_db`. Let's fetch the whole region and do a dot plot to examine the repeat structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "49\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0, 828999.0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD4CAYAAAAtrdtxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABG4ElEQVR4nO29e5QU1bX4/9nyEFGHh4oIiCOR+MCoPESNMUaN4M1kRfNSzDUh6l0Y4zXh5pulgzFXfETH3PwSfyY3Br5RQjRBvSZBf+kocNUkJAEEQUHEAcRRGUA0wzDoODLg/v3RVVDdXd1dVV3VfbrnfNbqNT2nz2OfU6dq1zlnn31EVbFYLBaLJQoHVFoAi8VisVQvVolYLBaLJTJWiVgsFoslMlaJWCwWiyUyVolYLBaLJTK9Ky1A3Bx++OFaX19fcj57VeklEio+ECpNJQhbr1rFrx2qsW26935In172XTBuyt0XsssLWn6QeMXi7FXl5S0d7Nm5nb2dO8NXWlVr6jN+/HgtlZ3v79bvz1+jO9/fHTj+jY+9qDc+9mLgNJUgbL1qFb92qMa22bzjPZ304z/r5h3vVVqUmqLcfSG7vKDlB4lXLM7Hvv9HPebG9Kfv0ONUIzxzRWtsn8iECRN0xYoVJefT0dVNXb8+oeLv6upm+MD+JZedJGHrVav4tUM1tk1re6fxfa4aKXdfyC4vaPlB4vnFWfDSFq55aFVG2JY539qz+61XQ1fajoPzEKUDzfrLJlrbO4H0hTORantIJoVfO1Rj21gFkgzl7gvZ5QUtP0i87Dj1jakcBdLS1ED39k0vBio0C6tEYqKuXx+uOXfUPkXyowXNxioSi8XS87jt8TXUN6YywkYO6EtLU0NJ+QZSIiLyHyKyVkReEpF5ItJPRAaLyCIR2eD8HeSJP0NENopIs4hM9oSPF5E1zm/3iqRXe0TkQBF5xAlfJiL1njRTnTI2iMjUkmqbMMMH9ue7k4/f97ca32wtFkvtUd+Y4oElb2SEtTQ18NcZF5acd1ElIiLDgW8BE1T1ZKAXMAVoBJ5W1dHA087/iMhJzu9jgIuAn4tILye7+4BpwGjnc5ETfjWwQ1WPA34C3O3kNRi4BTgDmAjc4lVWQSj3aMBVHKYpEDsqysSvPWwbWVzK3Reyy4ur/BseWZUz+jhl2CH+ow85oFduYHGCTmf1Bg4Skd5Af2ALcDEw1/l9LnCJ8/1i4GFV/UBVXwM2AhNF5CigTlWXaHo1/9dZady8HgMucEYpk4FFqtqmqjuARexXPEXp6Oq200rYdsjGrz1sG5VOrbRduftCdnlxlV/fmOLRVVsywlqaGnjiW+f6ytCr7vBhUcopqkRUtRX4EfAGsBXYqaoLgSNVdasTZyswxEkyHHjTk8VmJ2y48z07PCONqu4BdgKHFcgrAxGZJiIrRGTF22+/vS+8rl+fkqaVku5E5eqkpbZDMart4eHXHnG0UTnaIVvxmUItKeGk75di5fmVH6ZdP3fvX3JGH+eMGlhw7aOuXx/2dryzJW+EAgSZzhpEeqRwLDAMOFhEriiUxCdMC4RHTbM/QHW2qk5Q1QlHHHFExm+lKJAkb4py33RJKpBqfHjEbZ1VjnbwlmFau5f7wZs0JllnhbnW9Y0pVm95NyOspamBB6edXVwI/XBvMGkzCTKd9WngNVV9W1W7gd8DHwfecqaocP5ud+JvBo72pB9Bevprs/M9OzwjjTNlNgBoK5BX4gS9KaKa9NbKTVcr9SiVcrSDt4zs8tx+WEl6eh/IvgZxKfggfev8Hz6dM/qYdMLhJVteBSGIEnkDOFNE+jvrFBcA64AnANdaairwuPP9CWCKY3F1LOkF9OecKa9dInKmk8/XstK4eX0JeMZZN1kATBKRQc6IaJITVhaCKJCr5iyneVtHpLfCWrnpaqUepVKOdsiegoP9/dAERdJTyb4GcY8UC/Wt+sYUm9q6MsJamhqY/fUzYim7GIF2rIvIrcBlwB5gFfBvwCHAo8BI0ormy6ra5sT/HnCVE3+6qj7phE8AfgUcBDwJXK+qKiL9gAeBsaRHIFNUdZOT5irgJkeUH6jqnEKyxrVjPSjujuFq3O1sqR3szvXyke9ez74GST8TzrzjKba9mzkDNf28UUyffGLovDq6uhlwUN/nVXVC2LTW7YnFYrEExB1hVHIKt6Orm1NmLswIO0hg3V3Rpq7cOt3++VNf0A/3jg2b3ioRi8ViCUElZx0m3v4k29/7MCPs743nlTwKLWUkYt2eWHokplg1+WGabKbJU2kqoUDmLH6V+sZUhgIZfFB67SOOacxS6mSVSAFqZZ9I0lSqHlF3+SZlHhtHfmFkK9e+FJNMiXsi9Y0pbk29khHW0tTAylsKT1+V65pZJZIHE/aJVNr1QtA0lXjIlLLLN45NqH4KLI52CCpbudq91I1vlujc/LsXc8x2hx7SK5DZbjnvS7smUoCk5z4L5V/uBbxSyqvUHHHUMxjyxQ16NsOdqXUA3NRwYuTy46AS7W7CwnI5qLS1ZbbyAELv+XDrELQuckCvSAvrdiRSgLCHUsWZf6VdL4RNawJhFEhU31l1/fpwU8OJOQokTPlB5AuCW15cb5tJj+KSHNXHnV/Ut/hSZZn2q2WxuWt3FUiQuiTqO8tSnKSGjpV2vWAypbR5qb6z3B3jSRC2XnH1vbDTgUnmX+l8oyrKUmWpb0yx8JV3MsJKddcetC6l+M6q+emsqMPSsOkqPfztidRqm1eq71Vy+tbEfKMQRZaGe/7M2m3vZYSdN3owc64+q2wyWRPfPER9M4iSzpRO3JOo1TYPW6+42iHp9kwqf5P6QVhZ6htTOQqkpakhVgUSxIDnRwuaI58nYkcieTDp7cZisdQWn7xrEW/s3J0R1jBmCP/91dNjLyvpkUjvyJJVCVEVgVUgFoslCbIXzgV4LUFvu0HX+aJS80rEYrFYTOCkm1J0Znos4a5LxnD5mfUVkScurBKpcap9Wi5J+au9bSzVQWt7J2c3PZsRdiDQXIazPspBTS+s93Sq3WVFkvJXe9tYqoMx30vlKJDfX3tWzSgQ6AEL6z2dan/btiORYNRSXWqBOYtfzfF3dUhveOmO5JRHqX1AROzCehIUc01i+o1runzFSPq42VrA64oEaqde1UocLkvCUkl3NHY6y0MYp3rub83bOnzTJiFPT8DPuWGlqIRzzCj139XVvU+BlDJFl3S7Z+dfynG+Jh4F7OcwcWC/aD6vwlLXrw//esbIirxAFFUiInK8iLzg+XSIyHQRGSwii0Rkg/N3kCfNDBHZKCLNIjLZEz5eRNY4v93rnLWOcx77I074MhGp96SZ6pSxQUSmkhB+CqOQywD3on1r3qrIZ6yHkSfuvE3EdW54Z2pdxWUM8gKRhAuPsPV3z/be5YyKS/FrVU6v1aWcC2/imfL1jSkeWr45I6ylqYEXZoZXIFGuQ2t7J9+at6oibRJqTUREegGtwBnAdUCbqjaJSCMwSFVvFJGTgHnARGAY8L/AR1V1r4g8B3wbWAr8CbhXVZ8UkW8Cp6jqN0RkCvB5Vb1MRAYDK4AJgALPA+NVdUc+GUtZE4kyPZXkGeteL5xxDVVN98Lq3jwmyFaJqcwo9Y/rfPVyuz0pRW5TzpS/8v4lPLuhLSNs5IC+Jfm7inodSm2TqGsiYZXIJOAWVT1bRJqBT6nqVhE5Cvizqh4vIjMAVPUuJ80CYCbQAjyrqic44Zc76a9x46jqEhHpDWwDjgCmuHGcNLOccublk7FWF9bjvMGrYS3HYjGdSqx9JEm5FtankB5lABypqlsBHEUyxAkfTnqk4bLZCet2vmeHu2nedPLaIyI7gcO84T5p9iEi04BpACNHjgxZpeogzoe+VSAWS3Q+1fS/tLR/kBEWp8PEaiOwEhGRvsDngBnFovqEaYHwqGn2B6jOBmZDeiRSRD6LxWKJhGmjDxNmFcJYZ/0LsFJV33L+f8uZxsL5u90J3wwc7Uk3AtjihI/wCc9I40xnDQDaCuRlsVgsZeOiHz+bo0AuOWVoxRWICRtmwyiRy9k/lQXwBOBaS00FHveET3Esro4FRgPPOVNfu0TkTMcq62tZady8vgQ8o+nFmgXAJBEZ5Fh/TXLCLBaLJTYKPYjrG1O8sj3T6qmlqYF7vjI+abEKUu7TT/MRaDpLRPoDFwLXeIKbgEdF5GrgDeDLAKq6VkQeBV4G9gDXqepeJ821wK+Ag4AnnQ/A/cCDIrKR9AhkipNXm4jcDix34t2mqpmmEBaLxVIC+SwWq8FhYqUVCFi3J0UxYc6xXPT0uvak+oehJ7SLt47zV77J9EdXZ/x+eH9hxX9+piLylCsvOaDXC/rh3rFh87c71gtgypxjOejpde1J9Q9DT2kX9yFb35jKUSCrZ04quwKJq82D5tXR1U2vusOHRSnDjkSK0BPewlx6el17Uv3D0BPaZd7SFmbMX5sRVtcXVt9WmYXzahqJWAeMRaj1m8dLT69rT6p/GGq9XYqZ7VZCiVZkX5h+uLd4pFzsdJbFYumR+DlMrB94YI4C6QnTeaVgRyIWi6XHEXTToClmtCZjlYjFYukxnP/Dp9nU1pURVj/wQP7c+Om8aawCKYxVIhaLpUdgmsuSWsEqEYvFUtN8/AcL2LJrT0bYOaMG8uC0syskUW1hlYil7PQEk9FqpdaujR19JI+1zrKUFWvtYi61dG3OvOOpHAVyxekjrAJJALvZsIYx9a3SVLl6OiadKlkK2cpDgNes8ihK1EOp7EikRjH5rbLaH1K1iNtfqpkTGlM5CuSeS0+xCiRh7EikhrFv/JYwVGt/Wbx+O199YHlG2MG9YO0PrPIIQ7mOx7VUEdX4QLBUB6YoHL+F89UzJxkhW1RMadug2Oksi8USyttra3tnxadK5yx+NUeBDDn4AFqaGqrqAZyNydPQ+ajJ6axn/rYkUkeqtjeAStLT2iq7vknU3/vgSKptC8ldrE4dXd3cmVoHwPUXHMfwgf0Tk6UQJpvttrZ3VqxdSsUurDvsVY2kyb1vANX0FlAJqvFtqRSy65tE/d0H9MzH13Jnal0ibVtM7mIPrrp+fbip4URuajgxlgdl2Da8O7U2R4EMPaSXUQrkqjnLaW3vLB65ANX2chZoJCIiA4FfAicDClwFNAOPAPVAC3Cpqu5w4s8Argb2At9S1QVO+Hj2H4/7J+DbqqoiciDwa2A88E/gMlVtcdJMBW52RLlDVecWkrXUkQjge1SmJRM7Eqm9kUi5CSOLyaMPL3GMRCpF1JFIUCUyF1isqr8Ukb5Af+AmoE1Vm0SkERikqjeKyEnAPGAiMAz4X+CjqrpXRJ4Dvg0sJa1E7lXVJ0Xkm8ApqvoNEZkCfF5VLxORwcAKYAJp5fU8MN5VVn7EYZ1l0o1msVQCU+6By3/xN5a07MwIGzviUP7w75+skES1S2LTWSJSB3wSuB9AVXerajtwMeCOCuYClzjfLwYeVtUPVPU1YCMwUUSOAupUdYmmNdevs9K4eT0GXCAiAkwGFqlqm6M4FgEXha1kWCp185QyhdFTppYsyWPKdGV9YypHgbQ0NVgFYhhB1kRGAW8Dc0RklYj8UkQOBo5U1a0Azt8hTvzhwJue9JudsOHO9+zwjDSqugfYCRxWIK8MRGSaiKwQkRVvv/12gCqZRyk3rik3vaV6KDRvX+kzND5516Kc6avzRg82cvrKEkyJ9AbGAfep6ljgPaCxQHzxCdMC4VHT7A9Qna2qE1R1whFHHFFANHMp5cat9E1vqS6at3UUXQCuVF+qb0zxxs7dGWEtTQ3MufqsishjKU4QJbIZ2Kyqy5z/HyOtVN5ypqhw/m73xD/ak34EsMUJH+ETnpFGRHoDA4C2AnnVJKXcuFaBWILQ0dXNb5a9wb2XjzVqAXjCbX/KGX1cOnaYHX1UAUWViKpuA94UkeOdoAuAl4EngKlO2FTgcef7E8AUETlQRI4FRgPPOVNeu0TkTGe942tZady8vgQ846ybLAAmicggERkETHLCLBZLBNxR6/FD6yotyj7qG1O807l/guFA0qOPH142tnJCWQIT1O3J9cBvHMusTcCVpBXQoyJyNfAG8GUAVV0rIo+SVjR7gOtUda+Tz7XsN/F90vlAetH+QRHZSHoEMsXJq01Ebgdcxzi3qWpbxLpaLD2Gu1NrubFhjO9vpoxaP9KYYm9W2KwrxjL55GEVkccSDbtj3YMpZo0upslTDsrljjxI2xaKE3Q/QDmuYbYsd6fWct/iFq49pz6vIqkkC17awjUPrcoIG9gPXpgZfOrK1HujJ+4TsTvWHUyzcDJNnnLg7tpOase2t5xibVsoTtCdyeW4hn6y3NgwxlgFUt+YylEgLU0NoRWIifdGXDvWqw07EvGQ9NtN2LdXU9+2ksBbZyjfSKTQNamWkUjzto6CaxwmyOqOjrwMOfgAnvv+v0TKz9R7o9i1CIL1nWUAUS9A0gok6Nur+yZu4k2SBN43y7p+fcpSb1eBFLomheQIOmVRDmX4m2Vv5H0rN2HUVN+YylEgLU0NkRUImLOu46XYtQiah4mjrELU5EjE1EOpgrwRukrkpoYTjbxRkqJSSrOa57BdirVdpUYi0361jIWvvJMRdsKQ/jz1nfNiK8M04mjDahuJWCViID1pFGKpTarFYaJlP/ZkQ8MoRRFYBZJLLYwYqpmg/fnCHz3Dhnfezwg7fWQd//PNc5ISzVJhanJNpNJU47ymyfRUqxdTCNqf6xtTOQqkpanBKpAKUo5nkJ3OSgg7JRUvdiRSWQr1Z7/RxyWnDOWer4wvh2iWPLjKP6hfPTudZRimK5B8D4W4ld/K19sYd8zgkvOxCqSy5OsTdu3DXMrlmNVOZ1UxUYeq+aYn4p6GW/l6G1NmL2Xl69ZTjZdamOY8bWYqR4FccfoIq0AMoxwvs1aJGEzUHdXFKPSGkh1eygNv3DGDeXjambGMREyio6s7NgVebes8re2d1DemaO/aH9aX9Ojjji+eWhF5apVqedmwSsRQiimJUoeq2enc8sLIEIRaVCCluGbxXrdqMxg4rjHF2U3PZoTNumIs6ys0+ghyLkq1Uk3GOXZhPSGqcdORX3nWQCCXOF2zVIPBQBiHieXqL+5D9l/PGGmUW/s4Kfe9ZxfWDSKsVUQ+THh4myBDEpiyj8dPgZikuP3ctedb94ir3wch6UVjE65BtbzQ2emsBKjG42qrafhcKibX1RTZbv7di9RnKZCRA/oWXDhPst/7tUeSCsSEa+DFRJlc7HSWZR+mvukkgcl1rbRsppntlnOE4y3TtP6RtEzWi28VYeLbBKTf7Jq3dSReTpSzXuLGtAeEl0KyJdl3PnfvX3IUyJihB1fcbLeUEU7U9jKxf5goEwRUIiLSIiJrROQFEVnhhA0WkUUissH5O8gTf4aIbBSRZhGZ7Akf7+SzUUTudc5axzmP/REnfJmI1HvSTHXK2CAiU6lyyjksDVtG87YOLpu1JFFFErb+Jg/jy02h/T3F0hWjvjHF6i3vZoS1NDWQmv6p0HKWIkc+wipW1wy73H0njrKqra+HGYmcp6qneYY7jcDTqjoaeNr5HxE5ifQZ6WOAi4Cfi0gvJ819wDRgtPO5yAm/GtihqscBPwHudvIaDNwCnAFMBG7xKqtqpFzrJVFuoOOH1vHINWclau0Stv7VuL6UFH5tUew6F/v9E3cuzBl9nDNqYOyjj6Qe6H75es3Vyz0FVmodq/GlKdCaiIi0ABNU9R1PWDPwKVXdKiJHAX9W1eNFZAaAqt7lxFsAzARagGdV9QQn/HIn/TVuHFVdIiK9gW3AEaSV0adU9RonzSynnHn5ZLVrIvsxcV7XEj/FrnO+38u99pFUfzTJNL0aTftdkjbxVWChiCgwS1VnA0eq6lYAR5EMceIOB5Z60m52wrqd79nhbpo3nbz2iMhO4DBvuE8aSxGsAukZFLvO2b9//AcL2LJrT0bYVWeN5D8v/ljsshWSw4+Orm52dXWH2jvjl2+l+n4c5VbbfRtUiZytqlscRbFIRF4pEFd8wrRAeNQ0+wsUmUZ6moyRI0cWEM1i6dlkjz4EeM0Qf1cdXd3MfHwtq1vbmXvVxLyKxI6wzSLQmoiqbnH+bgf+QHp94i1nGgvn73Yn+mbgaE/yEcAWJ3yET3hGGmc6awDQViCvbPlmq+oEVZ1wxBFHBKmSpQaoRXcXYQk6dz7x9idzFMiDV51ujAKB9Bv4zIvHFFUg1bZmUOsUVSIicrCIHOp+ByYBLwFPAK611FTgcef7E8AUx+LqWNIL6M85U1+7RORMxyrra1lp3Ly+BDyj6cWaBcAkERnkLKhPcsKqkuZtHbS2d9LR1V3wAViph2MxuaLmGTVe9mKpF9fv1H3PrA+Uf756LXgp552kagjyQF35ehv1jSm2v/fhvrBD+6TXPs756JCMuEk/mIP0rbp+fQpOZVlDC/MIMhI5EvibiLwIPAekVPUpoAm4UEQ2ABc6/6Oqa4FHgZeBp4DrVNXd+Hot8EtgI/Aq8KQTfj9wmIhsBL6DY+mlqm3A7cBy53ObE1Z1NG/r4Is//zuXz17KjY+9yNT7n/O9qSrllM+dSpj6gL9cUfMM8tZYyMImn6nm8IH9ueTUo7h74QbmLW0pmH++NnV9QlWrIin2QB3VmOIL9y3JCPt743msud3f51WSb/hx9murQMzC7lgvE63tnXxl1lJ+MuU0jjvy0IKLh5VyyhdlUTNInkEXVAtZ2OTLZ97SFi4/s75o/vnadMFLW5h88rCi6auJ+SvfZPqjqzPChh3am398b3KeFGmSXmvwuwZ3p9ZyY8OYxMq0BCeqdZZVIkUodGOFvemqwWOrpbo5rjHFnqywSu84z8fdqbXct7iFa8+prylFUq0L/9btSQIUGuJHGf5bBWJJivueWU99lgIZffhBxioQgBsbxtSkAulpC/92JFKEOEciQaZnSsk/DuIus5Q2SiJ+qelMLCeuTYNJy5qdfy2OzFvbOzm0Xx/q+vUpe98stT3tSCQBil3MsA+77IXiQouMQd5oyuFCImi6OPLLjh9lQT5KOaVQyIIsqq+roHz+Z3/NUSBDD+lV8LyPfCT9Bp2df7Wd6hgEt067SvDbFbXP+LVn6GvZq0+0NwhVranP+PHjNQ52vr9bvz9/je58f3cs+bl5umze8V7R/Iv9Frd8rlxhKCZHWPnc+EHrF7X+cbSbV8Z88vr9H8d1O+bGP+Z8gspaKE6SZOcftq9VA946xdU3g/aZ7LLD9LPNO97TPocf06kRnrl2OstD9sij3MP7cqf3yy/KuQ3l9IlkGmGnKEuNd/4Pn2ZTW1dG2NgRh/KHf/9kbOVazCPMtXPjhr3e0rvvat2z+9SwstnpLAe/YWQ5jvisZHq//K45d1TofJNqp2p44HllDOOZuBh+/bG+MZWjQFqaGgIpkDDlWswjjAJx+03o+2dvd6SLb5WIg90Jm+6As/6yKfCahCU5vP1x3K2pnLWPK04fUTXu2i3loxLPMTudZcnAOxQu95GkSVKtUzm14q7dYj5Ju4LvcdSi+WEQ3AdILY3MqlEhfuz7KXZlDQjuumQMDaclexJCtbSPxRzsdJYPUcwPa3EKwPQHStA2z3cioInMX/km9Y2ZCuTw/kJLUwMNpw2vyHSTiW1lokwQn1zFHJKahFUiPgwf2J8Hrjw98EikJ84l3/DIqoqWH7bNwxwpm0259jLUN6ZyfF6tnjmJFf/5GaAyo0MT+3ZHVzd3ptYZJRPEJ1cxh6SmYZVIHsKerFZNUyWlcsMjq3h01ZaKKpJS2jxMWndU2rytI4qYgZi3tCVn7aOub3rtI1vGcvcxE/v2LgMfpHHi1+YmXgcXq0RiwsSLmxQXjDky42+l8O7PWLx+e5HY/mmL8W7XHqZMGMFvlr2RyFtgfWOKGfPXZoStnjmJ1beZ4/PKpL7tWhBef8FxRskF6Xa6qeHExI7INa2+LlaJWAri9+CcfPIwZl0x1ggX6h1d3fzHw6uY+sDyfYqk1Okn1zXNytfb+PxPF3Nr6hXGHT0g8MglOy8/bv7di4FHH6ZRSVcl7hu5qUYvpl+7JLAmvlXM4vXbc06ni5NKWzWF2dn94hs7OOejQ/ZNP4VZ0/Ky8vU2Hl2xmd17PqRv7wO4dMII3t7VVVRhutNd35q3al/Z+dqvWhwm+uFtX9fRoKU2sOeJOPQUJbJ4/XamPrCcuVednrgiqZQCiarAophnt7Z38tLmdm54bA1zrjyd4448FAj2Ztm8rYNLf/EPzj1+CNeddxzHD63LqIebx+W/+BtLWnZmpB01uB/P3HBBKFndfCul4F1PtdVmNm0pjFUiDj1FiUDyI5FKU6oCa97WkfFALxRv2twVtHXuZvwxA7n3K+ND+7Zq3tbBUQMPypsuiU2Dld4YWOnyLfGSuCt4EeklIqtE5I/O/4NFZJGIbHD+DvLEnSEiG0WkWUQme8LHi8ga57d7RUSc8ANF5BEnfJmI1HvSTHXK2CAiU8NWsJapZQUC0eeXF7y0JX2m/X3/KGhV1dreSUdXN79Z9gazp07gsWs/nqFAVr7eljdttsnl8UPrfOX9+A8W5CiQS04ZSktTQ8kL9ZV+gFe6fIsZhFlY/zawzvN/I/C0qo4Gnnb+R0ROAqYAY4CLgJ+LSC8nzX3ANGC087nICb8a2KGqxwE/Ae528hoM3AKcAUwEbvEqK4vFS0dXN1fev4RrHlrFMy9v4/0P9vLeB9mHxaZp3tax7+yH704+nuOH1mUogpWvtzFl9tK8iiSIyWV9Y4otuzLLb2lq4J6vjE/sTJNqolrltmQSSImIyAigAfilJ/hiYK7zfS5wiSf8YVX9QFVfAzYCE0XkKKBOVZdoeg7t11lp3LweAy5wRimTgUWq2qaqO4BF7Fc8Fss+Orq6+cLP/sazG9IP/X59evE/157FuGMG+8b9zbI3uPfysQwf2N9XEYw7ZjAPTzvTN71LPgVy5h1P5R19eNPGsZ5g8ia0QlSr3JZcgvrOuge4ATjUE3akqm4FUNWtIuLOrQwHlnribXbCup3v2eFumjedvPaIyE7gMG+4T5p9iMg00iMcRo4cGbBKllpjwrGHceaxg/nIkEO48pyP5I1X7AG+8vU2jhzQr6ACyUe28hDgNZ+1j7geniZvQitEtcptyaWoEhGRzwLbVfV5EflUgDzFJ0wLhEdNsz9AdTYwG9IL6wFkDITph1KZggn1CLvRK1+8+Svf5Lv/s5phgw5i3rQzA1t5ndCYoisr7ME8lnMrX2/joaVv0Lf3ASVtTnPbvdJtH5VqlbsSRLnHynVfBpnOOhv4nIi0AA8D54vIQ8BbzhQVzl93y/Bm4GhP+hHAFid8hE94RhoR6Q0MANoK5JU45T5zuloxqR6lnhLZvK2DW55Yy5mjBvF/vzYhkAJxXZZ4FYi7adBPgTRv6+DKOc/xfveeknZdm9TulmSJcq3L2T9Cmfg6I5HvqupnReS/gH+qapOINAKDVfUGERkD/Jb0Qvgw0ovuo1V1r4gsB64HlgF/An6qqn8SkeuAj6nqN0RkCvAFVb3UWVh/HhjniLASGK+qec1m4jTxtSORYJhajzCbFd09D1vb3y9oquvlxBkp3s+6ffKZ7YYxBw6K+4Awse0t8VKOkUjiJr4+NAEXisgG4ELnf1R1LfAo8DLwFHCdqu510lxLenF+I/Aq8KQTfj9wmIhsBL6DY+nlKIvbgeXO57ZCCsQlqvbNTlfp43Er5fI7bLlxt1NcFkvF3sRc9x3e+XmvhVa+tHen1lLfmKlAhhx8QEEFEsQcOArV7Bq+FPcppo7A4pSr1KO6y/Zyoao19Tlt3Dj9/vw1uvP93RqGne/vjpQuKXa+v1tvfOzFssnjlrd5x3uJtkOxfOO8Dt48svPbvOM9nfTjP+vmHe+FkuOYG/+Y8wkrS5yUu7/G1S+LtX8xGUy6V13i7rvlvP9VVYEVGuGZW3MOGHuJRLL6sNYiaQ5NsB2CjA7ivA5uHq3tnTnlFjszJluOab9almN5dcKQ/oF3nSfVr6q1v4Y9s8eLqfeqqXIljXV7YjDlXmsIU15UlytJ1WnBS1t8nSS6iuuac0dF9vzq57Jk9cxJsdVj5ettjDtmcGA3LZWm0mtglS6/XJS7npVYE6k5krTEikIlDiAKguv8MewZHkHLCNte85a2cM1Dq1jwUq7hXimuwz9516IcBTJ2xKEl+7zyMm9pC5fNXsr8lW9y2awliR5+FRdx9cso94XpVmlxylUtitIqEYekOqfpnT4K53x0SGLeg6McXTt78WsMOaQPJ48YuC/MS5Sbsb4xxRs7d2eEtTQ18Id//2RBGf0Wi/O5Tlm8fjvff3wtg/r35vRRh/HINWdVxUgkDqLeFyZPGdXivR4EO53lIanhY08ZfsdF2PZyH9zDB/Yv+TyRz//sr6zavCsj7NKxw/jhZWOLyrjy9TamP/JCxiZF1weX14WK+5C5M7WOho8NZdSQQ4w9ZClJavG+qOY6RZ3OCur2pEdgFz/9KfeNka+sfHJ4H8BBF2z98grjrt1PgXx73iq27OzirZ1d+8rP9sG1eP125r+whe9M+igAp44cVPX9Iyom1ztqn49zqq/seckBvYpHysVOZ1kKYsoQvaOrmztT6wLJEUSBeOs07tZUjgK5peGEwGsfzds6mPH7Nfy/l4/lf76R6/Rx3DGD6ehKnwN/5ZzlLNv0T4DYzuO2xEul+3zcHp6D5NXR1U2vusMjnXdtp7PKiHviXrGdxqYNiU2Qx1Ui2Q/e+55Zz7XnfzRSfru6ujm76dmM8F7AqyEWzt2b9F/PGJl3PcOVffeeD9n5/m5uuOiEsq19tLZ38tbOrkjOJJOQpVqm7Srd5ysxEpEDer2gH+4dWzRiFnY6q0y4c/X3Xj6WOX9vAfzfRL0uOOLoRHF0RhNuJj8Hi/c9s567F24A8FUkhep+zl0L2flBZtisK8YWPUvdi/tQLHat6vr14cqz6zlq4EH7/i8Hre2dTPnFErbt+oBHiri1L4cspaxVVTuVVkpJYqezyoQ7V3/80Dpuajgx71RGnNYnlR6Wl0q2/Nltcu35H+XGSaPzKhC/ui94aQv1jakMBXJon/S+j7+/2hbKKuyqOctpbe8seK06urppbe/kW/NWsSuEL68w5Is/fGB/Hv7GWSUrkDj6j99alan9Mu77Jmx++eInaRJtp7M8mDydVQmq/Q2oFPmz0xZbOI9iFVbordo7qtzV1R3oDTzsSDTukWu58k9a7lKJ+74Jm192/FLaK+nprJpUIs/8bUmkDlCOB24554VrSYG434vVya9971mwjnue3ZQRdnh/YcV/fqZkueKMmy9NsTyClFFKv7Om72aQdHvZHesOe1UjDUXLMfXjnQJJmlqaynK/+/nA8uLXvvWNqRwF0tLUEEmBeGUJ2q6lmokGKS+IAiml31nTdzMwtb3sSMSDHYmYRSkjkem/fZ75q7dl/FY/8ED+3Php3/yDyOJOJwCBZHFp3tbBex/sibwuEcd1rCbLKEtliDoSqUklYtdEejZBNg0GmWMuNK0UdI56/so3+d7vV7P7Q3jkmsw9JPbBbjEJO53loVqncOKmp7XD5+79S44COat+AC1NDb7+tAopAL+pM2/cIFZ0N//uRf7j0dWMPWYQD3z99BwFUq6pTUv1Y/K9XHNKJOqaSK1R7WsiYalvTLF6y7sZYS1NDcz7xifyPrALjUBm/WUT15w7irp+fZi/8k3feIUUyD0L1vHQ8s0oMKzuoBxnlV6T12rw3GupHKbfy0WViIj0E5HnRORFEVkrIrc64YNFZJGIbHD+DvKkmSEiG0WkWUQme8LHi8ga57d7RUSc8ANF5BEnfJmI1HvSTHXK2CAiU4vJG/VQqnJQzoeFyd5O4+QTdy7MGX2cM2pgxvRVIX9afh52ve7j5yx+lemPrs6rSPzo6OqmrWsv088bxRWnj8hx3uiVa/H67Vz6i+RcwLvGAJbqxfR7OchI5APgfFU9FTgNuEhEziR9DvrTqjoaeNr5HxE5CZgCjAEuAn4uIq5jr/uAacBo53ORE341sENVjwN+Atzt5DUYuAU4A5gI3OJVVvmoRGMXu1Gbt3Xwpfv+QfO2Dpq3deR1Dx4kr6BEMS0N87AMQpLTNfWNKTZ3ZLZVS1MDD047OyduPgVy6awlvuei1PXrQ2t7J7957k0mjBzA+ScNDSyXe9NPn3wid3zx1LzxWts7ufX/e5nTjx20bzc7xLfZ0HW38s0HV+xTJlahVAfVdJ2KKhHn+F13nqCP81HgYmCuEz4XuMT5fjHwsKp+oKqvARuBiSJyFFCnqkuc83x/nZXGzesx4AJnlDIZWKSqbaq6A1jEfsVjDIGHm/ohLe+8yxf/ezGXzVriq0gqNXTt6Orm2l8vD/3WXYik5v39Dou69pz60IdFvffBHvZ8CL/6R4tvew8f2J+5V0/kgavOCK2Qg8R38//xZWNzFuxL3d3synDmsYP426ttPPPyNu5MrQvsxNJSObKvqenTWYGss5yRxPPAccB/q+qNItKuqgM9cXao6iAR+RmwVFUfcsLvB54EWoAmVf20E34OcKOqflZEXgIuUtXNzm+vkh59fB3op6p3OOHfB95X1R9lyTeN9AiHkSNHjn/99dejtkdkgphhusefFjP5rJRpbkdXN8+8vI1Lxh0dW55xWyBlK48wDhP92nXx+u3GuWMvdXdzNu7Rwfncx1jMI+yG0zhI9DwRVd0LnCYiA4E/iMjJhWTxy6JAeNQ0XvlmA7MhbeJbQLbECHKBXc+txTy4Vuomr+vXx1eBlNKB41IgJ9+c4t09mWEPhjhdMZ9Jrje9KSa3cY96XKeS5XT8aEI7VjPZ18pkxR/KOktV24E/k55SesuZosL5604sbwa8T6IRwBYnfIRPeEYaEekNDADaCuRVEFOHfWC2bH5Ueig9f+Wb1DdmKpABB6bXPsIez1vMpLeWTW7DTI+VQq23oyWXINZZRzgjEETkIODTwCvAE4BrLTUVeNz5/gQwxbG4Opb0AvpzqroV2CUiZzrrHV/LSuPm9SXgGWfdZAEwSUQGOQvqk5ywvJhs4usudMbpHTSOOIWopGVIfWOK6Y+uzghraWrgxVvDrX247e4X7hL0RMQ4KOf1d+MF6XdxvDCUsx2rCROfR3ERZCRyFPCsiKwGlpNe6P4j0ARcKCIbgAud/1HVtcCjwMvAU8B1znQYwLXAL0kvtr9Keq0E4H7gMBHZCHwHx9JLVduA251ylwO3OWF5ScLE18QOEOSGj/pQyI4ftwIpJs+V9y/JWfuoH3hg6IXzYjJkt025FEgcLzlJjBDjemGwCiSTSo/mk8a6PfGQz1eT31x61EXKOBfIguQV1j8UENoVedhF4EL5hznnPAx+cpbiHbeU6xhXHwh7bU2bVzdRplLJV6dqqKt1e1Iifl5j3Qvvp0CimkuGfeCWmldYB4NQeO3AL02YNsj3tjvtV8tyFMjYEYfGNvrwq08xBZKvbqW+Wcb1MPG+8JSrzLiIe2rXBArVybT2jxM7EvHgNxIpFBeS6xxBHfzFXWYlRlZJjT5KJamRSJwE6Scmyeq9v+5Mrct7wmc1Uu11sl58HWrJi68pN39SfP5nf2XV5l0ZYeeNHsycq8+qkETVSTFlZ8IJgn5y1GL/ruY6WSXiUEtKpJYxdfRRi5jyYDNFDos/iW427OmUezG80iQp48Tbn2T7ex9mhF06dpivk8JqaKtqIM7d7+WSIy5sH0oeu7BehDCLqOWwww9DlHKSlLG+MZWhQPqSHn3kUyAmmkWaJk+cmNrmUam1+piKnc4KQFBT2h8taOaac0cVtJMv15tRKXPhccv4iTsX5njbDeKyxLS3SFPWF5LEtDYvlVqrT5JYE98ECWpKe825o5j1l01533y8e0+SdgtRbONYPtPV5m0dsd10zds6cty1H9wruMuSIHKU072Gqec6BBkBByVK3Qoda1BpTLtWtYhVIjEyfGD/vA8Z9y22tb2TW59Yy9T7nyuLIvHDb5jf0dXN936/JrYDkk6ckWLyPYszwhZMP4e1P4hv8bwSfppMeygVm7JJekpn5ettTJm91GhFYkkWO53lIemhr3cksqurO3H3EGH3OXR0dbO1/f2iXoYLcd8z67l74YaMsLq+sPq2ZCyvrMfY0nbex8HK19vyHmtgqR6sdVaJlGO+2823rl+fxN9oi9UnX1jd0OhyVcJst6crECg+Okq6r1kF0rOxIxEPtbYIV676XPfgclJrM4+YHXpIL5bebNwhlBaLJQ92JBID+Rz0JfUwTvohXw4FYjcN7qfcLyG19tJjyU+lXBIFwS6s56G1vXPfQnghR3xRqXYb9useXJ6jQE4ZdkjJCsTE9giycB/H9QxjIFDt/ccSnCjXupz9w05n+ZC95yPfInSpayjV+iaZpLt20/ZhuBZgNzecmOi+ltb2Tib9+C8s/M65gdd5qrX/lEpPrHc5RiLWd5ZDXJsNg24w7Emd+cs/X8zyNzLNf88ZNZAHp50dWxkmtuni9du5ft4qHrnmrJIs1wpx5f1LeHZDm3VAWQQTXzRqBatEHKwDxmTo6Wsfzds6ElMgLlfev8QqkACY+KJRCyS2Y11EjhaRZ0VknYisFZFvO+GDRWSRiGxw/g7ypJkhIhtFpFlEJnvCx4vIGue3e52z1nHOY3/ECV8mIvWeNFOdMjaIyFSqkGqetz5tZipHgVx7Tn1FFEjS7Vgo/6QVCGAVSECsAjGLIAvre4D/o6onAmcC14nISaTPQX9aVUcDTzv/4/w2BRgDXAT8XER6OXndB0wDRjsf1wb0amCHqh4H/AS428lrMHALcAYwEbjFq6wqTVCnjCYugAaRp74xRXtXZlhLUwM3NoxJSKr8JN2OlbpOQcszqf+YJIul8hRVIqq6VVVXOt93AeuA4cDFwFwn2lzgEuf7xcDDqvqBqr4GbAQmishRQJ2qLtH0HNqvs9K4eT0GXOCMUiYDi1S1TVV3AIvYr3gqStCHjon+lorJ7jf6mH7eqIpOXyXdjpW4TkH7kEkvIibJYjGDUCa+zjTTWGAZcKSqboW0ogFc05XhwJueZJudsOHO9+zwjDSqugfYCRxWIK9suaaJyAoRWfH222+HqVIg/G6YMA8dN47rk2rl620FzTmTvkFdZ5HZsi94aUvO6OPQPunRx/TJJ0Yqy8+nUtS3b7+2nre0JZJcfoQ9dyOO8oL0IW+8bH9nQWQJ4wutWFw/ma1C6dkEViIicgjwO2C6qhbqaeITpgXCo6bZH6A6W1UnqOqEI444ooBo4cn35pW9uBfkRr5s1hLmr3yTy2YvZcovlvgqEm95hfIsdY9KtrfhYxtTXPPQqox4q2dOYs3t0Ucffs75gjgMDBIP0gpkxvy1sSqSILS2d3Jnal1siiRovOy+UWwfE+zvd0EUSb64xfphLY5MTKiPCTIEQlWLfoA+wALgO56wZuAo5/tRQLPzfQYwwxNvAXCWE+cVT/jlwCxvHOd7b+Ad0gpkXxznt1nA5YVkHT9+vMbNzvd35/z//flr9oVn/5+PV7buVFXV51v+qZt3vFewvEJ5Bi2vEG7anz/drMfc+MeMz1l3PBU532yeb/ln3rL9wrPbtRi/XfJaSfKFZef7u/XGx17U/3h4VUntX0r52d+D9rsgZMcN0tcr0Q5JEsf9VY0yACs0gD7I/hQ18XXWJuYCbao63RP+X8A/VbVJRBqBwap6g4iMAX5LeiF8GOlF99GquldElgPXk54O+xPwU1X9k4hcB3xMVb8hIlOAL6jqpc7C+vPAOKfYlcB4Vc3rd7pcJr5+I5G459PDeuENi4lmu9Vgvum+IZouZ1yUo6+bhgl1LLcMSR5KdTbwVeB8EXnB+XwGaAIuFJENwIXO/6jqWuBR4GXgKeA6Vd3r5HUt8EvSi+2vAk864fcDh4nIRuA7OJZejrK4HVjufG4rpECSxG9+fvH67Rn/x02hPEsp77bH1+QokBOG9KelqcF3zr0UwuZR6Rs3COXwwmwS2XX1895Qa5hwfU2QIQh2s2EA/HbJLl6/nakPLGeu55hXE95eilFo9OGtJxCLWxe7u7i2sdd4P9Vw/+ejo6ubAQf1tTvWIbnpLL8Osnj99gwFYvLN9NXZf2fxpvaMsDFDDyY1/VMZYd56xnFTVPONZQmGvcbm3/+FcGW//fOnvqAf7h0bNr1VIjFi6s1k4tqHxVJrmHr/B/UDGHUkYl3Be8heDwh7drdpHehz9/4lR4GcN3pw4gokzjnyJOfba3Eu31I5ih1RXAnCbIqOilUiDtk2+DMfX8vUB54LfcaDKdQ3pli95d2MsJamhsT9MwXdBxJHXmHy99vnU4v7G0qhltrCpLpUsq+VwxODnc7ykL0esKurO9TZDibMiY67NUXb+5lhtzScwJXnfKRsMuQbPkdpozBTBPnyLxRu2uixUpjSf+PAxLpUQ1+zruAdeuKayLRfLWP218/ImbrqDWw0bO0j6TYqpMBMv4krTS21US3VpVzYM9YNoFIKZOEr7+QokFlXjGXyycPKLk8xkm6jfPnbB0pxTGyjqMrAxLrUKnZNpIppbe9k4Svv5IRfOnaYkQrEYgmDXbeqDqwSiYlyd/TRjSnObno2I+wLY4dzcG/4trNZMCym3aymyRMXXk8HpVDp9imHt2mT1jUs/lglEgNJvTH55Xfz716kvjGF95eB/dKWVzMvHsPC754X2BgguyyT3vpMkycuXE8HpSqSSrdPucq3CsR87MJ6TMS9kOdnYZL0pkHTFiNNkycuvJ4OSqHS7VPp8i3xYq2zHCppnQXx3lhuXpf/4m8sadmZ8dvoww9i0XfPj6Uci8VSOtWuVJP04msJSNxD/Lp+fahvTOUokJamBqtALBaDqPT0YiWxSiRG4lwI/Orsv+dMX50zamCOu/Y4MeEGMEGGYrheDUyh3LLEVZ5JbVgq7r0fJ9XSPlaJxEwhBeI9drRQB6lvTOV43G1pauDBaWcnuogf5FjeJInTzUlS6Tu6urkztS6243FLJUibhfUBV2p5+dJl+6arxTf3uOpUTe1j10TKRPO2Di79xRJ+9pWxnDpykK9bhk/cuZDNHZmd5orTR3DHF0/NCEtq7tXtsJV0GRGHm5MwZUVJ77aTKfPfhdqstb2Tq+Ys54ErT49ktRe2vHzx70ytA+CmhhNjPWrANJJYEy0XdmHdwVQl0tHVzXceWUXLO53MvXoih2adjmeSu/ZqurlLlbWa6hqV1vbO2BRIVKIq3p5wfUwhsYV1EXlARLaLyEuesMEiskhENjh/B3l+myEiG0WkWUQme8LHi8ga57d7nbPbEZEDReQRJ3yZiNR70kx1ytggIlPDVs4k6vr14ceXjWXu1RMZPrD/vhvjzDueylEg088bVdHzPqrppi1V1mqqa1QqrUAg2pHC1TSl05MJsibyK+CirLBG4GlVHQ087fyPiJwETAHGOGl+LiK9nDT3AdOA0c7HzfNqYIeqHgf8BLjbyWswcAtwBjARuMWrrMpB3J23rl+fjBu6vjHFtnf37vv/ANKjj+mTT4y13KSxN3maal3gNhW7Y30/Ua51ufpHUSWiqn8F2rKCLwbmOt/nApd4wh9W1Q9U9TVgIzBRRI4C6lR1iabnz36dlcbN6zHgAmeUMhlYpKptqroDWESuMkuMJN6C3AXOk29O5Yw+fn/tWWzynHVeCTq6ukPvpC7WTl5jgrjp6OpmwUtbSsojrkXncr81hzGECHMNTFNMVoFENzgpV3+Map11pKpuBXD+uttvhwNveuJtdsKGO9+zwzPSqOoeYCdwWIG8EsPb4HG/BbW2d/Kln/+D+sYU7+7ZH35on/ToY9wxg4H0DV8Jyx93zSasS45C7dS8rYPLZi3JeYgFrVu2NU/2b9f/5nmueWhVIEXipyzcRec4FEm535q9JqXFlPhls5YEuqZ2+shMovStcvbHuE18xSdMC4RHTZNZqMg0EVkhIivefvvtQIJm43cDxXkBzm16lq0dH2SEtTQ1sOb2/Wsfre2dXPfblWzbGZ9JZlDcNZsfXDImtEuOfO10/NA6Hvj66Rw/tG5fWNAHVfabtt+1ufOLp/CxYYdy1nFHFMwrn7IYPrA/d33hY7GtGZT7rdldZ7jm3FF5yz6kX29unHw8d6TWFVWWST944jQ17mmY7A4/qhJ5y5miwvnrvuZsBo72xBsBbHHCR/iEZ6QRkd7AANLTZ/nyykFVZ6vqBFWdcMQRhR8o+UjqBrrvmfXUN6bwDD4YcvABHH/kIb4PtVs+exKr3mhna3vW8YRlYGv7+9y9oDnyyCGbjq5u5r+wJdIIzxsvX5pN29/lpS27ePGNHQXzGj6wv6+Ja2t7JzN+v6bgw830t/KOrm5m/WWTr5yt7Z386/9dxl1PruPmhhMDKcskFUhcoz6LWQQy8XUspv6oqic7//8X8E9VbRKRRmCwqt4gImOA35JeCB9GetF9tKruFZHlwPXAMuBPwE9V9U8ich3wMVX9hohMAb6gqpc6C+vPA+McMVYC41U1e30mA5NMfAuZ7RYyu2ze1pHx9l5OssuOYy9Gkm9EpTozLHQdSq17uSi2T+Tdrj0V60/ZsphgKWbxJ7F9IiIyD/gUcDjwFmmLqfnAo8BI4A3gy+7DXUS+B1wF7AGmq+qTTvgE0pZeBwFPAterqopIP+BBYCzpEcgUVd3kpLkKuMkR5QeqOqdYheJSIqV0+OseXE5qbeYc9MgBffnrjAtLlqsS9GRb/Z5cdy9B28G2V3Xgd53sZkOHUpSI27Cl7PI1adOgxRIHQUdk1TJy6+m4HgS83gPAevEtGe/ibb459HzpOrq6ufBHz+QokPNGD45NgfhZJ4VNUyxdkvP/bt5xlmH6ekVUTKtXlHUsS2mUow/sKmD9GAarRByyb4CgCmTm42s5ZeZCNryTuRDe0tTAnKvPikW2bOukIBZOfnEKpUvSvNPNu7W9s0c6qAuDqfUKqhh6sgKJ65qV4uQyCHX9+nD9BcftM8hwy6NXn0gXz05n5cFdEyk0x+vnMPG80YNjUx5esuUIMvfsF6dQuijz2WHnyovFD+uAsRYfWqXUq1bbxHQ6urq59Ym13PK5MUDpyjSKk0u3/Cj3T2t7J8fUH7t9T/tbR4aV1Y5EfHDXRJq3deTdFVzfmMpRIAumn8OWjt2JmDFmd4wgHcUvjhsWdPNZsd+DvjG55RZTIGE2Wtbqw7IUBWLiKKYnsKurmxffTJvmx3ENwvYBt/xdIcr1ljF8YH/2drwTyf2DHYnkwTsSgf3u0T/ZtJD2rsy4088bxVXnHrdvUd50M8bF67cz9YHlzL3q9H3msdmLomEWU+N0fe234GcJjh2JVI4gsxflKD8q1jrLIal9IovXb+erDyzPCOsLrK9Syyu//RVRpszixj4ELZbKEFWJ9E5CmFrjhMYUWYMPHrzqdE4dWVanwrHiJ3v2YUFJrU0Uil9pJWaxWMJh10QK4Los8SqQoYf0oqWpgVNHDjLmiNSwxGmllVR8O78fDdte8VBKO/a0a2CVSB7qG1PcvXBDRlhLUwNLb97vjX73ng9900bxUpsUfov8hez5o9j6F3IAGDZ/t02qcc9B3NezmIGG394hq3hLp5R2TMo8N6osScTNxiqRLG57fE3OpsEhBx8QeNNgUAujctzwre2dTH3gubyKJB9hprHuTK3jp09vDFWPQgrE2ybVpkDivJ7FHBb6lVeNijcs5VCQpbRjlLTF+k6UvhUmzb59InJAr6KRfbAL6x78XJb8vfE8X4uHfJZEre2d/PTpjXktjLzz/EnP+be2d/L/LFjPzIuD246Hxe+Bb+I+h3Ksr8Rdhp+1TTn7j2kEtRisRuLcPxUlTUdXNwMO6mvdnkTl8z/7a44CGTmgb14FAumHZrai6OhKu+W+/oLj8iqQcr1pu7L8n8kfTbSc7AX4Ut/IXfPiOEliisHvt7jb2U+BVOtIzY+w18P7lm/qlF0p/b6U30tNU0pf6vEjkbgdJibxRpGULCaWm9TbZhQLsnxyVPKNuFZGH6W0oakjElPlCordJ+IQVIk03PNn1m57LyMsKZcllnCY8qAsJIcpMlYzJk57loqpcgXB7hMJgXXXbjam3IRxGB9Y8lNKG5ra/qbKlSQ9Sol8/AcL2LJrT0ZYw5gh/PdXT6+QRBaLxVLd9BglYkcfFovFEj81r0ROuilFZ9aewFsaTuDKcz5SGYEsFosliyTXUpJep6kKE18RuUhEmkVko4g0BknT2t5JfWOmAhnYLz36sArEYrGYQpIbj8uxqdl46ywR6QWsBy4ENgPLgctV9WW/+BMmTND2T9/KnqzwQns+LOYT9G0qX7y438aq2QrHFGrROisqJoxEavmM9YnARlXdpKq7gYeBi/NFXtO6M0OBuA4TrQKpXkp12hj325j1UVU6lfBPZTJJbwhOkmoYiXwJuEhV/835/6vAGar6754404BpAAf0O3R874FHoh9+uLd7+6YXKiFzAA4H3qm0EAUwTz45oBf64V5PiL+MufEKh8cnTzbmtWEmlZevlDaM+3pGo/JtWJwwMh6jqkeELaAaFtbFJyxD86nqbGA2gIis+GBrR+ghWTkRkRVRho3lwnT5wHwZrXylY7qMpssH5ZGxGqazNgNHe/4fAUQ6C9hisVgs8VINSmQ5MFpEjhWRvsAU4IkKy2SxWCwWqmA6S1X3iMi/AwuAXsADqrq2QJLZ5ZGsJEyX0XT5wHwZrXylY7qMpssHZZDR+IV1i8VisZhLNUxnWSwWi8VQrBKxWCwWS3RUtWY+wEVAM7ARaEwg/weA7cBLnrDBwCJgg/N3kOe3GY4szcBkT/h4YI3z273sn1Y8EHjECV8G1HvSTHXK2ABMzSPf0cCzwDpgLfBtk2QE+gHPAS868t1qknyeeL2AVcAfDZWvxcn7BWCFaTICA4HHgFdI98WzDJPveKft3E8HMN0wGf+D9D3yEjCP9L1jjHwZskZ5mJr4IX3jvwqMAvqSflCdFHMZnwTGkalEfoijsIBG4G7n+0mODAcCxzqy9XJ+e865sQR4EvgXJ/ybwC+c71OAR3T/A2KT83eQ832Qj3xHAeOc74eSdhdzkikyOnkd4nzv43TeM02RzyPnd4Dfsl+JmCZfC3B4VpgxMgJzgX9zvvclrVSMkc/nubENOMYUGYHhwGvAQc7/jwJfN0W+nDaM8yFbyY/TUAs8/88AZiRQTj2ZSqQZOMr5fhTQ7Fc+aeuys5w4r3jCLwdmeeM433uT3mkq3jjOb7NI+w8rJuvjpH2OGScj0B9YCZxhknyk9yE9DZzPfiVijHzOby3kKhEjZATqSD8AxUT5fNpyEvB3k2QkrUTeJP0g7w380ZHTCPmyP7W0JuI2vMtmJyxpjlTVrQDO3yFF5BnufPeTc18aVd0D7AQOK5BXXkSkHhhL+m3fGBlFpJeIvEB6WnCRqholH3APcAPgPUDAJPkg7bFhoYg877j8MUnGUcDbwBwRWSUivxSRgw2SL5sppKeLMEVGVW0FfgS8AWwFdqrqQlPky6aWlEhR9yhlJp88heSMkia3YJFDgN8B01W1wyQZVXWvqp5G+o1/ooicbIp8IvJZYLuqPl9AporJ5+FsVR0H/AtwnYh80iAZe5Oe8r1PVccC75GeejFFvv0Fpzcvfw74nwLylV1GERlE2snsscAw4GARucIU+bKpJSVSKfcob4nIUQDO3+1F5NnsfPeTc18aEekNDADaCuSVg4j0Ia1AfqOqvzdRRgBVbQf+TNoYwhT5zgY+JyItpL1Fny8iDxkkHwCqusX5ux34A2lP16bIuBnY7IwwIb3APs4g+bz8C7BSVd9y/jdFxk8Dr6nq26raDfwe+LhB8mVSaK6rmj6k34A2kdbe7sL6mATKqSdzTeS/yFzs+qHzfQyZi12b2L/YtZz0grK72PUZJ/w6Mhe7HnW+DyY9zzzI+bwGDPaRTYBfA/dkhRshI3AEMND5fhCwGPisKfJlyfop9q+JGCMfcDBwqOf7P0grYpNkXAwc73yf6chmjHweOR8GrjTwPjmDtGVWfyffucD1psiX045xP2Qr+QE+Q9oi6VXgewnkP4/0HGU3aY19Nel5xKdJm8M97W1w4HuOLM04VhFO+ATSpnuvAj9jv9ldP9JD642krSpGedJc5YRv9Hb8LPk+QXrouZr95oufMUVG4BTSprOrnbz/0wk3Qr4sWT/FfiVijHyk1xxeZL+Z9PcMlPE0YIVzneeTfhgZI58Trz/wT2CAJ8wYGYFbSZtIvwQ8SFpBGCOf92PdnlgsFoslMrW0JmKxWCyWMmOViMVisVgiY5WIxWKxWCJjlYjFYrFYImOViMVisVgiY5WIxWKxWCJjlYjFYrFYIvP/A9q2Ei5s9QoeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ref_file_name, roi_chr, roi_b, roi_e = 'hg19_tagged.fa', \"chr6_hg19\", 32130918, 32959917\n",
    "padding = 0\n",
    "\n",
    "#get a segment of a reference\n",
    "roi_seq = ref_db.get_sub_seq(ref_file_name, roi_chr, roi_b-padding, roi_e+padding)\n",
    "x, y = pgrtk.get_shmmr_dots(roi_seq, roi_seq, 32, 32, 2, 2)\n",
    "plt.plot(x, y, \".\", markersize=0.5)\n",
    "print(x[0])\n",
    "plt.xlim(0,len(roi_seq))\n",
    "plt.ylim(0,len(roi_seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of hits is 114\n",
      "The total aligned regions is 147\n"
     ]
    }
   ],
   "source": [
    "#find hits in the pangenomic reference\n",
    "aln_range = pgrtk.query_sdb(sdb, roi_seq, merge_range_tol=200000)\n",
    "\n",
    "print(\"The number of hits is\", len(aln_range))\n",
    "\n",
    "count = 0\n",
    "for sid, rgns in aln_range.items():\n",
    "    count += len(rgns)\n",
    "\n",
    "print(\"The total aligned regions is\", count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can examine, the target region sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('HG02818#1#JAHEOS010000026.1', 'HG02818.paternal.f1_assembly_v2_genbank.fa', 27438479) target region 3503365 - 4325934 len: 822569\n",
      "('HG02886#1#JAHAOU010000006.1', 'HG02886.paternal.f1_assembly_v2_genbank.fa', 55496981) target region 22522980 - 23317712 len: 794732\n",
      "('chr6_GL000250v2_alt_hg38', 'hg38_tagged.fa', 4672374) target region 3497464 - 3812023 len: 314559\n",
      "('HG03453#2#JAGYVV010000230.1', 'HG03453.maternal.f1_assembly_v2_genbank.fa', 5546504) target region 3402532 - 4193513 len: 790981\n",
      "('HG02145#1#JAHKSG010000017.1', 'HG02145.paternal.f1_assembly_v2_genbank.fa', 11219473) target region 6953506 - 7772155 len: 818649\n",
      "('HG01175#2#JAHALZ010000032.1', 'HG01175.maternal.f1_assembly_v2_genbank.fa', 31257284) target region 23360059 - 24179654 len: 819595\n",
      "('HG00733#1#JAHEPQ010000070.1', 'HG00733.paternal.f1_assembly_v2_genbank.fa', 39879370) target region 32136286 - 32923993 len: 787707\n",
      "('HG03453#1#JAGYVW010000148.1', 'HG03453.paternal.f1_assembly_v2_genbank.fa', 5669647) target region 1359413 - 1711637 len: 352224\n",
      "('HG00438#1#JAHBCB010000040.1', 'HG00438.paternal.f1_assembly_v2_genbank.fa', 31169400) target region 23356561 - 24160522 len: 803961\n",
      "('HG03540#1#JAGYVY010000082.1', 'HG03540.paternal.f1_assembly_v2_genbank.fa', 23258697) target region 16945983 - 17757447 len: 811464\n",
      "('HG00735#1#JAHBCH010000013.1', 'HG00735.paternal.f1_assembly_v2_genbank.fa', 39978829) target region 32215406 - 33034378 len: 818972\n",
      "('HG02559#2#JAGYVJ010000064.1', 'HG02559.maternal.f1_assembly_v2_genbank.fa', 59272483) target region 32122253 - 32942017 len: 819764\n",
      "('HG01243#1#JAHEOY010000117.1', 'HG01243.paternal.f1_assembly_v2_genbank.fa', 5690930) target region 3438255 - 4238724 len: 800469\n",
      "('HG03492#1#JAHEPI010000049.1', 'HG03492.paternal.f1_assembly_v2_genbank.fa', 20121680) target region 15861252 - 16689804 len: 828552\n",
      "('HG01071#1#JAHBCF010000017.1', 'HG01071.paternal.f1_assembly_v2_genbank.fa', 18736761) target region 299013 - 597352 len: 298339\n",
      "('HG02622#1#JAHAOO010000042.1', 'HG02622.paternal.f1_assembly_v2_genbank.fa', 33929730) target region 27501039 - 27783009 len: 281970\n",
      "('HG01123#1#JAGYYZ010000057.1', 'HG01123.paternal.f1_assembly_v2_genbank.fa', 32172171) target region 25895025 - 26182990 len: 287965\n",
      "('HG00741#2#JAHALX010000077.1', 'HG00741.maternal.f1_assembly_v2_genbank.fa', 33608122) target region 25651940 - 26443114 len: 791174\n",
      "('chr6_qbl_hap6_hg19', 'hg19_tagged.fa', 4611984) target region 3393290 - 4184339 len: 791049\n",
      "('HG01978#2#JAGYVR010000046.1', 'HG01978.maternal.f1_assembly_v2_genbank.fa', 9459559) target region 1449927 - 1807865 len: 357938\n",
      "('chr6_cox_hap2_hg19', 'hg19_tagged.fa', 4795371) target region 3602951 - 4403309 len: 800358\n",
      "('HG00741#1#JAHALY010000025.1', 'HG00741.paternal.f1_assembly_v2_genbank.fa', 56286211) target region 23365046 - 24159700 len: 794654\n",
      "('HG01891#2#JAGYVN010000027.1', 'HG01891.maternal.f1_assembly_v2_genbank.fa', 33395144) target region 27120084 - 27919956 len: 799872\n",
      "('NA21309#2#JAHEPB010000021.1', 'NA21309.maternal.f1_assembly_v2_genbank.fa', 26636088) target region 3340145 - 4140416 len: 800271\n",
      "('HG02572#2#JAHAOV010000201.1', 'HG02572.maternal.f1_assembly_v2_genbank.fa', 10833350) target region 5332858 - 6134375 len: 801517\n",
      "('HG01978#1#JAGYVS010000035.1', 'HG01978.paternal.f1_assembly_v2_genbank.fa', 60616846) target region 32160188 - 32981440 len: 821252\n",
      "('HG01928#1#JAGYVQ010000020.1', 'HG01928.paternal.f1_assembly_v2_genbank.fa', 40778302) target region 27374381 - 27732067 len: 357686\n",
      "('chr6_hg19', 'hg19_tagged.fa', 171115067) target region 32131290 - 32959865 len: 828575\n",
      "('HG01109#1#JAHEPA010000084.1', 'HG01109.paternal.f1_assembly_v2_genbank.fa', 34924654) target region 27049277 - 27790253 len: 740976\n",
      "('HG03516#1#JAGYYT010000073.1', 'HG03516.paternal.f1_assembly_v2_genbank.fa', 43524149) target region 23126159 - 23920647 len: 794488\n",
      "('HG01258#2#JAGYYU010000011.1', 'HG01258.maternal.f1_assembly_v2_genbank.fa', 39568798) target region 26230380 - 27044764 len: 814384\n",
      "('HG02622#2#JAHAON010000041.1', 'HG02622.maternal.f1_assembly_v2_genbank.fa', 60041455) target region 32133713 - 32447811 len: 314098\n",
      "('HG03098#1#JAHEPM010000086.1', 'HG03098.paternal.f1_assembly_v2_genbank.fa', 26959022) target region 22525026 - 22888592 len: 363566\n",
      "('HG01358#1#JAGYZB010000008.1', 'HG01358.paternal.f1_assembly_v2_genbank.fa', 39988763) target region 6939649 - 7296770 len: 357121\n",
      "('HG00621#1#JAHBCD010000020.1', 'HG00621.paternal.f1_assembly_v2_genbank.fa', 56350629) target region 23360091 - 23711816 len: 351725\n",
      "('HG03516#2#JAGYYS010000003.1', 'HG03516.maternal.f1_assembly_v2_genbank.fa', 59716798) target region 32273514 - 33096539 len: 823025\n",
      "('HG01358#2#JAGYZA010000082.1', 'HG01358.maternal.f1_assembly_v2_genbank.fa', 13116737) target region 6957238 - 7752650 len: 795412\n",
      "('HG02559#1#JAGYVK010000047.1', 'HG02559.paternal.f1_assembly_v2_genbank.fa', 59967020) target region 32167236 - 32465007 len: 297771\n",
      "('HG005#2#JAHEPN010000064.1', 'HG005.maternal.f1_assembly_v2_genbank.fa', 33628824) target region 7040698 - 7339378 len: 298680\n",
      "('HG03492#2#JAHEPH010000100.1', 'HG03492.maternal.f1_assembly_v2_genbank.fa', 11165420) target region 3433872 - 4223369 len: 789497\n",
      "('HG005#1#JAHEPO010000054.1', 'HG005.paternal.f1_assembly_v2_genbank.fa', 31959783) target region 5320458 - 5618569 len: 298111\n",
      "('NA18906#1#JAHEOO010000017.1', 'NA18906.paternal.f1_assembly_v2_genbank.fa', 34516894) target region 28384491 - 29178977 len: 794486\n",
      "('HG02257#1#JAGYVI010000022.1', 'HG02257.paternal.f1_assembly_v2_genbank.fa', 35352640) target region 27484494 - 28280971 len: 796477\n",
      "('HG00733#2#JAHEPP010000060.1', 'HG00733.maternal.f1_assembly_v2_genbank.fa', 27483925) target region 26893186 - 27481359 len: 588173\n",
      "('HG02055#2#JAHEPJ010000058.1', 'HG02055.maternal.f1_assembly_v2_genbank.fa', 27275867) target region 23120171 - 23914366 len: 794195\n",
      "('HG01123#2#JAGYYY010000050.1', 'HG01123.maternal.f1_assembly_v2_genbank.fa', 59396933) target region 32147422 - 32970028 len: 822606\n",
      "('HG02080#2#JAHEOV010000002.1', 'HG02080.maternal.f1_assembly_v2_genbank.fa', 31364656) target region 23360926 - 23649661 len: 288735\n",
      "('HG01243#2#JAHEOX010000097.1', 'HG01243.maternal.f1_assembly_v2_genbank.fa', 30954732) target region 26625704 - 27447351 len: 821647\n",
      "('HG00733#2#JAHEPP010000166.1', 'HG00733.maternal.f1_assembly_v2_genbank.fa', 551833) target region 5000 - 204370 len: 199370\n",
      "('HG00438#2#JAHBCA010000042.1', 'HG00438.maternal.f1_assembly_v2_genbank.fa', 56346258) target region 23361552 - 23719816 len: 358264\n",
      "('HG03540#2#JAGYVX010000013.1', 'HG03540.maternal.f1_assembly_v2_genbank.fa', 29520975) target region 23368683 - 24125419 len: 756736\n",
      "('HG01175#1#JAHAMA010000069.1', 'HG01175.paternal.f1_assembly_v2_genbank.fa', 11224566) target region 6945713 - 7774312 len: 828599\n",
      "('chr6_KI270758v1_alt_hg38', 'hg38_tagged.fa', 76752) target region 681 - 58584 len: 57903\n",
      "('chr6_dbb_hap3_hg19', 'hg19_tagged.fa', 4610396) target region 3481454 - 3709418 len: 227964\n",
      "('chr6_mann_hap4_hg19', 'hg19_tagged.fa', 4683263) target region 3475673 - 3773416 len: 297743\n",
      "('HG01952#1#JAHAME010000044.1', 'HG01952.paternal.f1_assembly_v2_genbank.fa', 33720779) target region 27450161 - 28271020 len: 820859\n",
      "('chr6_ssto_hap7_hg19', 'hg19_tagged.fa', 4928567) target region 3480367 - 3778125 len: 297758\n",
      "('HG03098#2#JAHEPL010000070.1', 'HG03098.maternal.f1_assembly_v2_genbank.fa', 28806995) target region 22510938 - 22862038 len: 351100\n",
      "('NA21309#1#JAHEPC010000026.1', 'NA21309.paternal.f1_assembly_v2_genbank.fa', 27383785) target region 3438116 - 4255912 len: 817796\n",
      "('chr6_hg38', 'hg38_tagged.fa', 170805979) target region 32163513 - 32992088 len: 828575\n",
      "('HG02486#1#JAGYVM010000005.1', 'HG02486.paternal.f1_assembly_v2_genbank.fa', 34332109) target region 26401396 - 27227220 len: 825824\n",
      "('HG02723#1#JAHEOU010000100.1', 'HG02723.paternal.f1_assembly_v2_genbank.fa', 28982488) target region 5052082 - 5841973 len: 789891\n",
      "('HG01106#1#JAHAMC010000024.1', 'HG01106.paternal.f1_assembly_v2_genbank.fa', 39876757) target region 6939124 - 7736742 len: 797618\n",
      "('HG01928#2#JAGYVP010000017.1', 'HG01928.maternal.f1_assembly_v2_genbank.fa', 59868708) target region 32100107 - 32398425 len: 298318\n",
      "('HG03486#2#JAHEOP010000002.1', 'HG03486.maternal.f1_assembly_v2_genbank.fa', 27376079) target region 3432636 - 4232522 len: 799886\n",
      "('HG002#2#JAHKSD010000045.1', 'HG002.maternal.f1_assembly_v2_genbank.fa', 32356768) target region 26179067 - 26986295 len: 807228\n",
      "('HG02145#2#JAHKSF010000005.1', 'HG02145.maternal.f1_assembly_v2_genbank.fa', 20185610) target region 3502917 - 4318859 len: 815942\n",
      "('HG02886#2#JAHAOT010000003.1', 'HG02886.maternal.f1_assembly_v2_genbank.fa', 33659657) target region 25612529 - 25968761 len: 356232\n",
      "('HG00673#1#JAHBBZ010000030.1', 'HG00673.paternal.f1_assembly_v2_genbank.fa', 56158886) target region 32028390 - 32824277 len: 795887\n",
      "('HG01071#2#JAHBCE010000076.1', 'HG01071.maternal.f1_assembly_v2_genbank.fa', 14785192) target region 6895908 - 7685015 len: 789107\n",
      "('HG02148#1#JAHAMG010000076.1', 'HG02148.paternal.f1_assembly_v2_genbank.fa', 56047404) target region 23140769 - 23938206 len: 797437\n",
      "('chr6_GL000251v2_alt_hg38', 'hg38_tagged.fa', 4795265) target region 3602845 - 4403203 len: 800358\n",
      "('NA20129#1#JAHEPE010000077.1', 'NA20129.paternal.f1_assembly_v2_genbank.fa', 27592746) target region 3426335 - 4225737 len: 799402\n",
      "('HG01891#1#JAGYVO010000024.1', 'HG01891.paternal.f1_assembly_v2_genbank.fa', 32274981) target region 25967087 - 26785111 len: 818024\n",
      "('HG002#1#JAHKSE010000066.1', 'HG002.paternal.f1_assembly_v2_genbank.fa', 32583635) target region 5444315 - 5742327 len: 298012\n",
      "('HG00621#2#JAHBCC010000005.1', 'HG00621.maternal.f1_assembly_v2_genbank.fa', 59876178) target region 32115696 - 32903256 len: 787560\n",
      "('HG02818#2#JAHEOR010000019.1', 'HG02818.maternal.f1_assembly_v2_genbank.fa', 18440780) target region 16915483 - 17713098 len: 797615\n",
      "('HG03579#1#JAGYVU010000035.1', 'HG03579.paternal.f1_assembly_v2_genbank.fa', 21125918) target region 16940123 - 17734844 len: 794721\n",
      "('HG00735#2#JAHBCG010000038.1', 'HG00735.maternal.f1_assembly_v2_genbank.fa', 31959602) target region 3503667 - 3800803 len: 297136\n",
      "('NA18906#2#JAHEON010000020.1', 'NA18906.maternal.f1_assembly_v2_genbank.fa', 29624882) target region 23346958 - 24162283 len: 815325\n",
      "('HG01952#2#JAHAMD010000016.1', 'HG01952.maternal.f1_assembly_v2_genbank.fa', 59464697) target region 32160890 - 32459032 len: 298142\n",
      "('HG02630#2#JAHAOP010000058.1', 'HG02630.maternal.f1_assembly_v2_genbank.fa', 24764972) target region 23131055 - 23500654 len: 369599\n",
      "('HG01361#1#JAGYYX010000108.1', 'HG01361.paternal.f1_assembly_v2_genbank.fa', 13154969) target region 6944105 - 7739347 len: 795242\n",
      "('HG02717#1#JAHAOS010000073.1', 'HG02717.paternal.f1_assembly_v2_genbank.fa', 32364301) target region 5416199 - 6229196 len: 812997\n",
      "('NA19240#2#JAHEOL010000047.1', 'NA19240.maternal.f1_assembly_v2_genbank.fa', 20513270) target region 3530359 - 4367914 len: 837555\n",
      "('chr6_GL000256v2_alt_hg38', 'hg38_tagged.fa', 4929269) target region 3481069 - 3778827 len: 297758\n",
      "('HG00673#2#JAHBBY010000031.1', 'HG00673.maternal.f1_assembly_v2_genbank.fa', 24834737) target region 734801 - 1474176 len: 739375\n",
      "('chr6_GL000255v2_alt_hg38', 'hg38_tagged.fa', 4606388) target region 3387694 - 4178743 len: 791049\n",
      "('NA19240#1#JAHEOM010000153.1', 'NA19240.paternal.f1_assembly_v2_genbank.fa', 5748104) target region 1451541 - 2246056 len: 794515\n",
      "('chr6_apd_hap1_hg19', 'hg19_tagged.fa', 4622290) target region 3447380 - 3761939 len: 314559\n",
      "('chr6_GL000254v2_alt_hg38', 'hg38_tagged.fa', 4827813) target region 3506514 - 4277282 len: 770768\n",
      "('HG01109#2#JAHEOZ010000001.1', 'HG01109.maternal.f1_assembly_v2_genbank.fa', 35424037) target region 27560048 - 28359916 len: 799868\n",
      "('chr6_mcf_hap5_hg19', 'hg19_tagged.fa', 4833398) target region 3512099 - 4282867 len: 770768\n",
      "('chr6_GL000253v2_alt_hg38', 'hg38_tagged.fa', 4677643) target region 3470053 - 3767796 len: 297743\n",
      "('HG02148#2#JAHAMF010000034.1', 'HG02148.maternal.f1_assembly_v2_genbank.fa', 27574774) target region 3417861 - 4208964 len: 791103\n",
      "('HG03486#1#JAHEOQ010000034.1', 'HG03486.paternal.f1_assembly_v2_genbank.fa', 27583779) target region 23355596 - 24150259 len: 794663\n",
      "('HG02109#1#JAHEPG010000124.1', 'HG02109.paternal.f1_assembly_v2_genbank.fa', 5653665) target region 3388333 - 4201417 len: 813084\n",
      "('HG01258#1#JAGYYV010000066.1', 'HG01258.paternal.f1_assembly_v2_genbank.fa', 40483040) target region 27171028 - 27969215 len: 798187\n",
      "('chr6_GL000252v2_alt_hg38', 'hg38_tagged.fa', 4604811) target region 3475869 - 3703833 len: 227964\n",
      "('HG02055#1#JAHEPK010000074.1', 'HG02055.paternal.f1_assembly_v2_genbank.fa', 27564723) target region 3492924 - 4219719 len: 726795\n",
      "('HG02486#2#JAGYVL010000026.1', 'HG02486.maternal.f1_assembly_v2_genbank.fa', 60776310) target region 32112628 - 32907289 len: 794661\n",
      "('HG02109#2#JAHEPF010000055.1', 'HG02109.maternal.f1_assembly_v2_genbank.fa', 9372866) target region 1436526 - 2191206 len: 754680\n",
      "('NA20129#2#JAHEPD010000054.1', 'NA20129.maternal.f1_assembly_v2_genbank.fa', 27332928) target region 23125343 - 23929432 len: 804089\n",
      "('HG02080#1#JAHEOW010000032.1', 'HG02080.paternal.f1_assembly_v2_genbank.fa', 56349618) target region 32203958 - 32997236 len: 793278\n",
      "('HG02572#1#JAHAOW010000052.1', 'HG02572.paternal.f1_assembly_v2_genbank.fa', 8124277) target region 189257 - 979210 len: 789953\n",
      "('HG01106#2#JAHAMB010000019.1', 'HG01106.maternal.f1_assembly_v2_genbank.fa', 123902156) target region 116016100 - 116806378 len: 790278\n",
      "('HG01071#1#JAHBCF010000092.1', 'HG01071.paternal.f1_assembly_v2_genbank.fa', 7272928) target region 6950885 - 7240236 len: 289351\n",
      "('chr6_chm13', 'chm13_tagged.fa', 172126628) target region 32018053 - 32812380 len: 794327\n",
      "('HG03579#2#JAGYVT010000002.1', 'HG03579.maternal.f1_assembly_v2_genbank.fa', 20167018) target region 3492400 - 4287049 len: 794649\n",
      "('HG02723#2#JAHEOT010000107.1', 'HG02723.maternal.f1_assembly_v2_genbank.fa', 27552632) target region 23141054 - 23422570 len: 281516\n",
      "('HG01361#2#JAGYYW010000059.1', 'HG01361.maternal.f1_assembly_v2_genbank.fa', 40042795) target region 32169967 - 32467732 len: 297765\n",
      "('HG02257#2#JAGYVH010000080.1', 'HG02257.maternal.f1_assembly_v2_genbank.fa', 32276632) target region 5412300 - 6222285 len: 809985\n",
      "('HG02717#2#JAHAOR010000061.1', 'HG02717.maternal.f1_assembly_v2_genbank.fa', 31157616) target region 23144624 - 23496198 len: 351574\n",
      "('HG02630#1#JAHAOQ010000088.1', 'HG02630.paternal.f1_assembly_v2_genbank.fa', 20219067) target region 15879189 - 16237163 len: 357974\n"
     ]
    }
   ],
   "source": [
    "rgn_lengths = []\n",
    "for k in list(aln_range.keys()):\n",
    "    b, e = aln_range[k][0][0:2]\n",
    "    print(sdb.seq_info[k], \"target region\", b, \"-\", e, \"len:\", e-b )\n",
    "    rgn_lengths.append(e-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'count')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATc0lEQVR4nO3dfbBkdX3n8feHGRREEMhcyAjo1V2KDbEikFmjYmUpFYOiQLk+QMXUJJiQrfUhZrNljZuH1bK2MptHjaZ2pRJxiEZlRWUWK1EyhqgxcRwQEFB2iE5c1nFmfEg0hiKC3/3jnHF67ty5NJd7uvvye7+qTvU5vz6nz7fvw6dP/87pX6eqkCS144hpFyBJmiyDX5IaY/BLUmMMfklqjMEvSY1ZO+0CxrFu3bqan5+fdhmStKrcdNNNX6+quYXtqyL45+fn2bFjx7TLkKRVJcnfL9ZuV48kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDVmVXxyV5JWk/lNH3nQdXZtvnAClSzOI35JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGjN48CdZk+RzSa7vl09MckOSnf3tCUPXIEk6YBJH/L8EfGFkeROwrapOB7b1y5KkCRk0+JOcClwI/NFI88XAln5+C3DJkDVIkg429BH/W4DXA98faTu5qnYD9LcnLbZhkiuS7EiyY9++fQOXKUntGCz4k7wQ2FtVNy1n+6q6sqo2VNWGubm5Fa5Oktq1dsDHPhe4KMkLgKOA45K8G9iTZH1V7U6yHtg7YA2SpAUGO+KvqjdU1alVNQ9cCny8ql4BbAU29qttBK4bqgZJ0qGmcR3/ZuD8JDuB8/tlSdKEDNnV8wNVdSNwYz//DeA5k9ivJOlQfnJXkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0ZLPiTHJVke5Jbk9yR5E19+4lJbkiys789YagaJEmHGvKI/z7g2VX1VOAs4IIkTwc2Aduq6nRgW78sSZqQwYK/Ov/ULx7ZTwVcDGzp27cAlwxVgyTpUIP28SdZk+QWYC9wQ1V9Bji5qnYD9LcnHWbbK5LsSLJj3759Q5YpSU0ZNPir6oGqOgs4FXhakqc8hG2vrKoNVbVhbm5usBolqTUTuaqnqv4BuBG4ANiTZD1Af7t3EjVIkjpDXtUzl+T4fv5o4LnAF4GtwMZ+tY3AdUPVIEk61NoBH3s9sCXJGroXmGuq6vokfwNck+SVwFeAlw5YgyRpgcGCv6puA85epP0bwHOG2q8kaWl+cleSGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhozVvAn2TZOmyRp9i05ZEOSo4DHAOv6r0hMf9dxwOMHrk2SNIAHG6vnF4HX0YX8TRwI/m8DfzhcWZKkoSwZ/FX1VuCtSV5TVW+bUE2SpAGNNTpnVb0tyTOB+dFtqurqgeqSJA1krOBP8ifAvwJuAR7omwsw+CVplRl3PP4NwJlVVUMWI0ka3rjX8d8O/PCQhUiSJmPcI/51wJ1JtgP37W+sqosGqUqSNJhxg/+NQxYhSZqcca/q+auhC5EkTca4V/V8h+4qHoBHAUcC362q44YqTJI0jHGP+I8dXU5yCfC0IQqSJA1rWaNzVtWHgWevbCmSpEkYt6vnxSOLR9Bd1+81/ZK0Co17Vc+LRubvB3YBF694NZKkwY3bx/9zQxciSZqMcb+I5dQkH0qyN8meJNcmOXXo4iRJK2/ck7tXAVvpxuU/BfjffZskaZUZN/jnquqqqrq/n94FzA1YlyRpIOMG/9eTvCLJmn56BfCNIQuTJA1j3OC/HHgZ8DVgN/ASwBO+krQKjXs555uBjVX1LYAkJwK/Q/eCIElaRcY94v+x/aEPUFXfBM4epiRJ0pDGDf4jkpywf6E/4h/33YIkaYaMG96/C3w6yQfohmp4GfDfBqtKkh7h5jd9ZKz1dm2+cMX3Pe4nd69OsoNuYLYAL66qO1e8GknS4MburumDfuywT3IacDXdd/V+H7iyqt7adxO9H5inG/PnZaPnDyRJw1rWsMxjuh/4lar6EeDpwKuSnAlsArZV1enAtn5ZkjQhgwV/Ve2uqpv7+e8AX6Ab7uFiYEu/2hbgkqFqkCQdasgj/h9IMk93+edngJOrajd0Lw7ASYfZ5ookO5Ls2Ldv3yTKlKQmDB78SR4LXAu8rqq+Pe52VXVlVW2oqg1zcw4LJEkrZdDgT3IkXei/p6o+2DfvSbK+v389sHfIGiRJBxss+JME+GPgC1X1eyN3bQU29vMbgeuGqkGSdKghP317LvAzwOeT3NK3/RdgM3BNklcCXwFeOmANkqQFBgv+qvoU3Ye9FvOcofYrSVraRK7qkSTNDoNfkhpj8EtSYxxaWcsyzsiCQ4wqKOnh84hfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSY9ZOuwBpHPObPjLWers2XzhwJdLq5xG/JDXG4Jekxhj8ktQYg1+SGmPwS1JjBgv+JO9MsjfJ7SNtJya5IcnO/vaEofYvSVrckEf87wIuWNC2CdhWVacD2/plSdIEDRb8VfUJ4JsLmi8GtvTzW4BLhtq/JGlxk+7jP7mqdgP0tycdbsUkVyTZkWTHvn37JlagJD3SzezJ3aq6sqo2VNWGubm5aZcjSY8Ykw7+PUnWA/S3eye8f0lq3qSDfyuwsZ/fCFw34f1LUvOGvJzzvcDfAGckuSfJK4HNwPlJdgLn98uSpAkabHTOqrrsMHc9Z6h9SpIe3Mye3JUkDcPgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSY9ZOuwBJWi3mN31k2iWsCI/4JakxBr8kNcbgl6TG2MevwYzbH7pr84UDVyIt7ZHSdz8uj/glqTEGvyQ1xuCXpMYY/JLUGE/uPgJ4ElU6vNZO3I7DI35JaozBL0mNMfglqTGP+D7+1d7/bf+khjTO39es/m9o+Tzil6TGTCX4k1yQ5K4kdyfZNI0aJKlVEw/+JGuAPwSeD5wJXJbkzEnXIUmtmsYR/9OAu6vqS1X1L8D7gIunUIckNSlVNdkdJi8BLqiqn++Xfwb4iap69YL1rgCu6BfPAO6aaKEHWwd8fYr7X4q1LY+1LY+1Lc+0antiVc0tbJzGVT1ZpO2QV5+quhK4cvhyHlySHVW1Ydp1LMbalsfalsfalmfWaptGV889wGkjy6cCX51CHZLUpGkE/2eB05M8KcmjgEuBrVOoQ5KaNPGunqq6P8mrgY8Ca4B3VtUdk67jIZqJLqfDsLblsbblsbblmanaJn5yV5I0XX5yV5IaY/BLUmuq6hE7Ae8E9gK3j7SdCNwA7OxvTxi57w3A3XSfGfipkfYfBz7f3/cHHOgiezTw/r79M8D8yDYb+33sBDYuqOs04C+BLwB3AL80Q7UdBWwHbu1re9Os1Dayzhrgc8D1s1QbsKt/zFuAHTNW2/HAB4Av0v3dPWOGajuj/5ntn74NvG4W6gN+me7/4HbgvXT/H1Ov62Fn40o+2KxNwE8C53Bw8P8WsKmf3wT8937+TLqwezTwJODvgDX9fdv7f5QAfwY8v2//j8D/7OcvBd4/8s/+pf72hH5+9I9jPXBOP38s8H/6/c9CbQEe288f2f8xPn0Wahup8T8Bf8qB4J+J2uiCf92CtlmpbQvw8/38o+heCGaitgV1rgG+Bjxx2vUBpwBfBo7ul68Bfnbada1INk4zmCcxAfMcHPx3Aev7+fXAXf38G4A3jKz30f4XtR744kj7ZcA7Rtfp59fSfTIvo+v0970DuGyJGq8Dzp+12oDHADcDPzErtdF97mMb8GwOBP+s1LaLQ4N/6rUBx9EFWGattkX+5p4H/PUs1EcX/P+XLnzXAtf39c3cz+2hTi328Z9cVbsB+tuT+vb9v+T97unbTunnF7YftE1V3Q/8I/BDSzzWIZLMA2fTHVnPRG1J1iS5ha6b7IaqmpnagLcArwe+P9I2K7UV8LEkN/VDjsxKbU8G9gFXJflckj9KcsyM1LbQpXRdKky7vqr6f8DvAF8BdgP/WFUfm3ZdK6HF4D+cww0lsdQQE8vZ5sAOk8cC1wKvq6pvz0ptVfVAVZ1Fd3T9tCRPmYXakrwQ2FtVNy1Rz1Rq651bVefQjTz7qiQ/OSO1raXr8vwfVXU28F26LopZqO3ATrsPdF4E/K8laptYfUlOoBtA8knA44Fjkrxi2nWthBaDf0+S9QD97d6+/XBDSdzTzy9sP2ibJGuBxwHfXOKxfiDJkXSh/56q+uAs1bZfVf0DcCNwwYzUdi5wUZJddKO6PjvJu2ekNqrqq/3tXuBDdCPRzkJt9wD39O/coDvJe86M1Dbq+cDNVbWnX552fc8FvlxV+6rqe8AHgWfOQF0P30r1Gc3qxKF9/L/NwSdmfquf/1EOPjHzJQ6cmPks3QnO/SdmXtC3v4qDT8xc08+fSNenekI/fRk4caSGAFcDb1lQ6yzUNgcc388fDXwSeOEs1LbgZ3UeB/r4p14bcAxw7Mj8p+leMKdeW7/OJ4Ez+vk39nXNRG0jNb4P+LlZ+X+gO7d1B925rtCdIH/NtOtakVycZigPPdH1Fe4Gvkf3CvpKuv6zbXSXSG3j4H/eX6U7E38X/Vn3vn0D3eVcfwe8nQOXYh1F97b0brqz9k8e2ebyvv3u0T/m/r5n0b1tu40Dl7C9YEZq+zG6SyVv6x/3N/r2qde2oM7zOBD8U6+Nrh/9Vg5cBvurs1Jbf/9ZwI7+9/phujCZidr6dR4DfAN43Ejb1OsD3kR3CeztwJ/QhfrU63q4k0M2SFJjWuzjl6SmGfyS1BiDX5IaY/BLUmMMfklqjMGvwST5pwXLP5vk7f38G5P858NsN5fke0l+cYnHPjLJ5iQ7k9yeZHuS5/f37UqybiWfS/+45yW5fuFzGWA/nx7icfvHvijJUp/aVQMMfs2ilwJ/SzdQ1eG8mW7wq6dU1VOAF9GNdDoz+k9iPmRV9cyVrmXksbdW1eahHl+rg8GvWXQZ8CvAqUkOGZgqyWOAXwBeU1X3AVTVnqq6ZpF1P9wPmnbH/oHT+kHo3tW/U/h8kl/u21+b5M4ktyV533IK79/JXJnkY8DV/buXa5N8tp/O7debS3JDkpuTvCPJ3+9/l7L/nVI6vz1S58v79vOS3JjkA0m+mOQ9SQ4Z22Wx57PgXdctI9O9Sf5dkmOSvLOv9XNJLl7Oz0GzbeJftq6mHN2P8rnficDWpTZIchrww1W1Pck1wMuB31uw2r8GvlJLD2y33+VV9c0kRwOfTXIt3TAep/TvFEhyfL/uJuBJVXXfSNty/DjwrKq6N8mfAr9fVZ9K8gS6YXh/BPivwMer6jeTXABcscjjvJjuE7dPBdb19X+iv+9suiECvgr8Nd04Rp9asP2Sz6e6gfhI8iK6EU8/TfdJ1Y9X1eX9NtuT/EVVfXdZPwnNJI/4NaR7q+qs/RPwG2NscyndF15AN3bLUt0943htklvpuo5OA06nG0PlyUne1ofu/heQ24D3pBuB8f6Hsc+tVXVvP/9c4O39C+BW4Lgkx9IN2/E+gKr6c+BbizzOs4D3Vjda6h7gr4B/29+3varuqarv0w35Mb/I9g/6fJKcTjf2zMurG4jsecCmvt4b6YYUeMLYz1yrgkf8mjWXAScn+el++fFJTq+qnSPr3A08IcmxVfWdwz1QkvPogvcZVfXPSW4EjqqqbyV5KvBTdINkvYxuXJQL6b617SLg15P8aHVjpD9Uo0fHR/T7v3d0hcW6ZhZ7Ckvcd9/I/AMs/r98yPNZUMMxdC+yv1D9yKL9Pv99Vd01Rn1apTzi18xIcgZwTFWdUlXzVTUP/Cbdu4AfqKp/Bv4Y+IN0Y7iTZH0OHSv9ccC3+tD/N3SjI9L3pR9RVdcCvw6ck+QI4LSq+ku6bo/jgceuwNP6GPDqked4Vj/7KboXHJI8j27QtIU+Aby8PycxRxfi28fZ6ZjP5yrgqqr65EjbR4HX7H9hSnL2OPvT6mLwa5p+Lck9+ye6o/0PLVjnWhbv7vk1um+VujPJ7XQjTu5bsM6fA2uT3EZ3FdDf9u2nADf23RnvovvKvDXAu5N8nm500t+v7vsIHq7XAhv6E6x3Av+hb38T8LwkN9ONQ78bWPju5UN03TW3Ah8HXl9VXxtzv0s+nyRPBF4CXD5ygncD3c/pSOC2/uf65of8jDXzHJ1TmoIkjwYeqKr7kzyD7tuxzppyWWqEffzSdDwBuKbvkvkXustTpYnwiF+SGmMfvyQ1xuCXpMYY/JLUGINfkhpj8EtSY/4/t1uuvKuFzHUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(rgn_lengths, bins=32);\n",
    "plt.xlabel(\"HLA Class II region size\")\n",
    "plt.ylabel(\"count\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some mapped regions are larger, let's check the repeat / CNV by quick dot plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABr0UlEQVR4nO29d3xc1Z33/z53+qjMqEtWl+Ui925jAzYG0zskQEJJICE92fIkm2w2m23P89vsZndTNpuElE0gJITQQzfFgAHbuHfLsqze26hMuXPvPb8/ZkbIxja2GVnF5/166aXRmbn3nnskzWe+59uElBKFQqFQKE6GNtYTUCgUCsX4RgmFQqFQKE6JEgqFQqFQnBIlFAqFQqE4JUooFAqFQnFK7GM9gWSTnZ0ty8rKxnoaCoVCMaHYtm1bl5Qy50TPTTqhKCsrY+vWrWM9DYVCoZhQCCHqT/ac2npSKBQKxSlRQqFQKBSKU6KEQqFQKBSnRAmFQqFQKE6JEgqFQqFQnBIlFAqFQqE4JUooFAqFQnFKlFAoFArFBMY6B60ilFAoFArFBEU3TRrC+qiLhRIKhUKhmGBYUqKbJlv6g0xx2tGEAMCIREblekooFAqFYgKhmya1wTBNkSjL0r04bTYABns6aX7n5VERCyUUCoVCMUEwLItNgSEsKSlxO3HabERDIfq7O6h+5ud0l1ahOZ1Jv+6kKwqoUCgUk4mR/ge7prHSn4omBJau093dSe2bz9GTmk6j30Xqew9TlffXeFJSkzoHJRQKhUIxTjEsi7pQbCtJE4Iyjwu7ptHd2kzbe6+x+fA+Wlx5FGutpIeiXHr1l5MuEnCaW09CiL8UQuwTQuwVQvxBCOEWQmQKIdYLIQ7Hv2eMeP23hBA1QohDQogrRowvFkLsiT/3IyFiHhghhEsI8cf4+GYhRNmIY+6JX+OwEOKeJN67QqFQjFtGikSZx0WZx4UmBI2Hq/nDE7/n1Q6Lo2EDEd7BmrJLuGTevaSYqUgr+RFQHyoUQohC4KvAEinlHMAG3A58E3hVSjkNeDX+M0KIWfHnZwNXAv8jhLDFT/dT4H5gWvzryvj4fUCvlLIS+C/ge/FzZQLfBZYDy4DvjhQkhUKhmGxYUmJJSUNYBxi2Iixdp6etlSdfeIyerkZsA9vId8FluRfgd3modm6n2d6KFGMgFHHsgEcIYQe8QAtwA/Db+PO/BW6MP74BeERKGZFSHgVqgGVCiAIgXUr5rpRSAg8ed0ziXI8Bl8atjSuA9VLKHillL7Ce98VFoVAoJhVhwzhGICq8bmQkQqCrk9rXn6bujafp767DZQuxevXNLPGXIvNS6S7wU5m5AtE9BKNgUXyoj0JK2SyE+D7QAISAl6WULwsh8qSUrfHXtAohcuOHFAKbRpyiKT4WjT8+fjxxTGP8XIYQIgBkjRw/wTHDCCHuJ2apUFJS8mG3pFAoFOMKw7IwLItnu/q5Njt9OC+ir6ON/S8/weEB0NMyKMyHjKw85rpzSesP03lhJTnubGz98Ooz/x/OZVdSPHUOGrYPueKZ8aFCEd/quQEoB/qAPwkh7jzVIScYk6cYP9tj3h+Q8gHgAYAlS5aMfj67QqFQJAHdNAHYFBhiisvBtdnpuO2xt+WO5kaefvxh2nUvqanNpAzWMyf9ajJSdFrzIphpXQR21WI072dHOAtH+gBFvQFk1AC7I6nzPJ2op8uAo1LKTgAhxBPASqBdCFEQtyYKgI7465uA4hHHFxHbqmqKPz5+fOQxTfHtLR/QEx9fc9wxG0735hQKhWK8EjYMnuroY1G6lxW+FOyaRnRoiCHLAuDPTz5Dz0AEu6hmVeFllC++CI+WSnX6EIM7NyA3BqhvaKdLpjN1Tjph70XMWLkOp8eT9Lmejo+iAVghhPDG/QaXAgeAZ4BEFNI9wNPxx88At8cjmcqJOa23xLepBoQQK+Lnufu4YxLnuhV4Le7HeAm4XAiREbdsLo+PKRQKxYTFsCxadIMFaR4qvG6cNhv64ABv/+qnPP2HPzA41E/Y2UqqPcDaeZ+gwJeLXltNK+20vPEWYtN+dtYGsNr7MGQ/M4pv4ONXXUNGVtaozPd0fBSbhRCPAdsBA9hBbJsnFXhUCHEfMTH5WPz1+4QQjwL746//kpTSjJ/uC8BvAA/wQvwL4FfAQ0KIGmKWxO3xc/UIIf4ZeC/+un+SUvZ8pDtWKBSKMcCwrFiinJQ0RaKUuJ1oQqAJQai/jwPbt7M+bJHvqsPt9nDp9BVYQ0H6ZQ2D5cvQo91Ea1voOtTJQDCHtBSDoYp5XDF/LktWLcTtdY3a3IU8ByVqzyVLliyRW7duHetpKBQKxTC6abIpMES+045d0yhxO7FrGr2tLWAT7Hn0f9gT9tITPsBFU1ex6vq76NjxNr32ITqsbjKGogzV93HQ6Kdv82FSfDNYeNNqKmbNJ9Ofid3x0XOnhRDbpJRLTvScysxWKBSKUSKRE9GiG8N+CIhlWXccreWZP30f0qfS25mHN6OWKe4SSnJLEHYbwakFHN3yJ8wdu9k6WMSQI51FC7Npmepl+YUrWbF67Tm7DyUUCoVCMQoYlkVtMHyMBQEgLYvejg52vfI8ewezKdDb0HwBLl17Jy0vbiBakYvUJO1bXuLI2420h9Pw0IHT6yB76jo+scTL1mCUSDSKy5Hc6KaToYRCoVAokkyioVBLJMpKfyp2TSM80E/UMGmtP8pb725Cd7sp0YKUlhSQa87A2TVA9fwIvrpddJPJxreOMhAyqJAlZC6ey/LLl+PPzufPNXUUp3qwaeeu+LcSCoVCoUgSlpQYlsWW/iDL0r2UeVzYgEBHC7sf/SmHZDbYbRieNhbnTKGzWzCjoILGjAF6nC7ynniZN8IruOaOeQitn4zsLBZdsZaqBSuxOWwc6OrhhunlOO127LYPJtVJy0KMgoAooVAoFIokEIxGaYvGAjwTDYXC/f001B+hbdOfeb0rl/T8RlJMyYr8ZWRmpXJosJPItBJCu9dTV/0e1cELcEc68GXnMGfRbcxeWkluURkHunpBSubkZZ90u8nQo/R3d5KRl590sVBCoVAoFB8BS0qC0SgPtfVyT0EmbrsdGTXoaW/n3d//ilZvJvqQICV1LxcWr6F0wQqM+jo8FTNp2/97jm7byOEtR2kLOigM2vDOuhB/Rh4XXZPFvp5eNh0+yvXTTm5FGHoUy5Lsq20iPTMLP+KEJS0+CkooFAqF4ixIRDTVhSK0RKLclZ+B1+FgsLuLAxu3sq32AEeDdkq0alJz/awp/BjpTjDq6hgsy2b/m4+StaeP/fp2+iODFHhziBQW8bFPfgy7y8G+tj5smsb108rxuj6YIyEti1AwxNade0n1Z5CWlUtpdhqalmyZUEKhUCgUZ8zIbaZErwgrFKKztYXdT/ycTWEf3pRWsrFYufIGfMEQ9f07sC25DNOwiOx/l+bNm2nuL8cpmyl2zWLRjZdQMGM27RJygPn5OQAntCLCQ0H6ujp4+1Ars2dOY2q+H5vdhpTWqNyvEgqFQqE4TSwpCRsGv23tGd5m0oQg0NbMW797gAGHYKv0kWM7wAWll2LkdZBnCALTiunbv4v0/e/SUVPL3hbIcpZT6kqlqbSSdddfTnZRCb/cU81n5804oTgkGOjv56mX32X58iVcc2EBTrcTTRNEwhGOtNUzo3gqtlMcfzYooVAoFIrTQDdNmiJRAO4pyMTrcCAti87mBjY89D02UUVFaCdl6TYWllxIXkomh/pa6cn3Mnj4AOHXd/KyVsqQ5cStd6BPW0zOghyuv/QKnG4X1T19fHbeDNJOUtRPWhaRcJT1m/ezaMkCKgoysNs1LMMiouvsOlTN/BnTky4SoIRCoVAoTolummhCDIe82jUNTQjCg0M01Bxi/Ya3OBqejSE2k1I0izWLLiOwZwfBxXk0D2wi+vrjHKmO0kghDtlCnpFN+aJrWHXlxVQPDHGgL0D9wNBJfRGWYWBZkr6uDgYd6Vy9aiFOtxMpLcLBMK11AYZcfcyfMR2Xe3TqPSmhUCgUipMQNgye7ern+hzfcAkOIxKhZt9Ojhw+QlugH5naSHq0ibmOMqpyZ2D0tnK4Ioh/80t4dr7HVmM+mugj15hN5aK1LL9oNpm5RZhC8m5NE5+bP4NFU/JPGPaqB0PUV+9j0OEjPSuX0owU7HaNwcFBalrrqGuvY/a06czIqsCR5B4UI1FCoVAoFMdhxYultkVNrs1Ox2mzYRkGDYcP0rn1DQ4c2EkgPQ0tx81FhRfgnGVx5MWdRK7KYe/h9WjbDvNeaxqD+gx8oo/UrDIuufVySqfOxOawsbejG6Tkc/NPvtUUCYbZtH0PIqOQ5VNzcTjtCCAQCPDHzU9z2/JryU/vIiVFO6VPIxkooVAoFIoRWFLSENYpcTuHS4FbhsGbTz3K1poDFLoK2O4YZK6/kiXF00gXTroON9C7NpvUuqM0bGinT2aihZuptM1hyiUruODSq7B7vMckztm0D77BW4aBjGdBHKxvpaC8krIcH06nDSNq0NsxyOa6d7h+3nzsjgA5qSvRNDtCjG45DyUUCoVCQayIX8L/kBAIgKPVB/B4Uni7upY24xB53hDzzWksLSnH70wlUJZJXWQT2san2TlwAX393aTbsimZfzOZRV5WX3kDUdPk2Zp6rq0sPWninKnrHNqzg6gniznTy5kzrQShaRhGFF03qWmpo7KghCV2L6kpKXjcxWjauXkLV0KhUCjOe3TTZEt/cNgPoQmBEYlwYN8eHnnlTa5aMhPD8S6LXOVcfNlt7P/9M4jpZWTkTKX93d8T3vIujZHlBDs6cfkzqZy6lLWfuAnsGrvbu0BKrq0sPaGzGmLbTF2dbbzcrPGZtYXY7DEhGQoO8erujRRnFjKrtBKny01uzoUIoY26FTESJRQKheK8JJET4bbbadGN4YgmgFB/gLo3X+StQw0YkZ1k5V1IpWclszy5eFwZLL3/Tlq6j7Bz30Y2bDhCiz6Daf0pDOWncPMn76CwYjb7u3vZ2NrF5+bPwON0njRxTiJ44Z1dzK2azmcvLSAlJSYm4UiYTXXvccmcFSC7MMxWnLL0nFkRI1FCoVAozjssKdnd18+fuwf569JcStxORNyBPdjVTd2GF+m0UhD+AAusqRRk5DNVuuhemok7fJi3n/sJQy0VhKIWkWAL0yI59GXB/V/6Cll5hcPbTPPyc05ZfuPhx19i4fxZXLVyPi6PGzNqYJqxjO+a5loWFuSjaT143KXn3IoYiRIKhUJx3pCIZtJNk+qwwVeKsvA6HERDIeoPH6SorJJtjz/OZtlHcWoTFxSspLv9EHaXh9olKWjbHqW+tYq68FScfa2INA8LVt1JZo6d0rmL6ZA2Ojp7TrnNZBkGfZ0d9NlSueiSi5ian4HDaSccDPPmph3kFqYxZ+oMynI9eL3FYyoQCZRQKBSK8wLDsmgI61hSognB9Tk+nIARibDjmYd5vbmLO25Kpd7XwQzpZdXqv2LoyCGqLy/G1AyCrz7KkaHF9Bu92CMwpXQuq6cvY9a6VehY/NuOA/ztkjmkut0n3GayDAOAvs4O/Dm5+DUbECviFwqFqG6txfR1U5qXhc1mIyWldMwFIoESCoVCMakxLAtNCJoi0WOimfSBft555UVmrVjNOwf70G01OBwO1s1cgt4rcfX20lWahXzpl+ztzePIwHy0rl1k58zGn17OVbdcTnFRGfu6e6kN9PM3i2bhT0n5wPWlZWHoUZqq95GWmc2Aw49fs6FpgsHBQWw2G6/ufovVs5aS6TpIWmrJuBGIBEooFArFpCVsGGwdCLHClzIsEtKyCPX38cwP/oW9rkJE2nsYmRuYlbaaDF8WBItoyIU93YdoeG47TfXp9Id2kTEwSEbuWq647yZ0h4uA3cPBnj7m5+cwPz/npJnVgZ4umvt1UrPLychNJ1PTEMBA/wA/fPlnrKu6hAsqi9BsfRTkX4Pd7j73C/UhKKFQKBSTjkRL0q0DIZakeYajmaKhEC11NdRseYtd0em02TdiT1nAZVPvxOrow67ZsIrLCe9/A2NvD0cPdCL6LOZeVsXWAQdX3XELOcXFPF19lGsrp5w0J8IyDPSIzusbNjJt/mJmVxZgs8esCCNq0N8dAim5a9lc/L40bHb7Oc2LOFPG56wUCoXiLEn4IuD9lqQAQ33d1L31EltaB+g1BYbnRZZGZlJckE+wtZ3GZUWYDjjwxsO8veEILSkevM4QKSVLKJt7IWXzV9FmSPI1jRuml5+0JWk0HKHu4B767el4SmdRkuvH4bSh67E5VTfX4k13UJxehBXwkppafk6yqz8KSigUCsWkwZLyGF+EJgRDgT7sNjuHnn6ALXo6ZrSWFClZahRidbbhcaXwbkE3ts1v0pc2lz9t6yJjqIncjCKmLr2DqmWzyMj0x62IUuw22wmtCGlZmIbJ7iNNDDkLWV6Rjc1uQwhJKBRi/a43Kc4spKqkAtNqx253kJm5ctxaESMZ/zNUKBSK0yARzTTSF9HZ3MjG537PqmvuYF/Ai8vdx7ILbiTDn8nrjz1A7XV+rvR60d/cx6EuJ8XTakjzNOOduYKbbrwRp8fHbw/W8vkpBae0IizDoKetBaFp+LOymZuZitNpY3BwkKMdjUgpWTtnBVJ2YVrteNzjI+z1dFFCoVAoJiQJYUg8ThTy04RgsLeH1vp6tr39Akd6gswNDzDoaeHKJXfgCfYSDNTRVhFlxoG3aMy+gd3dDlJ6e0nP9HHn2r8gvSQXb2oqezu6+Oy8GSfNiUhYEb0dbbT068yuLCXDbscyTQYHB/n9O09w+4rrgU4i+l78vqXjfpvpRCihUCgUEw5LSupCEco8ruEtphK3E0vX6e/r5eXHfkmnmUZzoItct0l6io915StJsfXjmV7Frl1/wvHeQeq05fS+sxmPMUjW1PnkF8/A6/EiRaxX9Zzc7JOW8Db0KL0dbTT366RlZjOrMgV7PHFu+8H9eFMd3DB/PjZ7Lx53OVA+IbaZTsTEnLVCoTgvSWRWJ76PxAyHqV7/OIGeAaq7g2TlDlCeJ1hcsIwMXyaN6RGi5gCdz/2MzuYeesLTSO/w45+axrwLbqFizXLS0tIwTJO9HV2nFAk9GOJAfRupGZnMrkzBZrdhRg3CwTA7qw+RWeAi32vHMDpwu6omrEAkmNizVygU5w0JK+J4X4QZjSIti5q9u3n2UADTPUCer5+1a+8lun0PfS4LA0mfPkTT+vfYMehgYBDyQ258+VlcsPYK/LkFvFjfzA0ZflwOxwlFwjIMDMPCpkHN1nepWnIBdrc7ljgXGOTNLTvw+VJYUFVJOLqHtNSlCDF1wosEKKFQKBTjnIQwaEJQ5nnfVyCkZCjQR/eezQhpsWlPLVOsbtKmZFMhL8TV3sPR+dl0H9nC9EAfwT/+mfc82fijdtLsLoqypzHztqVMKS1nb0cX11aWDjurTyQSB3a+x6FuuHb1YmasuBCb04m0JMGhIH/etJ4BrZul5ReArRu/dyk2m/OcrtNoooRCoVCMW3TTpEU3jgl3hXiTn907OVKzk8ZQNtddugDtwBssufhCPFE7nT4ngfQ0ut98hr49RziacoA9spAM+xSWXL6Q6VWziUobf2hsZ5phnHKbKTQUQrPZ2R7wcuPKSpzumAAYUYO+riHaQq1ctWwVXT2P4/eXY7e7J5yz+sNQQqFQKMYlhmWxpT/IsnTvsEBIy2KoP8DOdzayaXct2VoW3ilt+LNyuea6L1PXXAMlWTTvehXzjcMcqtewBj1gE7jTCrnp3pspLC0nqOt8f8tuZqR6TtiSdPhaQ0F+8+d3uO+Gi/j4qlm43A7MeERTQ3cz3nQHRV6BaR6mpPguHA7vuV6mc4ISCoVCMS6xa9pwxzmIVXk9tH0r7+17j6YOE7ttPxWX3s6CmdfjdbjoP3oYd14Kbds20fbqJpq1Qvy9uzEzVpFXMZ17/mo+qX4/ezq6QUq+s2weHqfzhLkRiYim1oEon7puJZ4UD9KS6LrOziN7eXnfej5z0VWketx4PeXA5PBFnIzJe2cKhWJCkajP5Bzx6d6uaUjLIhIKceSdl3l8VzNZLh271cZly9dRlDsVj8uNieSIvZnuh16nftBLn24j3TlIuOxirr/pKoqKSohEoxzo6mFObtYJrQhpWQhNQw/rHDjaTGpGJlVTU3A47UhL0tHSR6fezuySCvy2VLIyKyflNtOJUEKhUCjGnGA0SlMkys6BEDfm+ofFor+7g96uTprb2zjSBmKwnvKSBUxftwRbfRuHo/XsP3qEFfnL6d90mNeiqUwhAOkVXP3J6/BlF9JuacMiUZWdeUILwtR1ejrbycjJ51BdC1XlhdidDgQQCoWwO+zs7d3KwqJiEL2Ul3960m4znYjJL4UKhWJcE4xG+W1rD0Uux7BISEvS09HCpt/9E2+9sZHAQIiFqyu4ceEsZudkMXiojR1pPRSTg/OlDfR2dbK53SIz2E5qxQLu/OxdlM+Yz1vdA+imgU3TmJOb/QGRsAwDyzCo2beLrQ0BJIJZlcU43U4E0NrUxYNv/BE9EmZuvhtfegVeT+l5JRKgLAqFQjFGJBoKdRgW9xRk4nU4YrkKEZ1gRz/b9+/hteB0it2tXLH0RtJTPOwcehXbgrnUNh/Cs2krR7u2M1RVhRQWmXadWcuuZ8nlV1ATDNPX2cMN08tPWgo8PDDAkQP7mDZvEZ3eEtbMiLUkBTBNE6EJwim9XFHlwe124fFcMKn9EKfitCwKIYRfCPGYEOKgEOKAEOICIUSmEGK9EOJw/HvGiNd/SwhRI4Q4JIS4YsT4YiHEnvhzPxIiFsoghHAJIf4YH98shCgbccw98WscFkLck8R7VygUY4BhWRiWxabAEJaUlLideGw2TD3Ke688T/u2jeys2Y3UBcWuTq5YegP+zEzqw63smRJiqK2Wyn1HeKWliCYrRGpjOw5T8onPfZkLb7iJlPR0bJqNOXnZeF2uE/oiIkND/PHxl/jl1m4i0SgrpubgdsesDdM0OdhwGMs0yfM4KCq8AZvNed6KBJy+RfFD4EUp5a1CCCfgBf4WeFVK+a9CiG8C3wT+RggxC7gdmA1MAV4RQkyXUprAT4H7gU3A88CVwAvAfUCvlLJSCHE78D3gNiFEJvBdYAkggW1CiGeklL1JuXuFQnHOsKREN83hjnOJiCZpWVRveYeo6eLp3TWsXLqKFHc/pTOKmJF7CTn5DtoGm8gxUyncupO3uhoJN/VTkeOhYMnVLL/yWtzpPkLRKE9XH2WqP435+bkn9EVYhkFvexu6ZeKesZC/q8ohLS0Fy7Dobh8gKy/Ww7o0x4XNZsPrGT99q8eSDxUKIUQ6cDHwKQAppQ7oQogbgDXxl/0W2AD8DXAD8IiUMgIcFULUAMuEEHVAupTy3fh5HwRuJCYUNwD/ED/XY8B/x62NK4D1Usqe+DHriYnLH87+lhUKxbnGsCzqQhFaIlGWpXuHQ171YJCGo9U8tekRbL5sjMEWQoMFrFl2HU63k1qbZFvbfqpSSnj75Td490g2eC0KcrJYfM1qFl98GZYQmFJS3dN3yq0myzDo6WgnYEtD2DVuWOjBadeQlqSrI8Ceru1cnHMxDruDlBQlECM5HYuiAugE/lcIMR/YBnwNyJNStgJIKVuFELnx1xcSsxgSNMXHovHHx48njmmMn8sQQgSArJHjJzhmGCHE/cQsFUpKSk7jlhQKxbnAknK4mVCZx0WZxzUsEt2tjTS98QSH6lo4auh8LA+GZmdyyeorcXvc1PUdZdPmx8msDrJjig6ZKRQIjdKZF7Nw7RIyi0rY2tpO/cAQN82YetLsakOPYrPb6OloZ8CWRmlWKpom0IciHNrfwbS5UxhydzMn345NiyX2KZE4ltNZDTuwCPiplHIhMERsm+lkiBOMyVOMn+0x7w9I+YCUcomUcklOTs4ppqZQKM4VCSuiIaxT5HJg1zTsmkagq5OejlZeeOy/eaa7mzqvj4W2FPLKLsTfmILNsGjsrSda14144zDbu9OJ9uxlXt50Fq1bzaV3XIdvShFPVx8F4PppsYZCJ/JFhIeCbN66i6huMOhIpzgrFSEkkWCYp19+hQeP/pxgZIhct43srPPXWf1hnM6qNAFNUsrN8Z8fIyYU7UKIgrg1UQB0jHh98Yjji4CW+HjRCcZHHtMkhLADPqAnPr7muGM2nNadKRSKMSHhi2iLmsfUaIpVeN3Ju68+SlrxbPYMBblxYQ5hPYPg9l6KKqrY56inNtiIdvgwdYdrqI+2kuF0Urb8UxTOX45f19nb209toJ9rK0vxOJ0ntCJMXae7rZXqxg5yS8qxOx2UZDkwjCiHmo4QMSI4iw9zh9eJx+XAbp/4W03SkgjtRJ+tPzofujJSyjagUQgxIz50KbAfeAZIRCHdAzwdf/wMcHs8kqkcmAZsiW9TDQghVsT9D3cfd0ziXLcCr0kpJfAScLkQIiMeVXV5fEyhUIxDdNOkZijEM50B8h027JqGkBJ9KMi+t9/gpaceYn9wEFu4k8Ul6fhSF1KeNgX7jcto796D753f0fPaG0QsH2agk2mhqVx+86eYuXgVQdPkl/tqqMrO4KYZU0nzeE5oRUTDEeoP7mPQ6WPZ4vlUFmahaQLDiFLTUsf0wnLmlk5j0ZRCpld+DYfDO6FFQloSy7Awe8NI64N9OpLB6dpZXwEejkc81QKfJiYyjwoh7gMagI8BSCn3CSEeJSYmBvCleMQTwBeA3wAeYk7sF+LjvwIeiju+e4hFTSGl7BFC/DPwXvx1/5RwbCsUivGDbsb+xRNF/Mo8Lpw2G5Zh0HDoAIE929jY0Ue9PUql3c+85Yto33aE4pnzaAg2UvPsjzjUdRmtQwsoGKhlZl8L86++G+3mFPLLKjCl5EhfP5+dN4M0j+eEc7AMg67WZloGDFKzyyjLTsdu19B1HZvNxuHmWqYWlGAYrSCgMP967Hb3uVympDIsEH1hhBDYMtyjZlEIeYJOUROZJUuWyK1bt471NBSK84LENtMznQEWpHmGBQJiRfyOHtrHm+++R6/oIi3SjRAw01bMjMuvoUvrp79tN3u31dPe0EI02E2xPpOiqiwuvvdePGk+dNOM12fKBj7YJwJiAiE0ja6WZlr6I1RVlGB3OpDSIhKO8NrejaxbcDGRSCM2ux2PuxghtAltRZi6iRkIY/RGsKU7cWR70ewf7X6EENuklEtO9Jzy3CgUirNCN02aIlEArs/x4bTZhjvOdXW1ET6yn417jxAOtTPFl8GKRVdydMcmjJXleHKz2Pvcz6irHoLOCCIgmTrFjZFXxMWfuRtvWipGXCSqsjNPWgbc1HVaag4xZfoshlx+ZlV6cDjtRMIRDjUdYaB3iAumFuNw2DGMmEhMdIe1ETaINg7gLE3HnuFBaGLULIkEE3vFFAoF8H7l03NBosprYpvJrmnD/SK62lvo3ruNRw/XsbiqCoejiSsvuQx7SxRNdNE3W4OD29nZ1MaRLV3MWljJocHDZGuVrP70x0nPL8SblgrErIeThbyauk5PWwvBQB+F02bGndX2YSti16Fq5lSWMZi2A5c7AyE0vN6J7bCWlkRakmhjPyLVgWbXRl0gEiihUCgmONKyCHS048vNG1WxsOLb1HWhCADL0r3D20yh/gD9gwP87s/PMXvmAgqG9rJk+izshdMJd9SQtmAxW4+sp3NHNcG6ATpcHeSEByju8WKfewtVS6vIKSvjZzsP8PkFVXhdsZanJ3JWG3qUuoN7GLKnM3taFXZ37LWGEeVgwxEQkrnTKpB0kZW9Ak2zTwqBsAIRbBluXOX+c2JFjEQJhUIxwRGadk5Eoi4UGU6aA4atiP6eTg489SN60xaRHegnRxti1k0Xk5OfTygcpDoSoHPrn6h9twbRXkO/mcLq5Qs5UJNLxW13UuHx8lJjCzdp2jEi8YE5GAadzY00D5oMOgtZXpGNw+3AMiykkBxuPkqeKxdPei8WHZNim8kyLIye0DHO6nMpEAkm9ioqFAqAc7btBO8LRGigH8Ow2LlzD3sGVrJsehqzFmXhGrAIVR9Gz5/HntrXSD3SzLvvHmRgMBXDkU52ahX+6au4fm0JXp+PXW2dlKalYNO0E9ZnkpaFtCx6WluoaepkycJ52JxOhJDoIZ2OhkGmTPUzvagMXW/F4ymb8M5qeF8kAGwZ7o/srP4oKKFQKBQfiiYEZR4XmhAYEZ1ATydbH/0lgxmVTMkvZekCwazZM+na+hzZq27BmrEUy7QY2LyPZw8O4nOEKE218ETnsfhTN1AxYzqmlNhtNubnx6opnLAU+FCQob4eNE2QWTCF5XkF2J0OQqEQ1U1HCQ1EWVRVhbBBNNSGx1M04a2IRC6E2RcGwJ7pGVORACUUCsWkINlZuZaUw5ZDAk0Ihvr7aH57A5ub6kCmkentpGrOhTS8/FvCPkHBBTcStiQHdq8nz+5h64EQmVY/qdlVXHbHrbhcKdRHjWPCXk8W0RQJhXn8xY0sXbaYqXl+NLsNaZqEQiFe3b2RxQULyS7x4XDFrBCPp3jCWxFG2IChKLYMN/bMWL7IWGw1HY8SCoVigiMtidkbTlrClSUlDWF9uPxGgvbGOva/+TT7euwIOcSAabB68aXI/hZm3fJVcHnYseVZdlZ3cbi6lduvW01mWicX3/JpCktm4UnxohsGdPUOd5w7WV5EoKMdBNxy5YU4PZ5YZnXUYF9dNUKDNbOXIUQvdmfm8HETWSQsw8IyLCL7u3DNyh4X4jASJRQKxQRHaMnNytWEOEYkjEiEro5Wnv79g8hMH7MyM7D50siwRXD1R/HNXEJrcxPh0CBPP7eHa6dlEYm0ku71cvMXvoM/Kx/Tstje0k5toJ/rp5WfUCAAgoND6IP9w855oWmYpom0NAa6wqSFM8gs7Cei78XvWzqhxQFiAiEtiV4fwJbuwjUrG7t7/L0tj78ZKRSKMybZn0ATiXOhwUGaNr/C7rYQ0m5y5drLycstonPHO7SEO4ikezhcvY9fPf4yUzM8THMZ+KYvY1VRBeWzF2IIwa62TpCS+fk5zM/POaHDOjQUc9r++pm3uff6VWj22FtTOBKmtrWBGcVT8eV6cKZ24kmZCkyd8L4Iy7AI1/Zh97lwlvrOaV7EmTKxV1qhUIwKvR2tdO7awv7OCFZqJjMLHXRnlOFsa8VRVEFkRgn9ta1Miaaza+s7XJMbIXfqbEhdyLTFC5CaQNpsHOjoYk5uFjZNO6kVMTgwyG/+/A733XAR991wEZ6U2N581Ijyds07rJy6Alv8WG/qxE6aS5DwKdkzXNgzxt5Z/WEooVAoFMOE+gP09fXwwmMPEDamsmiam/beHRSvuJPM6n2EKgro7+un9+BOGl97jW2yhhtXryJqTaFy5WpsLheWEOzr6GJObvZJ/RAAphErJHi0uYtPXbcST4qHqBEdft6mCebk2LDMNqSc+CGvlmEhNHFM8pwjyzturYiRKKFQKBQAtNTXsPvPP6ZJTCMU8bN0RgoLLr8O07LwpKbSQz/dNQf54xvPkd9xgEaRw4wpufhLK8kpKsbujlVi1eCUAqGHdWx2G/trGplVWcysymJsdhtRI8qGAxtZU3UhDrsDTbOTk71qwgsExIr46fUBtDQnmk0b1Uqvo4ESCoXiPGewv4++3h5+9+hj1IX9LM9uIndKIbPXXok7JVZ3ydCjNLz3Du9uPYK7dStdabmUm24uuvka8svKP5Dwd7KQ12AwzPpNe7j6woXDAgEQCoVw2B24wx6EfP8NdKL7IRLlN+SAPuyHgPER8nomTOzfgkKh+Eh0tjWz/ulfE7YyyTP7ybAPsHrNp3B2B0jxxkRiqKefra9tYcubGzByh9D8KVROu4Zlly6nqKzstK5j6FF6O9rY3jzEuhVzcbqdw88NDg7y+3ce4xMrb6I8vxmhLRqNWz3nSEsS7QoOl98Y736IU6GEQqE4D9GDQVrbm3n7rfUcCAyysNjN0ks+RvXDT2FvayNr8YX09w0SDg6w7ek/0tRajxnp4raVnyd33jI8Xh/CceKtpQ9cK6xzqK6FqopC1uYKHM7Y244Zb3ZU21bLumluPB4PHs912GzOU51u3DOyy5zQBDb/xBYJUEKhUExaEtVeRybNhfoDGIbBvud+zYZmG4WZJtdOm8PMlZcSrj1M7pqZNGXpbNv0EM1va+iDRzFlO2kZNmZU3EHJ8kvxpHgJRiI8W32UG6aXnzDcNUFUNzhU18KMsinYnbHXJTrOHWw4zIziqVTkp+F23zjhBQKOLeJnz/Jgz/RMuG2mE6GEQqGYZFhSDmdXA8M1mnramzj055/Q7p7H1o5MCr3vsWbNF+h/dwuBA7vpKk7jN+8+zJr2O+l/czOZVZcw4OzGYZvJNffcSWp6Np4UL4ZpUt3Tx7WVpSct4mcaJprdTstAhJkVRTicdsx4+Y0NezeybuEaSnNc2Gw2UlImfsirtGSs/EYwFrU10ZzVH4YSCoViEjGy61wiu1oaBk3N9Wx6/RlqOuexsDREuWsbl6z5NI6eAPLi+bSE22h8bwfRXVVsYTOZfif3fGwdqWk3Y1kCT2oqezu6mON2fWhDod72NqobO1i2ZD7FGV4EMStif301DW21LK5IxeGw43ROfIGAeOJc+wDRmgDepQU43PZJJRKghEKhmBQkrIjju85FBgepfudFXt7VRlpqFp/8xAqyfX4O/mmIgW27CK1eRHXNa5T5LmTX1r3ktGwkc1o5ZWU34vPn4Pa6MUzzA+JwvEiYesx6aa4+QOH0KpbFq7waUYP+7hAtwWZmFldQ4OshK3PZhI9mgvcjmozeEAyZeJbk4/CefBtuIjPxf1sKxXmOYVk0RaKUuJ2s8KVg1zSkJQl0d3J040u80Swp9jSx+pL7yS0uQ2iCkhtv5o3aZ8k9cpiG3+6jqaCX6y+ew7vr07np3u+Qm1M+LBJ748lzp+pbfWjz2+RXTKNwehU2pxNLWpimSXVzLZUFJZSlunC53ORkr5wUIhENRmPbTBLsGZ4JkV39UZj4vzGF4jwlYUU0RaIUuRxoQsRKgQ/009/QxWsbX6HL0clN197GwaMZZJQXITRBKDLEM7t+hb3Zzba9BxjIPESOD2av+TpTL1hLVk4BRjwi6VTbTCOrvM5YvirWTEjTsKRF80AzhWmFTCssJRptw+stnjSJc5ZhEdraimdJAfZJuM10IpRQKBQTEMOyhp3VJW4n9njCW33NId5+6yXCUcnKVWsxdRcFucVoTfXYhaC3qZMoFl278pG71uOdlsKSOZ/ngtXXkJKeSQp8wIo4kUgYepTBnq5jqrwm0ITGlNQCBBCNtuF2F04KK0If0tEPduMoTMezpGDSbjOdiIn/21MozjN006RFN4ad1ZoQ9LS2YiL5w1NPUpiukZfppKggn+DhaiyKaMmOoB8+zAsPvEtagYV73x8Ipupcc8dPKCqdjWNE9NKprAiI1Wg6UNtMVUXhcJXXkUhpEQk34/EUDzcTioSCuDzeUVuT0cTUY9ZV9FAPzplZODyO88KKGMnEtgMVivMMw7LY0h9kitOOXdOwdJ2W2hqeeugXDPX3kOMa5LIrbmbN2jsJ1hwibdYcHC4PBAvZ8KdtZF3oZc3HV6OV5nHzX/6Q8soFx4hEghNGNMWL+NnsNmZVFg/nRRyPENqwQCRE4t03HyYSCiZ3MUaZRMhreG8nAO45OThTnOedSAAIKeWHv2oCsWTJErl169axnoZCkTSOb0tqWBZEo4SGBmne8ioH2yPUR9q49477iQ70o9cdwTANhiry6OiuZW7mEt576AEWfuxTpKdEcOZNJ9Dfg8+ffVrXTwhEoohfoj6Trus4naeXJDfRLArLsGI9qyWINOe4bCaUbIQQ26SUS0703OS/e4VigjLSWT2y41ykv5eaDS9S3WfhmTKVRRcXsCY9lRSXl459uwiUZtHRc5j0TouGjSHCmftpSitneXYWzvRU0LQzEonjq7xKSxKOhHl191tcvnDNaYnFRBMJKxDB5o8lzZ2PFsTxKKFQKMYZidIbdaEI8H7inB4MIi2LvU88wOahQm676RIyc/PRpKRz+7ukLlpBV44b0dtFU2M6tlaTy29cS0a2n2A4TLo//bTnYBomNrtteJspYUUk8iKahxpZNb0Uh2PyvIUk+kVYgQiazzWpw13PFLUSCsU4wrCsYYEo87go87iwaxrBQC8Hn3mAYDDIZms2H7vxEvKKinE4HVg2OOrvZ6C/m41//A1bXjxCZQiuvGE5OUW52N3OMxaJ/TWNx/gkILbVVN1cS1qmi/I8Dz7f1NMOd5WWdYYrcW4xwgZ6XQBpyQlf6XU0UKuhUIwTLCmpC0WOKeYnDYO+7k4OvvJnHqvxEdV17vv4WgqKi485NhKK0NRYT6jdwbyVM1h4zcWkZ/uHn0/kRZwOx1sRieS5w821VE4pw+F0kpJSetohr9KyCHS0jzuxkJYkGozGCvk1DWAvShvXfavHEiUUCsUYkvBDjBSHMo8LaRgEero5+taLPPHUYzSklLL24iwy8/NISY/1idAjESzDJNg9SMP6IGWls7ji03exeOVF2Ef0e0jkRZypWMTmZ9HU3whCUprjGt5qOpPEOaFpw/kW4wXLsBhqDjC4oQFTN3GW+c4Lh/XZoqKeFIoxwpKS2mAYeF8gACID/dRteYU3DjcQchSxZm4VM6vKcHhd2O2OYYF4/Z0NVGaW43S58Of48eVlnvRaiXpNp4NpmkghsWt2pLQIhurxekqBMxOI8UiiPlO0J4jeMIBjRgbuNPdYT2tcoKKeFIpxRMJ6SIS9jkyca6s/yq6nH2NH0EuWp4vbr7uFnKw8+vZsI3PBUvRIhFfffJWpGaVUOlIorypHczrQ7KcWgdMViVAoxOHWoww6e1k2ZTl2zY7XM3mqvCZCXh2ZXux+Dzbn6a3L+Y4SCoXiHDLSD5EQCbumEQkG6epo5fGHf0u7KOHWK+ZTXlyMPzcfS1pkLlhKKDiAx5vGNGcaGXlT8M+Zhs2dvGY/kXCEV/e8zdo5F2BKN7a4OEx0kRgZzaRCXs8OJRQKxTnCsCzsmja8xQTQ1VqP05XC/pcfp7ZfI+KMctMls5g/dwFdOzYTSk2lXe8kw0zlqV/9B7fe/w0qLliBsCe3GJ20JOEBk3XzLsLldiFl2YQXCIhFMxlNAzjLfJOumdC5RAmFQvEhHJ8ZfaYYljXcKyJRBhzg6IEdvPn8T5AZi7Flz+a6dXMIRwbILyzFsAzasuzUNb1FjiuLgz0B+jpCBPr7SPFnxhzT5ulvKSXyIk6G0ATpWe+37ZzoImEZVqwEx4hoJsXZo4RCoTgFiZaiIzOjz+RY3TRjtZlcjuGGQkYkQm93J2++soFD1iXcOmseVXOn4XLY0HfsZyA9ncaeo2x+4iFu/fS3Sc/MpbTYRnHW15hSUoFhmuxq6wQpmV+Q+6FiMTK7+sPEYqIjLYllWOj1AWzpLhwl6coPkQSUUCgUp2Cks/lM0E2ThrBOSyTKsnQvTpsNaVn0tbfSuv0tdvaEiPgruSfXxdT509HcsQJ79elRDu9+gDLKCbbC0MAQGTmxf9MpJRVAzIqYn58z/PjDOD4vYrIiLUm0K4gQAmepT+VEJJHTFgohhA3YCjRLKa8VQmQCfwTKgDrg41LK3vhrvwXcB5jAV6WUL8XHFwO/ATzA88DXpJRSCOECHgQWA93AbVLKuvgx9wB/F5/Gv0gpf/sR7lehOGPORiQSLUkTmdWRYJBt77xFw4G9RIsWsnBmERXTi5DRITS3g4ZAPWm6g9ef/xV33f4P5OaXUVKwkKKKaSe8xskEwrIk2gneHBMiIS05Jm+e0rJGLY8icU+J7za/yqxONmeyml8DDoz4+ZvAq1LKacCr8Z8RQswCbgdmA1cC/xMXGYCfAvcD0+JfV8bH7wN6pZSVwH8B34ufKxP4LrAcWAZ8VwiRcYb3qFB8JEKDQ6f9WktKWnRj2IqwaxptTfXse/Ex3theQ3d6JdevWUrV/KkIofPkI98kNNCLIxigt6sHe4OP0EAEm81+UpE46bUtSWNvEMs6NjfKkrGMaGlJ+rtDSOvc5k6NZma2qZsY3aFjCvmdryIhTRN5BkmVZ8JpragQogi4BvjliOEbgMSn+98CN44Yf0RKGZFSHgVqgGVCiAIgXUr5roxl+T143DGJcz0GXCqEEMAVwHopZU/cWlnP++KiUIw6ocEhXnzxmVOKxfGZ1SVuJ06bjb6uTg7u3clvH3uKXRRy5+3XcN/Hr8SXkYbNbqO1pZVDjQW0t3WQnzeHylnzuPELX6Ri5uyzmqumCYozvMdYFIZl0DjQiGkZH3BYnytGKzPbCBvo9QEsMxb+ej7XaLJ0nf7qvQzUHBgVsTjdracfAN8A0kaM5UkpWwGklK1CiNz4eCGwacTrmuJj0fjj48cTxzTGz2UIIQJA1sjxExyjUIw6Treb2dOX43R/MHv3+CqvZR7XcH2m3e9t4bUdGwnpLlZWLWHZRUvwpBx7jqKiMi4uyqaoqAybLfaveCYicaLtnJEikehdbUqDcLgxljg3RmU0knldacmYJdE0gKM45qw+X30R0jSxdJ2hhhoAUsunI04zEu5M+FChEEJcC3RIKbcJIdacxjlP9BuTpxg/22NGzvF+YltalJSUnMYUFYrTw2a3MXVW6QccwYZlDfeJSORFaEIw0NvFgdef5dnGNuZnzqJ8aglzF8/G4XIQNXQc9vcT5JxeD6s//TmcXs8ZzSmxhRPoaP/AJ3VpSaRIWDcaxemx4oGCiR/ymhAIayCC2afjKE47r+szWbpO/+G9RDqayF6xDs3pHBWRgNPbeloFXC+EqAMeAdYKIX4HtMe3k4h/74i/vgkYWdqyCGiJjxedYPyYY4QQdsAH9JziXMcgpXxASrlESrkkJyfnNG5JoTh9jheJRDOhIpdjuPSGtCx62lp44bFfs77aJFfOZN01a5m/bO6wSOyuf4OooR9zrjMRCWlZw/v9wAlForsjQEOggcaBRixpoQkNLd6SdCxJ+EnO+njDQu8aIryvCy3Nhav8/C7iZwwOMnj0EEKzkb1iHTaPZ9REAk5DKKSU35JSFkkpy4g5qV+TUt4JPAPcE3/ZPcDT8cfPALcLIVxCiHJiTust8W2qASHEirj/4e7jjkmc69b4NSTwEnC5ECIj7sS+PD6mUJxzEr6IkaU3AIIDA+zauZmNr7zJkG0m937ycu79zJWkprqGRcZhdzKvdPUxFgWcfvnvhBDB+wIxUiQsaWFh0qcdoDBtCsVpxWjjxIJIbIGdjViYujnsqHZkevHMzcHutp+3vggAMxSi46WH8RaWkz59DjbPmVmkZ8NHkeR/BR4VQtwHNAAfA5BS7hNCPArsBwzgS1LKxH/DF3g/PPaF+BfAr4CHhBA1xCyJ2+Pn6hFC/DPwXvx1/ySl7PkIc1YozopEjSZ43xcRDYUIDg6y7YWnadbtXHDRCgoLS/CkerH0KD073yNzwVI0ZyxH4kQisbejizm52R+aD2FZkuY+HX+u/ECoa+KNuDCtkPLspdhOs0/EuUITGoVphWcsXEbYILyvE0dhOo5sz3ktDiOxeTzkXX33ORGIBKrMuEJxCoy4P8CuacdENnW3NNH49noadTfeinksmV5IWkY62shP+XoUzek4ZQmQU5X/Pr7sRuLnRHhrf3doOIopsc00GUiUArcCEUhxnNfO6nOJKjOuUJwhlpQYlsU7fYMIYFVG2vBWU0dzAy/86TFcMy/kghm5FBQWoknrGJEAhkXiVCVATiUSx5fd0DRtOBciPctzTKjrZBEJI2wgB3WQnNfhruMN9VtQKI4jsc3UFImy0p86LBLRUIju1jY2b9hEJK+Kqy+YTXFpCZq02Lv5OYxI5APnOtMSICP7VI8UCWlJAp0hgGGBmEyfso2wgRE2iOzvQqQ6sWeprabT5Vy0mFUWhUIR5/h2pInHZjRKf083Na/9mX39Ni64aDmXFZfh8cbyIqTdRlZZMeIktZSkYYDDccLnDD2KPe7DON6KOD7aKuHqm0wCARANRgm+14JrWhauWdnndTTTmSAtC8swGOjpHg5w+KiVjk+G+o0oxoQzac15LggbBi26AbzvrDajUdo72ujZv429LYNQuIBr15bjy/Qds81ks9mZUrhwOGluJGY0Ssf2jeQuuhDbcWJh6FE2v/4myy+5GLvTcUIrYlgUhMSZ1g0ihROnF008Er4WETLwLp2C3Z3cHhuTGVPXCXS0M9DbTdGMWQhNw7AsGoZClKV6ky4WyrZTnHMS0T6nGxo6qnOxLMKGwbNd/Uxx2odFItDZxuE3nuV///g4j27v5oI1q7jh4nlkZGd8wBcBnFAkAKJRk1eaJdHoB+81EgnS1LuZSCT4/nlGiERv58D7b6ZCw+stGfN8iGSQ6BVh9sb6hdsy3Di8DiUSp4kRDtNcfQBfbh5FM2ZhxUWiurkB/7ZNMAr/VyrqSTEmjLVFYVgWhmWxKTBEkdtJkcuB2x57s2+s3s87Lz/GYOYS5pXkU1ZRSM6UvLO+VjgYxu39YAkQgKGBAClpvuGfE1aEaRkc7XpvXIa7fhRM3STa0I+zzDfp/CyjTSLhsungPgoqpiEcDrraWtjW0UWR14548Du4p82h4o5vY/eknPH5VdSTYtwx1iKxsXcAKSUFLsdwKfCB3h4CfT08/fQbuApW8bErlpLmS/nIfRxOJhLAB0Sit3OAjJw0bJp9UolEIuRVDug4StKVo/oMsIzYlmh/Vye+3DwKp1fR0dpC08AgjW2N5G74A1GvB/u8yym5+hNnJRIfxuT4K1QoPoShgQCe1HQglhNxYUasvmWiiN/+3bvYtvF5LGcBa9YtY+r0OXi8rlOdMik0Hj5I8bSZAFiY9Ih9+FiKDfukEQnLsDD7wirk9SwwdZ2mQ/tIy8zGn1eA0DSC0SjPvPM2s7ZsIq2hhUiqxvTv/B8yS2Zgd47O3+zk+EtUKE7B0ECA557/CTPXfIbUtPRhC8KMRmlvbaJ9/yb+dLCdsrRZ3HLtRfhzMk/oh0g2jYcP8tCv/i933fdtiqfNnHRWBHBMnwi11XRmJCKaimbMRmgaUgh0Xeex//4+le9twZxaQvb9yyidfQuZRZWjOpfJ8xepUJwAS0pcKWlcddUX6bS7h3MaettbObr1Vd441IwnYyYr8uax9urlpPiSb7afDF3quK1MdPl+ocDJIhIJJ7wViKD5XMqKOE0swxiu4dXf1UlaZhYmEGhvpTEwQGp7DbN2bKYnN5WlX/wKGUWV2GzOU580CShntmLS0tbewmBqBi3xxLlEpdemwwd48dlHaAoVk+JK4947LyPdn47DdeJch9HkSPVupk6fd86vO5okIppsGTHfjLIiTg9T12k8uJf0rBwyCgqxDIPO9la2NrdQ4LYhn3wSr6OV6IL5VKz4OGnZyW3No5zZivOKSDjMkaMH+NWGDXz95ttYmZOPjISImNDX3c4rzzyNyLuIz66pwp/uJyX93FkRxzPZRAIY7janBOL0MXU9ZkFkZOHLiUXY9XR1sKX6AFm/+19CRLHlZVNw7zfxFVeinWPLUwmFYlLR19fL48//mfquOu5bdx35eVOIhoZ478H/oFrLQeIlb9GlXLxkDqlp56765skY7OkmNTNrrKeRdJRInB5GOIzQNJqrD1A4vQpd1wl0d+HLzePIwR1k/vLn4Buk5P6/J7tkDt7s3A8/6SigNg4Vk4KoYVBTV80Tm16mfqiNO6++mVlVCzEiER77w+/Z0+1hsM/OVVddyroLF44bkXj6Nz9msKd7VK8T6O0a1fMrzhxpWUSDQba//By9rS0UTq8iEo3y7PqXqOtuJ9BSj3zkIagoZvrXf0jJorVjJhKgfBSKCU5XRzMWGm/UHublHdv5m6uuJtOXSWZGNr/5xU8YPHIQ75ATvdLH7Xd8EX9u9lhP+RhG26II9Hbx0B/+m7vu+DK+jOTce8JRDcpyOBssw4h1KRSQ6svA7nYTjkTYW1PN4f/6B8rsDmb9/fcA8PqzcMbDukcb5aNQTEraWhv52WP/H5ZWxPU33MXfF5dRXFjC7h1v86N//2cc0sLqH8Q+ZyGf/ewXPlBraTww2ttOvozspIqEZVgYPbEqtkII7CNKnStOTSKzOpE4JzSNoWAQMxrluaf/RMaBPcyN9KGvuRF/UflYT/cYlFAoJhzSsth3YDft3gIqK69m3cJF5OQVYOk6P/7R/+VQ3W5yPW6cfZ1c/ZnvMnvR8rGe8piSLJEwdROrP1ZK3Z45+UqdjxbH9DoX4MvJQwpBMBjkpz/4PsvmzSXz8WfxZAmCl17Jsru+MtZT/gBKKBQTBsswaGlpounADn65v5Hb5s7lltWX4fG62fLWyzz08/+P0pxpTLFpeDIq+Pw3/xtPWtpYT3vM6W5uJKuw+KyPl5bEMiyiDf3D5TeUQJweRjjMQE93TCBy84YbYh1qaUbseocLt7+Jfng7YqaPqfd8g5yps8dl4UclFIoJgWUYPPHQ71jf3MZUX5jb5y1lzcUX0VS9jQ0bN7Nvz2vkl8ymvfsA666+j8tv+MRYT3lc0N3cyJO/+z433fl/zkosRpbfcJSkY3OOn9Lw4xnLMLAMg+ot7zB92UpsTie6rnOwpQVDj5DSu4+uF55GpkP2PTdTvug6vL7xG/2mhEIxrunr7cafkcX2vTt4tr2HaxdXse6CC0lLz+CJhx6gbudPCbvmk26mE9R1/vV/Xjppye/zkazC4rMSiZGZ1ar8xplhhMM0Hz5AWmY205etRDgc6LrOu6+9xPSKclpefZb+3fvIuudKMnIXkDd1wbi0Ikai/qMU45ZDh/bw8Iuv8IWP3cbmeoM7V87lklVrwbL43Oc+SamnmC7/x8jPinD9nTcxa86isZ7yuORMRMIyYm01rUAEW4ZbJc6dAZZhIC2LlppDFE6rwuZ0EgyHqT5yhKF9OzEf+zO1jiHs/lRKvvxV8mcvTmrinB6NYkmJ25n8kh5KKBTjDtMw2bN3Gw+/8hJRYcNms/PpdQuxC4sHfv5T9u94mezMa9gUbeGmudP51M03KysiCViGRbi2D7vPpaKZzhBT12nYvxtfTqwMuKVphHSd3z39OAuMHuyv7iCcOUjZ57+BP7uU1LwpSbu2tCwGgkO8tGs3/Zbgk8uXJF0s1H+XYlxR39RAzZFmHnv7DW67aAVV0+eRm5dLNBzii5+7lsyyi3DkpDCoP85/fPmnzKiYOtZTnvBYhjW8tWTPcGHP8KgifqfJiSKawobBW0eO4GutZeajD2LmZSAvW8qiS24hs7Aiqdcf6mmlqaWV3U1HqMxMZ/qcC5VFoZi8BIcGqD5Szfdff5O/ufhqPrX6MhYvW0BH/WHIy+NzX/wUZSUf50jna/z1PX/J1DlzSUlJHetpT3hM3USvD2BLd+HI8eLI8ipL4jQZmTjnzy/Ane5DN03Wr3+B1G0bMZvr8eR4yb7z0xQvvQSH6+QNrM6UYH8XwYEeDj7yFxzxFrPu5u+Ql5WHzT46uUJKKBRjTl9vNz/+3YMMhAfIdpZRWlFISoqH3//6JzRuf5LsBbeQU7yAI+0v8v1//Dk5eWffllQRI9G3Wg7oOEt9KuT1DJCWhREOM9QfwJcb+1sMBAI89+xTFHld5PzhERw+Df9dnyG9cAr5M06Y7HzW9DQdYP/z36bbWcjRRZ/k1tlrmJKb3Eqyx6OEQjGmHDi4j427DtM/2Metl6xl1uwleN1Ovvu338Ct2RnMuA76B1g8p4yb/u7/4BiH2dUTjZgvohekwFXuUyGvp0ki5LWvrZWju7ez4LKrsKSko6mBp19/FbH5TRw2DcpyyLnlY5SvujrpcxjsaeXoE1+n31fMqqu+zjp/Ht5R6mo3EiUUinOOqesMhobo6u7i/z35EhcVVXLvjbdSVTWXf/3Hv8fnEoQcGVQH+5jmHeBvv/HdsZ7ypEKza7grMoYfKz4cyzCo27udYF8/M1dcyFz/pUi7nQNNjURa6yl8+WkqHUG6p69k0ee/QkpOflKvHw72Y5om1T2D1K/6Amsrl+P3nbu6ZUooFOcMaVlYhsHzv/kF+9wWn7/+k/zFlRcze+ZcbMLii5+/j+KMRewZ3IyRksvX77qNpYuWjvW0JyVKIE6PxN+sZrfjzy6gsGImpmmy9fWXEdnZ+M1WIr99nHy/i678eVz47X9MashrONjPQKCT6rd/DsCCK77F7NJSXI7R72o3EiUUinNC3dEaHH19dDrg8b4I969dgt+fycL5fiKD/XztG7fiybyEbsNg2bx1XLZ2DUUfoeyEQvFRsQyD3pZm2o7WMPOCi3Cn+ziwbTPh7Dze2FXNxfV/os+uUfzlr+JJ9ZMzfV5SE+d6OhvZ+fL3sYLNzL/2e6T5cnB7z00l2eNRQqEYVSzDoLG5ge/85lly0vu5Y+1a/vLa1cyftZjm2iP86Cf/Tf/gPjz564h0bWH5qpu46xN3jvW0Fec5lmHQ39VJxpTC4b7Vz+/YQUdHN1PW/5mbInX0LVjOrFs+kfRKr0ODfQC889J/EUwvYvWlf0lOfllSr3GmKKFQjBqmrnN445vkzp/LstKDXLH6K0wtm4k+NMCGZ/7IU3/8Bb6yWzEd+QTaDvKz//mjSpz7iEhLquils0RaFkLThkUiPTsHIxrl4LbNTF+ynOCLz7Kgrg5R5KenfDEXfOHr2GzJ3QLas/1Vmg+9xEXX/S2X3Pj3aHYnDuxYholmH7ugA/VfqRgVpGWxZ+Pr/HH/bv5yzlzuuuX/4vdl0d/dyne/di2LnJUstVXwbuBZLrr4Dm77+N+M9ZQnPNKSmL1hVXbjLDB1nf7uTnw5ee+LhGWxt6WZt4424CbMrLqD6EuXMv+ez+POyEyqL8IK9rPv4HsEX/gahVf9kJRUPwDhgUH6X92MoySXlHmVYyYWSigUo0JHRzuP7NnBneuuJDc3j21vvMArr71Dzb6XKJ52F2/VPMMX//rL3LpwKa4kJiKdzwhNKJE4QyzDAKDp0D5S/ZkITRsWiYONDYReeIRLNryG8enPkvLZj7NozR3JdVb3d2JGdRof+y4N3lRsl32fpQvXAqAPhQhv2I330oU4PZ4xtShUK1TFqNHe1kpGejp73nuTB5/8NV7PAgCaG17hV795WeVEKMYUfXCQ5sOHKJ07HwDTNLE7HHS1NnOkoZ607i76n30UFs9lxWdiFm8yndXhQDvbfnEPLfnzOZCawadXfIKctGycrlheRPhAHfaKKThTzk1/d9UKVTEmuITBl750HSmuCnKyLmPg6LNUrLmKf/onJRKKsSORWb39lRdwp6QiLYtINMruA3tZNG8hzT3t8Kv/wb5wFtO/8W0yy6uSXgZ8aLCP1tpt9BIgc84NfKFiFmnCTt+Dr+FZWkHaghm4q8rG1IoYiRIKxajw6uMP86Nf/g8XLbqfPY1PgtjG3/3XL8lW5TcUY8jI+kyLL7+GSDSKCbyybRveUD1meBqO118i7QvXUDz9Gpxp/qRePxTsZ3AwwGtvPEikYhFLb/w5M6bOxQiGiR5uwrWwhJQ5Y+eLOBlKKBRJ5dFH/pfqvTXofbWsmH8nL7z1AK+++fZYT2tSoCKaPhqmrjPQ040vNw+haXR1d/P0W29y+7p1pLyzAd/uzXTmlDLt01/GnpKSVF+E1MM0dzax982fYRgGa6/4a9Iy8iCsYwTDBF/bgXftQjwu57gTCVBCoUgiX77/BnyZs8CWhkjN5GDL00okkoSKaDp7Eg7r5uoDFE6vQgrBwMAADz7yEPLgHvpmFJN65BD26y+jeOGFSd1mknqYoaE+Ot9+kHc7D3Lx5X9Lhi8bjzsNPRSi56cv4V09A+/ahefMF3E2KKFQfGSeeuJ/2bvzAIWZizgS7GN6aSZ/+dUfJN0PYZkmAJpt/H3iGg1GWhAqounMkZaFqeu01ByiaOZsimbOxrAs9tceIdhST97LT7HQDNH8xptUfutvk+6LGOzroGn9D6kJD9Jcsoprlt5BUV4xRlgnfKAOd1UZmV+4Yswjmk6HD10VIUSxEOJ1IcQBIcQ+IcTX4uOZQoj1QojD8e8ZI475lhCiRghxSAhxxYjxxUKIPfHnfiSEEPFxlxDij/HxzUKIshHH3BO/xmEhxD1JvXvFRyIyOMCXPr2WXdtexp+fR50xxDVXXMI3/vrroyISA81HGGg5MiwYk5mEBZHoXQ0okTgDTF2nt6WZ5sMHmVI5A81ux7Asdu/bS0rtFoKPPMgin2RrqY8ln/sqWRWzkyYSUg9Tf3QvTzz9fY40vkHl0ju5a8X15PvyMMI6+uFGnNOK0ew23Gmp414k4DTCY4UQBUCBlHK7ECIN2AbcCHwK6JFS/qsQ4ptAhpTyb4QQs4A/AMuAKcArwHQppSmE2AJ8DdgEPA/8SEr5ghDii8A8KeXnhRC3AzdJKW8TQmQCW4ElgIxfe7GUsvdk81XhsaOPtCyajuzj//7L58ktvZRw0EXFzGxWLFnGggULk3otyzTRbLZjxGEyWxQjrQjlkzhzEkX8EttMACZgs9nY/NoL8MxzGHIAR1oKnnlLmPvxe5KaXR0e6Kb+hR9Q1/guM2/+AcLupKhgKkZEZ3D9VhwluXhmlWN3n9uifqfDRwqPlVK2Aq3xxwNCiANAIXADsCb+st8CG4C/iY8/IqWMAEeFEDXAMiFEHZAupXw3PqkHiQnOC/Fj/iF+rseA/45bG1cA66WUPfFj1gNXEhMixRgQDPTw9svr2fLuM2QU3UBzfwOl/mI+/5nPJf1ahq4z1F5HSl4ZoY4GUgrKJ51IHC8MI/0QSiTODMsw6GtvRWgaUypnIIWgo6mRt7e8w5rFi7H99hdY04qouvNbpPhy8GbnJu3aMqqzf/871NW+TUd3Oyuv/R6l5XMIDwwS2lsLQOq6JdjHqbP6wzgjH0V8S2ghsBnIi4sIUspWIURi1QuJWQwJmuJj0fjj48cTxzTGz2UIIQJA1sjxExwzcl73A/cDlJSUnMktKU4TaVkM9nbyg3+5mw7XpaS75xMaPMoX776PxYuT28HLMk0MXSfU2YglYxbFZBWJ44VB+SHOjuG+1UB6Vg5dbS00BkPsf/0VKja8SrCyCOZWMvuevyQ9vyip1x7s62D7n3+I3v4cztX/zrWrP0uG1x8rv/HoW6R//KIJ4Yc4FactFEKIVOBx4C+klP1x98IJX3qCMXmK8bM95v0BKR8AHoDY1tPJJqY4O0xd5+77rmZOzhzc3tVk2XQqphdy953fSPq19FCIofY69J4uMqqWYHc6J51AJDiRMCiROH0SfSIAhKbhzy/AiEaJSsk7u3eQ07ifqRu2MlTkJm/6AvJnLMKVmrwy3aGhAIODfbz17L9Cex3RFX/HLYsvxQpF6H8xVp8p/eMX4U6b+L3dT0sohBAOYiLxsJTyifhwuxCiIG5NFAAd8fEmYGQjgSKgJT5edILxkcc0CSHsgA/oiY+vOe6YDad1Z4qPjB4c4u///tvsOLidSxbcxdGBBqb6nfzD1/9P0q+VsCL6Dm7DP3MxaVOmYneOv33cZKOE4cwZ6Ycomjl7WCx0XWfTxjewejvIeewxbBlunB9fwoJLP4s7yYlztQffo3br7wlLyUXXfpOoJcn2F4BuYNS2TOhtphPxoUIR9xX8CjggpfzPEU89A9wD/Gv8+9Mjxn8vhPhPYs7sacCWuDN7QAixgtjW1d3Aj48717vArcBrUkophHgJ+H8jIqouB7511nerOG1eeOpRfvvEf1I15U7mz8lnf+3D/PP3HqS0OPlbe4auM9B6BE2z4Z+5GKdn/MaTK8aWkZnVhdOrhsuCSyHYtXMb1u9/R6otilFZyvTPfBF/SUVSndWdLUcIBgfZ8NC/UPbxr7O6YjZpaRmEBwYJPLcFd8X4dVZ/FE7HolgF3AXsEULsjI/9LTGBeFQIcR/QAHwMQEq5TwjxKLAfMIAvSSkTIStfAH4DeIg5sV+Ij/8KeCju+O4Bbo+fq0cI8c/Ae/HX/VPCsa0YPb7yxU/h8nkpL72epuYnmTn3Cv7xO8/hSUlJ+rX0UIhQZ8wNlZJXdl5YEYqzQ1oW/V2dw5nVRjRKX3Mjfc1NFM2Zh/3Nl3B4PKTffSWls65JuhWx570XaXr+3zBmXMOln/4niqfOA9NCHwoxuH4HliVxzSiZdCIBqnqsYgRbt2zkv3/6LxTlLaWtv4YMbwHf/ta38WdlJfU6eiiE0+PB0HV69m3BP3PxpPZFKD4aiYZCiccAkXCY93ZsZYpmob/6Nvl3fYyaB/6Tknu+RE5lVVLLbzTV7cWm2Tj8u3s4dMG3+OSKy/F6UjDCOpGD9QA4phWh2WwTWiROFR6rhELBO68/y9a9h2lsasEmchDBFq699QZWXXxJUq9jmSZ6OETXpg3krrx0WCyUFaE4GYlucwkrQg+FCHS1c6DmCIFnn6GQKN7SmVR+5jNEhvpJyclP6vWr971N+3NfpeTW/+VQbx8Xz7sAp9AIDw1CQ9dw4txk8EWoMuOKk/Lobx7g9W0bSEvJw2maBKP13Hff/cyZPz+p1wkPDRLqim0xuctLhsVBiYTiZCS2mtKzc4CYJbrxjVcpK7TD6y+Q3d9N5r13UzT/QuxeL3avN2nXNvUQgUAXR3Y9jfPS71NcPpcpRQbClAR27iP8ZjW+z62b0BbEmaCE4jyl/uB2Xn59I7U175Ltv4zWjje4/JLL+Pjtdyf1OoauY5kmXe9uwFNRiq9oBprNpraZFCfE1HVs8Q8PQtPw5eZhmiYDba2E+/twbt5IZ8Nh0hySlI/fTdkFlyd9DqGBHmo3/Io+vYfl676CP7sITIvw3ljARdqCmXhmV0yKsNfTRQnFecjXv/oZ2m1HKXAvYWrlBdQdrecXP/1t0kM19VCI3gNbyZqznNxVlyo/hOKkWIaBtCwObdrIjBUXYnM6kZZFKBSiurmJCq+Htj/8CV/vQby3fxp/WQVZ0+YmdQ7dHQ2k+XJ4640n6cuexhWzLsTny0YfCsVCXTUbrpml2N3O88aSSKCE4jxi40vP8O3//DcuXXg7GtNxamHu/8JXk15TKFGXaaijDkdmRszJp7aYFCfB1HWaDu2nuGrOsEhYhkFncyOvvvAM0wd6aHN7mPHFz9HXXMuURauSPocjB7dw5MVvs/gTv2LNuk9gCQ230zWcXe2/fQ2eORWTwhdxNiihOA8wo1E+e/+1FPtnc93i+3mz9bdcf+ldfObOTwHJTfoydJ1QRwOe3BI0YSOlYKqyIhQnJeGHSPVnxEJeLQthGAQ624m01jLt7fVEcjLJ/dp38OYV4M0rSN619TB1h7fhz69g9ws/IH3Nd8nIKUYTAsswY1917bESHOO4V8S5QEU9TXK6W+v567+9j8opN1M3cJD21v384ddPkToK+6t6KESkuxlPbsxZnaj8qlCMJBHuOvK7tCwMy+L5N95gaeEUvMFBOje+TGCojal3fpnM0qrkzkEP8+4D9+Nt3oTvzv/Fn19BWnoOdocdyzCH+0UA540VoaKezkNq929lw7tbWf/ew+RX3ELTkSf40le+x4KFi5LuizB0HWC4/EZim0mJhOJ4pGXR196KLydvOOwVoKO5kdT0dPRHfkVg5gxC/gyK7/gspR4XrvSMDznrmTE01E9tzXbM0A6ezr+Zr5XNJ9XmJLLvKFp8e8ldVXbeCMTpoIRiElK7fyvf/tE/UOpbxozsOzl09Bl+8YPHSfX7knqdRF5Ef/VOMmcvI3P2MuWLOEPOp54TCcsBGYtoSs/OQVoWwUiYPzz8Cy6ZM5OZsg8xZzZlF12DPTW5lQBM06Cnu5VXX/8lGQSZcvmP+XLpfNJdHkJ7a7GsET1PlEgcgxKKScbXvvJZhkJDTCm6jP0HnqU8dzp/+PVzSb+OoesMtBxB7+nCN2OhEoiz4Hzpg50QiESNpoQV0d3SRFftYQ5seovLdm1kINRNxsc/RdWVH0/q9UNDAZzuFA7uf5mG3Ru4dN1XSE31YYYNQq/uxbpmOZ45FYASiJOhhGKSEBkc4L6v3Epl8S1I6z2aG/fw5MMvY7MnrwdwAss0CXU0kDZlKpwnVV5Hg/Oh/4S0LPraWocFQmgaUcPAJgS7trxGxtYN5HQFCaxax9JPfTnp20w9LYfZ/tr/sPDyv6Rhz+ssuezLZGVOITw0SN/PXwOPA8s0z7tw1zNFCcUkYN+ubTz6wHeYmreCmo4n+PUPn8FmsydVJEY6pidrI6GxYLKKRKJXRKJPBIAlJeFQiKef/D3LygvIWf8kaXfcx9TK+fgKS5PWsxpiVsRQXzuHH/0relKm4PSkcfH138Zl9w53nMv84jrsdud5H9F0OiihmMD87Bc/J9QXoq+rF5l+IYIwv/7hMzhdyf10lGhJmjZl6jFioVAcTyJxbmSvCEtKTNNk9+Fqap94mLK9O+lfPI+0T36C8ouuT6pAAAwGunj7ie8StdkZWnEXl828hLS0DPShEPrhRlwzSydNfaZzhRKKCUhPWyN//4//gDurCp9hx5+dQUFZLrd/7BNJvc7xLUkVilMRS5zbR1pmNoXTq9DsdoxolN21Rwi3N6M/9RTT+44wdPGFzL7ryzjS0pIuEoGedmoObiEqbCy/4q9IzyzAZoIR1gm+tgPv2oVqm+ksUEIxgQgNBPjmN77MkL2RvOyrMTpeJ2/xNdx//xeTeh3LNLFMk/6WaqI9vRO+JallWGij4KsZT0hpJf1N97SvbVmYus5gXy9FM2JWhBQCaVn0tjRgrH8Mc8cu3BmpuG77FPMvv3VU5trR3siOF/6Znsp1XHHD3+FPy8KI6Ay+toPUdUtiXeeUSJwVSigmCBtffJL/+t/vs6z8DrrEXCLBRv7x3/9ASnryegBDTCQCTdUACGEjo2rJhO04ZxkW0pJEG/pxlvkmrVhIaREKNeLxFJ9TsUhEM/W2NNN2tIYZy1dhczoxLIujA4MUY9C35w20uloKP/9JckuXkZo7JalzsIb6AAgJOxu2PsOaq75DekYuWtQa9kUoK+Kjo4RinBMZHGDXtnd55uXXWTn1fh7veZ3PrLmce2//RlKvk6jPZJlmrPRGftmErfIqLYllWOj1AWzpLhwl6ZNWJACE0M65SFiGQV97K0LTyJhSiD+/AOFwEIqEaQiGqX/2QcKhLsLVjRTf+yXyqhYlfX49DXuJPP+vIAR5d/6Yay65B4/TgxEZsc00ifpWjyVKKMYxP/nZj+jZv560ijUY7lLeOvQgz/76KTJ9aUm9jqHrDLXVYUkTTYuJxEQNeTV1E7MvjBACZ2nMipiskUUjOdeWRKCjHQBfTh6a3U7YMDja28f+F/7EdGEj/+2Xyf78N0i7oZLUvORaEUODfXQ2V1P/0t8x98b/xJ9VhJbixxXWj7EiVDRT8lBCMU656bZLKa+4gnTftfR1t1FamMl//MsrSa/yapkmQ+11wxYETMyIJmlJjLCB0TyAlubEnu2Z1FbEWJCozTQy5FUKQdgw+P3GDSzraKZy49NEp84g7/6/pmDBxcm9vh6mrekgW17/Be6eo8y4/l/ILJkDgGWYKqJpFFFCMc744Y/+g7c3Ps/q6XezsfohZk+/iH/8l+8m/TojrQhgQm8zSUuidw1hNA/hqsrC7rafF1bEuSRhRSSS5oSmoes69eEI3TveYtHP/gN90VzyP/9Nssvn4UxPcrmYoT7aX/0ZTXueYtkdP8OdnoPPn49lmMPCoOozjR6qeuw4YfuWN/npz/+djJwLsQGbG3fx3AO/wZOS/JwIyzSHq7xOZCvCMizMvjDSiv0Na+ku7G712SfZHF/tFSASDvPu6y+RYwbQd+/AvbyK8iW34PZlJfXaengIPTSAtetpUpfexsBgL2lZJRgRnejhJoDzuk9EMlHVY8c5T/7pN+zf/Bozcq5nIHwYd1YKz//it7i9jqReRw+F6Nq7BZumkTFryYT1QwAYYQOGotj875fAUFZE8jnekgAwolG2vPoMricexpw3k6l3f5XUKcVoWnLfTiID3Wx84j8I2UKsueZvYn4Iy87Q7sNEGzpj4a7KWX1OUEIxhoSHBvjBf/w/9nXuZ6FrJfvbn+Suu77G6rVXJP1aeihEqLMRd04WqXlTJ2zIayKiKbK/C9esbOWHGCVG+iMSImEYUaK9PQzUHML25J+wZ5RRcfvXku6shthWU99bv8bZuZO5t/8XqRn56EMh+h99i/SPX4RnVoUKeT2HqK2nMeLTd1/GlBwHpc5VHBV2+nsO8ZOf/W/Sr5PoFdGzb8twr4iJuM0EsYgmOaBjy3BjGRY258S8j/HO8VaEaZpIabH5rafwv/QKeWtuQKuaRmr2FJwpyS0FXl+7m+K8Egbf+yOpS28jqIdJzcgffl4fCqloplFCbT2NI750/y0EIwMsKryDcKSaZjnIjFkz+NRd30zqdUYKxETuFZHo12DqJtGGfhwl6bGqq0okRo2RVoRpmuysPsT03DzEc+9S8KX/g7+4HM2e3LeOxqN7MPUQjU/ejznjRsou/wu0FD/Ht6RQIjE2KIviHBEaCHD/l25h1pQbaIhsIdzZyRe++E2Wrlg9Kh3nuvduJmvOcoAJKRAQd1bH+zVYgQgizakE4hxhmiY2m43w4CDb/ufHLP7iV7AMA6/fn/RrHdq7kY5nvkjRrJsRVVdTWj4P4XQn/TqKU6MsijEm0NHCN//5bgbyV9Ntt9HX1cPvH3wh6QKRyK7WbDacmdkTNuQVRkQ0SXle9G0YT5imyfZd21g0fzGuFC+z7r0BV4o36Ul9YT2CqYc4svMJMoqvpHjd18CWgpigH2wmM8oTOIqEBgL849/9Nf/8rbvxZ91AWbCelrbn+cMvnxsVkRhoPsJAyxEAfMXTJ6RIJOozWYEINr8bR7YXoYlj1suyLAzDGMNZTl4sw2CwtQnjNz8h2NmGEBr+rJlJFYme3g7CeoTnd7+Gzelh9Y1/z5Jb/h4LN/0vbsYI60m7liI5KItilPjy5+9kiHrys69nSs4gze3v8YOfPjhq19NsNtIKpw4/nohYhkXkaABXue+EFoRlWViWRW9vL83NzcyZMwd7kvfKz2eMcJi+DRtwpqZS/rXPkJoXy75OpkgcObqPLa/+O5dd94/M6TqMgzU43WkEth3A7rBjK8xU4a7jEGVRjAK33raWYv8Kcv1X09L6NosvuYr/GAWRSJTgSDBRt5oSIa9CE9j8zg9YEBATia6uLurr68nIyFAikWQsw6D/rTfp2PYyjhnTyKtYlVSBCOsRmpoOs/nVf6PCV0pmRi7lF9+LtGzooRD6rgYc04pIWzhTCcU4RDmzk4QRifCVr92LJ8WPy1PK1sFqrpk+n699/sujss1kmSZDbXUgOKbz3EQiIRBmIIzZp+Mq951wm0mLR9/09fXh8/mUQIwSlmEQHQjgykhudnVPbwcv73uTtPpNzF/9OQpyi5CWjYGdh7Ba+nCW5eKYVog7LTWp11WcGcqZPcp85fOfAKOLvLxbaDFrkYON/PEf/oNMf3Lr3SRCXhMCMZFLgVuGRbQ7iNUfxVmajj3jg0X8TNOkt7eXjIwMAoEAfr8f2wS814mCZrcnXSS6OxrY/PIPmJOeT8lVf01qei5GRCdyqBarpY+UdQtxejzKihjnKKH4iFxx2YWsXvlJWjveor37TdZdeS23Xn970q0IPRSid/9WXDk5pOZXTFiBgFj5DTmox0uBp38g5DXhiwgEArGoJyHw+/1omtopnShEI0EMI8reDT9h+eV/gd+fg7RsBLbuR7b1kbpuicqunkAooThLHvndL3nk8Ye4ZsnneLzmMZ7/5R9J8Y5O7HeiiF+iPtNEFIhEldeR5TfsTtsHBDUajRIIBAgEApSWlqJpmhKIJDLY1kJqfvJLbiQwLAtTD3Hg9V8ze+19LLvmb3G5UrEMk9D+EVaESpybUKj/wLPgs5++nOpd1UyfeRWb2MT9t94xKiIxMi8ipaA8ZqJPUJGItA8Sru3DCkRiInFcKXDLsohGo+zfvx+fz0dZWRl2u12JRBIZbGvh0L9/l8G2lqSf25IS3TR5p2+QRlOj6pJ7cbi8OGweQntriRysxzOrAt/Vy5UvYgKinNlnwA9/8gO6GuvIswrYqHcwO9vP7Pkzufm6jyX9WpZpMtR6lJSC8gkpDglM3URogsjRPhzFsW2m460I0zTp6elBCIHP58PhSG7VXMX7jIZFYVgWtcEwdk2jyOXArmlgWgCED9ThnFasmglNAJQz+yMSGgjwjW9/k3RfOcLp5wXs3H35ldx22bqkXifhrE74HyaySEhLogd1jOpe3HNycJX7T1jp1TAMAoEAABkZGcpZfQaM7A9xuiRTJIKRME6Hk7pQhJZIlJX+VOwSMC3CB+pwV5WpZkKTBCUUH8Lfff3LdEcO4stcTX/PUX78o5+ybd9eFs+ek9TrBAcDBA7uAgmevFzSC6dNSJEY7jjXPUTkYC/uxXkfcFYn6giZpklDQwMlJSXKF3GGWIZBf1fnMX0iziVdPe08eXATq2atYXp6GmUeF5ollUBMUiaEUAghrgR+CNiAX0op/3W0r/n5u69gEJP8KVeS7i0mFGjixz/6KUDSRUIPhejZvJHM5RfidHombESTqZtY/RGQ4MxKwb7cc0zHuUQ0U319PWVlZdhsNkpLS5UVcYZIy6K/q5P07JwxEYmhwT42PfevVNr9VCy8NLbVBKBxUoEwg4PYvMo3MVEZ9z4KIYQNqAbWAU3Ae8AdUsr9J3p9MnwU1122io/NvYe9zm5qIl0sys/j7775jY90zuM5PqNaD4UmbDMhy4jtR+t1ARwl6Wh27QN+iMQWk2VZSCnJzs5WFsQZMnKrSVoWlpTnVGT7uprwZxcBMbFw2B3Y7e4PtRzM4CB9z/8I/9VfVWIxjjmVj2Ii/KcuA2qklLVSSh14BLhhtC72+c/eyqKln+SnaZ0Mmi1876tfHhWRCDQeJNBUzUDLESzTnNAiodfFfAzOMt9JndUNDQ34fD6ysrKUSJwFiWZC0oqJsiUlB5qaMEd84BhN+rqa2PX7e+nrivWpTkn1o+EgfKAOyzj1HGzeVCUSE5yJsPVUCDSO+LkJWD7yBUKI+4H7AUpKSs76QqGBAOV5cMRoZkUkyn99/8dnfa5TYZkmenfPhM6LSKDZNZxlvlO2JFVbTB+dkc2EILamVUVF52xN/dlFzP/Er4ctCssw0Q83Dkc0fRhKJCY2E+Fj3YlSnI/ZL5NSPiClXCKlXJKTk3PWF/Kk+fjq3/yKC2aV8V/f+7ezPs+HYXc6yZq7fMLmRRzP6fStViLx0TneH3Gu1zQhEgCa3Ya7qkxlVp8nTAShaAKKR/xcBCQ/YyiOJ83Hp+/+7GidfpiJ2nVOoUigoprOHyaCULwHTBNClAshnMDtwDNjPCeFQqE4bxj3PgoppSGE+DLwErHw2F9LKfeN8bQUCoXivGHcCwWAlPJ54PmxnodCoVCcj0yErSeFQqFQjCFKKBQKhUJxSpRQKBQKheKUKKFQKBQKxSkZ97WezhQhRCdQf5aHZwNdSZzOZEWt0+mh1unDUWt0epyLdSqVUp4wY3nSCcVHQQix9WRFsRTvo9bp9FDr9OGoNTo9xnqd1NaTQqFQKE6JEgqFQqFQnBIlFMfywFhPYIKg1un0UOv04ag1Oj3GdJ2Uj0KhUCgUp0RZFAqFQqE4JUooFAqFQnFKlFDEEUJcKYQ4JISoEUJ8c6znMxoIIYqFEK8LIQ4IIfYJIb4WH88UQqwXQhyOf88Yccy34mtySAhxxYjxxUKIPfHnfiSEEPFxlxDij/HxzUKIshHH3BO/xmEhxD3n8NbPGCGETQixQwjxbPxntUbHIYTwCyEeE0IcjP9NXaDW6ViEEH8Z/1/bK4T4gxDCPSHXSEp53n8RK19+BKgAnMAuYNZYz2sU7rMAWBR/nAZUA7OAfwO+GR//JvC9+ONZ8bVwAeXxNbLFn9sCXECsA+ELwFXx8S8CP4s/vh34Y/xxJlAb/54Rf5wx1mtyirX6K+D3wLPxn9UafXCNfgt8Jv7YCfjVOh2zPoXAUcAT//lR4FMTcY3GfDHHw1f8F/DSiJ+/BXxrrOd1Du77aWAdcAgoiI8VAIdOtA7EeoJcEH/NwRHjdwA/H/ma+GM7sWxSMfI18ed+Dtwx1mtwknUpAl4F1vK+UKg1OnaN0uNvguK4cbVO78+rEGiMv1nbgWeByyfiGqmtpxiJX2iCpvjYpCVuoi4ENgN5UspWgPj33PjLTrYuhfHHx48fc4yU0gACQNYpzjUe+QHwDcAaMabW6FgqgE7gf+NbdL8UQqSg1mkYKWUz8H2gAWgFAlLKl5mAa6SEIoY4wdikjRsWQqQCjwN/IaXsP9VLTzAmTzF+tseMG4QQ1wIdUsptp3vICcYm9RrFsQOLgJ9KKRcCQ8S2UU7GebdOcd/DDcS2kaYAKUKIO091yAnGxsUaKaGI0QQUj/i5CGgZo7mMKkIIBzGReFhK+UR8uF0IURB/vgDoiI+fbF2a4o+PHz/mGCGEHfABPac413hjFXC9EKIOeARYK4T4HWqNjqcJaJJSbo7//Bgx4VDr9D6XAUellJ1SyijwBLCSibhGY72PNx6+iH06qiWm/Aln9uyxntco3KcAHgR+cNz4v3Osc+3f4o9nc6xzrZb3nWvvASt437l2dXz8SxzrXHs0/jiT2J52RvzrKJA51mvyIeu1hvd9FGqNPrg+bwEz4o//Ib5Gap3eX5/lwD7AG7+33wJfmYhrNOaLOV6+gKuJRQEdAb491vMZpXu8kJj5uRvYGf+6mtie5qvA4fj3zBHHfDu+JoeIR1rEx5cAe+PP/TfvZ/m7gT8BNcQiNSpGHHNvfLwG+PRYr8dprNca3hcKtUYfXJ8FwNb439NT8TcktU7HrtE/Agfj9/cQMRGYcGukSngoFAqF4pQoH4VCoVAoTokSCoVCoVCcEiUUCoVCoTglSigUCoVCcUqUUCgUCoXilCihUCgUCsUpUUKhUCgUilPy/wPcH+4/sVRaOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "seq_list = []\n",
    "i = 0\n",
    "seq_info = sdb.seq_info.copy()\n",
    "for k in list(aln_range.keys()):\n",
    "    ctg_name, source, _ = seq_info[k]\n",
    "    seq_id = k\n",
    "    rgns = aln_range[k].copy()\n",
    "   \n",
    "    for rgn in rgns:\n",
    "        b, e, length, orientation, aln = rgn\n",
    "\n",
    "        seq =  sdb.get_sub_seq(source, ctg_name, b, e)\n",
    "       \n",
    "        if orientation == 1:\n",
    "            seq = pgrtk.rc_byte_seq(seq)\n",
    "  \n",
    "        seq_list.append((i, \"{}_{}_{}_{}\".format(ctg_name, b, e, orientation), seq))\n",
    "        \n",
    "        x, y = pgrtk.get_shmmr_dots(roi_seq, seq, 80, 56, 12, 1)\n",
    "        plt.plot(x, y, \".\", markersize=0.1)\n",
    "        i += 1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new shimmer database for the detailed study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_list = []\n",
    "i = 0\n",
    "for k in list(aln_range.keys()):\n",
    "    ctg_name, source, _ = seq_info[k]\n",
    "    seq_id = k\n",
    "    rgns = aln_range[k].copy()\n",
    "    \n",
    "    rgns = pgrtk.merge_regions(rgns,tol=1000)\n",
    "\n",
    "    for rgn in rgns:\n",
    "        b, e, length, orientation, aln = rgn\n",
    " \n",
    "        seq =  sdb.get_sub_seq(source, ctg_name, b, e)\n",
    "        if orientation == 1:\n",
    "            seq = pgrtk.rc_byte_seq(seq)\n",
    "  \n",
    "        seq_list.append((\"{}_{}_{}_{}\".format(ctg_name, b, e, orientation), seq))\n",
    "        \n",
    "        i += 1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a shimmer indext database with smaller window (denser shimmers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_sdb = pgrtk.SeqIndexDB() \n",
    "new_sdb.load_from_seq_list(seq_list, w=80, k=56, r=12, min_span=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "147"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(new_sdb.seq_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of hits is 147\n",
      "The total aligned regions is 147\n"
     ]
    }
   ],
   "source": [
    "#find hits in the pangenomic reference\n",
    "aln_range0 = pgrtk.query_sdb(new_sdb, roi_seq, merge_range_tol=300000)\n",
    "\n",
    "aln_range = aln_range0\n",
    "\n",
    "print(\"The number of hits is\", len(aln_range))\n",
    "\n",
    "count = 0\n",
    "for sid, rgns in aln_range.items():\n",
    "    count += len(rgns)\n",
    "\n",
    "print(\"The total aligned regions is\", count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's loop through the hits (contigs that has alignments) and see the number of \"copies\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of copy: 1, number of haplotypes: 147\n"
     ]
    }
   ],
   "source": [
    "n_copy = {}\n",
    "for k in list(aln_range.keys()):\n",
    "    n_copy[k] = len(aln_range[k])\n",
    "copy_count = Counter(n_copy.values())\n",
    "for nc, nh in copy_count.items():\n",
    "    print(\"number of copy: {}, number of haplotypes: {}\".format(nc, nh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'count')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZFUlEQVR4nO3df7DddX3n8eeLgIAoAkugMYkGdzO2wFR+ZBHFcVlxJYoC44qGWUuqtOmy+KvbHSdsf+l0mLLtbqtoYWVUCPUHzYpIVgeVjVK1UuNFkd9ZUqGQJZL4q1LroOB7//h+Ug43N/d7EnPuPZc8HzPfOd/zPt/P97xPkntf+f48qSokSZrOPrPdgCRp/BkWkqRehoUkqZdhIUnqZVhIknrtO9sNjMrhhx9eS5Ysme02JGlOueWWW75bVfMn15+yYbFkyRImJiZmuw1JmlOS/P1UdXdDSZJ6GRaSpF6GhSSpl2EhSeplWEiSehkWkqRehoUkqZdhIUnqZVhIkno9Za/glqS5Ysnqzwy13P2XnDHiTnbOLQtJUi/DQpLUy7CQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb0MC0lSL8NCktTLsJAk9RppWCQ5JMknktyT5O4kL0pyWJIbk9zbHg8dWP6iJJuSbExy+kD9xCS3t9cuTZJR9i1JerJRb1m8F/hsVf0y8ALgbmA1sL6qlgLr23OSHA2sAI4BlgOXJZnX1nM5sApY2qblI+5bkjRgZGGR5GDgpcCHAKrqp1X1Q+AsYE1bbA1wdps/C7imqh6tqvuATcBJSRYAB1fVzVVVwNUDYyRJM2CUWxbPA7YBVyb5ZpIPJjkIOLKqtgC0xyPa8guBBwfGb261hW1+cn0HSVYlmUgysW3btj37aSRpLzbKsNgXOAG4vKqOB35M2+W0E1Mdh6hp6jsWq66oqmVVtWz+/Pm72q8kaSdGGRabgc1V9bX2/BN04fFw27VEe9w6sPzigfGLgIdafdEUdUnSDBlZWFTVd4AHkzy/lU4D7gLWAStbbSVwfZtfB6xIsn+So+gOZG9ou6oeSXJyOwvqvIExkqQZMOqvVX0r8NEkTwO+DbyJLqDWJjkfeAA4B6Cq7kyyli5QHgMurKrH23ouAK4CDgRuaJMkaYaMNCyq6lZg2RQvnbaT5S8GLp6iPgEcu0ebkyQNzSu4JUm9DAtJUi/DQpLUy7CQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb0MC0lSL8NCktTLsJAk9TIsJEm9DAtJUi/DQpLUy7CQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb32ne0GtPdYsvozvcvcf8kZM9CJpF3lloUkqddIwyLJ/UluT3JrkolWOyzJjUnubY+HDix/UZJNSTYmOX2gfmJbz6YklybJKPuWJD3ZTGxZ/NuqOq6qlrXnq4H1VbUUWN+ek+RoYAVwDLAcuCzJvDbmcmAVsLRNy2egb0lSMxu7oc4C1rT5NcDZA/VrqurRqroP2ASclGQBcHBV3VxVBVw9MEaSNANGHRYFfD7JLUlWtdqRVbUFoD0e0eoLgQcHxm5utYVtfnJ9B0lWJZlIMrFt27Y9+DEkae826rOhTqmqh5IcAdyY5J5plp3qOERNU9+xWHUFcAXAsmXLplxGkrTrRrplUVUPtcetwHXAScDDbdcS7XFrW3wzsHhg+CLgoVZfNEVdkjRDRhYWSQ5K8szt88ArgDuAdcDKtthK4Po2vw5YkWT/JEfRHcje0HZVPZLk5HYW1HkDYyRJM2CUu6GOBK5rZ7nuC3ysqj6b5OvA2iTnAw8A5wBU1Z1J1gJ3AY8BF1bV421dFwBXAQcCN7RJkjRDRhYWVfVt4AVT1L8HnLaTMRcDF09RnwCO3dM9SpKG4xXckqRehoUkqZdhIUnqZVhIknoZFpKkXoaFJKmXYSFJ6mVYSJJ6GRaSpF6GhSSpl2EhSeplWEiSehkWkqRehoUkqZdhIUnqZVhIknoZFpKkXoaFJKmXYSFJ6mVYSJJ6GRaSpF6GhSSpl2EhSeo18rBIMi/JN5N8uj0/LMmNSe5tj4cOLHtRkk1JNiY5faB+YpLb22uXJsmo+5YkPWEmtizeDtw98Hw1sL6qlgLr23OSHA2sAI4BlgOXJZnXxlwOrAKWtmn5DPQtSWpGGhZJFgFnAB8cKJ8FrGnza4CzB+rXVNWjVXUfsAk4KckC4OCqurmqCrh6YIwkaQaMesviPcA7gZ8P1I6sqi0A7fGIVl8IPDiw3OZWW9jmJ9d3kGRVkokkE9u2bdsjH0CSNMKwSPJqYGtV3TLskClqNU19x2LVFVW1rKqWzZ8/f8i3lST12XeE6z4FODPJq4ADgIOTfAR4OMmCqtrSdjFtbctvBhYPjF8EPNTqi6aoS5JmyMi2LKrqoqpaVFVL6A5cf6Gq3gisA1a2xVYC17f5dcCKJPsnOYruQPaGtqvqkSQnt7OgzhsYI0maAaPcstiZS4C1Sc4HHgDOAaiqO5OsBe4CHgMurKrH25gLgKuAA4Eb2iRJmiEzEhZVdRNwU5v/HnDaTpa7GLh4ivoEcOzoOpQkTccruCVJvQwLSVKvocIiyfphapKkp6Zpj1kkOQB4OnB4u4fT9mseDgaePeLeJEljou8A928B76ALhlt4Iix+BPzF6NqSJE22ZPVnepe5/5IzRvLe04ZFVb0XeG+St1bV+0bSgSRp7A116mxVvS/Ji4Elg2Oq6uoR9SVJGiNDhUWSvwT+JXArsP1Cue13gJUkPcUNe1HeMuDodotwSdJeZtjrLO4AfmmUjUiSxtewWxaHA3cl2QA8ur1YVWeOpCtJ0lgZNizeNcomJEnjbdizof561I1IksbXsGdDPcIT3073NGA/4MdVdfCoGpMkjY9htyyeOfg8ydnASaNoSJI0fnbrrrNV9SngZXu2FUnSuBp2N9RrB57uQ3fdhddcSNJeYtizoV4zMP8YcD9w1h7vRpI0loY9ZvGmUTciSRpfw3750aIk1yXZmuThJNcmWTTq5iRJ42HYA9xXAuvovtdiIfC/W02StBcYNizmV9WVVfVYm64C5o+wL0nSGBk2LL6b5I1J5rXpjcD3RtmYJGl8DBsWbwZeD3wH2AK8Dpj2oHeSA5JsSPKtJHcmeXerH5bkxiT3tsdDB8ZclGRTko1JTh+on5jk9vbapUky1XtKkkZj2LD4I2BlVc2vqiPowuNdPWMeBV5WVS8AjgOWJzkZWA2sr6qlwPr2nCRHAyuAY4DlwGVJ5rV1XQ6sApa2afmQfUuS9oBhw+JXq+oH259U1feB46cbUJ1/bE/3a1PRXZ+xptXXAGe3+bOAa6rq0aq6D9gEnJRkAXBwVd3cvnzp6oExkqQZMGxY7DNpd9FhDHGNRju+cSuwFbixqr4GHFlVWwDa4xFt8YXAgwPDN7fawjY/uT7V+61KMpFkYtu2bUN+NElSn2Gv4P4fwFeTfIJu6+D1wMV9g6rqceC4JIcA1yU5dprFpzoOUdPUp3q/K4ArAJYtW+btSCRpDxn2Cu6rk0zQ3TwwwGur6q5h36SqfpjkJrpjDQ8nWVBVW9oupq1tsc3A4oFhi4CHWn3RFHVJ0gwZ+q6zVXVXVb2/qt43TFAkmd+2KEhyIPBy4B66i/tWtsVWAte3+XXAiiT7JzmK7kD2hrar6pEkJ7ezoM4bGCNJmgHD7obaHQuANe2Mpn2AtVX16SQ3A2uTnA88AJwDUFV3JlkL3EV3s8IL224sgAuAq4ADgRvaJEmaISMLi6q6jSnOmKqq7wGn7WTMxUxxLKSqJoDpjndIkkZot778SJK0dzEsJEm9DAtJUi/DQpLUy7CQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb0MC0lSL8NCktTLsJAk9TIsJEm9DAtJUi/DQpLUy7CQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb0MC0lSr5GFRZLFSb6Y5O4kdyZ5e6sfluTGJPe2x0MHxlyUZFOSjUlOH6ifmOT29tqlSTKqviVJOxrllsVjwO9U1a8AJwMXJjkaWA2sr6qlwPr2nPbaCuAYYDlwWZJ5bV2XA6uApW1aPsK+JUmTjCwsqmpLVX2jzT8C3A0sBM4C1rTF1gBnt/mzgGuq6tGqug/YBJyUZAFwcFXdXFUFXD0wRpI0A2bkmEWSJcDxwNeAI6tqC3SBAhzRFlsIPDgwbHOrLWzzk+tTvc+qJBNJJrZt27ZHP4Mk7c1GHhZJngFcC7yjqn403aJT1Gqa+o7FqiuqallVLZs/f/6uNytJmtJIwyLJfnRB8dGq+mQrP9x2LdEet7b6ZmDxwPBFwEOtvmiKuiRphozybKgAHwLurqo/G3hpHbCyza8Erh+or0iyf5Kj6A5kb2i7qh5JcnJb53kDYyRJM2DfEa77FODXgNuT3Npq/xW4BFib5HzgAeAcgKq6M8la4C66M6kurKrH27gLgKuAA4Eb2iRJmiEjC4uq+gpTH28AOG0nYy4GLp6iPgEcu+e6kyTtCq/gliT1GuVuKM1xS1Z/Zqjl7r/kjBF3Imm2uWUhSeplWEiSehkWkqRehoUkqZdhIUnqZVhIknoZFpKkXoaFJKmXYSFJ6uUV3PqFDXult6S5yy0LSVIvw0KS1MuwkCT1MiwkSb0MC0lSL8NCktTLsJAk9TIsJEm9DAtJUi/DQpLUy7CQJPUaWVgk+XCSrUnuGKgdluTGJPe2x0MHXrsoyaYkG5OcPlA/Mcnt7bVLk2RUPUuSpjbKLYurgOWTaquB9VW1FFjfnpPkaGAFcEwbc1mSeW3M5cAqYGmbJq9TkjRiIwuLqvoS8P1J5bOANW1+DXD2QP2aqnq0qu4DNgEnJVkAHFxVN1dVAVcPjJEkzZCZPmZxZFVtAWiPR7T6QuDBgeU2t9rCNj+5PqUkq5JMJJnYtm3bHm1ckvZm43KAe6rjEDVNfUpVdUVVLauqZfPnz99jzUnS3m6mw+LhtmuJ9ri11TcDiweWWwQ81OqLpqhLkmbQTIfFOmBlm18JXD9QX5Fk/yRH0R3I3tB2VT2S5OR2FtR5A2MkSTNkZF+rmuTjwKnA4Uk2A38IXAKsTXI+8ABwDkBV3ZlkLXAX8BhwYVU93lZ1Ad2ZVQcCN7RJkjSDRhYWVXXuTl46bSfLXwxcPEV9Ajh2D7YmSdpF43KAW5I0xgwLSVKvke2GknbHktWfGWq5+y85Y8SdSBrkloUkqZdhIUnqZVhIknoZFpKkXh7g3ksNeyBZ0tT2tpMx3LKQJPUyLCRJvQwLSVIvw0KS1MuwkCT1MiwkSb0MC0lSL8NCktTLi/I0Jw1zQdRT5WIozW1PlQtg3bKQJPUyLCRJvQwLSVIvw0KS1MsD3FOYywdPnyoH0zS+5vLPx7D8OdqRWxaSpF5zJiySLE+yMcmmJKtnux9J2pvMibBIMg/4C+CVwNHAuUmOnt2uJGnvMSfCAjgJ2FRV366qnwLXAGfNck+StNdIVc12D72SvA5YXlW/0Z7/GvDCqnrLpOVWAava0+cDG2e00Z07HPjubDexE/a2e+xt99jb7pnJ3p5bVfMnF+fK2VCZorZDylXVFcAVo29n1ySZqKpls93HVOxt99jb7rG33TMOvc2V3VCbgcUDzxcBD81SL5K015krYfF1YGmSo5I8DVgBrJvlniRprzEndkNV1WNJ3gJ8DpgHfLiq7pzltnbF2O0aG2Bvu8fedo+97Z5Z721OHOCWJM2uubIbSpI0iwwLSVK/qnLayQR8GNgK3DFQOwy4Ebi3PR468NpFwCa66ztOH6ifCNzeXruUJ3b/7Q/8Vat/DVgyMGZle497gZVT9LYY+CJwN3An8PZx6Q84ANgAfKv19u5x6a29Pg/4JvDpceqrLXN/W++twMQ49QccAnwCuIfu392LxqE3umuqbh2YfgS8Yxx6a6//Nt3PwR3Ax+l+Psait136ffiLDH6qT8BLgRN4clj8CbC6za8G/lubP5rul+P+wFHA3wHz2msb2g9WgBuAV7b6fwL+Z5tfAfxVmz8M+HZ7PLTNHzqptwXACW3+mcD/bT3Men9tPc9o8/u1f8Anj0NvbZn/DHyMJ8JiLPpqy90PHD6pNhb9AWuA32jzT6MLj7HobaDHecB3gOeOQ2/AQuA+4MD2fC3w6+PQ2y7/PpztX8jjPgFLeHJYbAQWtPkFwMY2fxFw0cByn2t/sQuAewbq5wIfGFymze9Ld4VmBpdpr30AOLenz+uBfzdu/QFPB74BvHAceqO7Rmc98DKeCItZ72ugfj87hsWs9wccTPdLL+PW26R+XgH8zbj0RhcWD9L9wt4X+HTrcdZ729XJYxa77siq2gLQHo9o9e3/KLbb3GoL2/zk+pPGVNVjwD8A/2KadU0pyRLgeLr/wY9Ff0nmJbmVbjfejVU1Lr29B3gn8POB2jj0tV0Bn09yS7t9zbj09zxgG3Blkm8m+WCSg8akt0Er6Hb1MA69VdX/A/478ACwBfiHqvr8OPS2qwyLPWdntySZ7lYluzPmyW+aPAO4FnhHVf1oXPqrqser6ji6/8mflOTY2e4tyauBrVV1yzS9zHhfk5xSVSfQ3WH5wiQvHZP+9qXbJXt5VR0P/Jhu98k49NYN7i7YPRP4X9P0NaO9JTmU7qanRwHPBg5K8sZx6G1XGRa77uEkCwDa49ZW39ktSTa3+cn1J41Jsi/wLOD706zrSZLsRxcUH62qT45bfwBV9UPgJmD5GPR2CnBmkvvp7lz8siQfGYO+/llVPdQetwLX0d1xeRz62wxsbluI0B3oPmFMetvulcA3qurh9nwcens5cF9VbauqnwGfBF48Jr3tmt3df7W3TOx4zOJPefKBqT9p88fw5ANT3+aJA1NfpzvAu/3A1Kta/UKefGBqbZs/jG7/8KFtug84bFJfAa4G3jOpPuv9AfOBQ9r8gcCXgVePQ28DPZ7KE8csxqIv4CDgmQPzX6UL2XHp78vA89v8u1pfY9FbW+4a4E1j9rPwQrozoZ7e1rkGeOs49LbLvwtn+5fxOE90+z63AD+jS+nz6fYFrqc7FW39pH8Yv0t39sJG2pkKrb6M7rS5vwPezxOnvB1At8m8ie5Mh+cNjHlzq28a/AEYeP0ldJuUt/HEKYOvGof+gF+lOzX1trbeP2j1We9tYJlTeSIsxqIvuuMC3+KJU45/d8z6Ow6YaH+vn6L7BTQuvT0d+B7wrIHauPT2brrTje8A/pIuCMait12ZvN2HJKmXxywkSb0MC0lSL8NCktTLsJAk9TIsJEm9DAuNlST/OOn5ryd5f5t/V5L/spNx85P8LMlvTbPu/ZJckuTeJHck2ZDkle21+5Mcvic/S1vvqUk+PfmzjOB9vjqK9bZ1n5lkuqu1tRcwLPRUcQ7wt3Q3T9uZP6K7IduxVXUs8Bq6O/aOjXYF7i6rqhfv6V4G1r2uqi4Z1fo1NxgWeqo4F/gdYFGSHW6WluTpwG8Cb62qRwGq6uGqWjvFsp9qN/K7c/vN/NqNEa9qWyS3J/ntVn9bkruS3Jbkmt1pvG0xXZHk88DVbSvp2iRfb9Mpbbn5SW5M8o0kH0jy99u3hrZvkaXzpwN9vqHVT01yU5JPJLknyUeT7HDvoKk+z6Stu1sHpp8k+TdJDkry4dbrN5OctTt/Dhpvu/W/GGmEDmx3q93uMGDddAOSLAZ+qao2JFkLvAH4s0mL/SvggZr+Zovbvbmqvp/kQODrSa6lu+3LwrZFQpJD2rKrgaOq6tGB2u44EXhJVf0kyceAP6+qryR5Dt0tqH8F+EPgC1X1x0mWA6umWM9r6a60fgFweOv/S+214+luJ/EQ8Dd098r6yqTx036e6m4OSZLX0N2996t0Vyh/oare3MZsSPJ/qurHu/UnobHkloXGzU+q6rjtE/AHQ4xZQfelMtDdH2i6XVHDeFuSb9Ht1loMLKW7R8/zkryv/aLeHjq3AR9NdyfRx36B91xXVT9p8y8H3t9Ccx1wcJJn0t3i5RqAqvos8IMp1vMS4OPV3fX3YeCvgX/dXttQVZur6ud0t4dZMsX43s+TZCndvY3eUN3N8V4BrG793kR3+4nnDP3JNSe4ZaGngnOBI5P8h/b82UmWVtW9A8tsAp6T5JlV9cjOVpTkVLpf1i+qqn9KchNwQFX9IMkLgNPpbtz2err77pxB942KZwK/n+SY6r5TYFcN/i98n/b+PxlcYKrdRlN9hGlee3Rg/nGm/vnf4fNM6uEgumD+zWp3yG3v+e+rauMQ/WmOcstCc1qS5wMHVdXCqlpSVUuAP6bb2vhnVfVPwIeAS9N97wFJFmTH7xZ4FvCDFhS/THeXT9qxgX2q6lrg94ETkuwDLK6qL9LtkjkEeMYe+FifB94y8BmPa7NfoQspkryC7kZ+k30JeEM7xjKf7hf/hmHedMjPcyVwZVV9eaD2OeCt28MsyfHDvJ/mFsNCc83vJdm8faLbqrhu0jLXMvWuqN+j+7a3u5LcQXfn1G2TlvkssG+S2+jOnvrbVl8I3NR2tVxF9/WX84CPJLmd7i67f17d93f8ot4GLGsHme8C/mOrvxt4RZJv0H13wxZg8lbSdXS7kr4FfAF4Z1V9Z8j3nfbzJHku8DrgzQMHuZfR/TntB9zW/lz/aJc/scaed52V5ogk+wOPV9VjSV5E9611x81yW9pLeMxCmjueA6xtu4t+SncqsDQj3LKQJPXymIUkqZdhIUnqZVhIknoZFpKkXoaFJKnX/wcrgahHu34EFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rgn_lengths = []\n",
    "for k in list(aln_range.keys()):\n",
    "    for sid, rgns in aln_range.items():\n",
    "        for rgn in rgns:\n",
    "            b, e = rgn[0:2]\n",
    "            #print(seq_index[k], \"target region\", b, \"-\", e, \"len:\", e-b )\n",
    "            rgn_lengths.append(e-b)\n",
    "plt.hist(rgn_lengths, bins=32);\n",
    "plt.xlabel(\"HLA Class II region size\")\n",
    "plt.ylabel(\"count\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HG02818#1#JAHEOS010000026.1_3503365_4325934_0 Memory\n",
      "HG02148#2#JAHAMF010000034.1_3417861_4208964_0 Memory\n",
      "HG02257#1#JAGYVI010000022.1_27484494_28280971_1 Memory\n",
      "chr6_GL000256v2_alt_hg38_3481069_3778827_0 Memory\n",
      "HG00621#1#JAHBCD010000020.1_23360091_23711816_1 Memory\n",
      "HG02818#2#JAHEOR010000019.1_16915483_17713098_1 Memory\n",
      "HG01928#1#JAGYVQ010000020.1_27374381_27732067_1 Memory\n",
      "HG02145#2#JAHKSF010000005.1_3502917_4318859_0 Memory\n",
      "HG01123#1#JAGYYZ010000057.1_25895025_26182990_1 Memory\n",
      "HG01361#2#JAGYYW010000059.1_32800632_33087679_0 Memory\n",
      "chr6_ssto_hap7_hg19_4032775_4439545_0 Memory\n",
      "HG00733#1#JAHEPQ010000070.1_32136286_32923993_0 Memory\n",
      "HG02486#2#JAGYVL010000026.1_32112628_32907289_0 Memory\n",
      "HG00733#2#JAHEPP010000166.1_5000_204370_1 Memory\n",
      "chr6_GL000254v2_alt_hg38_3506514_4277282_0 Memory\n",
      "HG005#2#JAHEPN010000064.1_7595432_7953375_0 Memory\n",
      "HG02630#2#JAHAOP010000058.1_23131055_23500654_1 Memory\n",
      "HG02622#2#JAHAON010000041.1_32704698_33060395_0 Memory\n",
      "NA20129#1#JAHEPE010000077.1_3426335_4225737_0 Memory\n",
      "HG00741#1#JAHALY010000025.1_23365046_24159700_1 Memory\n",
      "HG01106#1#JAHAMC010000024.1_6939124_7736742_1 Memory\n",
      "HG01243#1#JAHEOY010000117.1_3438255_4238724_0 Memory\n",
      "chr6_chm13_32018053_32812380_0 Memory\n",
      "chr6_dbb_hap3_hg19_3951412_4240123_0 Memory\n",
      "chr6_GL000250v2_alt_hg38_3497464_3812023_0 Memory\n",
      "HG02109#1#JAHEPG010000124.1_3388333_4201417_0 Memory\n",
      "HG02055#2#JAHEPJ010000058.1_23120171_23914366_1 Memory\n",
      "HG00673#2#JAHBBY010000031.1_734801_1474176_0 Memory\n",
      "HG03516#2#JAGYYS010000003.1_32273514_33096539_0 Memory\n",
      "HG00735#2#JAHBCG010000038.1_3503667_3800803_0 Memory\n",
      "chr6_hg19_32131290_32959865_0 Memory\n",
      "HG02886#2#JAHAOT010000003.1_26225747_26525711_1 Memory\n",
      "HG00741#2#JAHALX010000077.1_25651940_26443114_1 Memory\n",
      "HG02717#2#JAHAOR010000061.1_23144624_23496198_1 Memory\n",
      "HG03098#2#JAHEPL010000070.1_23116735_23415017_1 Memory\n",
      "HG03453#1#JAGYVW010000148.1_1968603_2266612_1 Memory\n",
      "NA20129#2#JAHEPD010000054.1_23125343_23929432_1 Memory\n",
      "HG00438#2#JAHBCA010000042.1_23974572_24289464_1 Memory\n",
      "chr6_mcf_hap5_hg19_3512099_4282867_0 Memory\n",
      "HG005#1#JAHEPO010000054.1_5320458_5618569_0 Memory\n",
      "HG01361#1#JAGYYX010000108.1_6944105_7739347_1 Memory\n",
      "HG03098#1#JAHEPM010000086.1_23141672_23457242_1 Memory\n",
      "HG002#1#JAHKSE010000066.1_5444315_5742327_0 Memory\n",
      "NA21309#2#JAHEPB010000021.1_3340145_4140416_0 Memory\n",
      "HG01928#2#JAGYVP010000017.1_32654512_33012769_0 Memory\n",
      "HG01071#1#JAHBCF010000017.1_299013_597352_1 Memory\n",
      "HG02723#2#JAHEOT010000107.1_23141054_23422570_1 Memory\n",
      "chr6_mann_hap4_hg19_4126658_4409419_0 Memory\n",
      "HG03453#2#JAGYVV010000230.1_3402532_4193513_0 Memory\n",
      "chr6_GL000252v2_alt_hg38_3475869_3703833_0 Memory\n",
      "HG02080#2#JAHEOV010000002.1_23360926_23649661_1 Memory\n",
      "NA19240#1#JAHEOM010000153.1_1451541_2246056_1 Memory\n",
      "HG02559#1#JAGYVK010000047.1_32167236_32465007_0 Memory\n",
      "NA18906#2#JAHEON010000020.1_23346958_24162283_1 Memory\n",
      "HG03516#1#JAGYYT010000073.1_23126159_23920647_1 Memory\n",
      "HG01071#2#JAHBCE010000076.1_6895908_7685015_1 Memory\n",
      "HG01978#2#JAGYVR010000046.1_1449927_1807865_1 Memory\n",
      "HG02630#1#JAHAOQ010000088.1_15879189_16237163_1 Memory\n",
      "chr6_hg38_32163513_32992088_0 Memory\n",
      "HG03540#1#JAGYVY010000082.1_16945983_17757447_1 Memory\n",
      "HG02572#1#JAHAOW010000052.1_189257_979210_1 Memory\n",
      "HG01175#1#JAHAMA010000069.1_6945713_7774312_1 Memory\n",
      "chr6_GL000253v2_alt_hg38_4121038_4403799_0 Memory\n",
      "NA18906#1#JAHEOO010000017.1_28384491_29178977_1 Memory\n",
      "NA19240#2#JAHEOL010000047.1_3530359_4367914_0 Memory\n",
      "HG01358#1#JAGYZB010000008.1_7551447_7849561_1 Memory\n",
      "HG00621#2#JAHBCC010000005.1_32115696_32903256_0 Memory\n",
      "HG01978#1#JAGYVS010000035.1_32160188_32981440_0 Memory\n",
      "HG002#2#JAHKSD010000045.1_26179067_26986295_1 Memory\n",
      "HG02622#1#JAHAOO010000042.1_28115560_28413635_1 Memory\n",
      "HG01361#2#JAGYYW010000059.1_32169967_32467732_0 Memory\n",
      "chr6_ssto_hap7_hg19_3480367_3778125_0 Memory\n",
      "HG01175#2#JAHALZ010000032.1_23360059_24179654_1 Memory\n",
      "HG02055#1#JAHEPK010000074.1_3492924_4219719_0 Memory\n",
      "HG01243#2#JAHEOX010000097.1_26625704_27447351_1 Memory\n",
      "chr6_apd_hap1_hg19_3971218_4220350_0 Memory\n",
      "HG005#2#JAHEPN010000064.1_7040698_7339378_0 Memory\n",
      "HG01952#2#JAHAMD010000016.1_32713709_33064809_0 Memory\n",
      "HG02622#2#JAHAON010000041.1_32133713_32447811_0 Memory\n",
      "chr6_GL000251v2_alt_hg38_3602845_4403203_0 Memory\n",
      "chr6_cox_hap2_hg19_3602951_4403309_0 Memory\n",
      "HG02723#1#JAHEOU010000100.1_5052082_5841973_0 Memory\n",
      "HG02559#2#JAGYVJ010000064.1_32122253_32942017_0 Memory\n",
      "HG01071#1#JAHBCF010000092.1_6950885_7240236_1 Memory\n",
      "chr6_dbb_hap3_hg19_3481454_3709418_0 Memory\n",
      "HG02886#1#JAHAOU010000006.1_22522980_23317712_1 Memory\n",
      "HG03486#1#JAHEOQ010000034.1_23355596_24150259_1 Memory\n",
      "HG00733#2#JAHEPP010000060.1_26893186_27481359_1 Memory\n",
      "chr6_GL000256v2_alt_hg38_4033477_4440247_0 Memory\n",
      "HG00621#1#JAHBCD010000020.1_23966447_24265104_1 Memory\n",
      "HG03579#1#JAGYVU010000035.1_16940123_17734844_1 Memory\n",
      "HG01928#1#JAGYVQ010000020.1_27986731_28284846_1 Memory\n",
      "HG02886#2#JAHAOT010000003.1_25612529_25968761_1 Memory\n",
      "HG01123#1#JAGYYZ010000057.1_26515553_26813255_1 Memory\n",
      "HG02257#2#JAGYVH010000080.1_5412300_6222285_0 Memory\n",
      "HG03098#2#JAHEPL010000070.1_22510938_22862038_1 Memory\n",
      "HG03453#1#JAGYVW010000148.1_1359413_1711637_1 Memory\n",
      "HG02109#2#JAHEPF010000055.1_1436526_2191206_1 Memory\n",
      "HG00438#2#JAHBCA010000042.1_23361552_23719816_1 Memory\n",
      "HG01109#2#JAHEOZ010000001.1_27560048_28359916_1 Memory\n",
      "HG03492#2#JAHEPH010000100.1_3433872_4223369_0 Memory\n",
      "HG02630#2#JAHAOP010000058.1_23740053_24038216_1 Memory\n",
      "HG03098#1#JAHEPM010000086.1_22525026_22888592_1 Memory\n",
      "HG01891#1#JAGYVO010000024.1_25967087_26785111_1 Memory\n",
      "HG01891#2#JAGYVN010000027.1_27120084_27919956_1 Memory\n",
      "HG01928#2#JAGYVP010000017.1_32100107_32398425_0 Memory\n",
      "HG03492#1#JAHEPI010000049.1_15861252_16689804_1 Memory\n",
      "HG03579#2#JAGYVT010000002.1_3492400_4287049_0 Memory\n",
      "chr6_mann_hap4_hg19_3475673_3773416_0 Memory\n",
      "chr6_GL000250v2_alt_hg38_4021302_4270434_0 Memory\n",
      "HG01258#1#JAGYYV010000066.1_27171028_27969215_1 Memory\n",
      "HG01123#2#JAGYYY010000050.1_32147422_32970028_0 Memory\n",
      "chr6_GL000255v2_alt_hg38_3387694_4178743_0 Memory\n",
      "HG01358#2#JAGYZA010000082.1_6957238_7752650_1 Memory\n",
      "HG00735#2#JAHBCG010000038.1_4055466_4414241_0 Memory\n",
      "HG01109#1#JAHEPA010000084.1_27049277_27790253_1 Memory\n",
      "HG00673#1#JAHBBZ010000030.1_32028390_32824277_0 Memory\n",
      "chr6_qbl_hap6_hg19_3393290_4184339_0 Memory\n",
      "HG02717#2#JAHAOR010000061.1_23753078_24051118_1 Memory\n",
      "NA21309#1#JAHEPC010000026.1_3438116_4255912_0 Memory\n",
      "HG00438#1#JAHBCB010000040.1_23356561_24160522_1 Memory\n",
      "HG02080#1#JAHEOW010000032.1_32203958_32997236_0 Memory\n",
      "HG03540#2#JAGYVX010000013.1_23368683_24125419_1 Memory\n",
      "chr6_GL000253v2_alt_hg38_3470053_3767796_0 Memory\n",
      "HG005#1#JAHEPO010000054.1_5873324_6223697_0 Memory\n",
      "HG02717#1#JAHAOS010000073.1_5416199_6229196_0 Memory\n",
      "HG01358#1#JAGYZB010000008.1_6939649_7296770_1 Memory\n",
      "HG002#1#JAHKSE010000066.1_5997057_6348138_0 Memory\n",
      "HG02572#2#JAHAOV010000201.1_5332858_6134375_0 Memory\n",
      "HG03486#2#JAHEOP010000002.1_3432636_4232522_0 Memory\n",
      "HG02622#1#JAHAOO010000042.1_27501039_27783009_1 Memory\n",
      "HG02723#2#JAHEOT010000107.1_23755078_24052596_1 Memory\n",
      "HG01952#1#JAHAME010000044.1_27450161_28271020_1 Memory\n",
      "HG02145#1#JAHKSG010000017.1_6953506_7772155_1 Memory\n",
      "chr6_GL000252v2_alt_hg38_3945827_4234538_0 Memory\n",
      "HG02080#2#JAHEOV010000002.1_23982528_24280437_1 Memory\n",
      "chr6_apd_hap1_hg19_3447380_3761939_0 Memory\n",
      "HG02559#1#JAGYVK010000047.1_32719678_33069411_0 Memory\n",
      "HG01952#2#JAHAMD010000016.1_32160890_32459032_0 Memory\n",
      "HG01258#2#JAGYYU010000011.1_26230380_27044764_1 Memory\n",
      "HG02148#1#JAHAMG010000076.1_23140769_23938206_1 Memory\n",
      "HG01978#2#JAGYVR010000046.1_2063894_2361641_1 Memory\n",
      "HG02630#1#JAHAOQ010000088.1_16490203_16788417_1 Memory\n",
      "HG02486#1#JAGYVM010000005.1_26401396_27227220_1 Memory\n",
      "HG00735#1#JAHBCH010000013.1_32215406_33034378_0 Memory\n",
      "HG01106#2#JAHAMB010000019.1_116016100_116806378_1 Memory\n",
      "chr6_KI270758v1_alt_hg38_681_58584_1 Memory\n"
     ]
    }
   ],
   "source": [
    "f0 = open(\"/results/HLA-Class2_seq.fa\", \"w\")\n",
    "for sid in aln_range:\n",
    "    rgns = aln_range[sid]\n",
    "    c, s, _ = new_sdb.seq_info[sid]\n",
    "    print(c, s)\n",
    "    count = 0\n",
    "    for rgn in rgns:\n",
    "        b, e, l, orientation, alns = rgn\n",
    "        seq = new_sdb.get_sub_seq(s, c, b, e)\n",
    "        if orientation == 1:\n",
    "            seq = pgrtk.rc_byte_seq(seq)\n",
    "        \n",
    "        print(\">{}_{}_{} {}\".format(c, b, e, s), file=f0)\n",
    "        print(pgrtk.u8_to_string(seq), file=f0)\n",
    "       \n",
    "f0.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on built-in function generate_mapg_gfa:\n",
      "\n",
      "generate_mapg_gfa(...) method of builtins.SeqIndexDB instance\n",
      "    Convert the adjecent list of the shimmer graph shimmer_pair -> GFA\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    min_count : int\n",
      "        the minimum number of times a pair of shimmers must be observed to be included in the graph\n",
      "    \n",
      "    filenpath : string\n",
      "        the path to the output file\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    \n",
      "    None\n",
      "        The data is written into the file at filepath\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(new_sdb.generate_mapg_gfa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_sdb.generate_mapg_gfa(0, \"/results/HLA-ClassII.gfa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct the genome graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "from networkx.drawing import nx_pydot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_graph(seq_list, w=128, k=56, r=4, min_span=1, highlite=0):\n",
    "    new_sdb = pgrtk.SeqIndexDB() \n",
    "    new_sdb.load_from_seq_list(seq_list, w=w, k=k, r=r, min_span=8)\n",
    "    new_sdb.generate_mapg_gfa(0, f\"/results/{name}-{w}-{k}-{r}.gfa\")\n",
    "    sid_to_path = {}\n",
    "    spl=new_sdb.get_shmmr_pair_list()\n",
    "    for (k0, k1, sid, p0, p1, o) in spl:\n",
    "        sid_to_path.setdefault(sid,[])\n",
    "        sid_to_path[sid].append((p0, p1, o, k0, k1))\n",
    "    G = nx.DiGraph()\n",
    "    m = new_sdb.get_shmmr_map()\n",
    "    for sid, path in sid_to_path.items():\n",
    "        #print(sid)\n",
    "        path.sort()\n",
    "        v = path[0]\n",
    "        #if sid != 45:\n",
    "        #    continue\n",
    "        for w in path[1:]:\n",
    "            v = tuple(v[-2:])\n",
    "            w = tuple(w[-2:])\n",
    "            color = \"yellow\"\n",
    "            penwidth =  2 + 0.1 * (len(m[v]) * len(m[w])) ** 0.5\n",
    "            G.add_edge(v, w, penwidth = penwidth, color=\"orange\" )\n",
    "            G.add_node(v, fillcolor=color, style=\"filled\")\n",
    "            #G.add_node(w, size =  len(m[v]))\n",
    "            v = w\n",
    "    for sid, path in sid_to_path.items():\n",
    "        #print(sid)\n",
    "        path.sort()\n",
    "        v = path[0]\n",
    "        if sid != highlite:\n",
    "            continue\n",
    "        for w in path[1:]:\n",
    "            v = tuple(v[-2:])\n",
    "            w = tuple(w[-2:])\n",
    "            color = \"blue\"\n",
    "            penwidth =  2 + 0.1 * (len(m[v]) * len(m[w])) ** 0.5\n",
    "            G.add_node(v, fillcolor=color, style=\"filled\")\n",
    "            G.add_edge(v, w, penwidth = penwidth, color=color )\n",
    "            v = w\n",
    "    return G\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_graph2(seq_list, w=128, k=56, r=4, min_span=1, highlite=0):\n",
    "    new_sdb = pgrtk.SeqIndexDB() \n",
    "    new_sdb.load_from_seq_list(seq_list, w=w, k=k, r=r, min_span=8)\n",
    "    new_sdb.generate_mapg_gfa(0, f\"/results/{name}-{w}-{k}-{r}.gfa\")\n",
    "    sid_to_path = {}\n",
    "    spl=new_sdb.get_shmmr_pair_list()\n",
    "    for (k0, k1, sid, p0, p1, o) in spl:\n",
    "        sid_to_path.setdefault(sid,[])\n",
    "        sid_to_path[sid].append((p0, p1, o, k0, k1))\n",
    "    G = nx.DiGraph()\n",
    "    m = new_sdb.get_shmmr_map()\n",
    "    for sid, path in sid_to_path.items():\n",
    "        #print(sid)\n",
    "        path.sort()\n",
    "        v = path[0]\n",
    "        #if sid != 45:\n",
    "        #    continue\n",
    "        for w in path[1:]:\n",
    "            v = tuple(v[-2:])\n",
    "            w = tuple(w[-2:])\n",
    "            color = \"yellow\"\n",
    "            penwidth =  1 + 0.03 * (len(m[v]) * len(m[w])) ** 0.5\n",
    "            G.add_edge(v, w, penwidth = penwidth, color=\"orange\" )\n",
    "            G.add_edge(w, v, penwidth = penwidth, color=\"orange\" )\n",
    "            G.add_node(v, fillcolor=color, style=\"filled\")\n",
    "            #G.add_node(w, size =  len(m[v]))\n",
    "            v = w\n",
    "            \n",
    "    for sid, path in sid_to_path.items():\n",
    "        #print(sid)\n",
    "        path.sort()\n",
    "        v = path[0]\n",
    "        if sid != highlite:\n",
    "            continue\n",
    "        for w in path[1:]:\n",
    "            v = tuple(v[-2:])\n",
    "            w = tuple(w[-2:])\n",
    "            color = \"blue\"\n",
    "            penwidth =  1 + 0.05 * (len(m[v]) * len(m[w])) ** 0.5\n",
    "            G.add_node(v, fillcolor=color, style=\"filled\")\n",
    "            G.add_edge(v, w, penwidth = penwidth, color=color )\n",
    "            G.add_edge(w, v, penwidth = penwidth, color=color )\n",
    "            v = w\n",
    "    return G\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ht in (0, 1, 2, 5,):\n",
    "    w=128; k=56; r=12;\n",
    "    name=\"HLA-ClassII\"\n",
    "    G = get_graph2(seq_list, w=w, k=k, r=r, highlite=ht)\n",
    "    \n",
    "    #nx.write_gexf(G, f\"/results/{name}-{w}-{k}-{r}.gexf\")\n",
    "    \n",
    "    ## very slow for graphviz to layout, remove comment to execute\n",
    "    \n",
    "    #nx_pydot.write_dot(G, f\"{name}-{w}-{k}-{r}.dot\")\n",
    "    \n",
    "    ## for this version of graphviz, we need to create PDF to get the bounding box right\n",
    "    \n",
    "    #os.system(f\"\"\"neato -Tpdf {name}-{w}-{k}-{r}.dot -x -Grankdir=LR -Gsize=64,64 \\\n",
    "    #-Nfixedsize=true -Nwidth=.5 -Nheight=.5 -Nlabel=\"\" > {name}-{w}-{k}-{r}-h{ht}.pdf\"\"\")\n",
    "    #os.system(f\"\"\"pdf2svg {name}-{w}-{k}-{r}-h{ht}.pdf {name}-{w}-{k}-{r}-h{ht}.svg\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for r in (2, 8, 12):\n",
    "    w=128; k=32;\n",
    "    name=\"HLA-ClassII\"\n",
    "    \n",
    "    G = get_graph(seq_list, w=w, k=k, r=r, highlite=83)\n",
    "    \n",
    "    ## very slow for graphviz to layout, remove comment to execute\n",
    "    \n",
    "    #nx_pydot.write_dot(G, f\"{name}-{w}-{k}-{r}.dot\")\n",
    "    \n",
    "    ## for this version of graphviz, we need to create PDF to get the bounding box right\n",
    "    #os.system(f\"\"\"dot -Tpdf {name}-{w}-{k}-{r}.dot -x -Grankdir=LR -Gsize=64,64 \\\n",
    "    \n",
    "    #-Nfixedsize=true -Nwidth=.5 -Nheight=.5 -Nlabel=\"\" > {name}-{w}-{k}-{r}.pdf\"\"\")\n",
    "    #os.system(f\"\"\"pdf2svg {name}-{w}-{k}-{r}.pdf {name}-{w}-{k}-{r}.svg\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}