dev.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e02e2be9-1fad-41db-be3a-77997fbadf29",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap\n",
    "import matplotlib as mpl\n",
    "import upsetplot\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import yaml\n",
    "\n",
    "from snakemake.io import expand\n",
    "\n",
    "\n",
    "# from utils import * "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f17a0d5e-6e8b-4597-88fc-a328693827c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "config_file = 'config.yml'\n",
    "with open(config_file) as f:\n",
    "    config = yaml.safe_load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05eff97e-2064-4620-b8fe-3747f8a02483",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f5679f49-4c15-4bf2-9fd0-71919bde787d",
   "metadata": {},
   "source": [
    "## 230601_2 WHAT THE f"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f9c40e-29f9-47e5-9cb0-8e7c33ee6e08",
   "metadata": {},
   "source": [
    "## 230601 source_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5e5946f0-15e0-4557-88e4-f0af2e3dfc05",
   "metadata": {},
   "outputs": [],
   "source": [
    "studies = ['ad002', 'ad003', 'ad004']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "790d1ed7-5af7-4a08-96df-d9be861f140f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sources = ['vM21']+studies\n",
    "source_df = pd.DataFrame()\n",
    "source_df['source'] = sources\n",
    "source_df.index.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "32b215fe-d3ff-408a-bf2c-c223c2cbc7aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'data/230516/cerberus/ca_vM21_0_annot.h5'"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(config['data']['ca_annot'],\n",
    "               zip,\n",
    "               batch=batch,\n",
    "               source=source_df.loc[0, 'source'],\n",
    "               cerb_run=0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d8ffc3c5-2ac8-41be-87fd-b498df5daba7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'data/230516/cerberus/ca_vM21_0_annot.h5'"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cerb_run = 1\n",
    "expand(config['data']['ca_annot'],\n",
    "                               zip,\n",
    "                               batch=batch,\n",
    "                               source=source_df.loc[cerb_run-1, 'source'],\n",
    "                               cerb_run=cerb_run-1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "11404b41-9338-4da3-92d4-fd6cc8680106",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ad004'"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source_df.loc[source_df.index.max(), 'source']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2e991673-9e89-425c-a1c0-0db56801bec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "batch = '230516'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "3e3ec914-0dcc-467f-91d0-3503c9acde42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['data/230516/cerberus/ca_vM21_0_annot.h5']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(config['data']['ca_annot'],\n",
    "               zip,\n",
    "               batch=batch,\n",
    "               source=source_df.loc[0, 'source'],\n",
    "               cerb_run=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92e7f58-9805-49e4-bd89-f2949c7cd81e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "210c6bf0-0842-4bb2-9892-abd8e27aacdf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7d3dc6da-0111-4706-8e69-a3f2f87fb51e",
   "metadata": {},
   "source": [
    "## 230531 Fixing lapa filter stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3865c566-a3a3-4d10-8512-bbfa92d7d0e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def filt_lapa_gtf(gtf, filt_list):\n",
    "#     \"\"\"\n",
    "#     Filter LAPA GTF using a TALON-style pass list\n",
    "#     \"\"\"\n",
    "#     gtf = pr.read_gtf(gtf).as_df()\n",
    "#     gids, tids = get_ids_from_pass_list(filt_list)\n",
    "\n",
    "#     # first filter on tids\n",
    "#     gtf = gtf.loc[(gtf.transcript_id.isin(tids))|(gtf.Feature=='gene')]\n",
    "\n",
    "#     # then filter on gids\n",
    "#     gtf = gtf.loc[(gtf.gene_id.isin(gids))]\n",
    "\n",
    "#     gtf = pr.PyRanges(gtf)\n",
    "#     return gtf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e9cd10d4-af4c-4531-92bd-964a89d2af4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf = 'data/230516/lapa/ad002/lapa_corrected.gtf'\n",
    "filt_list = 'data/230516/lapa/ad002/lapa_pass_list.tsv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c4cc15a0-e7b5-4517-a5ed-cdac2d4fb43f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtf = pr.read_gtf(gtf).as_df()\n",
    "gids, tids = get_ids_from_pass_list(filt_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a52cfa3c-faa3-468f-ae32-a12840bac64b",
   "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>Chromosome</th>\n",
       "      <th>Source</th>\n",
       "      <th>Feature</th>\n",
       "      <th>Start</th>\n",
       "      <th>End</th>\n",
       "      <th>Score</th>\n",
       "      <th>Strand</th>\n",
       "      <th>Frame</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>...</th>\n",
       "      <th>NNC_transcript</th>\n",
       "      <th>ISM-suffix_to_IDs</th>\n",
       "      <th>ISM-suffix_transcript</th>\n",
       "      <th>NIC_transcript</th>\n",
       "      <th>ont</th>\n",
       "      <th>intergenic_novel</th>\n",
       "      <th>intergenic_transcript</th>\n",
       "      <th>antisense_gene</th>\n",
       "      <th>gene_antisense_to_IDs</th>\n",
       "      <th>antisense_transcript</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>gene</td>\n",
       "      <td>36511866</td>\n",
       "      <td>36528237</td>\n",
       "      <td>.</td>\n",
       "      <td>+</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSMUSG00000001138.13</td>\n",
       "      <td>Cnnm3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>transcript</td>\n",
       "      <td>36511866</td>\n",
       "      <td>36528237</td>\n",
       "      <td>.</td>\n",
       "      <td>+</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSMUSG00000001138.13</td>\n",
       "      <td>Cnnm3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "      <td>36511866</td>\n",
       "      <td>36513153</td>\n",
       "      <td>.</td>\n",
       "      <td>+</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSMUSG00000001138.13</td>\n",
       "      <td>Cnnm3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "      <td>36518951</td>\n",
       "      <td>36519095</td>\n",
       "      <td>.</td>\n",
       "      <td>+</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSMUSG00000001138.13</td>\n",
       "      <td>Cnnm3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "      <td>36520077</td>\n",
       "      <td>36520227</td>\n",
       "      <td>.</td>\n",
       "      <td>+</td>\n",
       "      <td>.</td>\n",
       "      <td>ENSMUSG00000001138.13</td>\n",
       "      <td>Cnnm3</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 44 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Chromosome  Source     Feature     Start       End Score Strand Frame  \\\n",
       "0       chr1  HAVANA        gene  36511866  36528237     .      +     .   \n",
       "1       chr1  HAVANA  transcript  36511866  36528237     .      +     .   \n",
       "2       chr1  HAVANA        exon  36511866  36513153     .      +     .   \n",
       "3       chr1  HAVANA        exon  36518951  36519095     .      +     .   \n",
       "4       chr1  HAVANA        exon  36520077  36520227     .      +     .   \n",
       "\n",
       "                 gene_id gene_name  ... NNC_transcript ISM-suffix_to_IDs  \\\n",
       "0  ENSMUSG00000001138.13     Cnnm3  ...            NaN               NaN   \n",
       "1  ENSMUSG00000001138.13     Cnnm3  ...            NaN               NaN   \n",
       "2  ENSMUSG00000001138.13     Cnnm3  ...            NaN               NaN   \n",
       "3  ENSMUSG00000001138.13     Cnnm3  ...            NaN               NaN   \n",
       "4  ENSMUSG00000001138.13     Cnnm3  ...            NaN               NaN   \n",
       "\n",
       "  ISM-suffix_transcript NIC_transcript  ont intergenic_novel  \\\n",
       "0                   NaN            NaN  NaN              NaN   \n",
       "1                   NaN            NaN  NaN              NaN   \n",
       "2                   NaN            NaN  NaN              NaN   \n",
       "3                   NaN            NaN  NaN              NaN   \n",
       "4                   NaN            NaN  NaN              NaN   \n",
       "\n",
       "  intergenic_transcript antisense_gene gene_antisense_to_IDs  \\\n",
       "0                   NaN            NaN                   NaN   \n",
       "1                   NaN            NaN                   NaN   \n",
       "2                   NaN            NaN                   NaN   \n",
       "3                   NaN            NaN                   NaN   \n",
       "4                   NaN            NaN                   NaN   \n",
       "\n",
       "  antisense_transcript  \n",
       "0                  NaN  \n",
       "1                  NaN  \n",
       "2                  NaN  \n",
       "3                  NaN  \n",
       "4                  NaN  \n",
       "\n",
       "[5 rows x 44 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gtf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1da2844b-61a9-4b54-b3bc-282060c5e67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gtf = pr.PyRanges(gtf)\n",
    "gtf.to_gtf('test')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c637802b-2539-40c8-9f9e-6b0cae100e1a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Chromosome', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand',\n",
       "       'Frame', 'gene_id', 'gene_name', 'gene_status', 'gene_type',\n",
       "       'talon_gene', 'havana_gene', 'level', 'transcript_id',\n",
       "       'transcript_type', 'transcript_status', 'transcript_name',\n",
       "       'talon_transcript', 'havana_transcript', 'source_attr',\n",
       "       'transcript_support_level', 'tag', 'protein_id', 'ccdsid',\n",
       "       'exon_number', 'exon_id', 'talon_exon', 'exon_status', 'ISM_to_IDs',\n",
       "       'ISM_transcript', 'ISM-prefix_to_IDs', 'ISM-prefix_transcript',\n",
       "       'NNC_transcript', 'ISM-suffix_to_IDs', 'ISM-suffix_transcript',\n",
       "       'NIC_transcript', 'ont', 'intergenic_novel', 'intergenic_transcript',\n",
       "       'antisense_gene', 'gene_antisense_to_IDs', 'antisense_transcript'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gtf.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cabf3b7-f5db-474d-b22e-07137c640b5b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6915832-9fcf-4661-9c5a-6a12a7e0e2e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first filter on tids\n",
    "gtf = gtf.loc[(gtf.transcript_id.isin(tids))|(gtf.Feature=='gene')]\n",
    "\n",
    "# then filter on gids\n",
    "gtf = gtf.loc[(gtf.gene_id.isin(gids))]\n",
    "\n",
    "gtf = pr.PyRanges(gtf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc303ba1-c621-46e3-a15d-bf78d60cf98c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89990ed6-7da9-4cbc-89c9-49b7343637ef",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2603c86-7372-4165-989e-22d4849dd6d6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf5a2d2d-2be4-4c26-bfda-cfd72184349c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "aa6f3abc-cb92-4e19-9d99-b9204ff0abf8",
   "metadata": {},
   "source": [
    "## 230517 getting pydeseq2 to work "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7f9465b9-0a71-49b2-b5f8-beebd4243fa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # save swan adata\n",
    "# swan_file = '/Users/fairliereese/Documents/programming/mortazavi_lab/bin/modelad_pipeline/data/230429/swan/swan.p'\n",
    "# ofile = '/Users/fairliereese/Documents/programming/mortazavi_lab/bin/modelad_pipeline/data/230429/swan/swan_gene_adata.h5ad'\n",
    "# save_swan_adata(swan_file, ofile, how='gene')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cf01cc54-d25c-49e0-bd20-db9958c15335",
   "metadata": {},
   "outputs": [],
   "source": [
    "# actual deseq2\n",
    "from pydeseq2.dds import DeseqDataSet\n",
    "from pydeseq2.ds import DeseqStats\n",
    "from pydeseq2.utils import load_example_data\n",
    "import scanpy as sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0bf370e7-8ba0-4bac-abef-86599a1e2707",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_file = '/Users/fairliereese/Documents/programming/mortazavi_lab/bin/modelad_pipeline/data/230429/swan/swan_gene_adata.h5ad'\n",
    "obs_col = 'genotype'\n",
    "obs_conditions = ['5xFADWT', '5xFADHEMI']\n",
    "threads = 4\n",
    "how = 'gene'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "39be5b81-d49f-4c75-bdb1-5a74735e4b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read(adata_file)\n",
    "\n",
    "# todo - rm this\n",
    "adata.obs['genotype'] = adata.obs['dataset'].str.rsplit('_', expand=True)[0]\n",
    "\n",
    "# remove novel genes\n",
    "if how == 'gene':\n",
    "    adata = adata[:, adata.var.loc[adata.var.gid.str.contains('ENSMUS')].index].copy()\n",
    "    \n",
    "# remove entries that don't belong to the groups we want to test\n",
    "adata = adata[adata.obs.loc[adata.obs[obs_col].isin(obs_conditions)].index,\n",
    "              :].copy()\n",
    "\n",
    "# remove unexpressed stuff\n",
    "sc.pp.filter_genes(adata, min_counts=1)\n",
    "\n",
    "# densify matrix \n",
    "adata.X = np.array(adata.X.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9103092a-cb56-4676-b09c-96f1f39f79f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting size factors...\n",
      "... done in 0.01 seconds.\n",
      "\n",
      "Fitting dispersions...\n",
      "... done in 56.09 seconds.\n",
      "\n",
      "Fitting dispersion trend curve...\n",
      "... done in 11.87 seconds.\n",
      "\n",
      "Fitting MAP dispersions...\n",
      "... done in 46.42 seconds.\n",
      "\n",
      "Fitting LFCs...\n",
      "... done in 6.42 seconds.\n",
      "\n",
      "Refitting 15 outliers.\n",
      "\n",
      "Fitting dispersions...\n",
      "... done in 0.06 seconds.\n",
      "\n",
      "Fitting MAP dispersions...\n",
      "... done in 0.03 seconds.\n",
      "\n",
      "Fitting LFCs...\n",
      "... done in 0.01 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# with adata\n",
    "dds = DeseqDataSet(adata=adata,\n",
    "                   design_factors=obs_col,\n",
    "                   n_cpus=threads,\n",
    "                   refit_cooks=True)\n",
    "dds.deseq2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2f1a7317-ee6d-4474-93d1-191252b03276",
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_res = DeseqStats(dds,\n",
    "                      n_cpus=threads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a087fd27-9814-4daf-9aec-8b293b876abd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pydeseq2.ds.DeseqStats at 0x7fa108383f10>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stat_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b06ed7e2-a531-4dfe-9bdb-c325daeea58e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running Wald tests...\n",
      "... done in 2.89 seconds.\n",
      "\n",
      "Log2 fold change & Wald test p-value: genotype 5xFADWT vs 5xFADHEMI\n"
     ]
    },
    {
     "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>baseMean</th>\n",
       "      <th>log2FoldChange</th>\n",
       "      <th>lfcSE</th>\n",
       "      <th>stat</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>padj</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gid</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000001.4</th>\n",
       "      <td>75.208923</td>\n",
       "      <td>-0.301879</td>\n",
       "      <td>0.081467</td>\n",
       "      <td>-3.705529</td>\n",
       "      <td>0.000211</td>\n",
       "      <td>0.007123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000028.15</th>\n",
       "      <td>6.509473</td>\n",
       "      <td>0.241976</td>\n",
       "      <td>0.372722</td>\n",
       "      <td>0.649213</td>\n",
       "      <td>0.516201</td>\n",
       "      <td>0.776239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000031.16</th>\n",
       "      <td>0.484195</td>\n",
       "      <td>0.609119</td>\n",
       "      <td>1.567579</td>\n",
       "      <td>0.388573</td>\n",
       "      <td>0.697592</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000037.16</th>\n",
       "      <td>1.336114</td>\n",
       "      <td>-0.282314</td>\n",
       "      <td>0.710100</td>\n",
       "      <td>-0.397569</td>\n",
       "      <td>0.690948</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000049.11</th>\n",
       "      <td>5.782117</td>\n",
       "      <td>-0.403003</td>\n",
       "      <td>0.391095</td>\n",
       "      <td>-1.030446</td>\n",
       "      <td>0.302801</td>\n",
       "      <td>0.619153</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000118383.1</th>\n",
       "      <td>2.859580</td>\n",
       "      <td>-0.271767</td>\n",
       "      <td>0.544775</td>\n",
       "      <td>-0.498862</td>\n",
       "      <td>0.617877</td>\n",
       "      <td>0.839385</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000118385.1</th>\n",
       "      <td>1.606901</td>\n",
       "      <td>-0.822110</td>\n",
       "      <td>0.720619</td>\n",
       "      <td>-1.140838</td>\n",
       "      <td>0.253937</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000118386.1</th>\n",
       "      <td>0.234249</td>\n",
       "      <td>-0.136794</td>\n",
       "      <td>1.958223</td>\n",
       "      <td>-0.069856</td>\n",
       "      <td>0.944308</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000118391.1</th>\n",
       "      <td>0.587604</td>\n",
       "      <td>-0.042947</td>\n",
       "      <td>1.112612</td>\n",
       "      <td>-0.038600</td>\n",
       "      <td>0.969209</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000118392.1</th>\n",
       "      <td>1.559768</td>\n",
       "      <td>-0.525404</td>\n",
       "      <td>0.700771</td>\n",
       "      <td>-0.749751</td>\n",
       "      <td>0.453405</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>30384 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                        baseMean  log2FoldChange     lfcSE      stat   \n",
       "gid                                                                    \n",
       "ENSMUSG00000000001.4   75.208923       -0.301879  0.081467 -3.705529  \\\n",
       "ENSMUSG00000000028.15   6.509473        0.241976  0.372722  0.649213   \n",
       "ENSMUSG00000000031.16   0.484195        0.609119  1.567579  0.388573   \n",
       "ENSMUSG00000000037.16   1.336114       -0.282314  0.710100 -0.397569   \n",
       "ENSMUSG00000000049.11   5.782117       -0.403003  0.391095 -1.030446   \n",
       "...                          ...             ...       ...       ...   \n",
       "ENSMUSG00000118383.1    2.859580       -0.271767  0.544775 -0.498862   \n",
       "ENSMUSG00000118385.1    1.606901       -0.822110  0.720619 -1.140838   \n",
       "ENSMUSG00000118386.1    0.234249       -0.136794  1.958223 -0.069856   \n",
       "ENSMUSG00000118391.1    0.587604       -0.042947  1.112612 -0.038600   \n",
       "ENSMUSG00000118392.1    1.559768       -0.525404  0.700771 -0.749751   \n",
       "\n",
       "                         pvalue      padj  \n",
       "gid                                        \n",
       "ENSMUSG00000000001.4   0.000211  0.007123  \n",
       "ENSMUSG00000000028.15  0.516201  0.776239  \n",
       "ENSMUSG00000000031.16  0.697592       NaN  \n",
       "ENSMUSG00000000037.16  0.690948       NaN  \n",
       "ENSMUSG00000000049.11  0.302801  0.619153  \n",
       "...                         ...       ...  \n",
       "ENSMUSG00000118383.1   0.617877  0.839385  \n",
       "ENSMUSG00000118385.1   0.253937       NaN  \n",
       "ENSMUSG00000118386.1   0.944308       NaN  \n",
       "ENSMUSG00000118391.1   0.969209       NaN  \n",
       "ENSMUSG00000118392.1   0.453405       NaN  \n",
       "\n",
       "[30384 rows x 6 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stat_res.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d2ee7f8b-1222-4bb1-ad3b-d503a8d77b43",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = stat_res.results_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7ccf943b-afdf-40ac-addc-ace807d0511a",
   "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>baseMean</th>\n",
       "      <th>log2FoldChange</th>\n",
       "      <th>lfcSE</th>\n",
       "      <th>stat</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>padj</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gid</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000001.4</th>\n",
       "      <td>75.208923</td>\n",
       "      <td>-0.301879</td>\n",
       "      <td>0.081467</td>\n",
       "      <td>-3.705529</td>\n",
       "      <td>0.000211</td>\n",
       "      <td>0.007123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000028.15</th>\n",
       "      <td>6.509473</td>\n",
       "      <td>0.241976</td>\n",
       "      <td>0.372722</td>\n",
       "      <td>0.649213</td>\n",
       "      <td>0.516201</td>\n",
       "      <td>0.776239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000031.16</th>\n",
       "      <td>0.484195</td>\n",
       "      <td>0.609119</td>\n",
       "      <td>1.567579</td>\n",
       "      <td>0.388573</td>\n",
       "      <td>0.697592</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000037.16</th>\n",
       "      <td>1.336114</td>\n",
       "      <td>-0.282314</td>\n",
       "      <td>0.710100</td>\n",
       "      <td>-0.397569</td>\n",
       "      <td>0.690948</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000000049.11</th>\n",
       "      <td>5.782117</td>\n",
       "      <td>-0.403003</td>\n",
       "      <td>0.391095</td>\n",
       "      <td>-1.030446</td>\n",
       "      <td>0.302801</td>\n",
       "      <td>0.619153</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        baseMean  log2FoldChange     lfcSE      stat   \n",
       "gid                                                                    \n",
       "ENSMUSG00000000001.4   75.208923       -0.301879  0.081467 -3.705529  \\\n",
       "ENSMUSG00000000028.15   6.509473        0.241976  0.372722  0.649213   \n",
       "ENSMUSG00000000031.16   0.484195        0.609119  1.567579  0.388573   \n",
       "ENSMUSG00000000037.16   1.336114       -0.282314  0.710100 -0.397569   \n",
       "ENSMUSG00000000049.11   5.782117       -0.403003  0.391095 -1.030446   \n",
       "\n",
       "                         pvalue      padj  \n",
       "gid                                        \n",
       "ENSMUSG00000000001.4   0.000211  0.007123  \n",
       "ENSMUSG00000000028.15  0.516201  0.776239  \n",
       "ENSMUSG00000000031.16  0.697592       NaN  \n",
       "ENSMUSG00000000037.16  0.690948       NaN  \n",
       "ENSMUSG00000000049.11  0.302801  0.619153  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "390375b8-76de-489d-88b9-8d020a7e302d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8752780f-43f5-4797-b9e6-ee0a70a1587a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cb31fbea-d416-45d0-8fa1-4a76a4f78e89",
   "metadata": {},
   "source": [
    "## idk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b8fad2d8-5939-4674-a303-7594aaa01e43",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_talon_nov_colors():\n",
    "    c_dict = {'Known': '#009E73',\n",
    "              'ISM': '#0072B2',\n",
    "              'NIC': '#D55E00',\n",
    "              'NNC': '#E69F00',\n",
    "              'Antisense': '#000000',\n",
    "              'Intergenic': '#CC79A7',\n",
    "              'Genomic': '#F0E442'}\n",
    "    order = ['Known', 'ISM', 'NIC', 'NNC', 'Antisense', 'Intergenic', 'Genomic']\n",
    "    \n",
    "    return c_dict, order\n",
    "\n",
    "def add_perc(ax, data, feature):\n",
    "    total = data[feature].sum()\n",
    "    ylim = ax.get_ylim()[1]\n",
    "    for p in ax.patches:\n",
    "        percentage = '{:.1f}%'.format(100 * p.get_height()/total)\n",
    "        x = p.get_x() + p.get_width() / 2 - 0.45\n",
    "        y = p.get_y() + p.get_height() + ylim*0.00625\n",
    "        ax.annotate(percentage, (x, y), size = 12)\n",
    "        \n",
    "def plot_read_novelty(df, opref, c_dict, order,\n",
    "                      ylim=None, title=None, \n",
    "                      datasets='all'):\n",
    "    sns.set_context(\"paper\", font_scale=1.6)\n",
    "    \n",
    "    temp = df.copy(deep=True)\n",
    "    \n",
    "    # filter on datasets\n",
    "    if datasets != 'all':\n",
    "        temp = temp.loc[temp.dataset.isin(datasets)]        \n",
    "    \n",
    "    # count number of reads per cat\n",
    "    temp = temp[['transcript_novelty', 'read_name']].groupby('transcript_novelty').count()\n",
    "    temp.reset_index(inplace=True)\n",
    "    temp.rename({'read_name':'counts'}, axis=1, inplace=True)\n",
    "    print(temp)\n",
    "    \n",
    "    # actual plotting\n",
    "    g = sns.catplot(data=temp, x='transcript_novelty',\n",
    "                y='counts', kind='bar', \n",
    "                palette=c_dict, order=order)\n",
    "    [plt.setp(ax.get_xticklabels(), rotation=90) for ax in g.axes.flat]\n",
    "    g.set_ylabels('Reads')\n",
    "    g.set_xlabels('Transcript novelty')\n",
    "    \n",
    "    # add percentage labels\n",
    "    ax = g.axes[0,0]\n",
    "    add_perc(ax, temp, 'counts')\n",
    "    \n",
    "    if ylim:\n",
    "        g.set(ylim=(0,ylim))\n",
    "    \n",
    "    # add title\n",
    "    if not title:\n",
    "        g.fig.suptitle('Reads per novelty category')\n",
    "    else:\n",
    "        g.fig.suptitle('{} reads per novelty category'.format(title))\n",
    "        \n",
    "    # save figure\n",
    "    fname = '{}_read_novelty'.format(opref)\n",
    "    g.savefig(fname+'.pdf', dpi=300)\n",
    "    \n",
    "def plot_transcript_novelty(df, oprefix, c_dict, order, \\\n",
    "                            ylim=None, title=None,\n",
    "                            whitelist=None, datasets='all', save_type='pdf'):\n",
    "    sns.set_context('paper', font_scale=1.6)\n",
    "    \n",
    "    temp = df.copy(deep=True)\n",
    "    \n",
    "    # remove transcripts that are not on whitelist\n",
    "    if whitelist:\n",
    "        temp = temp.loc[temp.transcript_ID.isin(whitelist)]\n",
    "    \n",
    "    # filter on datasets\n",
    "    if datasets != 'all':\n",
    "        temp = temp.loc[temp.dataset.isin(datasets)]        \n",
    "    \n",
    "    # count number of isoforms per cat\n",
    "    temp = temp[['transcript_ID', 'transcript_novelty', 'read_name']].groupby(['transcript_ID', 'transcript_novelty']).count()\n",
    "    temp.reset_index(inplace=True)\n",
    "    temp.drop('read_name', axis=1, inplace=True)\n",
    "    temp = temp.groupby('transcript_novelty').count()\n",
    "    temp.reset_index(inplace=True)\n",
    "    temp.rename({'transcript_ID': 'counts'}, axis=1, inplace=True)\n",
    "    print(temp)\n",
    "    \n",
    "    # actual plotting\n",
    "    g = sns.catplot(data=temp, x='transcript_novelty',\n",
    "                y='counts', kind='bar', \n",
    "                palette=c_dict, order=order)\n",
    "    [plt.setp(ax.get_xticklabels(), rotation=90) for ax in g.axes.flat]\n",
    "    g.set_ylabels('Isoforms')\n",
    "    g.set_xlabels('Transcript novelty')\n",
    "    \n",
    "    # add percentage labels\n",
    "    ax = g.axes[0,0]\n",
    "    add_perc(ax, temp, 'counts')\n",
    "    \n",
    "    if ylim:\n",
    "        g.set(ylim=(0,ylim))\n",
    "    \n",
    "    # add title\n",
    "    if not title:\n",
    "        g.fig.suptitle('Transcript models per novelty category')\n",
    "    else:\n",
    "        g.fig.suptitle('{} transcript models per novelty category'.format(title))\n",
    "        \n",
    "    # save figure\n",
    "    fname = '{}_isoform_novelty'.format(oprefix)\n",
    "    if save_type == 'png':\n",
    "        g.savefig(fname+'.png', dpi=300)        \n",
    "    elif save_type == 'pdf':\n",
    "        g.savefig(fname+'.pdf', dpi=300)        \n",
    "    \n",
    "    plt.show()\n",
    "    plt.clf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b24f418-26d6-4bfa-b53e-5931241fc518",
   "metadata": {},
   "outputs": [],
   "source": [
    "# with just bam tag change + reversing fastq reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a3c52e66-42e4-4c4f-a9fa-dd2cc21e6a7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "read_annot = 'data/230418/talon/annot_talon_read_annot.tsv'\n",
    "filt_ab = 'data/230418/talon/annot_talon_abundance_filtered.tsv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d6b3d1c0-35c3-44f3-90cd-e8cc477ff232",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  transcript_novelty   counts\n",
      "0          Antisense     9725\n",
      "1            Genomic    13959\n",
      "2                ISM  1406787\n",
      "3         Intergenic     4699\n",
      "4              Known  7080683\n",
      "5                NIC   676050\n",
      "6                NNC   538026\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAGXCAYAAADViRn0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABDUklEQVR4nO3dd5xU1fnH8c+Xji4CKlZUUGMDCxGMDRRjBOxiiQYVlMQGEtEoKBYsIdgFjChJFKzEghoSxUqLvQJKfjaaBREUkF6f3x/nzjLMzu7dXZi5s+zzfr3mtTv3nrn3mTuzz957zrnnyMxwzjlXuhpJB+Ccc4XOE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBEmSOSjpRkkrolHYtbR1L/6HNplnQsruqokokyLQmlP5ZKmiLpekn1k47RVR2SDogSaLOkY6mIKOaTk46jOqiSiTLN08A50aMfsBi4EXg2yaBclXMAcAPQLNkwKuwG4OSkg6gOaiUdwAaaZGaPpp5IGgy8C3SQdKCZfZBcaFWfpC3M7Oek43CbluiKb5WZrU46lvKq6meU6zGzNcDY6Okv0tdJaiDpz5I+k7RC0k+SnpO0X0a5GpKukTRO0mxJKyV9K2mEpJ2z7VdSr7TtTpd0HVn+CUmqK+k6SVMlLZH0c/S6B8tTXRBVMQyX1F7SG9E25kXLtslSvo6kqyRNlrQs2t+rktpllGsWbbu/pFMlvStpKfCvmHiGR6/bQtKQ6HitkPShpA6lvKaLpHei2JdIelvSmRllHpW0RtJOWV6/maSFkiZmLG8l6WlJP0Sf2TRJAyVtFvcegIeip2PTqnKGS9pe0ipJT5by2sujsseXtY+obFF0fD+JPov5kt6T1DOtTANJN0fHZG70PmZIulfSlmnljpSUuve4a3oVVGWPiaSDJL0efSbzJT0paedo/+OylI/9HKNy46Jt7CJppKR5wFJgt+i7+24px6tz9J56Zlufd2ZW5R7AkYAB12ZZ92y07jdpy7YAJhM+oKHARcDVwDRgCfDLtLL1gIWEP54rgQuAe6PXfgNsmbG/gdH+3geuAK4BvoqeG9AtrezfomWPRjFcBNwKfApsXY73bcAkQhXDEODiKM61wP+AorSytYBXgVXAw8AlwJ+Aj4HVwPFpZZtF2/4IWBC9pz8AF8TEMzx63VvAi0AvoC/wPbAC2Dmj/E1R+cnAVUAfYEq07Jq0ckdHy/pl2efZ0bruacs6AsuBL4Dros/sr8BKYCJQK61s/+j1zaLnhwAPRMv+HG3/bOCQaP0z0Xsp8fkAU4FZQM2Y49Qwes9G+OdzOXApMAx4Pa3cXsAcwne0d9rnuyr63GpH5bZNOw4T0mI+u5LHpA3h+70AGBB9V54BpgNzgXGV+RyjsuOAedFxegboEX1HtgbujF6zX5Zj9mIUU6Ok842ZVflEOSA64FsDe6d9gDOAOmnl746+IL/K2E4j4GtgbNoyAZtl2Wfqj/fKtGW7A2uAd4C6acu3BL6jZKL8CXhxA963RY/TMpb3jpb3T1v2x2jZKRllawMfAtPSljWLyq4CWlYgnuHR6x7IWH5w6vNJW/aL6Fh9nH58gc2jP7LVrEteNYCZwOdZ9vkq4Z9bg+h5PWB25mcQrTs1iqNr2rL+pCXKaFm3aNmRWfb3m2jdFRnLD8885mUcp3ujspdnWVcj7fc6RMkwo8zvo9efnuX7MDxL+Yoekzei479vRtnBUdlxlfkco+Xjom0MzBLnHoR/8kMylu8c7aPEe0vqkXgAlQp6XaLM9ngJ2D2trAj/FcezLqmmP/4Rfbj1s32JCck0VXYB8HTa+iuzfYGjdf0omSinEf6z7l/J923AZ1mW1wHmA1PSlr1POCPI9p77R9vaIyrbLHr+XAXjGR697hdZ1i0q5VidlaXsudG6P6Ytuzladljasp2iP6CH05YdH5W7OMv7bEJIqo+nlU+992Zpy7pReqIU4azsfxnLR0Tfm51ijlENwj/Ir0hLiuU4trXSvnu7RvHdkeX7MDzLa8t9TIBtorKjs2xne0omyop+juOiZY1LeZ+vRd/d+mnLboxec2hl/k5y8ajqdZTDCf/xOwGXEf6LNgWWpZVJfUHaERJm5uN8oGZUBgBJJ0t6M9rO/LSyDQlniym7RT+nZont0yzL/kioBvhY0syoLu4cSfXK/Y6z7MvMVhL+EHdPW7w3IQFme883RGW2zdjU5xWII920LMt+BLZKe75r9HNKlrKpZbulLRse/eyWtqwrIfEMT1u2d/TzPkq+zx+AzSj5PsvNwl/uMGAvSW0BJDUCTgfGmNnXMZvYGmhMaHhcG7c/SX+Q9BHrf/e+ilZvWeoL11eRY5L6XP4vcyNmNptQDZWuop8jwFwzm19KrPcT/iGcBiCpJuFv8hMze7OU1+RdVW/1/srMXo1+HyPpZUI920hJ7aIveeqfwQTCWUpp5gJIOolQz/k+oS5pFusS70g2oAHMzEYr9NXrQDgrPhLoAtwg6RAzm1vZbWdRA/gMKKsy/JOM50srsyMLjWjZqDLbi7b5VdRgc4akXma2jHDGMoN1DXaw7vPoR+jxkE1pf6Tl9RDhu/MHQv1eF6A+oW5zo5H0R+AeQvXCJYTqmxWEv9MXKf93Lx/HpCLK+l49R6jT/gPwCKFutSlwW+7DKr+qnijXY2b/kzSIUMF8FvA4IQEuIJz6v1rGy1O6EirBjzCz4g9Y0uaEM4N0qf/0+1DyDLJFKTEuAP4ZPZB0EaHyvgfhsjDOPpkLJNUh/Bf/Mm3x54RL1XFWGN0wUseqBSUTdMuMMinDCVUjp0iaSagfuyn6B5iSOgteXs7PNxsrc6XZPEnPAKdJ6kX4o/4GeKEc255HSEr7S6oRc1bZlfCPoEN6OUl7l/qK7CpyTFJXA3tlrpC0PeEqKl1lPsdSmdkqSQ8C10jai3BslxGSZsGoqpfeqS49V6nkbYK3EVqF+0uqFX3hHgX2ldQ1fSOSjpb036h7w0JJE6JV6WeiKddlWfZcVPZPkuqmbXdLQuJL31dNSZmJFiDV13OrLOuy2UPSaRnLehAuX0alLXuYkNj7ZduIpEpfjlbSc4SK+z+lVzVEXVWuJNQ9Pp/xmicJ9WndWFePODyjzEuEluIrJW2XuVNJtdK71pRicfSzrHIPEM4iBwH7Aw+WcSZdLPr+PU64ZL00S3zp36nU9mqkrRdwfRlxZ4u53MfEzH4g9FroJGnfjKJ9s2z7OSr+Ocb5W7TN64DjgKeiE4qCUVXPKFN9Dv9NaMUrZmY/SrqX8CGfCzxISBaHAsMVbvmaSDgDO59w2v9JVHZ/QqvwKcD4qI+dCJfK+xDODtL39YWkOwndbt6Q9AShYeX3hMum7dOKNwBmSxpNaDGcDexA+A+6GnisnO99SvQ+2hG6BLUhJJHPgTvSyg0Cfk34h9EOeJnQqLBTdCx2ZV19U86Z2ZeS/kz4Y3hb0mOEY3s2sC+hK9CMjNcsjs7kziZcvk0ws+kZZZZKOofwx/k/SQ8R6tsaED7jzoTvwvAywnuP8IfaL/pntgSYbmbvpO1ngqSphO/JWuDvFXj71xKqWe6R1J7QsLiScFa2J+FzAniK0F3sJUlPE+oSTyF8p7J5GzhaUh9CFZGZ2chKHJPehEaXiZLuI5wt/xpoRfjOF59xV+ZzjGNmMyS9BPwuWjSsIq/Pi6RbkyrzIK0fJeEyuVvG+u0IdTurCX9gHxMuxa+Ofl9K+LLPIZxtHpPx+vMJCWkZofL7MUK9yQxK9ikToSHpC8KXfzrhS5TqTtQtKleH0J3pbUJ1wApC16SngIPK+b5TZ1TtCV06lhKS38PAtlnK1yS0fL5NaIVeFsX3DHBGWrlmlLOrS8b2hxO1d2RZV+JYRcvPJnRbWRo93iZLC2qWz3q9HgRZyu0VxfN19DnMJdQzDyCtZZosrd7R8q6EhrKVqeOcZR+9onX/qcR3dgtCP83Pos9+fnQcLk4rU4NQbfQ54Xv9LaFBpnG2mAhVES8DP6eOUWWOSVT2YELd79IoticJ/1R/BF6o7OdISMAzynF8Tozewye5yBkb+kg8gA1+A9kT5cPRF6I94aype/TlbB+tbxN9KOcQLn3nRF+Sw/IdfwXfa9Y/YH/k7fhfFH0GJyYdS57eb5Po/Q7Nw746RvvqlfT7zvaoqpfepZK0K+G/3e5mlqqonibpUEIL8FjWdV8YQLib5kvgPOB1SfuZ2Wd5DtsVuKgusQfhsvQ/CYez0Umqb6FnQbpro58v5SGEXoQz04JqxEnZ5BIlcCDhcnhyqAcvVodQpwfruq38xcxS9/F+KOlI4EJCtyDnkNSccJtjJ0Krbk8rRyNOVSKpFvB1VMc+lXCXzTGEPspjibnnfwP2uw2hLvQgwvG9w0rvb5moTTFRploMDyZclqdbGf2cHf3M7Lw9FdglR3G5qukIQj/KHwn3Jg9NNpycWENoze5AuLKqRahjvoVwG2psR/lK2ofQI+Bnwp1OpbXuJ25TTJQfRj93MLOXSynzASGJ7kmobE7Zk3BLVUEys0p34HaVY2bDKbvFvMqzUEn4+wT2O44NuCkhn6pkopRUxLrb9QTsLOkA4CcLXXYeAR6U9CfCnQmNCF1ilprZg2a2SNJfgeslfU1osT6f0Er42/y+G+dcoVPU4lSlRHWJY7OsGmFm3aI6l76E/oU7E7o7fEQYwWRctI1ahNGGuhG6bkwC+prZxJKbdc5VZ1UyUTrnXD5V1VsYnXMubzxROudcjCrXmCPJvLrAOZcDpbbA+xmlc87FqLKJsqioaL1HzZo1ufTSMIrV22+/zW9+8xu23HJLmjRpwumnn87s2bNL3daMGTM49thjady4Mdtttx09e/Zk9eowhOPChQvp0KEDjRo1okuXLqxZs+6mjAsuuIBRo0aVtlnn3CaiyibKxYsXFz++//576tevz+mnnw7A/PnzueCCC5gxYwYzZ86kQYMGnHfeeaVu65JLLmGbbbZh9uzZfPzxx4wfP5777rsPgAceeIBWrVoxZ84cZsyYwbPPPgvAW2+9xXfffUfnzp1z/2adc4mqcnWU2TzzzDNss802tG3bFoBOnTqtt75nz54cccQRpb5++vTp9OzZk3r16rHddtvRsWNHPv300+J1J598MnXr1qVt27ZMmzaNNWvW0Lt3bx5//PHcvSnnXMGosmeU6UaMGMG5555LxiAYxSZMmECLFllnZgDgsssuY+TIkSxdupRvv/2WF198kY4dOwLQsmVLXn31VZYtW8bEiRNp0aIFgwcPplOnTuy6a97GvXXOJSnpcd4qMW6dpZsxY4bVqFHDpk2bZtlMmjTJGjdubBMmTMi63sxs6tSp9stf/tJq1qxpgHXt2tXWrl1rZmbLli2zP/zhD7bvvvtanz597Ouvv7ZWrVrZggUL7MILL7S2bdtav379St22c67KKDXvVPkzykceeYTDDz+c5s2bl1j35Zdf0qlTJwYNGlR8WZ5p7dq1dOzYkc6dO7NkyRLmzZvH/Pnz6dOnDwD16tVj2LBhTJ48mYEDB9K7d28GDBjAY489xtq1axk/fjzvvPMOY8aMyen7dM4lp8onyocffpiuXbuWWD5z5kyOPvporrvuOs4555xSX//TTz8xa9YsevbsSd26ddlqq60477zzeOGFkhPsjRkzBjOjY8eOTJkyhdatWyOJ1q1bM3ny5I36vpxzhaNKJ8o333yTb7/9tri1O+Xbb7/lqKOOomfPnlx00UVlbmPrrbemefPmDB06lNWrV7NgwQJGjBjBfvvtt1655cuX07dvX+655x4Amjdvzrhx41i5ciVvvPGG11c6twmr0olyxIgRdO7cmQYNGqy3/O9//zvTpk2jf//+6/W1TBkwYMB6LeOjRo1izJgxNGnShN13353atWtz9913r7fNAQMG0KVLF5o2bQrAhRdeyLx582jSpAlNmzbllFNOyeE7dc4lKa+jB0maQfYRxF8ws+PKuQ3LZ8zOuWqj1FsY892Psg1hCtWU7QmjjT+ZvXjZmtzeY2PEVGlzr/xrovt3zuVHXhOlmc1Nfy6pO2G+jEolSuecy4fE6igVeod3Bx61ktNkOudcwUiyMec3QHPgb6UVkNRfkqU/8heec84FSSbKPwDvmdmk0gqYWX8zU/ojj/E55xyQUKKMJj4/iTLOJp1zrlAkdUbZDVgBPJHQ/p1zrtzyniijRpzfAyPNbHG+9++ccxWVxHiURwK/AM5OYN/OOVdheU+UZjaWMnrAO+dcoanS93o751w+eKJ0zrkYniidcy6GJ0rnnIvhidI552J4onTOuRieKJ1zLoYnSueci+GJ0jnnYniidM65GJ4onXMuhidK55yL4YnSOedieKJ0zrkYniidcy6GJ0rnnIvhidI552J4onTOuRieKJ1zLoYnSueci+GJ0jnnYniidM65GJ4onXMuhidK55yL4YnSOedi5D1RStpW0oOSfpC0XNJnkk7IdxzOOVdetfK5M0lbAP8F3gFOBL4HmgEL8xmHc85VRF4TJdAH+MHMzk5bNiPPMTjnXIXk+9L7ZOBtSU9El95TJPWVVDPPcTjnXLnlO1HuBvQAvgY6AAOBq4Dr8xyHc86VW74TpYDJZnaVmX1kZo8BA4BLsxaW+kuy9Edeo3XOOfKfKGcDUzOWfQo0ltQgs7CZ9TczpT/yEqVzzqXJd6J8A9gzY9mewE9mtijPsTjnXLnkO1HeCRwo6XpJu0s6EegHDMlzHM45V255TZRm9iFwAqH1+xPgLuBu4JZ8xuGccxWR736UmNlLwEv53q9zzlWW3+vtnHMxPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzjkXwxOlc87F8ETpnHMxPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzjkXwxOlc87F8ETpnHMxPFE651yMvCZKSf0lWcZjdT5jcM65iqqVwD6/BNqmPbcEYnDOuXJLIlGuMbPvE9ivc85VShJ1lLtI+lbSTEmjJO2TQAzOOVdu+U6U7wDdgE7AeUB94C1Ju2QrnK1OM3+hOudckNdEaWYvmtk/zWyymb0OnATMAy4ppXx/M1P6I5/xOuccJNw9yMxWAh8CeyQZh3POlSXRRCmpJrAvMDvJOJxzriz57kd5p6QjJTWT1Bp4AmgGPJDPOJxzriLy3T1oB+BRoAnwE/A+cKiZTcpzHM45V255TZRmdlY+9+eccxuD3+vtnHMxPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzjkXwxOlc87F8ETpnHMxPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxahQopR0sqTz0p43l/SmpAWSnpTUYOOH6JxzyaroGeV1QOO050OA7YBBwKHALRspLuecKxgVTZS7AZMBJDUEjgEuN7MbgKuBkzZueM45l7yKJsqagEW/t4t+fzl6PhPYdiPF5ZxzBaOiiXIq0EXS5sAfgDfNbGm0bgdg3sYMzjnnCkFFE+XNQBfgZ6Aj8Je0dZ2ADyuyMUlnSTJJYyoYh3PO5U2tihQ2s39L2htoBUw2sy/SVr9BVH9ZHpJ2BW4HJlYkBuecy7cKJUoAM5sGTMuyfFh5tyGpNvAE0A9oT2g5d865ghSbKCV1rsgGzWxUOYr9GZhhZiMkta/I9p1zLt/Kc0b5dMbzVKu3siyD0DJeKknHAGcAB8TtWFJ/4IbYCJ1zLofK05jTPO3RFvgauBc4Atgr+vnXaHm7sjYkqQkwHOhqZgvidmxm/c1M6Y9yxOuccxtV7Bmlmc1M/S5pCPCwmV2fVuRzYKKkBYRO58eXsbkWwPbAa1JxzqsRbXs1cKCZTarIG3DOuVyraPego4AJpaybCBwZ8/r3gH0Jl92px7+AN6PfP6tgPM45l3MVbfVeSOg/+WqWdR2j9aUysyXAJ+nLojPRemb2SdYXOedcwiqaKIcAAyRtAzwD/ABsA5xG6Ijeb+OG55xzyatoh/OB0RlgP+BsQmu3gO+AHmY2tKIBmFm3ir7GOefyqTIdzu+X9ACwE6FhZjbwtZlZ2a90zrmqqcKJEiBKirOih3PObdIqnCgl1SUMgLEHUC9jtZnZzRsjMOecKxQVSpSSmhIGv9iRUDe5CqgTrV4RPfdE6ZzbpFS0H+XdwAzCAL0iTP+wNXAZ8C3wq40Ym3POFYSKXnofCvQE5qdeb2Y/AYMl1Sd0Hzp6I8bnnHOJq+gZ5RbAj2a2ltC5PH3qh/eAgzdWYM45VygqmiinEboEAXwKnJu27lTgx40RlHPOFZKKXnr/G/g18E/CmJLPS5pLaMTZFrhq44bnnHPJq+idOf3Sfn9R0mHAyUB94BUze3Hjhuecc8mrVIfzFDN7j1A36Zxzm6yK1lGi4HeS7pE0KpokDEknStpt44fonHPJqlCilLQjYabFBwmTgp0ENIpWHw/03ZjBOedcIajoGeUgwohBewAHsv68OWOJH7jXOeeqnIrWUR5DmO9mlqTMScS+I9za6Jxzm5QK11ESugJlsyWwbANicc65glTRRPkW8PtS1v2OMG+Oc85tUip66X0dMF7SG8CThPrKUyT1Iwy9duhGjs855xJXoTNKM3uXMI/3KuAOQmPONYR5c44k3OLonHOblArXUZrZu2Z2JNAAaBr9PB3ojI947pzbBJXr0lvSwUBXYGfgS2CQmU2TtAq4DTgfqA2MzFWgzjmXlNhEKakTMJpwmT0X+A3wO0m/Ax4HGgNPADeb2ec5jNU55xJRnkvva4APgKZmth0hMb5MGEloKfArMzvHk6RzblNVnkS5FzDAzGYDmNkSwq2KtYG+ZvZBDuNzzrnElSdRbkWYuzvdd9HPLzZuOM45V3jK2+ptpSxfU5GdSeou6WNJCyUtlvSRpHPjX+mcc8kpb4fzsZLWZlk+MWO5mVnDMrYzm9Bp/XNgNWHEoQclzTez0eWMxTnn8qo8ifLGjbUzM3shY9EgSV2BdoSWdeecKzixidLMNlqiTCepBmH+nT0JZ5nOOVeQNmgqiMqQtDMwFahLuPzuYWb/yXcczjlXXpUZZm1DfQccABwE3ADcLenX2QpK6i/J0h95jNM554AEEqWZrTazL83sIzO7DXiKUi69zay/mSn9kd9onXMumTPKTDWAekkH4ZxzpclrHaWkW4BXgZmEucCPBc4B/pjPOJxzriLy3ZizJWEGxx2AxcBnwDlm5qMOOecKVl4TpZldks/9OefcxlAIdZTOOVfQPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzjkXwxOlc87F8ETpnHMxPFE651wMT5TOORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzjkXwxOlc87F8ETpnHMxPFE651yMvCZKSX0lvStpoaQfJb0sqU0+Y3DOuYrK9xnlUcADQNvoMRt4VdIueY7DOefKrVY+d2Zmx6Q/l3Q+cApwDPC3fMbinHPllXQd5WZAbWBJwnE451ypkk6UdwDfAc9mWympvyRLf+Q3POecSzBRSroBOA04xcyWZStjZv3NTOmP/EbpnHN5rqNMkfRn4CLgaDObnEQMzjlXXnlPlJLuBLoCvzazj/O9f+ecq6i8JkpJ9xKS5KnA95K2i1YtNrPF+YzFOefKK991lD2AIuAlQh/K1ONPeY7DOefKLd/9KL0xxjlX5STdPcg55wqeJ0rnnIvhidI552J4onTOuRieKHPo3nvvpXXr1tStW5du3bplLXPTTTchiVdffbXU7TRr1oz69etTVFREUVERxxyzbmyR1157jebNm7PddtsxcuTI4uULFizgl7/8JYsWLdpo78e56iqRO3Oqix122IFrr72Wl156iWXLSt6l+dVXX/HUU0+x/fbbx25r9OjRHH300SWWX3bZZYwePZo1a9bQvn17Tj/9dGrWrMnVV19N3759adCgwUZ5L85VZ35GmUOdO3fm5JNPZquttsq6vkePHtx6663UqVOn0vtYsmQJLVu2ZP/996dOnTr8+OOPvPvuu0yfPp0zzjij0tt1zq3jZ5QJeeqpp6hbty7HHntsucp36dKFtWvX0qpVK26//Xb2339/ALbZZhsmTZoEQI0aNWjcuDEnnXQSI0aMyFnszlU3fkaZgEWLFnHNNdcwaNCgcpV/7LHHmDFjBjNnzqR9+/Z06NCBBQsWAHD//ffzxz/+kQsuuIBHHnmEoUOHcvTRR7N8+XI6dOhA+/btGT9+fA7fjXObPj+jTED//v0555xzaNasWbnKH3bYYcW/X3311YwYMYKJEydywgkncMABBzBu3DgAZs+ezRVXXMFbb73FEUccwT333MMOO+xAu3btmDlzJpLfGOVcZfgZZQJee+01Bg8ezHbbbcd2223H119/zRlnnMGtt95artdLwqzkGMa9e/fmlltuoX79+kyZMoXWrVvTrFkzVq1axdy5czf223Cu2vAzyhxavXo1q1evZs2aNaxZs4bly5dTq1YtXnvtNVatWlVcrk2bNtx111106tSpxDZmzZrF119/TZs2bVi7di1Dhgxh3rx5651lArzyyissX76c448/HoDmzZvz+uuvs9NOO7FixYpSG5Scc/E8UebQLbfcwo033lj8/NFHH+WGG26gf//+65WrWbMmjRs3pqioCICLLroICPWPixYt4uKLL+arr76iXr16HHDAAbz44ovrJb4VK1Zw5ZVX8vzzzxcvGzJkCN27d2fFihXcd9991KxZM4fv1LlNm7JdwhUySZaKucntPRKNZe6Vf010/865jarUSnyvo3TOuRh+6Z1D2176j0T3P2dI90T379ymws8onXMuhidK55yL4YnSOedieKJ0zrkYniidcy6GJ0rnnIvhidI552J4onTOuRieKJ1zLkbeE6WkdpL+JelbSSapW75jcM65ikjijLIImAL0AlYmsH/nnKuQvN/rbWYvAC9AGAko3/t3zrmK8jpK55yL4YnSOediFHSilNQ/avApfiQdk3Ou+inoRGlm/c1M6Y+kY9rUjBw5kr333pvNN9+c3XbbjYkTJ5Yos2LFCnr37s0OO+xA48aNueSSS9ab8+eyyy6jcePGHHLIIXzzzTfFyx9//HF69eqVl/fhXC4VdKJ0ufXKK6/Qp08fHnroIRYtWsSECRPYddddS5QbOHAg77//Pp988gmff/45H374IbfccgsA7777Lh988AHff/89hx9+OAMHDgRg4cKF3H777cXlnKvKkuhHWSTpAEkHEOao2Dl6vnO+Y6nubrjhBq6//noOPvhgatSowY477siOO+5Yotzo0aPp1asXW265JU2aNKFXr148+OCDAEyfPp3DDz+cunXr8utf/5pp06YB0K9fP6688kq22GKLvL4n53IhiTPK1sBH0aMOcGP0+00JxFJtrVmzhvfff5+5c+ey++6707RpU3r27MmyZcuylk+fhM7M+Oabb1i4cCEtWrRg4sSJLFu2jNdee40WLVrw/vvv89lnn/G73/0uX2/HuZzKe6I0s3GZ9Y7Ro1u+Y6nO5syZw6pVq3j66aeZOHEiH3/8MR999FHWS+WOHTsyaNAg5s6dy/fff8/gwYMBWLp0KS1btuTUU0/l4IMPZtasWVx11VX06tWLwYMHM3jwYNq1a0eXLl1YsGBBnt+hcxuP11FWU/Xr1wfg0ksvZfvtt2frrbfm8ssv54UXXihRtl+/frRq1YoDDjiAQw89lJNPPpnatWuz7bbbAtC7d28mTZrEP//5T5588knatWvH2rVrGTZsGK+99hp77713cd1lRRx55JHUq1ePoqIiioqK2HPPPbOW69+/P7Vr1y4uV1RUVFwFsHDhQjp06ECjRo3o0qULa9asKX7dBRdcwKhRoyocl6t+PFFWU40bN6Zp06ZI6zoSpP+ern79+tx77718++23TJs2ja222ooDDzyQGjXW//rMmTOHYcOGcf311/PJJ5+w3377Ubt2bdq0acPkyZMrFee9997L4sWLWbx4MZ999lmp5X77298Wl1u8eHFxo9QDDzxAq1atmDNnDjNmzODZZ58F4K233uK7776jc+fOlYrLVS+eKKux8847jyFDhvDDDz8wf/587r77bo4//vgS5b799lu+++47zIy3336bm2++mRtvvLFEucsvv5z+/fuz2Wab0bx5c9577z0WL17MuHHjsram58P06dNp3749devWpW3btkybNo01a9bQu3fv4ioE5+J4oqzGrrvuOtq0acMee+zB3nvvTatWrejXrx+zZs2iqKiIWbNmAfDVV19x6KGHsvnmm9O1a1cGDhzIMcccs962Xn/9dRYsWMApp5wCwEEHHcRxxx3HTjvtxNixY+nbt2+lYrz66qvZeuutOeywwxg3blyp5UaPHs2WW25JixYtGDp0aPHyli1b8uqrr7Js2TImTpxIixYtGDx4MJ06dUosebuqR+mtmVWBJEvF3OT2HonGMvfKv5a5fttL/5GnSLKbM6R7ovvfUO+88w777LMPderUYeTIkfTs2ZOPP/6Y3Xbbbb1yU6dOpVGjRmy77ba88847nHrqqdx1112cddZZLF++nF69evH2229z7LHH0rNnT0488UTGjh1Lnz59mDp1Ku3atfP+ng5Cd8XsKzxRVl5VT5SvnN80T5Fk95sHv4kvlKZjx44cd9xxXHrppWWWGzhwIO+99x7PPPNMiXWnn3463bt3Z9q0aXz88cc88MADHHPMMVxxxRV07NixQvG4TU6pidIvvV2VIYny/GMvrdyYMWMwMzp27MiUKVNo3bo1kmjdunWlG5tc9eCJ0hWkBQsW8NJLL7F8+XJWr17NY489xoQJE7Ke9T3//PPMnz8fM+Pdd99l8ODBnHTSSeuVWb58OX379uWee+4BoHnz5owbN46VK1fyxhtveH2lK5MnSleQVq1axbXXXkuTJk3YeuutGTJkCM899xx77LEHEydOpKioqLjsyJEj2X333WnQoAHnnnsuffr0oWvXruttb8CAAXTp0oWmTUN1w4UXXsi8efNo0qQJTZs2LW6Eci6bvI9w7lx5NGnShPfeey/rurZt27J48eLi50888UTs9m66af07ZBs2bMjLL7+8YUG6asMTpStY7/412XFSDuoxK9H9u8Lhl97OORfDE6VzzsXwROmcczE8UTrnXAxPlM45F8MTpXPOxfBE6ZxzMTxROudcDE+UzlXSTz/9xCmnnMLmm2/OLrvswuOPP5613NixY2nfvj0NGzakWbNm661bvXo1Z555Jo0aNaJjx478/PPPxesGDBjAXXfdlfP4zIw+ffqw1VZbsdVWW9GnT5/iQUVyOZVGoR+/dJ4onaukHj16UKdOHebMmcNjjz3GxRdfzKefflqi3Oabb87555/P7bffXmLdqFGjkMS8efNo2LAhw4YNA8LI7P/617/o1atXzuMbNmwYzz33HJMmTWLy5MmMHj2aBx54AMjtVBqFfvzSeaJ0rhKWLFnCM888w80330xRURGHH344J554Io888kiJsgcddBDnnHNO1hGKpk+fzpFHHkmtWrVo37598aRovXr14s4776RWrcrdZVyR+EaMGMEVV1xB06ZN2XHHHbniiisYPnx4cXy5mEqj0I9fJk+UzlXC559/Tq1atdhjjz2Kl+2///5Zz4jK0rJlS15//XVWrFjB2LFjadGiBc8++2zx9Bf5iO/TTz9l//33z1ouV1NpFPrxy+SJ0rlKWLx4MVtsscV6yxo2bMiiRYsqtJ1jjz2W5s2b06ZNGxo2bMiZZ57JjTfeyG233Ua/fv1o164dl1xyCStXrsxZfIsXL6Zhw4brlVu8eDFmRvfu3Vm4cCG/+tWvaNu2Lfvvvz+PPPIIl112GRdddBHt2rXj2muvrVBsFY2vLLk6fpk8UTpXCUVFRes1HAD8/PPPNGjQoELbkcTAgQOZPHkyw4YNY+DAgVx00UW89957vP/++4wfP56VK1fy4IMP5iy+zLI///wzRUVFSKJevXoMGzaMyZMnM3DgQHr37s2AAQN47LHHWLt2LePHj+edd95hzJgxOYuvLLk6fpkSSZSS/iDpS0krJE2S5JOVuCpljz32YPXq1XzxxRfFyyZNmkSLFi0qvc0pU6bw5ptvcsEFFzBlyhQOPPBAJFVqXvSKxNeiRQsmTZoUW25jTqVR6McvU94TpaRTgaHAncABwAvA85L2y3cszlXW5ptvTufOnbn++utZsmQJb7zxBs8//zznnHNOibJr165l+fLlrFq1CjNj+fLlJS4FzYyePXsyePBgatSoQfPmzfnvf//LypUrGT9+fIXrAysS37nnnstdd91VPH/7nXfeSbdu3dYrs7Gn0ij045cpiTPKK4FHzGyomf3PzK4GJgN/TCAW5yrtvvvuY9myZWyzzTacddZZDB06lBYtWpSYqmLChAnUr1+fY489llmzZlG/fv0S86I/9NBDtGzZkgMPPBCAzp07s8MOO9CkSRN+/PFHLrjggpzFd+GFF3LCCSew77770rJlS4477jguvPDC9baVi6k0Cv34pcvrdLWS6gBLgHPN7Im05TcDp5vZXuXYhk9XW05VfbpaH+Hc5Vmp09XmeyqIraN9zslY/j2wfZ5jcc4l6IvJZyW6/1/sFz/XUkq+zyh3AL4Ffm1mr6ct7wEMMLOGGeX7AzfkLUDnXLVmZlnPKvN9RjkPWA1sm7F8W2B2ZmEz6w/0z1Uw0WV8qafbSfP4NozHt2E8vnXy2phjZiuBD4DfZKzqCLyRz1icc668kpiu9g7gCUnvAuOAc4H9gd8nEItzzsXKe6I0s6clbQn0AQYBnwEnm9mG9Qh1zrkcSeKMEjMbBgxLYt8Zbkw6gBge34bx+DaMxxfJa6u3c85VRT4ohnPOxfBE6ZxzMTxROudcDE+UzjkXwxOl22RI2kJSie+0pJqStsj2mnyS1FXS+5J+ktQsWnaFpIoPvVMNSRokqcRsYZJ6Sbo7l/uutolSUh1JTSXtnP5IMJ7Xy/NIML5fSHpCUsMs6xpF6zZs0L8NIOlkwl1f9bKsrgd8IOn4vAaVRtKFhDFYnwPqs+5v70cKZIhBSX0llRhySlJ3SVclEVOGU4G3six/AzgtlzuudolS0u6SxgNLgZnA9OgxI/qZlCOB5sAXhE74pT2ScgXwnZktzFxhZguAb4A/5TuoNBcDd5jZ0swVZrYEuA3omfeo1ukFXGBmtxDGO0j5AKj8sN4b1wXA/2VZPhW4KM+xZLM18FOW5QuAJrnccSIdzhP2D8IZxhnAd0ChdCTtD3QDjgVGAA+Z2VdJBpThCEJ8pXkSKDnXaP60AC4tY/04cjjASjnsSkiKmZYBRVmWJyE1ulem76N1SZtB+B5m/l0cCXydyx1Xx0R5IHCQmU1NOpB0ZnYTcJOkXwPdgSmS3iEk9qfNbHmiAcIuZBnhKc0cIMmRdreKWS9gy3wEUoqvgX0IVzHp2hOuIgrBbMLfx4yM5a2BH/IeTUnDgDsl1QVS1VC/Bv4M3JLLHVfHRDmNwvkPXoKZvQa8JqkxcDZwFTBY0i5mVrG5PDeuxYSzotKG/d6VMHp9Ur4hDK7yeSnrD4jKJOV+YJCkZdHz5pKOAm4FrkkurPU8TohxKTA2WnYUcDfwcGJRRczsLklNgLuAOtHilcBdZnZnLvdd7W5hlHQEcD1waaGdVaaTtDdhRKWzgZ+BA6K6tqTieR5YYGZdS1n/CNDQzE7Mb2TF+78b6AQcmHmcJG1OuOx90cx6JxFfFMdNhLre+tGiFcBtZlYQg1NLqk1IiL9l/SqpkUBXM1ud9YV5Fn2e+0RPp+bj76I6JsplQG3CpdhqYE36ejPbLIm4ACQVAWcSLr33A54GHjSz8UnFlCKpLeEs4z7gL2Y2O1q+PeGM6GLgKDObkFB8WwMfERoohwD/i1btQ2jEWQv80szmJhFfiqT6hPrUGsCnSf7zK42kXxDOwAE+MrMvEwynIFTHRNmNMhpwzGxE/qJZR9JDhC4OnxLqJUcmfKldgqTfE5JQHcJZLsAWwCqgVzQqVGIk7USYCrkT6yaKMuBFoIeZZdYPJibq77kvMDPqNeCykPQCcJaZLYx+L5WZHZuzOKpboixUktYS6v/KrA7I5ZehPKJ5j84AfkFIRp8RGpuytZYmIqrf3Z0Q3xdmNj/hkJB0D/CJmf09SpJjgbaEet3jEjwTvw+4yswWR7+XyswuyVNYxaITiF5mtij6vVRmdl7O4qhuiTKqaxsLjDWzSUnHkyJpOOXoqpTLL4PLHUmzgNPM7N2o4/vfgWOAc4CDzaxtQnGNBU4xswXR76UxMzsqX3EVmuqYKO8h9Lval9BRdQKhj91YH2W9dJIOKk85M3s317FkI+nBchQzMyt7svMckbQc2N3MvpH0V2CtmV0qaTfgw8wZSF1hqXaJMkXSVkA7Qj+29sDewHwzy2kP/4qKbqtsAnxsZmviyucwjrWEM96yZr0zM6uZp5DWI+nFslYTLnPrJRjfN8CZZvZfSV8QLnefjXo3vO2JMl7UEHY5oe/kNmTcWWhm+2R73cZQHftRpiwn3Ma4LHqsJcF+dpJ+C2xpZkPTlg0l3FYGMF1SezPL6R0IZWie0H7Lxcw6ZVsedeC/lfD55rRTcoxngMckfQ40Al6Olh8AFEyrsqT2wNGEKaQzE9H5iQS1zt+BDoTeIBPI41111S5RShpAuPT+JaEhYhzwF2C8mWW7jzRfLiXtFkBJRwMXAtcRurr8mXXdcJJQpS49JO1HSJBHEXoRHGdmcxIM6QpCY93OhLPJVLeg7Qgt9YmTdA3hn8knFNbtvSnHAyck0fBV7S69o0vIuYRpc58slC4jkuYCR6camCQNAX5hZh2j58cC95pZIiP0pF16l1iVttzMLNF/vlEXoVuALsC/gL5mVtrdOi5NVD1wU9LdvEoTVVmclMSNItXujBI4hFAneTRwvaQfiFrBCQ063yUUVxGQ3o3lMMIdESmfEs4+knJIGetOIgwVtiJPsZQgqRHhjLsn4S6cdmb2ZlLxZBN1C9qN7Je1iXQPylAEvJp0EGW4ARgg6Vwz+zm29EZU7c4o00mqBRxEuBPmbKBmUmdE0X/LP5rZC9GYj/OAI83sjWj9L4ExZrZNEvFlI+lXwO2EQRP+CgxIqs+ipB8Jo0INIoz5mFWCrfJtCP/4mlGyQSyxRrB0UT/Kz8xsUNKxZCNpM2AU60YLWpW+3htzNjJJOxLOKo+Kfu5CGF5qXIJhPUUYkGBHoCOhjujttPWtyT5WYN5J2h0YSDiTfAI428xKGywjXxpHP/sCfcjeOm9AUgnpfsJVwekUZv0fhOHUrpd0OPAxJRPRbUkElWY44e9gOCHWvB3DandGGZ257Uo40ONY14cy0ZbHqOvDA8AJhNguMLOJaevHEs4ob00oRCRtQxjT8feEqoqrCqXTvqRdylMuqTppSUsIA5sUypBqJUgqa+BqS6p+PEXSYuDYJKopquMZ5R2ExFhQFfxmtgw4t4z17fMYTgmSbiC03H4OdIqGgyskhf4ffzLQlMIZe7IEMyvoLmCEMU/nJbHjandGmS4aVgozWxVXNg+xlHnDf8TM7LicB5NF1Oq9jJj+a0ndi17orfJRd6+BhCH+JlHysrYQBsYtFtWT/2wFlCAknUUYAq6rZZmSJJeq4xllagShqwmX4Ej6ijB0WCIjB0WS7ONXHg9T2GdtBd0qz7oO5v9m/eOYSuSF0JhTE7iW0Ke3EbAHME3Sn4EZZva3BMOD8E+mKTBHkjfm5JKk3oTO238F/hstbgvcJ6lRUi1+hT7YhZl1SzqGspjZO5nLsrXK5zuuNIlWnZTT1UBX4DLCtAspk4HeQNKJcmR8kdyodpfekqYBt5jZgxnLuwPXmNluyUTmNpYsrfLXFkCrfMGLbq/sYWavSFoE7G9m06L70d8wsyTnHEpUtZuuFtiRUM+WaXy0zlVRkraJ+gJOBRoArc3s3EJJkgpzo98u6VlJ20bLjpO0b9KxRXYi+5xDa4G6eY6lVJIOUZhrvHt01ZBz1TFRzgCyNYgcS8nZ51wVEbXKf0m4gaCTmXUolK5LAJLaERpxWhG+a5tHq1qS7DS66b4gHL9MJxDu/06UpCaSxgFvEHqv3AG8Jen1aCqQnKl2dZTAbcD90Z0uqTrKw4CzKIxJ3l3l3EBolZ8LXCHpimyFEhwh/i/ADWZ2e3RZm/I6Zc9Hnk+3A0MUJu8S0FbS+YShzc5JNLJgMNAQ2NfMPgWQ1JLQ0DiIcH9/TlS7RGlm/4gGoLiSMBoJhEu1U83s38lF5jZQobfK70e4TTbTXCCnZ0PlZWaPRF3m+gObAQ8R7ljrYWbPJBlbpBPQMZUkAczsE0k9gPJ0r6u0apcoAczsX4SRZdwmotBb5Qnzom8NZN790oJwJ1ZBiBo5H4wuZWsUWP/O2oTjmGlJtC5nqmMdpXNJeBa4KbpVFcAUpoW9lXCff0Exs3kFliQhVJX9RVKD1ILo95uBiaW+aiOojt2DmhAqgVOjOK83eEIhjOLiNj3RH/S/gQOB+oRb8bYm9LY4LrqFNVGS/kf26gsjzAjwBfB3M0tkKLaom9JLhHrKVONSC2Ah0MHMcjZoTHVMlKMJ8+PcR5ZRXMzsn0nE5aqHaKqFAwlXcx8U0j3z0ej/lxBGOUp14D+I0DI/PPrZnjBrYyJVV9EZeRfC3zCE9oXHzGx5TvdbDRPlz8BRZvZ+0rE4V0gk3Q/8YGbXZyy/AdjOzC6Obmc8xszaJBJkQqpjovwC6GxmU5KOxVUfkroC88zsP9HzvxC6o00FziqETvGS5gNtMoccjOpS3zWzxtHl7ztmtkVCMTYknOVmGyX+4Vzttzq2el8D/FnSOfkegcRVa1cTpqlIjXZ+GWGA4eOAOwkD+iZtLdlnhTwgWgehqmp1/kJaR1Inwi2pWwArWb/azAhdxHKiOibKm0hoBBJXre3EurEoTwZGmdng6E6TQpmnZjjwN0l7sq6O8leEcUhTYyMcASR1NXYX8CRhwri8zphaHRNlWSOQ7Jy3KFx1s5xw/zmEKUhSI/EsYt3tjEm7ktABvhehyw2EPp63EnqKQGi5T6oP8s7ArUlMK11tEqWke8zsMjO7sZT1OxC6ajiXCxOAOyX9lzCn/IvR8j2BbxKLKhJNtHcS8KCZDZS0BUDmbIdm9m0S8UX+S+gO9FW+d1xtEiVwvqTvzWxg5opoLpjXgB/zH5arJi4FhgKnAReb2exoeUdC38BEmdlqSY8Ce0XP8zodbDn9HbhDUlOyjxKfsxk2q02rdzQU/2jgEjN7KG351oQJxlYQug15A4+rliS9TRi4I/HEnU003Udpcjrlb7U5ozSzV6MuGg9L+tHM/iVpS0JF+hrgaE+Srpq7mXDG1gD4gHAPdbECuKUxscnPqs0ZZYqkSwjDSZ1JGCWlHnCkmc1NMi636ZG0FNjFzOZKWkbZk7Jtlr/Isss4Yysxr091vr232pxRppjZfdHl9nOE7hpHeJJ0OXIx8HPa74V+VlLw8/pIOoLQH3U34EQz+yaaLHBaLuf7rjZnlFmmgz2KcE/rerMfJjiwq3OuDJI6A49Fj7OBfaI5fXoQBhbJ2d9udRpmbU7G4wnC7HKZy53b6CRNk7RVluWNognvCkKBz+tzLaEx9ves3+L9JmGKjZypNpfehT4drNvkNSP73N11CXeKJS6a12cMIfG0Zf15fc4HTk0otJQ9gbFZli8AGudyx9UmUTqXhCj5pBwSDTyRUhM4Bkh8QIxIoc/rM49wd86MjOUHkeNj6InSudwaR2jEMcIo55kWEsaALASFPq/Po4TuS2cQjmfdqH/07cD9udyxJ0rncmt7Qvea7wiJKL2HxUozm5/1Vcko9Hl9biAMLvIV4Zh+Ev18BChxx93GVG1avZ1zZZN0H6FTd2fgB0Jir0U4E37RzK5MMLxikpqzbpT4DzPHz8zJPj1ROpcfSQ06W16FOq+PpMbA2WY2JHr+TBRfymqgay7Pzj1ROpcH0aCzIwlDrZUYdLYQ7sxJKbR5fSRdCexvZmdHzxcRBhJJNTgdAow0s/45i8ETpXO5F81wOJEEBp0tL0nnAv80sxUZy+sAZyZ11hsN1nGdmb0SPV9ESJzTouedgX5mdmDOYvBE6VzuSVoC7GdmeR9LsbwkrQG2zxz8Iuoo/0NS93pL+gFonZpXSNL/ESY4Sz3fDfjYzBqUsZkNUp3uzHEuSalBZwuZyH4/+rasu8xNQgOgKPXEzPbKmIytHtk782803j3IufxIbNDZOJJeT4UBPCtpZdrqmoTBfHM24EQ5zAT2JcxYmc1+5LjDuV96O5cHSQ46G0fS0OjXCwkzGaa3bq8k9Kt8KKnxWiXdARwPtMpseZe0OWHszH+b2Z9yFoMnSudyT9IupayqC5xiZrfmM55sJN0A3GFmS2IL51E0LOJH0dMhwGfR73sRTQFMSKLzchaDJ0rn8k/SPsDvgXOAhmZWJ+GQCpqknQi3KXYk1KVCqCoYQxhRaGZO9++J0rn8kFSfMLL+HwjzZU8gDPc3KpdnQ+UlaTOgD3A02TvF75pEXOmizue/iJ5+ma+uVt6Y41yOSWpNSI5nAl8TBnc4COhhZqU1UCThfuBY4HHCvekFdxYV3X2T94YvP6N0LockfUQYGGMk8IiZfRAtX0XoNF0wiTIaAq6zmWUb87Fa836UzuXWvoT5mR5PJckCthiYHVuqGvJE6Vxu7UOYYGy0pP+TdE0ZLeBJuxW4SlK1nW2xNH7p7VweSKoNnExo6T6KcJLSD/ibmf2YYGjFJL1IGGBiCaFzd2an+Go78Z4nSufyTFIzoDvQjdC6PN7MfpNkTACSHiprfXWed8oTpXMJkVQDOA74vZmdlHQ8rnSeKJ2r5rLMeZ+NmdlxOQ+mQHk/Suecz2cfw88onXMuhncPcs65GJ4onXMuhifKTYQkK8ejf9JxbghJMyTdW4nXXSapyvYBlNQs+vxOS1tWpd9TVeONOZuOQzKev0UYu+/xtGXf5C+cnDgFqMyUpJcRpmEtT+tuVXEZm957KlieKDcRZvZ2+nNJALMyl2eUqZ/UXM0VkYrTzD6KL+3cxueX3tWEpG7R5duvJL0YzQp4b7TuMknvSVooaa6kMZJaZLx+uKRPJLWT9KGkpZI+ktQuo9zxkt6VtDja3ofRdKLpZc6NXrtc0jxJL6Tuf46Jc71L77SYOkiaEm3vA0mHpJWZAewC9EirguhWxnEySVdKul7SbEk/SRoZjYOYXm5nSU9KWhAdi7GS2qStfyiaLTBz++1T7y9tWUdJb0bb+UnSg5IalRFj1vckaZCkWVFH9vTyh0ZlDi1tm65sniirn5HAG8AJwD+iZTsBfyXci3wesAJ4U9L2Ga/dLip3N3AqYT6V5yQ1gOJpQ0cB/yNcJp9OuPQvTjIKk9mPAD6MttEd+AJoUo44s9meMI7iHcBvCfcnvyRpm2j9KcD3wNOE6olDgP+UsT0I0wu0AM4HribM13J72ntoAIwnjCnZEziLcHU2TtJeUbEngD0ltcrY9lnAV2b2TrStk6N4PiMcjz8SBs79ZxnxlfaehhE+y8zbIbsD/zOzN2PetyuNmfljE3wQBl39U9rzbtGya2JeVxOoA/wAXJ62fDiwFtg3bdkB0TZPjp6fFj1vUMq2GxIGXHigjP2XGicwA7g3IyYDjkpb1pgwXNhfSntdOY7bexnL7gEWpD3vFR2LlmnLNgfmAsPTjuP3wG1pZWoDPwI3R89FmLjrnxn7OziK4/DoebPo+Wlx74nwz+XJtOdFhKlmr0j6O1mVH35GWf38K3OBpIMlvSxpHrCacEbZBNgjo+hsM5uS9jw16GzT6OdkYA3wuKQTs1w+HgJsRtlniKXGWYqFZpaabhULI2C/TphqobJezng+FWgoKTW3dFvgUzP7JG2/SwiNK22j52uAJ4EzFVUYAx2ALVnXwPYLQhIcKalW6gG8Txia7aBKxD4MOEnSVtHzMwkTmD1ciW25iCfK6me929Uk7UxIDHWAi4HDgTaEuZTrZbx2vRZnM0vN/1wvev454TJ1C+AZYK6k/0hKzbWS+uP9rqJxlmFuKa/NrDaoiMyW9fXeJ+Gs9ftS9rtl2vPHCZfCh0XPzwI+NrP/Rc9T1Q2jCFUG6Y8tgJ0rEfuTwFLCpGUQhnX7l5llO06unLzVu/rJvGe1I+Hy7JTobAwASVtSCWY2BhgTnX0dA9xJSBgHEy47AXYgvqtSee+tzazbhDB0WS5H6v6JMFVqtv0WT3ZlZm9Lmg6cJelD4CTgpoztQKjnfCfL9n6oaGBmtkzSo0B3Sa8QzqxvrOh23Pr8jNLVJySl4kFaJZ0INNiQjZrZYjMbBTxCGOUbQt/OpYQGo42loaSjUk+i1umjWD/xrKTk2fGG+C/QUmHK2dR+NyMMmTYxo+wThEatUwjVDk+krfs/wmRju5vZ+1kes8qIoaz3NAxoCdwXbf+l8r81l42fUbpU/d5wSfcT6iWvJvslbZkkXUi4zHyRcHm9M6HF9RUAM1so6UbgVoXpBp4j/LNuDzxhZu9XIv6fgH8o3HU0H7iGkPjvSSvzP+BoScdE5afbho0q/hDQG/i3pGsJjUd/IjToDMwo+3gU023Af83s69QKMzNJlwH/lLQ5oY5zEeG4dQAGWdQ6nkWp78nMpkh6G2hHaDhauwHv1eFnlNVe1DjTldCCPZrQ6nwWlbjsIzTmNCZ01XkZ+DPwLKGbTWp/t0XPD47WDSck58rsD8Il9iXAVcBThLrWjmaWXsd5DTCL0J3mPUKXo0ozs0XAEcAHwFBCV6a1QHsz+7+Msp8SjssOrH+XVGr9KEJS3A14jJAs+xKqKWaUEUbcexpF+IfxYIXenMvKh1lzVZak4UBrM2uZdCyFRtLrwGozOybpWDYFfunt3CZEUmtCz4X2gA+asZF4onRu0/IeoQ/mX8zsxaSD2VT4pbdzzsXwxhznnIvhidI552J4onTOuRieKJ1zLoYnSueci+GJ0jnnYvw/PkI+xY3eFuAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv(read_annot,sep='\\t')\n",
    "c_dict, order = get_talon_nov_colors()\n",
    "opref = 'test'\n",
    "plot_read_novelty(df, opref, c_dict, order,\n",
    "                      ylim=None, title=None, \n",
    "                      datasets='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a0dfa8b2-7e50-4c0c-9adb-7555e983c460",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16888510"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post minimap reads\n",
    "3633544+4507237+3836377+4911352"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9016cc7e-c2a5-4af3-813d-739c6846abe2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10267670"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post tc reads\n",
    "2242793+2700697+2299489+3024691"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e3618feb-92d0-47fc-8893-84581c4d64a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9729929"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post talon reads\n",
    "len(df.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b9490ac-ce57-4ac8-89a4-c311eb5ff408",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "febb1e69-1d8f-4468-865e-6017e68ff01c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3667217-8e4d-4737-a68d-b26b63f5e1f8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "312815af-5107-48d5-a5b0-2ceb907bcf6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# with just bam tag change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7c45965f-2bfa-42f7-b08f-61512b2c4f00",
   "metadata": {},
   "outputs": [],
   "source": [
    "read_annot = 'data/230418/talon/annot_talon_read_annot.tsv'\n",
    "filt_ab = 'data/230418/talon/annot_talon_abundance_filtered.tsv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "04ea9b20-d736-4930-a43a-c016b418ead4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(read_annot,sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "32c0da51-107a-4180-a80e-150d0f0415bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  transcript_novelty   counts\n",
      "0          Antisense     9726\n",
      "1            Genomic    13939\n",
      "2                ISM  1406832\n",
      "3         Intergenic     4692\n",
      "4              Known  7080673\n",
      "5                NIC   676048\n",
      "6                NNC   538010\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAGcCAYAAACC3IhRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/4UlEQVR4nO3dd5wV1fnH8c9DRxYBFRVEBTUSBcUCdhCMBbBjNCoqqAk2JKBRUCxrI6ixgZUkKjEKsWDLzxILzd4BxagRKdIEFQSk8/z+OLPL3cvdZXede+cu+32/XvPaO+XOPHN397lnzpw5x9wdERGJT42kAxAR2dQosYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWKsBM3Mz2yXpOJJmZi2jz6JW0rHIpk2JNSFmNt3MlpvZUjObZ2YPm1lB0nFVJ2Y2zsx+n8Bx9UW3iVNiTdax7l4A7AXsDVyRbDjZoRJiftDvIXeUWPOAu88DXiYkWADM7AAze8vMFpnZJDPrnLLubDP73MyWmNk0MzsvdX9mdpmZzTWzOWZ2Ttq67mY2NXrvbDP7U6aYzKy3mb1pZneb2WIz+6+Z/SZlfSMz+3t0nNlmdqOZ1Ux77x1m9j1QmGH/hWb2uJn9I4rlMzNrn7J+t6hEuShad1y0fP+ohF8zZdsTzWxy9LqGmQ0ys6/N7PvoGFtkOP5NQEfg7uiq4W4zu8fMbkvb7jkzG1DKZ9TGzF4xsx/MbL6ZXRkt38/M3o5inxvtu060bkL09knRcX8XLT/GzD6J3vOWme2Zcpx9zOzj6HN6wsz+ZWY3pqz/g5n9L4rjOTNrnrLOzewiM/sK+Kqi5yiV5O6aEpiA6cDh0esWwBTgrmh+O+B7oDvhy++IaL5ptP5oYGfAgEOBn4F9onVdgflAW6AB8BjgwC7R+rlAx+h1k6L3ZYivN7AGGADUBn4HLAa2iNY/DTwQHWNr4D3gvLT3XgzUAupn2H8hsCI6x5rAn4F3onW1gf8BVwJ1gMOAJUDraP3XwBEp+3oCGBS9/iPwTvSZ1o1iHBWtaxl9FrWi+XHA71P2sx8wB6gRzW8VfbbbZIi/YfRZXgrUi+b3j9btCxwQnXtL4HOgf8p7i38f0fzewHfA/tFn0Yvw91E3Ov8Z0XnVBnoAq4Abo/ceBiwE9om2Hw5MSDvWK8AWQP2KnKOmX/D/nXQAlQoaOgHPAbOjP5zeldzP4cAbwLIoaUwACnJ0DtOBpVHCcOA1oHG0biDwSNr2LwO9StnXM8Afo9cPAkNT1u1KycQ6EzgP2Hwj8fWO/gEtZdl7wJnANsBKUhImcBowNuW9Mzey/0Lg1ZT53YHl0euOwLyif/5o2SigMHp9I/Bg9Lph9PvbMZr/HPhNyvuaAatZn+RKTawp7z8iet0XeKGU+E8DPi7n77o/8HTKfHpivQ+4Ie09XxC+NDtFf+epv4c3WJ9Y/w7ckrKuIDrflinHOqwy56ip8lNVrQooIJTw+hG+vSvMzI4FHgdGE0oY+wF3A2tjirE8TnD3hkBn4NeE0gPAjsDJ0WXhIjNbBBxCSBKYWTczeye69FtEKPUVvbc5MCvlGDPSjnlStP0MMxtvZgeWEd9sj/77UvbVPIqvNjA3Jb4HCCXXIqkxlGZeyuufgXpRPWBzYJa7r0s79nbR68eAHmZWl1CC+8jdi85zR+DplLg+J/xOtylHPAAjgTOi12cAj5Sy3faEkvMGzGxXM/t3VGXxEzCE9b+fTHYELk37fW9P+Byas+HvIfWzbU7K79jdlxKubrYrZXso/zlKJVXJxOruL7j7YHd/ivCNXIKZ1Tazm8xsppn9HNVdnZSyviYhiQ5097vd/b/u/oW7P+7uy3N4KgC4+3jgYeAv0aJZhBJr45SpgbsPjZLJU9G227h7Y+AFQrUAhMvT7VN2v0Pasd539+MJSfAZwpdLabYzM0uZ34FQip1FKLFulRLf5u7eJvVQ5Tz9TOYA25tZ6t/nDoSSG+4+lZBMugGnExJtkVlAt7TPrp67z85wnEwx/hM43szaAbsRPqNMZgE7lbLuPuC/wK/cfXNClYaVsm3Rvm5Ki3kzdx9F+H2m/x5Sf79zCIkZADNrAGxJ9FlF0s+zvOcolVQlE2s5/B04ilBX1ZZQ7/SYmXWJ1u9D+EddYWYfRjcexprZwcmEC8CdwBHRH/s/gWPN7Cgzq2lm9cyss5m1INS51QUWAGvMrBtwZMp+Hgd6m9nuZrYZcG3RCjOrY2Y9zayRu68GfgJSS4Xptgb6RV9UJxP+CV9w97nAf4DbzGzz6IbRzmZ2aEyfxbuEEuzl0bE7A8cSri6KPEaod+xEqGMtcj9wk5ntCGBmTc3s+FKOM5+05Oju3wLvE0pxT5XxRftvoJmZ9TezumbW0Mz2j9Y1JHy2S83s18AFGznuX4HzLdyYMzNrYGZHm1lD4G1CibuvmdWKzmW/lPeOAs42s72iL90hwLvuPr2UuCtyjlJZSddF/NKJcAOkd8r8ToRksVPadn8n/BEBnEr4Fp8FnEJItMMJpbDWOYp7OtHNq5Rl96XEuD8wHviBkET/D9ghWncR4Z9zEeGfYzRRnVu0fhDhMnsOcE50rrsQkvJLwI+Ef/z3gUNKia838CahZL8Y+BI4MmV9oyjeb6P1HwOnprz3jY2cfyHwz5T5lpSs/2wTnf9iYCpwYtr7d4h+z/+XtrwGcAmhjnIJ4XJ9SCnHODA6rx+BYSn7OCParstGzqEtoW78x+jzLrqB1olQYl0KTASuT/08gPMJJdFFwCnRsq7R72NRtO4JoGG0rj3wSbS/J4AxwNVp+/s6+lv5N9AiZV2J+tyKnqOmyk0WfchVlpmtAM5394ej+ZMJpbZlaZvWAT5393ZmdhqhxHORu9+bsq8pwCvufklOgs9jZtabcGPnkKRjyTUz60S4atjR8/AfxMzeBe5394d+wT7y+hyruk2xwXBR9cYBhNJsqqIbXXOjn1PT1k8lpb5Kqh8zq02oYvhbviScqIrlC0Kzqp7AnoQrj8ruL+/OcVOzKdaxfhT9bO7u/0ubZkbrPiQk3dZp721NuESXasjMdiNcijcj1Hnni9bAJEJslwK/9VDPXWF5fI6blCpZFWDhmfqiZ63fBW4itGv9wd1nmtk/CA2n/0Roe9kYOAj42d0fjPbxF0JbxD8AXxHqIgcA7dz9i9ydjYhsaqpqYu0MjM2waqS7947aQg4i3ETZgXBz4WNCw/lx0T5qEW4q9AY2J5QIBrn7xOxGLyKbuiqZWEVE8tmmWMcqIpKoKtcqwMx0I1NE8kGpT9OpxCoiEjMlVhGRmFXpxFpQUFBiqlmzJhdffDEA77zzDkcccQRbbLEFTZs25eSTT2bu3NKb/k2fPp3u3bvTpEkTtt12W/r27cuaNWsAWLx4MUcddRSNGzemZ8+erF27vgOsPn36MGbMmOyeqIhUKVU6sS5durR4mjdvHvXr1+fkk08G4Mcff6RPnz5Mnz6dGTNm0LBhQ84+++xS93XhhRey9dZbM3fuXD755BPGjx/PvfeGp10feOAB9t57b+bPn8/06dN5+umnAXj77beZM2cOPXr0yP7JikiVUeVuXpXmqaeeYuutt6Zjx44AdOvWrcT6vn37cuihpXe+9M0339C3b1/q1avHtttuS9euXfnss8+K151wwgnUrVuXjh07Mm3aNNauXcuAAQN47LHHSt2niFRPVbrEmmrkyJGcddZZlOy2cr0JEybQpk2bjOsA+vfvz+jRo/n555+ZPXs2L774Il27dgWgbdu2vPrqqyxfvpyJEyfSpk0bhg0bRrdu3dhpp9K65BSRaivp7rUqOoWQS5o+fbrXqFHDp02btsE6d/dJkyZ5kyZNfMKECRnXu7tPnTrV99lnH69Zs6YD3qtXL1+3bp27uy9fvtz/8Ic/+B577OEDBw70WbNm+d577+2LFi3y8847zzt27OiDBw8udd8iskkqNU9tEiXWRx55hEMOOYRWrVptsO5///sf3bp146677iquJki3bt06unbtSo8ePVi2bBkLFy7kxx9/ZODAgQDUq1ePESNGMHnyZIYOHcqAAQMYMmQIjz76KOvWrWP8+PG8++67vPRSpTscEpFNyCaRWP/xj3/Qq1evDZbPmDGDww8/nKuvvpozzzyz1Pf/8MMPzJw5k759+1K3bl223HJLzj77bF544YUNtn3ppZdwd7p27cqUKVNo3749Zkb79u2ZPHlyrOclIlVTlU+sb731FrNnzy5uDVBk9uzZHHbYYfTt25fzzz+/zH1stdVWtGrVivvuu481a9awaNEiRo4cyZ577lliuxUrVjBo0CDuvPNOAFq1asW4ceNYtWoVb775pupbRQTYBBLryJEj6dGjBw0bNiyx/G9/+xvTpk2jsLCwRFvXIkOGDCnRcmDMmDG89NJLNG3alF122YXatWtzxx13lNjnkCFD6NmzJy1atADgvPPOY+HChTRt2pQWLVpw4oknZvFMRaSqqHK9W6mvABHJE6X2FZDTdqxmNp3MQ5+84O5HV2afTW+96BfFFIcFl92TdAgikkdy/YBAB6BmynwzwjApZY1tLyJSpeQ0sbr7gtR5MzuXMAyzEquIbDISu3ll4RGpcwljyy8vZZtCM/PUKbdRiohUXJKtAo4AWgF/LW0Ddy90d0udcheeiEjlJJlY/wC87+6TEoxBRCR2iSRWM9saOJ4ySqsiIlVVUiXW3sBKYFRCxxcRyZqcJ9boptXvgdHuvjTXxxcRybYkOrruDPwKOCOBY4uIZF3OE6u7j6WMR8FERKq6Kt8Ji4hIvlFiFRGJmRKriEjMlFhFRGKmxCoiEjMlVhGRmCmxiojETIlVRCRmSqwiIjFTYhURiZkSq4hIzJRYRURipsQqIhIzJVYRkZgpsYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWEVEYqbEKiISMyVWEZGYKbGKiMRMiVVEJGZKrCIiMVNiFRGJmRKriEjMcp5YzWwbM3vQzL4zsxVm9oWZHZvrOEREsqVWLg9mZpsDbwDvAscB84CWwOJcxiEikk05TazAQOA7dz8jZdn0HMcgIpJVua4KOAF4x8xGRVUBU8xskJnVzLSxmRWamadOuQ1XRKTicp1YdwYuAmYBRwFDgcuBazJt7O6F7m6pU+5CFRGpnFxXBRgwyd0vj+Y/NrNmwJXAtTmORUQkK3JdYp0LTE1b9hnQxMwa5jgWEZGsyHVifRNonbasNfCDuy/JcSwiIlmR68R6G7CvmV1jZruY2XHAYGB4juMQEcmanCZWd/8IOJbQOuBT4HbgDuDGXMYhIpJNub55hbu/DLyc6+OKiOSK+goQEYmZEquISMyUWEVEYqbEKiISMyVWEZGYKbGKiMRMiVVEJGZKrCIiMVNiFRGJmRKriEjMlFhFRGKmxCoiEjMlVhGRmCmxiojETIlVRCRmSqwiIjFTYhURiZkSq4hIzJRYRURipsQqIhIzJVYRkZgpsYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWEVEYpbTxGpmhWbmadOaXMYgIpJttRI45v+AjinznkAMIiJZk0RiXevu8xI4rohITiRRx7qjmc02sxlmNsbMdi9tw0xVB7kMVESkMnKdWN8FegPdgLOB+sDbZrZjpo3dvdDdLXXKXagiIpWT06oAd38xdd7M3gA+By4EBuYyFhGRbEm0uZW7rwI+AnZNMg4RkTglmljNrCawBzA3yThEROKU63ast5lZZzNraWbtgVFAS+CBXMYhIpJNuW5u1Rz4J9AU+AH4ADjI3SflOA4RkazJ9c2r03J5PBGRJKivABGRmCmxiojETIlVRCRmSqwiIjFTYhURiZkSq4hIzJRYRURipsQqIhIzJVYRkZgpsYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWEVEYqbEKiISMyVWEZGYKbGKiMRMiVVEJGZKrCIiMVNiFRGJmRKriEjMKpRYzewEMzs7Zb6Vmb1lZovM7HEzaxh/iCIiVUtFS6xXA01S5ocD2wJ3AQcBN8YUl4hIlVXRxLozMBnAzBoBRwKXuPu1wBXA8fGGJyJS9VQ0sdYEPHrdKXr9n2h+BrBNTHGJiFRZFU2sU4GeZtYA+APwlrv/HK1rDiyMMzgRkaqooon1BqAn8BPQFfhzyrpuwEcV2ZmZnWZmbmYvVTAOEZG8VasiG7v7v81sN2BvYLK7f5Wy+k2i+tfyMLOdgFuBiRWJQUQk31UosQK4+zRgWoblI8q7DzOrDYwCBgNdCC0LREQ2CRtNrGbWoyI7dPcx5djsJmC6u480sy5lHLsQuLYixxcRSVp5SqxPps0XtQqwDMsgtBwolZkdCZwC7LWxA7t7IVCY9n7PuLGISJ4oT2JtlfK6BfAo8Cwh4c4nNLE6GTiOcGOrVGbWFHgYOM3dF1U8XBGR/LfRxOruM4pem9lw4B/ufk3KJl8CE81sEeEhgWPK2F0boBnwmllxgbdGtO81wL7uPqkiJyAikm8q2tzqMGBCKesmAp038v73gT0I1QBF03PAW9HrLyoYj4hI3qloq4DFhParr2ZY1zVaXyp3XwZ8mrosKunWc/dPM75JRKSKqWhiHQ4MMbOtgaeA74Ctgd8S6lcHxxueiEjVU9EHBIZGJczBwBmE1gAGzAEucvf7KhqAu/eu6HtERPJZZR4QuN/MHgC2J9yImgvMcnc1gxIRoRKJFSBKojOjSUREUlQ4sZpZXUKHK7sC9dJWu7vfEEdgIiJVVYUSq5m1IHS2sh2hbnU1UCdavTKaV2IVkWqtou1Y7wCmE562MsJwLFsB/YHZwP4xxiYiUiVVtCrgIKAv8GPR+939B2CYmdUnNMc6PMb4RESqnIqWWDcHvnf3dYSHAVKHYnkfOCCuwEREqqqKJtZphCZWAJ8BZ6WsOwn4Po6gRESqsopWBfwb+A3wL0Kfqs+a2QLCTattgMvjDU9EpOqp6JNXg1Nev2hmBwMnAPWBV9z9xXjDExGpeir1gEARd3+fULcqIiKRitaxYsHpZnanmY2JBgXEzI4zs53jD1FEpGqpUGI1s+0II7E+SBgE8HigcbT6GGBQnMGJiFRFFS2x3kXo0WpXYF9Kjns1lo13dC0issmraB3rkUAvd59pZumDBs4hPOoqIlKtVbiOldC0KpMtgOW/IBYRkU1CRRPr28DvS1l3OmHcKxGRaq2iVQFXA+PN7E3gcUJ964lmNpjQleBBMccnIlLlVKjE6u7vAYcSqgP+Qrh5dSVh3KvOhEdeRUSqtQrXsbr7e+7eGWgItIh+ngz0QCMKiIiUryrAzA4AegE7AP8D7nL3aWa2GrgFOAeoDYzOVqAiIlXFRhOrmXUDnidc9i8AjgBON7PTgceAJsAo4AZ3/zKLsYqIVAnlqQq4EvgQaOHu2xIS6X8IPV39DOzv7mcqqYqIBOVJrL8Ghrj7XAB3X0Z4dLU2MMjdP8xifCIiVU55EuuWwNy0ZXOin1/FG46ISNVX3lYBXsrytXEFIiKyqShvYh1rZj8VTawfTHBi6nIzW1zWTszsXDP7xMwWm9lSM/vYzM4q6z0iIlVNeZpbXRfj8eYSnt76ElhD6GrwQTP70d2fj/E4IiKJ2WhidffYEqu7v5C26C4z6wV0IjTpEhGp8irTu1UszKyGmR0BtAbGJRWHiEjccp5YzWwHM1sKrASeAy529/8rZdtCM/PUKafBiohUQhIl1jnAXsB+wLXAHWb2m0wbunuhu1vqlMM4RUQq5ReN0loZ7r6G0N8AwMdmtivhhtZruY5FRCQbEqtjTVEDqJd0ECIicclpidXMbgReBWYA9YHuwJnAH3MZh4hINuW6KmALwtDZzYGlwBfAme6u7gZFZJOR08Tq7hfm8ngiIknIhzpWEZFNihKriEjMlFhFRGKmxCoiEjMlVhGRmCmxiojETIlVRCRmSqwiIjFTYhURiZkSq4hIzJRYRURipsQqIhIzJVYRkZgpsYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWEVEYqbEKiISMyVWEZGYKbGKiMRMiVVEJGZKrCIiMVNiFRGJmRKriEjMlFhFRGKW08RqZoPM7D0zW2xm35vZf8ysQy5jEBHJtlyXWA8DHgA6RtNc4FUz2zHHcYiIZE2tXB7M3Y9MnTezc4ATgSOBv+YyFhGRbEm6jnUzoDawLOE4RERik3Ri/QswB3g600ozKzQzT51yG56ISMUllljN7Frgt8CJ7r480zbuXujuljrlNkoRkYrLaR1rETO7CTgfONzdJycRg4hItuQ8sZrZbUAv4Dfu/kmujy8ikm05TaxmdjchqZ4EzDOzbaNVS919aS5jERHJllzXsV4EFAAvE9qwFk1/ynEcIiJZk+t2rLr5JCKbvKSbW4mIbHKUWEVEYqbEKiISMyVWEZGYKbGKiMRMiVVEJGZKrDlw99130759e+rWrUvv3r0zbnP99ddjZrz66qul7qdly5bUr1+fgoICCgoKOPLI9b0wvvbaa7Rq1Yptt92W0aNHFy9ftGgR++yzD0uWLIntfESkbIn0FVDdNG/enKuuuoqXX36Z5cs37G/m66+/5oknnqBZs2Yb3dfzzz/P4YcfvsHy/v378/zzz7N27Vq6dOnCySefTM2aNbniiisYNGgQDRs2jOVcRGTjVGLNgR49enDCCSew5ZZbZlx/0UUXcfPNN1OnTp1KH2PZsmW0bduWdu3aUadOHb7//nvee+89vvnmG0455ZRK71dEKk4l1oQ98cQT1K1bl+7du5dr+549e7Ju3Tr23ntvbr31Vtq1awfA1ltvzaRJkwCoUaMGTZo04fjjj2fkyJFZi11EMlOJNUFLlizhyiuv5K677irX9o8++ijTp09nxowZdOnShaOOOopFixYBcP/99/PHP/6RPn368Mgjj3Dfffdx+OGHs2LFCo466ii6dOnC+PHjs3g2IlJEJdYEFRYWcuaZZ9KyZctybX/wwQcXv77iiisYOXIkEydO5Nhjj2WvvfZi3LhxAMydO5dLL72Ut99+m0MPPZQ777yT5s2b06lTJ2bMmIGZumwQySaVWBP02muvMWzYMLbddlu23XZbZs2axSmnnMLNN99crvebGe4bjlYzYMAAbrzxRurXr8+UKVNo3749LVu2ZPXq1SxYsCDu0xCRNCqx5sCaNWtYs2YNa9euZe3ataxYsYJatWrx2muvsXr16uLtOnTowO233063bt022MfMmTOZNWsWHTp0YN26dQwfPpyFCxeWKMUCvPLKK6xYsYJjjjkGgFatWvH666+z/fbbs3LlylJvoIlIfJRYc+DGG2/kuuuuK57/5z//ybXXXkthYWGJ7WrWrEmTJk0oKCgA4PzzzwdC/emSJUu44IIL+Prrr6lXrx577bUXL774YolEuXLlSi677DKeffbZ4mXDhw/n3HPPZeXKldx7773UrFkzi2cqIgCW6VIyn5mZp8bc9NaLEowmWHDZPUmHICK5V+rNCpVYc2Cbi/+edAgAzB9+btIhiFQLunklIhIzJVYRkZgpsYqIxEyJVUQkZkqsIiIxU2IVEYmZEquISMyUWEVEYqbEKiISMyVWEZGY5TyxmlknM3vOzGabmZtZ71zHICKSTUmUWAuAKUA/YFUCxxcRyaqcd8Li7i8AL0DoqSrXxxcRyTbVsYqIxCyvE6uZFUb1sMVT0jGJiGxMXidWdy90d0udko5JRGRj8jqxSu6NHj2a3XbbjQYNGrDzzjszceLEDbZZuXIlAwYMoHnz5jRp0oQLL7ywxNhd/fv3p0mTJhx44IF8++23xcsfe+wx+vXrl5PzEEmSEqsUe+WVVxg4cCAPPfQQS5YsYcKECey0004bbDd06FA++OADPv30U7788ks++ugjbrzxRgDee+89PvzwQ+bNm8chhxzC0KFDAVi8eDG33npr8XYim7Ik2rEWmNleZrYXYcyYHaL5HXIdi5R07bXXcs0113DAAQdQo0YNtttuO7bbbrsNtnv++efp168fW2yxBU2bNqVfv348+OCDAHzzzTcccsgh1K1bl9/85jdMmzYNgMGDB3PZZZex+eab5/ScRJKQRIm1PfBxNNUBroteX59ALBJZu3YtH3zwAQsWLGCXXXahRYsW9O3bl+XLl2fcPnVAR3fn22+/ZfHixbRp04aJEyeyfPlyXnvtNdq0acMHH3zAF198wemnn56r0xFJVM4Tq7uPS78hFU29cx2LrDd//nxWr17Nk08+ycSJE/nkk0/4+OOPM166d+3albvuuosFCxYwb948hg0bBsDPP/9M27ZtOemkkzjggAOYOXMml19+Of369WPYsGEMGzaMTp060bNnTxYtWpTjMxTJHdWxCgD169cH4OKLL6ZZs2ZstdVWXHLJJbzwwgsbbDt48GD23ntv9tprLw466CBOOOEEateuzTbbbAPAgAEDmDRpEv/61794/PHH6dSpE+vWrWPEiBG89tpr7LbbbsV1ryKbIiVWAaBJkya0aNECs/Ut2lJfp6pfvz533303s2fPZtq0aWy55Zbsu+++1KhR8s9p/vz5jBgxgmuuuYZPP/2UPffck9q1a9OhQwcmT55c6Vg7d+5MvXr1KCgooKCggNatW2fcrrCwkNq1axdvV1BQUFznu3jxYo466igaN25Mz549Wbt2bfH7+vTpw5gxYyodn4gSqxQ7++yzGT58ON999x0//vgjd9xxB8ccc8wG282ePZs5c+bg7rzzzjvccMMNXHfddRtsd8kll1BYWMhmm21Gq1ateP/991m6dCnjxo3L2NqgIu6++26WLl3K0qVL+eKLL0rd7ne/+13xdkuXLi0+7gMPPMDee+/N/PnzmT59Ok8//TQAb7/9NnPmzKFHjx6/KD6p3nLeV4Dkr6uvvpqFCxey6667Uq9ePU455RQGDx7MzJkz2X333Zk6dSo77LADX3/9NWeddRbfffcd22+/PUOHDuXII48ssa/XX3+dRYsWceKJJwKw3377cfTRR7P99tvTunVrnnzyySROsdg333zDCSecQN26denYsSPTpk1j7dq1DBgwgMceeyzR2KTqs9S7u1WBmXlqzE1vvSjBaIIFl91T5vptLv57jiIp2/zh5yYdQiw6d+7MZ599hrvTunVrbrrpJjp37rzBdoWFhdxxxx3UrFmTZs2a0bdvXy644AIA7rnnHqZPn87111/P4YcfzpVXXsmXX37JTz/9xLXXXpvjM5IqqtQnQVVilWKvnNMi6RA44sFvN7rNzTffzO67706dOnUYPXo0xx57LJ988gk777xzie1OOeUU+vTpwzbbbMO7777LSSedROPGjTnttNM499xz6devH/vvvz/du3enXbt2XH311YwdO5bzzz+fqVOn0qlTJz3QIJWiOlapcvbff38aNmxI3bp16dWrFwcffHDG1gu77747zZs3p2bNmhx00EH88Y9/LK6CqFevHiNGjGDy5MkMHTqUAQMGMGTIEB599FHWrVvH+PHjeffdd3nppZdyfXqyCVBilSrPzChPlVZp27300ku4O127dmXKlCm0b98eM6N9+/a/qPWCVF9KrFKlLFq0iJdffpkVK1awZs0aHn30USZMmEDXrl032PbZZ5/lxx9/xN157733GDZsGMcff3yJbVasWMGgQYO48847AWjVqhXjxo1j1apVvPnmm7+49YJUT0qsUqWsXr2aq666iqZNm7LVVlsxfPhwnnnmGXbddVcmTpxIQUFB8bajR49ml112oWHDhpx11lkMHDiQXr16ldjfkCFD6NmzJy1ahPrl8847j4ULF9K0aVNatGhR3KpBpCLUKiAGm0qrgKpy8+q9e5Lvr2e/i2YmHYIkr9RWASqxiojETIlVRCRmSqwiIjFTYhURiZkSq4hIzJRYRURipsQqIhIzJVYRkZgpsYqIxEyJVUTy3g8//MCJJ55IgwYN2HHHHUvtjHzs2LF06dKFRo0a0bJlyxLr1qxZw6mnnkrjxo3p2rUrP/30U/G6IUOGcPvtt8cWrxKrSBZUlURQ3jjdnYEDB7Lllluy5ZZbMnDgwOKewnIxfthFF11EnTp1mD9/Po8++igXXHABn3322QbbNWjQgHPOOYdbb711g3VjxozBzFi4cCGNGjVixIgRQBhN4rnnnqNfv36/OM4iSqwiWVBVEkF54xwxYgTPPPMMkyZNYvLkyTz//PM88MADQPbHD1u2bBlPPfUUN9xwAwUFBRxyyCEcd9xxPPLIIxtsu99++3HmmWdm7JXsm2++oXPnztSqVYsuXboUDyzZr18/brvtNmrViq/ffyVWkZhVlURQkThHjhzJpZdeSosWLdhuu+249NJLefjhh4vj7NKlS8bxw4YNG/aLYgT48ssvqVWrFrvuumvxsnbt2mX8AihL27Ztef3111m5ciVjx46lTZs2PP3002y11VYcfPDBvzjOVEqsIjGrKomgInF+9tlntGvXLuN2bdu25dVXX2X58uVMnDiRNm3aMGzYMLp16xZLf7ZLly5l8803L7GsUaNGLFmypEL76d69O61ataJDhw40atSIU089leuuu45bbrmFwYMH06lTJy688EJWrVr1i2NWYhWJWVVJBBWJc+nSpTRq1KjEdkuXLsXdOffcc1m8eDH7778/HTt2pF27djzyyCP079+f888/n06dOnHVVVdVKkaAgoKCEvXLAD/99BMNGzas0H7MjKFDhzJ58mRGjBjB0KFDOf/883n//ff54IMPGD9+PKtWreLBBx+sdKxFlFhFYlZVEkFF4kzf9qeffqKgoAAzy/r4Ybvuuitr1qzhq6++Kl42adIk2rRpU6n9AUyZMoW33nqLPn36MGXKFPbdd1/MjA4dOsQyHE8iidXM/mBm/zOzlWY2ycw2HFdDpIqqKomgInG2adOGSZMmbXS7bIwf1qBBA3r06ME111zDsmXLePPNN3n22Wc588wzN9h23bp1rFixgtWrV+PurFixYoMSvbvTt29fhg0bRo0aNWjVqhVvvPEGq1atYvz48bFUX+Q8sZrZScB9wG3AXsALwLNmtmeuYxHJhqqSCCoS51lnncXtt9/O7NmzmTNnDrfddhu9e/cusU02xw+79957Wb58OVtvvTWnnXYa9913H23atNlgOJ4JEyZQv359unfvzsyZM6lfvz5HHnlkiX099NBDtG3bln333ReAHj160Lx5c5o2bcr3339Pnz59Kh1nkSRKrJcBj7j7fe7+ubtfAUwG/phALCJZUVUSQXnjPO+88zj22GPZY489aNu2LUcffTTnnXdeiX1lc/ywLbbYgmeeeYZly5Yxc+ZMTj/9dAA6duzI0qVLi7fr3Lkz7l5iGjduXIl9nXPOOdxzz/rhlGrVqsXo0aNZvHgxL7/88gb1zpWR0zGvzKwOsAw4y91HpSy/ATjZ3X9djn1ozKtK0phX8dGYV7n11eTTkg6BX+05Kn1RqWNe5TqxNgdmA79x99dTll8EDHH3RmnbFwLX5ixAEZEKcPeMyTW+Rw2ywN0LgcJsHiMqAZf6zZMvFGd8qkKMoDjjlss4c13HuhBYA2yTtnwbYG6OYxERyYqcJlZ3XwV8CByRtqor8GYuYxERyZYkqgL+Aowys/eAccBZQDvg9wnEIiISu5wnVnd/0sy2AAYCdwFfACe4+y9/3KFyrkvouBWlOONTFWIExRm3nMWZ01YBIiLVgfoKEBGJmRKriEjMlFhFRGKmxCoiEjMlVhGRmFXrxGpmdcyshZntkDolHVdVYGa/MrNRZtYow7rG0bpf3rFlDMxsczPb4G/dzGqa2S/vyigmZtbLzD4wsx/MrGW07FIzq3y3UNWUmd1lZhuMtmhm/czsjmwfv1omVjPbxczGAz8DM4Bvoml69DMxZvZ6eaYkY4xcCsxx98XpK9x9EfAt8KdcB5XOzE4gPO1XL8PqesCHZnZMToPKwMzOI/RR/AxQn/X/m9+TZ11qmtkgM9ugqzQzO9fMLk8ipgxOAt7OsPxN4LfZPnhed8KSRX8n/FOdAswB8qkxb2dCsv8PsC7ZUMp0KNC7jPWPAxsO95l7FwB/cfef01e4+zIzuwXoC/w755GV1A/o4+5jzGxgyvIPgQ3Hxk5WH2DD3rBhKvAocEtuw8loK+CHDMsXAU2zffDqmlj3BfZz96lJB5JBISFhdQdGAg+5+9dJBlSKHSm745z5QD5Uq7QBLi5j/Tiy3INaOe1ESKLplgMFGZYnqaj7z3TzonX5YDrhyz/9f6czMCvbB6+WVQHANPLvjxUAd7/e3XciJNedgClmNtbMzjCzTJezSVlKiK80OxE6NU/alhtZb8AWuQhkI2YBu2dY3gX4KsPyJM0lFE7StQe+y3EspRkB3GZmF5hZ62i6kNBXyf3ZPnh1LbFeDPzZzC7O01Ir7v4a8JqZNQHOAC4HhpnZju5esXGUs+Nt4GxCiS+Tc8lcx5Vr3xI6+fmylPV7Rdsk7X7gLjNbHs23MrPDgJuBK5MLK6PHCLH+DIyNlh0G3AH8I7GoUrj77WbWFLgdqBMtXgXc7u635SKAajcRLq/WAGuBlYSbWMVT0vGlxbob4abGfELJpUHSMUVxdYw+w2FAs5TlzYDh0bpOeRDnHcB/M31uQINo3R1JxxnFcz2hlL8umpYD1yUdV4Y4awOjohjXpkyPArWSji/D77hDNOXsf6dadsJiZr0p44aVu4/MXTQbMrMC4FRCqW9P4EngQXcfn2Rc6czs94QkWgcoGnR+c2A10M/dRyQVWxEz2wr4mFDtNRz4PFq1O+Gm1TpgH3dfkEyEJZlZfUK9cA3gM3fPh+qUjMzsV4QSP8DH7v6/BMPJK9UyseYzM3uI0BzkM0LrhdGeH5f+GUXjmJ0C/IpQX/kF8KS7Z7q5kQgz254w5Ho31g8A58CLwEXuPiOp2EoTtbvdA5jhofmabISZvQCc5u6Lo9elcvfuWY2lOiZWM3uWUDc01t0nJR1PKjNbB8wkNF0pVbb/MDZFUX31LoTk+pW7/5hwSMXM7E7gU3f/W5RUxxKqW5YBR7v7hITjuxe43N2XRq9L5e4X5iisEqJCST93XxK9LpW7n53NWKrrzatvCHfdbzOzRcAEwk2YsZ5ch9tF/kF+tavNyMz2K8927v5etmMpryiRvp90HKXoQbgpBKGpXWvCZfaZwE2EJJuk3VifL3YrY7vE/nZTk2W2E+fGVMsSaxEz2xLoRGjS0oXwB/Oju2e9AXFVF5WsnTLGVgfc3WvmKKSMzOzBcmzm7r7Bk0S5ZGYrgF3c/VszuwdY5+4Xm9nOwEeeNjS85LfqWmItsoLQEmB5NK0jP5rebCDqw6Ap8Im7r006HqBV0gGUU7My1hmhJFiPcKMwSQuBloS/vyMJzetgfVMhqYDoJuAlwG+ArUlrs+/umdoMx6ZaJlYzG0J4AmMfws2WccCfgfHunukxuJwxs98BW7j7fSnL7iM8RgjwjZl1cfesPz2yEVXiUsfdu2Vabma/IbQRXQfcmNOgMnsKeNTMvgQaEx5phlAdkHd3282sC3A4Yej69KR1TiJBlfQ34ChCi5oJ5PjvtVpWBUSXsQsIT2E8nk93hc3sDeARd38gmj+c8E92NaGp0E3AOHe/ILkoS1QFbLAqZbm7e159eZvZnoSEehih1cV17j4/2ajAzGoROlvZAXjY3T+Olg8Alrj735KML5WZXUn4MvqUDH1tlPZllktmthg4NqmbftU1se7P+nrVgwiP4Y1lfUuBOQnGtgA4vKi1gpkNB37l7l2j+e7A3R4ee01M9BmW5nhCkljp7vnwuGhRk6sbgZ7Ac8Agdy/taSwpg5l9C1yfD+2US2NmXwHHe0JPVlbLxJoqKinsR6hjOwOomWQpK3qksbW7z4zmPyK0Zb0lmt8R+NzdN0sqxtJEyfZWwjPj9wBDkm7SZGaNCY+E9iV0cjLQ3d9KMqbSRM2sdibz5XWiza1SRS1p9nH3aUnHUhozO53Qvvosd/9pY9vHLa8u03LJzLYjlFgPi37uSOixZ1yCYUG4edEWmBl1Ir0HJXtn2pLQAUreMLNdgKGEkuoo4IyiL4Y88DXh5tSdhL5OMzYVS7pZmJl1AEYTbmClt7RwINHWFWkeA44F7ko6kDI8A5wFfGdmswhPAxbTzassiC4TdiJ0czYOGEKoAsiHmwRPEDq42A7oSqjDeidlfXvC8+2JM7OtCV3u/Z5QjdI+3x64AJpEPwcBA8ncPCwfEtf9hKftTib/+ghONw+4xswOAT5hw6SVD/2xPkz4X3mYEK9uXmVb1Fv72HysY4uaiTxAKBHMI3R+PDFl/VjgJXe/OaEQi+K4ljCKwJeEy+vXkoynNFHVyUYlfQPTzJYBe7l7vnURuAEzK2uUDU+6/h/AzJYC3XXzKiFmVhvA3VdvbFtZL2oVsJyNNGVJ+tFbK+cYZklXXZjZ28CV7j52oxvLRpnZ14RWAbp5lUtRD1dXsL6z5q+BP+dBz1Zldh4RcXc/OuvBlMHMHqYcl1eJP1pYRZqFRc3qhgLXAJPY8PI6XzqQLiG6D/CT51kiMbPTgN8BvTzDuGxZP36efR45EbUNvIlw5/qNaHFHwvhIV7p7YpXyG+s8okjSCauqqCrNwqIvgCKp/5RGHjwanMrMagJXEW6qNgZ2dfdpZnYTMN3d/5pkfABm9jnQgtB3rG5e5cjFQF93T32O/Nnol3ElCd7tVMKMl7u/m74sU7OwXMeVQZekA6iAK4BeQH/CEChFJgMDgMQTK6GFRWKqa4l1JdAmvRVA1GzoU3fPp7GlJCYZmoVdlXTdalUUPXZ7kbu/YmZLgHZRiXU34M2kS//5oLoOJjgdyFRH2T1aJ5sQM9s66kN0KtCQ0CzsrHxLqmb2KzO71cyeNrNtomVHm9keSceWZnsyjyG2Dqib41jKZGYHmtm50VRWtVCsqmtVwC3A/Wa2D+vrWA8GTgPOTywqiV1as7BuedwsrBPwEvAWob6/QbSqLXAOcFJCoWXyFeFpxfQmascS+g9IXDSQ4BOEbkGLbl41MrNxwCnuvjCrx6+OVQEAZnYccBnrhxyeCtzs7v9OLiqJWxVqFvYm8Iy735p2ed0BeNrdWyQZXyozO5NQRz0IuBu4iDA0zyXAme7+VILhAWBmo4BfE54C/Cxa1pbQkfzn7t4zq8evrolVqocq1CxsCbCnu3+TllhbAv/Nt3p/MzuH0DSsqJ3wbOAady9Xq5Zsi/oz6Oru76QtPxB4wd2bZHxjTKprVYBUE+7eO+kYymkpsBVh2KBUbQhP4OWVqEXNgxZGwa2Rh+1sa5O5T41l0bqsqpY3r8ysqZmNNLPZZrbGzNamTknHJ9XS08D10SPNAG5heOmbCXWFecndF+ZhUoVw7+TPZtawaEH0+gZgYqnvikm1rAows+cJ41vdS+aOev+VRFxSfUX/9P8G9gXqE4Zq2QoYTxildXmC4ZUQtffOlDicMNzRV8Df3P3VnAaWImr69TLQiPU31NoQbmQd5e5Z7ciouibWn4DD3P2DpGMRSRUNebIv4Wryw3xsxWBhaKMLCb1xFT2AsR+hBcPD0c8uwInu/lwSMUJxh0Y9WT+q7FTgUXdfkfVjV9PE+hXQw92nJB2LSFVjZvcD37n7NWnLrwW2dfcLosdbj3T3DokEmbDqmlhPJozXfmYSHTSIpDOzXsBCd/+/aP7PhDbVU4HT8ulhBjP7EeiQ4cnFXwHvuXuT6FL8XXffPJEgKe4gZj8yj8jwj2weu7q2Crie0EHD/CR6FxfJ4ArC8DFFown0J3TMfTRwG6ED7Hyxjsyjx+4VrYNQ37omdyGVZGbdCI8tbw6somSdsBPas2ZNdU2sZXXQUK7+O0Vitj3hpg/ACcAYdx8WPSmU2E2gUjwM/NXMWrO+jnV/whNuRR0bHQokWdV2O/A4YdDInA9pX62qAszsTnfvX8b65sB4d/9V7qISATP7HjjU3T+NOr3+q7s/aGatCB0DNdjILnImGvTwcqAfsG20eB6hV7i/uPvaaGihde4+N6EYlxEeuPg6ieNXtxLrOWY2z92Hpq+Ixm96Dfg+92GJMAG4zczeAPYBXoyWtyYMMJkXLIxqfDzwoLsPNbPNAdJHQnX32UnEl+INQvMqJdYc6AE8b2bzUx+9i54eeR34GTgqqeCkWrsYuA/4LXBBSkmvK6E9Zl5w9zVm9k/Cc/gbJNQ88jfgL2bWgswjMmR1VN5qVRUAYGanECquT3H358xsC0JSNaCzu/+YaIAiec7M3gGudfe8Sfjp0kZkSJf1ERmqW4kVd388KqGOMrNTCcM310ZJVaS8biCUBhsCHxKevy+WJ4+4tkry4NWuxFrEzK4BriXciT3U3ecnHJJUM2b2M7Cjuy8ws+WU3a3hZrmLrGxVaXyupFSrEmuGEVBXE75tHzKz4oVJ980p1cYFwE8pr6tKKadKjM9lZocS2gbvDBzn7t9GozNPc/cJ2Tx2tUqsQHqpdFQiUYgAqUOtu/vDCYZSIe4+PukYNsbMegCPRtPuQJ1oVQNCB91ZTazVtipAJJ+Y2TTCY6Lfpy1vDHzk7jslElgposdX+wC7AOe7+3wzOxqYmQ99cJjZR8Bwd38orePwvQkdXTfL5vGrZX+sInmoJZCpbrIu4fHrvBGNzzUJ2JswAGfq+FyFCYWVrjUwNsPyRUBWRw+A6lcVIJJXoiRV5MCog5MiNYEjgbzpgCXyZ0Jzq6LxuYq8TmiPmw8WEh5Pn562fD9y8HkqsYokaxzhppUTRhFIt5jQ92k+2RM4I8PyBYTOufPBPwlNwk4hfLZ1zexwwiCI92f74EqsIslqRmimNIeQsBakrFuVp22rq8L4XNcSOrb5mvD5fhr9fATY4JH2uOnmlYhUiJndS2iA3wP4jvCFUItQ4n7R3S9LMLwSok5sikZk+Ci9D9msHVeJVSQ/JNkxc0Xk8/hcZtYEOMPdh0fzT0UxFlkD9Mr2lYASq0geiDpmHg00JEPHzPn05FWRfByfy8wuIzStOiOaX0LoxKboJtuBwGh3L8xqHEqsIsmLRj6dSEIdM1eEmZ0F/MvdV6YtrwOcmmTpOuog5mp3fyWaL27DGs33AAa7+75ZjUOJVSR5SXfMXBFmthZolt7ZipltSRhkMLG+AszsO6B90RhhZvZfwqCGRfM7A5+4e8NsxqEHBETyQ1HHzFWBkblfg21Yf8mdlIZAQdGMu/86bSDGemR+ECNWam4lkh8S7Zi5PMzs9eilA0+b2aqU1TUJnV9n9Rn8cpgB7EEY3TaTPcnBAwKqChDJA0l3zFweZnZf9PI8QmfxqXf/VxHatT6U5JDyZvYX4Bhg7/TWCWbWgNB/7L/d/U9ZjUOJVSR5ZrZjKavqAie6+825jKcsZnYtYdDAZRvdOMeiTuw/jmaHA19Er39NNLw4IekuzGocSqwi+cfMdgd+D5wJNHL3Oht5i0TMbHvCY6tdCfXBEKovXgIudPcZWY9BiVUkP5hZfeBU4A/A/oT6ylHAmGyXsCrCzDYDBgKHk/lhhrzo4jB6WKBoKPv/5bIZm25eiSTMzNoTkumpwCxCByL7ARe5e2k3YZJ0P6G7wMcIfRzkZekseroqkZt+KrGKJMjMPiZ0xDIaeMTdP4yWryY0bM+7xBp1bdjD3TP1dyqoHatI0vYAngEeK0qqVcBSYG7SQeQzJVaRZO1OGFDweTP7r5ldWUYLgXxxM3C5mSXeBCxfqSpAJA+YWW3gBEJLgMMIhZ7BwF/Tx8FKmpm9SOjMZBmhIX76wwzVfpRjJVaRPGNmLYFzgd6Eu+7j3f2IJGNKZWYPlbXe3c/OVSz5SolVJE+ZWQ3gaOD37n580vFI+Smxiki5mNkL5djM3f3orAeT59SOVUTKa37SAVQVKrGKiMRMza1ERGKmxCoiEjMl1mrAzLwcU2HScf4SZjbdzO6uxPv6m1mVbXdpZi2j399vU5ZV6XPaFOjmVfVwYNr824S+Kh9LWfZt7sLJihOBygxp3J8wlHN57nhXFf3Z9M6pSlFirQbc/Z3UeTMDmJm+PG2b+kmOD19eRXG6+8cb31okN1QVIJhZ7+hycn8zezEaMfTuaF1/M3vfzBab2QIze8nM2qS9/2Ez+9TMOpnZR2b2s5l9bGad0rY7xszeM7Ol0f4+ioYjTt3mrOi9K8xsoZm9UPTs/EbiLFEVkBLTUWY2Jdrfh2Z2YMo204EdgYtSqkR6l/E5uZldZmbXmNlcM/vBzEZH/X6mbreDmT1uZouiz2KsmXVIWf9QNHpo+v67FJ1fyrKuZvZWtJ8fzOxBM2tcRowZz8nM7jKzmdFDB6nbHxRtc1Bp+5SKU2KVVKOBN4Fjgb9Hy7YH7iE8x342sBJ4y8yapb1322i7O4CTCGMgPWNmDaF42OExwOeEy/aTCVURxUnJzC4DRgIfRfs4F/gKaFqOODNpRug79C/A7wjPtL9sZltH608E5gFPEqpLDgT+r4z9QRjeow1wDnAFYXylW1POoSEwntCfal/gNMKV4Tgz+3W02SigtZntnbbv04Cv3f3daF8nRPF8Qfg8/kjoXPpfZcRX2jmNIPwu0x+NPRf43N3f2sh5S0W4u6ZqNhE6Jv5TynzvaNmVG3lfTaAO8B1wScryh4F1wB4py/aK9nlCNP/baL5hKftuROjU44Eyjl9qnMB04O60mBw4LGVZE0KXd38u7X3l+NzeT1t2J7AoZb5f9Fm0TVnWAFgAPJzyOc4DbknZpjbwPXBDNG+Ewfn+lXa8A6I4DonmW0bzv93YORG+jB5PmS8gDFd9adJ/k5vapBKrpHoufYGZHWBm/zGzhcAaQom1KbBr2qZz3X1KynxRB80top+TgbXAY2Z2XIbL2QOBzSi7BFpqnKVY7O5FQzbjoUf51wnDnlTWf9LmpwKNzKxoLPuOwGfu/mnKcZcRbiZ1jObXAo8Dp1pU4Q0cBWzB+huKvyIkzdFmVqtoAj4gdDO4XyViHwEcb2ZbRvOnEgYr/Ecl9iVlUGKVVCUeWTSzHQiJpA5wAXAI0IEwdnu9tPeWuCPv7kVjzteL5r8kXDZvDjwFLDCz/zOzovGRiv7Z51Q0zjIsKOW96dUYFZHe8qDEeRJKxfNKOe4WKfOPES7ND47mTwM+cffPo/mi6o8xhCqM1GlzYIdKxP448DNhgEIIXRQ+5+6ZPif5BdQqQFKlP9/clXC5eGJU2gPAzLagEtz9JeClqHR3JHAbIcEcQLgMBmjOxpt+lfc57PS6WQjd8GWz9/sfCEMtZzpu8WB27v6OmX0DnGZmHwHHA9en7QdCPe27Gfb3XUUDc/flZvZP4Fwze4VQcr+uovuRjVOJVcpSn5DEijsyNrPjgIa/ZKfuvtTdxwCPEHrQh9C29mfCDbK4NDKzw4pmorv3h1EyUa1iw9L3L/EG0NbC8NVFx92M0P3fxLRtRxFu4p1IqAYZlbLuv4SBBXdx9w8yTDPLiKGscxoBtAXujfb/cvlPTcpLJVYpS1H95MNmdj+hXvUKMl9il8nMziNc9r5IuNzfgXBH+hUAd19sZtcBN1sY8uMZwhd/F2CUu39Qifh/AP5u4amyH4ErCV8Ud6Zs8zlwuJkdGW3/jf+yHvsfAgYA/zazqwg3y/5EuIE1NG3bx6KYbgHecPdZRSvc3c2sP/AvM2tAqKNdQvjcjgLu8qj1QAalnpO7TzGzd4BOhBtl637BuUopVGKVUkU3o3oR7vA/T7grfxqVuAwl3LxqQmj69B/gJuBpQrOlouPdEs0fEK17mJDMK3M8CJf8FwKXA08Q6oq7untqHe2VwExC86T3CU24Ks3dlwCHAh8C9xGahq0Durj7f9O2/YzwuTSn5FNwRevHEJLozsCjhOQ6iFBtMr2MMDZ2TmMIXzAPVujkpNzUbaBskszsYaC9u7dNOpZ8Y2avA2vc/cikY9lUqSpApJows/aElh1dAHXSkkVKrCLVx/uENrB/dvcXkw5mU6aqABGRmOnmlYhIzJRYRURipsQqIhIzJVYRkZgpsYqIxOz/Ab2JkzxzO2nVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_dict, order = get_talon_nov_colors()\n",
    "opref = 'test'\n",
    "plot_read_novelty(df, opref, c_dict, order,\n",
    "                      ylim=None, title=None, \n",
    "                      datasets='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e14d42b9-e4a0-419e-92a5-8f014ef89a3c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16888524"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post minimap reads\n",
    "3633558+4507243+3836381+4911342"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b634e6ed-d553-4176-a292-b1951eca3d62",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10267664"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post tc reads\n",
    "2242821+2700669+2299503+3024671"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "33137a5e-4b75-472c-a5b0-2b26e11a73fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9729920"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# post talon reads\n",
    "len(df.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f9a0bf5b-5368-40e7-9de5-b8e31ee51e6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.375"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "6/16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71df0c98-525e-4820-a311-31cf2049d3d9",
   "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}