macroH2A_IPS_StemCells_2021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was implemented to\n",
    "- compute the *correlation networks* for the three examined experimental conditions: **epi** (Episome), **m11** (macroH2A1.1), and **m12** (macroH2A1.2);\n",
    "- summarize the obtained networks and the counts of genes belonging to **NHEJ**, **DDR not NHEJ**, **Reprogramming**, and **Senescence** pathways, which we represented as bar plots;\n",
    "- export the biggest m11 and m12 networks as SIF and attribute files for easy representation in Cytoscape (sizes **62** and **8**, respectively).\n",
    "- finally compute a one-way ANOVA among all conditions (m11, m12, episome) and select the significantly most varying genes to be represented by heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os.path\n",
    "import math\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from IPython.display import display\n",
    "\n",
    "from multiprocessing import Pool\n",
    "from functools import partial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Importing gene expression data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hidden": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Gene Symbol</th>\n",
       "      <th>gene_biotype</th>\n",
       "      <th>EPI</th>\n",
       "      <th>3_I</th>\n",
       "      <th>3_II</th>\n",
       "      <th>1_1_1</th>\n",
       "      <th>1_1_2</th>\n",
       "      <th>1_1_3</th>\n",
       "      <th>2_I</th>\n",
       "      <th>2_II</th>\n",
       "      <th>2_III</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TBC1D3L</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>7.148288</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>36.919539</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DNAH5</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>21.165576</td>\n",
       "      <td>11.080615</td>\n",
       "      <td>12.798294</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.395609</td>\n",
       "      <td>0.751568</td>\n",
       "      <td>4.215694</td>\n",
       "      <td>2.008978</td>\n",
       "      <td>2.188068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SAP25</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>7.201708</td>\n",
       "      <td>7.178799</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.870301</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>EYA1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>25.060854</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.741551</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.318718</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>CYP39A1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.528252</td>\n",
       "      <td>15.767867</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>15.017370</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Gene Symbol    gene_biotype        EPI        3_I       3_II  1_1_1  \\\n",
       "0      TBC1D3L  protein_coding   7.148288   0.000000  36.919539    0.0   \n",
       "4        DNAH5  protein_coding  21.165576  11.080615  12.798294    0.0   \n",
       "7        SAP25  protein_coding   0.000000   7.201708   7.178799    0.0   \n",
       "10        EYA1  protein_coding   0.000000  25.060854   0.000000    0.0   \n",
       "11     CYP39A1  protein_coding   0.000000   0.528252  15.767867    0.0   \n",
       "\n",
       "       1_1_2     1_1_3       2_I       2_II     2_III  \n",
       "0   0.000000  0.000000  0.000000   0.000000  0.000000  \n",
       "4   0.395609  0.751568  4.215694   2.008978  2.188068  \n",
       "7   0.000000  0.000000  8.870301   0.000000  0.000000  \n",
       "10  0.741551  0.000000  0.000000   0.318718  0.000000  \n",
       "11  0.000000  0.000000  0.000000  15.017370  0.000000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if os.path.exists(\"./playgrounds/data/NAR_2021/CTL_vs_macroH2A11_epi.xlsx\"):\n",
    "    m11_file = \"./playgrounds/data/NAR_2021/CTL_vs_macroH2A11_epi.xlsx\"\n",
    "    m11_file_huv = './playgrounds/data/NAR_2021/HUVEC_vs_macroH2A11_epi.xlsx'\n",
    "    m12_file = \"./playgrounds/data/NAR_2021/CTL_vs_macroH2A12_epi.xlsx\"\n",
    "    m11_vs_m12_file = \"./playgrounds/data/NAR_2021/macroH2A11_epi_vs_macroH2A12_epi.xlsx\"\n",
    "    gene_lists =  pd.read_excel(\"./playgrounds/data/NAR_2021/gene_lists_17032021.xls\", sheet_name=\"lists\")\n",
    "    \n",
    "    import macroH2A_IPS_NAR_2021_workers\n",
    "elif not os.path.exists(\"./data/NAR_2021/CTL_vs_macroH2A11_epi.xlsx\"):\n",
    "    ## Conditional branch for Google Colab\n",
    "    !git clone https://github.com/mazzalab/playgrounds.git\n",
    "    m11_file = \"./playgrounds/data/NAR_2021/CTL_vs_macroH2A11_epi.xlsx\"\n",
    "    m11_file_huv = './playgrounds/data/NAR_2021/HUVEC_vs_macroH2A11_epi.xlsx'\n",
    "    m12_file = \"./playgrounds/data/NAR_2021/CTL_vs_macroH2A12_epi.xlsx\"\n",
    "    m11_vs_m12_file = \"./playgrounds/data/NAR_2021/macroH2A11_epi_vs_macroH2A12_epi.xlsx\"\n",
    "    gene_lists =  pd.read_excel(\"./playgrounds/data/NAR_2021/gene_lists_17032021.xls\", sheet_name=\"lists\")\n",
    "    \n",
    "    ## Import the worker function for parallel computing\n",
    "    import playgrounds.macroH2A_IPS_NAR_2021_workers\n",
    "    ## Install Pyntacle\n",
    "    !pip install pyntacle==1.3.2\n",
    "else:\n",
    "    m11_file = \"data/NAR_2021/CTL_vs_macroH2A11_epi.xlsx\"\n",
    "    m11_file_huv = 'data/NAR_2021/HUVEC_vs_macroH2A11_epi.xlsx'\n",
    "    m12_file = \"data/NAR_2021/CTL_vs_macroH2A12_epi.xlsx\"\n",
    "    m11_vs_m12_file = \"data/NAR_2021/macroH2A11_epi_vs_macroH2A12_epi.xlsx\"\n",
    "    gene_lists =  pd.read_excel(\"data/NAR_2021/gene_lists_17032021.xls\", sheet_name=\"lists\")\n",
    "    \n",
    "    import macroH2A_IPS_NAR_2021_workers\n",
    "    \n",
    "m11_m12_epi = pd.read_excel(m11_file, sheet_name=\"CTL_vs_macroH2A11_epi\", \n",
    "                            index_col=None, na_values=['NA'], usecols = \"A,C,J:R\", engine='openpyxl')\n",
    "m11_m12_epi = m11_m12_epi.loc[m11_m12_epi['gene_biotype'] == \"protein_coding\"]\n",
    "\n",
    "m11_m12_huv = pd.read_excel(m11_file_huv, sheet_name=\"HUVEC_vs_macroH2A11_epi\", \n",
    "                            index_col=None, na_values=['NA'], usecols = \"A,C,M:U\", engine='openpyxl')\n",
    "m11_m12_huv = m11_m12_huv.loc[m11_m12_huv['gene_biotype'] == \"protein_coding\"]\n",
    "\n",
    "epi_vs_m11 = pd.read_excel(m11_file, sheet_name=\"CTL_vs_macroH2A11_epi\", \n",
    "                           index_col=None, na_values=['NA'], usecols = \"A,E,H\", engine='openpyxl')\n",
    "epi_vs_m12 = pd.read_excel(m12_file, sheet_name=\"CTL_vs_macroH2A12_epi\", \n",
    "                           index_col=None, na_values=['NA'], usecols = \"A,E,H\", engine='openpyxl')\n",
    "m11_vs_m12 = pd.read_excel(m11_vs_m12_file, sheet_name=\"macroH2A11+epi vs macroH2A12+ep\", \n",
    "                           index_col=None, na_values=['NA'], usecols = \"A,E,H\", engine='openpyxl')\n",
    "\n",
    "display(m11_m12_epi.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Analysis of relationships"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## macroH2A1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "temp = m11_m12_epi.iloc[:,5:8]\n",
    "temp.index = m11_m12_epi.iloc[:,0]\n",
    "\n",
    "# remove genes with more than one zero\n",
    "temp = temp[(temp == 0).sum(1) < 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "num_processors = 5\n",
    "\n",
    "#Create a pool of processors\n",
    "p=Pool(processes = num_processors)\n",
    "func = partial(macroH2A_IPS_NAR_2021_workers.range_calc_Cor, temp)\n",
    "\n",
    "chunklen = math.ceil(temp.shape[0] / num_processors)\n",
    "chunks = [range(0, temp.shape[0])[i * chunklen:(i + 1) * chunklen] for i in range(num_processors)]\n",
    "\n",
    "#get them to work in parallel - highly intensive !!\n",
    "output = p.map(func, [chunk for chunk in chunks])\n",
    "# output = p.map(func, [i for i in range(0, temp.shape[0])])\n",
    "p.close()\n",
    "\n",
    "m11_cor = np.vstack(output)\n",
    "print(m11_cor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "m11_cor_df = pd.DataFrame (m11_cor)\n",
    "m11_cor_df.columns = temp.index\n",
    "m11_cor_df = m11_cor_df.set_index(temp.index)\n",
    "m11_cor_df.to_csv(\"m11_cor_df.txt\", sep=\"\\t\", header=True, index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## macroH2A1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "temp = m11_m12_huv.iloc[:,8:11]\n",
    "temp.index = m11_m12_huv.iloc[:,0]\n",
    "\n",
    "# remove genes more than one zero\n",
    "temp = temp[(temp == 0).sum(1) < 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 -0.9975208980091912 -0.9976873881590403]\n",
      " [0 0 0 ... 0 0 0]\n",
      " ...\n",
      " [0 0 0 ... 0 0.9998810590275212 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]]\n"
     ]
    }
   ],
   "source": [
    "num_processors = 5\n",
    "\n",
    "#Create a pool of processors\n",
    "p=Pool(processes = num_processors)\n",
    "func = partial(macroH2A_IPS_NAR_2021_workers.range_calc_Cor, temp)\n",
    "\n",
    "chunklen = math.ceil(temp.shape[0] / num_processors)\n",
    "chunks = [range(0, temp.shape[0])[i * chunklen:(i + 1) * chunklen] for i in range(num_processors)]\n",
    "\n",
    "#get them to work in parallel - highly intensive !!\n",
    "output = p.map(func, [chunk for chunk in chunks])\n",
    "p.close()\n",
    "\n",
    "m12_cor = np.vstack(output)\n",
    "print(m12_cor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "m12_cor_df = pd.DataFrame (m12_cor)\n",
    "m12_cor_df.columns = temp.index\n",
    "m12_cor_df = m12_cor_df.set_index(temp.index)\n",
    "m12_cor_df.to_csv(\"m12_cor_df.txt\", sep=\"\\t\", header=True, index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Episome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>EPI</th>\n",
       "      <th>3_I</th>\n",
       "      <th>3_II</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gene Symbol</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TBC1D3L</th>\n",
       "      <td>7.148288</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>36.919539</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DNAH5</th>\n",
       "      <td>21.165576</td>\n",
       "      <td>11.080615</td>\n",
       "      <td>12.798294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAP25</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>7.201708</td>\n",
       "      <td>7.178799</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CYP39A1</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.528252</td>\n",
       "      <td>15.767867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>KCNMA1</th>\n",
       "      <td>2.125832</td>\n",
       "      <td>6.211489</td>\n",
       "      <td>12.575072</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   EPI        3_I       3_II\n",
       "Gene Symbol                                 \n",
       "TBC1D3L       7.148288   0.000000  36.919539\n",
       "DNAH5        21.165576  11.080615  12.798294\n",
       "SAP25         0.000000   7.201708   7.178799\n",
       "CYP39A1       0.000000   0.528252  15.767867\n",
       "KCNMA1        2.125832   6.211489  12.575072"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp = m11_m12_epi.iloc[:,2:5]\n",
    "temp.index = m11_m12_epi.iloc[:,0]\n",
    "\n",
    "# remove genes that contain more than one zero\n",
    "temp = temp[(temp == 0).sum(1) < 2]\n",
    "temp.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " ...\n",
      " [0 0 0 ... 0 (-0.9979656578155047, 0.04061446866363006) 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]]\n"
     ]
    }
   ],
   "source": [
    "num_processors = 5\n",
    "\n",
    "#Create a pool of processors\n",
    "p=Pool(processes = num_processors)\n",
    "func = partial(macroH2A_IPS_NAR_2021_workers.range_calc_Cor, temp)\n",
    "\n",
    "chunklen = math.ceil(temp.shape[0] / num_processors)\n",
    "chunks = [range(0, temp.shape[0])[i * chunklen:(i + 1) * chunklen] for i in range(num_processors)]\n",
    "\n",
    "#get them to work in parallel - highly intensive !!\n",
    "output = p.map(func, [chunk for chunk in chunks])\n",
    "p.close()\n",
    "\n",
    "epi_cor = np.vstack(output)\n",
    "print(epi_cor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "epi_cor_df = pd.DataFrame(epi_cor)\n",
    "epi_cor_df.columns = temp.index\n",
    "epi_cor_df = epi_cor_df.set_index(temp.index)\n",
    "epi_cor_df.to_csv(\"epi_cor_df.txt\", sep=\"\\t\", header=True, index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Format and filter datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "## Load from file if needed\n",
    "# m11_cor_df = pd.read_csv(\"m11_cor_df.txt\", sep=\"\\t\", dtype=object, index_col='Genes')\n",
    "# m12_cor_df = pd.read_csv(\"m12_cor_df.txt\", sep=\"\\t\", dtype=object, index_col='Genes')\n",
    "# epi_cor_df = pd.read_csv(\"epi_cor_df.txt\", sep=\"\\t\", dtype=object, index_col='Genes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Episome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "epi_cor_conn_df = epi_cor_df.reset_index().melt(id_vars='Genes').query('value != 0')\n",
    "epi_cor_conn_df = epi_cor_conn_df[epi_cor_conn_df['value']!='0']\n",
    "\n",
    "epi_cor_conn_df.set_index(epi_cor_conn_df['Genes'], inplace=True)\n",
    "epi_cor_conn_df['value2'] = [eval(x) for x in epi_cor_conn_df['value'].to_numpy()]\n",
    "epi_cor_conn_df[['b1', 'b2']] = pd.DataFrame(epi_cor_conn_df['value2'].tolist(), index=epi_cor_conn_df.index)\n",
    "epi_cor_conn_df = epi_cor_conn_df.drop(['value', 'value2'], axis=1)\n",
    "epi_cor_conn_df.columns = [\"Gene Symbol A\", \"Gene Symbol B\", \"Corr\", \"p_value\"]\n",
    "epi_cor_conn_df[['adj_p_value']] = macroH2A_IPS_NAR_2021_workers.fdr(epi_cor_conn_df['p_value'])\n",
    "\n",
    "epi_cor_conn_df = epi_cor_conn_df[(epi_cor_conn_df['adj_p_value'] <0.01)]\n",
    "\n",
    "epi_cor_conn_df.to_csv(\"epi_cor_conn_filtered.txt\", sep=\"\\t\", header=True, index=True)\n",
    "epi_cor_conn_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### macroH2A1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "m11_cor_conn_df = m11_cor_df.reset_index().melt(id_vars='Genes').query('value != 0')\n",
    "m11_cor_conn_df = m11_cor_conn_df[m11_cor_conn_df['value']!='0']\n",
    "\n",
    "m11_cor_conn_df.set_index(m11_cor_conn_df['Genes'], inplace=True)\n",
    "m11_cor_conn_df['value2'] = [eval(x) for x in m11_cor_conn_df['value'].to_numpy()]\n",
    "m11_cor_conn_df[['b1', 'b2']] = pd.DataFrame(m11_cor_conn_df['value2'].tolist(), index=m11_cor_conn_df.index)\n",
    "m11_cor_conn_df = m11_cor_conn_df.drop(['value', 'value2'], axis=1)\n",
    "m11_cor_conn_df.columns = [\"Gene Symbol A\", \"Gene Symbol B\", \"Corr\", \"p_value\"]\n",
    "m11_cor_conn_df[['adj_p_value']] = macroH2A_IPS_NAR_2021_workers.fdr(m11_cor_conn_df['p_value'])\n",
    "\n",
    "m11_cor_conn_df = m11_cor_conn_df[(m11_cor_conn_df['adj_p_value'] <0.01)]\n",
    "\n",
    "m11_cor_conn_df.to_csv(\"m11_cor_conn_filtered.txt\", sep=\"\\t\", header=True, index=True)\n",
    "m11_cor_conn_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### macroH2A1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "m12_cor_conn_df = m12_cor_df.reset_index().melt(id_vars='Genes').query('value != 0')\n",
    "m12_cor_conn_df = m12_cor_conn_df[m12_cor_conn_df['value']!='0']\n",
    "\n",
    "m12_cor_conn_df.set_index(m12_cor_conn_df['Genes'], inplace=True)\n",
    "m12_cor_conn_df['value2'] = [eval(x) for x in m12_cor_conn_df['value'].to_numpy()]\n",
    "m12_cor_conn_df[['b1', 'b2']] = pd.DataFrame(m12_cor_conn_df['value2'].tolist(), index=m12_cor_conn_df.index)\n",
    "m12_cor_conn_df = m12_cor_conn_df.drop(['value', 'value2'], axis=1)\n",
    "m12_cor_conn_df.columns = [\"Gene Symbol A\", \"Gene Symbol B\", \"Corr\", \"p_value\"]\n",
    "m12_cor_conn_df[['adj_p_value']] = macroH2A_IPS_NAR_2021_workers.fdr(m12_cor_conn_df['p_value'])\n",
    "\n",
    "m12_cor_conn_df = m12_cor_conn_df[(m12_cor_conn_df['adj_p_value'] <0.01)]\n",
    "\n",
    "m12_cor_conn_df.to_csv(\"m12_cor_conn_filtered.txt\", sep=\"\\t\", header=True, index=True)\n",
    "m12_cor_conn_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SIF from data/NAR_2021/epi_cor_conn_filtered_df.sif imported\n",
      "SIF from data/NAR_2021/m11_cor_conn_filtered_df.sif imported\n",
      "SIF from data/NAR_2021/m12_cor_conn_filtered_df.sif imported\n"
     ]
    }
   ],
   "source": [
    "from pyntacle.io_stream.importer import PyntacleImporter\n",
    "g_epi = PyntacleImporter.Sif(\"data/NAR_2021/epi_cor_conn_filtered_df.sif\")\n",
    "g_m11 = PyntacleImporter.Sif(\"data/NAR_2021/m11_cor_conn_filtered_df.sif\")\n",
    "g_m12 = PyntacleImporter.Sif(\"data/NAR_2021/m12_cor_conn_filtered_df.sif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select components that contain at least one gene of the lists above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "cc_epi = g_epi.clusters()\n",
    "cc_epi = [c.vs[\"name\"] for c in cc_epi.subgraphs() if c.vcount() >2 and \n",
    "          (set(c.vs[\"name\"]).intersection(gene_lists['NHEJ'].tolist()) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Senescence'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Reprogramming'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['DDR_not_NHEJ'].tolist()))\n",
    "          )]\n",
    "\n",
    "cc_m11 = g_m11.clusters()\n",
    "cc_m11 = [c.vs[\"name\"] for c in cc_m11.subgraphs() if c.vcount() >2 and \n",
    "          (set(c.vs[\"name\"]).intersection(gene_lists['NHEJ'].tolist()) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Senescence'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Reprogramming'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['DDR_not_NHEJ'].tolist()))\n",
    "          )]\n",
    "\n",
    "cc_m12 = g_m12.clusters()\n",
    "cc_m12 = [c.vs[\"name\"] for c in cc_m12.subgraphs() if c.vcount() >2 and \n",
    "          (set(c.vs[\"name\"]).intersection(gene_lists['NHEJ'].tolist()) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Senescence'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['Reprogramming'].tolist())) or \n",
    "           (set(c.vs[\"name\"]).intersection(gene_lists['DDR_not_NHEJ'].tolist()))\n",
    "          )]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 10, 5, 3, 3, 3, 3, 3]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[len(c) for c in cc_m12]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams.update({'figure.figsize':(2.8,1.5), 'figure.dpi':300})\n",
    "plt.rcParams.update({'font.size': 7})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAIACAYAAAC7LEJGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAC4jAAAuIwF4pT92AABzd0lEQVR4nO3dd5wV1f3/8ddHkCooCCiIuhRBigqCxtgoGkQxIAaQCEoxYu/GGhUTvxpLflFUkKiADQugIpqIgix2I1giTaVLbyJtZVk4vz9m7nJ3uffu3t27OzPwfj4e87hz75yZ+cye3bv3c08Zc84hIiIiIiISVfsFHYCIiIiIiEhpKKkREREREZFIU1IjIiIiIiKRpqRGREREREQiTUmNiIiIiIhEmpIaERERERGJNCU1IiIiIiISaUpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpSmpERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiIiIiKRpqRGREREREQiTUmNiIiIiIhEmpIaERERERGJtIpBBxBVZlYDOB5o7y/tgKaA+UUaOecWF/NYBwOXAGcDLYFawA5gFTADeAV40znnMngJIiIiIiJ7BdPn5JIxs6+BNimKFCupMbNzgOeBg4so+hFwvnNuXXFjFBERERHZF6j7WclZ3PovQDZey0rxD2DWFniT3QnNNKAf8Fvgd8AtwGp/22nA22amOhMRERERiaOWmhIys2uBtXjdw+Y755yZZQMd/CJFttSY2RvAef7TB51ztyUoUwv4AjjKf+k859zEUl+AiIiIiMheQt/6l5Bzbphz7mXn3I+lGOtymv+YC9yX5Dw/A/9IsI+IiIiIiKCkJmg1/cd1zrktKcotiluvUYbxiIiIiIhEjpKaYP3gP9Yxs+opymUl2EdERERERNCUzkEbCQwDKgF3AncULmBmBwI3+U83AS9m4sRmVg+om+ZuB+BNX70Jb3KEn/C6zomIiIiIVAIOj3s+3Tn3S3mcWElNsJ4EWgBXALeb2QnAKGAhXgLRFrgRqA/8DPRyzq1Ocqx0XQnck6FjiYiIiIgU1gN4qzxOpKQmQM65XcCVZjYJuAFvGuczCxX7Ffg7MNw591M5hygiIiIiEnoaUxMwM2sI/BE4OUmRKv72/mZWodwCExERERGJCLXUBMjvbvY2UA/YDNwKvA4sBSoDxwHXA38A7gc6mlkP59yvGTj9cGBcmvscDYyPPXnzzTdp2rRpBkIRERERkaibP38+5513XvxL5dbLSElNQMxsf+A1vIRmO9DBOfd1XJFc4GPgYzP7B97Ymi7A3SSYUCBdzrk1wJo0Yy7wvGnTprRq1aq0oYiIiIjI3qncJpRS97PgdGX3VM0vFUpoCrsL2OqvX2qFswsRERERkX2YkprgxDdxfJmqoHNuGzDbf1oHr3VHRERERERQUhOkHXHr+xejfHyZHUlLiYiIiIjsY5TUBGdB3HrHVAXNrC7Q2n/6s7+IiIiIiAhKaoI0FW/GM4DzzKx7okL+hAIj2d1SM9E558ohPhERERGRSNDsZyVkZk2BUwu9fGjcei8zWxf3fItzLn86ZOfcZjO7G/gnXnL5hpm9ALzB7imd2wBXAsf4u/0MDM3gZYiIiIiIRJ6SmpI7FRidYvvDhZ4vIe4eLwDOuUfNrBpeorI/MMBfElkM9HHOLSlJsCIiIiIieyslNQFzzt1vZq8Bg/HG1jQDagJ5wFrgW+AtYKw/C5qIiIiIiMRRUlNCzrkxwJgMHWs+GbihpoiIhJ9zjl27dqHhkSISBmbGfvvtt8dN1qNGSY2IiEgZy8vLY/PmzWzatIlt29ToLiLhU61aNWrWrEmNGjWoWDF6KUL0IhYREYmQLVu2sGzZMrXMiEiobdu2jW3btrF69WoaNmzIAQccEHRIadGUziIiImVECY2IRI1zjmXLlrFly5agQ0mLkhoREZEykJeXp4RGRCIpltjk5eUFHUqxqfuZiIhIGdi8efMeCU3VqlWpXbs2VatWjfygXBHZOzjnyMnJYcOGDeTk5BR4ffPmzdSqVSvA6IpPSY2IiEgZ2LRpU4HnVatW5cgjj1QyIyKhs//++1OjRg2WLFlSILHZtGlTZJIadT8TERHJMOfcHrOc1a5dWwmNiISWmVG7du0Cr23bti0yXWiV1IiIiGTYrl279nitatWqAUQiIlJ8VapU2eO1RO9nYaSkRkREJMMSfbOpVhoRCbv99tszNVBLjYiIiIiISDlQUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpFYMOQERERHYbmj006BDKxdCOQ4MOQTJs6NCh3HvvvQCMHj2agQMHFrnP4sWLadSoEQAdOnQgOzt7jzI7d+7ko48+YurUqXz22WfMnTuXdevWUbFiRerVq0f79u254IIL6NmzJxUqVEg77i5duvD+++8DcMghh7Bs2TIqVizeR+ScnBy++eYbZs6cyYwZM5gxYwbz5s1j586dAEybNo2OHTumHVMyGzZsYMaMGQXOt3Tp0vztUZl+uSwoqSkhM6sBHA+095d2QFMgdiOCRs65xWkesykwAOgKHAkcCGwAVgCfApOdc29nIn4RERGRsFu7di0tW7Zk3bp1e2zLzc1l8eLFLF68mPHjx9O+fXteffVVGjduXOzjL126lKlTp+Y/X716Ne+88w49evQo1v4XXnghb775ZrHPV1pNmjRh48aN5Xa+KFFSU3IfAm0ycSAz2w/4G3AzUKnQ5kP95Xjg94CSGhEREdknbN++PT+hady4MT179uTkk0+mYcOGAHzzzTc8/vjjzJo1ixkzZtCpUydmzpxJnTp1inX80aNHs2vXLgCqV6/O1q1befbZZ4ud1MS3jFSrVo02bdqwevVqFixYkM5lFlv8+erUqUO7du345JNP2LJlS5mcL0qU1JRc/K2hfwG+Bo7GS0CKfxCzCsCLQF//pQXAs8AXwEbgECAL6AI0L03AIiIiIlFiZnTs2JG77rqLzp0777H9xBNPZODAgfTp04eJEyeydOlS7r77boYPH17ksZ1zjB49GoBTTjmF4447juHDh/Of//yHlStXUr9+/SKP0bNnT3r27Em7du1o0aIFFSpUYODAgWWW1Pz5z3+mRYsWtGvXjiOPPBKArKwsJTUoqSmNUcBaYAYw3znnzCybNJMa4DZ2JzTPAFc653YkKDfCzAq34oiIiIjstQ477DCmTZuWskylSpV4+umn+c9//kNubi6vvvoqTz75JGaWcr8pU6awZMkSAAYOHEibNm0YPnw4eXl5PP/889x6661FxjdgwIDiX0wG3HnnneV6vijR7Gcl5Jwb5px72Tn3oyvhqCwzawzc7T/9EBiSJKGJnTO3JOcRERHZV40ZMwYzw8zyB6GPHz+erl27Ur9+fapWrUrz5s25+eabWbNmTYF9FyxYwHXXXcfRRx9NtWrVqF27Nl27di0wBqOwnTt3MmXKFG655RZOP/10Dj30UCpVqsQBBxxA06ZN6devH1OmTCl2/Hl5ebz44ov07t2bRo0aUb16dWrWrEmLFi3o3r07zzzzTMLxJrFrjg3Wnzt3LldffTXNmzenRo0amNkeY0Hy8vIYPXo03bp1o0GDBlSuXJnatWvTrl077rjjDlasWFHsuMtb3bp1OeaYYwBvMP369euL3OfZZ58FvG5jffr0oX379rRu3RqAUaNGlV2wUibUUhOsq9g9huaOkiZHIiIiUrSdO3fSr18/xo4dW+D1H374gX/84x+MHz+e7OxssrKymDhxIv379y/QrScnJ4fJkyczefJkRowYweWXX77HOW666SYee+yxPV7fsWMHCxYsYMGCBYwdO5a+ffsyZswYKleunDTeb7/9lj59+vDDDz/ssW3evHnMmzePSZMm8fHHHzNmzJikx3nuuee47LLL2L59e9IyixYtonv37syaNavA67m5ufz888989dVXPProo4wcOZKLLroo6XGClJu7+7vfomZB27BhQ35Sd/7551OzZk0ABg0axE033cQPP/zAxx9/zKmnnlpm8UpmKakJ1oX+40/OuU9iL5pZHaAWsME5V/RXDSIiIlKku+++m08//ZRu3boxaNAgGjVqxNq1axk5ciRvvPEGS5Ys4dJLL+Xvf/87vXv3pn79+tx3332ceOKJVKhQgalTp3Lfffexbds2rrvuOs4880yaNm1a4Bx5eXk0aNCA7t27c9JJJ9GkSROqV6/O6tWrmTt3Lk888QQLFy7klVdeoU6dOjz++OMJY/3222855ZRT2Lp1K+BNd3zRRRfRsmVLKleuzIoVK/jiiy94/fXXU17zjBkzGDt2LLVr1+b666/n1FNPpXLlysyZM4esrCzA+4DfoUMHfvrpJwBOO+00rrjiCpo1a8bGjRt56623GDFiBDk5OVx88cVUqVKF3r17l7I2MmvVqlXMnTsXgPr161OrVq2U5V988cX8JG/QoEH5r/fv359bb72VvLw8nn32WSU1EaKkJiBm1ojd428+8cfL3AQMwZsYIFZuCfAq8FAmExwzqwfUTXO3Jpk6v4iISHn79NNPueeeexg6dGiB17t06UL37t15++23mTJlCueccw6tWrXigw8+KPDh+MQTT6Rp06b06dOH3NxcRo4cycMPP1zgWDfccAOPPvpowvucdO3alWuvvZb+/fvzyiuvMGLECG6++eb8Ad8xeXl59OrVKz+heeyxx7j22msLlDn++OM599xz+dvf/lbgPiWFzZ49myZNmvDxxx9z6KG7h/2ecMIJ+es333xzfkIzZMgQnnrqqQLjUc444wx69OjB2WefTW5uLkOGDKFLly4ceOCBSc+7fPnyPVp9EslUl7b77ruPvLw8APr161dk+Vj3siOPPJJOnTrlv16vXj26devGxIkTGTduHMOGDaNGjRoZiVHKlsbUBKd13Ppa4GPgfuISGt+RwC3At2bWNoPnvxKYleYyMYPnFxERKVdt27blnnvu2eN1M+O6667Lf75mzRpGjx6d8Nv+Xr165U8nnGgAe5MmTVLeuLFChQo8/vjjVKhQgZ07dya8x8nYsWOZP38+AIMHD94joSnsiCOOSLl9+PDhBRKaeGvXruWll14CoFGjRgwbNizhAPvOnTtz4403ArBx48Yix5z85S9/4ZhjjilyOeuss1Iepzhi3QEBDj74YG677baU5WfOnMm3334LeAP9C19vrOVm69atvPrqq6WOT8qHkprg1I5bvxQ4AZgP9MbrelYVOJHdicRhwNtmlm7rioiIiOB9g59sRqy2bXd/b9i6dWvatGmTsJyZ5W8rzrS9W7duZcmSJcyZM4dZs2Yxa9YsVq1aRe3a3seAr776ao994hOdW265pchzpNKgQQN+97vfJd0+bdq0/LEogwYNSjnG56qrrspfnzx5cqniypR58+bRt29fdu3ahZnx3HPPcfDBB6fcJzZBgJklnL2sW7du1KtXr0BZCT8lNcGJb8usAiwDfuucG++c2+ic+9U59yXQE3jZL9cAKHp+QREREdnD0UcfnXRbfKtMixYtUh4nVnbTpk0Jty9atIjrrruOJk2aUKNGDbKysmjVqlWBFoq1a9cCJJy5bObMmQAceuihNG9eulvUtWnTJuXUxv/73//y108++eSUx2rYsGF+q1CspSOZ0aNH45wrclm0aFEaV1PQkiVL6NKlCxs3bgTggQceoFu3bin3ycnJ4eWXvY9Vp59+Oo0bN96jTMWKFenfvz8An3/+OXPmzClxjFJ+lNQEJ6fQ83udc3u8s/kzot0IxKZ6ztSUI8PxusClsxTv9roiIiIhVL169aTb9ttv90eiatWqpTxOrGzsTvTxXnvtNVq0aMGwYcNYuHAhRU1sum3btj1eiyU8hx12WMp9iyPWIpRM/NTHybqoxYvdkLI4UyaXpWXLltG5c+f8sUB33XVXse4rM2HChPwkKDbddSLxkwdoeudoUFITnM2Fnr+drKBzbhXeTT4B6vmTDJSKc26Nc252OgtQNrfHFRER2QvMnz+fiy++mO3bt1OtWjVuv/12pk+fzvLly/n111/ZtWtXfgvF4YcfDpAy6Snq5pHFUdTUxvHnz8T5ysOKFSvo3LkzCxcuBOC2227jr3/9a7H2je9ONmjQoPz7+RReYve8AXj++efZsSPpbQQlJJTUBGdJ3PoOP3Epbvl6ZRCPiIiIlMKoUaPypwmeMGEC999/P6effnr+jSzjk4YNGzYkPU7dut7w2WXLlpVtwECdOnXy11euXFlk+ViZosatlJUVK1bQqVMnfvzxR8Abc/TAAw8Ua9+FCxcyffr0tM+5du1aJk2alPZ+Ur40pXNwZgEOMGA/M7Mibr4Z/1XLzjKNTERERNL23XffAV6Xr65duyYtN3v27PzpmhNp3749S5cuZdWqVXz//felHleTyrHHHpu//tlnn3HGGWckLbt8+fL86aOPO+64MospmZUrV9KpU6f8m5H++c9/5sEHHyz2/qNGjcpvmfrb3/62xz2GCtu5cyeDBw8mNzeXUaNGcf7555c8eClzSmoC4pzbamZf4s1wVgFoBCxMsctRcetl/9WNiIiIpCV2n5ScnBx27tyZtOvXsGHDUh6nZ8+e+TfVfOihh8p0Bq5OnTpRqVIlcnNzGT16NLfccguVKlVKWHb48OH565mYijkdhROam2++mYceeqjY++/cuZMxY8YAXkvY7bffXmTXPIBXXnmFt99+m3fffZfly5dnZJyTlA11PwvWK3HrfZIVMrPmQOwrkXnF6KomIiIi5SzWopKTk8OLL76YsMyLL77I008/nfI4ffv2pVmzZoDXuvD444+nLB8bLF8SdevWzb9Z5cKFC7nhhhsSlsvOzuaRRx4B4KCDDmLw4MElPme6Vq1aRadOnfj+++8BL6EpfNPTokyePJnly5cD3r2GipPQgFcX4CVFzz33XFrnTEd2dnb+eJ6OHTuW2Xn2ZmqpCdbTwJ+B+sDtZjbJH5Cfz8yqAs/gdVMDeKJ8QxQREZHiGDhwIMOGDcM5x5AhQ/jqq68455xzqFu3LkuWLOHll19m/PjxNGvWjA0bNuTPclZYxYoVee211zjllFPYunUr1157La+//joXXXQRLVu2pHLlyqxcuZIvv/yScePG0b59+/xWiJJ45JFHmDJlCj/99BPDhw9nzpw5XHHFFTRt2pSNGzcyadIknnzyyfzB8v/617848MADS3y+dKxfv57OnTvnJzS///3vGTBgALNmzUq5X6NGjQrMdhff2nXBBRcU+/w9evSgatWq5OTkMHr0aG6//fYCY6NWrVrFu+++W2Cf2I1TAd59910WL15cYHuqWdeK8s033/DNN98UeG3Lli3564V/Dw499NCUXSH3JkpqSsjMmgKnFno5fi7EXmYWP0XzFufc+PjCzrktZnY58AZQE/jUzP4fMBVvyudjgZuAVv4u2cDIjF2EiIiIZEybNm34xz/+wU033URubi7Dhg3bo6vZUUcdxaRJk1LeEBO8MSsfffQRvXr1YuHChWRnZ5OdnZ2wbPv27UsVd+3atZk+fTrdu3dn1qxZSc9VtWpVnnrqKXr37l2q86Xju+++Y+7cufnPJ02aVKxB+9OmTctv8Ygf6N+gQQNOO+20Yp//gAMO4JxzzmHChAnMnz+fDz/8kA4dOuRvnzdvXoHpnwtLNOanNEnNm2++yb333pt0e+FYOnTooKRGinQqMDrF9sLtokuA8YULOefeMrMBwAi8xGaovxT2DtDPOZdXkmBFRESk7N1www20b9+eRx99lE8++YT169dz0EEH0bRpU84//3yuvPLKlPfLide2bVvmzZvHCy+8wJtvvslXX33FunXrqFKlCg0aNKB58+acd9559OhR+tvINWrUiK+//prnn3+e8ePH8/XXX7N+/XqqV69Oo0aNOOuss7jmmmto0KBBqc9V3uKnZO7Tp0+BexIVR9++fZkwYQLgtfjEJzUSHlbUTaEkMTMbSOqkprAlzrmsFMc7ArgC6AYcAVQB1gCfA2Occ/8ucbAZYmat8GZtA2DWrFm0atUqxR4iIvumvLy8/ClnY4466igqVtR3iSISXqV975o9ezatW7eOf6l14aEVZUXvriXknBsDjMng8ZYCt/uLiIiIiIgUk2Y/ExERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpSmpERERERCTSlNSIiIiIiEikKakpITOrYWYdzOwmM3vZzH4ws11m5vwlq5THHx53LGdmAzMTuYiIiIjI3qVi0AFE2IdAm7I4sJl1Bi4vi2OLiIiIiOxt1FJTcha3/guQDawq9UHNDgCe9Y+/urTHExERERHZ2ympKblRwIVAM6CWc64T8H0GjvswkAV8AzyTgeOJiIiIiOzV1P2shJxzwzJ9TDM7A7gMyAMGAz0yfQ4REQm3oUODjqB87CvXuS8ZOnQo9957LwCjR49m4MCBRe6zePFiGjVqBECHDh3Izs7eo8zOnTv56KOPmDp1Kp999hlz585l3bp1VKxYkXr16tG+fXsuuOACevbsSYUKFdKOu0uXLrz//vsAHHLIISxbtoyKFYv3ETknJ4dvvvmGmTNnMmPGDGbMmMG8efPYuXMnANOmTaNjx45px5TMrFmzmDx5Mh9//DHfffcdK1euJC8vj9q1a3PsscfSrVs3BgwYwIEHHpixc0aFkpqQKNTt7GHn3NdmpqRGRERE9llr166lZcuWrFu3bo9tubm5LF68mMWLFzN+/Hjat2/Pq6++SuPGjYt9/KVLlzJ16tT856tXr+add96hR4/ifQS78MILefPNN4t9vtLo2LEj06dPT7ht1apVrFq1ivfee48HHniA5557ji5dupRLXGGh7mfh8QhwJF4XtnsDjkVEREQkcNu3b89PaBo3bsxNN93EhAkT+OKLL/jiiy8YOXIkrVu3BmDGjBl06tQpYQKUzOjRo9m1axcA1atXB+DZZ58t9v7Oufz1atWqcfLJJ9OkSZNi75+OZcuWAVCrVi0GDx7Mc889x0cffcTMmTMZN24c3bp1A7wEp3v37nz00UdlEkdYqaUmBMzsTLxuZw64xDm3vRzOWQ+om+ZuZfNXKiIiIpKAmdGxY0fuuusuOnfuvMf2E088kYEDB9KnTx8mTpzI0qVLufvuuxk+fHiRx3bOMXr0aABOOeUUjjvuOIYPH85//vMfVq5cSf369Ys8Rs+ePenZsyft2rWjRYsWVKhQgYEDB7JgwYL0L7YIzZo146677qJv375Urly5wLbjjz+eXr168eCDD3Lbbbexfft2Lr/8cmbPnp3xOMJKLTUBM7Ma7J4Q4Ann3CfldOorgVlpLhPLKTYRERERDjvsMKZNm5YwoYmpVKkSTz/9NJUqVQLg1VdfLdCCksyUKVNYsmQJAAMHDmTQoEEA5OXl8fzzzxcrvgEDBjBgwABat25dovE86fj3v//NgAED9kho4t166620bdsWgDlz5vC///2vTGMKEyU1wYt1O1sC3BFwLCIiInuVMWPGYGaYWf4g9PHjx9O1a1fq169P1apVad68OTfffDNr1qwpsO+CBQu47rrrOProo6lWrRq1a9ema9euBcZgFLZz506mTJnCLbfcwumnn86hhx5KpUqVOOCAA2jatCn9+vVjypQpxY4/Ly+PF198kd69e9OoUSOqV69OzZo1adGiBd27d+eZZ55J2N0qds2xwfpz587l6quvpnnz5tSoUQMz22MsSF5eHqNHj6Zbt240aNCAypUrU7t2bdq1a8cdd9zBihUrih13eatbty7HHHMMABs2bGD9+vVF7hPrZlatWjX69OlD+/bt87uyjRo1quyCLWPxCeAPP/wQYCTlS93PAmRmvwOG+E+HOOe2BBmPiIjI3mznzp3069ePsWPHFnj9hx9+4B//+Afjx48nOzubrKwsJk6cSP/+/dmyZfe/5pycHCZPnszkyZMZMWIEl1++532yb7rpJh577LE9Xt+xYwcLFixgwYIFjB07lr59+zJmzJiU37p/++239OnTJ+EH03nz5jFv3jwmTZrExx9/zJgxY5Ie57nnnuOyyy5j+/bkvdsXLVpE9+7dmTVrVoHXc3Nz+fnnn/nqq6949NFHGTlyJBdddFHS4wQpNzc3f72oVpMNGzbkJ3Xnn38+NWvWBGDQoEHcdNNN/PDDD3z88ceceuqpZRZvWUnn57A3UVITkELdzkY7594r5xCGA+PS3KcJ6oImIiIRdffdd/Ppp5/SrVs3Bg0aRKNGjVi7di0jR47kjTfeYMmSJVx66aX8/e9/p3fv3tSvX5/77ruPE088kQoVKjB16lTuu+8+tm3bxnXXXceZZ55J06ZNC5wjLy+PBg0a0L17d0466SSaNGlC9erVWb16NXPnzuWJJ55g4cKFvPLKK9SpU4fHH388Yazffvstp5xyClu3bgW86Y4vuugiWrZsSeXKlVmxYgVffPEFr7/+esprnjFjBmPHjqV27dpcf/31nHrqqVSuXJk5c+aQlZUFeB/wO3TowE8//QTAaaedxhVXXEGzZs3YuHEjb731FiNGjCAnJ4eLL76YKlWq0Lt371LWRmatWrWKuXPnAlC/fn1q1aqVsvyLL76Yn+TFup0B9O/fn1tvvZW8vDyeffbZSCY106ZNy19v1apVgJGULyU1wfkHcASwCripvE/unFsDrCmyYBwzK6NoREREyt6nn37KPffcw9BCN8np0qUL3bt35+2332bKlCmcc845tGrVig8++KDAh+MTTzyRpk2b0qdPH3Jzcxk5ciQPP/xwgWPdcMMNPProownvc9K1a1euvfZa+vfvzyuvvMKIESO4+eabOfLIIwuUy8vLo1evXvkJzWOPPca1115boMzxxx/Pueeey9/+9jeWLl2a9Jpnz55NkyZN+Pjjjzn00EPzXz/hhBPy12+++eb8hGbIkCE89dRTBf7nn3HGGfTo0YOzzz6b3NxchgwZQpcuXVLeC2X58uV7tPokkqkubffddx95eXkA9OvXr8jyse5lRx55JJ06dcp/vV69enTr1o2JEycybtw4hg0bRo0aNTISY3kYP358/s+9ffv2NGvWLOCIyo/G1ATAzOoDl/pPJwNnmVnfwgvQOm6338Rta1TuQYuIiERc27Ztueeee/Z43cy47rrr8p+vWbOG0aNHJ/y2v1evXjRs2BAo+I14TJMmTVLeuLFChQo8/vjjVKhQgZ07dya8x8nYsWOZP38+AIMHD94joSnsiCOOSLl9+PDhBRKaeGvXruWll14CoFGjRgwbNizhl5idO3fmxhtvBGDjxo1Fjjn5y1/+wjHHHFPkctZZZ6U8TnHEugMCHHzwwdx2220py8+cOZNvv/0W8Ab6F77eWMvN1q1befXVV0sdX3lZtmwZV199NeD9ThdOuPd2SmqCEd+BdgDwcpLlD3HlLo97vUP5hCkiIrL36NevX9JeB7EZowBat25NmzZtEpYzs/xtxZm2d+vWrSxZsoQ5c+Ywa9YsZs2axapVq6hduzYAX3311R77xCc6t9xyS5HnSKVBgwb87ne/S7p92rRp+WMwBg0alHKMz1VXXZW/Pnny5FLFlSnz5s2jb9++7Nq1CzPjueee4+CDD065T2yCADNjwIABe2zv1q0b9erVK1A27LZu3UqPHj1YvXo14LW+dezYMdigypmSGhEREdknHH300Um3xbfKtGjRIuVxYmU3bdqUcPuiRYu47rrraNKkCTVq1CArK4tWrVoVaKFYu3YtQMKZy2bOnAnAoYceSvPmzVNfVBHatGmTsvt4/JS/J598cspjNWzYML9VKNbSkczo0aNxzhW5LFq0KI2rKWjJkiV06dKFjRs3AvDAAw/k34AymZycHF5++WUATj/9dBo3brxHmYoVK9K/f38APv/8c+bMmVPiGMvD9u3b6dmzZ36CfO6553L//fcHHFX5U1ITAOfcYuecFbUA98btNihu25iAQhcREYms2B3jE9lvv90fiapVq5byOLGysTvRx3vttddo0aIFw4YNY+HChUXeL2Xbtm17vBZLeA477LCU+xZHrEUomfipj5N1UYsXuyFlcaZMLkvLli2jc+fO+WOB7rrrLm699dYi95swYUJ+EhSb7jqR+MkDwjy9c25uLr169eL9998H4He/+x3jxo1L2QVyb6WkRkRERCQD5s+fz8UXX8z27dupVq0at99+O9OnT2f58uX8+uuv7Nq1K7+F4vDDDwdImfRkYoKeoqb0jT9/VCYEWrFiBZ07d2bhwoUA3Hbbbfz1r38t1r7x3ckGDRqUfz+fwkvsnjcAzz//PDt27MjsRWTAjh076NOnD2+//TYAnTp1YuLEiVSpUiXgyIKx76VxGWJmTYHC8/zFf8XRy8zi25S3OOfGl31kIiIiEoRRo0blTxM8YcIEunbtmrTshg0bkm6rW7cuS5cuZdmyZRmPsbA6derkr69cuZKWLVumLL9y5UqAIsetlJUVK1bQqVMnfvzxR8Abc/TAAw8Ua9+FCxcyffr0tM+5du1aJk2axPnnn5/2vmVlx44dXHDBBUyc6N1po2PHjrz99ttUrVo14MiCo6Sm5E4FRqfYXnjKiSWAkhoREZG91HfffQd4Xb5SJTSzZ8/On645kfbt27N06VJWrVrF999/X+pxNakce+yx+eufffYZZ5xxRtKyy5cvz58++rjjjiuzmJJZuXIlnTp1yr8Z6Z///GcefPDBYu8/atSo/Japv/3tb3vcY6iwnTt3MnjwYHJzcxk1alRokpq8vDwuuOAC3njjDcC7h9E777xTZLfJvZ2SGhEREZEMiN0nJScnh507dybt+jVs2LCUx+nZs2f+TTUfeuihMp2Bq1OnTlSqVInc3FxGjx7NLbfcQqVKlRKWHT58eP56JqZiTkfhhObmm2/moYceKvb+O3fuZMyYMYDXEnb77bcX2TUP4JVXXuHtt9/m3XffZfny5RkZ51QaSmiS05iaEnLOjSnOYP+4JasE5xiqyQFERESiIdaikpOTw4svvpiwzIsvvsjTTz+d8jh9+/bNv2niqFGjePzxx1OWjw2WL4m6devm36xy4cKF3HDDDQnLZWdn88gjjwBw0EEHMXjw4BKfM12rVq2iU6dOfP/994CX0KR7D5bJkyezfPlywLvXUHESGvDqAryk6LnnnkvrnOnIzs7OH8+TbCrmvLw8/vjHP+YnvLGEJtUEGPsStdSIiIiIZMDAgQMZNmwYzjmGDBnCV199xTnnnEPdunVZsmQJL7/8MuPHj6dZs2Zs2LAhf5azwipWrMhrr73GKaecwtatW7n22mt5/fXXueiii2jZsiWVK1dm5cqVfPnll4wbN4727dvnt0KUxCOPPMKUKVP46aefGD58OHPmzOGKK66gadOmbNy4kUmTJvHkk0/mD5b/17/+xYEHHlji86Vj/fr1dO7cOT+h+f3vf8+AAQOYNWtWyv0aNWpU4MN+fGvXBRdcUOzz9+jRg6pVq5KTk8Po0aO5/fbbC0yosGrVKt59990C+8RunArw7rvvsnjx4gLbU826lkq/fv0YP94bydCkSRMefPDBIqfErlevXv49d/Z2SmpEREREMqBNmzb84x//4KabbiI3N5dhw4bt0dXsqKOOYtKkSSlviAnemJWPPvqIXr16sXDhQrKzs8nOzk5Ytn379qWKu3bt2kyfPp3u3bsza9aspOeqWrUqTz31FL179y7V+dLx3XffMXfu3PznkyZNYtKkSUXuN23atPwWj9hAf/BuRnraaacV+/wHHHAA55xzDhMmTGD+/Pl8+OGHdOiw+x7o8+bNKzD9c2GJxvyUNKl57bXX8tcXLFjASSedVOQ+99xzD0OHDi3R+aJGSY2IiEiI7COfP/ZaN9xwA+3bt+fRRx/lk08+Yf369Rx00EE0bdqU888/nyuvvLLY3YXatm3LvHnzeOGFF3jzzTf56quvWLduHVWqVKFBgwY0b96c8847jx49epQ67kaNGvH111/z/PPPM378eL7++mvWr19P9erVadSoEWeddRbXXHMNDRo0KPW5ylv8lMx9+vQpcE+i4ujbty8TJkwAvBaf+KRGwsOKuimUSIyZtQLy23tnzZpFq1atAoxIRCSc8vLy8qecjTnqqKP2yRviiUh0lPa9a/bs2bRu3Tr+pdbOudmZizA5TRQgIiIiIiKRpqRGREREREQiTe3gEglDs4dm9ngdM3s8EREREQmOWmpERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiIiIiKRFrqkxsxqBB2DiIiIiIhER+iSGmCFmY00s3ZBByIiIiIiIuEXxqSmOvAn4L9mNsPM/mRm1YMOSkREREREwimMSc2ngPlLW2AkXuvNcDNrE2RgIiIiIiISPqFLapxzpwItgWHAz3jJTQ3gMmCmmX1hZoPMrGqAYWJmNcysg5ndZGYvm9kPZrbLzJy/ZBXjGJXN7Gwze9jMppvZKjPLNbPN/vFeMLOzyuFyREREREQiq2LQASTinJsHXG9mtwK9gUuB0/zN7YFngP9nZi8C/3LOfRdAmB8CbUq6s5kdC3wE1EyweX/gKH/pb2bvAf2cc+tKej4RERERkb1V6Fpq4jnntjvnXnTOdQBaAI8CG/Babw4ErgS+MbNPzOwiM6tSjuFZ3PovQDawKo39a7I7ofkOuB/oDpwAnAzcACzxt3cB3jezyqWIV0RERERkrxTqpCaec+5759yNwGFAf7yWktjYm5OAMcByM/unmbUoh5BGARcCzYBazrlOwPdp7L8LeBM4wTl3rHPuTufcJOfcDOfcZ865R4Fjgf/65dsA12UqeBERERGRvUVkkpoY51yuc26sc64jcCqw0t9kQC3gWmCWmb1nZqeUYRzDnHMvO+d+dM65Euz/qXOup3NuRooym4AhcS/1LUmsIiIiIiJ7s8glNQD+AP2XgKnAoezuCraF3a03ZwAfmtkTZmaJjxR+zrlvgfX+02ZBxiIiIiIiEkaRSWrM7GB/prF5wAd4rRaV8RKYT4GLgDpAK+AJIMffdgVweSBBZ04l/3FnoFGIiIiIiIRQ6JMaM+tkZi8Dy4CH8ForDNgKPAUc55w71Tn3kt81ba5z7lq8mcNm+mWHJDl86JlZe7wprQFmBxmLiIiIiEgYhXJKZzOrAwwC/gQ0jb3sP36Hl8y84JzbkuwYzrmVZnYn8G7cMaLonrj1lzJ1UDOrB9RNc7cmmTq/iIjIvigrK4slS5bQoUMHsrOzgw5HZK8RuqTGzF4BzsO7V0sskckFxgMjnHOfpHG4hf5jtYwFWI7M7FLgXP/pD8DTGTz8lRRMmEREJAyGDg06gvKxr1yniJSL0CU1QJ+49UXASGBUCW88mQssxZs+OVLM7HTgSf/pdqCvcy43wJBEREREREIpjEmNA94GRgCTSzJdcv6BnFsKZGUornJjZm2Bt/Baq3YBg51zXwcblYiIiJTW4sWLgw5BZK8UxqQmyzn3U9BBBMXMjgHeBw7ES/Auc86NLYNTDQfGpblPE2BiGcQiIiIiIlJioUtq9vGEphXevXcOxktornTOPVMW53LOrQHWpLNPhG/3IyIiIiJ7sVBO6WxmR/hL5WKUrRIrXx6xlRUza4l3/53YjGRXO+eeCjAkERGRvcqyZcu48847OfHEE6lTpw6VKlXikEMO4Xe/+x3Dhw9n+/btCfcbOnQoZoaZsXjxYvLy8hgxYgSnnHIKdevWpWrVqjRv3pwbbriBlStXpowhKysLM6Njx45Jy8yYMYMhQ4bQqlUratSokR9n69atOf/883nyySdZtmxZ0v3z8vIYPXo03bp1o0GDBlSuXJnatWvTrl077rjjDlasWJEyxo4dO2JmZGVlAbBx40buvfdejjvuOA488EBq1arFb3/7W1544QV27So4bHnSpEmce+65NGzYkCpVqpCVlcWVV17JqlWrUp4zZsOGDdx///2cdtppHHLIIVSqVIk6depw6qmn8uCDD7Jp06ZiHUf2PaFrqTGzk4GP8e5Dk4U3SD6V6sBcoLKZ/cY5N7NsI8w8P6GZBtTzX7raOTc8wJBERET2Ko899hi33nrrHonLmjVrmDJlClOmTOHRRx/lrbfe4uijj056nE2bNtGhQwc+/fTTAq//8MMP/PDDD4wePZo33niDTp06lSjOoUOH8te//pXCQ4rXrFnDmjVrmD17Nm+88QZr165laIIZ5BYtWkT37t2ZNWtWgddzc3P5+eef+eqrr3j00UcZOXIkF110UZHxzJ07l3POOWePsUCff/45n3/+OdOmTePZZ58lLy+Pyy67jNGjRxcot2TJEkaMGMHEiRP56KOPaNy4cdJzvfrqqwwZMmSPxGX9+vV88sknfPLJJ/zzn/9kwoQJnHLKKUXGLvuW0CU1wPn+4zvOufVFFXbOrTezt4ALgF54N9yMDDNrgddCE0tornHOPZliFxEREUnDfffdx1133QXA4YcfzlVXXUWrVq1o0KABa9eu5d///jdPPfUUP/74I2eeeSYzZ87kkEMOSXisIUOG8MUXX9C1a1eGDBnCkUceyapVq3jppZcYO3Ysv/zyC+eeey4zZ85MmRwl8t5773HvvfcC0LBhQ6644grat29P3bp1yc3NZfHixXz55ZdMmjQp4f4bNmygQ4cO/PST15P/tNNO44orrqBZs2Zs3LiRt956ixEjRpCTk8PFF19MlSpV6N27d9J4tm3bRvfu3Vm7di233XYbZ511FjVq1OC7775j6NChLFmyhNGjR3PGGWcwc+bM/NahgQMH0rhxY9atW8fIkSN5/fXXWbFiBZdeeilTp05NeK4XXniBAQMG4Jzj4IMP5qqrrqJt27Ycfvjh/PLLL0ydOpUnnniC1atXc/bZZ/PFF1/QokWLtH6+sncLY1JzGt54kvfS2OddvKSmQ5lEVEbM7Gi8hCb2znmtc+6JAEMSERHZq3z22Wfcc493W7ZLLrmEESNGsP/++xcoc9ZZZ9G3b1/OOOMMli9fzl/+8heefjrxreG++OILbr75Zh5++OECr59zzjmcfPLJXH311Wzbto3LL7887ZtrvvSSd4/tatWq8fnnn3PYYYcV2P6b3/yGCy64gEceeYT16/f83vfmm2/OT2iGDBnCU089VWA87BlnnEGPHj04++yzyc3NZciQIXTp0oUDDzwwYTxr165l+/btfPrppxx77LH5r7dr145OnTrRsmVLtm3bxg033MC6deu455579mg96tKlC+eeey7vvPMOH3zwAbNmzaJ169YFyixevJjLLrsM5xxnn30248aNo3r16gXKdO7cmYEDB3LKKaewdu1arr32Wt5///0ifqKyLwnjmJrYX/DcNPb5odC+Zc7MmprZwPgFODSuSK9C23sV2r8xXkIT2+cpYJqZtS5iKfhOLCIiIkn93//9H7t27aJx48YJE5qY3/72t1x55ZWA12qQk5OTsFzz5s154IEHEm676qqr6Ny5MwDTp0/nm2++SSvW2Hic5s2b75HQFHbwwQcXeL527dr8pKhRo0YMGzYs4QQ/nTt35sYbbwS8sTKjRo1KeZ577723QEITc+SRR9KzZ8/8c7dp0yY/eSzsuuuuy19PlOg98sgj5OTkULNmTV566aU9EpqYo446irvvvhuAKVOmsGjRopSxy74ljElNbKD8tjT2iZWtl7JUZp0KjC60NI/b/nChbY8U2v90oH7c88uB74qxlFviJiIiEmVbtmxh8uTJAPTq1StpQhMTG7y/fft2ZsyYkbDM4MGDqVgxeUeXIUOG5K+/++67acUbS2Rmz57NJ598kta+06ZNIzfXu0f3oEGDqFw5+VxLV111Vf567OeTTP/+/ZNua9u2bf56v379ks6Sevzxx+evL1iwYI/tr7/+OgBdu3alVq1aKeOJn2Dh448/TllW9i1hTGp+8R8Td2ZNLFY2nURIRERE9mJfffUVeXl5ADz00EP5M5glW37/+9/n75tsFrOTTjop5Tnjt3/77bdpxXvJJZdgZuTm5nL66afTpUsX/vnPf/LFF1/kJyzJ/O9//8tfP/nkk1OWbdiwIUcccUSRMdapU4c6deok3R6fgKQa3xJfrvAkAEuXLs3/Wb/22mtF1tExxxyTv29RM83JviWMSU2sLTGdaUM6+49LMxxLUs65Mc45S2PJKuX+sWVxeV2jiIhIlK1Zk9bt2ArYti3x96TJJhCIOfTQ3T3R161bl9Y5Tz31VJ577jkOOuggdu3axfvvv8+NN97ISSedxEEHHcSZZ57J008/za+//rrHvvFjbOJjSKZ+/fp77FdYsm5gMfvtt/tjZLVq1YpVbufOnQW2lUUdyb4pjBMFTANOAC4zs8eccyknNjez+sBleJMLfFAO8YmIiEgExFppAO666y769OlT7H0bNmyY8PV0bkRdkptWX3TRRfTo0YPx48fz/vvv8+mnn7J06VJycnKYOnUqU6dO5f/+7/948803adOmTf5+8VNAR+lm2fF1NHjwYG644YZi71uvXnmOOpCwC2NS8zRwE3Ag8K6Z9XTOJRwJ5g+2f90vuxP4V7lFKSIiIqFWt27d/PVt27btMetWSaxatYpmzZql3B5TeDB/cdWsWZPBgwczePBgwLtp6Pvvv8/TTz/NZ599xpIlS+jevTs//vhj/tiZ+G5iK1eupGXLlinPEeu6VdIYMyW+jn755ZeM1JHsm0LX/cw5twD4f4ABxwBzzOwlM/uTmZ1jZmf76y8Bs/0yDnjUOfd9cJGLiIhImLRt2za/69O///3vjBzz888/L/b24447LiPnbNiwIYMGDeKTTz7hrLPOAuCnn34qMFA+foayzz77LOXxli9fztKlSzMaY0k1atSI2rVrA/DBBx/scXNUkeIKXVLjux14CS+xqQz0BUYCk4C3/fW+/jYDXgRuCSRSERERCaXatWvnz5Y1d+5cxo8fX+pjjh49ukCXqcKeeeaZ/PUuXbqU+nzxzKzAMdeuXZu/3qlTJypVqpQfY6qJBYYPH56/HkuSgrLffvtx3nnnAfDzzz/z+OOPBxqPRFcokxrn3C7n3EXAYGABXuKSaJkPDHDODXDxnUlFREREgKFDh+aPMbn00kuLnCp52bJlPPvss0m3z5s3j7/85S8Jtz311FNMmTIFgNNOO63AVMbFMXbsWDZv3px0+65duwpMwdykSZP89bp169KvXz8AFi5cmHRsSnZ2No884t1l4qCDDsrv4hakO+64I78b3Z133skbb7yRsvz69et54gndq1wKCuOYmnzOuTHAGDM7BmiHdw8bA9YAM5xzswIMT0RERELutNNO4/777+f2229n48aNdOjQgT59+tCjRw8aN27Mfvvtx7p16/juu++YPHky06ZN46STTuKSSy5JeLzf/OY3PPjgg3z33XdceumlHHnkkaxatYqxY8fy4osvAlC1alWeeuqptGO94447uOyyyzj77LPp0KEDLVq0oFatWmzdupX58+czatQoPvroI8CbOvqEE04osP8jjzzClClT+Omnnxg+fDhz5szhiiuuoGnTpmzcuJFJkybx5JNPsmPHDgD+9a9/ceCBB6YdZ6Y1adKEUaNG0b9/f3Jzczn//PPp1q0bvXv3pnnz5lSuXJkNGzYwe/ZsPvjgA959913q1KnD1VdfHXToEiKhTmpinHOxG0+KiIiIpOW2227jkEMO4brrrmPz5s28/PLLvPzyy0nLp/qgP3LkSK644gr+/e9/JxynU7NmTV5//fUiB+ons2XLFsaNG8e4ceOSljnhhBPyb1gZr3bt2kyfPp3u3bsza9YssrOzyc7O3qNcLOnq3bt3iWIsCxdeeCEHHXQQgwYNYs2aNbzzzju88847ScuHIRmTcIlEUiMiIiJSGoMGDaJnz548++yzTJ48mVmzZrFhwwbASwaaNm3KSSedlN9KksyBBx7I9OnTeeqpp3j55Zf5/vvv2bp1K4cffjjdunXjlltuoUGDBiWK8eOPP+aDDz5g2rRp/O9//2P16tWsXbuWChUqcMghh9CuXTt69+5N7969C9z7JV6jRo34+uuvef755xk/fjxff/0169evp3r16jRq1IizzjqLa665psQxlqVzzjmHxYsX89xzz/Hvf/+br7/+mnXr1rFz504OOuggGjduzAknnMBZZ50V+FggCR/TUBQpLjNrBeR3+Zs1axatWrUql3MPzR6a2eN1zOzxRETi5eXl8eOPPxZ47aijjqJiRX2XGEVDhw7l3nvvBWDRokVkZWUFG5BIGSnte9fs2bMLT8vd2jk3O3MRJhfqd1fzRvYdDTQFagIVitrHOfd8WcclIiIiIiLhEcqkxswq4E3RfDVwaBq7OkBJjYiIiIjIPiR0SY2f0LwNxCZitwDDERERERGRkAtdUgMMAWKjvzYAo4D/AuuBXUEFJSIiIiIi4RTGpKa//7gQOMU5tzrIYEREREREJNwSzwcYrJZ4Y2PuU0IjIiIiQRs6dCjOOZxzmvlMJKTCmNRU8h+/DTQKERERERGJhDAmNT/5j1UDjaIIZlbDzDqY2U1m9rKZ/WBmu8zM+UtWmsc7zcyeN7NFZpZjZuvM7HMzu8XMapTRZYiIiIiIRF4Yx9T8G7gO+C3wacCxpPIh0Ka0B/HvxfMYcE2hTVWAg4HfAFea2fnOua9Kez4RERERkb1NGFtqHgU2AzeYWe2AY0klfqrpX4BsYFUJjvMguxOahcDleInMWXgzvwEcCfzHzI4sUaQiIiIiInux0CU1zrmlwB+Bg4APzKxNoAElNwq4EGgG1HLOdQK+T+cAZnYccJP/9EfgBOfcSOfcf51z7znnLgGu9bfXAx7JTOgiIiIiInuP0HU/M7NY68Q84HhgppnNBuYA24rY3fmJQJlzzg3LwGFuZ3dieY1zbkOC8zxuZhcCJwG9zKy5cy6t5ElERMqX17O4IOdcAJGIiBTfrl173hIy0ftZGIUuqQEG4k3pjP9oQCt/KY5ySWpKy8yqAuf6TxcD76Uo/i+8pAagN3Bf2UUmIiKltd9+e3aEyMnJYf/99w8gGhGR4vn111/3eC3R+1kYhTVKi1sKPy9qiYr2QHV/fZpL/RXe1Lj1jmUWkYiIZISZUa1atQKvbdiwQa01IhJazjk2bCjYaahatWpqqSkp51xYE61Mi295mp2qoHNuqZltAQ7AuzlpqZlZPaBumrs1ycS5RUT2BTVr1mTbtt29pnNycliyZAm1a9emSpUqkfn2U0T2brt27eLXX39lw4YN5OTkFNhWs2bNgKJKX+iSmn1I/ExmS4tR/iegBVDfzPZ3zu0o5fmvBO4p5TFERCSJGjVqsHr16gKtMzk5OSxfvjzAqEREisfMqFEjOrdK1NdEwYn/LdlcjPLxZaLzGyYiso+qWLEiDRs2jEzXDRGRGDOjYcOGVKwYnfYPJTXBqRq3nluM8tvj1qslLSUiIqFxwAEHKLERkUiJJTQHHHBA0KGkJdTpl5kdAQwCTgXq4yUCZznn5seVORpoAGx1zn0RSKAlE99psVIxyleOWy9qauviGA6MS3OfJsDEDJxbRGSfccABB9C0aVM2b97Mpk2bCoyzEREJi2rVqlGzZk1q1KgRqRaamNBGbGZ3AHcDsfkvDW+K58IJQEtgPLDDzBo659aWX5Slkm53svh0uTjd1VJyzq0B1qSzj75pFBEpmYoVK1KrVi1q1aqFc45du3ZpJjQRCQUzY7/99ov857xQJjVm9lfgTrxEJhf4DmiXpPgbwAq8lpyeePd0iYL4yQGOKEb5WJlVGZgkQEREAmJmVKhQIegwRET2KqEbU2Nmx+AlNOB1j2rgnDshWXn//i6v4yVAncs+woyZFbee8saifje8WEtNyumfRURERET2NaFLaoCr8BKUb4A/Ouc2pC4OwOf+Y+uyCqoMzAC2+uudLHWb3xlx69PLLiQRERERkegJY1LTAW/szJPOuV3F3Gex/9iwTCIqA865HOAd/2kW0CVF8SFx6+kO7hcRERER2auFMamJJSb/S2OfWItH1KY6/jteAgfwuJnVLlzAzK4BTvKfvuGcm1dewYmIiIiIREEYJwqIJVr7pyxVUC3/cUuGY0nKzJriTTUd79C49V5mti7u+Rbn3Pj4ws65r83sn8CNwFHAl2b2d7yudwcCfwQG+8XXATdl7gpERERERPYOYUxq1uDN9NUI+KyY+/zGf1xeJhEldiowOsX2hws9X4I39XRhfwaqAFcCjUk8e9tPwPnOuUUliFNEREREZK8Wxu5nsRtonlucwmZWCbgUrxvXJ2UVVFlxzu1yzl2FN5boJbzkZzuwAfgvcDvQ2jk3I7goRURERETCK4wtNeOAPkBvM3vYOfd1soJmVgF4Bq+Fw+ElBeXCOTcGGJPB430IfJip44mIiIiI7CtC11LjnJuAN91xBeA9MxtoZvETADgzO8DMeuG16vTDS2imOOc+Kv+IRUREREQkSGFsqQH4A954mgbAs8DT7J4lbDpQG+9eNviPi4GLyjdEEREREREJg9C11AA4534C2gHv4iUtFdidxNTBi9v8ZQpwknNuTQChioiIiIhIwMLaUoNzbjVwjpm1w2u5+Q1QDy/mtcBMYIJz7uPgohQRERERkaCFNqmJcc7NxEtgRERERERE9hDK7mciIiIiIiLFpaRGREREREQiTUmNiIiIiIhEWujG1JjZwlLs7pxzTTIWjIiIiIiIhF7okhogq5jlYvetsQSviYiIiIjIPiKMSc2HFJ2cVAeaATX9sj8Aq8o4LhERERERCaHQJTXOuY7FKWdmBpwHDAMOAvr70z+LiIiIiMg+JLITBTjPG8DpQGVgopkdHHBYIiIiIiJSziKb1MQ45xYBjwMNgOuDjUZERERERMpb5JMa3wf+Y89AoxARERERkXK3tyQ1v/iPRwYahYiIiIiIlLu9Jak5yn+M7JTOZtbZzMaY2Twz22xmO8xsvZl9YmZ/NbPDg45RRERERCSMQjf7WbrMrApwq//0hyBjKQkzqwq8CJyfYHNt4GR/ucnMrnHOjSrP+EREREREwi50SY2ZHVGMYvvhfeA/AbgWaIHXSvNaGYZWVp5nd0LzC/Ao8DmwHu9GpP2AHkA14BkzW+Oce7v8wxQRERERCafQJTXAYkrWjexrvHvWRIaZHQf08p+uBdo5536KK/IlMM7Mrsab4c2AvwFKakREREREfGEdU2NpLDuAUcAZzrlfA4m25E6LW3+6UEIT70lghb/exswOKNuwRERERESiI4wtNfcWo8wuYDOwAPjEObehbEMqMzXj1hcnK+Scc2a2BO9ePAA1gC1lGJeIiIiISGSELqlxzhUnqdlbxE9skJWskJkZu6er/gVYU4YxiYiIiIhESuiSmn3M28By4DDgUjMb7pxbnqDc5exupXnKObeztCc2s3pA3TR3a1La84qIiIiIZJqSmgA55341sx7AW3hJy3dm9ijwBbABr3WmL/AHf5dXgXsydPorM3gsEREREZHAKKkJmHNuppm1xWuNuZ7EY4q+AB52zk0oz9hERERERKIgdEmNmV1cFsd1zj1fFsfNkG54rTG1kmw/HhhgZvOdc9+WX1giIiIiIuEXuqQGGEPJ7lOTisO7yWWo+BMAPAMM9l+aCjwE/BfYCtTHS3iGAr8HzjCzvs65SRk4/XBgXJr7NAEmZuDcIiIiIiIZE8akBrz7z+wL/sTuhOY1oK9zLj6hWwqMMLP/ADOAg4GXzKyxc25daU7snFtDmrOoeTmYiIiIiEi4hDGp6YQX133Ab/AGzL8GfA6swkt4DgFOAnrjfdD/HPgLUOpZwcrZ5XHrfy6U0ORzzi02s2F4421qAH8EHi+H+EREREREQi90SY1zbrqZjQdOxJvt6zLn3KYERZ83s1uBkcAFwJXOuV7lGGomtPIfVzvnlhZR9su49ZZlFI+IiIiISOTsF3QAhZlZX+B84HPn3B+TJDQAOOc2A/3wWmp6mtmF5RRmpuzwH/cvRtn4MjuSlhIRERER2ceELqnBG2PiKGb3Kr/L1jC8bmmDyjCusrDAf6xtZscUUbZzgv1ERERERPZ5YUxqYh/uv09jnx8K7RsVb8atP2lm1RIVMrNTgMv8pzuBt8s4LhERERGRyAjdmBrgIP8x2T1bEomVPTCzoZS5/4fXunQEcBrwrZk9iTd+Jjal87nApezufjbMOaeWGhERERERXxiTmlV4H/J7Ah8Uc5+ecftGhnNuk5mdAUwAjgWaAv9MVhx4Eri5nMITEREREYmEMHY/ex9vfMxlZtatqMJmdjZe1yzn7xspzrn5QHugLzAeWIjXSpMH/AzMxBsz1NY5d41zbldQsYqIiIiIhFEYW2oeBi4CKgFvmtlo4FlgpnMuD8DMKgLtgEvwum9VALb7+0aOc24H3vTVrwYdi4iIiIhI1IQuqXHO/WhmlwDP4SUrl/hLnpltxGuRqcXu2A3YBfzJOfdj+UcsIiIiIiJBCmP3M5xzY4EueF2xzF/2B+oC9fz12OuLgC7OuZeCiVZERERERIIUupaaGOfcNDNrBnQDuuINpD8YL5FZD3wLTAbe1jgTEREREZF9V2iTGgA/WZnkLyIiIiIiInsIZfczERERERGR4lJSIyIiIiIikRbq7mdm9ltgCHAqUB+oAhzrnJsTV+Zk4Ghgk3NufCCBSpnLHtMxo8cbml3o+dCMHl5EREREylEokxr/PjTD8aZyBm9yAPCmcy6sCvAMsMvMmjjnlpRDiCIiIiIiEhJh7X42Ai+hMWAVMCFZQefcB8ACv+x55RGciIiIiIiER+iSGjM7ld0tNP8EjnTO9S5it9fxkpoOZRmbiIiIiIiETxi7n13qP37onLupmPt86T+2LIN4REREREQkxELXUgOcgjd25l9p7LPMf2yQ+XBERERERCTMwpjU1Pcf56WxT47/WCXDsYiIiIiISMiFManZ5T/WSGOfev7jLxmORUREREREQi6MSc1K/zErjX1O9R+XZjYUEREREREJuzAmNZ/gzWT2x+IUNrOaeDfodMCHZRhXmTOzCmbW18zGmdlCM9tmZpvMbL6ZTTKzm8xM44ZEREREROKEcfazl4ABQBcz6+Kcey9ZQTM7ABgPHIKX1IwplwjLgJm1xYv/2ASbawBNgHOB9UT4OkVEREREMi10SY1zboqZvQt0Bd40s4eAl+OKVDKzlsA5wHV4M5454GXn3LflHnAGmNnJwH+AmsBO4BVgErAYr44Ox0t2zg8oRBERERGR0ApdUuO7EJgOHAPc5S/O3zazUFnzX7us3KLLIDM7GJiAl9CsBc51zv03QdFXgDvMrFJ5xiciIiIiEnZhHFODc24jcBLwFLADL3FJtOwCngE6OOe2BhJs6T0IHOqvX5gkocnnnMst+5BERERERKIjrC01OOdygCvN7G/A74Hf4E3dXBGvRWMmMNE5tziwIEvJzOoA/f2nU5xzU4KMR0REREQkikKX1JhZbKD8BufcMufcSuBf/rK3+QNQ2V/PHzdkZpXxxgrtBFapdUZEREREJLkwdj/7Bvga6BZwHOXh5Lj1T8zseDN7C9gMLASWAJvM7H0z6xpIhCIiIiIiIRe6lhogB6gCfBV0IOWgtf/o8G4g+hR71kll4EzgTDN7ErjGOecoJTOrB9RNc7cmpT2viIiIiEimhTGpWQk0IpytSJlW2390wAi8+njCX58P1ALOA+73y16F9/P5vwyc+0rgngwcR0REREQkUGFMaj7CS2raAV8EHEtZq+E/7ucvtznnHozbvhoYaWaf4/0sKgN3mtnTzrk15RtqsAZ+M6ZY5ca0GVimcYiIiIhI+ISxNeQpvJaL682sWtDBlLGcuPWlwCOJCvk3FY1NlFAV6FXGcYmIiIiIREboWmqcc1+Y2d3AfcB7ZjbQOTc/6LjKyOa49f8453amKPs2cI2//htgeCnPPRwYl+Y+TYCJpTyviIiIiEhGhS6p8RMa8GZBOxmY53e/+hrYgDfNcVLOub+WaYCZtQRo4a8vLUbZmHqlPbHffS2tLmxmVtrTioiIiIhkXOiSGmAoXvcz/Mf9gN/6S3FEKan5DohN1VyhiLLx21MmdiIiIiIi+5IwjqkBsLil8PNUS9RMi1s/qoiy8duXlUEsIiIiIiKRFMaWmkZBB1COpgDrgYOBbmZWzTm3LUnZC+LWpyUpIyIiIiKyzwksqTGzD/C6lw12zsWPF4l1PVtexMD5yHPO7TCzB/BmPasNDDOzSwvfXNPMugN9/afLgDfLNVARERERkRALsqWmI14CU73Q64uBXcCxwJzyDSkQjwPn402KcAnQ2MxGAAvYffPNy/G61+0E/uSc2x5MqCIiIiIi4ROG7meJxsJEcXxMiTjncs3sXGA80Bno5C+FbQYGOOcml2d8IiIiIiJhF+REAVv8x7oBxhAKzrmfnXNnAH2At4DlQC6wEfgSb0a3Js65NwILUkREREQkpIJsqVmA18XsEjP70Dm3K8BYQsE5N470b4gpIiIiIrJPCzKpmQQcB1wInGlmC/BaJ2JGm9nWNI/p/BYPERERERHZRwSZ1DwE9ARaAYcA9eK2GdA+zeMZu2dOExERERGRfURgSY1zbouZnQRcBXQBGgKVgSPxkpOVwI6g4hMRERERkWgIdPYz59xWvBabh2KvmVlsbE0X59y+MKWziIiIiIiUQpCzn4mIiIiIiJRaGO5TU9gg/3FZoFGIiIiIiEgkhC6pcc49F3QMIiIiIiISHep+JiIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpCbEzKyRmW0xM+cvi4OOSUREREQkbJTUhJSZGTAKqB50LCIiIiIiYaakJryuAjoCqwOOQ0REREQk1JTUhJCZNQb+7j+9KshYRERERETCTklNyBTqdva6c25CwCGJiIiIiISakprwuRroAPyMWmlERERERIqkpCZEzKwJ8ID/9Ebn3Kog4xERERERiQIlNSFRqNvZe865McFGJCIiIiISDRWDDkDyXQOcDmwBhpT1ycysHlA3zd2alEUsIiIiIiKloaQmBAp1O7vdObekHE57JXBPOZxHRERERKRMKakJmN/tbDRQDfgYeDLYiPYN2YuzCzwfmp2dsFxxDe04tFT7i4iIiEjJaUxN8K4FTgN+Bf7knHMBxyMiIiIiEilqqQmQmTUF7vef3uuc+74cTz8cGJfmPk2AiWUQi4iIiIhIiSmpCdaf8bqd/QIsM7O+RZSvHldmq3NuUklP7JxbA6xJZx+vp5yIiIiISLgoqQlWZf/xQOCFYpSvA7zsry8BSpzUiIiIiIjsLTSmRkREREREIk1JTYCccwOdc1bUErfLkrjXs4KKW0REREQkTJTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1TOkdAockCREREREQkjlpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpFYMOYF9nZk2ArsDpwLFAQ6AysBGYA7wHPOucWx1UjCIiIiIiYaakJkBmNgYYkGRzXaCDv9xqZlc5514sr9hERERERKJCSU2wGvqP24C3gA+A74HNwOFAb6AfUBN43sx2OOdeDSJQEREREZGwUlITrBXA9Xjdy7YU2vY18JaZvQO8DBjwhJlNcs5tK98wRURERETCSxMFBMg5d7Fz7rEECU18mVeAN/2ndYAzyyM2EREREZGoUFITDVPj1psFFoWIiIiISAgpqYmGSnHrOwOLQkREREQkhDSmJho6x63PzsQBzawe3gxr6WiSiXOLiIiIiGSSkpqQM7MTgbP9pyuA7Awd+krgngwdS0REREQkMEpqQszMagJj2N1N8A7nXG5wEe09Bn4zpsDzrDGLk5bNHtixTGMRERERkdLRmJqQMrMKeFM5t/Bfes0591yAIYmIiIiIhJJaakLIzAx4FjjHf2kGcEmGTzMcGJfmPk2AiRmOQ0RERESkVJTUhIyf0IwABvgvfQuclepeNiXhnFsDrEkztkyGICIiIiKSEep+Fj6PA5f567OAM51zGwKMR0REREQk1JTUhIiZPQ5c5T+dDXR2zq0LMCQRERERkdBTUhMSZjYMuNp/Gkto1gYYkoiIiIhIJCipCQEzewy4xn8aS2jSGu8iIiIiIrKvUlITMDP7J3Ct/3QOSmhERERERNKi2c8CZGYPAtf7T9fhtdbUM7N6KXb72Tm3vKxjExERERGJCiU1wbogbr0OMLUY+zwHDCyTaEREREREIkjdz0REREREJNLUUhMg51xW0DGIiIiIiESdWmpERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiIiIiKRpqRGREREREQiTUmNiIiIiIhEmqZ0lr3KwG/GBB3CHoYOjfbxRUREpPwMzR6auWN1zNyxwk4tNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpSmpCxMxqmNktZva5ma0zsxwzW2Rmz5vZ6UHHJyIiIiISRhWDDkA8ZnY88DpwZKFNWf5ykZk9DlznnHPlG52IiIiISHippSYEzCwL+A+7E5pRwFnAb4DLgYX+69cAfy/v+EREREREwkxJTTj8A6jnr1/rnLvEOfeec+6/zrmRwAnAj/72m83s2ECiFBEREREJISU1ATOzo4Hz/aefOeceL1zGObcBr5UGvDq7rZzCExEREREJPSU1wesdt/6vFOXeA5b46783s6plF5KIiIiISHQoqQlex7j1D5IV8icHiG0/AGhXhjGJiIiIiESGZj8LXiv/cYtzbmkRZefErbcEPi7pSc2sHlA3zd2Ojn8yf/78kp4+bct//blMj+9+3pp025pFa4rcf3bd2cn3L3r3Upmd/NQiIiISMcX53FFcqT6flIUEnw0rlde5TbMDB8fMKgO/+k/nOOdaFVG+D/Cq//R+59ydpTj3UOCeku4vIiIiIlKEHs65t8rjROp+Fqwaceubi1E+vkyNpKVERERERPYhSmqCFT/YP7cY5bfHrVfLcCwiIiIiIpGkMTXByolbL06fw8px69tKee7hwLg09zkAaA9sAn4BfqJ4yVhpNAEmxj3vASwo43OWlb3lWnQd4aLrCJe95Tpg77kWXUe46DrCJdPXUQk4PO759FIcKy1KaoKVbneyA5Lsmzbn3BqgJCPRvijNedNlZoVfWuCci+TQ+L3lWnQd4aLrCJe95Tpg77kWXUe46DrCpYyu4+tS7l8i6n4WIOfcdnYnFoenKus7Mm69qJnSRERERET2CUpqgjfLf6xhZkcUUbZl3Hrkvg0QERERESkLSmqCF9/XsHOyQua1D3byn24FZpZlUCIiIiIiUaGkJnjj49aHpCjXBcjy1992zuWkKCsiIiIiss9QUhMw59wc4E3/6W/N7OrCZcysNvC4/3QX8PfyiU5EREREJPw0+1k43AicCtQBHjeztsDLeNMmtwFuAxr7Zf/pnPsmgBhFREREREJJSU0IOOcWmdnZwOt4s6AN9pfChgN/Ls/YRERERETCTklNSDjnZphZa+AK4HygKVAdWAV8DDztnCu3GxiJiIiIiESFkpoQcc5tAh70FxERERERKQZNFCAiIiIiIpGmpEZERERERCJN3c8k7NYC9xZ6HlV7y7XoOsJF1xEue8t1wN5zLbqOcNF1hMvech2Ycy7oGEREREREREpM3c9ERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiIiIiKRpqRGREREREQiTUmNiIiIiIhEmpIaERERERGJNCU1IiIiIiISaUpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRFrFoAMQiTGz/YGzgXb+0hg4GKgFbAdWADOAV4FJzjkXUKjFZmadgYuBk4DDgCrAJmAeMBV42jn3U0Cx1QCOB9r7SzugKWB+kUbOucVFHKMZ0Nk/znFAPbw6qwJsBObiXeco59yyjF/E7jhKfS1xxzoYuATvd7El3u/fDmAV3u/fK8CbZfH7Z2ZNgK7A6cCxQEOgMt7Pcg7wHvCsc251imMEXieZuI64YwVSH5l6Pwq6PjL9vhrk30cRcTUCvgOq+y8tcc5lJSkb+N9IMulcR9w+Qb5nLQaOLGbxG5xzjyY4RuD1kYnriDtWKP5GzKwC0Bv4A97f/qFAHrAG72eaDbzsnFtRaL/A66PEnHNatIRiwfvg44q5fALUDzrmFNdSFZhQjOvYCgwOKMavi4gtqxjHeKaY9bUVuCLM1+If5xxgXTGu50OgToavYUwxf5a/AP3DWieZuo4Q1EdG3o9CUB8Ze18Nsj6KuEYDphWKYXFY/0YydR1hqBNgcRq/X9eHtT4ycR1hqI+4ONoC3xYjjoFhrI+SLmqpkbDZAEzH+zZjEd43Gz8DB+J9Y/An4BjgZGCqmR3vnPs1oFhTeR4431//BXgU+BxYD2QB/YAeQDXgGTNb45x7u5xjtLj1X/ASg6Pxvs0prp3+fp8C/wNW4tWZAw7H+7Z+AN51Djezrc6550sf+h5KfS1m1hZ4E9jff2ka3pv7QuAAvH8SNwGHAKcBb5vZyc65XaUN3tfQf9wGvAV8AHwPbMb7WfbG+72pCTxvZjucc68mOE7QdZKR6whBfUBm3o+Cro+MXEdI6iOZq4COwGr//EUJQ50kktZ1hKxOZgCDiiizMsnrYaqPEl9HWOrDzE4G/oP3HrsTr2VoEl7iVhHvZ3osuz+jFBam+khP0FmVFi2xBW+M135FlKmI90Ep9k3BlUHHnSDG4+LiWwMcnqTc1XHlvg4gzmuBPwJHAea/lh0XU1YxjlGxGGVOwWt6d3hvjinrOMBreSOu/N+TlKkF/BBXrkcGr+F54DrggBRl+sadey1QLWx1ksHrCLo+MvJ+FIL6yNR1BFofKWJvDGzxz/eHuHMvTnW9QdZJBq8j8DphdwtHdimOEXh9ZOg6wlAfB/s/n9jnjxOLKF8pjPVR4usPOgAtWtJdgN/GvSGMCzqeBPHFJyv/l6KcAcvjyib9EFiOsWfHxZOVweNOjjtuqzBeC7u7DGxPVRfAZXHHfSSAOor/x9k9SnWSznVEqD4y8n4UgvpIeR1hrA//PTT2dz7Bf63IZCBsdVLS6whDnZCBZCAM9ZGJ6whJfcR3HTszqvVR0kWzn0kUbYpbrxlYFMnFx7Q4WSHnvSssiXupRlkFFAJhrzPYHdc659yWFOUWxa0HUWdT49ableI4QddJUdcRlfrI1M8x6Poo6vxhrI+rgQ54XemuKoPjl1edlPQ6wlgnZSnov5GiBFofZlYH6O8/neKcm5KpYycRuvpQUiNR1D9ufV5gUST3Q9x6VrJCZmbsnm3lF7ym4r2OmR0KnOk/3QX8GGA4qcTqrY6ZVU9RLivBPuWpUtz6zpIcICR1UtR1RKU+Sv1+FJL6KOo6QlUf/gx7D/hPb3TOrcrw8culTkp5HaGqk7IUkr+RogRdH3/Am2US4OXYi2ZW2cwamdkRZlYp8a7pCWt9KKmR0DOz/czsUDPraGYvAbf5m3KBpwIMLZm38bqVAVxqZoclKXc50MBff8o5V6IPqGFkZtXN7CgzuxZvgoSD/E0vOefWBRdZSiP9x0rAnYkKmNmBeAM9wfuW6sVyiKuwznHrs4u7UwjrpKjrCGV9ZOr9KOj6KMF1hKY+/C+ERuFNe/yec25Mho5brnWSgesITZ0AR5vZp2b2s5nlmtlqM8s2s7+YWb2SHDCgv5HSXEfQ9XFy3PonZna8mb2FN0nLQryeIZvM7H0z65ruwYN+zyqWoPu/adGSaAHqkHoqwY3AOUHHmSL+duweL7MBuBs4CzgB6AWMj7uWV4DKQcfsx50dF1dWmvs+UUSdvU45jhtK91rwvuQZHrfP+3iTD/wGOAO4Ge+eHrE67RxA/ZyI16rh/N+vPQZ5hrlO0rmOMNVHpt6Pgq6P0lxHyOrjWv88m4EjC21La0xNkHVS2usIQ51QvKmQtwFDwlwfmbiOoOsDmOkfexfefXJ2FHE9T+BPrBO2+ijxzyDoALRoSbSQ/J/vLuAfwCFBx1iMa6iHl8xsSHItnwN/CDrOQjFnx8WXlea+yd78vsdL6FK+eYblWvBumPZekmvJwesqknBGuzK+npp4N66MxTIganVSkusIQ31k6v0o6PrIxHUEXR9AE7z7Yzjg6gTbY7EsDnOdZPI6gqwTvG5H/wFuwOuO1BZv0onBFBxI7vBuWhnW+sjYdQRVH3hjdRzeF0a5/vrjeDcArYQ3lfRleLeWiMVzZxjro8Q/g6AD0KIl0QJUAFr7y3F4XVXuxGs+zcPr4tUw6DiLuIZBpL75VS7eNKrHBR1rXMzZcfFlpblv/bg6OxHvG6pxeB+YlgGXh/1a8O6x8jy7p1ZNtCwGbgcqlOO1VADeiYvh1YjWSVrXEZb6yNT7UdD1UdrrCLo+8GYJ+9A/z0ck+EAVH0dY/0YyeR0hqJNaRWy/0P/dcngtB83DVh+ZvI4g64M9b/p5a5JyxwG/+mW2AfXCVh8l/hkEHYAWLekseDewin1rsgpoEXRMCWI04Nm4N5YpQBe8/qf7A0cAV+DdYM3hfVv3+6Dj9mPPjos7K0PHPBdviksHDAvrteB1DYzVySbgFqAp3jdcNYBTKdhtcDJQpZx+n8bEnfdLStncH0SdpHsdYa2PQjFm5P0oqL+RdK4jDPWBd/8jh/dtd7IPlWklNUHUSaauIwx1Uszr/VtcDE+GrT4ydR1B1wfwU9yxl5AiaQKGxZVN+35/YaiPhHEFHYAWLekueN26tvl/TNlBx5Mgvkvj3ixeJUnzLN4MKLFvVjYBdUIQe3Zc7FkZPO4jccftGLZrwUs2F/llfwXapij7j7jj3l/G12B4g7Zj5/sGqB21Okn3OsJaH0nOn5H3oyD+Rop7HWGoD7wPh7HuWrelKJeRpKas6iRT1xGGOknjmuvHnX9emOojU9cRhvqgYLfep4oo2yWu7HNRrI9Ei2Y/k8hxzq0BPvafdjCz+kHGk8Dlcet/dv5ff2HOucV435aA9y3OH8s4riC9HrfeN7AokuvK7mk2X3LOfZ2i7F14H0rAm93OyjCux/H6QAPMwruZ2oYMHbs86yTd6whrfewhg+9Hgf6NFHEdYaiPPwPV8Ka/X2ZmfRMtceWrx73++xKesyzqJFPXEYY6KRbn3Eq8L/DA655VUkH/jaS6jjDUR/x975amUbZEs9MRwv/rFYMOQKSE4qcPzAJWBhRHIq38x9XOuaLeWL6MW29ZRvGEQXx9NQosiuRaxa1/mbQU4JzbZmaz8foX18H7h7A60wGZ2ePsvhHfbLyZcjI5bWa51EkJryN09VGETLwfheFvJNl1hKE+YvffOBB4oRjl67D7Xh1LgEklOGdZ1EmmriMMdZKOhF/upSkMfyPJriMM9fEdXnIF3vi5VOK3l/R2EmGojwLUUiNRFX/vl82BRZHYDv9x/2KUjS+zI2mp6AtzfUHBn33g9WZmw/DuMg67E4G1GT5NmddJKa4jVPVRDJn4WYbhbyRZDFGrj0wJQ50kE5k68W/UWMd/ujxV2SIEWh9FXEcY6mNa3PpRRZSN376shOcL3d+HkhqJHDPLwptqEbw+4AuCiyahWDy1zeyYIsrG34AwbNeRSRfGrf8vsCiSi//Zd0xV0Mzq4s0EA/Czv2SMmT0GXOM/jSUCazJ5Dl+Z1kkpryM09VGUDL4fBfo3UsR1BF4fzrmBzjkraonbZUnc61klPG3G6ySD1xF4naThSrxxdVDwg3e6gv4/kuo6wlAfU/CmawboZmbVUpS9IG69pHUSdH3sKehBPVq0xBagH3BwEWXqAv9l9+C0Z4KOO0GMQ+Pi+xColqTcKXiz3zi8qSKbhCD27LjYs4oo2ww4qxjHHMzuGy1uL6/rTPNaauBN1uD8WLsnKbc/Xj/i2HFHZzjmf8YdezYpptoMc51k4DoCr49MvB+FoT4ydB2B10ca11vUAPvA6yRD1xF4neDdj6V6EWXip0LOA1qHrT4ydB2B14d//Jvi/45JPF14d7wpmR3ejGmV47YFXh+lWcwPTiRwZpaNd+fdd/A+kM7B+wajIl4zZwdgAFDL32U+8FuX2XEGpWZmNfH6th7hvzQfeBKvn+1WvBlUzsWbJS3WBP1P59yN5RxnU7wpJuPdBjT31/9MwT6zW5xz4+P274j3Dc/3wJvADLwm+W14/cVb430bdHrcMW5wzj2aoUvIV9pr8Y9xPd6HcfDe8F8A3sAbcFkZaIP3TV2s9e1nvBlu4gdcluYaHsSbAhQ/1guAolo2fnbO5XeDCEOdZOI6/ONcT7D1kU0p349CUh+lvg7/ONcTYH0Ul5nFPtQscQlaaMJQJ8VR1HX4Za4n+L+RNnj3N/oI72e6EagCtMD7OZ4Vt8sdzrkHCh2jI+H4GynVdfjHuZ6A/0bMrBLez/Nk/6VpwAi8lqRawHl4kxlVxEtKujnnJsft35EI/H0kFXRWpUVLbKHgN+tFLe8AhwYdc4praUrqG2/Gll14M0PtF0CMA9P4ee/xjSFeE3tx910PDAzrtcQd5w5234k51bIIOCHD17A4zWtwwJiw1UkmriMk9ZGdRvwJ349CUh+lvo4w1Eca11vU33jgdZKJ6whDnaTxu7UFuCKs9ZGJ6whDfcTFUAuYWsT5NwE9w1gfpVnUUiOhYWaH4f1BdcD7RuMQvG4RFfCmv5wPfAG84pz7Ipgoi8/M9gfOB3oBx+NdT2W8AXULgU+AUc65bwOKbyAwOo1dCnxj6F/fiXj1dRpwJN4sLgfhfauzGq+f7WS8O8f/kom4EynttRQ6VlO8pvWOeE3xNfG6G6zFS1TfAsY657aVPOKE512M9zNMx3POuYFxxwi8TjJxHYWOF1R9lPr9KCT1kdH31aDqo7iK0VITeJ0UR3FaauLKBvU30g7ohNcSeDRwsL/sxPvA+x3eh+sxzrn1SY4ReH1k4joKHS8UfyNm1hvoD7TD+5vfBvwI/Ad4wiWYtCUM9VEaSmpERERERCTSNPuZiIiIiIhEmpIaERERERGJNCU1IiIiIiISaUpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpSmpERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiJSjsws28ycmS0OOpZkzKyjH6Mzs4FBxyMiUhQlNSIiIiIiEmlKakREpNzEffs/JuhYRERk71Ex6ABEREQkXJxz2YAFHYeISHGppUZERERERCJNSY2IiIiIiESakhoRkb2UmQ2MG8PS0X/tbDObaGYrzGy7/zjOzH5TzGOamfUys1fNbLGZbTOzLWY2z8xGmtmxSfZbbGYu7qUBcbHFL1l++c/859+niGVh3H63JylzelyZi5OUqWhmg8zsnbifywYzm2lm95tZgyJ+JgVmMzOz2mZ2t7//en/bo6mOkeCY9czsv3Gx/yWd/eOOc5iZ3ecf62cz2+Ff2/dm9r6Z3WlmbRLsl3T2s0Lbirt0TBJfiX6fREQK05gaEZF9g5nZU8BlhV6vD/QCzjezPznnRqc4wOHABOCEBJub+8ulZvZ/wN3OOZegXHFNBU4CmplZQ+fcskKxZAGN4l46A3ggwXE6FzpmAWbWCHgLaF1oUyWgFnA8cL2ZXeace6GooM3seGASkDIRKuIYRwHvAo2BHcClzrnnSnCcbsCrQPVCm2r5SzPgTOB3QMeSxltS5fz7JCJ7OSU1IiL7hnuB04APgGeAH4GqwLnAjXj/D4ab2YfOuQWFdzazQ4HPgMOAXcB4vA/vi/A+eB8LXAW0Af4CbAfuiztEF7xE4Tv/+US/XGHL/cepwJ3++hlA4Q/1ZxR6foqZVXHO/Zqk3PfOueXxG8ysNjAdONx/6SNgBPADcBDQHbgC7+f0vJn96pwblyDmmAOAN4GDgWHAO8A6oCFQIcV+8TGdhPdzrQNsBv7gnHu/OPsWOk4dYCxeQpODV+fvAav9Iofg1VVXIN1k4UvgmCLKNMD7HakB5AIrC8VX2t8nEZGCnHNatGjRomUvXICBeB9YY8vfk5S7NK7Mw0nK/Nvfvh5on6RMReA1v1wukJWgTOw8Y4qIvTLeh3EHvJBg+1h/21v+uRxwRqEy1eO2PZngGKPi4hkJWIIynfE+UDvgZ+DABGWy447zK3BKEdcWK7+40Os9gG3+thXAcaWo+0FxMfUsouzBCV7rGLf/wDTPXQP42t93F9CvrH6ftGjRoiW2aEyNiMi+4X/AHUm2jWL3N/gdC280s/bA2f7TW5xzMxIdxDmXh9eykQvsj/fBukScc9uBT/ynhVtlYHe3sonAf5OUO92PA7wWqnxmVhfo5z9dBFzrnNujxcI59wHw//ynBwGDiwj9/znnPimizB7M7ErgdbxWoTnASc65b9M9Tpz6cet7dLuL55xbX4rzFGBmsUSkjf/SX5xzLxUqU+6/TyKy91NSIyKyb3jBObcr0Qbn3E5gpv+0aYIif/Afd+GN0UjK/4Ac62J2agnijBf7MF7fzFrEXjSz1njdpwCmxJUrnNTEnu8CphXa1gmvOxzAaD+JSubJuPWzioh5TBHbC/AHyj/gn2M/4EPgVOfc0nSOk0B8V7shpTxWOkbgdWkDeMY5d3+CMkH9PonIXkxJjYjIvmFOEdtj39bXTLAtNjPafsDmoma6Atr55esnOFY64lsYzkiwvsA5tySuXDszOzBBuW+ccxsKHTt+Vq1PUwXhvEkKYknGcSmKbnXO/ZDqWIVUBF4AbvOfvwZ0cc79nMYxknkTWOuvP2xm35jZvWb2OzNLVMel5s/Q9if/6WS8VpZEgvp9EpG9mJIaEZF9w7YitsdacRL9X6hXwnNWK+F+MTOBX/z1M+NejyUrU/zHz/GurwJeC0xsoHwsAUnU/erguPVVxYglNtD94BRl0k1GDmN3F7g5wIVFtBgVm3PuF7wuXrEk6zjgbrzJAn72p5u+yx+wX2pm1h/4m//0W6C3330skaB+n0RkL6bZz0REpCix/xXrSW/q39zSnNQ5t9PMpuPNQtbRzGIziHXwH6f65XLN7CO8rmFn4LVSdALML1dgPI3P4tYzNVXwzjTLrwK+x7uelsDT/rTaCbsJpss5N9PMWuIlN+cCp/jn2Q9vqurjgdv86apfLOl5zKwT8Kz/dDnQzTm3OcUugfw+icjeTUmNiIgUZS3ePUNqA0udc5vK8dxT8ZKaA/G6Ie2H10XOUXCczBR2JzXEPe7Am6q5sHVx6/UpunterOtTxgbV482qdjZeEtYFbyB8ZTO72B/nVGr+cd72F8zsILwJFC4E+uC1fowxs2+cc7PSPb6fNL2ONz5pM3COKzR1dgJB/j6JyF5K3c9ERKQosUkEjN2DwMtLfNexMyk4TmZdgnItzKxBXLnPnXNbExz3f3Hrv00VgJkdBhzhPy3NjGR7cM7l4CVtk/yXLgReNbP9k+9VqvNtdM695Zzry+77AFUALkj3WGZWH/gP3qxweUAv59z/Uu7kCfL3SUT2UkpqRESkKBPi1m8r5Qfu2M0xKxensHNuNrunmz6D3clK4XEy37C7FWUQu2dxSzad8TR2d2caZGaVkpQDuDJufXIRIafNH0fzByB2Y88/AK+bWbF+RqUQfy1109nRzA7Aa/2JJXuXO+feK+bumfx9EhEBlNSIiEgRnHMfsXtQfltgdKoP3Ga2n5n9we+aVFisa1KzNEKIjYk52V+IiycWY3x3tD8n2LcA59xaIHb/lMbAPxOVM7OOwM3+04149/TJOOfcDuCPQGxsy7nAJDOrWpLjmdnZZnZEEcXOiVtfkMaxK+BNxXy8/9J9zrlnU+xSQIZ/n0REAI2pERGR4ukPfAEciTdj12/N7Gn/tZ+B6kAWcBLQE29mr9+x51iV6UAT4Hgz+yvezTN/idu+OMGsWVPxPvBX8Z/nkniczFSgF974G4CteDOjJXMzXpe2w4Er/Q/NI4D5eF2qfg9cxe4beA7xZxUrE/7ECAOAHOBSvJ/ff8zsXOfcljQPdwHQ359AYTJed7s1eF9mNsSro9jMaxvZnUwVx13sTog+xusu17qIfRYV6gaYqd8nERFASY2IiBSDc261mZ2Ed1+VM/FaNx5IsctOvKSisIeA3kANvA/HdxXa3ghYXOi1wl3IPnPOJZqiunC5j/wWkISccxvMrAPwFtAabyaujgmK5uB1rxqXYFtGOed2mdll/jmvxZsZ7T0zO7sECVUFkl9TzGq8sTArU5QpLCtu/VR23xwzlU5AduxJBn+fREQAdT8TEZFics6tcs79Du9D8tN435r/gveBczMwF29cyBXAYc65zxIc43u8WcxGAvPw7i+Tckpl59xiYGHcSwnHyTjnfmT3TTIhSdezQvsswusCdQneoPdVeDOmbQS+Bv4ONHXOPV/UsTLFea4DHvRf+i0w1cxqp3GY64AewKPAZ3iJ4ja8Vq5VeD/DG4BmzrmPMxN5ejLx+yQiEmNeN2QREREREZFoUkuNiIiIiIhEmpIaERERERGJNCU1IiIiIiISaUpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERiTQlNSIiIiIiEmlKakREREREJNKU1IiIiIiISKQpqRERERERkUhTUiMiIiIiIpGmpEZERERERCJNSY2IiIiIiESakhoREREREYk0JTUiIiIiIhJpSmpERERERCTSlNSIiIiIiEikKakREREREZFIU1IjIiIiIiKRpqRGREREREQiTUmNiIiIiIhEmpIaERERERGJNCU1IiIiIiISaUpqREREREQk0pTUiIiIiIhIpCmpERERERGRSFNSIyIiIiIikaakRkREREREIk1JjYiIiIiIRJqSGhERERERibT/D3wD2P7VPWeuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 840x450 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "kwargs_m11 = dict(alpha=0.5, bins=30)\n",
    "kwargs_m12 = dict(alpha=0.5, bins=3)\n",
    "kwargs_epi = dict(alpha=0.5, bins=3)\n",
    "\n",
    "temp_cc_m11 = [len(c) for c in cc_m11]\n",
    "temp_cc_m12 = [len(c) for c in cc_m12]\n",
    "temp_cc_epi = [len(c) for c in cc_epi]\n",
    "\n",
    "plt.hist(temp_cc_m11, **kwargs_m11, color='g', label='macroH2A1.1')\n",
    "plt.hist(temp_cc_m12, **kwargs_m12, color='b', label='macroH2A1.2')\n",
    "plt.hist(temp_cc_epi, **kwargs_epi, color='r', label='episome')\n",
    "\n",
    "plt.gca().set(ylabel='frequency', xlabel='network size')  # title='Frequency Histogram', \n",
    "plt.xticks(np.arange(min(min(temp_cc_m11), min(temp_cc_epi), min(temp_cc_m12)), \n",
    "                     max(max(temp_cc_m11), max(temp_cc_epi), max(temp_cc_m12))+2, \n",
    "                     5))\n",
    "plt.yticks(np.arange(0,20,2))\n",
    "\n",
    "# plt.xlim(50,75)\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig(\"FigureXX_A.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Count genes in gene lists for the clusters of sizes > 2 of mH2A1.1 co-expression network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cc_m11 = g_m11.clusters()\n",
    "cc_m11_ddr_not_nhej = [(set(c.vs[\"name\"]).intersection(gene_lists['DDR_not_NHEJ'].tolist()), len(c.vs[\"name\"])) \n",
    "               for c \n",
    "               in cc_m11.subgraphs() if c.vcount() >2 and set(c.vs[\"name\"]).intersection(gene_lists['DDR_not_NHEJ'].tolist())]\n",
    "\n",
    "cc_m11_nhej = [(set(c.vs[\"name\"]).intersection(gene_lists['NHEJ'].tolist()), len(c.vs[\"name\"])) \n",
    "               for c \n",
    "               in cc_m11.subgraphs() if c.vcount() >2 and set(c.vs[\"name\"]).intersection(gene_lists['NHEJ'].tolist())]\n",
    "\n",
    "cc_m11_reprogramming = [(set(c.vs[\"name\"]).intersection(gene_lists['Reprogramming'].tolist()), len(c.vs[\"name\"])) \n",
    "               for c \n",
    "               in cc_m11.subgraphs() if c.vcount() >2 and set(c.vs[\"name\"]).intersection(gene_lists['Reprogramming'].tolist())]\n",
    "\n",
    "cc_m11_senescence = [(set(c.vs[\"name\"]).intersection(gene_lists['Senescence'].tolist()), len(c.vs[\"name\"])) \n",
    "               for c \n",
    "               in cc_m11.subgraphs() if c.vcount() >2 and set(c.vs[\"name\"]).intersection(gene_lists['Senescence'].tolist())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyMAAAH1CAYAAADlKyHeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAC4jAAAuIwF4pT92AABtdUlEQVR4nO3debyUc//H8fen7bRoLy2KIltFol0UEW4hlCXa7LvbToWyC78bt7sQKpWlIokWSwtSokQLSlRKqbRo387398d1zTTnNDNn5mzXnNPr+XjMY665ru91XZ+Zc51zrs98N3POCQAAAADyW5GgAwAAAABwYCIZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABAIkhEAAAAAgSAZAQAAABCIYkEHgLxlZuUltYlY9YekXQGFAwAAgNRSQlLtiNfTnHOb8uvkJCOFXxtJY4MOAgAAAAXCBZI+zK+T0UwLAAAAQCBIRgAAAAAEgmZahd8fkS8++OAD1atXL6hYAAAAkEJ+/fVXdezYMXLVHzGK5gmSkcIvQ2f1evXqqUGDBkHFAgAAgNSWrwMd0UwLAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgkkPAQApxTkn55zS09ODDgUAUpaZqUiRIjKzoEPJEZIRAEDgnHPavn27Nm7cqM2bN5OIAECCSpcurXLlyqls2bIqVqzg3doXvIgBAIWKc06rVq3Spk2bgg4FAAqcbdu2adu2bfrrr79Uq1YtHXTQQUGHlBT6jAAAAvX333+TiABADjnntGLFCm3ZsiXoUJJCzUghYWY9JPWIsql0/kYCAInbvn271q5dG3QYAFAohBKSevXqFZgmWwUjSiSijqQ2QQeR2cLuwXWqqj/UBXZuBIPrreCJ9g1epUqVVKFCBRUtWjSAiACgYAj1tVu/fr22b9+eYf3mzZtVsWLFAKNLHMlI4bFU0rQo60tLapq/oQBAYrZu3Zrhdfny5VWtWrWAogGAgqV48eIqW7asli1bliEh+eeff0hGkL+cc0MkDcm83swaSJqf3/EAQFZC3+pFKl++fEDRAEDBZGaqVKmSVq5cGV63bds2OecKxLC/dGAHAAQi2vC9JUqUCCASACjYSpYsud+6gjJEOskIACAQzu3fz6YgfIsHAKmmSJH9b+mj/Y1NRSQjAAAAAAJBMgIAAAAgECQjAAAAAAJBMgIAAAAgECQjAAAAAAJBMgIAAAAgEEx6CAAoUBZ2P3CG/60/NO+H5ly6dKnq1q0bdVvp0qVVrlw5VahQQfXr11fjxo3Vrl07tWzZMqFj9+jRQ0OHDt1vfdGiRVWuXDmVK1dONWvWVOPGjdW0aVN17NhRFSpUyFHMRYsWVdmyZVW7dm01bdpUV1xxhU4//fSE4oXUt29f9evXL/y6e/fuGjJkSNx9fv31Vx155JFxy0deC1OmTFHbtm2zjGXq1Kk67bTTEj5uoho1aqS5c+cmfS7kDWpGAABAVNu2bdPq1av1888/6/3339eDDz6oVq1aqX79+ho+fHi2j7t3715t2LBBy5Yt04wZMzRgwAD17NlTNWrUUPfu3bVq1aocHXvjxo2aN2+e3njjDbVr104dO3bUtm3bsn3MVNW3b1+ZmcxMS5cuzZNzDBs2TPPmzcuTYwMSNSMAAMDXpEkTDR48OPx6z5492rhxo1atWqVvv/1WH330kRYvXqyffvpJXbt21ejRo/X222+rVKlSWR77jTfeUNOmTcOvt2zZoo0bN2rBggWaPn26PvroI+3YsUNvvvmmxo0bp2HDhuncc89NOua9e/dqzZo1+vTTT/XSSy9p+/btGjt2rHr27Kl33303yU8E6enpuv/++/Xxxx8HHUqWMl9jsSRyvSL/kIwAAABJUpkyZdSwYcOo2y6//HI999xz+uCDD3T99ddr7dq1Gjt2rK644gqNHj066gzQkerWrRv12GeffbbuuusurV69WnfddZfeeustbdiwQZ06ddLkyZOzbBIWK+YzzzxTnTt3VuvWrbVr1y6NHDlSvXr1UqNGjeIeD/tUrVpVa9eu1fjx4zV16tSEmlYFKdY1htRGMy0AAJAQM9OFF16oWbNmqXLlypKkMWPG6M0338zxsatXr64RI0borrvukiTt2LFDXbt21Z49e7J9zKZNm+qyyy4Lv540aVKO4zyQ9OrVK5xk3nfffQFHg8KKZAQAACSlTp06GjBgQPj1Y489pvT09Fw59tNPP63jjjtOkrRkyRKNGDEiR8eLbLazfPnyHB2rTp06MrNwDcG6dev00EMP6bjjjlPZsmVVtmxZnXjiiXryyScT7qPy+eefq2vXrqpbt65Kly6tsmXL6uijj9Z1112n2bNnR91nyJAhMrMMHc3r1q0b7j8SevTo0SNH77dhw4bq1q2bJGnWrFkaNWpUjo4HREMyAgAAktapU6fwiFZLlizRnDlzcuW4RYsW1R133BF+ndMb4GLFikVdzqlvv/1WjRo10qOPPqr58+dry5Yt2rJli77//nv16tVLrVu31saNG2Puv337dnXu3FlnnHGGhg8frqVLl2r79u3asmWLFi1apEGDBqlp06a6/fbbtXfv3lyLO1mPPPKISpYsKUnq3bt3jmqqgGhIRgAAQNKKFCmi9u3bh19PmzYt1459zjnnhJe/+uqrHN2MR44Edfjhh+corpBVq1apQ4cO2rZtmx588EFNmTJFs2fP1ttvv60TTjhBkvT999/rnnvuibq/c06dO3fW6NGjJXk1Li+99JJmzpyp6dOn6+mnn1aVKlXknNOLL76o2267LcP+HTt21Lx583TjjTeG102aNEnz5s3L8Hj88cdz/F5r166tW2+9VZK0ePFivfrqqzk+JhCJDuwAACBbTjrppPDyL7/8kmvHrV69umrWrKk///xTmzZt0l9//aWaNWsmfZxly5aFhyAuUaKELrroolyJb9GiRapZs6a+/vprHXHEEeH1J554os477zw1btxYixcv1rBhw/T000+rUqVKGfYfNmxYeHSqk046SZMnT1a5cuXC21u1aqUuXbqodevWWrZsmQYMGKBLLrlEbdq0kSRVqFBBFSpU0MEHHxze56ijjlKdOnVy5f1l9sADD+i1117Thg0b9Mgjj6hbt2466KCDsn2833//XVWqVEmoXF4ct1atWgnNZ4P8QTICAACypWrVquHlv//+O9eP/eeff4aPnWgykp6err/++kuffPKJevfurX/++UeS1KdPH9WqVSvX4nvxxRczJCIhZcqU0e23365bbrlFO3fu1IwZM/Ybovg///mPJK92afjw4RkSkZBatWrplVde0dlnny1J+r//+79wMpLfKlasqAceeED33nuv/vrrLz333HN6+OGHs328q666KhejS/64gwcPznF/GuQemmkBAIBsifx2PHTTn9/HnjZtWoZO20WLFlXNmjXVo0cPrVy5UrVq1dKAAQP04IMP5lps5cqVU8eOHWNub9asWXj5119/zbBtzZo14dm/27Ztq2OOOSbmcc466yzVq1dPktfRPcj+Grfeeqtq164tSXr22We1Zs2awGJB4UIyAgAAsmXz5s3h5fLly6fksTt06KCePXvmRkhhRx99tIoWLRpze2jYY2n/ROrHH38ML7dq1SrLc4XKbN26VUuWLEk21FxTsmRJPfLII5K8CStDy9kxZcoUOeeyfEyZMiVPjkutSGohGQEAANmydu3a8HLmfhH5dewmTZpk6LQ9bdo0DRw4MDz53csvv6wOHTrkaq1C6dKl426PnAAyc+f7yOZs1atXz/JcNWrUiLpvELp16xYedvnVV1/dr9YHyA6SEQAAkC3fffddePnYY4/NteOuWrVKq1atkuR11o530x6agT30OPXUU3XDDTdo9uzZ4VG5Pv/8c/Xt2zfX4ssJ51x42cwCjCR5RYoU0ZNPPilJ2r17t3r37h1wRCgMSEYAAEDS9u7dq08++ST8OjQRYG4YP358ePnUU0/NUNOQqBIlSujNN98MN5l65plnAm3mFBI52lMo4Yonskxk86+gnHvuueGO9KNGjcqQkALZQTICAACSNnr0aC1btkySN6xso0aNcuW4e/fu1fPPPx9+3blz52wfq0qVKurTp48kadeuXerVq1dOw8ux448/Prw8Y8aMLMuHypQpU2a/0buCqlnp37+/JK+W59577w0kBhQeJCMAACApS5cu1c033xx+3adPn2zVXkRz7733av78+ZK8JOfyyy/P0fFuvPFGHXLIIZK8b/IjJ0EMwsEHHxyeGHHKlClatGhRzLKffvqpFi9eLElq167dfjPIh2ZGl6SdO3fmfrAxNGvWTJ06dZLkvYeJEyfm27lR+JCMAACAhDjnNGbMGDVr1izcmbpTp0668sorc3zs1atX64orrtD//d//SZJKlSqlYcOGxR21KhFpaWnhvg3OuZToO3LHHXdI8uZE6dq1q7Zs2bJfmT///FPXX399+PWdd965X5lQkiUpblKTF5544olwcpSTkbUAJj0EAACSvOFjQ7USkrRnzx5t2rRJq1at0rfffqtx48aFv6mXpIsuukjDhg1LqLlQ5tmxt27dqo0bN2rBggWaPn26PvroI+3atUuSN3rW8OHDM8zXkRNXX321nn76aS1btkxjxozR3Llzw7UTQejatatGjhypjz/+WLNmzVKjRo10991366STTtLevXs1ffp0PfPMM+G5PG666aaoEx6ecsopMjM559SrVy+ZmerVqxdOEsqVK5dhlvbcdOSRR+raa6/VwIEDM4x8FqREZ2CXpPr16+dabR5yhmQEAABI8kbHCg3dGk/9+vXVu3dvdenSJeFjJzI7dqlSpdS5c2c9/fTTCQ17m6gSJUqoT58+uvbaa+Wc08MPP6yxY8fm2vGTZWYaNWqUunXrptGjR+u3337TTTfdFLXcLbfcEp6xPbPDDjtM11xzjQYNGqT58+frvPPOy7C9e/fuGjJkSF68BUnSww8/rDfffFNbt27Ns3MkI5mZ3Tds2KAKFSpE3VbQRjkr6EhGAABAVKVKlVLZsmVVsWJF1a9fX40bN1b79u3VvHnzHB23aNGiOuigg1SuXDnVqlVLjRs3VtOmTdWxY8eYN4g51aNHDz311FNasmSJPvzwQ82ePVsnnXRSnpwrEaVKldKoUaP02WefaciQIZo+fbpWr14dnkG+TZs2uuGGG7KM8eWXX1azZs301ltvaf78+dqwYUO+zdRerVo13XXXXYWimVZkn5us5pFB7rLI8a5R+JhZA0nhOvf58+erQYMG+Xb+hd2D+3ah/lCu7QMN11vBsmfPngxNfiSv6UfmTroAkNfeeeed8GAJffr00aOPPhpwRMnJyd/TBQsWhCcJ9TV0zi3I3Qhjo7EcAAAADmg//vhjeLl+/foBRnLgIRkBAADAAWvbtm166623JHlNCE8//fSAIzqwUBcOAACAA8rKlSu1Zs0a/frrr3rmmWfCE3hefvnlqlatWsDRHVhIRgAAAHBA6d27t4YOHZphXePGjfX8888HE9ABjGQEAAAAB5wiRYqoQoUKql+/vjp16qTrr78+w6z2yB8kIwAAADigDBkyJE/nYEHi6MAOAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBDFgg4AAIBk2OC7gw4h37iez+bLeZYuXaq6deuGX1eoUEG//fabKlasGHe/1q1ba/r06ZIk51yGbVOnTtVpp50mSerevbuGDBmSUCx16tTRsmXLoh4z83GT8f333+uEE05I6lwA8h41I7nMzOqa2RYzc/5jadAxAQCQjI0bN+qJJ54IOgwABwBqRnKRmZmkNySVCToWAABy4qWXXtKtt96qQw89NOhQYrrgggv02GOPJVS2Xr16eRwNgOwgGcldN0tqK+kvSdWCDQUAgORVrVpVa9eu1Y4dO/TQQw8l3LwqCBUqVFDDhg2DDgNADtBMK5eY2eGSnvJf3hxkLAAAZFf79u3VvHlzSdKwYcM0b968gCMCUJiRjOSCTM2z3nfOvRdwSAAAZFv//v0lSenp6br//vsDjgZAYUYykjtukdRG0gZRKwIAKOBOPfVUnXvuuZKk8ePHa9q0aQFHBKCwIhnJITM7QtKT/ss7nXOrg4wHAIDc8NRTT6lIEe824d577w04GgCFFR3YcyBT86xPnHND8vh8B0uqmuRuR+RFLACAwq1hw4bq1q2bhgwZolmzZmn06NHq1KlTto61ceNGzZ8/P6Gyu3fvzvXjlilTJsM8KgBSB8lIztwq6VRJWyRdlw/nu0nSw/lwHgAA9Mgjj+idd97Rjh071KtXL3Xs2FHFiiV/6zB27FiNHTs21+NL9Lht2rTR1KlTc/38AHKOZlrZlKl51gPOuWVBxgMAQG6rXbu2br31VknS4sWLNWjQoIAjAlDYkIxkg988a7Ck0pK+kvS/YCMCACBvPPDAA6pYsaIkqV+/ftq6dWvSx+jevbuccwk9DjvssFw/LrUiQOoiGcme2ySdImmHpGuccy6fzjtAUsMkHxfkU2wAgEKoYsWKeuCBByRJf/31l5577rmAIwJQmJCMJMnM6kl6wn/Zzzn3S36d2zm3xjm3IJmHpCX5FR8AoHC69dZbVbt2bUnSs88+qzVr1gQcEYDCgmQkeffIa561SdIKM7ss2iOifJmI9ecFEzIAANlXsmRJPfLII5KkzZs369FHHw04IgCFBaNpJS/Nfy4vaVgC5atIettfXiZpXF4EBQBAXurWrZuee+45zZ8/X6+88or+/e9/Bx0SgEKAmhEAAJClIkWK6KmnnpLkzQXSu3fvgCMCUBiQjCTJOdfDOWdZPSJ2WRaxvk5QcQMAkFPnnnuu2rRpI0kaOXKkfvkl37pNAiikaKYFAAAS9vTTT6tFixZyzmndunWBxpLMzO41atRQ5cqV8zgiAMkiGQEAAAlr3ry5Lr74Yr333ntBh5LUzO7/+c9/6OcCpCCSEQBAgeJ6Pht0CAe8J554QmPHjtWePXuCDiVHdu7cKUkqXbp0wJEABy6SEQAADnB16tRRMvP3HnXUUdq9e3fcMm3btk3qmCFLly7Nk+Nm5pzTxo0bJUmVKlXK8fEAZA/JSB7J1IkdAACkkMWLF2vHjh2SpPr16wccDXDgYjQtAABwwHnttdfCy+3btw8wEuDARs0IAAAo9LZu3arffvtN69at03vvvacBAwZIkipWrKhu3boFHB1w4CIZAQAAhd63336r0047LcO6kiVLavjw4apatWpAUQEgGQEAAAeM0qVLq3bt2mrXrp3uuOMO1atXL+iQgAMayQgAACj0cmsULgC5iw7sAAAAAAKRcsmImZUNOgYAAAAAeS/lkhFJf5rZK2Z2UtCBAAAAAMg7qZiMlJF0jaRZZvadmV1jZmWCDgoAAABA7krFZORrSeY/Gkt6RV5tyQAzOyHIwAAAAADknpRLRpxzrSXVl/SipA3ykpKykq6XNNvMvjGznmZWKsAwAQAAAORQyiUjkuSc+9k5929Jh0jqJulL7astaSLpNXm1Jf81s+MCCxQAAABAtqVkMhLinNvpnBvunGsj6VhJz0taLy8pKS/pJklzzWy6mXU1s5LBRQsAAAAgGSmdjERyzv3inLtTXm3JlZK+0L7akhaShkhaaWb/MbNjAwsUAAAAQEIKTDIS4pzb5Zx7yznXVlJrSav8TSapoqTbJM03s0/M7OSAwgQAAACQhQKXjEiSmbUxsxGSPpdUXV4iIklbtK+2pJ2kL8zsJTOz6EcCAAAAEJQCk4yYWWUzu8vMfpY0WdJlktLkJR5fS+oqqYqkBpJekrTd33ajpBsCCRoAAABATCmfjJjZaWb2tqQVkvpLOkpekrFV0suSGjnnWjvnRvhNuH5yzt0m6UhJs/2y1wUUPgAAAIAYUjIZMbMqZnaPmf0i6TNJl2hfLcg8STdLqumcu8k5Ny/aMZxzqyT19l/Wy4ewAQAAUABMnTpVZiYz05AhQ4IO54BWLOgAMjOzdyR1lFRc+/qC7JI0WtJA59z0JA73m/9cOtcCBAAEyu4aF3QI+cY9d16+nGfp0qWqW7du1G1FixZV2bJlVbt2bTVt2lRXXHGFTj/99HyJC0Dhl3LJiLxakJDfJb0i6Q3n3LpsHGuXpOWS0nMjMAAADjR79+7Vxo0btXHjRs2bN09vvPGGLrjgAr311lsqXZrv+gDkTComI07SR5IGSprknHPZPpBzyyXVyaW4AAAo9Jo0aaLBgweHX+/du1dr1qzRp59+qpdeeknbt2/X2LFj1bNnT7377rsBRgpkX9u2bZWDW0zkolRMRuo45/4IOggAAA5EZcqUUcOGDfdbf+aZZ6pz585q3bq1du3apZEjR6pXr15q1KhRAFECKCxSrgM7iQgAAKmpadOmuuyyy8KvJ02aFGA0AAqDlEtGJMnMDvUfaQmULRkqnx+xAQBwIGvatGl4efny5THLrVixQr1791azZs1UpUoVlShRQtWqVdOZZ56pAQMGaOfOnTH37du3b3iko6VLl2rPnj0aOHCgTj75ZFWtWlWlSpXS0UcfrTvuuEOrVq2KG2/oOD169JAk/fTTT7rlllt09NFHq2zZsjIzffDBBxn22bNnjwYPHqxzzz1XNWvWVFpamipVqqSTTjpJvXr10p9//pn1ByXpjz/+0G233aYjjzxSpUqVUtWqVdW6dWu9+uqr2rt3r5YuXRqOr2/fvvvtH23Ep08++USXXHKJDjvsMKWlpcnMtHHjxvA+f//9twYPHqxu3brp+OOPV/ny5VW8eHFVrlxZzZo10/33368//oj/ve+QIUPC5506daokafTo0Tr77LNVo0aN8Od/9913a82aNRn2XbJkiW6//XYdc8wxKl26tCpVqqSzzz5bn3/+edxz1qlTR2amtm3bSpJWrlype+65R8cee6wOOuggVa1aVaeffrrGjcs4gEV6erpGjBih008/XTVq1FDJkiV11FFH6YEHHtA///wT83xZjaYVbfuMGTPUpUsXHXrooUpLS1PVqlV17rnnauLEiXHfW0hOr4fCKuWaaZlZK0lfyZtHpI6k2H+tPGUk/SQpzcyaO+dm522EAAAcuIoVKxZ1OdILL7yg++67b7+EY82aNfrss8/02Wef6fnnn9eHH36oY445Ju75/vnnH7Vp00Zff/11hvWLFi3SokWLNHjwYI0ZM0annXZalrEPHTpU119/fdxE6Pfff9f555+v+fPnZ1i/a9cubdiwQXPmzNHzzz+vV155RV27do15nAkTJqhz587aunVreN2OHTu0bt06TZ8+XSNGjNBLL72UZcyRbrrpJg0cODBumSOOOEKbNm3ab/369eu1fv16ffvtt3rxxRf1xhtvZKjlimXv3r264oor9NZbb2VYv2jRIj333HMaPXq0pk6dqjp16mjs2LG68sortWXLlnC57du3a9KkSZo0aZIGDhyoG27Ieh7qr776Sh07dtTff/8dXrd161ZNmTJFU6ZM0cMPP6y+ffvqn3/+0WWXXaYJEyZk2H/x4sV66qmnNG7cOH3xxReqVKlSlufMyqOPPqq+ffsqPX3fmEjr1q3T+PHjNX78ePXr108PPfRQzP3z4nooLFKxZuQi//lj59zfcUtK8st8KO+9dMrLwAAAONDNm7dveq/DDz98v+2PPfaY/v3vf2vnzp2qXbt2+KZw9uzZmjhxom677TaVKFFCixcv1hlnnKG//vor7vmuu+46ff311zr77LP1/vvva/bs2fr444/VpUsXSdKmTZvUoUMH/fzzz3GP89133+naa69VhQoV9OSTT+rLL7/UrFmzNGTIENWpU0eSd8Pepk2bcCJyyimn6K233tJ3332nzz77TLfddpuKFy+u7du3q1u3bho1alTUc/3www+68MILtXXrVpmZunXrpvHjx2v27NkaM2aMzj//fH3xxRe6/vrr48Yc6YUXXtDAgQN1wgkn6LXXXtM333yjr776Ss8++6xKlCgRLrd37161aNFC/fr10wcffKCZM2dq1qxZGj16tG666SaVLFlS27dvV9euXTVz5swsz/vQQw/prbfe0rnnnqvRo0eHf44XXnihJGnZsmW69tprNXv2bHXu3FmVKlXS888/r6+//lrffPONnnjiifCoa7fffrt+/fXXuOdbtWqVLrzwQpmZHn/8cX355ZfhBKpy5cqSpH79+umrr75Sz549NXHiRHXp0kUffvihZs+erQ8//FBt2rSRJC1YsED33HNPwp9xLG+88YYeeughNW7cWK+//rq++eYbff3113rsscd00EEHSfJq87766quo++fF9VCYpFzNiKRT5I2o9UkS+0yUdKmkNnkSEQAA0LJlyzR8+HBJUokSJXTRRRdl2D5jxgw9/PDDkqSrr75aAwcOVPHixTOUOeuss3TZZZepXbt2Wrlypfr06aNBgwbFPOc333yju+++W88880yG9f/617/UqlUr3XLLLdq2bZtuuOGGcJOiaBYsWKAjjjhCX331lapXrx5eH9ns7O677w43Ybruuuv08ssvy8zC29u1a6cLLrhA55xzjnbt2qXrrrtO7du3V/ny5TOc64YbbgjXvrz++uvq2bNneNuJJ56ojh076uGHH9YjjzwSM97M5s6dq44dO2rUqFEZaqROPvnkDOVmz56to446ar/9mzZtqosvvlh33nmnTj75ZP3111/q06ePPvvss7jn/frrr8M1EZHat2+v888/Xx999JE+++wz/etf/1KDBg00efJkVaxYMVyuWbNmqlevni655BLt2rVLr7zyyn4/y0iLFi1S7dq1NX36dNWuXTu8vkmTJmrevLlatGgh55wuvvhirVmzRoMHDw43wZO8z/fss89Ws2bNNHfuXA0fPlzPPPNMjmpHvvzyS3Xp0kVvvvmmihYtGl7fsmVLNW3aVGeddZacc3rhhRfUunXr/fbPi+uhMEnFmpFD/OefkthnUaZ9AQBALkhPT9eqVas0dOhQnXzyyeF2+H369FGtWrUylH388ceVnp6uww8/PGoiEtKyZUvddNNNkqRhw4Zp+/btMc9/9NFH68knn4y67eabbw5PwDht2jTNnTs37nsZMGBAhkQk0tq1azVixAhJUt26dfXiiy9mSERCTj/9dN15552SpI0bN+qNN97IsH3OnDnhGodzzjknw41npIcffjjqqGWxHHTQQXr99ddjNo0LiZaIRDriiCN03333SZImT56sDRs2xC3fuHHjcIIZycx0++23h1+HEoPIRCSkU6dO4WtlypQpcc8nSS+++GKGRCSkWbNmatWqVfh8HTt2zJCIhBQvXjx8fe3atWu/Jn7JOvjggzVo0KAMiUhI+/bt1aRJE0mKmgzn1fVQmKRiMlLVf96WxD6hsgfnciwAABxQpk2bFu5Ea2YqWrSoatasqR49emjlypWqVauWBgwYoAcffDDDflu2bAmPrtWpU6eYiUhIqKPyzp079d1338Usd9VVV8W9Ab/uuuvCy/E6EtesWVNnnnlmzO1TpkzRrl27JEk9e/ZUWlrsMXRuvvnm8HLmEcU++WRfw45YN56SVKRIkag30rGcd9552fp2f+3atfr111+1YMECzZ8/X/Pnzw83LXLOZZnAXXHFFVGTMslLVEIaNmyoE044IWo5MwtvW7JkSdzzlStXTueff37M7ZHnvPLKK2OWO/HEE8PLWZ0zK507d447wWezZs0keX1IMvfXyavroTBJxWZamyRVllQtiX1CZZNJYAAAQJI6dOgQ9aZqzpw52rNnjySpf//+6t+/f8LHjDciVosWLeLuG7n9hx9+iFnuhBNOiHlTLUk//vhjeDn07XsstWrV0qGHHqrly5fvd87IPjWRTcCiyWp7pMib8KyMGzdOgwYN0pdffplhpK1o1q1bF3d7vAEGImtBjj322LjHCZWNN8KV5NXsFCkS+7vyRM8ZWS6rc2alfv36cbeH+rKEzhXZbC+vrofCJBWTkd/lJSOnKfF+I6f7z7HHGAQAAFnKPAP7+vXrtXDhQv3vf//T/Pnz9fLLL2vx4sWaOHFihhqLzEO8JmPbttjfJVarFv+7ychmV/FurLOqVYgcuSlWU65INWrU0PLlyzPsl/k4Bx8cv8FGVu8tUiK1Irt379aVV16pkSNHJnzceJ+95E2CGUtk0hCv5iCybORoVMmeL5lzRpbbu3dv3GNmJdH3Fu1ceXU9FCapmIxMkdRU0vVm9oJzbnW8wmZWQ9L18jq9T86H+AAAKLSizcB+6qmn6qqrrlLHjh01YcIEff755+rbt68ee+yxcJlQrYgkPfjgg7rkkksSPmfmvieR4tVmJFM2Wnv/SM65bJ0zJyLPmZWs4pekp59+OpyINGjQQLfffrtatGih2rVrq0yZMuGmc5MnT1a7du2SjgF560D9WaRiMjJI0l2SykuaaGYXOud+j1bQzA6X9L5fdq+kV/MtSgAADiAlSpTQm2++qWOOOUZ///23nnnmGfXs2VNHHHGEJKlq1arhstu2bcu1zrirV6+O2yl79ep931lGNpdJVpUqVcLLq1atyrJpTqhpWeZzRtZgrFmzJjxscDQ5qU2KJjQHyeGHH65vvvkmZi1D5toc5J0gr4eCIuU6sDvnlkj6P0km6ThJC81shJldY2b/MrNz/OURkhb4ZZyk551zvwQXOQAAhVuVKlXUp08fSd4oRb169Qpva9y4cbi5yvjx43PtnFnNhRG5vVGjRtk+z/HHHx9enjFjRtyyK1euDM8+n/mckUnYrFmz4h4nq+3J+Pvvv8Mzw19wwQVxmzvl5nkRX1DXQ0GScsmI7wFJI+QlJGmSLpP0iqRxkj7yly/zt5mk4ZLuDSRSAAAOIDfeeKMOOcQbSX/UqFHhDrqVKlUKj5D1008/afTo0blyvsGDB2doApbZa6+9Fl5u3759ts9z2mmnhScPHDx4cHhkrWgGDBgQXj7rrLMybIuMYciQITGPkZ6erqFDh2Yz2v1FfkaRs3xn9s8//+TqeRFfUNdDQZKSyYhzLt0511XSVZKWyEs4oj1+ldTdOdfdHagN7QAAyEdpaWnq3bu3JK+Ne+RkeH379g33t7j22ms1ffr0uMdasWKFXn/99bhlfv7553BtTGYvv/xyeNK+U045JcNwrsmqWrWqrrjiCknSb7/9pjvuuCNqualTp+rZZ5+VJFWoUEFXXXVVhu1NmjQJD/U6YcKEDIMBROrXr194pvfcULVq1XCToA8//DDqzPY7duzQFVdcobVr1+baeRFfUNdDQZKSyUiIc26Ic+5ISY3kJSb3SbrfXz7eOXeUc25YkDECAHCgufrqq3XYYYdJksaMGROeq+KUU07RE088IcmbELBNmzbq0qWL3n33XX377beaPXu2Jk2apGeffVZnnnmm6tSpE/PmLKR58+Z6+umnde655+qDDz7Q999/rwkTJqhr16668cYbJUmlSpXSyy+/nOP39eyzz4Yn2xswYIBOO+00jRw5UnPmzNHkyZN1xx13qH379uFak1dffXW/2dclL0kK1bJcffXV6tGjhyZOnKjvv/9eY8eOVceOHfXII49kGJY4p53mixQpou7du0vy+tG0aNFCAwcO1Ndff62ZM2dqwIABOuGEE/TRRx/plFNOydG5kJwgroeCJBU7sO/HOTdP0rwsCwIAgDxXokQJ9enTR9dee62cc3r44Yc1duxYSdL999+vatWq6fbbb9fmzZv19ttv6+233455rGg385FeeeUV3XjjjRo/fnzUvijlypXT+++/n2WH80RUqlRJ06ZN0/nnn6/58+dr6tSpUWfVDiU/nTt3jnqcxo0ba/To0br88su1detWDR06dL8mOC1bttT//ve/cG1OyZIlcxz/o48+qlmzZmn69OlaunRpeBbySD179tSVV14ZHk0LeS+o66GgSOmaEQAAkJp69OgRHknrww8/1OzZs8PbevbsqeXLl4drQGrUqKG0tDSlpaWpRo0aOuWUU3TPPfdo8uTJGjduXNzzlC9fXtOmTdOLL76oli1bqlKlSkpLS1O9evV0++2366effsrVG+u6devq+++/1+uvv65zzjlH1atXV/HixVWhQgU1btxY999/v3799Vd169Yt7nHOO+88LVy4UDfffLOOOOIIpaWlqVKlSmrZsqVefPFFTZs2LcNQrlklZYkoU6aMpkyZohdeeEHNmzdX2bJllZaWptq1a6tTp076+OOP9cYbb8SdVBB5I4jroaAwuloUbmbWQFK4EeL8+fPVoEGDfDv/wu7BVTPWH8q1faDheitY9uzZo8WLF2dYd+SRR2aYSA8Hpr59+6pfv36SpN9//z3ucKgF2dChQ9WjRw9J3twfp512WrABIVA5uR5y8vd0wYIFmYfibuicW5DwyXMopf/im9dg7hhJ9SSVk5TljD/OuTfzOi4AAICcevNN75YlLS1NTZo0CTgaBO1AvR5SMhkxs6Lyhuq9RVL1JHZ1kkhGAABAoBYtWhR3ssbnnntOkydPliR17txZZcuWza/QEACuh9hSLhnxE5GPJIUGZj5whhMAAACFQrdu3bR7925dfPHFatKkiapWraqdO3fqp59+0ogRI/T5559L8oYHDo1AhsKL6yG2lEtGJF0nKTSD0HpJb0iaJelvSelBBQUAAJCMOXPmaM6cOTG3V69eXWPHjg0PJ4zCjeshulRMRq70n3+TdLJzbv9ZewAAAFLYSy+9pHHjxmnKlCn6448/tG7dOu3atUsVK1ZUw4YN1aFDB1177bUqU6ZM0KEiH3A9xJZyo2mZ2QZ5ndWvds4NCTicAo/RtHAg4XorWBhNCwByR0EeTSsVB5ou4T//EGgUAAAAAPJUKiYjf/jPpQKNAgAAAECeSsVkZLz/3DLQKAAAAADkqVRMRp6XtFnSHWZWKeBYAAAAAOSRlOsl6JxbbmaXSxolabKZ9XDOzQ04rJRnZj0k9YiyqXT+RgIAAAAkJuWSETN7w1/8WdKJkmab2QJJCyVty2J355y7Oi/jS2F1JLUJOggAAAAgUSmXjMj7dj80RqaTNwN7A/+RiAM1GVkqaVqU9aUlNc3fUAAAAICspWIyInkJSLzXyMSfk2VI5vWZ5xkBAAAAUkXKJSPOuVTsVA8AAAAgl3HjDwAAACAQJCMAAAAAAkEyAgAAACAQKddnJJKZHSqpp6TWkmpIKiXpLOfcrxFljpFUU9JW59w3gQQKAAAAIGkpm4yYWS9JD0kqHlolb6jfEpmK1pc0WtJuM6vlnFubf1ECAAAAyK6UbKZlZo9IelRe4rFb0pw4xcdI+lNeYnVh3kcHAAAAIDekXM2ImR0nqbf/cpSkG51z680sPVp555wzs/cl3SLpdEmv5k+kAIAgbFrZLugQ8k35Qz7P1/Pt2bNHH3zwgUaPHq05c+Zo9erV2rFjhypXrqyqVauqVq1aatmypVq3bq3mzZurdOnS+RofgMIn5ZIRSTfLa5L1vaTLnXNRk5BMZspLRhrmZWAAABRWCxcu1OWXX64ff/xxv22rV6/W6tWrNW/ePE2YMEGSdP311+vll1/O7zABFDKpmIy0kdc35H8JJiKStNR/rpUnEQEAUIj98ccfatOmjdatWydJOumkk9StWzc1atRI5cuX19atW7Vs2TLNmjVL48eP1+LFiwOOGEBhkYrJSCih2P+rmdi2+s/UFwMAkKQHHnggnIjcddddeuaZZ2RmGcqcfPLJ6tKli55//nnNmjVLK1euDCJUAIVMKiYjoU71xeOWyqii/7wll2MBAKBQ27t3rz744ANJUtWqVfXUU0/tl4hk1qxZs3yIDMCBIBVH01rjP9dNYp/m/jNf0wAAkIS1a9dq61avgcHhhx+uYsVy/j2lc06jR4/WpZdeqjp16qh06dI66KCDdMwxx+j666+P2i8lkpnJzNSjRw9J0rJly3TnnXfq6KOPVunSpVWhQgW1atVKAwYM0J49exKK6eeff9Ydd9yhRo0aqVKlSkpLS1PNmjV13nnnafjw4dq7d2/c/VeuXKk+ffqoWbNmqlixoooXL65KlSrp6KOP1plnnqnHH39cc+fOjXuMVatWqV+/fjrllFNUrVo1lShRQtWrV9dJJ52kG2+8UZ9++mnc97NixQr17t1bzZo1U5UqVVSiRAlVq1ZNZ555pgYMGKCdO3fG3Ldv377hz3Xp0qWSpLfffltnnnmmqlWrprS0NB166KHq0aOHfvrpp7jvI2TPnj0aPny4OnfurLp166pMmTIqV66cjj32WJ1//vl67bXXwjVu0ezevVuDBw/W+eefr1q1aqlkyZIqX768jj/+eN1xxx36/fffE4oDBZs554KOIQMze0dSZ0nvOue6RKxPl9eX5Djn3MKI9SUkLZSXvAxyzt2QzyGnNDNrIGl+6PX8+fPVoEGDfDv/wu7xv13LS/WHpta1jbzH9Vaw7NmzZ7++B0ceeWSWN8OMppW7NmzYoEqVKkmSKlWqpNWrV6t48WQaJ2T0xx9/6OKLL9a3334bs4yZqXfv3nrkkUei1sKE1nXv3l2XXHKJunTpok2bNkU9Vvv27fXRRx/FjDk9PV29evXSs88+GzfhaNasmcaMGaOaNWvut+3jjz/WpZdeGk7aYmnTpo2mTp0addvzzz+vBx54QDt27Ih7jClTpqht27b7rX/hhRd03333xU04jjzySH344Yc65phj9tvWt29f9evXT5L0008/6d5779W4ceOiHqdkyZL64IMPdNZZZ8U81w8//KBLLrlEixYtivt+unfvriFDhuy3ft68ebr44ovj9j8qXry4XnzxRd1wA7d2Wcnu31NJWrBggRo2zDAGVEPn3ILcjTC2VGymNUrSJZI6m9kzzrnvYxU0s6KSXpN0uLxEZUT+hAgAQOFQsWJFHX744frtt9+0fv163XDDDXrppZdUqlSppI+1evVqtWzZUitXrlSRIkXUqVMnnXfeeapbt66KFy+uH3/8Uf/73/80d+5cPfbYY0pLS1OfPn1iHm/evHkaNWqUypcvr169eqlVq1YqVaqUfvjhBz322GP6/fff9cknn6h///7q3bt31GNcd911ev311yVJxx57rG644QYdddRROvjgg7Vy5Uq9//77Gjp0qGbNmqVzzjlHM2bMyDBk8bp169SlSxdt3bpVpUqV0jXXXKP27durWrVqkqS//vpLc+fO1cSJE2M2b+vTp48ef/xxSd6Nfs+ePXX22WerVq1a2rVrlxYtWqQpU6ZozJgxUfd/7LHH9OCDD0qSateurZtvvlkNGjRQzZo1tXbtWo0fP14vv/yyFi9erDPOOEOzZ88OxxfrM/nyyy918cUXq0uXLqpTp442bNigd999V4MGDdKOHTvUtWtXLVq0SBUqVNhv/x9++EEnn3xyODlr06aNunbtqvr16ystLU1//vmnvvnmG73//vtRz79w4UKdfPLJ2rx5s0qUKKFu3bqpXbt2qlOnjpxz+vbbb/Xiiy9qyZIluvHGG3XQQQfpyiuvjPl+ULClXM2IJJnZLElNJP0t6R5JI+X1B3Hyhu/9Q9LZku6X1Njf7TPnXOwU/gBFzQgOJFxvBQs1I1nLr3lGXn75Zd14443h1xUrVtS5556rVq1aqUmTJjr++OOVlpaW5XH+9a9/acKECapUqZImTZqkJk2a7Fdmz5496tKli0aNGqXixYtr0aJFqlOnToYykTf1DRs21OTJk1W1atUMZVavXq0GDRpo/fr1ql69ulasWKGiRYtmKDNy5EhdeumlkqSHH35YDz30kIoU2b+F+pgxY3TxxRfLOafHHnssQ2IzePBgXXXVVZKk999/XxdeGHt+5b///luVK1fOsG7q1Kk67bTTJEm1atXSpEmTVL9+/aj7b968WXv27FHFihXD62bMmKHWrVsrPT1dV199tQYOHBi1FmjGjBlq166dtm/frmuuuUaDBg3KsD2yZkSSBg4cGLXGoXfv3nriiSckSS+99JJuvvnmDNv37NmjY489Vr/++qskr8bmtttui/mZLF++XIceemj4dXp6uo4//ngtWLBAdevW1aeffqojjjhiv/22bt2qs846S9OnT1elSpW0dOlSlS1bNuZ5DnQFuWYkFfuMSNLF8mZVryzpdUmb5CUikjRN0kZJ78pLREzSMkld8z1KAAAKgRtuuEH33HNP+PWGDRs0fPhw3XTTTWrWrJnKli2rFi1aqG/fvjHb8X/33XfhOUj69+8fNRGRpGLFimngwIEqUaJEuM9APEOGDNkvEZGk6tWrh5OE1atX6+eff96vzCOPPCJJOuWUU9S3b9+oiYgkXXjhhbroooskSa++mnHu5FWrVoWX27WLnwhnTkQiY5CkESNGxExEJKls2bIZEhFJevzxx5Wenq7DDz88ZiIiSS1bttRNN90kSRo2bJi2b98e8zwdOnSI2fTp7rvvDp8jWpOzt956K5yIXHXVVXETEUkZEhHJS+gWLPDucwcNGhQ1EZGkMmXKhOexWb9+vUaPHh33PCi4UjIZcc79IekkSRPlJRtF/WdJqiIvbvMfn0lq4ZxbE+VQAAAgAf3799esWbN06aWX7jez+u7du/XNN9+oX79+Ouqoo3TnnXdq165dGcq89957kqQiRYqEayNiqVy5so477jhJ0ldffRWzXIMGDXTSSSfF3B45qlfoBjnkl19+Cd/0Xn755XHjkRTup7F8+XL98ccf4fWHHHJIeDlzopKVjRs3atq0aZKkFi1a6NRTT01q/y1btmjSpEmSpE6dOmXZlyf0Hnbu3KnvvvsuZrnu3bvH3FaxYkUdeeSRkvb/TCWFR16TpHvvvTduPNGErpODDz44y+SuYcOG4QQv3nWCgi0V+4xIkpxzf0n6l5mdJK+mpLmkg+XFvFbSbEnvOee4OgEAyAVNmzbVO++8E76Z/fbbbzVv3jzNnDlTCxd6Y8fs2bNH//nPf/Tbb79pzJgx4SZV33zzjSSvGU4yzWkiax4yi1eLIGWsifjnn38ybAvFI0k33XRTuNYg0Zhq164tSerYsaOqVq2qtWvX6p577tHw4cN1wQUXqHXr1mrevLnKlSsX8zjff/+90tO9+ZtDTbWSMWfOnPDoWv3791f//v2Teg+xJPq5Zv5MJWn27NmSvJqpo48+OuF4QkI/lzVr1mQ5hHSkeO8HBVvKJiMhzrnZ8hIPAACQD9LS0nTyySfr5JNPDq9btGiRHnroIb377ruSpLFjx2r06NHq3LmzJO/mMju2bdsWc1vmGprMIptdZR4pK7vxZI6pfPnymjBhgrp06aJFixbphx9+0A8//BA+/wknnKCOHTvq2muvVfXq1TMcZ+3ateHlaKN0ZSW33kNmiX6u0UYfC72nyBqjZOTFdYKCLeWTEQAAELyjjjoqXGsSaqrz9ttvh5OR0Df4lStXjjm8bTQlSpTI7VAzxCN5nbVbt26d8L5162ac6uykk07SwoULNWHCBH300UeaPn26Fi5cqPT0dM2ZM0dz5szRU089pVdeeSXDqE+RgwQlUwsQ7T08+OCDuuSSSxLet1atWkmfLxnZeT/Svvd03HHH6a233kp4vzJlymTrfEh9JCMAACBhN998czgZ+eWXX8Lrq1atql9++UXr16/XoYceGrf5Un6I7PS+Z8+ezKMFJa1o0aLq0KGDOnToIMnrD/LFF1/orbfe0siRI7Vt2zb16NFDJ5xwQvhcBx98cHj/lSuTn5c58j1s27Ytx+8hN1StWlXLly/XihUrsr3/ihUrtGrVKjVo0CDbSQ0Kj5TswA4AAFJTZPOcyGZSoY7mzjlNnDgx3+PKLLLj+/jx43P9+BUqVND555+vd955JzyHyN69e8PN2CTpxBNPDA83PGXKlKTP0bhx4/BnnBfvITtCo6StXr06QzKaqNDPZd26dXEnxsSBI+WSETP7LQePJUHHDwBAYTZr1qzwcuSwrBdffHF4+amnntLu3bvzNa7MGjVqFI5v4sSJeXrjGzlTeWQ/kfLly4dHuJo5c6a++OKLpI5bqVKl8P4//fRTSgxvGznPSjId6kMir5NHH300V2JCwZZyyYikOgk+DvMfmdcDAIAE7dixQyeeeKJGjhy533C9mS1evFi9evUKv+7UqVN4+ZRTTtEZZ5whyRtFqmfPntq5c2fMY6Wnp+u9994Lj9KV28wsPMmfc06dOnXSTz/9FHefn3/+OUPNhiRNmDBBy5cvj7tfZK1F5nkzHnzwwXBTpCuuuCJuDFu2bNGGDRsyrOvbt294/2uvvVbTp0+PG8uKFSvCM87nhcsuu0xHHXWUJOmNN97Qf//737jlI4dJlrxhlo855hhJ0kcffaT77rtP8Sbg3rVrl9544w399ddfOYwcqSoV+4x8oX0THMZSRtJRksr5ZRdJWp3HcQEAUCh9//33uvTSS1WxYkWdd955atGihY4++mhVrFhRe/bs0dKlS/XZZ59lmEyvbdu26tKlS4bjDB8+XM2bN9eyZcs0YsQIzZgxQ9dee62aN2+uihUrauvWrVq6dKlmzpypMWPGaOXKlfr000+zHGo2u6644gp9+eWXeuWVV7R8+XKdeOKJ6tatm84++2wdeuihSk9P15o1azR37lyNHz9eM2bMUJcuXTLMk/Luu+9q+PDhOuWUU3TWWWfp+OOP18EHH6z09HStWLFCY8aM0YgRIyR5TbciO7BLUps2bdSrVy89/vjjWrFihU488URdddVVOvvss3XIIYdo9+7dWrJkiSZPnqz33ntPY8aMCdeGSF6S98QTT+iBBx7Qxo0b1aZNG11yySW64IILdPjhh6tIkSJat26d5s2bp0mTJmnKlClq0aKFrr766jz5TIsVK6aRI0fq5JNP1tatW3Xbbbfp/fffV9euXVW/fn2lpaVp1apV+vbbbzVq1Cg1adJEQ4YMybD/+++/r1atWmnjxo3q37+/Jk6cqKuuukonnniiypYtq82bN2vx4sX6+uuv9cEHH+jvv//W4sWLVa1atTx5TwhWyiUjzrm2iZQz72uCjpJelFRB0pX+MMAAACBBRYsWVc2aNfXnn39qw4YNevPNN/Xmm2/G3adLly569dVX95vRvFq1apo5c6a6du2qzz77TL/99pseeOCBuOfO61GSBg4cqMMOO0z9+vXTjh079Oqrr8advLB8+fL7rdu7d6+mTp0ad5SwatWqafTo0apRo8Z+2x577DGVL19effr00Y4dOzRgwAANGDAg4fdw//33q1q1arr99tu1efNmvf3223r77beTeg+5qVGjRvryyy/VqVMn/fbbb3E/m1Afk0jHHnusZs2apcsuu0xz5szRjz/+qH//+98xz5eWlqa0tLRcih6pJuWSkUQ5r05vjJnNlTcPyVgza+Sc+zvYyAAAean8IZ8HHUKhUrx4ca1YsUKzZ8/W5MmTNXPmTP3yyy9auXKltmzZohIlSqhChQo66qij1LJlS11++eU6/vjjYx6vevXq+vTTTzVt2jSNGDFC06dPDx+rdOnSqlWrlho2bKjTTz9dF154YZ5/221meuCBB9SjRw8NGjRIn3/+eXjUr6JFi6pKlSo6+uij1apVK3Xo0CHDrO6S9MILL+iiiy7SlClT9M0332jVqlVas2aN9uzZo0qVKqlBgwbq0KGDrrrqqrgjiN1zzz267LLL9PLLL+vTTz/VkiVLtHnzZlWpUkU1a9ZUixYt1KlTp5iztPfs2VMXXnihXn/9dU2aNEnz58/X+vXrJXl9S+rVq6cWLVronHPOUZs2bXLvA4yhcePG+vnnnzVs2DB98MEHmjNnjtatW6eSJUuqZs2aOvroo9WxY0ddcMEFUfc/8sgj9d1332ncuHEaPXq0Zs6cqdWrV2vbtm0qW7asDj30UDVq1Ejt2rXTBRdcoAoVKuT5e0IwLF47vYLCzPpJelDS4865B4OOJ5WYWQNJ80Ov58+frwYNGuTb+Rd2D27IvvpDC/61jeRwvRUse/bs0eLFizOsO/LII1WsWIH9ngwAApGTv6cLFizIPGx0Q+fcgtyNMLZU7MCeHZP95wvjlgIAAACQMgpLMrLJfz4s0CgAAAAAJKywJCNH+s+0kwAAAAAKiAKfjJhZSUn3+S8XBRkLAAAAgMSlXC9BMzs0gWJFJFWS1FTSbZKOlVcrMjIPQwMAAACQi1IuGZG0VNlrbvW9vDlHAAAAABQAqdpMy5J47Jb0hqR2zrkdgUQLAAAAIGmpWDPSL4Ey6ZI2S1oiabpzbn3ehgQAAAAgt6VcMuKcSyQZAQAAAFDApWozLQAAAACFHMkIACAQZrbfOueYLgoAkpWenr7fumh/Y1MRyQgAIBBFiuz/L2jXrl0BRAIABduOHfuP4RTtb2wqSrk+I2bWLS+O65x7My+OCwDIHjNTqVKltH379vC6TZs2qUyZMgFGBQAFi3NO69dnHMupdOnSBaZmJOWSEUlDlL15RuJxkkhGACDFlClTZr9kpGjRoqpQoYKKFi0aYGQAkNrS09O1Y8cOrV+/PsPfUUkqV65cQFElLxWTEcmbPwQAUMgddNBBWrduXYZ169ev3+9bPgBAYsxMZcuWDTqMhKViMnKavLgek9Rc0npJIyXNlLRaXqJSTVILSZ0lVfa39ZG0N4B4AQDZVKpUKVWtWlVr164NOhQAKPDMTLVq1VKxYql4ix9dykXqnJtmZqMlNZP0rqTrnXP/RCn6ppndJ+kVSZdKusk51ykfQwUA5ILKlStr165d2rRpU9ChAECBFUpEDjrooKBDSUrKJSNmdpmkiyTNcM5dHq+sc26zmV0h6TBJF5pZF+fcW/kRJwAgd5iZatSooYoVK2rDhg3avHlz1GEqAQD7K126tMqVK6eyZcsWqBqRkFSM+Cp5Hc7/m0hh55wzsxcltZTUUxLJCAAUMKGRtUqVKiXnnJxzJCQAEIeZqUiRIgVm1KxYUjEZOc5//iWJfRZl2hcAUECZWfifLACgcEvFv/QV/OeKSewTKls+d0MBAAAAkFdSMRlZ7T9fmMQ+obKr45YCAAAAkDJSMRn5VN7wvdeb2blZFTazcyRdL6+fyad5HBsAAACAXJKKycgzknZKKirpAzN71cyam1m4f4uZFfPXvSrpQ3l9X3b5++YLMzvCzG42s3fN7Ccz22xmu8xsjZlNNbNeZlYtv+IBAAAACpqU68DunFtsZldLGiovIbnaf+wxs43yakAqal/sJild0jXOucX5EaOZDZHUPcbmqpLa+I/7zOxm59zw/IgLAAAAKEhSLhmRJOfcW2a2StIgSYf7q4vLu9HP7HdJ1zrnJudXfJJq+c/b5NXMTJY3+tdmSbXlzQx/haRy8iZn3O2cezcf4wMAAABSXkomI5LknJtiZkdJOlfS2ZKOl1RZXk3I35J+kDRJ0kfOufwejP5PSf+W9Lpzbkumbd9L+tDMPpb0trx4XzKzcc65bfkbJgAAAJC6UjYZkSQ/yRjnP1KGc65bAmXeMbNLJXWUVEXSGfJqUQAAAAAoNTuwFyafRywfFVgUAAAAQAoiGclbJSKW9wYWBQAAAJCCUjoZMbOWZjbYzBab2RYz22Nm9TOVaWVmV5lZp6DijOP0iOUFgUUBAAAApKCU7DPizykyQN6QvpLXCVzyhvXNrKSk1ySlm9kRzrll+RBilsysmaRz/Jd/SpqaC8c8WNFHFIvniJyeFwAAAMgLKZmMSBoo6Sp5ScgqSV9LujhaQefcZDNbIm8I4I6SXsinGGMys3KShmhfzVMv59yuXDj0TZIezoXjAAAAAIFLuWTEzFrLqxFxkv4j6T7n3B4zizd87/uS7pE30WCgyYiZFZU3pO+x/qqRzrmhAYYEAMhlC7tb1oXyUP2h0RoKAEDBk3LJiKRr/ecvnHN3JbjPt/5z/bil8piZmaTXJf3LX/Wd9jU1AwAAABAhFZORk+XViryaxD4r/OeauR9OYvxEZKCk7v6qHySdFWVSxJwYIGlUkvscIWlsLsYAAAAA5IpUTEZq+M8/J7HPdv+5ZC7Hkoz/SrreX54v6Qzn3PrcPIFzbo2kNcns4+VIAAAAQOpJxaF9Q31Dyiaxz8H+86ZcjiUhZvZfSTf7LxdIOt05ty6IWAAAAICCIhWTkVX+c50k9mntPy/P3VCyZmYvSrrFfxlKRNbmdxwAAABAQZOKych0eUP6Xp5IYX8Y3evk9TP5Ig/jinbuFyTd6r8MJSJJNaMCAAAADlSpmIyM8J/bm1n7eAXN7CBJoyVV81cNycO4Mp/7P5Ju818uFIkIAAAAkJSU68DunPvMzCZKOlvSB2bWX968HSElzKy+vOFzb5c3gpaT9LZz7of8iNHMnpb0b//lOnm1Iwf7M6THssE5tzKvYwMAAAAKipRLRnxdJE2TdJykB/1HaIan2ZnKmr/ueuWfSyOWq0j6PIF9hkrqkSfRAAAAAAVQKjbTknNuo6QWkl6WtFtewhHtkS7pNUltnHNbAwkWAAAAQLakas2InHPbJd1kZo9KOk9Sc3lD+BaTtFZebchY59zSAGKrk9/nBAAAAAqblEtGzOx4f3G9c26Fc26VvNnYk5mRHQAAAECKS8VmWnMlfS/p3IDjAAAAAJCHUjEZ2e4/zwk0CgAAAAB5KhWTkdAM7KkYGwAAAIBckoo3/F/6zycFGgUAAACAPJWKycjL8uYU+beZlQ46GAAAAAB5I+WSEefcN5IeklRP0idmVi/gkAAAAADkgVQc2vchf3GupFaSfjazmfJG2FovaW+8/Z1zj+RpgAAAAAByRcolI5L6ymumJf+5iKSW/iMRJCMAAABAAZCKyYgkWRavY3FZFwEAAACQClIxGakbdAAAAAAA8l5gyYiZTZZXk3GVc25ZxKZQ7cZK51zc/iEAAAAACq4ga0bayks8ymRav1RSuqTjJS3M35AAAAAA5JdUGNo3Wn+QRPuIAAAAACiggkxGtvjPVQOMAQAAAEBAgkxGlvjPV5tZKtTQAAAAAMhHQfYZGSepkaQuks4wsyWSdkVsH2xmW5M8pnPOtcutAAEAAADknSCTkf6SLpTUQFI1SQdHbDNJTZI8nol5RgAAAIACI7BkxDm3xcxaSLpZUntJtSSlSTpMXlKxStLuoOIDAAAAkLcCnfTQObdVXg1J/9A6M0v3F9s75xjaFwAAACik6DgOAAAAIBCB1ozE0NN/XhFoFAAAAADyVMolI865oUHHAAAAACDv0UwLAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEgmQEAAAAQCBIRgAAAAAEoljQASB3mFkPST2ibCqdv5EAAAAAiSEZKTzqSGoTdBDw2F3jAju3e+68wM6NYAR5vW288/nAzl3+kM8DO/eBjOvtwGKD7w7s3K7ns4GdG/mHZKTwWCppWpT1pSU1zd9QAAAAgKyRjBQSzrkhkoZkXm9mDSTNz+94AAAAgKzQgR0AAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGAAAAAASCZAQAAABAIEhGcoGZlTWze81sppmtM7PtZva7mb1pZqcGHR8AAACQiooFHUBBZ2YnSnpf0mGZNtXxH13N7L+SbnfOufyNDgAAAEhd1IzkgJnVkTRB+xKRNySdJam5pBsk/eavv1XSU/kdHwAAAJDKSEZy5jlJB/vLtznnrnbOfeKcm+Wce0VSU0mL/e13m9nxgUQJAAAApCCSkWwys2MkXeS/nOGc+2/mMs659fJqRSTvs74/n8IDAAAAUh7JSPZ1jlh+NU65TyQt85fPM7NSeRcSAAAAUHCQjGRf24jlybEK+Z3WQ9sPknRSHsYEAAAAFBiMppV9DfznLc655VmUXRixXF/SV9k5oZkdLKlqkrsdE/ni119/zc6ps23Jxnw9XQZuwYLgTr4uq0si7ywI8n0HjOst//30y5bAzl12Y3CfeZDXmsT1FoQgr7dArVwX2KkP5P9n+SnKvWGJ/Dy/Mdps8swsTdIO/+VC51yDLMpfIuld/+UTzrne2TxvX0kPZ2dfAAAAIAEXOOc+zK+T0Uwre8pGLG9OoHxkmbIxSwEAAAAHEJKR7InshL4rgfI7I5ZL53IsAAAAQIFEn5Hs2R6xnEi7urSI5W05OO8ASaOS3OcgSU0k/SNpk6Q/lFgChew5QtLYiNcXSFoSUCwo/LjekJ+43pBfuNbyVwlJtSNeT8vPk5OMZE+yza4OirFvUpxzayStycau32T3nEiOmWVetcQ5Rw885AmuN+QnrjfkF661QHwf1IlpppUNzrmd2pcU1I5X1ndYxHJww5AAAAAAKYRkJPvm+89lzezQLMrWj1gmswcAAABEMpITke3pTo9VyLy6xtP8l1slzc7LoAAAAICCgmQk+0ZHLF8Xp1x7SXX85Y+cc9vjlAUAAAAOGCQj2eScWyjpA/9lSzO7JXMZM6sk6b/+y3RJT+VPdAAAAEDqYzStnLlTUmtJVST918waS3pb3hC6J0i6X9Lhftn/OOfmBhAjAAAAkJJIRnLAOfe7mZ0j6X15o2pd5T8yGyDpnvyMDQAAAEh1JCM55Jz7zswaSrpR0kWS6kkqI2m1pK8kDXLO5evkMQAAAEBBQDKSC5xz/0h62n8AAAAASAAd2AEAAAAEgmQEAAAAQCBopgXkrrWS+mV6DeQVrjfkJ6435BeutQOIOeeCjgEAAADAAYhmWgAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkIwAAAAACQTICAAAAIBAkI0A2mVlxMzvfzPqZ2UdmttDM/jKzXWa22cx+MbMRfhkLOl4UXmZW18y2mJnzH0uDjgkFn5ktjbimsnr8O+h4kfrMrKiZXWZmo8zsNzPbZmb/mNmvZjbOzO4ys5pR9kszs3PM7Bkzm2ZmqyP+1y4ys2FmdlYQ7wk5Z865oGMACiQzqyXpjwSLfy2pk3NuVR6GhAOQn+hOltQ2YvUy51ydQAJCoeEntYclWPwO59zzeRcNCjozayxpiKTjsyja0zk3JGK/4yV9KalcAqf5RNIVzrl12QwTASgWdABAAbde0jRJ30n6XdJqSRsklZfUSNI1ko6T1ErS52Z2onNuR0CxonC6WV4i8pekasGGgkLqO0k9syjDFy2IycxaSZogL6HYK+kdSeMkLZV3L1pbXpJyUZTdy2lfIjLP32+mvGuuuKTmkv4tL3FuL+lTM2vhnNuZN+8GuY2aESCbzKyIJDnn0uOUKSbpfUnn+atuds4NyIfwcAAws8Ml/SipjKROkkb7m6gZQY5F1IxMc861DTYaFFRmVlnSfEnVJa2V1ME5NytO+RLOuV0Rr1tJukfS486572LsU07Sp5Ka+avuc871z6W3gDxGnxEgm5xz6fESEb/MHklPRqw6LW+jwoHCb571hrxE5H3n3HsBhwQA0TwtLxGRpC7xEhFJikxE/NdfO+cujJWI+GX+kXRdxKrLshss8h/JCJD3/olYTqTNK5CIWyS1kdcs8OaAYwGA/ZhZFUlX+i8/c859llfncs79IOlv/+VReXUe5D76jAB578qI5Z8DiwKFhpkdoX01bnc651YHGQ8AxHCxpDR/+e3QSjNLk1RTXv+R1ZlrQ3KghP+8N5eOh3xAzQiQy8ysiJlVN7O2ZjZC0v3+pl2SXg4wNBQCmZpnfRI56gyQR44xs6/NbIM/nOpfZjbVzPqY2cFBB4eU1ipiebqZnWhmH0raLOk3Scsk/WNmn5rZ2Tk5kZk1kVTWf7kgJ8dC/qJmBMgFflX02jhFNslrK/tTPoWEwutWSadK2qKMbaSBvFJNGUdqO9h/tJHUy8z+7Zx7NZDIkOoa+s9OUmt5X8hlvvdMk3SGpDPM7H+SbnXZG13p4YjlEdnYHwGhZgTIW07S/0k62jk3PuhgULBlap71gHNuWZDxoNDbLWmipDslnSnpRHnfdF8tbz4HSSol6RUzuyOQCJHqKvnPTtJAeYnIS5IayEtCqku6Qd4w+ZLX/61Xsicxs2sldfBfLpI0KPshI78xtC+QC8ysqKRj/ZdFJVWW1FLeN9eHyPuHfoNzbkUwEaKg85tnTZN0iqSvJJ2a+dtDMwu9Zmhf5JiZVXTObYizvYukN+X9zdsjqaFz7pf8ig+pz8zWyft/GHK/c+7pKOUaSfpGXoKyXVId59yaBM9xqqTP5M05slNSS+fc9zmNHfmHmhEgFzjn9jrn5vuPH5xzk51zj8v79udzSedK+s7Mjo1/JCCm2+QlIjskXZPNZgxAwuIlIv72t7Svpq6YvGsUiLQ9Ynm5pGejFfJHwgo19Sslb96kLPmzun8oLxFJl3QViUjBQzIC5CHn3BZJXeX9Qa4mr5oaSIqZ1ZP0hP+yH98+I4VETuLaLrAokKo2RyxPcM7FG+Xqo4jl5lkd2MyOkzfRYXl5zcCu9xNkFDB0YAfymHNujZl9Ja/NdRszq+GcWxV0XChQ7pFUWt5ACCvMLKsJvcpElNnqnBuXp9HhgOWcW+U3xakiqVbQ8SDlLNO+JszLEygbEneUNjMLtTqoLC8Ruck591p2g0SwSEaA/LEuYrmOJJIRJCM0Tn95ScMSKF9F+8b0XyaJZAR5iSaDiGWepNCQvUWzKBu5PWYNipnVlzRZUlV/1S3OOYbNL8BopgXkj0MiljfHLAUABYiZVZeX/ErSyiBjQUqaErF8ZBZlI7dHHezFT0SmaF/NyS3OuQHRyqLgIBkB8piZ1ZE3spYkbZO0JLhoUBA553o45yyrR8QuyyLW1wkqbhwQbpIUuvamxCuIA9Jnkv72l881s9Jxyl4asbzfteQPADNZ+xKRW51z/8uVKBEokhEgm8zsCjOrnEWZqpJGyhvpQ5Leds5tj7MLAATOzM4xszJZlOmifXNC7JU3fwQQ5pzbrX0jrlWS9KI/THkGZna+pFA/txWSPsi0/Rh5iUho8s3bnHNcb4UE84wA2WRmU+WN+PGxpKmSFkraIK8v1iHyZifuLqmiv8uv8sY/X5f5WEBuYJ4R5Bb/79sJ8kY4+lLSL5I2Siopr0PypZLOitill3PuSQGZmFkJeTUdrfxVU+SNLLlE3v/HjvImPiwmL6k91zk3KWL/w+XNrVTDX/WypERqRH7xkyGkOJIRIJv8f9ZtEiw+XtLVzrnVeRcRDnQkI8gtSfx92yrpHuccw5YjJjOrKGm0pNPjFNssqbtzbkymfXtIGpyN09Z1zi3Nxn7IZ4ymBWTfFZLayvuHfYK86uOq8kYE2SSvJuQbSe84574JJkQAyJa7JJ0mr/b3GHlDqFaW98313/JGSfpc0hDn3N+xDgJI4Qk025lZZ0lXSjpJ3v/LbZIWS5og6SXn3NrgokRQqBkBAAAAEAg6sAMAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjAAAAAAIBMkIAAAAgECQjABAHjGzqWbmzGxp0LHEYmZt/RidmfUIOp7cYmZlzGy5/76uCToeoCAws4pmttH/vekYdDw4MJCMAAAKowcl1Za0WNKQYEOBmfVNNuk1szoR+0yNUaaon1A/amafmdlKM9tpZlvN7HczG2VmncysaDbj/iQihtVmViyJfUuZWUszu8XMhpjZfDPbE3G8ttmJKc75KplZezN7wMzeM7NlEedyiRzDObdB0nP+y+fNrFRuxghEk/AvFQAkI+Kf31DnXI8gY8GBxcxqS/q3//Jh59yeAMNBHjGzqpIWSqoSZXMJSXX8RydJ35nZpc6535I4/qGS2kWsqibpXEljEzzEW5I6Jnq+XLBEUoVcOM7zkm6TdJikWyX1z4VjAjFRMwIABzDn3FTnnPmPIUHHk0t6S0qT9JukdwKOBXknTfsSkd/kfaN/saTm/uN6SfP97U0kTTGzaIlLLD217z5pq/98dRL7W8TyNklfy0sY8krk+dZJmiRpS7IHcc5tlvSS//I+MyuTC7EBMZGMAAAKDTM7WFJ3/+Uw51xCzVNQIDlJUyW1c84d4Zy72zn3vnNulv94VdJJ2leTcaikRxI5sJmZvGREkqZLGuovn2NmNRKMb4ykHpKOk1TOOXeypK8S3Dc7npGXjNVxzlV1zp0t6e9sHutNeZ9vJSWXgAFJIxkBABQmPSWVlHcj9WbAsSAPOedWOudOc85NjlNml6RrJe3yV13qJxpZOUNeMyXJ63M02F8uJqlbgvEN9R/znXN7E9knJ5xzj/vJ2LJcONbv2pc4XZ/T4wHxkIwAhYCZ9cjcKdLMzjGzsWb2p9+h80+/M2fzBI9pfsfPd81sqZltM7MtZvazmb1iZsfH2G9pps6S3SM7UUY86vjlZ/ivf4kTy28R+z0Qo8ypEWWi3iyYWTEz62lmH0d8LuvNbLaZPWFmNbP4TDKMjuV3GH3I3/9vf9vz8Y4R5ZgHm9msiNj7JLN/xHEOMbPH/GNtMLPd/nv7xcw+NbPeZnZClP1ijqaVaVuij7Yx4svW9ZQNXf3n7+L1D4jxO9PJzCaa2Soz2+5/ds/6tS2R+x5hZi/4sW/zP+eJZtYu2rki9itqZmeYWX8z+8K8DtG7/M/hVzMbYWZnJPpG/ev5Sv/3+nfzOm3/Y2Y/mdmHZnaNRWmWFPG+h/ivjzWzl/z3u9mijKSU09+dIDnn1kqa57+sJKlyAruFagO2SRrpnPtO+5p8XZW7EaasUBPH+mZ2UqCRoHBzzvHgwaOAP+Q1BXD+4zRJL0e8zvzYK6lnFserLWlWnGM4SemSHpVkmfZdmsV+oUcdv/xjEetqRYmlTqb9PosRc9+IModE2V5X3g1JvJi2Seoa53OZ6pdbKulESSujHOP5aOVjHO9IeW3Inbxvbrtn8+d/rry24Vl95lOj7Ns2YnuPONsSfbTNzespyc+hbrSfQwK/M+0kjYgT29KI6/UCSZvjlL0hzjmfT/AzfFtSWhbxN5L0SwLHGhJl3/A2eU3adkTZr2Mu/+70jXWdxdmnTrxrN8lr48eIY1XMomyliM9kWMT6OyOO0TqbcQyJ97uS2w9F/D3Oxr6NImLtl9ex8jhwH4ymBRQ+/SSdImmypNfkDW1aSlIHef9Mi0kaYGZfOOf260xpZtUlzZB0iLwbxNGSxkn6XdJuScdLulnSCZL6SNopL6EIaS9vJJvQN5Fj/XKZrfSfP5fX4VjybgqHZiqX+dvmk82spHNuR4xyvzjnVkZuMLNKkqbJuymWpC8lDZS0SN7oM+dLulHe5/Smme1wzo2KEnPIQZI+kPcN64uSPpbXYbSWpISGEDWzFvI+1yrybm4vds59msi+mY5TRd6oPWUkbZf3M/9E0l9+kWryflZny7upSMa38tq7x1NT3jVSVl5CtSpTfDm9npJxWsTyzCT2e0RSK3k/x8F+bFXlNU+5UF5znUFmdr+kUfLeYx95CdZeeddeH0mlJb1gZp85536Ncp5ikv6U9KEf3xJ5HaOrSTpW0i2SDpd0mbzr6dZowZpZI3n9GEIdi6dJGiZvZKmd8n4mzSVdlMX7biKpi6T18hKlr/z968u7ic2L351851+Dx/ovVzlv+Np4rpTXOV7a1zxLkoZLelrez/Fq5W3/j1QwX96XHAdJOl3Sw8GGg0Ir6GyIBw8eOX8o47e8TtJTMcpdG1HmmRhlxvvb/5bUJEaZYpJGat83+nWilIn5zWymcmnybqIzfAsZsf0tf9uH/rmcvA6rkWXKRGz7X5RjvBERzyuK8u27vH+2O/0yGySVj1JmasRxdkg6OYv3Fiq/NNP6C+R9k+zk3Zw2ysHPvmdETBdmUbZylHVtI/bvkeS5y0r6XvtqNq7Iq+spwXj+F/Fe6iX5O9M3ShmTlziFyvzlv9/9vlmX1DmB360jJBWLE1NRebUiTtIeSYfF+KwWR5zrtize56FR1kW+718lVY+zf2797vSNOE5vSQ0TeLSP2GdqDn5HXsrqZ5Op/NzQ723m9yvvSwgn7ya9bDZiGRIRS9vsvqckzrc0dL5s7j/N33+rpCJ5HS+PA/NBnxGg8PlRUq8Y297Qvm/M22beaGZNJJ3jv7zXee2k9+O8eRtulHfjWFz7Rp1JmnNup7xveaX9a0Ek70ZH8mpYZsUod6ofh+TVCIWZNxfBFf7L3+XdvO1XQ+C8TrD/57+soKzbhf+fc256FmX2Y2Y3SXpf3jfJCyW1cM79kOxxIkSO7PN5vILOueyOrLMf8yZ/GymvRkOS+jjnRmQqk9/X05ERy6uT2O97eTWKmeNykl6IWHWwvCaO0b5ZHy1phb98WpTtcs4tcXHmPHFeJ+db5dW2FFX0OSq6SKrnL7/hnHsx1vH8Yy6Pt13STc65qJ9VHv7uPCav5jSrx6QsjpMlMztL3rUleQnxU1mUP0le8yTJmyMp8/sN1ZSUkXRpTuMrAEL/L0rLq3EDch3JCFD4DHPOpUfb4N/szPZf1otS5GL/OV3Su/FO4t/Yhppitc5GnJFCN9E1zCzUnEJm1lBeExZJ+iyiXOZkJPQ6XdKUTNtOk9dsTJIG+8lPLP+LWD4ri5iHZLE9A78D95P+OYpI+kJeu/OsbhazEtkk7bocHisZA+U1/ZKk15xzT0Qpk9/XUygx2+2cS2Z+hRHRbrJ930csz3fOzY1WyN8/tO2IRE5qZmXM7DAzq29mDf3rvbq8ZlOS1y8ps44RyzmdjO5PSfGaBubV706+MLNj5HXCLiLv2/3uCSTkoY7rTvs3GZW8pnxrMpUtzCI/r0SHNAaSQjICFD4Ls9ge+udSLsq20EhbRSSFRtWJ+ZA3hr+U839Skd/ot4uyvMR5w1WGyp1kZuWjlJvrnFuvjCJHafo6XhDOuRWSQslBozhFtzrnFsU7VibF5LXpv99/PVJS+xjfsCfrA0lr/eVnzGyumfUzszPNLNrPOMfMG/HrGv/lJO375jmz/L6eSvvPyX6uP8fZFnmsn7I4TqhszM/dzOr6I3EtkddXaKmkBcpYI1DVLx5tgr7QZ7TaORdzBLoEzY2ThEl587sjebVLltVDXsf5bDGzw+T1nargr3rAOfdxFvuUknS5//ILF2U0Nr9ma7j/soWZ1c9ujAVE5N9TJj9EniAZAQqfbVlsD9WaRPv9PzjKukSUzrpIXLMlbfKXI4c2DSUZn/nPM+W9v6Lym8L4HbhDNz/RmilFDuOZSNOdUAfseMN/Jnuze4j2NXdZKKlLFt8yJ8w5t0leU6hQctRI0kPybsQ2+EOvPuh34s0xM7tS3qhXkvSDpM5xmh7l9/W0238umeR+W2NtyFTLmJPfLZnZJfISmtvkdVTPar6LaJ9DKFFZGWVbsjIn7pnlxe9OnjOzWvKaa4Y63T/qnHs6gV0v1r7kZUiccoMjlgv7ML+lIpZ3xywF5ACjaQGIFPqb8Lei9CmJY1fWRWJzzu01s2nyRuZpa2ahEana+M+f++V2mdmX8pqBtJNXK3Ca9t3URZv8LPKGL963wMlIdgKz1fKGYW0jb6SiQWZ2TazmdMlyzs32v6E9R96oaSf75ykir6nPiZLuN7PrnXPDYx8pPjM7TdLr/suVks51zm2Os0t+X0/hmgkzKxavf0Z+M7N68iZhTJOX1LwgaaK8DuR/S9oVqqUws+XybqTjJSu5cS1ndR3nxe9OnvLnO5ksL9mTvME8Hkpw98hmV4PNbHDMkvt0M7MHnHOF9UY9MrHcGFQQKNxIRgBEWivpaHnj7C93zv2Tj+f+XF4yUl5eU5Qi8pq7OGXsB/KZ9iUjinjeLW/Y0czWRSzXUNbN2EJNhHKts7e8kYbOkZc8tZfXQTvNzLq5XJqZ2T/OR/5DZlZBXsf+LpIukfct+xAzm+ucmx/rOLH4yc778voQbJb0L5dpCOUo8vt6Wiaphb9cSfva9qeCq7RvuNiLnXMT45StFGfbWkmHyhtGOq+lwu9OwvxEZIr2DWTQ3zkXdZLUKPsern1ffiSjqqTz5P1uFEaR12KOZ3YHoqGZFoBIoc7tpn2dk/NLZBOrM5SxH8i6KOWO9W8+QuVmOueiNbf5MWK5ZbwAzOwQeTd6ktcEKdc457bLS7bG+au6SHrXzIrH3itH59vonPvQOXeZ9s3jUlTZGAHIzGpImiCvCcseSZ2ccz/G3cmT39fTvIjlo/PhfMkIzdeyPl4iYmYNFL9tfmhEsupmltfvMSV+dxLhX6NTJB3lr3rGOXdfEoe4Svtqgh6U13ck3uNK7avBK8xNtY7xn39PclAIIGEkIwAivRexfH8Ob5RDkxKmxS3lc84t0L5hJNtpX5KRuR/IXO375rWn9o0KFmtY2ynad9PQ08xKxCgnSTdFLOd4WNHM/H4iF8ubOE/+8vtmltBnlAOR76VqzFJRmNlB8mpbQjeaNzjnPklw99y8nhIR2cm6ecxSwQi1RCgV0QwxmtuyOM6YiOV7cxZSllLmdyeeKInIs865hD8b/+fRw3+5VtKTzrl3sniMkNcvS5LO9pOxQsXvjxf6+xp3AAMgJ0hGAIQ5577Uvs7ijeW1m455o2xmRczs4hgjyoSa8BwVZVssoT4frfyHIuIJxRjZbOueKPtm4JxbKyk0/8Xhkv4TrZyZtZV0t/9yo7w5WXKd37b8cu0bkaeDpHH+SD5JM7NzzOzQLIr9K2J5SRLHLipvSN7QELOPOedej7NLBrl8PSXiK3lNyKR9zbVSRWjkq1LyvlXfjz84wLVZHOcd7Rus4CozizpLe8Qxa8fbHk+q/e7EOHd1eX8PQrVEzzrn7omzSzRnyRtkQpJGJ9F08h3/uaik7kmeM2Fm1jZi1LmpeXWeKCJ/hybk43lxgKHPCIDMrpT0jaTD5I0A1dLMBvnrNshrQlJH3j+qC+X9Ez9T+7cnnyZvvoUTzewReZMWborYvjRKB+PP5d2oh0ZD2qXo/UA+l9RJXv8SyRsNaWac93S3vKZftSXd5N/sDpTXebiCvDbfN2vfxInX+aNU5Qm/w353eTPPXyvv85tgZh2y0RTiUklX+h37J8lrWrNG3pdNteT9jEIjeW3UviQoEQ9qXyLzlbxmZQ2z2Of3TM3lcut6ypJzbreZjfXP2dbMiqdQx+Ih8mo9TNKrZnaivNnp18r7bC6Xd00vktdOP2oNlnNujz8q13R5n92LZnaRvKGjF8rrn1RDUlN5s8J/p33f+mdHSv3uRDKzyvK+hAglIuMkDc3GNRrZcT3ufDiZjJX3O1xKXs3Rk5FDJfuJUubmiZHzO51tZnUiNzrnhiRx/gzM7ATtm4Q05KCI7T0ybVudRd+l9v7zDvl90YA84VJgGngePHjk7CHvZsP5j7ZZlB0SKhunTHV5k6G5BB57JLWMcoyjJf0TZ786Ufapk6nM1BjxHZmp3IQEPqO68voUxHsv2yR1i3OMqX65pQn+XOKWl3dj+kLE+b+WVD7Jn/2QLN5T6LFa3iSLmfdvG1GmRzaPHfnY7/rLjespic/j9IhjXZhLvzOhckNy8rsl6Q55w//Geu+L/Gt7abzr3z9WY3m1XFl9nvvFnOj7yeXfnb6xrrM4+9SJ2Ge/zyLTtZuta1Re0rfLX79SUpEkr7fREcdtk9P4Ypwj8jjxrom+SZ4v3rGKy0uUnbxJQbP1+8iDRyIPmmkB2I9zbrVz7kx5/wQHyfvGdZO8oUA3y5srYZS8ye4Occ7NiHKMX+SNivWKvEnltsn7xxbvvEslRU40FrUfiHNusfZNsCbFaKKVaZ/f5d3AXS2vycFqeSNwbZQ3y/ZTkuo5597M6li5xXlulxSaA6GlpM/NLN5oSpndLukCSc9LmiHvRnabvBus1fI+wzskHeWc+yp3Ik9OblxPSZxrsvZ1vO6eg7BznXPuP/JGbHpfXv+oPfJGrJopr/9HY//aTuRY38vrXHy1vBqBlfJqRTbJ+zw/kJdw/TsX4k65351c1E37anVGuuSH234nYrkwzcj+L+2bdPP5AOPAAcCci3tvAABAgWJml8q7Sdwt7yZ5eRa7AIhgZuPlDUf+qXOufVblgZygZgQAUNiMlPeNfXF5s9EDSJCZtZKXiDhJvQIOBwcAkhEAQKHivCr/2+T3TzCzZEZ0Aw50j/vPQ51z38UtCeQCRtMCABQ6zrmvzOxaeaNA1da+4XABxGBmFeWNhDhN0v8CDgcHCPqMAAAAAAgEzbQAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAgSEYAAAAABIJkBAAAAEAg/h/23adigiT8/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 840x450 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "\n",
    "z = [x[1] for x in cc_m11_ddr_not_nhej]\n",
    "counter_cc_m11_ddr_not_nhej = Counter(z)\n",
    "\n",
    "z = [x[1] for x in cc_m11_nhej]\n",
    "counter_cc_m11_nhej = Counter(z)\n",
    "\n",
    "z = [x[1] for x in cc_m11_reprogramming]\n",
    "counter_cc_m11_reprogramming = Counter(z)\n",
    "\n",
    "z = [x[1] for x in cc_m11_senescence]\n",
    "counter_cc_m11_senescence = Counter(z)\n",
    "\n",
    "labels = counter_cc_m11_ddr_not_nhej + counter_cc_m11_reprogramming + counter_cc_m11_nhej + counter_cc_m11_senescence\n",
    "# labels.most_common()\n",
    "labels = [x[0] for x in sorted(labels.items())]\n",
    "zero_dict = dict(zip(labels, [0] * len(labels)))\n",
    "\n",
    "dict_cc_m11_ddr_not_nhej = {**zero_dict, **counter_cc_m11_ddr_not_nhej}\n",
    "dict_cc_m11_nhej = {**zero_dict, **counter_cc_m11_nhej}\n",
    "dict_cc_m11_reprogramming = {**zero_dict, **counter_cc_m11_reprogramming}\n",
    "dict_cc_m11_senescence = {**zero_dict, **counter_cc_m11_senescence}\n",
    "\n",
    "counts_m11_ddr_not_nhej = list(dict_cc_m11_ddr_not_nhej.values())\n",
    "counts_cc_m11_nhej = list(dict_cc_m11_nhej.values())\n",
    "counts_cc_m11_reprogramming = list(dict_cc_m11_reprogramming.values())\n",
    "counts_cc_m11_senescence = list(dict_cc_m11_senescence.values())\n",
    "\n",
    "x = np.arange(len(labels))  # the label locations\n",
    "width = 0.20\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "rects1 = ax.bar(x - width*2, counts_m11_ddr_not_nhej, width, label='DDR not NHEJ', align='edge', color='#d55e00')\n",
    "rects2 = ax.bar(x - width, counts_cc_m11_nhej, width, label='NHEJ', align='edge', color='#009e73')\n",
    "rects3 = ax.bar(x, counts_cc_m11_reprogramming, width, label='Reprogramming', align='edge', color='#0072b2')\n",
    "rects4 = ax.bar(x + width, counts_cc_m11_senescence, width, label='Senescence', align='edge', color='#f0e442')\n",
    "\n",
    "ax.set_ylabel('frequency')\n",
    "ax.set_xlabel('network size (macroH2A1.1)')\n",
    "# ax.set_title('Scores by group and gender')\n",
    "ax.set_xticks(x)\n",
    "ax.set_xticklabels(labels)\n",
    "ax.yaxis.set_minor_locator(MultipleLocator(1))\n",
    "ax.legend()\n",
    "\n",
    "plt.savefig(\"FigureXX_B.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Create network (sif) and attribute files for cluster (size 62) of macroH2A1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Graph successfully exported to SIF at full path:\n",
      "/scratch/tom/playgrounds/g_m11_62.sif\n"
     ]
    }
   ],
   "source": [
    "cc_m11 = g_m11.clusters()\n",
    "cc_m11_62 = [c for c in cc_m11.subgraphs() if c.vcount() == 62][0]\n",
    "cc_m11_62_nodes = pd.DataFrame(cc_m11_62.vs[\"name\"])\n",
    "cc_m11_62_nodes = cc_m11_62_nodes.rename(columns={0: 'Node'})\n",
    "cc_m11_62_nodes[[\"NHEJ\"]] = [c in gene_lists['NHEJ'].tolist() for c in cc_m11_62_nodes['Node']]\n",
    "cc_m11_62_nodes[[\"DDR_not_NHEJ\"]] = [c in gene_lists['DDR_not_NHEJ'].tolist() for c in cc_m11_62_nodes['Node']]\n",
    "cc_m11_62_nodes[[\"Reprogramming\"]] = [c in gene_lists['Reprogramming'].tolist() for c in cc_m11_62_nodes['Node']]\n",
    "cc_m11_62_nodes[[\"Senescence\"]] = [c in gene_lists['Senescence'].tolist() for c in cc_m11_62_nodes['Node']]\n",
    "\n",
    "diff_expr = [x['Gene Symbol'] for index, x in epi_vs_m11.iterrows() if \n",
    "             x['pvalue']<0.05 and (x['log2FoldChange'] > 1 or x['log2FoldChange'] < -1)]\n",
    "cc_m11_62_nodes[[\"DEG\"]] = [c['Node'] in diff_expr for index2, c in cc_m11_62_nodes.iterrows()]\n",
    "\n",
    "cc_m11_62_nodes.to_csv(\"g_m11_62_nodes_attributes.tsv\", index = False, header=True, sep='\\t')\n",
    "cc_m11_62_nodes.head()\n",
    "\n",
    "from pyntacle.io_stream.exporter import PyntacleExporter\n",
    "PyntacleExporter.Sif(file='g_m11_62.sif', graph=cc_m11_62, header=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create network (sif) and attribute files for cluster (size 10) of macroH2A1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Graph successfully exported to SIF at full path:\n",
      "/scratch/tom/playgrounds/g_m12_10.sif\n"
     ]
    }
   ],
   "source": [
    "cc_m12 = g_m12.clusters()\n",
    "cc_m12_10 = [c for c in cc_m12.subgraphs() if c.vcount() == 10][0]\n",
    "cc_m12_10_nodes = pd.DataFrame(cc_m12_10.vs[\"name\"])\n",
    "cc_m12_10_nodes = cc_m12_10_nodes.rename(columns={0: 'Node'})\n",
    "cc_m12_10_nodes[[\"NHEJ\"]] = [c in gene_lists['NHEJ'].tolist() for c in cc_m12_10_nodes['Node']]\n",
    "cc_m12_10_nodes[[\"DDR_not_NHEJ\"]] = [c in gene_lists['DDR_not_NHEJ'].tolist() for c in cc_m12_10_nodes['Node']]\n",
    "cc_m12_10_nodes[[\"Reprogramming\"]] = [c in gene_lists['Reprogramming'].tolist() for c in cc_m12_10_nodes['Node']]\n",
    "cc_m12_10_nodes[[\"Senescence\"]] = [c in gene_lists['Senescence'].tolist() for c in cc_m12_10_nodes['Node']]\n",
    "\n",
    "diff_expr = [x['Gene Symbol'] for index, x in epi_vs_m12.iterrows() if \n",
    "             x['pvalue']<0.05 and (x['log2FoldChange'] > 1 or x['log2FoldChange'] < -1)]\n",
    "cc_m12_10_nodes[[\"DEG\"]] = [c['Node'] in diff_expr for index2, c in cc_m12_10_nodes.iterrows()]\n",
    "\n",
    "cc_m12_10_nodes.to_csv(\"cc_m12_10_nodes_attributes.tsv\", index = False, header=True, sep='\\t')\n",
    "cc_m12_10_nodes.head()\n",
    "\n",
    "from pyntacle.io_stream.exporter import PyntacleExporter\n",
    "PyntacleExporter.Sif(file='g_m12_10.sif', graph=cc_m12_10, header=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select largest components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl = g_epi.components()\n",
    "cl_sizes = cl.sizes()\n",
    "giant_component_index = cl_sizes.index(max(cl_sizes))\n",
    "lc_epi_names = [g_epi.vs[cl[x]][\"name\"] for x in set(cl.membership) if giant_component_index == x][0]\n",
    "\n",
    "cl = g_m11.components()\n",
    "cl_sizes = cl.sizes()\n",
    "giant_component_index = cl_sizes.index(max(cl_sizes))\n",
    "lc_m11_names = [g_m11.vs[cl[x]][\"name\"] for x in set(cl.membership) if x == giant_component_index][0]\n",
    "\n",
    "cl = g_m12.components()\n",
    "cl_sizes = cl.sizes()\n",
    "giant_component_index = cl_sizes.index(max(cl_sizes))\n",
    "lc_m12_names = [g_m12.vs[cl[x]][\"name\"] for x in set(cl.membership) if giant_component_index == x][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "## m11\n",
    "cc_m11_nodes = pd.DataFrame(index=lc_m11_names, columns=[\"Node\"])\n",
    "cc_m11_nodes[\"Node\"] = lc_m11_names\n",
    "cc_m11_nodes[[\"NHEJ\"]] = [c in gene_lists['NHEJ'].tolist() for c in cc_m11_nodes['Node']]\n",
    "cc_m11_nodes[[\"DDR_not_NHEJ\"]] = [c in gene_lists['DDR_not_NHEJ'].tolist() for c in cc_m11_nodes['Node']]\n",
    "cc_m11_nodes[[\"Reprogramming\"]] = [c in gene_lists['Reprogramming'].tolist() for c in cc_m11_nodes['Node']]\n",
    "cc_m11_nodes[[\"Senescence\"]] = [c in gene_lists['Senescence'].tolist() for c in cc_m11_nodes['Node']]\n",
    "\n",
    "diff_expr = [x['Gene Symbol'] for index, x in epi_vs_m11.iterrows() if \n",
    "             x['pvalue']<0.05 and (x['log2FoldChange'] > 1 or x['log2FoldChange'] < -1)]\n",
    "cc_m11_nodes[[\"DEG\"]] = [c['Node'] in diff_expr for index2, c in cc_m11_nodes.iterrows()]\n",
    "\n",
    "cc_m11_nodes.to_csv(\"cc_m11_nodes_attributes.tsv\", index = False, header=True, sep='\\t')\n",
    "\n",
    "\n",
    "## m12\n",
    "cc_m12_nodes = pd.DataFrame(index=lc_m12_names, columns=[\"Node\"])\n",
    "cc_m12_nodes[\"Node\"] = lc_m12_names\n",
    "cc_m12_nodes[[\"NHEJ\"]] = [c in gene_lists['NHEJ'].tolist() for c in cc_m12_nodes['Node']]\n",
    "cc_m12_nodes[[\"DDR_not_NHEJ\"]] = [c in gene_lists['DDR_not_NHEJ'].tolist() for c in cc_m12_nodes['Node']]\n",
    "cc_m12_nodes[[\"Reprogramming\"]] = [c in gene_lists['Reprogramming'].tolist() for c in cc_m12_nodes['Node']]\n",
    "cc_m12_nodes[[\"Senescence\"]] = [c in gene_lists['Senescence'].tolist() for c in cc_m12_nodes['Node']]\n",
    "\n",
    "diff_expr = [x['Gene Symbol'] for index, x in epi_vs_m12.iterrows() if \n",
    "             x['pvalue']<0.05 and (x['log2FoldChange'] > 1 or x['log2FoldChange'] < -1)]\n",
    "cc_m12_nodes[[\"DEG\"]] = [c['Node'] in diff_expr for index2, c in cc_m12_nodes.iterrows()]\n",
    "\n",
    "cc_m12_nodes.to_csv(\"cc_m12_nodes_attributes.tsv\", index = False, header=True, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Compute ANOVA among all conditions (m11, m12, episome)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Gene Symbol</th>\n",
       "      <th>gene_biotype</th>\n",
       "      <th>EPI</th>\n",
       "      <th>3_I</th>\n",
       "      <th>3_II</th>\n",
       "      <th>1_1_1</th>\n",
       "      <th>1_1_2</th>\n",
       "      <th>1_1_3</th>\n",
       "      <th>2_I</th>\n",
       "      <th>2_II</th>\n",
       "      <th>2_III</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TBC1D3L</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>7.148288</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>36.919539</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DNAH5</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>21.165576</td>\n",
       "      <td>11.080615</td>\n",
       "      <td>12.798294</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.395609</td>\n",
       "      <td>0.751568</td>\n",
       "      <td>4.215694</td>\n",
       "      <td>2.008978</td>\n",
       "      <td>2.188068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SAP25</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>7.201708</td>\n",
       "      <td>7.178799</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.870301</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>EYA1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>25.060854</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.741551</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.318718</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>CYP39A1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.528252</td>\n",
       "      <td>15.767867</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>15.017370</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Gene Symbol    gene_biotype        EPI        3_I       3_II  1_1_1  \\\n",
       "0      TBC1D3L  protein_coding   7.148288   0.000000  36.919539    0.0   \n",
       "4        DNAH5  protein_coding  21.165576  11.080615  12.798294    0.0   \n",
       "7        SAP25  protein_coding   0.000000   7.201708   7.178799    0.0   \n",
       "10        EYA1  protein_coding   0.000000  25.060854   0.000000    0.0   \n",
       "11     CYP39A1  protein_coding   0.000000   0.528252  15.767867    0.0   \n",
       "\n",
       "       1_1_2     1_1_3       2_I       2_II     2_III  \n",
       "0   0.000000  0.000000  0.000000   0.000000  0.000000  \n",
       "4   0.395609  0.751568  4.215694   2.008978  2.188068  \n",
       "7   0.000000  0.000000  8.870301   0.000000  0.000000  \n",
       "10  0.741551  0.000000  0.000000   0.318718  0.000000  \n",
       "11  0.000000  0.000000  0.000000  15.017370  0.000000  "
      ]
     },
     "execution_count": 324,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m11_m12_epi.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 357,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/t.mazza/miniconda3/envs/playgrounds/lib/python3.7/site-packages/scipy/stats/stats.py:3709: F_onewayConstantInputWarning: Each of the input arrays is constant;the F statistic is not defined or infinite\n",
      "  warnings.warn(F_onewayConstantInputWarning())\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>Gene Symbol</th>\n",
       "      <th>gene_biotype</th>\n",
       "      <th>EPI</th>\n",
       "      <th>3_I</th>\n",
       "      <th>3_II</th>\n",
       "      <th>1_1_1</th>\n",
       "      <th>1_1_2</th>\n",
       "      <th>1_1_3</th>\n",
       "      <th>2_I</th>\n",
       "      <th>2_II</th>\n",
       "      <th>2_III</th>\n",
       "      <th>fvalue</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>adj_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TBC1D3L</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>7.148288</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>36.919539</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.06819</td>\n",
       "      <td>0.196896</td>\n",
       "      <td>0.655542</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DNAH5</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>21.165576</td>\n",
       "      <td>11.080615</td>\n",
       "      <td>12.798294</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.395609</td>\n",
       "      <td>0.751568</td>\n",
       "      <td>4.215694</td>\n",
       "      <td>2.008978</td>\n",
       "      <td>2.188068</td>\n",
       "      <td>21.1798</td>\n",
       "      <td>0.00107409</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SAP25</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>7.201708</td>\n",
       "      <td>7.178799</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.870301</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.29337</td>\n",
       "      <td>0.332655</td>\n",
       "      <td>0.732052</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>EYA1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>25.060854</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.741551</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.318718</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.18013</td>\n",
       "      <td>0.361688</td>\n",
       "      <td>0.705583</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>CYP39A1</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.528252</td>\n",
       "      <td>15.767867</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>15.017370</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.49759</td>\n",
       "      <td>0.627978</td>\n",
       "      <td>0.826801</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Gene Symbol    gene_biotype        EPI        3_I       3_II  1_1_1  \\\n",
       "0      TBC1D3L  protein_coding   7.148288   0.000000  36.919539    0.0   \n",
       "4        DNAH5  protein_coding  21.165576  11.080615  12.798294    0.0   \n",
       "7        SAP25  protein_coding   0.000000   7.201708   7.178799    0.0   \n",
       "10        EYA1  protein_coding   0.000000  25.060854   0.000000    0.0   \n",
       "11     CYP39A1  protein_coding   0.000000   0.528252  15.767867    0.0   \n",
       "\n",
       "       1_1_2     1_1_3       2_I       2_II     2_III   fvalue      pvalue  \\\n",
       "0   0.000000  0.000000  0.000000   0.000000  0.000000  2.06819    0.196896   \n",
       "4   0.395609  0.751568  4.215694   2.008978  2.188068  21.1798  0.00107409   \n",
       "7   0.000000  0.000000  8.870301   0.000000  0.000000  1.29337    0.332655   \n",
       "10  0.741551  0.000000  0.000000   0.318718  0.000000  1.18013    0.361688   \n",
       "11  0.000000  0.000000  0.000000  15.017370  0.000000  0.49759    0.627978   \n",
       "\n",
       "   adj_pvalue  \n",
       "0    0.655542  \n",
       "4           1  \n",
       "7    0.732052  \n",
       "10   0.705583  \n",
       "11   0.826801  "
      ]
     },
     "execution_count": 348,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp = pd.DataFrame(index = m11_m12_epi.index, columns = ['fvalue', 'pvalue'])\n",
    "for idx, row in m11_m12_epi.iterrows():\n",
    "    temp.at[idx,'fvalue'], temp.at[idx,'pvalue'] = stats.f_oneway(row[2:5], row[5:8], row[7:11])\n",
    "    \n",
    "m11_m12_epi_anova = pd.concat([m11_m12_epi, temp], axis=1)\n",
    "m11_m12_epi_anova[['adj_pvalue']] = macroH2A_IPS_NAR_2021_workers.fdr(m11_m12_epi_anova['pvalue'])\n",
    "m11_m12_epi_anova.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Print data for Heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "temp4hm = m11_m12_epi_anova[m11_m12_epi_anova['pvalue']<= 0.05]\n",
    "temp4hm = temp4hm.drop(['gene_biotype', 'fvalue', 'pvalue', 'adj_pvalue'], axis=1)\n",
    "temp4hm.head()\n",
    "\n",
    "temp4hm.to_csv(\"m11_m12_epi_4hm.txt\", sep=\"\\t\", header=True, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Compute ANOVA among all conditions (m11, m12, episome)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/t.mazza/miniconda3/envs/playgrounds/lib/python3.7/site-packages/scipy/stats/stats.py:3709: F_onewayConstantInputWarning: Each of the input arrays is constant;the F statistic is not defined or infinite\n",
      "  warnings.warn(F_onewayConstantInputWarning())\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>Gene Symbol</th>\n",
       "      <th>gene_biotype</th>\n",
       "      <th>1_1_1</th>\n",
       "      <th>1_1_2</th>\n",
       "      <th>1_1_3</th>\n",
       "      <th>2_I</th>\n",
       "      <th>2_II</th>\n",
       "      <th>2_III</th>\n",
       "      <th>HU-1</th>\n",
       "      <th>HU-2</th>\n",
       "      <th>HU-3</th>\n",
       "      <th>fvalue</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>adj_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DLL4</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>1.967663</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.888776</td>\n",
       "      <td>3.875864</td>\n",
       "      <td>1.055344</td>\n",
       "      <td>171.317048</td>\n",
       "      <td>159.895481</td>\n",
       "      <td>131.157147</td>\n",
       "      <td>5.91105</td>\n",
       "      <td>0.0313693</td>\n",
       "      <td>0.140032</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>A2M</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>50.713592</td>\n",
       "      <td>42.028808</td>\n",
       "      <td>35.046734</td>\n",
       "      <td>5.77696</td>\n",
       "      <td>0.0329851</td>\n",
       "      <td>0.144026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>FABP4</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>8.925116</td>\n",
       "      <td>5.769943</td>\n",
       "      <td>4.567326</td>\n",
       "      <td>10.841104</td>\n",
       "      <td>24.417390</td>\n",
       "      <td>22.339012</td>\n",
       "      <td>1600.273607</td>\n",
       "      <td>1421.770685</td>\n",
       "      <td>1531.110873</td>\n",
       "      <td>6.34627</td>\n",
       "      <td>0.0267785</td>\n",
       "      <td>0.131111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CXCR4</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>0.991714</td>\n",
       "      <td>8.655214</td>\n",
       "      <td>8.221472</td>\n",
       "      <td>16.755013</td>\n",
       "      <td>102.556588</td>\n",
       "      <td>37.232976</td>\n",
       "      <td>1526.319178</td>\n",
       "      <td>1319.388751</td>\n",
       "      <td>1293.156542</td>\n",
       "      <td>6.24151</td>\n",
       "      <td>0.0278001</td>\n",
       "      <td>0.132937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>IL33</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>1.917317</td>\n",
       "      <td>4.845609</td>\n",
       "      <td>0.883048</td>\n",
       "      <td>0.952738</td>\n",
       "      <td>11.779960</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>422.037827</td>\n",
       "      <td>405.282691</td>\n",
       "      <td>365.291858</td>\n",
       "      <td>6.0656</td>\n",
       "      <td>0.029631</td>\n",
       "      <td>0.136465</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Gene Symbol    gene_biotype     1_1_1     1_1_2     1_1_3        2_I  \\\n",
       "0        DLL4  protein_coding  1.967663  0.000000  0.000000   4.888776   \n",
       "1         A2M  protein_coding  0.000000  0.000000  0.000000   0.000000   \n",
       "2       FABP4  protein_coding  8.925116  5.769943  4.567326  10.841104   \n",
       "3       CXCR4  protein_coding  0.991714  8.655214  8.221472  16.755013   \n",
       "4        IL33  protein_coding  1.917317  4.845609  0.883048   0.952738   \n",
       "\n",
       "         2_II      2_III         HU-1         HU-2         HU-3   fvalue  \\\n",
       "0    3.875864   1.055344   171.317048   159.895481   131.157147  5.91105   \n",
       "1    0.000000   0.000000    50.713592    42.028808    35.046734  5.77696   \n",
       "2   24.417390  22.339012  1600.273607  1421.770685  1531.110873  6.34627   \n",
       "3  102.556588  37.232976  1526.319178  1319.388751  1293.156542  6.24151   \n",
       "4   11.779960   0.000000   422.037827   405.282691   365.291858   6.0656   \n",
       "\n",
       "      pvalue adj_pvalue  \n",
       "0  0.0313693   0.140032  \n",
       "1  0.0329851   0.144026  \n",
       "2  0.0267785   0.131111  \n",
       "3  0.0278001   0.132937  \n",
       "4   0.029631   0.136465  "
      ]
     },
     "execution_count": 349,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp = pd.DataFrame(index = m11_m12_huv.index, columns = ['fvalue', 'pvalue'])\n",
    "for idx, row in m11_m12_huv.iterrows():\n",
    "    temp.at[idx,'fvalue'], temp.at[idx,'pvalue'] = stats.f_oneway(row[2:5], row[5:8], row[7:11])\n",
    "    \n",
    "m11_m12_huv_anova = pd.concat([m11_m12_huv, temp], axis=1)\n",
    "m11_m12_huv_anova[['adj_pvalue']] = macroH2A_IPS_NAR_2021_workers.fdr(m11_m12_huv_anova['pvalue'])\n",
    "del temp\n",
    "m11_m12_huv_anova.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Compute ANOVA among all conditions (m11, m12, HUVEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 366,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "temp4hm = m11_m12_huv_anova[m11_m12_huv_anova['adj_pvalue']<= 0.05]\n",
    "temp4hm = temp4hm.drop(['gene_biotype', 'fvalue', 'pvalue', 'adj_pvalue'], axis=1)\n",
    "temp4hm.head()\n",
    "\n",
    "temp4hm.to_csv(\"m11_m12_huv_4hm.txt\", sep=\"\\t\", header=True, index=False)\n",
    "del temp4hm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Print system and required packages information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPython 3.7.6\n",
      "IPython 7.12.0\n",
      "\n",
      "numpy 1.18.1\n",
      "pandas 1.1.4\n",
      "matplotlib 3.1.3\n",
      "sklearn 0.24.0\n",
      "traitlets 4.3.3\n",
      "IPython 7.12.0\n",
      "ipywidgets 7.5.1\n",
      "openpyxl 3.0.3\n",
      "pyntacle 1.3.2\n",
      "\n",
      "compiler   : MSC v.1916 64 bit (AMD64)\n",
      "system     : Windows\n",
      "release    : 10\n",
      "machine    : AMD64\n",
      "processor  : Intel64 Family 6 Model 60 Stepping 3, GenuineIntel\n",
      "CPU cores  : 8\n",
      "interpreter: 64bit\n",
      " \n",
      "last updated: Mon Mar 08 2021 17:36:33 W. Europe Standard Time\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -v -m -p numpy,pandas,matplotlib,sklearn,traitlets,IPython,ipywidgets,openpyxl,pyntacle\n",
    "\n",
    "# date\n",
    "print(\" \")\n",
    "%watermark -u -n -t -z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}