{ "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 }