{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Protein Domains\n",
    "Isotools facilitates integration of functional data on protein isoforms, in the form of protein domains. \n",
    "\n",
    "Three different types of sources for protein domains can be used:\n",
    "\n",
    "* **annotation based**: These domains are based on UniProt UCSC genome browser tracks. Annotation tables can be retrieved from [UCSC Table Browser](https://genome.ucsc.edu/cgi-bin/hgTables). Select group: = “Genes and Gene Predictions”, track: “UniProt”, and chose from the available tables (e.g. “domains”). Generally, this is the fastest approach, but it does not cover novel exonic regions nor novel splice junctions. \n",
    "* **alignment based**: Second, domain models are aligned to the protein sequences, using hmmer. This depends on model and metadata files for the domains. For Pfam domains, these files are available for download on the interpro website, at \"https://www.ebi.ac.uk/interpro/download/Pfam/\". \n",
    "* **web request**: Last, domains are obtained with the ebi interpro REST API. This approach is quite slow, and cannot be applied genome wide. \n",
    "\n",
    "In this tutorial we learn how to use the different resources to functionally annotate the identified transcript isoforms, to depict the domain structure, and to integrate the domain information with differential splicing analysis. \n",
    "\n",
    "To run this tutorial, download the transcriptome object file 'PacBio_isotools_substantial_isotools.pkl' from [here](https://oc-molgen.gnz.mpg.de/owncloud/s/gjG9EPiQwpRAyg3) to a subfolder 'demonstration_dataset'. Further, the domain resources for annotation based and alignment based domains need to be downloaded, as described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "This is isotools version 0.3.2rc6, but data has been pickled with version 0.3.2rc2, which may be incompatible\n"
     ]
    }
   ],
   "source": [
    "from isotools import Transcriptome\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "path='demonstration_dataset'\n",
    "isoseq=Transcriptome.load(f'{path}/PacBio_isotools_substantial_isotools.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add annotation based domains.\n",
    "As described above, annotation tables need to be downloaded [here.](https://genome.ucsc.edu/cgi-bin/hgTables) I downloaded the \"Domains\", as well as \"Transmembrane\", and \"Interest\" annotations as tsv for chr8 in a subfolder \"demonstration_dataset/annotation_domains\", and named them as follows:\n",
    "\n",
    "``` bash\n",
    "$ ls demonstration_dataset/annotation_domains/\n",
    "ucsc_domains.tsv  ucsc_transmembrane.tsv\n",
    "```\n",
    "The following code adds the annotation from each of the tables in the subfolter "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "domains\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2523/2523 [00:13<00:00, 190.60domains/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "transmembrane\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1103/1103 [00:02<00:00, 424.64domains/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "interest\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1988/1988 [00:08<00:00, 246.38domains/s]\n"
     ]
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "from glob import glob\n",
    "for file in glob(f'{path}/annotation_domains/ucsc_*.tsv'):\n",
    "    what=Path(file).stem.split('_',2)[1]\n",
    "    print(what)\n",
    "    isoseq.add_annotation_domains(file,category=what, name_col='name', progress_bar=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the domain structure of genes\n",
    "Isotools implements specific visualization functionality to depict the domain structure of genes of interest, to explore different isoforms in the context of protein domains. Here, we want to visualize the top expressed transcripts of the differentially spliced gene RIPK2, which has an exon skipping event between positions 89780160 and 89786592.\n",
    "While GM12878 exclusively expresses transcripts including the exon, K562 mostly skipps this exon. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGdCAYAAAACMjetAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhRUlEQVR4nO3de3xU1b3//1cukgvkIgRIyIWLUEzLpQQBwSsSIVpbFEW/SFpirVVP4GhBiqAUhSqCCuhR8bRAsD/0aAUifSAEEQ3o0XoJgkaGayPhCGiDJkEDCST798eeGWbCZCdDZjIT8n4+HnlkZu+1115r9syez6y19tohhmEYiIiIiIhHoYEugIiIiEgwU7AkIiIiYkHBkoiIiIgFBUsiIiIiFhQsiYiIiFhQsCQiIiJiQcGSiIiIiAUFSyIiIiIWwgNdgPNBXV0dhw8fJiYmhpCQkEAXR0RERJrAMAyOHz9Ot27dCA1tuP1IwZIPHD58mNTU1EAXQ0RERM7BoUOHSElJaXC9giUfiImJAcwXOzY2NsClERFpITt2wFVX+S//rVvh5z/3X/7S5lVWVpKamur8Hm+IgiUfcHS9xcbGKlgSkbajQwf/569zqrSAxobQaIC3iIiIiAUFSyIiIiIWFCyJiIiIWFCwJCIiImJBwZKIiIiIBQVLIiIiIhYULImIiIhYULAkIiIiYkHBkoiIiIgFBUsiIiIiFhQsiYiIiFhQsCQiIiJiQcGSiIiIiAUFSyIiIiIWFCyJiIiIWFCwJCIiImJBwZKIiIiIBQVLIiIiIhYULImIiIhYULAkIiIiYkHBkoiIiIgFBUsiIiIiFhQsiYiIiFhQsCQiIiJiQcGSiIiIiAUFSyIiIiIWFCyJiIiIWFCwJCIiImJBwZKIiIiIhfBAF0BEzg+lpaWUlZU1O5/q6moiIiJ8UCLf8bZMwVgHB5+WzWajGvBHTauBCJvN9/m2lWPTBiUkJJCWluaXvBUsiUizlZaWkp6eTlVVVbPzCgsLo7a21gel8h1vyxSMdXDwddnCQkOpravzWX5u+WZn+z7fNnRs2pro6GhsNptfAiYFSyLSbGVlZVRVVbFq1SrS09PPOZ8NGzYwe/bsZufjS96WKRjr4ODrsvmrrq0tX18I5rK1BjabjezsbMrKyhQsiUhwS09PJyMj45y3t9m7XZqbjy95W6ZgrIODr8vmr7q2tnx9IZjLJgqWREQkwOqPdyspKQHOBBDN4c9xLP5w5MgRtm/fHuhiBIVgOnYKlkREJGDM8W59qao6eda6bB+MWYqOjsRm2xM0X7qNGXfzOGqqawJdjKAQGRXJnt3BcewULImISMCY491OsmoV+Hqojs0G2dkn/TaOxR9qqmtgHJAQ6JIEWBmcXBs8x07BkoiIBFx6Omiojl0C0C3QhRBXmpRSRERExIKCJRERERELCpZERERELChYEhGRZgsJCeGNN94IwH4hALsFICcnhxtvvDEwO5cWpQHeIiLSoJycHF566SUAwsPD6dixIwMGDGDChAlccMEFznRHjhzhwgsvDFQxxdVx4D1gr/1xeyARuBTo5ZLuPeAdIBO4rF4enwHrXJ53ALoD1wLxLsvzgIP2x2FANJAE/Bz4aXMrEjwULImIiKWsrCzy8vKora3lm2++oaCggPvuu4+LLrrImSYxMdFv+6+pgXbt/JV3De38lXkgfA+sACKB0UAXoA7YD7wJTHFJ+xlmkOT4X18EMNkl3zeB14G76qXLAEba91MJ7AZWYwZMv2pmfYKEz4Ol+r9CUlJSGD9+PHPnziUyMtKZLiQkhPz8fGcTZkhIiHNdbGws/fr1Y968eVxzzTXOfMvLy92aeVevXk12djaPPfYY06ZNY/78+axdu5bdu3cTFRXFiBEjWLBgAX379rUs8+uvv87s2bP56quv6NOnDwsWLOD6668/p/r76s7r4j/BNCtsIPnyveqLmZYleEVERDiDoeTkZDIyMrj00ksZNWqUM43rOb2mpoapU6eyZs0avv/+e7p27co999zDzJkzAfO9N2XKFLZs2eLc/tixM/t75BGza23yZHjsMTh4EOrqYN8+uPNO+Phj6NULnnnm7LIeOgTTpsFbb4FhmMsOHz7sdguRqVOncu211/L8888TERFBSUkJhw4dYtq0abz11luEhoZyxRVX8Mwzz9CjRw8AamtrmT59OitWrCAsLIw777wTw7GDYPKm/f9dgGsM2AUY5PL8K+AUZpCzEygFPJ0WY1z+ZwAbgZOYwZjDBS7p4oBUzOkP1gE/Ay6ixTT1XOTt94BfWpYcv0JOnTpFUVERkyZNIiQkhAULFlhul5eXR1ZWFmVlZTz00EPccMMNFBcX06tXr7PSLlu2jNzcXF588UXuuOMOALZu3Upubi5Dhgzh9OnTzJo1i9GjR7Nr1y7at2/vcZ8ffPABEyZMYP78+dxwww288sor3HjjjWzfvp1+/fp5Ve9Dhw4xdOhQn9x5XfzHn3embi3MWZPT9V6Vc3bNNdeQlpZGaWnpWeueffZZ/vGPf/D3v/+dtLQ0Dh06xKFDhwCoq6tj7NixdOjQga1bt1JcXExOTg4zZ8K1157JY/9+WLMG1q6FsDAzWBo3Drp2hY8+gooKuP9+9/2eOgVjxsDw4fDee7B3L9xyC0yePJnRo0c7033yySekpaWxefNm+3anGDNmDMOHD+e9994jPDycP//5z2RlZfH555/Trl07nn76aVauXMmKFStIT0/n6aefJj8/3/mDPihUYbYgjcI9UHKIcnm8HeiP2XXWD7N1yeqU+ANgA0Jo2mjngcAm+zYtGCw1ddZ3x/dAfHx8k9L7JVhy/RWSmppKZmYmmzdvbjRYio+PJzExkcTERJYuXUpycjKbN2/m7rvvdku3cOFC5syZw6uvvspNN93kXF5QUOCWbuXKlXTp0oWioiKuvPJKj/t85plnyMrKYvr06QDMmzePzZs389xzz/Hiiy96Ve9jx4755M7r4j/+vjN1a2HOmuy796rjjunStnTr1s1jsFRaWkqfPn24/PLLCQkJoXv37s51W7Zs4YsvvqCkpITU1FRnr0JREXzyCQwZYqarqYG//Q06dzafv/UW7N4NmzZBN/uEjY8/Dtddd2a/r71mBlXLlpkDv0+dMpcfPXqUwsJCZ7qoqCiWLVvm7H5btWoVdXV1LFu2zFmevLw84uPjKSwsZPTo0SxZsoSZM2cybtw4AF588UU2bdrU3JfQt76z/29s9u+TwC7gd/bnAzDHHmVhdr05VAOP2R/bX0uG4TkQqy8U6ASUNyGtD82bN6/RniHX74GABkuuiouL+eCDD9w+LE0RFWWGwDU17vfImTFjBi+88ALr1693awL2pKKiAoCOHTs2mObDDz9k6tSpbsvGjBljeVVHdXU11dXVzueVlZVu63XX6ODX1ruN/HX3eWlbGuqGysnJ4dprr6Vv375kZWVxww03OFt2bDYbqamppKamum0TE2PensQRLHXvfiZQMreD1NQzgRKYLUiudu40W6Ri7F1CdXXm/5qaGg4cOEBsbCwAvXv3dhuntHPnTvbv30+MY0O7kydPcuDAASoqKjhy5AjDhg1zrgsPD+eSSy4Jzq64xhQDHTEHfYM5IDsO+BKzq82hHXA35likfcAXgLcNaSGNJ/Glnj17+uX71y/B0vr16+nQoQOnT5+murqa0NBQnnvuuSZvX1VVxcMPP0xYWBhXXXWVc/nGjRtZt24dW7ZsabTps66ujvvvv5/LLrvMsjvt6NGjdO3a1W1Z165dOXr0aIPbzJ8/n0cffbSJtZFg5IsbdIq0dYcPH/a4PCMjg5KSEjZu3Mjbb7/NrbfeSmZmJqtXr25y3g2MnLD0ww8weDC8/LL5/Msv4cYbIT8/n6uvvpr169cDuI2fNbf7gcGDB/OyY0MXnV0jtmDnaBdobCjiduBbwPVrzMDsinONM0IwW4cAOnNmkPe4JpSlDjjGeXPbFr8ESyNHjmTp0qX8+OOPLF68mPDwcG6++eZGt5swYQJhYWGcOHGCzp07s3z5cgYMGOBcP2DAAMrKypgzZw5Dhw6lQ4cODeaVm5tLcXEx77//vk/q5GrmzJlurVGVlZVn/UqS4NbWu0odzdAi5+qdd95xjkPyJDY2lttuu43bbruNW265haysLL777jvS09OdY5hcz5vHj8NPLS41T083B28fOQJJSeayf/7TPU1GhtkV16ULxMaCo9E/NTWVuLi4BvPOyMjgtddeo0uXLs7Wp/qSkpL46KOPnEM6Tp8+TVFRUXD1IkQDvYGP8dxddgLzarXDQA7uY5hOACuBf2MGRp5cDjyDOQVBY0HQDszuvvNk+gC/BEvt27end+/eAKxYsYKBAweyfPly7rzzTsvtFi9eTGZmJnFxcR6j+eTkZFavXs3IkSPJyspi48aNZzWbgjmYb/369Wzbto2UlBTLfSYmJvLNN9+4Lfvmm28sL4ONiIggIiKiwfXqkghevu5+au189V4tKSnxST4SnKqrqzl69Kjb1AHz589n0KBBfPbZZ2elX7RoEUlJSQwaNIjQ0FBef/11EhMTiY+PJzMzk/79+zNx4kSWLFlCcXExYLYIXXJJw2XIzISf/AQmTYInnzQDoYceck8zcaK5buxYmDsXysvN5QsXLuSpp55qMO+JEyfy5JNPMnbsWObOnUtKSgoHDx5k7dq1/PGPfyQlJYX77ruPJ554gj59+nDxxRezaNEiyh07CCbXY04d8FfMK926YrbyHAA+xRxsnQz08LBtMmbr0mgP68DsqksH3gUmuiw/hTmfk+vUAR8ClwA9m1MZ75WUlLB9+3bLNOdy3vP7mKXQ0FBmzZrF1KlTuf32251jkTxJTEx0BlkN6d69O1u3bnUGTAUFBc6AyTAMpkyZQn5+PoWFhfTs2fhRGj58OFu2bOF+l8sqNm/ezPD6neFN0KlTJ6Kjo/WLPchFR0eTkNDYCMjzW0JCgt6r0mQFBQUkJSURHh7OhRdeyMCBA3n22WcJDw/nN7/5zVnpY2JiWLhwIfv27SMsLIwhQ4awYcMGQkPNy6jWrVvHlClT3C68mT/fugyhoZCfb04dMHQo9OgBzz4LWVln0kRHw7ZtMGOGeeWco2WppqamwRYjc7totm3bxowZMxg3bhzHjx8nOTmZUaNGObebNm0aR44cYdKkSYSGhvLb3/6Wm266yTk2Nmh0xBxntA3zarQfMFucumEO4F6D2ULkSTpmkGM1HPhSYDnwf4CjLWK7/S8Ms7WqGzDenl8Lmz17dpMuNvH2e6BFJqUcP34806dP5/nnn+eBBx5odn6pqakUFhYycuRIxowZQ0FBAbGxseTm5vLKK6+wbt06YmJinOOO4uLiGgzS7rvvPq666iqefvppfvGLX/Dqq6/y6aef8pe//OWcymWz2TTPUpDTPEuQlpbm0/equvXOXytXrmTlypUe17mO8XEd6HzXXXdx1131Zy48Iy0tjXXrzOmht2/fzuDBg+nU6cz6Rx4x/+r7yU/MKQFc1R9fnZgI9qn+2L7dbLF6+OGH3YKlRYsWnZV3YmKic45AT8LDw1myZAlLlixpME3QiAF+Yf+rb4bFdpdzJpAahPu8TA6pwCMuz+84h/L5UVOHWDi+B+pfoNWQFgmWwsPDmTx5MgsXLuTee+9tcM4jb6SkpLgFTJs2bWLp0qUAXH311W5p8/LyyMnJ8ZjPiBEjeOWVV3j44YeZNWsWffr04Y033vB6jiWHtLS0Nv9FLK2D3qsicr7x1xALnwdLDf0CefDBB3nwwQedz+tfbtnY5Zee8k1OTmbv3r1NzqMh48ePZ/z48ee0rYiIiJzfmjIPp4iIiEibpWBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEQovMsyQiImLFH3eJarV3ntK8xkH3GihYEhGRgDFvvRNJdvZJv+QfHR3Zqm5v1C6iHTVrawJdjKAQGRU8x07BkoiIBIx56509brfe2bBhA7Nnz27yrSustLbbG61ds5akpKRAFyMoBNOxU7AkIj5zLnfzdlVSUuKTfHzJ2zIFYx0cfF22YK6rQ1lZmTMQC+byOspWXl6uYMnO9dg1xt/HNMQ413uEiFNlZSVxcXFUVFRY3tla5HxVWlpKeno6VVVVzc4rLCyM2tpaH5TKd7wtUzDWwcHXZQsLDaW2rs5n+fk93zZ0bNqa6OhobDabV61RTf3+VsuSiDSb2ZVia/KvQCvV1dVERET4oFS+422ZgrEODj4tm81GdXY2/qhpdV0dEatWQTO74c7Kt60cmzbIn912CpZExCfS0tKCZnyBnCfS08EPd5AX8ZbmWRIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbEQHugCyPmltLSUsrKyQBfDUkJCAmlpaS2+32B7baqrq4mIiAh0Mc4SjOXytkzBWAcHn5bNZqMa8EdNq4EIm833+baVY9MG+fPcrmBJfKa0tJT09HSqqqoCXRRL0dHR2Gy2Fg2YgvG1CQsLo7a2NtDFOEswlsvbMgVjHRx8Xbaw0FBq6+p8lp9bvtnZvs+3DR2btsaf53YFS+IzZWVlVFVVsWrVKtLT0wNdHI9sNhvZ2dmUlZW1aLAUbK/Nhg0bmD17dtCUxyEYy+VtmYKxDg6+Lpu/6tra8vWFYC5ba+Dvc7uCJfG59PR0MjIyAl2MoBQsr43N3r0RLOVxCMZyeVumYKyDg6/L5q+6trZ8fSGYyyYKliTAWnocj80PYyBa0/4dSkpKgOApj0MwlsvbMgVjHRwaK1uwjOfz5WsYqDqdqyNHjrB9+/ZAFyMoBNOxU7AkAVNaWkrfi/ty8sTJlt1xiHlCaklHjhyBEMj2wxiM5gi28jgEY7m8LVMw1sGhobJFRkWyZ/eeAIzn60tV1dnnAV+8htHRkdhsLVun5hh38zhqqmsCXYygEIj3Y0MULEnAlJWVmYHSOCChpXYKrIXy8vIW2qGpvLwcDFq2riLeKIOTa08GaDzfSVatAl8P1bHZIDu75evUHDXVNTpPQMDejw1RsCSBlwB0C3QhWkhbqquIF9LTQUN17HSeCDqalFJERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSUTOlg/8j5/zeBdY2sx9SNAICQnhjTfeCMB+IQC7BSAnJ4cbb7wxMDuXFqWr4USCQT6w0/44FIgDBgJXAGHNzPckMMHL7a7DnOrAn0YAw/y8D2m2nJwcXnrpJQDCw8Pp2LEjAwYMYMKECVxwwQXOdEeOHOHCCy8MVDHF1XHgPWCv/XF7IBG4FOjlku494B0gE7isXh6fAetcnncAugPXAvEuy/OAg/bHYUA0kAT8HPhpcysSPHweLNX/YKWkpDB+/Hjmzp1LZGSkM11ISAj5+fnOqDwkJMS5LjY2ln79+jFv3jyuueYaZ77l5eVuv1xWr15NdnY2jz32GNOmTWP+/PmsXbuW3bt3ExUVxYgRI1iwYAF9+/ZtsLx//etf+dvf/kZxcTEAgwcP5vHHH2fo0KHnVP9gu7N8SwrGGYtbld7AWKAW2Ae8iXnyucJD2tP496dOZONJmk03V281srKyyMvLo7a2lm+++YaCggLuu+8+LrroImeaxMREv+2/pgbatfNX3jW081fmgfA9sALzMzwa6ALUAfsxzylTXNJ+hhkkOf7XFwFMdsn3TeB14K566TKAkfb9VAK7gdWYAdOvmlkfLzX1e8jb2cH9crp1fLBOnTpFUVERkyZNIiQkhAULFlhul5eXR1ZWFmVlZTz00EPccMMNFBcX06tXr7PSLlu2jNzcXF588UXuuOMOALZu3Upubi5Dhgzh9OnTzJo1i9GjR7Nr1y7at2/vcZ+FhYVMmDCBESNGEBkZyYIFCxg9ejRffvklycnJXtX70KFDDB06NKjuLC+tSBgQY388BLABezCDJUcLUTLwMeYn937gG2Aj8H/ABUA6MAbzJPcuZ1qrHrH/nwT0BCqATcABIATzF2MW4GgYqN8ilQd0te93u72sl2CeIJvqa+BlzBaly+3l2w3cW2+facCHmEFjP3u5HK1rO4F/Asfs9e1pX9/Bvv4EsMFerxogFvP1G2Rf31i9xaOIiAhnMJScnExGRgaXXnopo0aNcqZx/QFcU1PD1KlTWbNmDd9//z1du3blnnvuYebMmYD5o3LKlCls2bLFuf2xY2f298gjZtfa5Mnw2GNw8CDU1cG+fXDnnfDxx9CrFzzzzNllPXQIpk2Dt94Cw946evjwYbf7rU2dOpVrr72W559/noiICEpKSjh06BDTpk3jrbfeIjQ0lCuuuIJnnnmGHj16AFBbW8v06dNZsWIFYWFh3HnnnRiGv5tfz8Gb9v93Aa4xYBfOfA4AvgJOYX6GdwKlmJ+9+mJc/mdgnm9O4v6D6gKXdHFAKuZcUeuAnwEX0WKaOut7dHQ0NpuN+Pj4JqX3S7Dk+sFKTU0lMzOTzZs3NxosxcfHk5iYSGJiIkuXLiU5OZnNmzdz9913u6VbuHAhc+bM4dVXX+Wmm25yLi8oKHBLt3LlSrp06UJRURFXXnmlx32+/PLLbs+XLVvGmjVr2LJlC7/5zW+aXGeAY8eOBdWd5Vua467Z4iMXYH75O5RgBkGOt2UN8P9hnpjuAn4E/oEZLNyEGZSUAdWYLVYAUZhBiGO732J2+20DVmEGLg2dFXYCw+37OgS8gXlybcqJ8F/Aa5hN+JdYpPsK86Q7CfgO89dpIjDYvr4WuAboZK/vJns5HOfHd4B/AxMxuwO+w2yFc2x7LvUWj6655hrS0tIoLS09a92zzz7LP/7xD/7+97+TlpbGoUOHOHToEAB1dXWMHTuWDh06sHXrVoqLi8nJyWHmTLj22jN57N8Pa9bA2rUQFmYGS+PGQdeu8NFHUFEB99/vvt9Tp2DMGBg+HN57D/buhVtugcmTJzN69Ghnuk8++YS0tDQ2b95s3+4UY8aMYfjw4bz33nuEh4fz5z//maysLD7//HPatWvH008/zcqVK1mxYgXp6ek8/fTT5OfnO3s/gkIVZgvSKNwDJYcol8fbgf6YP0T6YbYuWTW0/ID5Ay6Epo12Hoj5+bTRosHSvHnzuP766y3T2Gw2srOzKSsrC2yw5Kq4uJgPPviA7t27e7VdVJR5VGtq3O+RM2PGDF544QXWr1/v9qvGk4qKCgA6duzY5P1WVVVx6tQpy22qq6uprq52Pq+srHRb31bvGt2auuFKSkpa9GaVjhuDNomBGVzsx31MzwWYTdqOT20RZiBwE2dOjNdjDqq+FrO1JdyeJsYln532ffwK88QHZjD1BGaw0ruBcnUFrrY/7oTZwvUvGj8R2jBbjX6FeVK2EmmvQyjQGehj34cjWHL9WHXEHFv1V8yAMAKz5SgRswUO3FuMijm3ekuDunXr5jFYKi0tpU+fPlx++eWEhIS4nf+3bNnCF198QUlJCampqc4hGEVF8MknMGSIma6mBv72N+jc2Xz+1luwezds2gTd7LNbP/44XHfdmf2+9poZVC1bZg78PnXKXH706FEKCwud6aKioli2bJmz+23VqlXU1dWxbNkyZ3ny8vKIj4+nsLCQ0aNHs2TJEmbOnMm4ceMAePHFF9m0aVNzX0Lf+s7+v7FbpZwEdgG/sz8fgNl6nIV793g18Jj9sf21ZBieA7H6QjHPE+VNSOtDPXv29Mv3r1+CpfXr19OhQwdOnz5NdXU1oaGhPPfcc03evqqqiocffpiwsDCuuuoq5/KNGzeybt06tmzZ0mg0X1dXx/33389ll11Gv36NnaHPmDFjBt26dSMzM7PBNPPnz+fRRx9tcp4SfGbPnh18rWB7MU9MdZhf6v05E5zAmW4wh39jBgauJ640+7ZlnOmaqu8bzJPq4/WWn8Ycl9CQrvWex2C27lj5GrNet2J2ETamC+6/WmMwy+twGCgEjmKe8B29IBX2bYdgtmAdwQziLubMr+Vzrbc0qKFuqJycHK699lr69u1LVlYWN9xwg7Nlx2azkZqaSmpqqts2MTHmvdwcwVL37mcCJXM7SE09EyiB2YLkaudOs0Uqxv7joK7O/F9TU8OBAweIjY0FoHfv3m7jlHbu3Mn+/fuJiYlxy+/kyZMcOHCAiooKjhw5wrBhZ369hIeHc8kllwRnV1xjijF/bDiGmSVhdp99ifsPknbA3ZjnpH3AF5gtu94IaTxJa+CXYGnkyJEsXbqUH3/8kcWLFxMeHs7NN9/c6HYTJkwgLCyMEydO0LlzZ5YvX86AAQOc6wcMGEBZWRlz5sxh6NChdOjQ0LcB5ObmUlxczPvvv9/kcj/xxBO8+uqrFBYWug1Gr2/mzJlMnTrV+byysvKsD74Et6Y01fpSk7ooewK/4MzYpfpXwV1w1hbnpgbzvlPjPKzzPLTP5KnpvbHviQsxm/4/A35C41f2We3D0e3YG7gZs5utArMbrdaepg/wB8wT+wHgb5gB1BjOvd7SoMOHD3tcnpGRQUlJCRs3buTtt9/m1ltvJTMzk9WrVzc57waGmVr64QcYPBgcoyu+/BJuvBHy8/O5+uqrWb9+PcBZ5/cffviBwYMHnzUsA6Cza8QW7BwdIo1dY7Qd+BZw/c1vYH5OXYOlEMzWITBbeh2DvD19huqrwxxbeJ7c484vwVL79u3p3dts016xYgUDBw5k+fLl3HnnnZbbLV68mMzMTOLi4jy+QZOTk1m9ejUjR44kKyuLjRs3nvVLAMz+6fXr17Nt2zZSUlKaVOannnqKJ554grffftstQPMkIiKCiAhdytOa+auptiFN6qK8gDMnpqboDOzADAIcP5JLMU9wjmb4MM4OaJIwf1m2x/9XvUUDtwErMa+iGc+5T4VQhjmGKxPzVzCYLU31tce8CufnwKfAW5jBUkvWuw145513nOOQPImNjeW2227jtttu45ZbbiErK4vvvvuO9PR05xgm1x+Zx4/DTy0uNU9PNwdvHzkCSUnmsn/+0z1NRobZFdelC8TGgmOERGpqKnFxcTQkIyOD1157jS5dujhbn+pLSkrio48+co5/PX36NEVFRcE15CIa88fEx3juLjuBebXaYSAH9zFMJzA/p//GPLd4cjnwDOYUBI0FQTswW3/Pk+kD/D5mKTQ0lFmzZjF16lRuv/1251gkTxITE51BVkO6d+/O1q1bnQFTQUGBM2AyDIMpU6aQn59PYWEhPXv2bFIZFy5cyGOPPcamTZu45BKr0adN05rG7viSV+NypPn6Y15Rlo/ZXVeFeaXKAM50wcVjtrCUYZ4YI+3b/S/wKuaVMLGY4wpsmJcPN/ydcm46YA7YfglzwPYtnFvAFGff7iPMQeLfYg7QdvUO5km8M2Zr017OnPhbut7nkerqao4ePeo2dcD8+fMZNGgQn3322VnpFy1aRFJSEoMGDSI0NJTXX3+dxMRE4uPjyczMpH///kycOJElS5a4TNsCVqffzEz4yU9g0iR48kkzEHroIfc0Eyea68aOhblzobzcXL5w4UKeeuqpBvOeOHEiTz75JGPHjmXu3LmkpKRw8OBB1q5dyx//+EdSUlK47777eOKJJ+jTpw8XX3wxixYtotyxg2ByPebUAX/FfJ93xWzlOYD54+EizDF9PTxsm4zZujTawzowPyPpmOediS7LT2HO5+Q6dcCHmJ/Tpn0N+0xTxqOey3d0i1z/MX78eKZPn87zzz/PAw880Oz8UlNTKSwsZOTIkYwZM4aCggJiY2PJzc3llVdeYd26dcTExHD06FEA4uLiGgzSFixYwJ/+9CdeeeUVevTo4dymQ4cOlt18nnTq1Ino6OgmX7oo0iztgF9jBkh/xX3qAIfBmIOX/4LZAuWYOuAO4G3M8T3VmIFDT/w395HjCreVwFrMbjRvtQduBLZgBkxJmCd111nCwzDrVY75eqRhBmdgvl4tXe/zREFBAUlJSYSHh3PhhRcycOBAnn32WcLDwz1eNRwTE8PChQvZt28fYWFhDBkyhA0bNhAaavazrlu3jilTprhdpTx/vnUZQkMhP9+cOmDoUOjRA559FrKyzqSJjoZt22DGDPPKOUfLUk1NTYMtRuZ20Wzbto0ZM2Ywbtw4jh8/TnJyMqNGjXJuN23aNI4cOcKkSZMIDQ3lt7/9LTfddJPzQqKg0RFznNE2zKvRfsBsceqGOYB7DWYLkSfpmEGO1bVTlwLLMacrcXTcbOfMlCJR9n2Np2njFH2sqeNRo6OjSUhobCT8GSGGj0eneZo8EszxQIsWLaKkpIT27dt7nJTS9XlT8v36668ZOXIknTp1YtOmTQ02s+bl5ZGTk+NxXY8ePTh48OBZy+fMmcMjjzzSSG1NlZWVxMXFUVFRQXl5eZuelDI7O7vJTdPbt29n8ODB8Htarl/7MPAX8+qXiRMnNprcV15++WUziG7Juop4w/7Z8KZryfG+bk53lOM8UFRkdqP50vbtZouVo3y+KK+/OM8RoPMEeP1+dLyPmjp1j2NSStfvb6uA2uctSytXrvS4/MEHH+TBBx90Pq8fozUWs3nKNzk5mb179zY5D0+++uorr7exkpaW5tWsoCIiIuIb/pq6RzfSFREREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMRCi0xKKWKpJaelCvQUWIHev0hDAvze9MeND1rtzRR0ngi610DBkgRMQkICkVGRnFx7smV3HALx8fEtusv4+Hjznm1rW3S3Il6JjIr0alZjX0hISCA6OpLsbP+cB6KjW75OzdEuoh01a2sCXYygEIj3Y0MULEnApKWlsWf3nhad8dwxy3iS406cLSQpKQkMmjy7rL9t2LCB2bNnB015HIKxXN6WKRjr4NBY2RyzGrektLQ0bDb384AvX8NA1Kk51q5Z2+Lnp2AVTMdOwZL4XFu9kbBIa1dWVubVjxfHzbOD+TPvWqdgLq+jbOXl5QqW7Lx5P/r7mPr83nBtUVPvLXO+Ky0tJT09naqqqkAXxVJ0dDQ2m61Ff7EE42sTFhZGbW1toItxlmAsl7dlCsY6OPi6bGGhodTW1fksP7/n24aOTVtzLuf2gN0bTtouszndFvQ3Eg5cV0NwvTbV1dVEREQEuhhnCcZyeVumYKyDg0/LZrNRnZ2NP2paXVdHxKpV4OOuzDZzbNogf57bFSyJT+lGwg3TayPipfR08MNNUUW8pXmWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbEQHugCyPmltLSUsrKyQBfDUkJCAmlpaS2+32B7baqrq4mIiAh0Mc4SjOXytkzBWAcHn5bNZiMB8MenqRQos9l8nm+bOTZtkD/P7QqWxGdKS0tJT0+nqqoq0EWxFB0djc1ma9GAKRhfm7CwMGprawNdjLMEY7m8LVMw1sHB12WLDg3FVlfn04CpFEgPDaUqO9uHuZra0rFpa/x5blewJD5TVlZGVVUVq1atIj09PdDF8chms5GdnU1ZWVmLBkvB9tps2LCB2bNnB015HIKxXN6WKRjr4ODrsjk/T/i2dakMqKqr8/lr2JaOTVvj73O7giXxufT0dDIyMpqUNti6pkTaskB1UTfEm3NJU9js3Xq+ztcXHGWLj48PbEGCSDC9HxUsScCUlpbS9+K+nDxxsmV3HAJHjhxp0V0eOXIEQiDbD90KzRFs5XEIxnJ5W6ZgrINDQ2WLjIpkz+49QfMF1RaNu3kcNdU1gS5GUAim96OCJQmYsrIyM1AaByS01E6BtVBeXt5COzSVl5eDQcvWVcQbZXBy7ckW76IWdzXVNTpPQNC9HxUsSeAlAN0CXYgW0pbqKiLnRueJoKN5lkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRORs+cD/+DmPd4GlzdyHiEgL0NVwIsEgH9hpfxwKxAEDgSuAsGbmexKY4OV212FOdeBPI4Bhft6HSFt0HHgP2Gt/3B5IBC4Fermkew94B8gELquXx2fAOpfnHYDuwLVAvMvyPOCg/XEYEA0kAT8HftrcigQPnwdLOTk5vPTSS2bm4eGkpKQwfvx45s6dS2RkpDNdSEgI+fn53Hjjjc7nDrGxsfTr14958+ZxzTXXOPMtLy/njTfecKZbvXo12dnZPPbYY0ybNo358+ezdu1adu/eTVRUFCNGjGDBggX07dvXssxLlixh6dKllJaWkpCQwC233ML8+fPdyttUbXlGapsfbnrZpvQGxgK1wD7gTcyTzxUe0p7Gvz91vH/re0/3CxXxve+BFZif4dFAF6AO2I95TpnikvYzzCDJ8b++CGCyS75vAq8Dd9VLlwGMtO+nEtgNrMYMmH7VzPp4qanfQ97ODu6X021WVhZ5eXmcOnWKoqIiJk2aREhICAsWLLDcLi8vj6ysLMrKynjooYe44YYbKC4uplevXmelXbZsGbm5ubz44ovccccdAGzdupXc3FyGDBnC6dOnmTVrFqNHj2bXrl20b9/e4z5feeUVHnzwQVasWMGIESPYu3cvOTk5hISEsGjRIq/qfejQIYYOHRpUN0uVViQMiLE/HgLYgD2YwZKjhSgZ+Bjzk3s/8A2wEfg/4AIgHRiDeZJ7lzOtVY/Y/08CegIVwCbgABCC+YsxC7jQnq5+i1Qe0NW+3+32sl6CeYJsqq+BlzFblC63l283cG+9faYBH2IGjf3s5XK0ru0E/gkcs9e3p319B/v6E8AGe71qgFjM12+QfX1j9RZp7d60/78LaOeyvAtnPgcAXwGnMD/DOzHvXuwpdohx+Z+Beb45ifsPqgtc0sUBqZhzRa0DfgZcdE41OSdNnTnfcdPdpt5exi/BUkREBImJiQCkpqaSmZnJ5s2bGw2W4uPjSUxMJDExkaVLl5KcnMzmzZu5++673dItXLiQOXPm8Oqrr3LTTTc5lxcUFLilW7lyJV26dKGoqIgrr7zS4z4/+OADLrvsMm6//XYAevTowYQJE/joo4+8rvexY8eC6mapLc1xI0jxkQswv/wdSjCDoN/Yn9cA/x/mieku4EfgH5jBwk2YQUkZUI3ZYgUQhRmEOLb7LWa33zZgFWbg0tBZYScw3L6vQ8AbmCfXppwI/wW8htmEf4lFuq8wT7qTgO8wf50mAoPt62uBa4BO9vpuspfDcX58B/g3MBGzO+A7zFY4x7bnUm+R1qIKswVpFO6BkkOUy+PtQH/MHyL9MFuXrBpafsD8ARdC00Y7D8T8fNpo0WBp3rx5XH/99ZZpXG+6G9BgyVVxcTEffPAB3bt392q7qCjzqNbUuN8jZ8aMGbzwwgusX7+eUaNGWeZRUVEBQMeOHRtMM2LECFatWsXHH3/M0KFD+de//sWGDRv49a9/3eA21dXVVFdXO59XVla6rQ/GmzS2hNbUDVdSUsL27dtbdH9NZmAGF/txH9NzAWaTtuNTW4QZCNzEmRPj9ZiDqq/FbG0Jt6eJcclnp30fv8I88YEZTD2BGaz0bqBcXYGr7Y87YbZw/YvGT4Q2zFajX2GelK1E2usQCnQG+tj34QiWXD9WHTHHVv0VMyCMwGw5SsRsgQP3FqNizq3e0ihff/Jbz5kkyHxn/9/YrVJOAruA39mfD8BsPc7CvXu8GnjM/viU/f8wPAdi9YVinifKm5DWh3r27OmX71+/BEvr16+nQ4cOnD59murqakJDQ3nuueeavH1VVRUPP/wwYWFhXHXVVc7lGzduZN26dWzZssU5lqkhdXV13H///Vx22WX069fwGfr222+nrKyMyy+/HMMwOH36NPfccw+zZs1qcJv58+fz6KOPNrk+Enxmz54dfK1gezFPTHWYX+r9OROcwJluMId/YwYGrieuNPu2ZZzpmqrvG8yT6uP1lp/GHJfQkK71nsdgtu5Y+RqzXrdidhE2pgvuv1pjMMvrcBgoBI5invAdg9Ar7NsOwWzBOoIZxF3MmV/L51pvaVTw3jJYPCrG/LGRaH+ehNl99iXuP0jaAXdjnpP2AV9gtux6I6TxJK2BX4KlkSNHsnTpUn788UcWL15MeHg4N998c6PbTZgwgbCwME6cOEHnzp1Zvnw5AwYMcK4fMGAAZWVlzJkzh6FDh9KhQ0PfBpCbm0txcTHvv/++5T4LCwt5/PHHeeGFFxg2bBj79+/nvvvuY968eQ1+mc6cOZOpU6c6n1dWVpKamtpo/SR4NKWp1pea1EXZE/gFZ8Yu1b8K7gIfFaYG875T4zys8zy0z+Sp6b2xK+YuxGz6/wz4CY1f2We1D0e3Y2/gZsxutgrMbrRae5o+wB8wT+wHgL9hBlBjOPd6S6NW0bRYuKlsKAA7J45OlMauMdoOfAu4/uY3MD+nrsFSCGbrEJgtvY5B3p4+Q/XVYY4tPE/uceeXYKl9+/b07m22aa9YsYKBAweyfPly7rzzTsvtFi9eTGZmJnFxcXTu3Pms9cnJyaxevZqRI0eSlZXFxo0biYmJOSvd5MmTWb9+Pdu2bSMlJcVyn7Nnz+bXv/41v/ud2R7Zv39/fvzxR37/+9/z0EMPERp69tk7IiKCiAhdytOa+auptiFN6qK8gDMnpqboDOzADAIcrUulmCc4RzN8GGcHNEmYvyzb4/+r3qKB24CVmFfRjOfcp0IowxzDlYn5KxjMlqb62mNehfNz4FPgLcxgqSXr3cak4/4dKwESjflj4mM8d5edwLxa7TCQg/sYphOYn9N/Y55bPLkceAZzCoLGgqAdmK2/58n0AX4fsxQaGsqsWbOYOnUqt99+u3MskieJiYnOIKsh3bt3Z+vWrc6AqaCgwBkwGYbBlClTyM/Pp7CwkJ49ezZavqqqqrMCorCwMGd+56I1jd3xJa/G5Ujz9ce8oiwfs7uuCvNKlQGc6YKLx2xhKcM8MUbat/tf4FXMK2FiMccV2DAvH3YEIr7SAXPA9kuYA7Zv4dwCpjj7dh9hDhL/FnOAtqt3ME/inTFbm/Zy5sTf0vUWCYTrMacO+Cvm+7wrZivPAcwfDxdhjunr4WHbZMzWpdEN5B2HGRm/i3kRhcMpzPmcXKcO+BDzc9r417BPNWU86rl8R7fI9R/jx49n+vTpPP/88zzwwAPNzi81NZXCwkJGjhzJmDFjKCgoIDY2ltzcXF555RXWrVtHTEwMR48eBSAuLq7BIO2Xv/wlixYtYtCgQc5uuNmzZ/PLX/7SGTQ1VadOnYiOjm7ypYsizdIO+DVmgPRX3KcOcBiMOXj5L5gtUI6pA+4A3sYc31ONGTj0xH9zHzmucFsJrMXsRvNWe+BGYAtmwJSEeVJ3nSU8DLNe5ZivRxpmcAbm69XS9RZpaR0xxxltw7wa7QfMFqdumAO412C2EHmSjhnkWF07dSmwHHO6EkfHzXbOTCkSZd/XeHzbN9tETR2PGh0dTUJCYyPhz2iRYCk8PJzJkyezcOFC7r333gbnPPJGSkqKW8C0adMmli41751w9dVXu6XNy8sjJyfHYz4PP/wwISEhPPzww3z99dd07tyZX/7ylzz22GMe01tJTU3FZrO16UkpFSieo5vOcX1XzOb0hrTnzFQDrmIa2Wf9dXd4SNPYrOD184jBfUK8kbjP0+SpPNfVe97f/ufqEZfHV9n/GtJYvUXOBzGY4x9/4WHdDIvtLudMIDUI93mZHFJx/8x5OjcEUFOn7nFMSln/avaG+DxYWrlypcflDz74IA8++KDzef0ursa6vDzlm5yczN69e5uchyfh4eHMmTOHOXPmeL2tJ2lpaV7NCioiIiK+4a+pe3QjXRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQstMg8SyKWWnJaqkBPgRXo/Ys0RO/N4KFjEXSvgYIlCZiEhAQioyI5ufZky+44BOLj41t0l/Hx8eY929a26G5FvBIZFenVrMbie+0i2lGztibQxQgKwfR+VLAkAZOWlsae3XtadMZzxyzjSUlJLbZPwNyf0fTZZf1tw4YNzJ49O2jK4xCM5fK2TMFYB4fGyuaY1VgCZ+2atS1+fgpWwfR+VLAkPtdWbyQs0tqVlZV59ePF3591X+fvuNl3MJ6jHGUrLy9XsGTnzfvR38c0xDiXe4SIm8rKSuLi4qioqCA2NjbQxQmY0tJS0tPTqaqqCnRRLEVHR2Oz2Vr0F0swvjZhYWHU1tYGuhhnCcZyeVumYKyDg6/LFh0aiq2uDl9+mkqB9NBQqurqfJirqS0dm7bmXM7tTf3+VsuS+ExaWlqruJFwIJp2g/G1qa6uJiIiItDFOEswlsvbMgVjHRx8WjabjYTsbJ8GSgBpgK2ujrJVq8DHXZlt5ti0Qf48tytYEp/SjYQbptdGpOnSgLT0dPDDTVFFvKV5lkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxEB7oAogEs9LSUsrKynySV3V1NRERET7JyxeCrTwOwVgub8sUjHVw8GnZbDaqAX/UtBqIsNl8n29bOTZtUEJCAmlpaX7JW8GSSANKS0tJT0+nqqrKJ/mFhYVRW1vrk7x8IdjK4xCM5fK2TMFYBwdfly0sNJTaujqf5eeWb3a27/NtQ8emrYmOjsZms/klYFKwJNKAsrIyqqqqWLVqFenp6c3Ka8OGDcyePdsneflCsJXHIRjL5W2ZgrEODr4um7/q2try9YVgLltrYLPZyM7OpqysTMGSSCCkp6eTkZHRrDxs9u4EX+TlC8FWHodgLJe3ZQrGOjj4umy+ys+X3d31uXbNtIZjEx8fH9iCBBF/dqt5S8GSiIgEjNnd3ZeqqpNnrcv2QTdcdHQkNtueoPnSbcy4m8dRU10T6GIEhcioSPbsDo5jp2BJREQCxuzuPsmqVeDr3iebDbKzT/qta8YfaqprYByQEOiSBFgZnFwbPMdOwZKIiARcejoEWc9Y4CQA3QJdCHGleZZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREpNlCQkJ44403ArBfCMBuAcjJyeHGG28MzM6lRelqOBERaVBOTg4vvfQSAOHh4XTs2JEBAwYwYcIELrjgAme6I0eOcOGFFwaqmOLqOPAesNf+uD2QCFwK9HJJ9x7wDpAJXFYvj8+AdS7POwDdgWuBeJflecBB++MwIBpIAn4O/LS5FQkeXgVL9T80KSkpjB8/nrlz5xIZGelMFxISQn5+vjPiDgkJca6LjY2lX79+zJs3j2uuucaZb3l5uduvktWrV5Odnc1jjz3GtGnTmD9/PmvXrmX37t1ERUUxYsQIFixYQN++fRss75dffsmf/vQnioqKOHjwIIsXL+b+++93S3Mu+Vrx50y04htNvVmlzQ838RRpjbKyssjLy6O2tpZvvvmGgoIC7rvvPi666CJnmsTERL/tv6YG2rXzV941tPNX5oHwPbACiARGA12AOmA/8CYwxSXtZ5hBkuN/fRHAZJd83wReB+6qly4DGGnfTyWwG1iNGTD9qpn18VJTz9vezg7udcuS40Nz6tQpioqKmDRpEiEhISxYsMByu7y8PLKysigrK+Ohhx7ihhtuoLi4mF69ep2VdtmyZeTm5vLiiy9yxx13ALB161Zyc3MZMmQIp0+fZtasWYwePZpdu3bRvn17j/usqqqiV69ejB8/nj/84Q8e05xLvg05dOgQQ4cO9dmNV8U/dLNKEe9EREQ4g6Hk5GQyMjK49NJLGTVqlDON64/kmpoapk6dypo1a/j+++/p2rUr99xzDzNnzgTMH5VTpkxhy5Ytzu2PHTuzv0ceMbvWJk+Gxx6Dgwehrg727YM774SPP4ZeveCZZ84u66FDMG0avPUWGIa57PDhw263N5k6dSrXXnstzz//PBEREZSUlHDo0CGmTZvGW2+9RWhoKFdccQXPPPMMPXr0AKC2tpbp06ezYsUKwsLCuPPOOzEcOwgmb9r/3wW4xoBdgEEuz78CTmEGOTuBUsBT7BDj8j8D2AicxAzGHC5wSRcHpGLOFbUO+BlwES2mqbO+O26629Tby3gdLLl+aFJTU8nMzGTz5s2NBkvx8fEkJiaSmJjI0qVLSU5OZvPmzdx9991u6RYuXMicOXN49dVXuemmm5zLCwoK3NKtXLmSLl26UFRUxJVXXulxn0OGDGHIkCEAPPjggx7TnEu+DTl27JjPbrwq/uHNzSodaUXkbNdccw1paWmUlpaete7ZZ5/lH//4B3//+99JS0vj0KFDHDp0CIC6ujrGjh1Lhw4d2Lp1K8XFxeTk5DBzJlx77Zk89u+HNWtg7VoICzODpXHjoGtX+OgjqKiAeh0FnDoFY8bA8OHw3nuwdy/ccgtMnjyZ0aNHO9N98sknpKWlsXnzZvt2pxgzZgzDhw/nvffeIzw8nD//+c9kZWXx+eef065dO55++mlWrlzJihUrSE9P5+mnnyY/P9/ZQxIUqjBbkEbhHig5RLk83g70x+w664fZumTV0PIDYANCaNpo54HAJvs2LRgszZs3j+uvv94yjetNd/0WLLkqLi7mgw8+oHv37l5tFxVlHrGaGvf738yYMYMXXniB9evXu/1i8aSiogKAjh07erXvxvgi32C8SaOYvLmRprrhRKx169bNY7BUWlpKnz59uPzyywkJCXH7jtiyZQtffPEFJSUlpKamOodpFBXBJ5+A/fctNTXwt79B587m87fegt27YdMm6Gaf3frxx+G6687s97XXzKBq2TJz4PepU+byo0ePUlhY6EwXFRXFsmXLnN1vq1atoq6ujmXLljnLk5eXR3x8PIWFhYwePZolS5Ywc+ZMxo0bB8CLL77Ipk2bmvsS+tZ39v+N3SrlJLAL+J39+QDMsUdZmF1vDtXAY/bH9teSYXgOxOoLBToB5U1I60M9e/b0y/ev18HS+vXr6dChA6dPn6a6uprQ0FCee+65Jm9fVVXFww8/TFhYGFdddZVz+caNG1m3bh1btmxpNFKvq6vj/vvv57LLLqNfv37eVqHZ+VZXV1NdXe18XllZ6bMyiIi0Fg11Q+Xk5HDttdfSt29fsrKyuOGGG5wtOzabjdTUVFJTU922iYkx7+XmCJa6dz8TKJnbQWrqmUAJzBYkVzt3mi1SMfYuobo6839NTQ0HDhwgNjYWgN69e7uNU9q5cyf79+8nxrGh3cmTJzlw4AAVFRUcOXKEYcOGOdeFh4dzySWXBGdXXGOKgY6Yg77BHJAdB3yJ2dXm0A64G3Ms0j7gC8DbhrSQxpO0Bl4HSyNHjmTp0qX8+OOPLF68mPDwcG6++eZGt5swYQJhYWGcOHGCzp07s3z5cgYMGOBcP2DAAMrKypgzZw5Dhw6lQ4cODeaVm5tLcXEx77//vrfFt9TUfOfPn8+jjz7q032LiLQ2hw8f9rg8IyODkpISNm7cyNtvv82tt95KZmYmq1evbnLeXg4ZBeCHH2DwYHj5ZfP5l1/CjTdCfn4+V199NevXrwdwuyDJ3O4HBg8ezMuODV10do3Ygp2jQ6Sxa4y2A98Crl9jBmZXnGuwFILZOgTQmTODvMc1oSx1wDHOm3vceT3PUvv27enduzcDBw5kxYoVfPTRRyxfvrzR7RYvXsyOHTs4evQoR48eZdKkSW7rk5OTKSws5OuvvyYrK4vjx497zGfy5MmsX7+ed999l5SUFG+L3yBv8p05cyYVFRXOP0dfvIhIW/HOO+9YnvtiY2O57bbb+Otf/8prr73GmjVr+O6770hPT3cbw+Rw/Dj81OJS8/R0c/D2kSNnlv3zn+5pMjLMQeBdukDv3mZLFJjja+Pi4hrMOyMjg3379tGlSxd69+7t9hcXF0dcXBxJSUl89NFHzm1Onz5NUVFRwwUOhGigN/AxUONh/QngG+AwkAPc4/KXAxwC/m2R/+WYrVKeY2R3OzC7+86T6QOaNWYpNDSUWbNmMXXqVG6//XbnWCRPEhMT6d27t2V+3bt3Z+vWrYwcOZKsrCwKCgqczaKGYTBlyhTy8/MpLCykZ8+ezSm607nkGxERYXnpuca6BK+SkhKgacfIkVakrauurubo0aNuUwfMnz+fQYMG8dlnn52VftGiRSQlJTFo0CBCQ0N5/fXXSUxMJD4+nszMTPr378/EiRNZsmQJxcXFgNkidMklDZchMxN+8hOYNAmefBIqK+Ghh9zTTJxorhs7FubOhfJyc/nChQt56qmnGsx74sSJPPnkk4wdO5a5c+eSkpLCwYMHWbt2LX/84x9JSUnhvvvu44knnqBPnz5cfPHFLFq0iHLHDoLJ9ZhTB/wV80q3rpitPAeATzEHWycDPTxsm4zZujTawzowu+rSgXeBiS7LT2HO5+Q6dcCHwCWAb76qm6ykpITt27dbpjmX7+hmT0o5fvx4pk+fzvPPP88DDzzQ3OxITU2lsLCQkSNHMmbMGAoKCoiNjSU3N5dXXnmFdevWERMTw9GjRwGIi4trMEirqalh165dzsdff/01O3bsoEOHDs7A7VzybUinTp2Ijo5u8qWLEhhhYWE6RiJeKCgoICkpifDwcC688EIGDhzIs88+S3h4OL/5zW/OSh8TE8PChQvZt28fYWFhDBkyhA0bNhAaanZmrFu3jilTprhdcTx/vnUZQkMhP9+cOmDoUOjRA559FrKyzqSJjoZt22DGDPPKOcdw0pqaGud4JU+io6PZtm0bM2bMYNy4cRw/fpzk5GRGjRrl3G7atGkcOXKESZMmERoaym9/+1tuuukm50VBQaMj5jijbZhXo/2A2eLUDXMA9xrMFiJP0jGDHKvrqy4FlgP/Bzg6Ybbb/8Iwr7jrBoy359fCZs+e3aSrmKOjo0lIaGwk/BnNDpbCw8OZPHkyCxcu5N577/V6biJPUlJS3AKmTZs2sXTpUgCuvvpqt7R5eXnk5OR4zOfw4cMMGnRmYomnnnqKp556iquuusp5ZcS55NuQ1NRUbDabJqUMct5MSqmgStq6lStXsnLlSo/rXMf4uA50vuuuu7jrrvozF56RlpbGunXm9NDbt29n8ODBdOp0Zv0jj5h/9f3kJ+aUAK7qj69OTAT73Mls3262WD388MNuwdKiRYvOyjsxMdE56bIn4eHhLFmyhCVLljSYJmjEAL+w/9U3w2K7yzkTSA3CfV4mh1TgEZfnd5xD+fyoqVP3OCalbOoFWl4FSw19YB588EG3eYzqXx3Q2NUCnvJNTk5m7969Tc7Dkx49ejS6na+vZEhLS/NqVlARERHxDX9N3aMb6YqIiIhYULAkIiIiYkHBkoiIiIgFBUsiIiIiFhQsiYiIiFhQsCQiIiJiodnzLImIiDSXP2580GpvpqCp+oLuNVCwJCIiAZOQkEB0dCTZ2Sf9kn90dKRXMzUHWruIdtSs9XRjt7YnMip4jp2CJRERCZi0tDRstj1udz7YsGEDs2fPbvJszFYcMzW3FmvXrCUpKSnQxQgKwXTsFCyJNMIXN0b25ga+LSHYyuMQjOXytkzBWAcHX5ctmOvqUFZW5gzEgrm8jrKVl5crWLJzPXaN8fcxDTF8fb+PNqiyspK4uDgqKiosb9YorUtpaSnp6elUVVX5JL+wsDBqa2t9kpcvBFt5HIKxXN6WKRjr4ODrsoWFhlJbV+ez/Pyebxs6Nm1NdHQ0NpvNq9aopn5/q2VJpAFm94Dvbozc1Bv4tpRgK49DMJbL2zIFYx0cfFo2m43q7Gz8UdPqujoiVq2CZnbDnZVvWzk2bZA/u+0ULIlY0I2RRQIoPR38cFNUEW9pniURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELIQHugAicn4oLS2lrKys2flUV1cTERHhgxL5jrdlCsY6OPi0bDYbCUCab3JzUwqU2Ww+z7fNHJs2KCEhgbQ0f7wbFSyJiA+UlpaSnp5OVVVVs/MKCwujtrbWB6XyHW/LFIx1cPB12aJDQ7HV1fk0YCoF0kNDqcrO9mGuprZ0bNqa6OhobDabXwImBUsi0mxlZWVUVVWxatUq0tPTzzmfDRs2MHv27Gbn40velikY6+Dg67LZbDays7Mpw7etS2VAVV2dz1/DtnRs2hrne7GsTMGSiAS39PR0MjIyznl7m73bpbn5+JK3ZQrGOjgEc9k88XU5g7n+jrLFx8cHtiBBxJ/dat5SsCQiIhIkxt08jprqmkAXIyhERkWyZ/eeoAiYFCyJiIgEiZrqGhgHJAS6JAFWBifXnvRbt5q3FCyJiIgEkwSgW6ALIa40z5KIiIiIBQVLIiIiIhYULImIiIhYULAkIiIiYkHBkoiIiIgFXQ0nIiJyPjkOvAfstT9uDyQClwK9XNK9B7wDZAKX1cvjM2Cdy/MOQHfgWiDeZXkecND+OAyIBpKAnwM/bW5FgoeCJRERkfPF98AKIBIYDXQB6oD9wJvAFJe0n2EGSY7/9UUAk13yfRN4HbirXroMYKR9P5XAbmA1ZsD0q2bWJ0h4FSzl5OTw0ksvmRuGh5OSksL48eOZO3cukZGRznQhISHk5+dz4403Op87xMbG0q9fP+bNm8c111zjzLe8vJw33njDmW716tVkZ2fz2GOPMW3aNObPn8/atWvZvXs3UVFRjBgxggULFtC3b98Gy/vll1/ypz/9iaKiIg4ePMjixYu5//77z0r39ddfM2PGDDZu3EhVVRW9e/cmLy+PSy65xJuXB/DdndfFf4JpCv1A8uV71eaHu8OLyDl40/7/LqCdy/IuwCCX518BpzCDnJ2Ydy/2dFqMcfmfAWwETmIGYw4XuKSLA1Ix54paB/wMuOicanJOmnou8vZ7wOuWpaysLPLy8jh16hRFRUVMmjSJkJAQFixYYLldXl4eWVlZlJWV8dBDD3HDDTdQXFxMr169zkq7bNkycnNzefHFF7njjjsA2Lp1K7m5uQwZMoTTp08za9YsRo8eza5du2jfvr3HfVZVVdGrVy/Gjx/PH/7wB49pvv/+ey677DJGjhzJxo0b6dy5M/v27ePCCy/08pWBQ4cOMXToUJ/ceV38x593pm4tSktLSU9P13tV5HxShdmCNAr3QMkhyuXxdqA/ZtdZP8zWJatT4g+ADQihaaOdBwKb7Nu0YLCUnZ3dpHSO74Gm3ovP62ApIiKCxMREAFJTU8nMzGTz5s2NBkvx8fEkJiaSmJjI0qVLSU5OZvPmzdx9991u6RYuXMicOXN49dVXuemmm5zLCwoK3NKtXLmSLl26UFRUxJVXXulxn0OGDGHIkCEAPPjggx7TLFiwgNTUVPLy8pzLevbsaVmXhhw7dswnd14X//H3nalbi7KyMp++Vx13TBeRAPrO/r+xW6WcBHYBv7M/H4A59igLs+vNoRp4zP74lP3/MDwHYvWFAp2A8iak9aF58+Zx/fXXW6Zx/R7wW7Dkqri4mA8++IDu3bt7tV1UlBne1tS43yxwxowZvPDCC6xfv55Ro0ZZ5lFRUQFAx44dvdp3ff/4xz8YM2YM48ePZ+vWrSQnJ/Mf//Ef3HVX/U7ZM6qrq6murnY+r6ysdFsfjHe0FndtvdvI13dfb+uvZ1vn66Ovd5OfFQMdMQd9gzkgOw74ErOrzaEdcDfmWKR9wBfANV7uK6TxJL7Us2dPv3z/eh0srV+/ng4dOnD69Gmqq6sJDQ3lueeea/L2VVVVPPzww4SFhXHVVVc5l2/cuJF169axZcsW51imhtTV1XH//fdz2WWX0a9fP2+r4OZf//oXS5cuZerUqcyaNYtPPvmE//zP/6Rdu3ZMmjTJ4zbz58/n0UcfbdZ+JbCa2lQrIo3TpylIONoOGhuKuB34FnD9GjMwu+Jc44wQzNYhgM6cGeQ9rgllqQOOcd7c487rYGnkyJEsXbqUH3/8kcWLFxMeHs7NN9/c6HYTJkwgLCyMEydO0LlzZ5YvX86AAQOc6wcMGEBZWRlz5sxh6NChdOjQocG8cnNzKS4u5v333/e2+Gepq6vjkksu4fHHHwdg0KBBFBcX8+KLLzYYLM2cOZOpU6c6n1dWVpKamtrsskjLaetdpY5maBFfWAX48tNkQwHYOYkGegMf47m77ATm1WqHgRzcxzCdAFYC/8YMjDy5HHgGcwqCxoKgHZjdfefJ9AFeB0vt27end+/eAKxYsYKBAweyfPly7rzzTsvtFi9eTGZmJnFxcXTufPaRSE5OZvXq1YwcOZKsrCw2btxITEzMWekmT57M+vXr2bZtGykpKd4W/yxJSUn89KfuRzM9PZ01a9Y0uE1ERAQRERENrleXRPDydfdTa+er92pJSYlP8pHWKR33BgkJoOsxpw74K+aVbl0xW3kOAJ9iDrZOBnp42DYZs3VpdAN5x2Ee7HeBiS7LT2HO5+Q6dcCHwCXAuQ0BPmclJSVs377dMs25nPeaNWYpNDSUWbNmMXXqVG6//XbnWCRPEhMTnUFWQ7p3787WrVudAVNBQYEzYDIMgylTppCfn09hYeE5D8Ku77LLLmPPnj1uy/bu3ev1OCyATp06ER0drV/sQS46OpqEhMZGQJ7fEhIS9F4VOR91xBxntA3zarQfMFucumEO4F6D2ULkSTpmkGM1ZPhSYDnwf4CjvWK7/S8Ms7WqGzAe3zY3NtHs2bObdLGJt98DzZ6Ucvz48UyfPp3nn3+eBx54oLnZkZqaSmFhISNHjmTMmDEUFBQQGxtLbm4ur7zyCuvWrSMmJoajR48CEBcX12CQVlNTw65du5yPv/76a3bs2EGHDh2cgdsf/vAHRowYweOPP86tt97Kxx9/zF/+8hf+8pe/nFPZbTab5lkKcppnCdLS0nz6XlW3nkgQiQF+Yf+rb4bFdpdzJpAahPu8TA6pwCMuz+84h/L5UVOHWDi+B+pfoNWQZgdL4eHhTJ48mYULF3Lvvfc2OOeRN1JSUtwCpk2bNrF06VIArr76are0eXl55OTkeMzn8OHDDBp05mg/9dRTPPXUU1x11VUUFhYC5vQC+fn5zJw5k7lz59KzZ0+WLFnCxIkTPebZmLS0tDb/RSytg96rInK+8dcQC6+CpZUrV3pc/uCDD7rNY2QYhtv6+s+bkm9ycjJ79+5tch6e9OjRo0nb3XDDDdxwww1e5y8iIiLnv6bMwykiIiLSZilYEhEREbGgYElERETEgoIlEREREQsKlkREREQsNHvqABEREfEhTdUXdK+BgiUREZEg0S6iHTVrawJdjKAQGRUZNHdbULAkIiISJNauWUtSUlKgixEUguluCwqWREREgkRSUpJu8h2EFCyJiM+cy928XZWUlPgkH1/ytkzBWAcHX5fN33X0df5t6di0Nf5+3UKMc7mPiLiprKwkLi6OiooKYmNjA10ckRZXWlpKeno6VVVVzc4rLCyM2tpaH5TKd7wtUzDWwcHXZYsODcVWV4cvO0tKgfTQUKrq6nyYq6ktHZu2Jjo6GpvN5lXXXVO/v9WyJCLNlpaWhs1mo6ys+ZewVFdXExER4YNS+Y63ZQrGOjj4tGw2GwnZ2T4NlADSAFtdHWWrVkET7iDvjTZzbNogf45xUrAkIj6RlpYWNIMxpfVLA9LS00HjdyQIaFJKEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRETk3CQkQGSkf/KOjDTzFwkC4YEugIiItFJpabBnD5SV+T7vhAQzf5EgoGBJRETOXVqagho576kbTkRERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJRERERELCpZERERELChYEhEREbGgYElERETEgmbw9gHDMACorKwMcElERESkqRzf247v8YYoWPKB48ePA5CamhrgkoiIiIi3jh8/TlxcXIPrQ4zGwilpVF1dHYcPH8YwDNLS0jh06BCxsbGBLpbfVVZWkpqa2ibq25bqCm2rvm2prtC26tuW6gptq76+qqthGBw/fpxu3boRGtrwyCS1LPlAaGgoKSkpzua82NjY8/6N6qot1bct1RXaVn3bUl2hbdW3LdUV2lZ9fVFXqxYlBw3wFhEREbGgYElERETEgoIlH4qIiGDOnDlEREQEuigtoi3Vty3VFdpWfdtSXaFt1bct1RXaVn1buq4a4C0iIiJiQS1LIiIiIhYULImIiIhYULAkIiIiYkHBkoiIiIgFBUs+8vzzz9OjRw8iIyMZNmwYH3/8caCL5LX58+czZMgQYmJi6NKlCzfeeCN79uxxS3P11VcTEhLi9nfPPfe4pSktLeUXv/gF0dHRdOnShenTp3P69OmWrEqTPPLII2fV5eKLL3auP3nyJLm5uXTq1IkOHTpw8803880337jl0VrqCtCjR4+z6hsSEkJubi7Quo/ttm3b+OUvf0m3bt0ICQnhjTfecFtvGAZ/+tOfSEpKIioqiszMTPbt2+eW5rvvvmPixInExsYSHx/PnXfeyQ8//OCW5vPPP+eKK64gMjKS1NRUFi5c6O+qeWRV31OnTjFjxgz69+9P+/bt6datG7/5zW84fPiwWx6e3g9PPPGEW5pgqG9jxzYnJ+esemRlZbmlOV+OLeDxMxwSEsKTTz7pTNNajm1TvnN8dR4uLCwkIyODiIgIevfuzcqVK70rrCHN9uqrrxrt2rUzVqxYYXz55ZfGXXfdZcTHxxvffPNNoIvmlTFjxhh5eXlGcXGxsWPHDuP666830tLSjB9++MGZ5qqrrjLuuusu48iRI86/iooK5/rTp08b/fr1MzIzM43PPvvM2LBhg5GQkGDMnDkzEFWyNGfOHONnP/uZW13+/e9/O9ffc889RmpqqrFlyxbj008/NS699FJjxIgRzvWtqa6GYRjffvutW103b95sAMa7775rGEbrPrYbNmwwHnroIWPt2rUGYOTn57utf+KJJ4y4uDjjjTfeMHbu3Gn86le/Mnr27GmcOHHCmSYrK8sYOHCg8c9//tN47733jN69exsTJkxwrq+oqDC6du1qTJw40SguLjb+53/+x4iKijL++7//u6Wq6WRV3/LyciMzM9N47bXXjN27dxsffvihMXToUGPw4MFueXTv3t2YO3eu2/F2/awHS30bO7aTJk0ysrKy3Orx3XffuaU5X46tYRhu9Txy5IixYsUKIyQkxDhw4IAzTWs5tk35zvHFefhf//qXER0dbUydOtXYtWuX8V//9V9GWFiYUVBQ0OSyKljygaFDhxq5ubnO57W1tUa3bt2M+fPnB7BUzfftt98agLF161bnsquuusq47777Gtxmw4YNRmhoqHH06FHnsqVLlxqxsbFGdXW1P4vrtTlz5hgDBw70uK68vNy44IILjNdff925zGazGYDx4YcfGobRuurqyX333WdcdNFFRl1dnWEY58+xrf8FU1dXZyQmJhpPPvmkc1l5ebkRERFh/M///I9hGIaxa9cuAzA++eQTZ5qNGzcaISEhxtdff20YhmG88MILxoUXXuhW1xkzZhh9+/b1c42sefpCre/jjz82AOPgwYPOZd27dzcWL17c4DbBWN+GgqWxY8c2uM35fmzHjh1rXHPNNW7LWuOxNYyzv3N8dR7+4x//aPzsZz9z29dtt91mjBkzpsllUzdcM9XU1FBUVERmZqZzWWhoKJmZmXz44YcBLFnzVVRUANCxY0e35S+//DIJCQn069ePmTNnUlVV5Vz34Ycf0r9/f7p27epcNmbMGCorK/nyyy9bpuBe2LdvH926daNXr15MnDiR0tJSAIqKijh16pTbcb344otJS0tzHtfWVldXNTU1rFq1it/+9reEhIQ4l59Px9ahpKSEo0ePuh3LuLg4hg0b5nYs4+PjueSSS5xpMjMzCQ0N5aOPPnKmufLKK2nXrp0zzZgxY9izZw/ff/99C9Xm3FRUVBASEkJ8fLzb8ieeeIJOnToxaNAgnnzySbeui9ZU38LCQrp06ULfvn259957OXbsmHPd+Xxsv/nmG958803uvPPOs9a1xmNb/zvHV+fhDz/80C0PRxpvvqN1I91mKisro7a21u1AAXTt2pXdu3cHqFTNV1dXx/33389ll11Gv379nMtvv/12unfvTrdu3fj888+ZMWMGe/bsYe3atQAcPXrU42vhWBdMhg0bxsqVK+nbty9Hjhzh0Ucf5YorrqC4uJijR4/Srl27s75cunbt6qxHa6prfW+88Qbl5eXk5OQ4l51Px9aVo2yeyu56LLt06eK2Pjw8nI4dO7ql6dmz51l5ONZdeOGFfil/c508eZIZM2YwYcIEtxuO/ud//icZGRl07NiRDz74gJkzZ3LkyBEWLVoEtJ76ZmVlMW7cOHr27MmBAweYNWsW1113HR9++CFhYWHn9bF96aWXiImJYdy4cW7LW+Ox9fSd46vzcENpKisrOXHiBFFRUY2WT8GSeJSbm0txcTHvv/++2/Lf//73zsf9+/cnKSmJUaNGceDAAS666KKWLmazXHfddc7HAwYMYNiwYXTv3p2///3vTfrwtGbLly/nuuuuo1u3bs5l59OxFdOpU6e49dZbMQyDpUuXuq2bOnWq8/GAAQNo164dd999N/Pnz29Vt8v4f//v/zkf9+/fnwEDBnDRRRdRWFjIqFGjAlgy/1uxYgUTJ04kMjLSbXlrPLYNfecEC3XDNVNCQgJhYWFnjc7/5ptvSExMDFCpmmfy5MmsX7+ed999l5SUFMu0w4YNA2D//v0AJCYmenwtHOuCWXx8PD/5yU/Yv38/iYmJ1NTUUF5e7pbG9bi21roePHiQt99+m9/97neW6c6XY+som9VnNDExkW+//dZt/enTp/nuu+9a7fF2BEoHDx5k8+bNbq1KngwbNozTp0/z1VdfAa2vvg69evUiISHB7X17vh1bgPfee489e/Y0+jmG4D+2DX3n+Oo83FCa2NjYJv8wVrDUTO3atWPw4MFs2bLFuayuro4tW7YwfPjwAJbMe4ZhMHnyZPLz83nnnXfOaqb1ZMeOHQAkJSUBMHz4cL744gu3k5PjRP3Tn/7UL+X2lR9++IEDBw6QlJTE4MGDueCCC9yO6549eygtLXUe19Za17y8PLp06cIvfvELy3Tny7Ht2bMniYmJbseysrKSjz76yO1YlpeXU1RU5EzzzjvvUFdX5wwahw8fzrZt2zh16pQzzebNm+nbt2/QddM4AqV9+/bx9ttv06lTp0a32bFjB6Ghoc4uq9ZUX1f/93//x7Fjx9zet+fTsXVYvnw5gwcPZuDAgY2mDdZj29h3jq/Ow8OHD3fLw5HGq+/ocxuzLq5effVVIyIiwli5cqWxa9cu4/e//70RHx/vNjq/Nbj33nuNuLg4o7Cw0O2S06qqKsMwDGP//v3G3LlzjU8//dQoKSkx1q1bZ/Tq1cu48sornXk4LuMcPXq0sWPHDqOgoMDo3LlzUFxeXt+0adOMwsJCo6SkxPjf//1fIzMz00hISDC+/fZbwzDMS1bT0tKMd955x/j000+N4cOHG8OHD3du35rq6lBbW2ukpaUZM2bMcFve2o/t8ePHjc8++8z47LPPDMBYtGiR8dlnnzmv/nriiSeM+Ph4Y926dcbnn39ujB071uPUAYMGDTI++ugj4/333zf69Onjdnl5eXm50bVrV+PXv/61UVxcbLz66qtGdHR0QC4vt6pvTU2N8atf/cpISUkxduzY4fZZdlwd9MEHHxiLFy82duzYYRw4cMBYtWqV0blzZ+M3v/lN0NXXqq7Hjx83HnjgAePDDz80SkpKjLffftvIyMgw+vTpY5w8edKZx/lybB0qKiqM6OhoY+nSpWdt35qObWPfOYbhm/OwY+qA6dOnGzabzXj++ec1dUCg/Nd//ZeRlpZmtGvXzhg6dKjxz3/+M9BF8hrg8S8vL88wDMMoLS01rrzySqNjx45GRESE0bt3b2P69Oluc/EYhmF89dVXxnXXXWdERUUZCQkJxrRp04xTp04FoEbWbrvtNiMpKclo166dkZycbNx2223G/v37netPnDhh/Md//Idx4YUXGtHR0cZNN91kHDlyxC2P1lJXh02bNhmAsWfPHrflrf3Yvvvuux7fu5MmTTIMw5w+YPbs2UbXrl2NiIgIY9SoUWe9BseOHTMmTJhgdOjQwYiNjTXuuOMO4/jx425pdu7caVx++eVGRESEkZycbDzxxBMtVUU3VvUtKSlp8LPsmFOrqKjIGDZsmBEXF2dERkYa6enpxuOPP+4WYBhGcNTXqq5VVVXG6NGjjc6dOxsXXHCB0b17d+Ouu+4664fq+XJsHf77v//biIqKMsrLy8/avjUd28a+cwzDd+fhd9991/j5z39utGvXzujVq5fbPpoixF5gEREREfFAY5ZERERELChYEhEREbGgYElERETEgoIlEREREQsKlkREREQsKFgSERERsaBgSURERMSCgiURERERCwqWRERERCwoWBIRERGxoGBJRERExIKCJREREREL/z/IEYRwHAXGQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "goi=isoseq['RIPK2']\n",
    "dom_cols={'domains':'green','transmembrane':'pink', 'interest':'yellow'}\n",
    "es_region=[[89780160, 89786592]] # region to highlight (red vertical stripe)\n",
    "ax, blocks=goi.plot_domains(source='annotation',\n",
    "                            categories=['domains','transmembrane', 'interest'], \n",
    "                            ref_trids=False, \n",
    "                            trids=goi.filter_transcripts('SUBSTANTIAL and not MONO_EXON'), \n",
    "                            label='name',\n",
    "                            domain_cols=dom_cols, \n",
    "                            include_utr=True, \n",
    "                            max_overlap=3, \n",
    "                            highlight=es_region)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see, the skipped exon does not disturb the ORF, and also does not cover a domain, but the disordered region between the Protein kinase and CARD domains.\n",
    "\n",
    "## Annotate tables with domains\n",
    "\n",
    "To annotate all identified differential spliced events, we first import the table prepared in the [previous tutorial](09_diff_splice.html). During export, the list format of the transcript id is lost, so we need to make sure they are imported correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ast\n",
    "\n",
    "diff_splice=pd.read_csv(f'{path}/demonstration_dataset_differential_events.csv', converters={col:ast.literal_eval for col in ('trA', 'trB')})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the function add_domains_to_table, to annotate the differential events with the domains. \n",
    "This function has a parameter \"modes\", specifying, how the domains are of the event are summarized. \n",
    "Domains of transcripts sets from different columns are aggregated with set operations. \n",
    "Additionally,  there is a parameter \"overlap_only\", to consider only domains overlapping the event region (as defined by the \"start\" end \"end\" columns of the table).\n",
    "For example, to get the union of all domains from the transcripts in both sets trA and trB, the \"mode\" would be \"trA|trB\". For the top hit, the exon skipping of RIPK2, this reveals the disordered region already identified in the plot above. \n",
    "Domains unique to either set can be retrieved with the modes 'trA-trB' and 'trB-trA' respectively, which are evaluated as set differences. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gene</th>\n",
       "      <th>trA|trB domains</th>\n",
       "      <th>trA-trB domains</th>\n",
       "      <th>trB-trA domains</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>RIPK2</td>\n",
       "      <td>{Disordered}</td>\n",
       "      <td>{}</td>\n",
       "      <td>{Disordered}</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ASAH1</td>\n",
       "      <td>{}</td>\n",
       "      <td>{CBAH}</td>\n",
       "      <td>{}</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>RECQL4</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SNHG6</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SMIM19</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "      <td>{}</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     gene trA|trB domains trA-trB domains trB-trA domains\n",
       "0   RIPK2    {Disordered}              {}    {Disordered}\n",
       "1   ASAH1              {}          {CBAH}              {}\n",
       "2  RECQL4              {}              {}              {}\n",
       "3   SNHG6              {}              {}              {}\n",
       "4  SMIM19              {}              {}              {}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from isotools.domains import add_domains_to_table\n",
    "\n",
    "# get domains of the union of sets trA and trB, directly overlapping the event\n",
    "diff_splice=add_domains_to_table( table=diff_splice, \n",
    "                                 transcriptome=isoseq, \n",
    "                                 source='annotation',\n",
    "                                 modes=['trA|trB'], \n",
    "                                 overlapping_only=True,\n",
    "                                 naming='name',\n",
    "                                 insert_after='pvalue' )\n",
    "\n",
    "# get domains unique to sets trA and trB, from the complete transcript\n",
    "diff_splice=add_domains_to_table( table=diff_splice, \n",
    "                                 transcriptome=isoseq, \n",
    "                                 source='annotation',\n",
    "                                 modes=['trA-trB', 'trB-trA'], \n",
    "                                 overlapping_only=False,\n",
    "                                 naming='name',\n",
    "                                 insert_after='pvalue' )\n",
    "pd.set_option('display.max_columns', None)\n",
    "                        \n",
    "diff_splice[['gene','trA|trB domains','trA-trB domains','trB-trA domains']].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}