Generate-common-TSS-TTS-HCT116_HAP1_step2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "aerial-climate",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
       "    function(sublime_keymap, cell, IPython) {\n",
       "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
       "        var cells = IPython.notebook.get_cells();\n",
       "        for(var cl=0; cl< cells.length ; cl++){\n",
       "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
       "        }\n",
       "    }\n",
       ");\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
    "    function(sublime_keymap, cell, IPython) {\n",
    "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
    "        var cells = IPython.notebook.get_cells();\n",
    "        for(var cl=0; cl< cells.length ; cl++){\n",
    "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
    "        }\n",
    "    }\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "corresponding-poultry",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:85% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# change the cell width\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:85% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "suited-vocabulary",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "%config InlineBackend.print_figure_kwargs={'facecolor' : \"w\"}\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from more_itertools import chunked\n",
    "# import all relevant func and modules from a companion .py file\n",
    "from supp_lib import *\n",
    "# import all samples as dictionary ...\n",
    "from samples import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "romance-booth",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's start with the super useful intro on how to extract a BED file with gene annotation from HiGlass docs ...\n",
    "# https://docs.higlass.io/data_preparation.html#gene-annotation-tracks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wired-hotel",
   "metadata": {},
   "source": [
    "### Further filter and preprare our RefSeq-based list of genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "tamil-train",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_df = pd.read_csv(\n",
    "    \"./gene_annot_data/hg19/geneAnnotationsExonUnions.bed\",\n",
    "    sep=\"\\t\",\n",
    "    header=None,\n",
    "    names = [\"chr\",\n",
    "            \"txStart\",\n",
    "            \"txEnd\",\n",
    "            \"geneName\",\n",
    "            \"citationCount\",\n",
    "            \"strand\",\n",
    "            \"union_geneId\", # \"refseqId\",\n",
    "            \"geneId\",\n",
    "            \"geneType\",\n",
    "            \"geneDesc\",\n",
    "            \"cdsStart\",\n",
    "            \"cdsEnd\",\n",
    "            \"exonStarts\",\n",
    "            \"exonEnds\"]\n",
    ")\n",
    "genes_df.head()\n",
    "\n",
    "### we'll be removing duplicated genes because they're likely unamppable and/or harder to deal with anyways ...\n",
    "### there are of course expections, but will do it anyways:\n",
    "geneId_grp = genes_df.groupby(\"geneId\")\n",
    "genes_df = geneId_grp.first()[geneId_grp.size() == 1].reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "noticed-comparative",
   "metadata": {},
   "source": [
    "### some stats ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "resident-oakland",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='geneType,strand'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAFdCAYAAADhUAFNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5bUlEQVR4nO3deVhV5fo//vdOckJFSEwCFWxzFEEExCkVRcM5NDXHAhG0nIfjSf00qA2GaeUQnSPlgKaikoojkpYplsMWo8xSNFBAQlQcUJnv3x982T9RnGCtBWzfr+vyuthrw3OvB2G/Wet+1to6EREQEdFT7Zny3gEiIip/DAMiImIYEBERw4CIiMAwICIiMAyIiAiPEQajRo1C/fr14eLiYtx29epV+Pj4wNHRET4+PsjIyAAAiAgmTZoEvV4PV1dXxMbGGr8mLCwMjo6OcHR0RFhYmHH78ePH0aJFC+j1ekyaNAlc6UpEpD3do64zOHDgAGrVqgU/Pz+cPHkSAPD222/DysoKM2fORHBwMDIyMjB//nzs2rULS5cuxa5du3DkyBFMnjwZR44cwdWrV+Hp6QmDwQCdTodWrVrh+PHjsLS0RJs2bbBkyRK0bdsWvXv3xqRJk9CrV69H7ni9evVgb2+vyDeBiOhpkJiYiMuXL5f4nNmjvtjLywuJiYnFtkVGRmL//v0AAH9/f3Tp0gXz589HZGQk/Pz8oNPp0K5dO1y7dg2pqanYv38/fHx8YGVlBQDw8fFBVFQUunTpghs3bqBdu3YAAD8/P2zduvWxwsDe3h4Gg+GRn0dERIU8PT0f+Nwjw6AkaWlpsLGxAQA0aNAAaWlpAICUlBQ0bNjQ+Hl2dnZISUl56HY7O7v7tj9IaGgoQkNDAQDp6eml2XUiIipBmRvIOp0OOp1OiX15pDFjxsBgMMBgMMDa2lqTmkRET4NShcHzzz+P1NRUAEBqairq168PALC1tUVSUpLx85KTk2Fra/vQ7cnJyfdtJyIibZUqDHx9fY0rgsLCwtCvXz/j9tWrV0NEcPjwYVhYWMDGxgY9evRAdHQ0MjIykJGRgejoaPTo0QM2NjaoU6cODh8+DBHB6tWrjWMREZF2HtkzGDZsGPbv34/Lly/Dzs4Oc+fOxcyZMzF48GAsX74cjRs3xsaNGwEAvXv3xq5du6DX61GzZk2sXLkSAGBlZYX33nsPrVu3BgC8//77xmbyV199hZEjR+LOnTvo1avXYzWPiYhIWY9cWlpRFS1VJSKix/Ow101egUxERAwDIiJiGBAREUp50Rlpx37mzif6/MTgPirtCRGZMh4ZEBERw4CIiBgGREQEhgEREYENZMKTN6kBNqqJTA2PDIiIiGFAREQMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIilDEMvvjiCzg7O8PFxQXDhg1DVlYWEhIS0LZtW+j1egwZMgQ5OTkAgOzsbAwZMgR6vR5t27ZFYmKicZxPPvkEer0eTZs2xZ49e8o0ISIienKlDoOUlBQsWbIEBoMBJ0+eRH5+PsLDwzFjxgxMnToVZ8+ehaWlJZYvXw4AWL58OSwtLXH27FlMnToVM2bMAACcOnUK4eHh+OOPPxAVFYVx48YhPz9fmdkREdFjKdORQV5eHu7cuYO8vDzcvn0bNjY2+OGHHzBo0CAAgL+/P7Zu3QoAiIyMhL+/PwBg0KBB2LdvH0QEkZGRGDp0KKpVqwYHBwfo9XocPXq0bLMiIqInUuowsLW1xfTp09GoUSPY2NjAwsICrVq1Qt26dWFmZgYAsLOzQ0pKCoDCI4mGDRsCAMzMzGBhYYErV64U237v1xARkTZKHQYZGRmIjIxEQkICLl68iFu3biEqKkrJfbtPaGgoPD094enpifT0dFVrERE9TUodBnv37oWDgwOsra3x7LPPYsCAATh06BCuXbuGvLw8AEBycjJsbW0BFB5JJCUlASg8vXT9+nU899xzxbbf+zX3GjNmDAwGAwwGA6ytrUu760REdI9Sh0GjRo1w+PBh3L59GyKCffv2oXnz5vD29kZERAQAICwsDP369QMA+Pr6IiwsDAAQERGBrl27QqfTwdfXF+Hh4cjOzkZCQgLi4+PRpk0bBaZGRESPy6y0X9i2bVsMGjQIHh4eMDMzg7u7O8aMGYM+ffpg6NChePfdd+Hu7o7AwEAAQGBgIN544w3o9XpYWVkhPDwcAODs7IzBgwejefPmMDMzQ0hICKpUqaLM7IiI6LHoRETKeydKw9PTEwaDobx3Q3X2M3c+0ecnBvdRvUZp6xBR+XrY6yavQCYiIoYBERExDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIpQxDK5du4ZBgwahWbNmcHJywi+//IKrV6/Cx8cHjo6O8PHxQUZGBgBARDBp0iTo9Xq4uroiNjbWOE5YWBgcHR3h6OiIsLCwss2IiIieWJnCYPLkyejZsyf++usvxMXFwcnJCcHBwejWrRvi4+PRrVs3BAcHAwB2796N+Ph4xMfHIzQ0FGPHjgUAXL16FXPnzsWRI0dw9OhRzJ071xggRESkjVKHwfXr13HgwAEEBgYCAKpWrYq6desiMjIS/v7+AAB/f39s3boVABAZGQk/Pz/odDq0a9cO165dQ2pqKvbs2QMfHx9YWVnB0tISPj4+iIqKKvvMiIjosZU6DBISEmBtbY2AgAC4u7sjKCgIt27dQlpaGmxsbAAADRo0QFpaGgAgJSUFDRs2NH69nZ0dUlJSHridiIi0U+owyMvLQ2xsLMaOHYsTJ07A3NzceEqoiE6ng06nK/NOFgkNDYWnpyc8PT2Rnp6u2LhERE+7UoeBnZ0d7Ozs0LZtWwDAoEGDEBsbi+effx6pqakAgNTUVNSvXx8AYGtri6SkJOPXJycnw9bW9oHbSzJmzBgYDAYYDAZYW1uXdteJiOgepQ6DBg0aoGHDhjh9+jQAYN++fWjevDl8fX2NK4LCwsLQr18/AICvry9Wr14NEcHhw4dhYWEBGxsb9OjRA9HR0cjIyEBGRgaio6PRo0cPBaZGRESPy6wsX7x06VKMGDECOTk5aNKkCVauXImCggIMHjwYy5cvR+PGjbFx40YAQO/evbFr1y7o9XrUrFkTK1euBABYWVnhvffeQ+vWrQEA77//PqysrMo4LSIiehJlCgM3NzcYDIb7tu/bt+++bTqdDiEhISWOM2rUKIwaNaosu0JERGXAK5CJiIhhQEREDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiggJhkJ+fD3d3d/Tt2xcAkJCQgLZt20Kv12PIkCHIyckBAGRnZ2PIkCHQ6/Vo27YtEhMTjWN88skn0Ov1aNq0Kfbs2VPWXSIioidU5jBYvHgxnJycjI9nzJiBqVOn4uzZs7C0tMTy5csBAMuXL4elpSXOnj2LqVOnYsaMGQCAU6dOITw8HH/88QeioqIwbtw45Ofnl3W3iIjoCZQpDJKTk7Fz504EBQUBAEQEP/zwAwYNGgQA8Pf3x9atWwEAkZGR8Pf3BwAMGjQI+/btg4ggMjISQ4cORbVq1eDg4AC9Xo+jR4+WZbeIiOgJlSkMpkyZgk8//RTPPFM4zJUrV1C3bl2YmZkBAOzs7JCSkgIASElJQcOGDQEAZmZmsLCwwJUrV4ptv/dr7hUaGgpPT094enoiPT29LLtORER3KXUY7NixA/Xr10erVq2U3J+HGjNmDAwGAwwGA6ytrTWrS0Rk6sxK+4WHDh3Ctm3bsGvXLmRlZeHGjRuYPHkyrl27hry8PJiZmSE5ORm2trYAAFtbWyQlJcHOzg55eXm4fv06nnvuOeP2Ind/DRERaaPURwaffPIJkpOTkZiYiPDwcHTt2hVr166Ft7c3IiIiAABhYWHo168fAMDX1xdhYWEAgIiICHTt2hU6nQ6+vr4IDw9HdnY2EhISEB8fjzZt2igwNSIielylPjJ4kPnz52Po0KF499134e7ujsDAQABAYGAg3njjDej1elhZWSE8PBwA4OzsjMGDB6N58+YwMzNDSEgIqlSpovRuERHRQ+hERMp7J0rD09MTBoOhvHdDdfYzdz7R5ycG91G9RmnrEFH5etjrJq9AJiIihgERETEMiIgIKjSQyxvPfxMRPTkeGRAREcOAiIhM8DSRVng6iohMCY8MiIiIYUBERAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiICw4CIiMAwICIiMAyIiAgMAyIiAsOAiIjAMCAiIjAMiIgIDAMiIgLDgIiIwDAgIiIwDIiICAwDIiJCGcIgKSkJ3t7eaN68OZydnbF48WIAwNWrV+Hj4wNHR0f4+PggIyMDACAimDRpEvR6PVxdXREbG2scKywsDI6OjnB0dERYWFgZp0RERE+q1GFgZmaGzz77DKdOncLhw4cREhKCU6dOITg4GN26dUN8fDy6deuG4OBgAMDu3bsRHx+P+Ph4hIaGYuzYsQAKw2Pu3Lk4cuQIjh49irlz5xoDhIiItFHqMLCxsYGHhwcAoHbt2nByckJKSgoiIyPh7+8PAPD398fWrVsBAJGRkfDz84NOp0O7du1w7do1pKamYs+ePfDx8YGVlRUsLS3h4+ODqKioss+MiIgem5kSgyQmJuLEiRNo27Yt0tLSYGNjAwBo0KAB0tLSAAApKSlo2LCh8Wvs7OyQkpLywO0lCQ0NRWhoKAAgPT1diV0nIiIo0EDOzMzEwIEDsWjRItSpU6fYczqdDjqdrqwljMaMGQODwQCDwQBra2vFxiUietqVKQxyc3MxcOBAjBgxAgMGDAAAPP/880hNTQUApKamon79+gAAW1tbJCUlGb82OTkZtra2D9xORETaKXUYiAgCAwPh5OSEadOmGbf7+voaVwSFhYWhX79+xu2rV6+GiODw4cOwsLCAjY0NevTogejoaGRkZCAjIwPR0dHo0aNHGadFRERPotQ9g0OHDmHNmjVo0aIF3NzcAADz5s3DzJkzMXjwYCxfvhyNGzfGxo0bAQC9e/fGrl27oNfrUbNmTaxcuRIAYGVlhffeew+tW7cGALz//vuwsrIq47SIiOhJlDoMOnbsCBEp8bl9+/bdt02n0yEkJKTEzx81ahRGjRpV2l0hIqIy4hXIRETEMCAiIoYBERGBYUBERGAYEBERGAZERASGARERgWFARERgGBARERgGREQEhgEREYFhQEREYBgQEREYBkREBIYBERGBYUBERGAYEBERGAZERASGARERgWFARERgGBARERgGREQEhgEREYFhQEREYBgQEREYBkREBMCsvHeAnh72M3c+8dckBvepsHWITAmPDIiIiGFAREQMAyIiAsOAiIjABjIRUamY2kIFhgERUQWmVejwNBEREVWcMIiKikLTpk2h1+sRHBxc3rtDRPRUqRBhkJ+fj/Hjx2P37t04deoU1q9fj1OnTpX3bhERPTUqRBgcPXoUer0eTZo0QdWqVTF06FBERkaW924RET01dCIi5b0TERERiIqKwjfffAMAWLNmDY4cOYIvv/yy2OeFhoYiNDQUAPDXX3+hWbNmj10jPT0d1tbWyu0061SqGqxTcWuwjnY1EhMTcfny5ZKflApg06ZNEhgYaHy8evVqGT9+vKI1WrVqpeh4rFO5arBOxa3BOhWjRoU4TWRra4ukpCTj4+TkZNja2pbjHhERPV0qRBi0bt0a8fHxSEhIQE5ODsLDw+Hr61veu0VE9NSoEBedmZmZ4csvv0SPHj2Qn5+PUaNGwdnZWdEaY8aMUXQ81qlcNVin4tZgnYpRo0I0kImIqHxViNNERERUvhgGRETEMCAiogrSQFZDVlYWduzYgYMHD+LixYuoUaMGXFxc0KdPH0Wb07/88gu+/fZbHDx4EKmpqcXqvP7667CwsFCsFgDcunUL1atXR5UqVRQdt0hBQQHi4uKKfc/q16+v2Phaf78yMjKMc7G3t8czz6jz948WdQwGw30/zz4+PrC0tFSsRnJyMsLDw0v8venVq5ei87p06RIOHTpUrI6np6diNbR6DQDUnwug/nxMsoE8e/Zs7NixA126dEGrVq1Qv359ZGVl4cyZM/jxxx+RlZWFzz77DK6urmWq06tXL7zwwgvo168fPD0976uzfft2TJs2rUzLZAsKChAeHo61a9fi2LFjqFatGrKzs1GvXj306dMHb775JvR6fZnmAQDnzp3D/PnzsXfvXjg6OsLa2to4l5o1a+LNN9+Ev79/mX64tfh+AcD169cREhKC9evXIycnxziXtLQ0tGvXDuPGjYO3t3eZamhZZ+XKlVi6dCkcHBzu+3k+dOgQXFxc8OGHH6JRo0ZlqhMQEICUlBT07du3xP+f48ePIzg4GF5eXmWq8+OPPyI4OBhXr16Fu7t7sTrnzp3DoEGD8O9//xt16tQpdQ2tXgO0mItm81H0ErYKYseOHQ99Pi0tTY4dO1bmOunp6Yp8zsN4eXnJBx98IHFxcZKfn2/cfuXKFYmIiJABAwbImjVrylRDRGTo0KHy008/SUFBwX3PpaWlyRdffCGrVq0qUw0tvl8iIi+//LKsXr1aMjIy7nvOYDDI5MmT5Ztvvqk0db788ku5ffv2A58/ceKE7N27t8x1fv/994c+n52dLfHx8WWuM336dDl//nyJz+Xm5sqWLVskIiKiTDW0eg3QYi4i2szHJI8MytOOHTvQt29fxcbLzc3Fs88+W+bPISqN2NhYeHh4lPduKOKff/5BgwYNyns3FKP0fJ6KBvInn3yiWa33339f0fEe50VerSBQ+8KZgQMHqjr+3ebMmWMydcLCwlSvUSQoKEj1Gt9//73qNQCgd+/eqtfQai6A8vN5KsJg06ZNmtVS+0Crffv2qo5/N4PBoOr4f//9t6rj323btm0mU2fx4sWq1yiixYmDGTNmqF4DMK25AMrPx2RXE5WXZcuWqTp+VlaWquPfTclVREUuXLgAoPAHOTc3F0lJScYf6rI2QB9Gq7OhpnbWdfbs2eW9C4oZPXp0ee+CopSej8n2DBwcHKDT6SAiSE1NxQsvvAARgU6n0/QvUiUcOHAAQOELzejRo43v+wCgzCs7tObt7W38fzEYDGjdurXx/+WHH35QrW5BQYFqy0q1qHP39+348ePw9PTU5PumloCAAON8tm/fXmwF2YoVK8pxz56cqczFZI8MEhISjB+7u7vjxIkT5bg3ZbNy5Urjx1euXMGqVauMLwSVLQx+/PFH48fu7u6avZBpEQRq1lm1ahWAwj8I+vTpU+xnojIaOXKk8eOYmBj4+/uX386UkanMxWTDwJTc/Yvv4eFRqf7aIGU0btzY+HG1atWKPa6MOnfubPy4du3axR5XNqYyl6eigdyhQwdN66WmpiI7O1uVsU3prN7kyZPLexcqJXt7+/LeBUVVrVq1vHdBMZV5LibbMyhPL7/8Ms6dO4eBAwdi4cKFio4dHR2N7t27Kzrmw3z11Vd47rnnMHDgQJiZVe4DycjISDRo0ABt27Y1iTpacXJyAgCMHz8eEyZMKOe9KZuXX34Zzz77LMaPH6/o9UDlRcn5VO7f7gpq7969EBGcOnVK8bG1DAKg8EgkJiYGa9eu1Wx5plqOHDmC33//HXl5edi9e3elr6OVP//8E1euXMHhw4fLe1fKbPXq1UhNTTWJuQDKzodHBgq7desWtmzZgvXr12Pnzp3lvTtUTkzpqvCYmBisX78eISEh5b0rZZaUlITw8HD85z//Ke9dUYSS8+GRgQJycnKwc+dOrFu3Dnv27MHAgQPx1ltvlfdulcrnn39+3zYLCwu0atUKbm5u2u9QGeTm5uK///2vcWlu586d8dZbb6n2Ii0i+OGHH7Bu3Trs2LEDaWlpqtTRwokTJ7Bu3Tps2rQJDg4OGDBgQHnvUqmlp6dj06ZNWL9+PS5evIhXX321vHepTNSaz1MVBqmpqbCyskK1atUUGS86Ohrr169HdHQ0vL294efnh2PHjlXqZX8GgwEGgwGvvPIKgMJ7Lbm6uuJ///sfXnvtNbz99tuK11TrHPvYsWORm5uLcePGAQDWrFmDsWPHFrtOQwmHDx/GunXrsHXrVly9ehUhISGK94q0cObMGaxfvx7r169HvXr1MGTIEIhIseXAlcXNmzexefNmrFu3DmfOnMGAAQOQkJCA5OTk8t61UtFkPmW6zV0l061bN7G3t5d///vfioyn0+nEy8tL/v77b+M2BwcHRcZ+lJCQEAkPD5fc3FxFx+3UqZPcvHnT+PjmzZvi5eUlt2/fFicnJ0VrFZk1a5b07dtXevbsqei4rq6uj7WttGbNmiV6vV66du0qX3/9tVy+fFns7e0VG/9RmjVrJs2aNZOlS5cqMl7Rz/PddybV6udZpPD3s2fPnrJ9+/Yyj1W9enXx8vKSAwcOGO/EW1nnIqLNfJ6qIwOlG7uxsbEIDw/Hyy+/jCZNmmDo0KHIz89XZOxHEZUau5cuXSp25PTss88iLS0NNWrUUOyI6l7z5s1TZdwqVarg3LlzePHFFwEU3gtJyTcF+uabb/Cvf/0LY8eOxSuvvIJq1apBp9MpNv6jKN3Y3bx5M8LDw+Ht7Y2ePXti6NChmi5lVrIZ+sknnyA8PBzjxo3DsGHDMGTIEAX28PEp3ajWZD6KRksFlpmZKWvWrJHevXurMv6hQ4dkwoQJYmNjIz179pRly5apUkdtH3zwgbi7u8ucOXNkzpw50qpVK5k7d65kZmbK8OHDFamRk5MjixcvloEDB8rAgQNlyZIlkpOTo8jYd9u7d680bNhQOnfuLF5eXtK4cWP54YcfFBs/Ly9Pdu/eLX5+fmJrayuvv/66NGjQQPGjtQc5ePCgjBs3TvFxMzMzZe3atdK3b1+pWbOmvPXWW7Jnzx7F69zrwoUL8umnnyo65rlz5+Tjjz8WFxcXqVatmgQHB8vp06cVrVESNeYiou58THo1UUmN3QEDBhjPh6uhoKAAe/fuRXh4uKJXCmvZ2D127Bh+/vlnAIUX7Hl6eio6flBQEHJzc42X7a9ZswZVqlRR/Fw+AGRnZ+P06dMAgKZNm6p2dJOdnY0dO3Zg/fr1OHjwILp164Z169YpXqekxu7EiRMVr1MkIyMDmzZtQnh4uCq3DimpGapWv+XkyZNYt24dNm7ciLNnzyo+vpZzAZSfj0mGwb2N3SFDhmDixIlITExUtE7RHTgfRMm7cA4fPrzExm5iYqLijd38/HykpaUhLy/PuE3JubRs2RJxcXGP3FZamzdvfujzaq+MuXnzJrZs2QI/Pz9Fxiupsbtw4UKcP39ekfG1VlIzdMOGDZWyuWtKczHJMHjmmWfQqVMnrFq1Cg4ODgCAJk2aKH630hYtWhjvVlhEp9MhPT0dly5dUrR/4OXlhV27dqFWrVoAgMzMTPTp0wdRUVFo1aqVYn2QpUuXYu7cuXj++edRpUoV4w3xfvvtN0XGBwrvr7Rp06Zi5/IHDRqE2NhYRcYPCAgAUNj/+Pnnn9GtWzfjqpiXXnoJO3bsUKTO6tWrH/q8UmFQ9PO8fPly4/tdq/HzXLt2bWPPo+hnWqfTIS8vDzk5OcX+OCiLGjVqoE2bNvjoo4/QsWNH6HQ6xedTdNfiIkU/x0DhnM6dO6dIHS3mAmgzH5NsIGvV2P3999+LPU5MTDS+qfz//d//KVpLq8bu4sWLcfr0aTz33HOKjXmvBQsWwNvbG02aNIGI4Pz584ouxy0aq3v37jh16hRsbGwAFC4tvvsOk2V17NixErdv27YNKSkpioWBVo3dmzdvFnucmZmJkJAQLFu2TNG1+Vo0Q+99Y6aCggJs3LgRCxcuhLu7u2J1tGpUazIfRToPFZgWjd0zZ86Iv7+/NGvWTL7++mtVmqFaNHZFRLp06aJJAzQrK0vi4uIkLi5OsrKyVKnRrFmzYo/z8/Pv26aUgoICWbNmjbi4uMjgwYMlLi5O8RpaNXYzMjJk9uzZ4uDgIO+8845cvnxZ8Roi2jR38/PzZdWqVeLs7CwjRoyQP/74Q9Hxi2jVqFZzPiZ5mqgkajR2T548iY8//hh//PEH3n77bQwbNkzRpYv3UruxCwCBgYE4ffo0+vTpU+yIY9q0aWUeW+tz+RMmTEB8fDyGDRsGANiwYQP0ej2WLl2qWI28vDysWrUKCxcuRLt27TBr1iw0bdpUsfEfRI3G7uXLl/HZZ59hw4YNGDVqFCZOnAgLCwtFxn4UpZuhubm5WLFiBb744gt07NgRM2fONJ5iU5sajWot5mOSYaBVY7dKlSpo2LAh+vTpU2IILFmyRJE6RdRu7ALA3LlzS9yuxNsfanUu/26bN2/GwYMHART2XZQ83RESEoLFixejW7dumDFjRqW/tbS5uTmsra0REBCA2rVr3/e8En8QaMXOzg5mZmaYMmVKib8jle32GlrMxyTDQKvGblhY2EOfV/Idj7Ro7Gqle/fuCAsLu+9c/p49e8p5z57MM888g/r168Pa2rrE5p5S/zdaNXbnzJnz0IvmlHo/ZC2aoSNHjnzgXHQ6nWJnB7RqVGsxH5MMg3vd3didNGmSquuyi1y4cEHRv9r1ej2OHDmiWmN3ypQpWLRoEV555ZUSf+iUvMrZyckJf/75p/FxQUEBnJ2di21Twt0vojk5OcjNzYW5uTlu3LihyPiPWtqp1ruR3dvY/eyzz1Spc7dbt27B3NxckbGuXLlS7PHdzVAPDw989913itTRginNxaQbyFo0dn/++WfZtGmTpKWliYhIXFycDBs2TOzs7BSto3Zj12AwiIjI/v37S/ynpPHjx0v37t1l5cqVsnLlSunZs6dMmDBB0Rr3KigokC1btsiMGTNUrSNS2OT79ttvFR9Xi8ZucnKyHDt2TLKzs0VEJC0tTWbNmiU2NjaK11K7uZuXlyfp6enGx9nZ2bJs2TJVFhFo0ahWez4mGQa///67DB06VFq0aCFr1qyRvLw8VepMnz5dmjVrJkOHDhVPT09555135Pnnn5dFixbJnTt3FK01atQo6dChg8ybN08+++wz47/K6rvvvpMpU6bIlClTZPPmzZrVdXNzU2ys69evy7x582T8+PGyZ88eKSgokCVLlkjjxo3F19dXsTrp6ekyc+ZMcXBwkA8//FCuXbum2Nh3++KLL6RevXrSrl07cXd3l6+//lqsrKxkypQpcvHiRcXq5OTkyP/+9z9p2rSpBAYGFrsxnlLWr18vderUERsbG/Hy8pI9e/aIra2t9O/fX44fP65YHS3mIqLNfEzyNJFWjd3mzZsjNjYW1atXR0ZGBho2bIiTJ0+q0khUs7EL/P99lgepjL2Ju1cvFRQUwGAw4KeffsIvv/yiyPj9+vWDpaUl2rdvj3379uHSpUsQESxevFjRW4Ro1dht3rw5YmJiYGVlhQsXLuBf//oXDh06hFatWikyfhEtmqEuLi7YunUr9Ho9YmNj0b59e0RERCh+KxqtGtVazMckLzpT8p5AD1O9enVUr14dAGBpaQlHR0fVVpQo9aL/IEUreYrezeqNN94AAHz77beK34lT7XP5RbZv32782MzMDPb29oiMjFRs/L///tt44WFQUBBsbGxw4cIF48+EUv7zn/8Yv1/3XhimpOrVq8PKygpA4Sq1pk2bKh4EQOH79up0OsTFxd13CxKdTqfIC2jVqlWNSy89PDzg6Oioyj3JtJgLoM18TPLI4GGUbOzWrVsXXl5eAApXERw8eND4GFCm6aplYxcA3N3dceLEiWLbPDw8FLtVxL1EBJGRkTh8+DCCg4NVqaGWe78van6fHkTJxm79+vUxdOhQ4+Pw8PBij5VeKq0mOzu7YkdMn3/+ebHHlWmZLKDNfEw2DH755RekpKTAy8sL9evXx2+//Ybg4GAcPHgQSUlJitT46aefHvp8586dy1zj+PHjaNWq1QNrKVHjbm5ubggJCUGHDh0AAD///DPGjRuHX3/9VdE69yophMrqzJkzGDt2LNLS0nDy5En89ttv2LZtG959911Fxq9SpQrMzc2Nyz3v3LmDmjVrGpcXKnmkk5KSgtTUVLi6uqJq1aq4dOkSFi1ahFWrVuHixYuK1NByqXR+fj4yMjJQr149AIVHiKtWrcIXX3yhyKqyB51WLaLkkbbacwE0mo8inYcKRsvGrkjhbQLy8/ONj/Py8uTWrVuK19GCwWAQV1dXady4sTRq1EhatmypaMNNpLB5XPRv06ZNMmPGDGnXrp2iNUREvLy85MiRI8Waxs7OzorXUZtWjV2taNXc1YIpzcUkw8DJycn4on/16lUxNzeXhIQE1eq1bdv2vreKbN++vSJju7i4SIsWLR74Ty3Xrl1TbdXKyJEjjf+CgoLko48+Mi7NVZKnp6eIFF9B1LJlS0Vr5OXlSdOmTRUd815OTk5y5coVERE5f/68VKtWzbgUWA0vv/yyZGRkGB9fvXpVunfvrtj4zs7OxlU3x48fl6pVq8q2bdsUG/9uly5dko8//lhGjx4tAQEBxn9K0XIuIiJJSUnSv39/sba2FmtraxkwYIAkJSUpMrZJNpC1bOwCQFZWlvHW0gBQq1Yt3L59W5GxtWzsAsD169cxd+5cHDhwAEDhaaj3339f0XvUKHmH0oepV68ezp07Z/w+RUREGK96VkqVKlXQtGlTxS8yvJtWjd0i6enpqFu3rvGxpaUlLl26pNj4WjV3gcIVX506dcLLL7+syn3DtJwLUHhLl+HDh2PTpk0ACl8HAgIC8P3335d5bJMMg7///hu+vr4AChuUCQkJxseA8k1Xc3NzxMbGwsPDA0Dhef4aNWooMnbRVazff/99sXPq8+fPh4eHh+JN11GjRsHFxQUbN24EUPguZAEBAY+8ydyTUPtcfpGQkBCMGTMGf/31F2xtbeHg4IBvv/1W0RpA4U3jnJ2d0aZNm2LNXKV+zpKTkzFp0iTj49TU1GKPlW7sVqlSpVi4nT9/XtE/PC5dulTsnfuuXbtW7LGSzd3bt29j/vz5io13Ly3nAhQGddE9voDC21QsWrRIkbFNMgzuXT44ffp0VestWrQIr732Gl544QWICP755x9s2LBB0RoigkOHDhVr7BYUFChaAwDOnTtX7BL62bNnK/62mqNHj8aCBQvw5ptvAgBcXV0xfPhwxcOgSZMm2Lt3L27duoWCgoIS1+gr4cMPP1Rl3CILFiwo9ljNowIA+Pjjj9GxY0d07tzZuEouNDRUsfFHjx5dbInsvY+V1LdvX+zatQu9e/dWZXwt5wIAzz33HL799lvjnXjXr1+v3C1qFDnZVEFp2djNycmR33//XX7//XdVbnuhRWNXRKRdu3Zy8OBB4+OYmBjFm7tanMsXEVm0aJFcv35dCgoKJDAwUNzd3RW9/39BQYEin1MRpaeny/bt22X79u3FboFQ2dSqVUt0Op1Ur15dateuLbVq1ZLatWuX926VWmJiorzyyitSr149sba2ln79+sn58+cVGdukw0DNxu69Dh06JGvXrpWwsDDjPzWo2dgVETlx4oQxdBo3bixubm7y66+/KlqjZ8+ecvbsWXF3dxcRkU2bNknPnj0VrSEi4urqKiIiUVFR0r9/fzl58qSxphI6d+4sS5Ysue+XMTs7W/bt2yd+fn6ycuVKxeqp3dgtEhMTI5mZmSIismbNGpk6daokJiYqXkft5q6W1GzsasUkTxMVUbOxe7c33ngD586dg5ubm7FJpdPpFHvbQ0Cbxi5QeJ1BXFyccY18nTp1FB0f0O5cvvy/9f+7du2Cn58fnJ2dFX27yKioKKxYsQLDhg1DQkIC6tatizt37qCgoADdu3fHlClTFH2LRbUbu0XGjh1rvKL2888/R2BgIPz8/B55Xc2TUru5CwCHDh2Cm5sbzM3N8e233yI2NvaBt44oCzUbu3dLT0/H119/jcTExGK3LlfkrgvlnUZqeumll4qdSjEYDKqsZ2/WrJnqpwMGDBgg77//vpw7d07OnTsnc+bMkVdffVXxOrNmzbrvr8933nlH8Toihafxbty4ocrYIoVLWH18fESv18utW7fkxo0b4uHhoUqtnJwcuXjxYrHvndI8PDyKHYUkJiYqeqRTpGjMuXPnyjfffFNsm5LUODV4rxYtWkhBQYH8+uuv4ubmJl9++aV4eXkpXqekuagxv/bt28vbb78tGzZskIiICOM/JZh0GBw9elSaNGkiHTt2lA4dOsiLL76oyvrsQYMGqX7xj1Y/bCXd1VPpFwK1z+UXyc/Pl+PHjxtfoC9fvqzKexNrZffu3dKwYUN5/fXXZcSIEdKoUSOJiopSvI6Xl5fMmzdPHB0dJTU1VfLz88XFxUXxOu+8847s3LlT8XHvplWwde3a1XiH5Ly8PFmzZo107dpV8TpqBqhJh4GI+o1dkcL3Gqhbt650795dXnnlFeM/JWnR2BUp/Evq7jeov337tjRv3lzRGmqfy7/bd999J1OnTpVp06ZpeqtstWjR2E1NTZXPPvtMDhw4ICKFF7qp0QPTormrVbCp2di9m5oBarL3Jiry888/33d+Tclz+cCD71Gk5H2Dfv31V/j7++P69esACs8Xr1q1Ci1btlSsBlB4/cL27duNa5lXrlwJX19fvP3224rVcHV1xW+//YbJkyejS5cuePXVV1W5N9G4ceNw9uxZ4zK8DRs24MUXXzRewFfZlHT+e/LkyYq9o5rc9ZaNZfmciuSff/7BunXr0Lp1a3Tq1AkXLlzA/v37FX8N0Ert2rVx69YtVKtWDc8++6yi98Ey6TB4UGNXqYt0yuOXR83GbpGoqCjs3bsXAODj44MePXooOn5AQABSUlKQkJCAuLg45Ofno0uXLjh+/LiidZo1a4Y///zT+P0vKChA8+bN8ddffylaRyuurq6Ii4vDb7/9hoCAAAQGBmLjxo2KNXa7dOmCgQMHol+/fsUarDk5OYiJiUFYWBi8vb0xcuRIReqp2dzV+ndT1cauVlQ53qgg1G7sarm0UMvGbpHt27erMq5W5/L79OlTbElkYmKi9O3bV/E6WlH7/PedO3ckJCREXnrpJbGxsREnJydxcHCQRo0aSVBQkMTGxipWS0Td5q7Wy37VbOzeraRlv7zO4DGo3dgt6ZfH3t5elV8eLRq7Wo6vxbl8Ly8vqVGjhnTu3Fm6dOkiNWvWFC8vL1V6OlrQ6vy3iDaro9QMN62DTYuVUSLqBqhJnyby9vbGr7/+ijZt2qBatWrG7UrfmwgAcnNzcfnyZdSoUaPYWnCluLq64tixY8Z53LlzB56envjjjz8Ur1VEjfP4gHbn8rV4vwktmdr5786dO6Nnz55YuXIlDhw4gPr166Nly5bGd49Titq/mwDw7rvv4qWXXlLtthdFit5A6YMPPoCtrS0CAwMVe1Mlkw4Drd4QRgtaNHbvdfToUbRp00bxcdU+ly8m1gg1tfkUMaVwU7Oxezc1A9Qkw8BUf3nUbuwWUXsFVt++fRESEmJcBXP+/HlMmDCh2HsWl4XWjVC1mdp8TPX3UwtqBqhJhoGp/fLca8eOHejbt68qY6u9Agso/Ovm2LFjaNOmDXQ6HY4ePQpPT0/jrTXKehovKysLK1aswNq1a0u8TcS4ceMUvU2E2kqaT1ZWFvLz8yvlfEzx91Pt215oEaAmGQam9mJwLzXfeN3JyQmnTp1S9a8yLc/la3G+WEumMB9TCzeg+LLfkSNHIigoqNIt+zXJMLibKfzy3Eutxi4AvPbaa1iyZIni7wgG8PQA3c9Ufj/VbOwC2gSoyYeBKVKrsQuouwLLFE8PEAHarYwC1AtQhkElUtlvrWHqp+/o6WUKK6MYBpWEFo1dLZnK6QF6upnSqU+GQSWhdmO3Y8eOiImJQe3atYvVUGu9NJEpMKVTnwyDSkLNxi4RlY4prYxiGFQSWt5aIyYmBvHx8QgICMDly5dx8+ZNODg4KF6HyJRU9lOfDINKQqtba8ydOxcGgwGnT5/GmTNncPHiRbz22ms4dOiQonWIqGIxK+8doMej1f2UtmzZghMnTsDDwwMA8MILL+DmzZua1Cai8vNMee8APVzHjh0BFN4Iq06dOsZ/RY+VVrVqVeh0OmMT+datW4rXIKKKh0cGFVxMTAwAaPbX+eDBg/Hmm2/i2rVr+Prrr7FixQoEBQVpUpuIyg97BpWIVo3d77//HtHR0RAR9OjRAz4+PorXIKKKhWFQSWjV2J0xYwbmz5//yG1EZFrYM6gktmzZgm3btsHc3ByAeo3d77///r5tu3fvVrwOEVUs7BlUEmo3dv/73//iq6++wt9//w1XV1cAhVcfZ2ZmokOHDorWIqKKh2FQSajd2B0+fDh69eqFWbNmITg42Li9du3asLKyUqwOEVVM7BlUIlo1duPi4nDw4EEAQKdOndCyZUtV6hBRxcGeQSUxY8YM+Pj4YMGCBVi4cCF8fHwwY8YMxessWbIEI0aMwKVLl3Dp0iW8/vrrWLp0qeJ1iKhi4ZFBJVHSuya5urrit99+U7SOq6srfvnlF2Oj+tatW2jfvr3idYioYmHPoILTurErIsb3SwCAKlWqgH8vEJk+hkEFp3VjNyAgAG3btsWrr74KANi6dSsCAwMVr0NEFQtPE1Uiajd2CwoKcPjwYVSvXt14G4xOnTpVmvuxE1HpMQwqiSVLliA0NBQDBgwAUHgR2pgxYzBx4kRF67i7u+PEiROKjklEFR/DoJLQqrE7ffp0tG/fHgMGDKgU79tKRMrg0tJKQqvG7rJly/Daa6+hatWqqF27tmq3yiaiioUN5EpCq8Yu38iG6OnE00SVgNaN3c2bNyMmJgY6nQ6dOnVC//79ValDRBUHw6CS0KqxO27cOJw9exbDhg0DAGzYsAEvvvgiQkJCVK9NROWHYVBJaNXYbdasGf78809jjYKCAjg7O+PPP/9UrSYRlT82kCsJrRq7er0eFy5cMD5OSkqCXq9XvA4RVSw8MqBiOnfujGPHjqFNmzbQ6XQ4evQoPD09YWFhAQDYtm1bOe8hEamBq4kqES0aux988IHiYxJRxccjg0pC7cauiDyyF/E4n0NElRPDoJJQu7HbpUsXDBw4EP369UOjRo2M23NychATE4OwsDB4e3tj5MiRitQjooqFp4kqiaLGbuPGjQEo39iNiorCihUrMGzYMCQkJKBu3bq4c+cOCgoK0L17d0yZMoU3rCMyYTwyqCS0bOzm5ubi8uXLqFGjBurWravYuERUcTEMKomffvrpoc937txZoz0hIlPEMKjg2NglIi3worMKztvbG0uXLi12IRhQ2Nj94Ycf4O/vj7CwsHLaOyIyFTwyqOCysrKwYsUKrF27tsTG7rhx49jYJaIyYxhUImzsEpFaGAZERMSeARERMQyIiAgMAyKjK1euwM3NDW5ubmjQoAFsbW2Nj3NycjTZh3nz5qk6vr29PS5fvqxqDaqc2DMgKsGcOXNQq1YtTJ8+XdO6tWrVQmZm5n3bRQQigmeeKdvfb/b29jAYDKhXr16ZxiHTwyMDqnQ+/PBDNG3aFB07dsSwYcOwcOFCnDt3Dj179kSrVq3QqVMn/PXXXwCAkSNHYtKkSXjppZfQpEkTREREGMdZsGABWrduDVdXV8yePbvEWjdv3oSDgwNyc3MBADdu3DA+7tKlCyZPngw3Nze4uLjg6NGjAIBbt25h1KhRaNOmDdzd3REZGXnfuKmpqfDy8jJ+7cGDBzFz5kzcuXMHbm5uGDFiBBITE9G0aVP4+fnBxcUFSUlJGDt2LDw9PeHs7Fxsn+3t7TF79mx4eHigRYsWxvlfuXIF3bt3h7OzM4KCgsC//eiBhKgSOXr0qLRs2VLu3LkjN27cEL1eLwsWLJCuXbvKmTNnRETk8OHD4u3tLSIi/v7+MmjQIMnPz5c//vhDXnzxRRER2bNnj4wePVoKCgokPz9f+vTpIz/99JOxzuzZs2XBggUiIjJy5EjZsmWLiIgsW7ZMpk2bJiIinTt3lqCgIBER+emnn8TZ2VlERGbNmiVr1qwREZGMjAxxdHSUzMzMYvNYuHChfPTRRyIikpeXJzdu3BAREXNzc+PnJCQkiE6nk19++cW47cqVK8av6dy5s8TFxYmISOPGjWXJkiUiIhISEiKBgYEiIjJx4kSZO3euiIjs2LFDAEh6evoTf9/J9PGupVSpHDp0CP369UP16tVRvXp1vPLKK8jKysLPP/+M1157zfh52dnZxo/79++PZ555Bs2bN0daWhoAIDo6GtHR0cYL9jIzMxEfHw8vL6/7agYFBeHTTz9F//79sXLlSnz99dfG54reX8LLyws3btzAtWvXEB0djW3btmHhwoUACi8cvHDhApycnIxf17p1a4waNQq5ubno378/3NzcSpxv48aN0a5dO+PjjRs3IjQ0FHl5eUhNTcWpU6fg6uoKABgwYAAAoFWrVti8eTMA4MCBA8aP+/TpA0tLy8f5NtNTiGFAlV5BQQHq1q2LX3/9tcTnq1WrZvxY/t9pEhHBrFmz8Oabbz5y/A4dOiAxMRH79+9Hfn4+XFxcjM/de08onU4HEcF3332Hpk2bPnBMLy8vHDhwADt37sTIkSMxbdo0+Pn53fd55ubmxo8TEhKwcOFCHDt2DJaWlhg5ciSysrLum2eVKlWQl5f3yHkR3Y09A6pUOnTogO3btyMrKwuZmZnYsWMHatasCQcHB2zatAlA4Qt9XFzcQ8fp0aMHVqxYYWzWpqSk4NKlSw/8fD8/PwwfPhwBAQHFtm/YsAEAEBMTAwsLC1hYWKBHjx5YunSpMXhOnDhhrNGtWzcAwPnz5/H8889j9OjRCAoKQmxsLADg2WefNfYn7nXjxg2Ym5vDwsICaWlp2L1790PnCBSGzrp16wAAu3fvRkZGxiO/hp5ODAOqVFq3bg1fX1+4urqiV69eaNGiBSwsLLB27VosX74cLVu2hLOzc4lN27t1794dw4cPR/v27dGiRQsMGjQIN2/efODnjxgxAhkZGcbTQkWqV68Od3d3vPXWW1i+fDkA4L333kNubi5cXV3h7OyM9957D0Bh09jMrPBgfP/+/WjZsiXc3d2xYcMGTJ48GQAwZswYuLq6YsSIEfftQ9HnN2vWDMOHD0eHDh0e+f2aPXs2Dhw4AGdnZ2zevLnYu9gR3Y1LS6nSyczMRK1atXD79m14eXkhNDQUHh4eqtaMiIhAZGQk1qxZY9zWpUsXLFy4EJ6eno81xpdffolGjRrB19dXrd0kKjX2DKjSGTNmDE6dOoWsrCz4+/urHgQTJ07E7t27sWvXrjKNM2HCBIX2iEh5PDIgIiL2DIiIiGFARERgGBARERgGREQEhgEREQH4/wDypw6eZQHMKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# by geneType and strand ...\n",
    "genes_df.groupby([\"geneType\",\"strand\"]).size().plot.bar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "yellow-magnet",
   "metadata": {},
   "source": [
    "### merging RNA-seq expression data with our gene annotation ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "practical-desire",
   "metadata": {},
   "outputs": [],
   "source": [
    "rsem_df = pd.read_csv(\"merged_TPM_genes.tsv\",sep=\"\\t\")\n",
    "# genes_Df  overlap   rsem\n",
    "# ~4_300     ~23_200   ~2_500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "flying-sentence",
   "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>Id</th>\n",
       "      <th>12</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>MIR941-2</td>\n",
       "      <td>MIR941-2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TTTY1B</td>\n",
       "      <td>TTTY1B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>F8A1</td>\n",
       "      <td>F8A1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>C9orf16</td>\n",
       "      <td>C9orf16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>CXorf56</td>\n",
       "      <td>CXorf56</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Id        12\n",
       "0  MIR941-2  MIR941-2\n",
       "3    TTTY1B    TTTY1B\n",
       "5      F8A1      F8A1\n",
       "8   C9orf16   C9orf16\n",
       "9   CXorf56   CXorf56"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAD4CAYAAAAdIcpQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZI0lEQVR4nO3df1BU1/3/8ReGkLGmphWlibumgKuIoIkRdJzpGJv6I2KxRhEhRkPUrEZaTVpNO5/O9NdYNZ1pm7aYWCwxzkSx6cQWkYI1Nc04xqqI00SJDFVsAROj4I8ERQXO94983YoKd1fZu3vl+ZjJDHvZe/d9uHFfc865954IY4wRAACd6BHqAgAA4Y+wAABYIiwAAJYICwCAJcICAGApMtQF3I6+ffsqNjY21GUAgKMcP35cp0+fDmgfR4ZFcXGxiouL9eUvf1nl5eWhLgcAHCUlJSXgfRw5DJWenq78/Hzdd999oS4FALoFR4YFAMBehAUAwBJhAQCw5OgJ7nPnzoW6FADoFhzZs2CCGwDs5ciwAADYi7AAAFhy5JxFl6gq/d/PCZNDVwcAOAA9CwCAJcICAGCJsAAAWHLknAX3WQCAvRzZs+A+CwCwlyPDAgBgL8ICAGCJsAAAWCIsAACWCAsAgCXCAgBgibAAAFgiLAAAlsIqLJqampSSkqJt27aFuhQAwDWCGhbz5s1TTEyMkpOT220vKytTQkKCPB6PVq9e7dv+0ksvKTMzM5glAQBuQVDDIicnR2VlZe22tba2Kjc3V6WlpaqsrFRhYaEqKyu1Y8cODR06VDExMcEsCQBwC4L6IMGxY8fq+PHj7bbt27dPHo9H8fHxkqSsrCwVFRXps88+U1NTkyorK9WzZ0+lpaWpR4+wGiUDgG7L9qfO1tfXa8CAAb7Xbrdbe/fuVV5eniTp9ddfV9++fTsMivz8fOXn50uSTp06FfyCAQDh94jynJycTn/v9Xrl9XolSSkpKTZUBACwPSxcLpdqa2t9r+vq6uRyuQI6ButZAIC9bJ8USE1NVXV1tWpqanT58mVt3rxZU6dODegYrGcBAPYKalhkZ2drzJgxqqqqktvtVkFBgSIjI5WXl6dJkyYpMTFRmZmZSkpKCui4xcXF8nq99CwAwCYRxhgT6iJuVUpKisrLy29t56rS//2cMLlrCgIAB7iV706uTQUAWAq7q6H8wQQ3ANjLkT0LJrgBwF6ODAsAgL0YhgIAWHJkz4JhKACwlyPDAgBgL4ahAACWHNmzYBgKAOzlyLAAANiLsAAAWCIsAACWmOAGAFhyZM+CCW4AsJcjwwIAYC/CAgBgibAAAFhighsAYMmRPQsmuAHAXo4MCwCAvQgLAIAlwgIAYKnbhsX8Dfs1f8P+UJcBAI7QbcMCAOA/wgIAYIn7LAAAlhzZs+A+CwCwlyPDAgBgL8ICAGCJsAAAWCIsAACWCAsAgCXCAgBgibAAAFgiLAAAlsImLD788EMtWrRIGRkZevXVV0NdDgDgGkENi3nz5ikmJkbJycnttpeVlSkhIUEej0erV6+WJCUmJmrt2rV68803tXv37mCWBQAIUFDDIicnR2VlZe22tba2Kjc3V6WlpaqsrFRhYaEqKyslSVu3btWUKVOUlpYWzLIAAAEKaliMHTtWffr0abdt37598ng8io+PV1RUlLKyslRUVCRJmjp1qkpLS7Vx48ZglgUACJDtT52tr6/XgAEDfK/dbrf27t2rf/zjH9qyZYsuXbrUac8iPz9f+fn5kqRTp04FvV4AQBg9onzcuHEaN26c5fu8Xq+8Xq8kKSUlJchVAQCkEISFy+VSbW2t73VdXZ1cLldAx2A9CwCwl+2Xzqampqq6ulo1NTW6fPmyNm/erKlTpwZ0DNazAAB7BTUssrOzNWbMGFVVVcntdqugoECRkZHKy8vTpEmTlJiYqMzMTCUlJQV03OLiYnm9XnoWAGCTCGOMCXURtyolJUXl5eW3tO/8//uZJKng6dT2v0iYfLtlAUBYu5XvzrC5gxsAEL7C5mqoQDDBDQD2cmTPggluALCXI8MCAGAvhqEAAJYc2bNgGAoA7OXIsAAA2MuvsPjggw+CXUdAuCkPAOzlV1gsXrxYo0aN0iuvvBIWX9AMQwGAvfwKi127dmnjxo2qra3VyJEj9eSTT2rHjh3Brg0AECb8nrMYNGiQVqxYoZdeeknvvvuulixZoiFDhmjLli3BrA8AEAb8Cov3339fL7zwghITE7Vz504VFxfrww8/1M6dO/XCCy8Eu0YAQIj5dZ/Fd77zHS1YsEArV65Uz549fdv79++vFStWBK24jnCfBQDYy6+wKCkpUc+ePXXXXXdJktra2tTc3KwvfOELmjNnTlALvJn09HSlp6ezUh4A2MSvYajx48fr4sWLvtcXLlzQ+PHjg1YUACC8+BUWzc3Nuvfee32v7733Xl24cCFoRQEAwotfYdGrVy9VVFT4Xh84cKDd3AUA4M7m15zFyy+/rJkzZ6p///4yxujjjz/WH//4x2DX1iEmuAHAXn6FRWpqqo4cOaKqqipJUkJCgu6+++6gFtYZJrgBwF5+P6J8//79On78uFpaWnxDUnPnzg1aYQCA8OFXWMyZM0dHjx7Vww8/7Lt8NiIi4s4Mi6rS//2cMDl0dQBAGPErLMrLy1VZWamIiIhg1wMACEN+XQ2VnJysjz/+ONi1AADClF89i9OnT2vo0KEaNWqU7rnnHt/2rVu3Bq0wAED48CssfvKTnwS5DABAOPMrLB599FH95z//UXV1tcaPH68LFy6otbU12LV1iPssAMBefs1ZrFu3ThkZGVq4cKEkqb6+XtOmTQtmXZ1ipTwAsJdfYbFmzRrt3r1bvXv3lvT5QkiffPJJUAsDAIQPv8LinnvuUVRUlO91S0sLl9ECQDfiV1g8+uijWrlypS5evKgdO3Zo5syZSk9PD3ZtAIAw4VdYrF69Wv369dOwYcP0+9//XmlpaSFZIQ8AEBp+XQ3Vo0cPPfvss3r22WeDXQ8AIAz5FRZxcXE3naM4duxYlxcEAAg/fj8b6qrm5mb96U9/UmNjY9CKAgCEF7/mLKKjo33/uVwuPf/88yopKQl2bQCAMOFXz+LaJVXb2tpUXl6ulpaWoBUFAAgvfoXF9773vf/tEBmp2NhYvfnmm11ezF/+8heVlJTo/Pnzmj9/viZOnNjlnxEQ1rYAAEl+hsU777xzyx8wb948bdu2TTExMTp06JBve1lZmZYuXarW1lYtWLBAP/jBDzRt2jRNmzZNZ86c0bJly0IfFgAASX6Gxa9+9atOf//d7363w9/l5OTo29/+drtV9VpbW5Wbm6sdO3bI7XYrNTVVU6dO1dChQyVJK1asUG5urj+lAQBs4NcEd3l5uV599VXV19ervr5ea9euVUVFhT799FN9+umnne47duxY9enTp922ffv2yePxKD4+XlFRUcrKylJRUZGMMfr+97+vyZMn65FHHrnp8fLz85WSkqKUlBSdOnXKz2YCAG6HXz2Luro6VVRU6Itf/KKkz9e3mDJlit54441b+tD6+noNGDDA99rtdmvv3r363e9+p7ffflvnzp3Tv//9by1atOiGfb1er7xeryQpJSXllj4fABAYv8Li5MmT7R4kGBUVpZMnT3Z5MUuWLNGSJUss38d6FgBgL7/CYu7cuRo1apSeeOIJSZ9ftfT000/f8oe6XC7V1tb6XtfV1cnlcvm9f3p6utLT0+lZAIBN/AqLH/7wh5o8ebJ27dolSVq/fr1GjBhxyx+ampqq6upq1dTUyOVyafPmzdq0aZPf+9OzAAB7+TXBLUkXLlxQ7969tXTpUrndbtXU1Pi1X3Z2tsaMGaOqqiq53W4VFBQoMjJSeXl5mjRpkhITE5WZmamkpCS/i+7KlfLmb9h/28cAgDudXz2Ln/70pyovL1dVVZWeeeYZXblyRU899ZR2795tuW9hYeFNt6elpSktLS2wagEAIeFXWPz5z3/WwYMHfZez9u/f3/KS2WAKyTAUd3MD6Mb8GoaKiopSRESE7zHlTU1NQS3KSlcOQwEArPkVFpmZmVq4cKHOnj2rdevWafz48SyEBADdiOUwlDFGs2bN0pEjR9S7d29VVVXpZz/7mSZMmGBHfTfF1VAAYC/LsIiIiFBaWpo++OCDkAbEtbjPAgDs5dcw1COPPKL9+7nEFAC6K7+uhtq7d6/eeOMNxcbGqlevXjLGKCIiQu+//36w67sphqEAwF6dhsV///tfPfjgg9q+fbtd9fgl5MNQ115GK3EpLYA7XqdhMW3aNFVUVOirX/2qZsyYobfeesuuugAAYaTTOQtjjO/nY8eOBb0YAEB46jQsrt6Ed/3PAIDupdNhqH/961/q3bu3jDG6ePGievfuLUm+Ce7z58/bUuT1mOAGAHt1Ghatra121RGQkE9wA0A34/cjygEA3RdhAQCwRFgAACz5dQd3uOnqCe6rq+UVPJ16awdgrQsAdzhH9ixYzwIA7OXIsAAA2IuwAABYIiwAAJYICwCAJcLiGlevigIAtEdYAAAscZ8FAMCSI3sW3GcBAPZyZM8irHE3N4A7kCN7FgAAexEWAABLhAUAwBJhAQCwRFgAACwRFgAAS4QFAMAS91kEE/dcALhDhE3P4tixY5o/f74yMjJCXQoA4DpBDYt58+YpJiZGycnJ7baXlZUpISFBHo9Hq1evliTFx8eroKAgmOWEVlXpzf8DAAcIaljk5OSorKys3bbW1lbl5uaqtLRUlZWVKiwsVGVlZTDLAADcpqCGxdixY9WnT5922/bt2yePx6P4+HhFRUUpKytLRUVFfh8zPz9fKSkpSklJ0alTp7q6ZADATdg+Z1FfX68BAwb4XrvdbtXX16uhoUGLFi3SwYMHtWrVqg7393q9Ki8vV3l5ufr162dHyQDQ7YXN1VDR0dFau3atX+9lPQsAsJftPQuXy6Xa2lrf67q6OrlcroCOwXoWAGAv23sWqampqq6uVk1NjVwulzZv3qxNmzYFdIxg9iyursNd8HRqlx/bEvdlAAhTQe1ZZGdna8yYMaqqqpLb7VZBQYEiIyOVl5enSZMmKTExUZmZmUpKSgrouPQsAMBeQe1ZFBYW3nR7Wlqa0tLSgvnRAIAuFDZ3cAeiuLhYXq83qBPcV4ejQoYb9wCEEUeGBcNQAGAvR4YFAMBeYXOfRSC63X0W/gxFcfUUgCByZM+CYSgAsJcjwwIAYC+GoTpx7Q16Ib1Z73YEeqNfsN8PwJEc2bNgGAoA7OXIsAAA2IuwAABYIiwAAJaY4IZ/QjmRzSQ6EHKO7FkwwQ0A9nJkWAAA7EVYAAAsERYAAEtMcPshqGtbOHG9Cn9r9mdiuqNjMZENhBVH9iyY4AYAezkyLAAA9iIsAACWCAsAgCXCAgBgibAAAFgiLAAAlrjP4hZcv2re/A37262gd6ur6l1/nIDczj0N/hyzK98fqntLgvFAQh5yiG7CkT0L7rMAAHs5MiwAAPYiLAAAlggLAIAlwgIAYImwAABYIiwAAJYICwCAJcICAGCJsAAAWAqbx300NTVp8eLFioqK0rhx4zR79uxQlwQA+P+C2rOYN2+eYmJilJyc3G57WVmZEhIS5PF4tHr1aknSli1blJGRoXXr1mnr1q3BLAsAEKCghkVOTo7KysrabWttbVVubq5KS0tVWVmpwsJCVVZWqq6uTgMGDJAk3XXXXcEsCwAQoKAOQ40dO1bHjx9vt23fvn3yeDyKj4+XJGVlZamoqEhut1t1dXV6+OGH1dbW1uEx8/PzlZ+fL0k6depU0GrvyNUnyvqzvaueRusElm0LxpNmb+eYge4bwBNlY39Qom/0ONDx3+LaYznlSbh2P12Xp/neXAj/LrZPcNfX1/t6EJLkdrtVX1+v6dOn66233tJzzz2n9PT0Dvf3er0qLy9XeXm5+vXrZ0fJANDthc0Ed69evbR+/Xq/3hvq9SwAoLuxvWfhcrlUW1vre11XVyeXyxXQMVjPAgDsZXtYpKamqrq6WjU1Nbp8+bI2b96sqVOnBnSM4uJieb1eehYAYJOghkV2drbGjBmjqqoqud1uFRQUKDIyUnl5eZo0aZISExOVmZmppKSkgI5LzwIA7BXUOYvCwsKbbk9LS1NaWlowPxoA0IXCZoI7EExwA4C9HPlsKIahAMBejgwLAIC9IowxJtRF3Kq+ffsqNjb2lvY9derUHXVTH+0Jb3dSe+6ktkjdsz3Hjx/X6dOnAzquo8PidqSkpKi8vDzUZXQZ2hPe7qT23EltkWiPvxiGAgBYIiwAAJa6bVh4vd5Ql9ClaE94u5Pacye1RaI9/uq2cxYAAP91254FAMB/hAUAwFK3DIubrQEebmpra/X1r39dQ4cOVVJSkn7zm99IkhobGzVhwgQNGjRIEyZM0JkzZyRJxhgtWbJEHo9Hw4cPV0VFhe9YGzZs0KBBgzRo0CBt2LAhJO25qrW1VSNGjNA3v/lNSVJNTY1Gjx4tj8ejWbNm6fLly5KkS5cuadasWfJ4PBo9enS7FRdXrVolj8ejhIQEbd++PRTNkCSdPXtWGRkZGjJkiBITE7Vnzx7Hnp9f//rXSkpKUnJysrKzs9Xc3Oy4czNv3jzFxMQoOTnZt60rz8eBAwc0bNgweTweLVmyRMEcwb9ZW5YvX64hQ4Zo+PDheuKJJ3T27Fnf7zr6u3f0XdfRue2U6WZaWlpMfHy8OXr0qLl06ZIZPny4OXz4cKjLusGJEyfMgQMHjDHGnD9/3gwaNMgcPnzYLF++3KxatcoYY8yqVavMiy++aIwxpqSkxDz++OOmra3N7Nmzx4waNcoYY0xDQ4OJi4szDQ0NprGx0cTFxZnGxsbQNMoY88tf/tJkZ2ebKVOmGGOMmTlzpiksLDTGGLNw4ULzyiuvGGOMWbNmjVm4cKExxpjCwkKTmZlpjDHm8OHDZvjw4aa5udkcO3bMxMfHm5aWlhC0xJi5c+eadevWGWOMuXTpkjlz5owjz09dXZ2JjY01Fy5cMMZ8fk7Wr1/vuHPz7rvvmgMHDpikpCTftq48H6mpqWbPnj2mra3NPP744+avf/2rrW3Zvn27uXLlijHGmBdffNHXlo7+7p1913V0bjvT7cLivffeMxMnTvS9XrlypVm5cmUIK/LP1KlTzd/+9jczePBgc+LECWPM54EyePBgY4wxXq/XbNq0yff+q+/btGmT8Xq9vu3Xv89OtbW15rHHHjN///vfzZQpU0xbW5uJjo72/QO49txMnDjRvPfee8YYY65cuWKio6NNW1vbDefr2vfZ6ezZsyY2Nta0tbW12+7E81NXV2fcbrdpaGgwV65cMVOmTDFlZWWOPDc1NTXtvmC76nycOHHCJCQk+LZf/z472nKtLVu2mCeffNIYc+N32NW/e0ffdZ39u+tMtxuG6mgN8HB2/PhxHTx4UKNHj9bJkyf1wAMPSJLuv/9+nTx5UlLH7Qqn9j7//PP6xS9+oR49Pv/frqGhQV/60pcUGRl5Q23X1h0ZGan77rtPDQ0NYdOempoa9evXT88884xGjBihBQsWqKmpyZHnx+VyadmyZXrwwQf1wAMP6L777tPIkSMde26u1VXno76+Xm63+4btofLaa69p8uTJkgJvS2f/7jrT7cLCaT777DPNmDFDL7/8snr37t3udxEREYqIiAhRZYHZtm2bYmJiNHLkyFCX0iVaWlpUUVGh5557TgcPHlSvXr1umP9yyvk5c+aMioqKVFNToxMnTqipqUllZWWhLqvLOeV8WPn5z3+uyMhIzZ4929bP7XZh0RVrgNvlypUrmjFjhmbPnq3p06dLkr7yla/oo48+kiR99NFHiomJkdRxu8Klvbt379bWrVsVGxurrKws7dy5U0uXLtXZs2fV0tJyQ23X1t3S0qJz584pOjo6bNrjdrvldrs1evRoSVJGRoYqKioceX7efvttxcXFqV+/frr77rs1ffp07d6927Hn5lpddT5cLpfq6upu2G63119/Xdu2bdPGjRt9wRdoW6Kjozs8t526nTE1J7py5YqJi4szx44d8036HDp0KNRl3aCtrc3MmTPHLF26tN32ZcuWtZuwW758uTHGmG3btrWbsEtNTTXGfD5hFxsbaxobG01jY6OJjY01DQ0Ntrbleu+8845vgjsjI6PdRNuaNWuMMcbk5eW1m0SdOXOmMcaYQ4cOtZvMi4uLC9kE99e+9jVz5MgRY4wxP/7xj82yZcsceX7++c9/mqFDh5qmpibT1tZm5s6da37729868txcP87flefj+gnukpISW9tSWlpqEhMTzSeffNLufR393Tv7ruvo3Ham24WFMZ9fCTFo0CATHx9vVqxYEepybmrXrl1Gkhk2bJh56KGHzEMPPWRKSkrM6dOnzWOPPWY8Ho/5xje+4fsfua2tzSxevNjEx8eb5ORks3//ft+xCgoKzMCBA83AgQPNa6+9Fqom+VwbFkePHjWpqalm4MCBJiMjwzQ3NxtjjLl48aLJyMgwAwcONKmpqebo0aO+/VesWGHi4+PN4MGDg3pFipWDBw+akSNHmmHDhplvfetbprGx0bHn50c/+pFJSEgwSUlJ5qmnnjLNzc2OOzdZWVnm/vvvN5GRkcblcpk//OEPXXo+9u/fb5KSkkx8fLzJzc294eKGYLdl4MCBxu12+74Prga2MR3/3Tv6ruvo3HaGx30AACx1uzkLAEDgCAsAgCXCAgBgibAAAFgiLAAAlggLAIAlwgIAYOn/AeW75bfCmgf2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "yyy = pd.merge(\n",
    "    rsem_df,\n",
    "    genes_df,\n",
    "    how = 'left',\n",
    "    left_on=\"Id\",\n",
    "    right_on=\"geneName\"\n",
    ")\n",
    "\n",
    "yyy[yyy.isna().any(axis=1)][\"rsem.out.S442_NT\"].plot.hist(log=True,bins=100)\n",
    "yyy[\"rsem.out.S442_NT\"].plot.hist(log=True,bins=100,alpha=0.3)\n",
    "\n",
    "display(pd.merge(\n",
    "    yyy[yyy.isna().any(axis=1)][\"Id\"].reset_index(drop=True),\n",
    "    pd.read_csv(\"./gene_annot_data/hg19/refGene_sorted\",sep=\"\\t\",header=None)[12],\n",
    "    how=\"inner\",\n",
    "    left_on=\"Id\",\n",
    "    right_on=12\n",
    ").drop_duplicates().head(5))\n",
    "# pd.read_csv(\"./gene_annot_data/hg19/refGene_sorted\",sep=\"\\t\",header=None)[12]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "laden-hamilton",
   "metadata": {},
   "source": [
    "### after \"filtering\" our list of genes we are left with ~ 27k of them, and let's get to extracting TSS/TTS and other useful info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "alive-variance",
   "metadata": {},
   "outputs": [],
   "source": [
    "# columns needed for TSS/TTS bed files output:\n",
    "tx_cols = [\"chr\", \"txStart\", \"txEnd\"]\n",
    "\n",
    "pgenes_df =  genes_df[ genes_df[\"strand\"] == \"+\" ][tx_cols].reset_index(drop=True)\n",
    "mgenes_df =  genes_df[ genes_df[\"strand\"] == \"-\" ][tx_cols].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genetic-crowd",
   "metadata": {},
   "source": [
    "#### carefully extract TSS from the genes, considering their orientation, and possible duplications ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "unavailable-wealth",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " there are 66 duplicated pTSSs ...\n",
      "finally there are 13921 +TSSs to work with...\n",
      " there are 58 duplicated mTSSs ...\n",
      "finally there are 13476 -TSSs to work with...\n",
      "there are still 1 duplicated TSSs ... drop em\n",
      "finally there are 27396 -TSSs to work with...\n"
     ]
    }
   ],
   "source": [
    "pTSS = pgenes_df.copy()\n",
    "mTSS = mgenes_df.copy()\n",
    "# after copying, pTSS and mTSS are dataframes that have chr, txStart, txEnd  of the\n",
    "# transcripts, that we could modify as we please ...\n",
    "\n",
    "# for genes on p-strand, we keep txStart, and modify txEnd to be ~several bp downstream:\n",
    "pTSS[\"txEnd\"] = pTSS[\"txStart\"] + 1\n",
    "pTSS[\"strand\"] = \"+\"\n",
    "# report how many TSS were duplicated here:\n",
    "print(f\" there are {pTSS.duplicated().sum()} duplicated pTSSs ...\")\n",
    "pTSS = pTSS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(pTSS)} +TSSs to work with...\")\n",
    "\n",
    "# for genes on m-strand, we keep txEnd, and modify txStart to be ~several bp downstream (relative to the gene):\n",
    "mTSS[\"txStart\"] = mTSS[\"txEnd\"] - 1\n",
    "mTSS[\"strand\"] = \"-\"\n",
    "# report how many TSS were duplicated here:\n",
    "print(f\" there are {mTSS.duplicated().sum()} duplicated mTSSs ...\")\n",
    "mTSS = mTSS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(mTSS)} -TSSs to work with...\")\n",
    "\n",
    "# after modifying pTSS and mTSS separately - let's bring them back together ...\n",
    "TSS = pd.concat([pTSS,mTSS]).reset_index(drop=True)\n",
    "print(f\"there are still {TSS.duplicated(subset=tx_cols).sum()} duplicated TSSs ... drop em\")\n",
    "TSS = TSS.drop_duplicates(subset=tx_cols).reset_index(drop=True)\n",
    "print(f\"finally there are {len(TSS)} -TSSs to work with...\")\n",
    "# we save them in the other notebook ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "italic-travel",
   "metadata": {},
   "source": [
    "#### carefully extract TTS (aka 3' end of the gene ) from the list of genes, considering their orientation, and possible duplications ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "beneficial-newspaper",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " there are 120 duplicated pTTSs ...\n",
      "finally there are 13867 +TTSs to work with...\n",
      " there are 68 duplicated mTTSs ...\n",
      "finally there are 13466 -TTSs to work with...\n",
      "there are still 2 duplicated TTSs ... drop em\n",
      "finally there are 27331 -TTSs to work with...\n"
     ]
    }
   ],
   "source": [
    "pTTS = pgenes_df.copy()\n",
    "mTTS = mgenes_df.copy()\n",
    "# after copying, pTTS and mTTS are dataframes that have chr, txStart, txEnd  of the\n",
    "# transcripts, that we could modify as we please ...\n",
    "\n",
    "# for genes on p-strand, we keep txEnd, and modify txStart to be ~several bp upstream:\n",
    "pTTS[\"txStart\"] = pTTS[\"txEnd\"] - 1\n",
    "pTTS[\"strand\"] = \"+\"\n",
    "# report how many TTS were duplicated here:\n",
    "print(f\" there are {pTTS.duplicated().sum()} duplicated pTTSs ...\")\n",
    "pTTS = pTTS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(pTTS)} +TTSs to work with...\")\n",
    "\n",
    "# for genes on m-strand, we keep txStart, and modify txEnd to be ~several bp upstream (relative to the gene):\n",
    "mTTS[\"txEnd\"] = mTTS[\"txStart\"] + 1\n",
    "mTTS[\"strand\"] = \"-\"\n",
    "# report how many TTS were duplicated here:\n",
    "print(f\" there are {mTTS.duplicated().sum()} duplicated mTTSs ...\")\n",
    "mTTS = mTTS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(mTTS)} -TTSs to work with...\")\n",
    "\n",
    "# after modifying pTTS and mTTS separately - let's bring them back together ...\n",
    "TTS = pd.concat([pTTS,mTTS]).reset_index(drop=True)\n",
    "print(f\"there are still {TTS.duplicated(subset=tx_cols).sum()} duplicated TTSs ... drop em\")\n",
    "TTS = TTS.drop_duplicates(subset=tx_cols).reset_index(drop=True)\n",
    "print(f\"finally there are {len(TTS)} -TTSs to work with...\")\n",
    "# we save them in the other notebook ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "herbal-simpson",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !scp allTTS_wstrand.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp allTSS_wstrand.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp ActiveTSS-13412-1bp.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp ActiveTSS-NO-CTCF-8849-1bp.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp activeTTS_wstrand.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp ActiveTTS-withoutH3K4-1bp.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/\n",
    "# !scp ActiveTTS-withoutH3K4-NO-CTCF-10933-1bp.bed ghpcc:/nl/umw_job_dekker/users/sv49w/ALV2021/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abandoned-telephone",
   "metadata": {},
   "source": [
    "### Let's get TTS corresponding to active TSSs\n",
    "using our `ActiveTSS-13412-1bp.bed` file , and our main `gene_df` \"database\"\n",
    "\n",
    "also we need to do this for `ActiveTSS-12216-1bp_HCT116_HAP1.bed` file - list of TSS active between HCT116 and HAP1 ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "retained-medicaid",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fname = \"ActiveTSS-13412-1bp.bed\"\n",
    "fname = \"ActiveTSS-12113-1bp_HCT116_HAP1.bed\"\n",
    "act_tss_df = bioframe.read_table(\n",
    "    fname,\n",
    "    names=[\"chrom\",\"start\",\"end\",\"strand\"]\n",
    ")\n",
    "act_tss_df = act_tss_df.sort_values([\"chrom\",\"start\"])\n",
    "act_tss_df = act_tss_df.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "built-romance",
   "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>chrom2</th>\n",
       "      <th>start2</th>\n",
       "      <th>end2</th>\n",
       "      <th>strand2</th>\n",
       "      <th>geneId</th>\n",
       "      <th>chr</th>\n",
       "      <th>txStart</th>\n",
       "      <th>txEnd</th>\n",
       "      <th>geneName</th>\n",
       "      <th>citationCount</th>\n",
       "      <th>strand</th>\n",
       "      <th>union_geneId</th>\n",
       "      <th>geneType</th>\n",
       "      <th>geneDesc</th>\n",
       "      <th>cdsStart</th>\n",
       "      <th>cdsEnd</th>\n",
       "      <th>exonStarts</th>\n",
       "      <th>exonEnds</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>29369</td>\n",
       "      <td>29370</td>\n",
       "      <td>-</td>\n",
       "      <td>653635.0</td>\n",
       "      <td>chr1</td>\n",
       "      <td>14361</td>\n",
       "      <td>29370</td>\n",
       "      <td>WASH7P</td>\n",
       "      <td>4</td>\n",
       "      <td>-</td>\n",
       "      <td>union_653635</td>\n",
       "      <td>pseudo</td>\n",
       "      <td>WASP family homolog 7, pseudogene</td>\n",
       "      <td>29370</td>\n",
       "      <td>29370</td>\n",
       "      <td>14361,14969,15795,16606,16857,17232,17605,1791...</td>\n",
       "      <td>14829,15038,15947,16765,17055,17368,17742,1806...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>714013</td>\n",
       "      <td>714014</td>\n",
       "      <td>-</td>\n",
       "      <td>100288069.0</td>\n",
       "      <td>chr1</td>\n",
       "      <td>661138</td>\n",
       "      <td>714014</td>\n",
       "      <td>LOC100288069</td>\n",
       "      <td>4</td>\n",
       "      <td>-</td>\n",
       "      <td>union_100288069</td>\n",
       "      <td>ncRNA</td>\n",
       "      <td>-</td>\n",
       "      <td>714014</td>\n",
       "      <td>714014</td>\n",
       "      <td>661138,665277,667396,668360,668686,670798,6742...</td>\n",
       "      <td>665184,665335,667587,668593,668744,670993,6744...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  chrom2  start2    end2 strand2       geneId   chr  txStart   txEnd  \\\n",
       "0   chr1   29369   29370       -     653635.0  chr1    14361   29370   \n",
       "1   chr1  714013  714014       -  100288069.0  chr1   661138  714014   \n",
       "\n",
       "       geneName  citationCount strand     union_geneId geneType  \\\n",
       "0        WASH7P              4      -     union_653635   pseudo   \n",
       "1  LOC100288069              4      -  union_100288069    ncRNA   \n",
       "\n",
       "                            geneDesc  cdsStart  cdsEnd  \\\n",
       "0  WASP family homolog 7, pseudogene     29370   29370   \n",
       "1                                  -    714014  714014   \n",
       "\n",
       "                                          exonStarts  \\\n",
       "0  14361,14969,15795,16606,16857,17232,17605,1791...   \n",
       "1  661138,665277,667396,668360,668686,670798,6742...   \n",
       "\n",
       "                                            exonEnds  \n",
       "0  14829,15038,15947,16765,17055,17368,17742,1806...  \n",
       "1  665184,665335,667587,668593,668744,670993,6744...  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of duplicated genes among the list of active ones 1865 - why?\n"
     ]
    }
   ],
   "source": [
    "# Find \"active genes\" corresponding to our active TSS - simply by overlapping TSS and with the genes\n",
    "active_genes_df = bioframe.overlap(\n",
    "    act_tss_df,\n",
    "    genes_df,\n",
    "    how=\"left\",\n",
    "    cols2=tx_cols,\n",
    "    suffixes=('2', '')\n",
    ")\n",
    "active_genes_df = active_genes_df.astype(\n",
    "    {'start2':int,\n",
    "     'end2':int,\n",
    "     'txStart':int,\n",
    "     'txEnd':int,\n",
    "     'citationCount':int,\n",
    "     'cdsStart':int,\n",
    "     'cdsEnd':int,\n",
    "    }\n",
    ")\n",
    "# active_genes_df - they are good, but there are plenty of duplicates - i.e. same transcripts\n",
    "# but it's ok - since we are going to drop_duplicates for TTS anyways ...\n",
    "display(active_genes_df.head(2))\n",
    "print(f'Number of duplicated genes among the list of active ones {active_genes_df.duplicated([\"chr\",\"txStart\",\"txEnd\"]).sum()} - why?')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "innovative-arrival",
   "metadata": {},
   "source": [
    "#### get to extracting TTS, corresponding to the active TSS (active genes) ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "owned-respondent",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " there are 1015 duplicated active pTTSs ...\n",
      "finally there are 6246 active +TTSs to work with...\n",
      " there are 924 duplicated active mTTSs ...\n",
      "finally there are 6091 active -TTSs to work with...\n",
      "there are still 2 duplicated active TTSs ... drop em\n",
      "finally there are 12335 active-TTSs to work with...\n"
     ]
    }
   ],
   "source": [
    "# columns needed for TSS/act_TTS bed files output:\n",
    "act_pgenes_df =  active_genes_df[ active_genes_df[\"strand\"] == \"+\" ][tx_cols].reset_index(drop=True)\n",
    "act_mgenes_df =  active_genes_df[ active_genes_df[\"strand\"] == \"-\" ][tx_cols].reset_index(drop=True)\n",
    "\n",
    "act_pTTS = act_pgenes_df.copy()\n",
    "act_mTTS = act_mgenes_df.copy()\n",
    "# after copying, act_pTTS and act_mTTS are dataframes that have chr, txStart, txEnd  of the\n",
    "# transcripts, that we could modify as we please ...\n",
    "\n",
    "# for genes on p-strand, we keep txEnd, and modify txStart to be ~several bp upstream:\n",
    "act_pTTS[\"txStart\"] = act_pTTS[\"txEnd\"] - 1\n",
    "act_pTTS[\"strand\"] = \"+\"\n",
    "# report how many act_TTS were duplicated here:\n",
    "print(f\" there are {act_pTTS.duplicated().sum()} duplicated active pTTSs ...\")\n",
    "act_pTTS = act_pTTS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(act_pTTS)} active +TTSs to work with...\")\n",
    "\n",
    "# for genes on m-strand, we keep txStart, and modify txEnd to be ~several bp upstream (relative to the gene):\n",
    "act_mTTS[\"txEnd\"] = act_mTTS[\"txStart\"] + 1\n",
    "act_mTTS[\"strand\"] = \"-\"\n",
    "# report how many act_TTS were duplicated here:\n",
    "print(f\" there are {act_mTTS.duplicated().sum()} duplicated active mTTSs ...\")\n",
    "act_mTTS = act_mTTS.drop_duplicates().reset_index(drop=True)\n",
    "print(f\"finally there are {len(act_mTTS)} active -TTSs to work with...\")\n",
    "\n",
    "# after modifying act_pTTS and act_mTTS separately - let's bring them back together ...\n",
    "act_TTS = pd.concat([act_pTTS,act_mTTS]).reset_index(drop=True)\n",
    "print(f\"there are still {act_TTS.duplicated(subset=tx_cols).sum()} duplicated active TTSs ... drop em\")\n",
    "act_TTS = act_TTS.drop_duplicates(subset=tx_cols).reset_index(drop=True)\n",
    "print(f\"finally there are {len(act_TTS)} active-TTSs to work with...\")\n",
    "\n",
    "# DROP IT !\n",
    "act_TTS.to_csv(\"activeTTS_HCT116-HAP1_wstrand.bed\",index=False,header=None,sep=\"\\t\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}