08_revision_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional analyses for manuscript revisions\n",
    "This notebook contains additional analyses performed for a revised version of the manuscript. In particular, two analyses are performed:\n",
    "1. Determining whether there is a bias in the linear arrangement of motifs in strong enhancers and silencers.\n",
    "2. Associating differentially expressed genes in Crx-/- vs. wildtype P21 retinas with the activity of nearby library members."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import itertools\n",
    "\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable\n",
    "from pybedtools import BedTool\n",
    "from IPython.display import display\n",
    "\n",
    "sys.path.insert(0, \"utils\")\n",
    "import fasta_seq_parse_manip, modeling, plot_utils, predicted_occupancy, sequence_annotation_processing\n",
    "\n",
    "data_dir = \"Data\"\n",
    "downloads_dir = os.path.join(data_dir, \"Downloaded\")\n",
    "figures_dir = \"Figures\"\n",
    "all_seqs = fasta_seq_parse_manip.read_fasta(os.path.join(data_dir, \"library1And2.fasta\"))\n",
    "# Drop scrambled sequences\n",
    "all_seqs = all_seqs[~(all_seqs.index.str.contains(\"scr\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_utils.set_manuscript_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load in MPRA data and other metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mapping activity class to a color\n",
    "color_mapping = {\n",
    "    \"Silencer\": \"#e31a1c\",\n",
    "    \"Inactive\": \"#33a02c\",\n",
    "    \"Weak enhancer\": \"#a6cee3\",\n",
    "    \"Strong enhancer\": \"#1f78b4\",\n",
    "    np.nan: \"grey\"\n",
    "}\n",
    "color_mapping = pd.Series(color_mapping)\n",
    "\n",
    "# Sort order for the four activity bins\n",
    "class_sort_order = [\"Silencer\", \"Inactive\", \"Weak enhancer\", \"Strong enhancer\"]\n",
    "\n",
    "# MPRA measurements\n",
    "activity_df = pd.read_csv(os.path.join(data_dir, \"wildtypeMutantPolylinkerActivityAnnotated.txt\"), sep=\"\\t\", index_col=0)\n",
    "activity_df[\"group_name_WT\"] = sequence_annotation_processing.to_categorical(activity_df[\"group_name_WT\"])\n",
    "activity_df[\"group_name_MUT\"] = sequence_annotation_processing.to_categorical(activity_df[\"group_name_MUT\"])\n",
    "# Only keep sequences that were measured in WT form\n",
    "activity_df = activity_df[activity_df[\"expression_log2_WT\"].notna()]\n",
    "\n",
    "# TF occupancy metrics, also separate out the WT sequences\n",
    "occupancy_df = pd.read_csv(os.path.join(data_dir, \"predictedOccupancies.txt\"), sep=\"\\t\", index_col=0)\n",
    "wt_occupancy_df = occupancy_df[occupancy_df.index.str.contains(\"WT$\")].copy()\n",
    "wt_occupancy_df = sequence_annotation_processing.remove_mutations_from_seq_id(wt_occupancy_df)\n",
    "wt_occupancy_df = wt_occupancy_df.loc[activity_df.index]\n",
    "n_tfs = len(wt_occupancy_df.columns)\n",
    "\n",
    "# PWMs\n",
    "pwms = predicted_occupancy.read_pwm_files(os.path.join(\"Data\", \"Downloaded\", \"Pwm\", \"photoreceptorAndEnrichedMotifs.meme\"))\n",
    "pwms = pwms.rename(lambda x: x.split(\"_\")[0])\n",
    "motif_len = pwms.apply(len)\n",
    "mu = 9\n",
    "ewms = pwms.apply(predicted_occupancy.ewm_from_letter_prob).apply(predicted_occupancy.ewm_to_dict)\n",
    "\n",
    "# WT sequences measured in the assay\n",
    "wt_seqs = all_seqs[all_seqs.index.str.contains(\"WT\")].copy()\n",
    "wt_seqs = sequence_annotation_processing.remove_mutations_from_seq_id(wt_seqs)\n",
    "wt_seqs = wt_seqs[activity_df.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis for linear arrangement bias\n",
    "For each TF besides CRX, identify strong enhancers with exactly one position occupied by CRX and exactly one position occupied by the other TF. Count the number of times the occupied position is 5' (left) or 3' (right) of the CRX position. Because all sequences are centered on CRX motifs in a forward orientation, we do not need to consider orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>left</th>\n",
       "      <th>right</th>\n",
       "      <th>binom_pval</th>\n",
       "      <th>binom_qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>GFI1</th>\n",
       "      <td>8</td>\n",
       "      <td>10</td>\n",
       "      <td>0.814529</td>\n",
       "      <td>0.814529</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MAZ</th>\n",
       "      <td>13</td>\n",
       "      <td>25</td>\n",
       "      <td>0.072951</td>\n",
       "      <td>0.255330</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MEF2D</th>\n",
       "      <td>7</td>\n",
       "      <td>3</td>\n",
       "      <td>0.343750</td>\n",
       "      <td>0.802083</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NDF1</th>\n",
       "      <td>17</td>\n",
       "      <td>22</td>\n",
       "      <td>0.522397</td>\n",
       "      <td>0.609464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NRL</th>\n",
       "      <td>14</td>\n",
       "      <td>31</td>\n",
       "      <td>0.016094</td>\n",
       "      <td>0.112661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RORB</th>\n",
       "      <td>14</td>\n",
       "      <td>20</td>\n",
       "      <td>0.391528</td>\n",
       "      <td>0.548140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RAX</th>\n",
       "      <td>8</td>\n",
       "      <td>4</td>\n",
       "      <td>0.387695</td>\n",
       "      <td>0.678467</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       left  right  binom_pval  binom_qval\n",
       "GFI1      8     10    0.814529    0.814529\n",
       "MAZ      13     25    0.072951    0.255330\n",
       "MEF2D     7      3    0.343750    0.802083\n",
       "NDF1     17     22    0.522397    0.609464\n",
       "NRL      14     31    0.016094    0.112661\n",
       "RORB     14     20    0.391528    0.548140\n",
       "RAX       8      4    0.387695    0.678467"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "occupied_cutoff = 0.5\n",
    "strong_mask = activity_df[\"group_name_WT\"].str.contains(\"Strong\")\n",
    "strong_mask = strong_mask & strong_mask.notna()\n",
    "# {tf name: [number of times TF is 5' of central CRX, number of times 3']}\n",
    "tf_order_counts = {}\n",
    "\n",
    "for tf in ewms.index.drop(\"CRX\"):\n",
    "    left_counts = right_counts = 0\n",
    "    # Get strong enhancers with a motif for this TF\n",
    "    has_tf_mask = (wt_occupancy_df[tf] > occupied_cutoff) & strong_mask\n",
    "    has_tf_seqs = wt_seqs[has_tf_mask]\n",
    "    # Get the predicted occupancy landscape only for CRX and the other TF\n",
    "    for seq in has_tf_seqs:\n",
    "        occupancy_landscape = predicted_occupancy.total_landscape(seq, ewms[[\"CRX\", tf]], mu) > occupied_cutoff\n",
    "        # Only consider the sequence if there is exactly one CRX and exactly one of the other TF\n",
    "        fwd_counts = occupancy_landscape[tf + \"_F\"].sum()\n",
    "        rev_counts = occupancy_landscape[tf + \"_R\"].sum()\n",
    "        if (occupancy_landscape[\"CRX_F\"].sum() == 1) and (occupancy_landscape[\"CRX_R\"].sum() == 0) and (fwd_counts + rev_counts == 1):\n",
    "            # By construction, the motif will only be in the F or R column\n",
    "            if fwd_counts == 1:\n",
    "                col = occupancy_landscape[tf + \"_F\"]\n",
    "            else:\n",
    "                col = occupancy_landscape[tf + \"_R\"]\n",
    "            \n",
    "            tf_start_pos = col[col].index[0]\n",
    "            # CRX start position should always be the same, but just in case it's not do this\n",
    "            crx_occ = occupancy_landscape[\"CRX_F\"]\n",
    "            crx_start_pos = crx_occ[crx_occ].index[0]\n",
    "            if tf_start_pos < crx_start_pos:\n",
    "                left_counts += 1\n",
    "            elif tf_start_pos > crx_start_pos:\n",
    "                right_counts += 1\n",
    "            # else they start at the same position, ignore it\n",
    "            \n",
    "    tf_order_counts[tf] = [left_counts, right_counts]\n",
    "    \n",
    "tf_order_counts = pd.DataFrame.from_dict(tf_order_counts, orient=\"index\", columns=[\"left\", \"right\"])\n",
    "tf_order_counts[\"binom_pval\"] = tf_order_counts.apply(stats.binom_test, axis=1)\n",
    "tf_order_counts[\"binom_qval\"] = modeling.fdr(tf_order_counts[\"binom_pval\"])\n",
    "display(tf_order_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now do the same analyses as above for silencers with NRL. Then create a contingency table to determine whether the left/right positioning of NRL motifs is independent of whether a sequence is a strong enhancer or silencer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>left</th>\n",
       "      <th>right</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Silencer</th>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Strong enhancer</th>\n",
       "      <td>14</td>\n",
       "      <td>31</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 left  right\n",
       "Silencer            6      5\n",
       "Strong enhancer    14     31"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The linear arrangement of CRX-NRL motifs is independent of whether a sequence is a strong enhancer or silencer, Fisher's exact test p=0.17, odds ratio=2.7\n"
     ]
    }
   ],
   "source": [
    "silencer_mask = activity_df[\"group_name_WT\"].str.contains(\"Silencer\")\n",
    "silencer_mask = silencer_mask & silencer_mask.notna()\n",
    "has_nrl_mask = (wt_occupancy_df[\"NRL\"] > occupied_cutoff) & silencer_mask\n",
    "has_nrl_seqs = wt_seqs[has_nrl_mask]\n",
    "\n",
    "silencer_left_counts = silencer_right_counts = 0\n",
    "for seq in has_nrl_seqs:\n",
    "    occupancy_landscape = predicted_occupancy.total_landscape(seq, ewms[[\"CRX\", \"NRL\"]], mu) > occupied_cutoff\n",
    "    # Only consider the sequence if there is exactly one CRX and exactly one of the other TF\n",
    "    fwd_counts = occupancy_landscape[\"NRL_F\"].sum()\n",
    "    rev_counts = occupancy_landscape[\"NRL_R\"].sum()\n",
    "    if (occupancy_landscape[\"CRX_F\"].sum() == 1) & (occupancy_landscape[\"CRX_R\"].sum() == 0) & (fwd_counts + rev_counts == 1):\n",
    "        # Samw as above\n",
    "        if fwd_counts == 1:\n",
    "            col = occupancy_landscape[\"NRL_F\"]\n",
    "        else:\n",
    "            col = occupancy_landscape[\"NRL_R\"]\n",
    "\n",
    "        tf_start_pos = col[col].index[0]\n",
    "        crx_occ = occupancy_landscape[\"CRX_F\"]\n",
    "        crx_start_pos = crx_occ[crx_occ].index[0]\n",
    "        if tf_start_pos < crx_start_pos:\n",
    "            silencer_left_counts += 1\n",
    "        elif tf_start_pos > crx_start_pos:\n",
    "            silencer_right_counts += 1\n",
    "\n",
    "silencer_order_counts = pd.Series([silencer_left_counts, silencer_right_counts], index=[\"left\", \"right\"])\n",
    "# Join with strong enhancer counts\n",
    "contingency_table = pd.DataFrame.from_dict({\"Silencer\": silencer_order_counts, \"Strong enhancer\": tf_order_counts.loc[\"NRL\", [\"left\", \"right\"]]}, orient=\"index\").astype(int)\n",
    "display(contingency_table)\n",
    "odds, pval = stats.fisher_exact(contingency_table)\n",
    "print(f\"The linear arrangement of CRX-NRL motifs is independent of whether a sequence is a strong enhancer or silencer, Fisher's exact test p={pval:.2f}, odds ratio={odds:.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of gene expression changes in Crx-/- retina\n",
    "First, read in the RNA-seq data, CPM normalize each replicate, compute mean expression across replicates, and then compute the fold changes between Crx-/- and WT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rnaseq_df = pd.read_csv(os.path.join(downloads_dir, \"rogerRnaseqRaw2014.txt\"), sep=\"\\t\", usecols=[\"Wta\", \"Wtb\", \"Crxa\", \"Crxb\"])\n",
    "# Add a pseudocount too when CPM normalizing\n",
    "rnaseq_df = (rnaseq_df + 1) / (rnaseq_df.sum() / 1e6)\n",
    "# Compute averages and then log2 FC\n",
    "ko_wt_fc_log2 = np.log2(rnaseq_df[[\"Crxa\", \"Crxb\"]].mean(axis=1) / rnaseq_df[[\"Wta\", \"Wtb\"]].mean(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, read in a BED file of the library and intersect it with [Supplementary file 4](https://doi.org/10.7554/eLife.48216.022) from Murphy *et al.*, 2019 to associate sequences to genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bed_columns = [\"chrom\", \"begin\", \"end\", \"label\", \"score\", \"strand\"]\n",
    "# Load in BED file\n",
    "library_bed = pd.read_csv(os.path.join(data_dir, \"library1And2.bed\"), sep=\"\\t\", header=None, names=bed_columns)\n",
    "# Pull out sequences that were measured\n",
    "library_bed = library_bed.set_index(\"label\").loc[activity_df.index].reset_index()[bed_columns]\n",
    "library_bed = BedTool.from_dataframe(library_bed).sort()\n",
    "\n",
    "# Read in ATAC data and intersect with the library\n",
    "atac_df = pd.read_excel(os.path.join(downloads_dir, \"murphyAtac2019.xlsx\"), sheet_name=\"peakUnion_counts\", skiprows=1)\n",
    "atac_bed = BedTool.from_dataframe(atac_df[[\"Chr\", \"Start\", \"End\", \"PeakID\"]]).sort()\n",
    "atac_bed = atac_bed.intersect(library_bed, wo=True).to_dataframe()\n",
    "atac_bed = atac_bed[[\"chrom\", \"start\", \"end\", \"name\", \"thickEnd\"]].rename(columns={\"thickEnd\": \"library_label\"})\n",
    "atac_df = atac_df.set_index(\"PeakID\").loc[atac_bed[\"name\"]].reset_index()\n",
    "atac_df[\"library_label\"] = atac_bed[\"library_label\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the gene closest to every library member, and then get the fold change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nearest_gene_fc(gene):\n",
    "    if gene in ko_wt_fc_log2.index:\n",
    "        return ko_wt_fc_log2[gene]\n",
    "    else:\n",
    "        return np.nan\n",
    "    \n",
    "library_nearest_gene = atac_df.set_index(\"library_label\")[\"Nearest Gene\"]\n",
    "library_gene_fc = library_nearest_gene.apply(get_nearest_gene_fc)\n",
    "library_gene_fc = library_gene_fc[library_gene_fc.notna()]\n",
    "activity_has_gene_df = activity_df.loc[library_gene_fc.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now get genes that have an absolute fold change of 2 (log2 = 1) or more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "library_gene_de = library_gene_fc[library_gene_fc.abs() >= 1]\n",
    "up_mask = (library_gene_de > 0).replace({False: \"Down-regulated\", True: \"Up-regulated\"})\n",
    "de_direction_grouper = library_gene_de.groupby(up_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine enrichment of silencers being near up-regulated genes and enhancers being near down-regulated genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direction of differential expression is independent of having a silencer nearby, Fisher's exact test p=0.001, OR=2.13\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>False</th>\n",
       "      <th>True</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nearest Gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Down-regulated</th>\n",
       "      <td>0.866279</td>\n",
       "      <td>0.133721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Up-regulated</th>\n",
       "      <td>0.752427</td>\n",
       "      <td>0.247573</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   False     True \n",
       "Nearest Gene                      \n",
       "Down-regulated  0.866279  0.133721\n",
       "Up-regulated    0.752427  0.247573"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direction of differential expression is independent of having an enhancer nearby, Fisher's exact test p=0.022, OR=1.52\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>True</th>\n",
       "      <th>False</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nearest Gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Down-regulated</th>\n",
       "      <td>0.680233</td>\n",
       "      <td>0.319767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Up-regulated</th>\n",
       "      <td>0.582524</td>\n",
       "      <td>0.417476</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   True      False\n",
       "Nearest Gene                      \n",
       "Down-regulated  0.680233  0.319767\n",
       "Up-regulated    0.582524  0.417476"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Fraction of genes near an enhancer')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAD/CAYAAAC6q7DNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHodJREFUeJzt3XmYXFW57/HvDwhTBiGQEwbtJEgCxyhBjSOG4YoKBxAURR5yQVSMwI3KkcMBBwRynfC5+jiAQFARmRWCMiigBwkgKoIQJB6MDAmTHKaQpEMgJLz3j7WK7BTV3auH6q5Qv8/z1NNVa6/a+62uzpu111p7bUUEZmZDbb2hDsDMDJyMzKxFOBmZWUtwMjKzluBkZGYtwcnIzFqCk5GZtQQnIzNrCU5GZtYSNhjqAAbLlltuGePHjx/qMMzazu233/5kRIzpqV7bJKPx48dz2223DXUYZm1H0qKSej5NM7OW4GRkZi3BycjMWsKgJSNJoyVdLmm5pEWSDumi3nGS7pa0TNIDko6r275Q0gpJnflx3eB8AjNrpsHswD4dWAmMBXYGrpY0LyLm19UTcBhwF/Ba4DpJD0XExZU6+0XEbwcjaDMbHIPSMpI0HDgQODEiOiPiZuAK4ND6uhHxzYj4S0Ssioi/A78EdhmMOM1s6AzWadokYFVELKiUzQMmd/cmSQKmAfWtpwskPSHpOklTBjZUMxsKg5WMRgBL68qWACN7eN/JpBjPqZRNB8YD44DfAddK2qzRmyXNkHSbpNueeOKJPoRtZoOlKBlJ6m/S6gRG1ZWNApZ1c8yZpL6jfSLi+Vp5RPw+IlZExLMR8XXgGVLr6WUiYnZETI2IqWPG9DgB1MyGUI8d2JLWBzolbVZNCr20ANhA0sSI+Ecum8LLT79qx/w4cAKwa0Q83MO+g9Tpbcb4E64e6hDaxsJv7DOg++uxxRMRq0nJZIu+HiQilgNzgFmShkvaBdgfOK++rqTpwNeA90TE/XXbOiTtImlDSRvnYf8tgd/3NTYzaw2lQ/sXAFdJ+i7wMKk1AkBEXF+4j6OBHwOPA08BR0XEfEnTgF9HxIhc7yukxPfn1H8NwPkRcSSpj+kM0pD/c8CdwN4R8VRhDGbWokqT0VH558l15QFsV7KDiHgaOKBB+U2kDu7a6wnd7GM+sFPJ8cxs3VKUjLpLEGZmA6F4lEzSMEnTJH0kvx6eJzOamfVb6dD+G0id2GcDP8rFu5H6gMzM+q20ZXQG8OWI2BF4IZfNBd7VlKjMrO2UJqPJwPn5ecBLw/WbNCMoM2s/pcloIfDmaoGktwL3DnRAZtaeSof2TyQt+XEmsKGkzwNHAp9sWmRm1laKWkYRcRWwFzCG1Fc0DvhgRHhhMzMbEMWLq0XEHaRZ1GZmA650aH9OvmyjWjZN0qXNCcvM2k1pB/ZuwC11ZX8A9hjYcMysXZUmo+eA+tnWI1gz58jMrF9Kk9G1wFmSRgHkn6cB1zQrMDNrL6XJ6FjSyoxPS3oceBp4FXBMswIzs/ZSetX+YmAfSVsBrwEeiojHmhqZmbWV3q5t/SJpYbRNJW0nqWgtIzOznhS1jCTtRbpaf+u6TQGsP9BBmVn7KW0ZnQ78X2B4RKxXeTgRmdmAKJ2BvTlwVkREjzXNzPqgtGX0I+BjzQzEzNpbacvo7cBnJJ0ArDWKFhG7DnhUZtZ2SpPRD/PDzKwpSucZndvsQMysvZVetS9Jn5R0vaS7ctmukg5qbnhm1i5KO7BnAZ8AZgMduexh4PhmBGVm7ac0GR0O7BsRF7Pm1tYPUHg3WTOznpQmo/WBzvy8loxGVMrMzPqlNBn9Cvi2pI0g9SGRZmRf2azAzKy9lCajz5GuS1tCWjqkk7Qov/uMzGxAlA7tLwU+IGksqQPbS4iY2YDqMhlJatRqeiI/XtoeES82JzQzayfdtYxWsaazuhHhJUTMbIB012c0gTR039Wjtr2IpNGSLpe0XNIiSYd0Ue84SXdLWibpAUnH1W0fL+l3kp6VdI+kPUtjMLPW1WXLKCIWDfCxTgdWAmOBnUm3y54XEfPr6gk4DLgLeC1wnaSH8hwngItIt0n6t/y4VNLEiHhigOM1s0HUXZ/R7IiYkZ+fRxenbBFxWE8HkTQcOBB4fUR0AjdLugI4FDihbn/frLz8u6RfArsAF0uaBLwJeG9ErAAuk3RM3veZPcVhZq2ruz6jByrP7+3ncSYBqyJiQaVsHunmkF3K85mmAWflosnA/RGxrG4/k/sZn5kNse5O075eeX5KP48zAlhaV7YEGNnD+04m9WudU9nPkgb72bbRmyXNAGYAdHR0NKpiZi2i9Kr9PSRNyM+3knSupHPyrYtKdJLuu1Y1CljWoG7tmDNJfUf7RMTzfdlPRMyOiKkRMXXMmDGFoZrZUCidgf0DYHV+/m1gGOm2RbML378A2EDSxErZFKC+8xoASR8n9SW9OyIermyaD2wnqdqi6nI/ZrbuKF3pcduIeFDSBsD7SJeCrAQeLXlzRCyXNAeYJekI0mja/sA76+tKmg58DdgjIu6v288CSXcCJ0n6ErA3sBOpA9vM1mGlLaOl+VKQ3YC/5RExSC2kUkcDmwCPk4bnj4qI+ZKmSape/f8VYAvgz5I686M6UnYwMBVYDHwD+JCH9c3WfaUto+8DfwY2BI7JZbsA95QeKCKeBg5oUH4TqWO69npCD/tZCOxeelwzWzeUXih7qqTLgdURcV8ufgQ4ommRmVlbKW0ZUTdH6GWvzcz6o7TPyMysqZyMzKwlOBmZWUso7jOqydeLqfbai6uZ2UAovRxkG0lzJD1FWnTthcrDzKzfSk/TziIlnneTrg97E3AFcGST4jKzNlN6mvZOoCNf1hERMU/SJ4BbgLObF56ZtYvSltFq0ukZwDOSxgDL6WLpDjOz3ipNRn8iLfEKcC1wCTAHuK0ZQZlZ+yk9TTuUNYnrGOBY0sJo32lGUGbWfkqvTXum8nwF6cp6M7MB40mPZtYSnIzMrCU4GZlZS+gxGUlaX9INkjYajIDMrD31mIwiYjXpVtZuRZlZ05QmmFOAMySNyy2l9WqPZgZnZu2jdJ7RD/PPQytlIt3yev0BjcjM2lJpMup2kXwzs/4qnfS4qNmBmFl7K15cTdL7SfdN25K1F1c7rAlxmVmbKV1c7STSmkbrAR8GniLdWfaZ7t5nZlaqdDTs48B7IuLfgZX5537A+GYFZmbtpTQZbRYRd+fnKyUNi4hbSadtZmb9VtpndJ+kyRExH7gbOErSYtL97s3M+q00GX0J2CI//zxwATACOLoZQZlZ+ykd2v9V5fmfgO2bFpGZtaXeDO3vSBpJGxsRMyXtAGwUEXc1LTozaxulQ/sfBm4kLcBfm1c0Evh2k+IyszZTOpo2izS0fyTpTiEA84ApTYnKzNpOaTL6F6B2OhaVn9G4+stJGi3pcknLJS2SdEgX9faQ9DtJSyQtbLB9oaQVkjrz47rSGMysdZUmo9tZ+4p9gIOBW3txrNOBlcBYYDppSZLJDeotB34MHNfNvvaLiBH58d5exGBmLaq0A/szwHX5LrLDJV0LTAKKEoGk4cCBwOsjohO4WdIVpAR3QrVunkx5q6Q9C2Mzs1eA0qH9e/Jo2r7AVcBDwFU5sZSYBKyKiAWVsnn0fQb3BXlhtzuA4yJiXqNKkmYAMwA6Ojr6eCgzGwzFQ/sR8Szwsz4eZwSwtK5sCWlErremA38hrRzwWeBaSTtW7+1WExGzgdkAU6dOLe7fMrPBVzq0P0HShZL+JunB6qPwOJ3AqLqyUcCy3gQLEBG/j4gVEfFsRHydtHLAtN7ux8xaS2nL6ELgPtJtrZ/tw3EWABtImhgR/8hlU4D5fdhXvaCyvpKZrZtKk9FkYJeIeLEvB4mI5ZLmALMkHQHsDOwPvLO+bu4L2hAYll5qY+DFiFgpqQN4DfBnUqvu06TF3n7fl7jMrHWUDu3fCLyxn8c6GtgEeBy4CDgqIuZLmiap2hG+K7AC+BXQkZ/X5hKNBM4grRbwCLAXsHdEPNXP2MxsiJW2jBYC10i6HHisuiEivlyyg4h4GjigQflNpA7u2usb6OK0Ky9hslNhzGa2DilNRsNJQ/rDSKdJNR6hMrMBUTrP6GPNDsTM2pvvCGtmLcHJyMxagpORmbUEJyMzawmll4PsIWlCfr61pHMlnSNpq+aGZ2btonRo/wekO8gCfCv/XEG6CPX9Ax3UUBp/wtVDHUJbWPiNfYY6BGsxpclo24h4UNIGpKQ0jrRQ2qNNi8zM2kppMloqaSzweuBvEdEpqXb9mJlZv5Umo++TLk7dEDgml+0C3NOMoMys/ZTOwD41X5e2OiLuy8WPAEc0LTIzayu9Gdp/ANhG0kfy60eA+wc+JDNrR6VD+28gLZB2NvCjXLwb6S4eZmb9VtoyOgP4ckTsCLyQy+YC72pKVGbWdkqT0WTg/Pw8IK3eSFoszcys30qT0ULgzdUCSW8F7h3ogMysPZUO7Z8IXC3pTGBDSZ8HjgQ+2bTIzKytFLWMIuIq0nrTY0h9ReOAD0aE73NvZgOiNzdxvIO0qL6Z2YArSkb50o/DSbcYGlHdFhGHDXxYZtZuSltG55Juungl8D/NC8fM2lVpMtoLmNDofvZmZgOhdGj/QWCjZgZiZu2ttGX0U+CXkr5L3WlaRFw/4FGZWdspTUYz88+v1ZUHsN3AhWNm7ap0CZEJzQ7EzNpb8RIikoZJmlZbQkTScEnDmxeambUTLyFiZi3BS4iYWUvwEiJm1hIGbQkRSaMlXS5puaRFkg7pot4ekn4naYmkhQ22j8/bn5V0j6Q9S2Mws9ZVmoxqS4icwpolRH4OfKkXxzqddK+1scB04AxJkxvUW07qizqui/1cBNwBbAF8EbhU0phexGFmLWhQlhDJo24HAidGRGdE3AxcARza4Fi3RsR5NFjsX9Ik4E3ASRGxIiIuA/6a921m67DBWkJkErAqIhZUyuaRRuR6YzJwf0Qsq9tPoxaWma1DSpcQmdXFpueBh4FrIqK7q/lHAEvrypYAI0uOX7efJQ32s22jypJmADMAOjo6enkoMxtMpX1Gk4DjgT2A7fPP44E3AkcB90vaq5v3dwKj6spGAcsa1O1Or/YTEbMjYmpETB0zxt1KZq2sNBmtBxwcEdMi4pCImAYcRLrD7NtJp2/f6Ob9C4ANJE2slE0B5vcy3vnAdpKqLaq+7MfMWkxpMnofqcO56ipg7/z8fLq5YDbPSZoDzMqXkewC7A+cV19X0nqSNgaGpZfaOK80Se5zuhM4KZd/ANgJuKzwc5hZiypNRveRTseqjszlAFsCz/awj6NJkyQfJw3PHxUR8/P1bp2VersCK4BfAR35eXXU7mBgKrCY1Br7UEQ8Ufg5zKxFlY6mHQHMkXQ88Aipw3g18MG8fQfSXKQuRcTTwAENym+isq52RNwAqJv9LAR2L4zbzNYRpUuI/CX397wd2Ab4J/CHiHghb78RuLFpUZrZK15v5hm9ANzUxFjMrI0Vr2dkZtZMTkZm1hK6TEaS3l95PmxwwjGzdtVdy+j8yvOnmh2ImbW37jqwH5M0E/gbafb0HjQYcvetisxsIHSXjA4HZgGfBTak8XrXvlWRmQ2ILpNRRNwC7Akg6d6I2H7QojKztlM66XF7AEkdpNnXD0fEQ80MzMzaS+mtiraSNJe05vUc4D5JN0rapqnRmVnbKJ1ndCZpRcXNI2JrYHPSOtRnNiswM2svpZeDvAvYunIt2nJJ/0m6aNbMrN9KW0aLgdfVle0APDOw4ZhZuyptGX0T+K2kHwGLSHcH+Rg9LBtiZlaqdDTtbEn3AYeQVlZ8FDgkIv6rmcGZWfvozRIi1wOebW1mTeGr9s2sJTgZmVlLcDIys5bgZGRmLaH09tajgf8AdqZyJw+AiNi1CXGZWZspHU27ENgI+Bk93x/NzKzXSpPRO4ExEfF8M4Mxs/ZV2md0F/DqZgZiZu2ttGV0PXCNpHOAx6obIqLRCpBmZr1SmoymAQ8D76krDxovR2tm1iul16bt0exAzKy9FV+bJmlzYD/SsrOPAFdGxOJmBWZm7aV02dl3APcBR5Ku2v8UaenZdzQxNjNrI6Uto+8AR0fExbUCSR8Bvge8pRmBmVl7KR3an0Sa8Fh1KeDbF5nZgChNRv8ADq4r+zDp1K2IpNGSLpe0XNIiSYd0UU+STpX0VH6cKkmV7ZH30ZkfPyyNwcxaV+lp2jHAVZI+Q1p2djwwEdi3F8c6HVgJjCVd43a1pHkRMb+u3gzgAGAKaerAb4AHWPtOJFMi4t5eHNvMWlxRyyjfXfa1wGnA7cD3ge1zeY8kDQcOBE6MiM6IuBm4Aji0QfWPAt+KiIcj4hHgW6RbbZvZK1hvlp1dDJzfx+NMAlZFxIJK2TxgtwZ1J+dt1XqT6+rcKGk94BbgcxGxsI9xmVmL6DIZSbomIvbKz28inTK9TOESIiOApXVlS4CRXdRdUldvhCRFRJAS2B+BTYGvkE4fd46IVQ0+wwzSaR8dHR0FYZrZUOmuZfTTyvP+dhJ3AqPqykYBywrqjgI6cyIiIm7M5SslfZaU5P4V+Gv9jiJiNjAbYOrUqQ2TqZm1hi6TUURcWHl5T0T8qb6OpLcWHmcBsIGkiRHxj1w2BajvvCaXTQFu7aHeS6EC6ma7ma0DSof2f9NF+TUlb46I5cAcYJak4ZJ2AfYHzmtQ/afA5yRtK2kb4FjgJwCSJkvaWdL6kkaQOrcfAf678HOYWYvqNhlJWk/S+rw0/UfrVR4TgZf103TjaGAT4HHgIuCoiJgvaZqkzkq9s4ArSadddwNX5zJI0wIuIZ2a3U+aYrBvRLzQizjMrAX1NJq2ijUd1/WJ50Xgq6UHioinSfOH6stvorKudu4b+s/8qK97PbBD6THNbN3RUzKaQOqPmQtUR80CeCIiVjQrMDNrL90mo4hYBCBpB2B19XRI0jBJG3ldbDMbCKUd2NcBb64rezNw7cCGY2btqjQZ7QTUD+3fShp2NzPrt9Jk9AxpJKtqLLB8YMMxs3ZVmowuAy6U9HpJm0p6A2k+UP0aR2ZmfVKajL5Imlh4K+kSjj8Cfwe+0KS4zKzNlN4d5Dng/0iaCWwJPFm7VszMbCAULyGSjciPkbXFFyPi/oEOyszaT1EykvQ64ALWrL4o1szMXr85oZlZOyntM/oB8DtgNOm6sM1J14t9tElxmVmbKT1NmwK8JyJeyIucLZF0HOlC1r6u/mhm9pLSltFzwLD8/ElJHfm9WzQlKjNrO6XJ6CbgoPz8UuDXpItnr29GUGbWfkqH9g+qvPwC6fRsJGsvTWtm1mc9JqO8uNp/Ae+LiOcj4kXcT2RmA6zH07SIWE1a16j0lM7MrNdKE8wpwBmSxuX1p19afraZwZlZ+ygd2q/dqqh6B9jaxEdPejSzfitNRhOaGoWZtb1uk5GkrSLisdrys2ZmzdJTn8+C6gtJc5oYi5m1sZ6SUf2dWndvUhxm1uZ6SkZes8jMBkVPHdgbSNqDNS2k+te1GyuamfVLT8noceDHlddP1b0OYLuBDsrM2k9PN3EcP0hxmFmb8wxqM2sJTkZm1hKcjMysJTgZmVlLcDIys5YwaMlI0mhJl0taLmmRpEO6qCdJp0p6Kj9OVe0mbWn7zpJul/Rs/rnzYH0GM2uewWwZnQ6sBMYC00nrI01uUG8GcADpjiQ7AfsBnwKQtCHwS9JKk5sD5wK/zOVmtg4blGQkaThwIHBiRHRGxM3AFay9PlLNR4FvRcTDEfEI8C3g8Lxtd9LcqO/kJXC/R5oN/r+a/BHMrMkGq2U0CVgVEdVVAOYBjVpGk/O2RvUmA3dFRPWaubu62I+ZrUNKF1frrxGkO9FWLSHdYaRR3SV19UbkfqP6bd3tB0kzSKd9AJ2S/t7LuNcVWwJPDnUQvaFThzqClvNK/g7HlVQarGTUCYyqKxsFLCuoOwrojIiQ1Jv9EBGzgdl9ingdIum2iJg61HFY3/k7HLzTtAWkK/4nVsqmAPMb1J2ftzWqNx/YqTq6RurkbrQfM1uHDEoyiojlwBxglqThknYB9gfOa1D9p8DnJG0raRvgWOAnedsNwGrgM5I2kjQzl3sZE7N13GAO7R8NbEJaluQi4KiImC9pWj79qjkLuBL4K+nOtVfnMiJiJWnY/zDgGeDjwAG5vJ294k9F20Dbf4dae2DKzGxo+HIQM2sJTkb2EkkLJe05SMf6iaSvDMaxXikG83cm6XBJNw/GsWqcjLL8D3GFpGWSnpF0i6QjfQvvxiSFpO2HOo7eaBSzpJMlnT9UMTWLpBskHTHUcfSG/6Gtbb+IGEmapPUN4HjgR0Mb0hqSBmtemPWBv5/+cTJqICKWRMQVwEeAj0p6vaRXSfqppCfyqgNfqrWa8us35+fT8//Ak/PrT0j6RX5+sqSf5f0skzRfUpcT3XL9SyWdL2kpcLik9SSdIOm+vKrBzySNrrznsBzPU5JOrJ561TfzJe0u6eEujv1WSX/IrcR/SjqtdkGypBtztXmSOiV9JJfvK+nOSstyp8r+3ijpL/lzXwJs3Osvpslqvw9JX5D0ZP7dTS+of7ykx4Bzcnl3v4c3Sboj/x5+LumS2nfS6NSoqxaopM0lXZX/Hhfn56/O274KTANOy9/Pabl8R0m/kfS0pL9LOqiyvy0kXSFpqaRbgdf253fZF05G3YiIW4GHSV/s94FXke6GshtpesHHctW5rLnB5W7A/cCulddzK7t9P3AxsBnpYuHTeghjf+DSXP8C4NOk6Q27AdsAi0krIiDpdcAPSKsibJ3j3bY3n7liNfDvpMsU3gG8mzQ9g4iofbYpETEiIi6R9EbSnWM+BWxBmo5xhdJ8sA2BX5DmlY0Gfk66cLoVbUX6zNuSLtqeLWmHHuqPJrWmZxT8Hi4nzZsbTZri8oE+xrkeKfmNAzqAFeS/pYj4InATMDN/PzOVLlb/DXAh8C/AwcAP8t8MpL+h50h/Nx/Pj0HlZNSzR0l/OAcDn4+IZRGxkLSaQG3Vgbmk5AApcX298ro+Gd0cEb+KiNWkf5zV2eaN/CEifhERL0bECuBI4It5VYPngZOBD+VThA8BV0bEzXnu1Zfp4404I+L2iPhjRKzKn/esymdqZAZwVkT8KSJWR8S5wPPA2/NjGGm1hRci4lLgz32Ja5CcmFeFmEua53ZQN3VfBE7K9VfQ8+9hA+B7+fcwB7i1LwFGxFMRcVlEPBsRy4Cv0v33sy+wMCLOyd/pHcBlwIclrU/6z+HLEbE8Iu4mLc8zqJyMerYt6Q9oGLCoUr6INa2OucA0SVsD6wM/A3aRNJ7UOrmz8r7HKs+fBTaWtEE+vevMj19X6jxUF8844PJ8CvAM8N+kVsxYUkvppfoR8SzpXne9JmlSbvo/lk8Rv0ZqMXRlHHBsLa4c22tyTNsAj9SttrCo0U6abDXpe6waBrxQeb04XzFQswjYRlJH5fupTtJ9IiKeq7zu7e+h/vstImlTSWflU/KlwI3AZjmxNDIOeFtdXNNJLbsxpL/xaiyD/v04GXVD0ltICecXpD/Y6tXHHcAjABFxLymxfBq4MSKWkpLODFJL6MWejhURF+Qm9YiI2Lu6qa7qQ8DeEbFZ5bFxXvvpn8CrK/FvQjpVqFkObFp5vVU3IZ0B3ANMjIhRwBeo3Em4gYeAr9bFtWlEXJTj2lZa65rCjm721SwPAuPryiaw9j+8zfMpTU0H8GhEPFj5fkZUtjf6fnrze3hN5fla34+k7r6fY4EdgLfl76d26lzbd6O45tbFNSIijgKeAFbVxTLo34+TUQOSRknal9S3c35EzCO1dr4qaaSkccDnSCtO1swFZrLmlOyGutcD5cwcx7gc6xhJ++dtlwL7SXpn7p84mbUTyJ3AvyktAbwVcEw3xxlJWvalU9KOwFF12/+Hte8mfDZwpKS3KRkuaR9JI4E/kP7YPyNpmKQPAm/tw2fvr0uAL0l6tdJAwJ6klUQvrat3iqQNJU0jnd78vBfH6On3sBqYmVvD+7P272EeMFlpaeWNSd9fV0aS+omeURrAOKlue/33cxUwSdKh+TsYJuktkv41dxnMAU7OLa7XkfrLBpWT0dqulLSM9L/IF4Fvs6aT+tOk/7nuB24mdQRWb/U9l/QHcmMXrwfKd0kd39flWP8IvA0gIubnOC8m/S/cSboW8Pn83vNIf/ALgetI/zi78h/AIaTlWc5uUPdk4Nzc5D8oIm4DPknqRF0M3EteoTP3X30wv36aNEo5p/cfvd9mAbeQvr/FwDeB6bmPpOaxvO1R0oDBkRFxT+kBCn8PnyBdW/m/SUni+bx9QY7xt8A/cpxd+Q7pWs8nSX8D19Rt/y6pL3GxpO/lfqX3kvo+H82f81Rgo1x/Jmm9sMdIHeznlH7mgeJr017BJI0g/dFPjIgHhjqeVidpd1JL+NU91R3AY/4JODMiBv0ff6txy+gVRtJ+uak9HPh/pNUPFg5tVFYjaTdJW+XTtI+S1uOqb9W0JSejV579Sc3wR4GJwMHh5m8r2YF0qvwMqRP6QxHxz6ENqTX4NM3MWoJbRmbWEpyMzKwlOBmZWUtwMjKzluBkZGYtwcnIzFrC/wcpM4N6QY0LCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARsAAAD/CAYAAADBj3IBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHwdJREFUeJzt3Xuc1VW9//HXWwEzLilJmNgAJehPTLqQmUbqo+xy0vBkXn6YZmUc4Wedrie1NOORnfR36tdFU+limpejGRhKXjrH5PKwND2FSRmmQSoSKNdBE5HP74/1nfjOZs/sNczsPcOe9/PxmMfsvfb67vWZ2cOH9V3f9V1LEYGZWb3t0tsBmFn/4GRjZg3hZGNmDeFkY2YN4WRjZg3hZGNmDeFkY2YN4WRjZg3hZGNmDTGgtwPoKXvttVeMGTOmt8Mw63ceeOCBpyNiRK16DUs2koYDPwDeCTwNnBMR11WpdxswuVQ0CPhTRLy2s/cfM2YM999/fw9GbGY5JC3PqVcz2UjaBTgSWBQRm7sR06XAZmAk8DpgnqTFEbGkXCki3lPR/t3AXd1o18z6gJpjNhGxFfhZdxKNpMHA8cB5EdEaEYuAucCpNY4bQ+rlXL2jbZtZ35A7QLxA0qHdaGc8sCUilpbKFgMTahx3GrAwIpZ1o20z6wNyx2yWA7dJ+hnwOPCPdSki4vyM44cAGyrK1gNDaxx3GvCVjl6UNA2YBtDS0pIRhpn1ltxksztwc/F43x1opxUYVlE2DNjY0QGS3grsDdzUUZ2ImAXMApg0aZIX5jHrw7KSTUR8uJvtLAUGSBoXEY8UZROBJZ0c8yFgdkS0drNtM+sDsif1STpA0nmSLime7y/p4JxjI2ITMBuYKWmwpMOBKcCPO2hrd+BE4Ee58ZlZ35aVbCSdACwERpHGUSCNt3yjC23NIJ2OrQKuB6ZHxBJJkyVV9l6OA9YBv+zC+5tZH6acNYgl/RE4OSIWS1obEXtKGgisyJk52AiTJk0KT+ozazxJD0TEpFr1cgeIXwE8WDyO0vedblB2zNnzejuEfmHZ197b2yFYH5M7ZvMA20/AOxm4r2fDMbNmlduz+QRwp6SPAoMl3UGaqPfOukVmZk0l99L3w5IOAI4BbiVN7LvVl6XNLFdWspE0Cng2Im4sle0paZ+IWFG36MysaeSO2dzM9jOH9wXm9Gw4ZtascpPN+Ij4fbmgeH5Az4dkZs0oN9mslrRfuaB4/kzPh2RmzSg32fwQ+KmkYyQdKOlY0g2S369faGbWTHIvfX8NeAH4D+BVpKtR36drtyuYWT+We+l7K/B/iy8zsy7LXvBc0v6kZSGGlMsj4oc9HZSZNZ/ceTbnAueTlvJ8tvRSkMZzzMw6lduz+SRwSEQ8WLOmmVkVuVejngMermcgZtbccpPNecB3JL1S0i7lr3oGZ2bNI/c06kfF9zNKZSKN2ezakwGZWXPKTTZj6xqFmTW93Hk2WXv5mpl1pCvzbN4HHAHsRTqFAiAiTuvwIDOzQu7uCl8Crijqn0C6AfNdpB0QskgaLmmOpE2Slkua2kndN0haIKlV0t8k/WtuO2bWN+VeTfoIcHREfArYXHw/FhjThbYuBTYDI4FTgMskbbfXt6S9gNtJye3lwH7AnV1ox8z6oNxks0dEPFQ83ixpYETcRzqtqknSYOB44LyIaI2IRcBctl9EHeDTwB0RcW1EPB8RGyPij5lxmlkflZtsHi31Qh4Cpks6FVibefx4YEtELC2VLQa269kAhwJrJN0jaZWkWyS1ZLZjZn1U7gDxF0mnNADnANeSbsickXn8EGBDRdl60q6alfYF3gAcDfweuJi0g+bhlRUlTQOmAbS0OB+Z9WW5l75/Xnp8L2kcpStagWEVZcOAjVXqPgfMiYjfAEj6MvC0pJdFxPqKuGYBsyDtiNnFmMysgbpy6ftlwP5sv8TEXRmHLwUGSBoXEY8UZROBJVXqPkj7nTadRMyaQO4SE6eTria1sv0SE6+udXxEbJI0G5gp6QzgdcAU4LAq1a8kLUH6bVIyOg9YVNmrMbOdS+4A8YXAByJiZESMLX3VTDQlM4DdgVWkMZjpEbFE0mRJ/9jsrugpnQvMK+ruB3Q4J8fMdg65p1ED6OZcl4hYAxxXpXwh25+aXQZc1p32rHmNOXteb4fQLyz72nt79P1yezYXAV/0khJmtqM67NlIepxtg7MC9gb+TVK7vaIiwteczaymzk6jPtiwKMys6XWYbCJifiMDMbPmlnvpexBwOumSdeVgrpeYMLOacq9GXUWahHcL8Lf6hWNmzSo32bwbGBsR2evXmJmV5V7K/iuwWz0DMbPmltuzuRr4maRvUXEalXlvlJn1c7nJ5qzi+1cryrPujTIzy11iwlu5mFm3ZN9+IGlgcdPkScXzwcVyn2ZmNeXurvBa0po03wN+UBQfAfywTnGZWZPJ7dlcBpwfEQcALxRl84G31iUqM2s6uclmAnBN8TggLYhFWp/GzKym3GSzDHhjuUDSIcCfezogM2tOuZe+zwPmSbocGCTpHOBM4GN1i8zMmkpWzyYibiXdsjCCNFYzGnh/RHinSjPLkr27QkT8lvx9oszM2vEyn2bWEE42ZtYQDUs2koZLmiNpk6TlkqpuzyLpAkkvSGotffn+K7OdXPaYTQ+4FNgMjCSt+DdP0uKIqLYr5g0R4TWQzZpI7rKgw4HPUn1Z0LdlHD8YOB44KCJagUWS5gKnAmd3NWgz2/nk9myuIy2edSPtt9/NNR7YEhFLS2WLSfdXVXOspDXAU8AlxaZ1ZrYTy002hwEjIuL5HWxnCLChomw9MLRK3RuBWaRFut5M2vd7XURcX1lR0jRgGkBLi7evMuvLcgeIHwT27UY7rcCwirJhwMbKihHxh4hYEREvRsQ9wLeAD1R704iYFRGTImLSiBEjuhGemdVbbs/mLuB2SVcCK8svRETOMhNLgQGSxkXEI0XZRKDa4HClIO3IaWY7sdxkMxl4Aji6ojzIWNMmIjZJmg3MlHQGaaB5Cun0rB1JU4AFwDrgTcAngHMz4zSzPip3WdCjeqCtGaTEtAp4BpgeEUskTQZui4i2q1wnF/V2IyW4iyLiqh5o38x6UZfn2UgSpdOaiNiac1xErAGOq1K+kNLl9Ij4312Nycz6vtxlQUcVs3+fAbaQVutr+zIzqyn3atTlpNm/byddWXoDMJe0po2ZWU1dmWfTUgz0RkQslvRR4B7SIuhmZp3K7dm8SDp9AlgnaQSwCRhVl6jMrOnkJpt7gX8qHt8B3ADMBu6vR1Bm1nxyT6NOZVti+iTwGdKtBt+sR1Bm1nxy59msKz1+DvhK3SIys6bklfrMrCGcbMysIZxszKwhaiYbSbtKulvSbo0IyMyaU81kExEvAmNz6pqZdSQ3gXwZuEzS6KKns0vbVz2DM7PmkTvP5vvF91NLZSKtZ7Nrj0ZkZk0pN9mMrWsUZtb0cif1La93IGbW3LIXz5L0PtLWK3vRfvGs0+oQl5k1mdzFs74EXFHUP4G0rOe7SOsEm5nVlHs16SPA0RHxKWBz8f1YYEy9AjOz5pKbbPaIiIeKx5slDYyI++h4R0szs3Zyk82jkiYUjx8Cpks6FVib25Ck4cU6xpskLZc0tUb9QZL+KOmJ3DbMrO/KHSD+IvDy4vE5wLWkHRFmdKGtS0nrGI8k7Rs1T9LiiOhoo7rPAaupvkWvme1kci99/7z0+F5gv640ImkwcDxwUES0AoskzSVNEjy7Sv2xwAeBT+M1js2aQvbtBpIOkHSepEuK5/tLOjjz8PHAlohYWipbDEzooP53SLtgPpcbn5n1bbmXvk8gbYk7CmibVzMU+EZmO0OADRVl66lyiiTpn4FdI2JORlzTJN0v6f7Vq1dnhmJmvSG3ZzOTdOn7TNJOC5B6JhMzj28FhlWUDQM2lguK062LSft71xQRsyJiUkRMGjFiRGYoZtYbcgeIXwE8WDyO0veoXn07S4EBksZFxCNF2USgcnB4HGnuzsK0yy+DgJdJWgkcGhHLMtszsz4mt2fzAO3v+AY4Gbgv5+CI2ETa+mWmpMGSDgemAD+uqPoQ8CrS1arXAWcAfyseP54Zq5n1Qbk9m08Adxa7YA6WdAdp0PedXWhrBvBDYBXpdofpEbFE0mTgtogYEhFbgJVtB0haA2yNiJVV39HMdhq5l74flnQAcAxwK6mXcWtxGTtLRKwBjqtSvpA0gFztmLuBfXPbMLO+K/uu74h4FrixjrGYWRPLSjbFJLsLSWMn7XohEdFSh7jMrMnk9myuAx4lbbv7bP3CMbNmlZtsJgCHR8TWegZjZs0r99L3AuD19QzEzJpbbs9mGXC7pDmULk0DRMT5PR2UmTWf3GQzmHTJeyBp0l2b3BnEZtbP5c6z+XC9AzGz5uYdLc2sIZxszKwhnGzMrCGcbMysIXJX6juquGUBSa+UdJWkKyXtXd/wzKxZ5PZsvsu2Ffq+TroEvhWYVY+gzKz55M6zGRURf5U0gLTt7mjStiwr6haZmTWV3GSzQdJI4CDgDxHRKmkQqYdjZlZTbrL5DvAb0prAnyzKDgcerkdQZtZ8cmcQX1TcF/ViRDxaFD9JWiPYzKymrlz6/guwj6STiudPAo/1fEhm1oxyL32/lrQdy/eAHxTFR5AWMDczqym3Z3MZcH5EHAC8UJTNB95al6jMrOnkJpsJwDXF44B/7AW1e25DkoZLmiNpk6TlkqZ2UO9Tkh6TtEHSCkn/r7jkbmY7sdxkswx4Y7lA0iHAn7vQ1qWkuTkjgVOAyyRNqFJvLvCGiBhGutQ+kczteM2s78rtMZwHzJN0OTBI0jnAmcDHcg4u9vA+Hjio2GtqkaS5pF02zy7XLV3tAhBppvJ+mXGaWR+V1bOJiFuBdwMjSGM1o4H3R8Sdme2MB7ZExNJS2WLS6dl2JE2VtAF4mtSzuaKDetMk3S/p/tWrV2eGYma9oSub1P2WtIXujhgCbKgoWw8M7aCt64DrJI0DTiPt912t3iyK+7MmTZrkJUrN+rDcTeoGAadTfZO60zLeohUYVlE2DNjY2UER8YikJaQbQd+fE6uZ9U25PZurSKczt9BBL6OGpcAASeMi4pGibCKwJOPYAcBrdqBNM+tDcpPNu4GxEbFuRxqJiE2SZgMzJZ1B6iFNAQ6rrFu8PjciVkk6EDgHuGNH2jWzviP30vdfgd262dYM0rycVcD1wPSIWCJpsqTWUr3Dgd9L2gT8vPg6t5ttm1kvy+3ZXA38TNK3qDiNioi7ct4gItYAx1UpX0hpHMjbxpg1p9xkc1bx/asV5QG8uufCMbNmlbvExNh6B2JmzS17iQlJA4vxlZOK54OLmcFmZjV5iQkzawgvMWFmDdGwJSbMrH9r5BITZtaPNWSJCTOzRi0xYWb9XKOWmDCzfi53iYmZHbz0PPAEcHtE7Mjd4GbWT+QOEI8HPg8cRVqi86ji+euB6cBjkt5dlwjNrCnkJptdgJMjYnJETI2IycCJpB0yDyWdXn2tXkGa2c4vN9m8i7TrQdmtwHuKx9fgGzLNrBO5yeZR0ulS2ZlFOcBewLM9FZSZNZ/cq1FnALMlfZ60x/co4EW2rQu8P2kujplZVblLTPxPsdPBocA+wFPAryLiheL1BcCCukVpZju9rsyzeQFYWMdYzKyJZa9nY2bWHU42ZtYQHSYbSe8rPR7Y3YYkDZc0R9ImScslTe2g3uckPSRpo6S/SPpcd9s2s97X2ZjNNWzbxfIZtt/RsqsuBTYDI0n7Rs2TtDgiKjeqE2nL3QdJm9PdKenxiPjPbrZvZr2os2SzUtJZwB9Iu1keRUoE7eRs5VKsVXw8cFBEtAKLJM0FTgXOrni/i0tP/yTpZ6S9pJxszHZinSWb04GZwL8Cg6i+3nDuVi7jgS0RsbRUtpi0jnGHJAmYDFyR0YaZ9WEdJpuIuAd4B4CkP0fEft1oZwiwoaJsPTC0xnEXkMaVrqz2oqRpwDSAlpaWboRnZvWWu3jWfgCSWiS9RdKruthOK9uP+QwDNnZ0QHEKdxrw3oh4voO4ZkXEpIiYNGLEiC6GZGaNlLuVy96S5pPWHJ4NPCppgaR9MttZShr3GVcqmwhUDg63tfcR0ljO2yPiicw2zKwPy51nczlpjGXPiHglsCfw26K8pmInhtnAzGJzu8OBKcCPK+tKOoW0ze/REfFYZnxm1sflJpu3Ap8pkkZb8vg34LAutDWDtPXLKuB6YHpELCl22Wwt1fsK8HLgN5Jai6+spGZmfVfuvVFrgQNJvZs2+wPrchuKiDXAcVXKF5IGkNuee19xsyaUm2wuBv5L0g+A5aTdFT6Ml5Uws0y5S0x8T9KjwFTgYGAFMDUi/ruewZlZ8+jKEhN3ATVnC5uZVeO7vs2sIZxszKwhnGzMrCGcbMysIXK33x0OfJa0Ds2Q8msR8bY6xGVmTSb3atR1wG7AjXh/KDPbAbnJ5jBgREd3X5uZ1ZI7ZvMgsG89AzGz5pbbs7kLuF3SlcDK8gsRUW0FPzOzdnKTzWTgCeDoivKg+nKhZmbt5N4bdVS9AzGz5pZ9b5SkPYFjgVHAk8AtEbG2XoGZWXPJXRb0LcCjwJmku77/hbQ06FvqGJuZNZHcns03gRnljeIknQR8G3hTPQIzs+aSe+l7PGlCX9lNQHe2dzGzfiQ32TwCnFxRdgLp1MrMrKbc06hPArdK+gRpWdAxwDjgmDrFZWZNJvfS9z2SXgO8F9gHuAX4ebGIuZlZTdlLTETE2oi4JiIuLr53KdFIGi5pjqRNkpZLmtpBvaMk/VLSeknLutKGmfVdHfZsJN0eEe8uHi8kzRbeTheWmLgU2AyMJC1VMU/S4oio3BVzE2lW8vXAuZnvbWZ9XGenUVeXHn+/O41IGgwcDxwUEa3AIklzgVNJ2+z+Q0TcB9wn6R3dadPM+pYOk01EXFd6+nBE3FtZR9Ihme2MB7ZExNJS2WLgiMzjq5I0DZgG0NLS0p23MrM6yx2z+UUH5bdnHj8E2FBRth4Ymnl8VRExKyImRcSkESNGdOetzKzOOr0aJWkXQOmhVDxu8xpgS2Y7rcCwirJhwMbM481sJ1fr0vcWtg0MVyaWrcCFme0sBQZIGhcRjxRlE4HKwWEza1K1ks1YUm9mPlC+6hTA6oh4LqeRiNgkaTYwU9IZpKtRU0jLjbZT9KYGAQPTU70E2BoRm3PaMrO+qdNkExHLASTtD7wYES+0vSZpoKTdurAu8QzSJe1VwDPA9IhYImkycFtEtO3a8Dbgl6XjniMluyMz2zGzPih3gPhO4I0VZW8E7shtKCLWRMRxETE4IlrarnZFxMJSoiEi7o4IVXwdmduOmfVNucnmYKDy0vd9pHEXM7OacpPNOtLM37KRpNm+ZmY15SabnwLXSTpI0kslvZY0w7hyjRszs6pyk80XgD+STp02Ar8G/oTvXTKzTLlLTPwd+D+SzgL2Ap6OiKo3ZpqZVZO9u0JhSPE1NE0ohoh4rKeDMrPmk5VsJB0IXEu6+hSkiX5tPZtd6xOamTWT3DGb75Im2g0n3VC5J3AF8KE6xWVmTSb3NGoicHREvCBJEbFe0ueAh4Br6heemTWL3J7N30n3KgE8LamlOPbldYnKzJpObrJZCJxYPL4JuI10v9Jd9QjKzJpP7qXvE0tPzyWdPg2l/dKhZmYdqplsJO0K/Dfwroh4PiK24nEaM+uimqdREfEiaV2b7G1fzMwq5SaQLwOXSRotaVdJu7R91TM4M2seuZe+27ZyObVU1jaxz5P6zKym3GQztq5RmFnTq7W7wt4RsbJteVAzsx1Va8ylvKkcxaLlZmZdVivZqOL5kXWKw8yaXK1k02Nr1kgaLmmOpE2Slkua2kE9SbpI0jPF10VqW8/CzHZatQaIB0g6im09nMrnRETuLQuXAptJaxe/DpgnaXFEVG5UNw04jm3LWfwC+AtweWY7ZtYH1Uo2q0h7PbV5puJ5AK+u1YikwcDxwEER0QoskjSXdCn97IrqHwK+HhFPFMd+HfgYTjZmO7Vam9SN6aF2xgNbIqI84LwYOKJK3QnFa+V6E3ooDjPrJV1dFnRHDSEtulW2nnQzZ7W66yvqDSnW0Wk3hiRpGum0C6BV0p96KN6+Zi/g6d4Ooit0UW9H0Oc082c4OqdSo5JNKzCsomwYaaeGWnWHAa3VFliPiFnArJ4Ksq+SdH9ETOrtOGzH+TNs3M2VS0mDy+NKZROBysFhirKJGfXMbCfSkGQTEZuA2cBMSYMlHQ5MAX5cpfrVwKcljZK0D/AZ4EeNiNPM6qeRd23PAHYnXeG6HpgeEUskTZbUWqp3BXAL8HvSIl3zirL+rOlPFfuBfv8ZynvNmVkjeD0aM2sIJ5t+RNIySe9oUFs/kvSVRrTVTBr5e5N0uqRFjWgL+lGyKf6hPSdpo6R1ku6RdKZXG6xOUkjar7fj6KpqcUu6QFLTrZst6W5JZ/R2HLn62z+0YyNiKGkS0teAzwM/6N2QtpHUqHlPtoP8Ge24/pZsAIiI9RExFzgJ+JCkgyS9TNLVklYXd6V/sa3XUzx/Y/H4lOJ/zwnF849Kurl4fIGkG4v32ShpiaQOJ3IV9W+SdI2kDcDpxdrOZ0t6tLjr/UZJw0vHnFbE84yk88qnRpVdcElHSnqig7YPkfSropf3lKRLJA0qXltQVFssqVXSSUX5MZJ+V+oZHlx6v9dL+p/i574BeEmXP5gGaPudSDpX0tPF7++UjPqfl7QSuLIo7+x38QZJvy1+Fz+RdEPb51Lt1KWjXqSkPSXdWvxNri0e71u8diEwGbik+IwuKcoPkPQLSWsk/UnSiaX3e7mkuZI2SLoPeE13fpdd1S+TTZuIuA94gvShfQd4GenG0iOA04APF1Xns20tnyOAx4C3lZ7PL73t+4D/BPYA5gKX1AhjCmnjvz2Aa4GPk+56PwLYB1hLumMeSQeS9l0/BXhlEe+orvzMJS8CnyJNo38L8HbS9AQiou1nmxgRQyLiBkmvJ92E+y+knVCvAOZK2q1IUjeT5k0NB35CuvG2r9qb9HOPIt34O0vS/jXqDyf1iKdl/C7mkOaGDSdN8/jnHYxzF1JyGw20AM9R/D1FxBdIm0eeVXxGZynd8PwL4DrgFcDJwHeLvxtIf0d/J/3tfKT4aph+nWwKK0h/FCcD50TExohYBnydbQu8z2fbTaOTgX8vPa9MNosi4ufFFjg/pv1s6Gp+FRE3R8TWiHgOOBP4QkQ8ERHPAxcAHyi67x8AbomIRRGxGTifHVxzKCIeiIhfR8SW4ue9guo3xraZBlwREfdGxIsRcRXwPHBo8TUQ+GZEvBARNwG/2ZG4Gui8Yh+0+aS5XCd2Uncr8KWi/nPU/l0MAL5d/C5mA/ftSIAR8UxE/DQino2IjcCFdP4ZHQMsi4gri8/1t8BPgROU9n87Hjg/IjZFxEPAVTsS145yskn/uw0g/WMpr7W8nG29hvnAZEmvJO0mcSNwuKQxpN7F70rHrSw9fhZ4iaQBxelXa/F1W6nO4xXxjAbmFN3zdcAfSb2QkaSezj/qR8SzpGU/ukzS+KJbvrI4hfsq6X/7jowGPtMWVxHbq4qY9gGerLh/rbfWrX6RbfvStxkIvFB6vraY1d5mObCPpJbSZ1SeaLo6Iv5eet7V30XlZ5xF0kslXVGcNm8AFgB7FImjmtHAmyviOoXUMxtB+jsvx9LQz6hfJxtJbyIllJtJf4zlu1dbgCcBIuLPpMTxcWBBRGwgJZVppJ7M1lptRcS1RXd3SES8p/xSRdXHgfdExB6lr5dExJPAU8C+pfh3J3Xj22wCXlp6vncnIV0GPAyMi4hhpG2VO1sR8XHgwoq4XhoR1xdxjZLarajY0sl71dNfgTEVZWNp/w9rz+KUo00LsCIi/lr6jIaUXq/2GXXld/Gq0uN2n5Gkzj6jzwD7A28uPqO209u2964W1/yKuIZExHRgNbClIpaGfkb9MtlIGibpGNLYyjURsZjUW7lQ0lBJo4FP036b4fnAWWw7Zbq74nlPubyIY3QR6whJU4rXbgKOlXRYMTZwAe0TxO+Af1JagnVv4JOdtDOUtOxHq6QDgOkVr/+N9gujfQ84U9KblQyW9F5JQ4Ffkf6QPyFpoKT3A4fswM/eE24AvihpX6XB9ncAx5J+d2VfljRI0mTS6cdPutBGrd/Fi8BZRY92Cu1/F4uBCZJeJ+klpM+wI0NJ4zTrlC4SfKni9crP6FZgvKRTi89hoKQ3SfpfxWn9bOCCosd0IGm8qmH6W7K5RdJG0v8AXwC+wbZB4I+T/td5DFhEGmQrr0o4n/ThL+jgeU/5Fmlg+c4i1l8DbwYollD9OClJPkVajmMVabwA0hjRYmAZcCfpH15HPgtMJS3z8b0qdS8Ariq64ydGxP2kFRMvIQ1a/xk4vYhrM/D+4vka0lW+3tqJYyZwD+kzXAtcDJxSjFG0WVm8toI0KH9mRDyc20Dm7+KjwDrgg6Qk8Hzx+tIixv8CHini7Mg3SfcTPk36O7i94vVvkcbz1kr6djGu807S+OOK4ue8CNitqH8Wab2olaQB7Ctzf+ae4HujdmKShpD+oMdFxF96O56dgaQjSb3ZfWvV7cE27wUuj4iG/uPua/pbz2anJ+nYohs8GPgP0t3xy3o3KiuTdISkvYvTqA8BB7N9r6TfcbLZ+UwhdZFXAOOAk6utYmi9an/S6ew60iDvByLiqd4Nqff5NMrMGsI9GzNrCCcbM2sIJxszawgnGzNrCCcbM2sIJxsza4j/D55RyIz2DDxEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Silencers\n",
    "silencer_count_contingency = de_direction_grouper.apply(lambda x: activity_has_gene_df.loc[x.index, \"group_name_WT\"].str.contains(\"Silencer\").value_counts()).unstack()\n",
    "oddsratio, pval = stats.fisher_exact(silencer_count_contingency)\n",
    "print(f\"Direction of differential expression is independent of having a silencer nearby, Fisher's exact test p={pval:.3f}, OR={oddsratio:.2f}\")\n",
    "\n",
    "silencer_count_contingency = silencer_count_contingency.div(silencer_count_contingency.sum(axis=1), axis=0)\n",
    "display(silencer_count_contingency)\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar([0, 1], silencer_count_contingency[True], tick_label=silencer_count_contingency.index)\n",
    "ax.set_ylabel(\"Fraction of genes near a silencer\")\n",
    "\n",
    "# Enhancers\n",
    "enhancer_count_contingency = de_direction_grouper.apply(lambda x: activity_has_gene_df.loc[x.index, \"group_name_WT\"].str.contains(\"enhancer\").value_counts()).unstack()\n",
    "oddsratio, pval = stats.fisher_exact(enhancer_count_contingency)\n",
    "print(f\"Direction of differential expression is independent of having an enhancer nearby, Fisher's exact test p={pval:.3f}, OR={oddsratio:.2f}\")\n",
    "\n",
    "enhancer_count_contingency = enhancer_count_contingency.div(enhancer_count_contingency.sum(axis=1), axis=0)\n",
    "display(enhancer_count_contingency)\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar([0, 1], enhancer_count_contingency[True], tick_label=enhancer_count_contingency.index)\n",
    "ax.set_ylabel(\"Fraction of genes near an enhancer\")"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}