Jupyter_notebooks/Supplementary_Figure_2/Motif_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Motif analyses on CLIP/RIP binding sites and DRNs in diffBUM-HMM and deltaSHAPE data.\n",
    "For this to work you need to have the pyCRAC package and MEME tool suite installed (refer to README.md).\n",
    "To install the MEME tool suite, please follow the instructions for quick install here: http://meme-suite.org/doc/install.html?man_type=web"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: numpy in /anaconda3/lib/python3.7/site-packages (1.18.2)\n",
      "Requirement already satisfied: pandas in /anaconda3/lib/python3.7/site-packages (0.25.1)\n",
      "Requirement already satisfied: seaborn in /anaconda3/lib/python3.7/site-packages (0.9.0)\n",
      "Requirement already satisfied: matplotlib in /anaconda3/lib/python3.7/site-packages (3.2.1)\n",
      "Requirement already satisfied: scipy in /anaconda3/lib/python3.7/site-packages (1.1.0)\n",
      "Requirement already satisfied: pyCRAC in /anaconda3/lib/python3.7/site-packages (1.4.7)\n",
      "Requirement already satisfied: python-dateutil>=2.6.1 in /anaconda3/lib/python3.7/site-packages (from pandas) (2.7.5)\n",
      "Requirement already satisfied: pytz>=2017.2 in /anaconda3/lib/python3.7/site-packages (from pandas) (2018.7)\n",
      "Requirement already satisfied: cycler>=0.10 in /anaconda3/lib/python3.7/site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /anaconda3/lib/python3.7/site-packages (from matplotlib) (2.3.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /anaconda3/lib/python3.7/site-packages (from matplotlib) (1.0.1)\n",
      "Requirement already satisfied: pysam>=0.6 in /anaconda3/lib/python3.7/site-packages (from pyCRAC) (0.15.4)\n",
      "Requirement already satisfied: six>=1.9.0 in /anaconda3/lib/python3.7/site-packages (from pyCRAC) (1.12.0)\n",
      "Requirement already satisfied: cython>=0.19 in /anaconda3/lib/python3.7/site-packages (from pyCRAC) (0.29.2)\n",
      "Requirement already satisfied: setuptools in /anaconda3/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib) (40.6.3)\n"
     ]
    }
   ],
   "source": [
    "! pip install numpy pandas seaborn matplotlib scipy pyCRAC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Bad key \"text.kerning_factor\" on line 4 in\n",
      "/anaconda3/lib/python3.7/site-packages/matplotlib/mpl-data/stylelib/_classic_test_patch.mplstyle.\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "http://github.com/matplotlib/matplotlib/blob/master/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import math\n",
    "import sys\n",
    "import glob\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy.special import comb\n",
    "from more_itertools import pairwise\n",
    "from scipy.stats import hypergeom,chisquare,fisher_exact\n",
    "from matplotlib import rcParams\n",
    "from collections import defaultdict\n",
    "from pyCRAC.Methods import sortbyvalue,contigousarray2Intervals\n",
    "\n",
    "rcParams['font.family'] = 'sans-serif'\n",
    "rcParams['font.sans-serif'] = ['Arial']\n",
    "rcParams['pdf.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalizeIntervalLength(start,end,chromosomelength,length=20):\n",
    "    \"\"\" Extends or trims the interval coordinates to a set length. Default = 20 bp. \"\"\"\n",
    "    newstart = int()\n",
    "    newend = int()\n",
    "    if start - length < 0: start = 1 # to make sure that the start is always a positive number \n",
    "    if end + length > chromosomelength: end = chromosomelength # to make sure that interval doesn't go beyond the chromosome boundaries.\n",
    "    actlength  = end - start\n",
    "    difference = length - actlength\n",
    "    if difference == 0:\n",
    "        return start,end\n",
    "    else:\n",
    "        newstart = round(float(start) - float(difference)/2.0)\n",
    "        if newstart < 0: newstart = 1\n",
    "        newend = round(float(end) + float(difference)/2.0)\n",
    "        if newend > chromosomelength: newend = chromosomelength\n",
    "        return int(newstart),int(newend) # convert back to integers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def intervalsFromClosePositions(positions,mindistance=5):\n",
    "    \"\"\" Merges positions that are within a specific distance\n",
    "    to their neighbours:\n",
    "    [1,3,5,15,20,30,35,36,37,69,70,80,90,91]\n",
    "    should become:\n",
    "    [(1, 5), (15, 20), (30, 37), (69, 70), (80, 80), (90, 91)]\n",
    "\n",
    "    \"\"\"\n",
    "    start = None\n",
    "    end\t  = None\n",
    "    intervallist = list()\n",
    "    for nr,i in enumerate(sorted(list(positions))):\n",
    "        if not start:\n",
    "            start = i\n",
    "        try:\n",
    "            if positions[nr+1] - positions[nr] <= mindistance:\n",
    "                continue\n",
    "            elif positions[nr+1] - positions[nr] > mindistance:\n",
    "                end = i\n",
    "                intervallist.append((start,end))\n",
    "                start = positions[nr+1]\n",
    "        except IndexError:\n",
    "            if start:\n",
    "                end = i\n",
    "                intervallist.append((start,end))\n",
    "            break\n",
    "    return intervallist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the big dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('../../Data/New_data_table_Xist.txt',\\\n",
    "                    sep=\"\\t\",\\\n",
    "                    header=0,\\\n",
    "                    index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking positions not considered by deltaSHAPE:\n",
    "Any positions that have -999 values should not be further considered as these were positions where there was insufficient coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2517\n"
     ]
    }
   ],
   "source": [
    "positionstomask = data[(data[\"SHAPE_reactivity_ex_vivo_1\"] < -900) | \n",
    "                       (data[\"SHAPE_reactivity_ex_vivo_2\"] < -900) |\n",
    "                       (data[\"SHAPE_reactivity_in_cell_1\"] < -900) | \n",
    "                       (data[\"SHAPE_reactivity_in_cell_2\"] < -900)].index\n",
    "print(len(positionstomask))\n",
    "data.loc[positionstomask,data.columns[11:]] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['nucleotide', 'CELF1', 'FUS', 'HuR', 'PTBP1', 'RBFOX2', 'TARDBP',\n",
       "       'SHAPE_reactivity_ex_vivo_1', 'SHAPE_reactivity_ex_vivo_2',\n",
       "       'SHAPE_reactivity_in_cell_1', 'SHAPE_reactivity_in_cell_2',\n",
       "       'deltaSHAPE_rep1', 'deltaSHAPE_rep2', 'dStruct', 'dStruct_scaled',\n",
       "       'unmodified', 'ex_vivo', 'in_vivo', 'both', 'scaled_unmodified',\n",
       "       'scaled_ex_vivo', 'scaled_in_vivo', 'scaled_both', 'base_paired_to'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>nucleotide</th>\n",
       "      <th>CELF1</th>\n",
       "      <th>FUS</th>\n",
       "      <th>HuR</th>\n",
       "      <th>PTBP1</th>\n",
       "      <th>RBFOX2</th>\n",
       "      <th>TARDBP</th>\n",
       "      <th>SHAPE_reactivity_ex_vivo_1</th>\n",
       "      <th>SHAPE_reactivity_ex_vivo_2</th>\n",
       "      <th>SHAPE_reactivity_in_cell_1</th>\n",
       "      <th>...</th>\n",
       "      <th>dStruct_scaled</th>\n",
       "      <th>unmodified</th>\n",
       "      <th>ex_vivo</th>\n",
       "      <th>in_vivo</th>\n",
       "      <th>both</th>\n",
       "      <th>scaled_unmodified</th>\n",
       "      <th>scaled_ex_vivo</th>\n",
       "      <th>scaled_in_vivo</th>\n",
       "      <th>scaled_both</th>\n",
       "      <th>base_paired_to</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>C</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>G</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>G</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>C</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5</td>\n",
       "      <td>T</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  nucleotide  CELF1  FUS  HuR  PTBP1  RBFOX2  TARDBP  \\\n",
       "1          C      0    0    0      0       0       0   \n",
       "2          G      0    0    0      0       0       0   \n",
       "3          G      0    0    0      0       0       0   \n",
       "4          C      0    0    0      0       0       0   \n",
       "5          T      0    0    0      0       0       0   \n",
       "\n",
       "   SHAPE_reactivity_ex_vivo_1  SHAPE_reactivity_ex_vivo_2  \\\n",
       "1                      -999.0                      -999.0   \n",
       "2                      -999.0                      -999.0   \n",
       "3                      -999.0                      -999.0   \n",
       "4                      -999.0                      -999.0   \n",
       "5                      -999.0                      -999.0   \n",
       "\n",
       "   SHAPE_reactivity_in_cell_1  ...  dStruct_scaled  unmodified  ex_vivo  \\\n",
       "1                      -999.0  ...             NaN         NaN      NaN   \n",
       "2                      -999.0  ...             NaN         NaN      NaN   \n",
       "3                      -999.0  ...             NaN         NaN      NaN   \n",
       "4                      -999.0  ...             NaN         NaN      NaN   \n",
       "5                      -999.0  ...             NaN         NaN      NaN   \n",
       "\n",
       "   in_vivo  both  scaled_unmodified  scaled_ex_vivo  scaled_in_vivo  \\\n",
       "1      NaN   NaN                NaN             NaN             NaN   \n",
       "2      NaN   NaN                NaN             NaN             NaN   \n",
       "3      NaN   NaN                NaN             NaN             NaN   \n",
       "4      NaN   NaN                NaN             NaN             NaN   \n",
       "5      NaN   NaN                NaN             NaN             NaN   \n",
       "\n",
       "   scaled_both  base_paired_to  \n",
       "1          NaN             NaN  \n",
       "2          NaN             NaN  \n",
       "3          NaN             NaN  \n",
       "4          NaN             NaN  \n",
       "5          NaN             NaN  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now doing Motif analyses:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First I need to make fasta files for all the protein binding sites in Xist and find their motifs:\n",
    "I will normalize the interval size to a fixed lenght so that MEME can deal better with it. I will group together those DRNs that are within 2nt distance from each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['nucleotide', 'CELF1', 'FUS', 'HuR', 'PTBP1', 'RBFOX2', 'TARDBP',\n",
       "       'SHAPE_reactivity_ex_vivo_1', 'SHAPE_reactivity_ex_vivo_2',\n",
       "       'SHAPE_reactivity_in_cell_1', 'SHAPE_reactivity_in_cell_2',\n",
       "       'deltaSHAPE_rep1', 'deltaSHAPE_rep2', 'dStruct', 'dStruct_scaled',\n",
       "       'unmodified', 'ex_vivo', 'in_vivo', 'both', 'scaled_unmodified',\n",
       "       'scaled_ex_vivo', 'scaled_in_vivo', 'scaled_both', 'base_paired_to'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>nucleotide</th>\n",
       "      <th>CELF1</th>\n",
       "      <th>FUS</th>\n",
       "      <th>HuR</th>\n",
       "      <th>PTBP1</th>\n",
       "      <th>RBFOX2</th>\n",
       "      <th>TARDBP</th>\n",
       "      <th>SHAPE_reactivity_ex_vivo_1</th>\n",
       "      <th>SHAPE_reactivity_ex_vivo_2</th>\n",
       "      <th>SHAPE_reactivity_in_cell_1</th>\n",
       "      <th>...</th>\n",
       "      <th>dStruct_scaled</th>\n",
       "      <th>unmodified</th>\n",
       "      <th>ex_vivo</th>\n",
       "      <th>in_vivo</th>\n",
       "      <th>both</th>\n",
       "      <th>scaled_unmodified</th>\n",
       "      <th>scaled_ex_vivo</th>\n",
       "      <th>scaled_in_vivo</th>\n",
       "      <th>scaled_both</th>\n",
       "      <th>base_paired_to</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>C</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>G</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>G</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>C</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5</td>\n",
       "      <td>T</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>-999.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  nucleotide  CELF1  FUS  HuR  PTBP1  RBFOX2  TARDBP  \\\n",
       "1          C      0    0    0      0       0       0   \n",
       "2          G      0    0    0      0       0       0   \n",
       "3          G      0    0    0      0       0       0   \n",
       "4          C      0    0    0      0       0       0   \n",
       "5          T      0    0    0      0       0       0   \n",
       "\n",
       "   SHAPE_reactivity_ex_vivo_1  SHAPE_reactivity_ex_vivo_2  \\\n",
       "1                      -999.0                      -999.0   \n",
       "2                      -999.0                      -999.0   \n",
       "3                      -999.0                      -999.0   \n",
       "4                      -999.0                      -999.0   \n",
       "5                      -999.0                      -999.0   \n",
       "\n",
       "   SHAPE_reactivity_in_cell_1  ...  dStruct_scaled  unmodified  ex_vivo  \\\n",
       "1                      -999.0  ...             NaN         NaN      NaN   \n",
       "2                      -999.0  ...             NaN         NaN      NaN   \n",
       "3                      -999.0  ...             NaN         NaN      NaN   \n",
       "4                      -999.0  ...             NaN         NaN      NaN   \n",
       "5                      -999.0  ...             NaN         NaN      NaN   \n",
       "\n",
       "   in_vivo  both  scaled_unmodified  scaled_ex_vivo  scaled_in_vivo  \\\n",
       "1      NaN   NaN                NaN             NaN             NaN   \n",
       "2      NaN   NaN                NaN             NaN             NaN   \n",
       "3      NaN   NaN                NaN             NaN             NaN   \n",
       "4      NaN   NaN                NaN             NaN             NaN   \n",
       "5      NaN   NaN                NaN             NaN             NaN   \n",
       "\n",
       "   scaled_both  base_paired_to  \n",
       "1          NaN             NaN  \n",
       "2          NaN             NaN  \n",
       "3          NaN             NaN  \n",
       "4          NaN             NaN  \n",
       "5          NaN             NaN  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "chromosomelength = len(data.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting the length of each interval for motif analyses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixedlength = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting the threshold for selecting DRNs in the diffBUM-HMM data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins = [\"CELF1\",\"FUS\",\"HuR\",\"PTBP1\",\"RBFOX2\",\"TARDBP\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CELF1:\t21\n",
      "FUS:\t107\n",
      "HuR:\t12\n",
      "PTBP1:\t22\n",
      "RBFOX2:\t3\n",
      "TARDBP:\t17\n"
     ]
    }
   ],
   "source": [
    "proteinfastas = list()\n",
    "for protein in proteins:\n",
    "    outfilename = \"%s_Xist_RNA_binding_sites.fa\" % protein\n",
    "    #print(outfilename)\n",
    "    outfile = open(outfilename,\"w\")\n",
    "    proteinfastas.append(outfilename)\n",
    "    indices = data[data[protein] > 0].index\n",
    "    intervals = intervalsFromClosePositions(indices)\n",
    "    print(\"%s:\\t%s\" % (protein,len(intervals)))\n",
    "    count = 1\n",
    "    for (i,j) in intervals:\n",
    "        start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "        sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "        sequence = \"\".join(sequence)\n",
    "        outfile.write(\">%s_%s\\n%s\\n\" % (protein,count,sequence))\n",
    "        count += 1\n",
    "    outfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteinfastas = \" \".join(proteinfastas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CELF1_Xist_RNA_binding_sites.fa FUS_Xist_RNA_binding_sites.fa HuR_Xist_RNA_binding_sites.fa PTBP1_Xist_RNA_binding_sites.fa RBFOX2_Xist_RNA_binding_sites.fa TARDBP_Xist_RNA_binding_sites.fa\n"
     ]
    }
   ],
   "source": [
    "print(proteinfastas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running meme on these sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$proteinfastas\"\n",
    "\n",
    "DIR=$HOME/meme/bin\n",
    "\n",
    "for i in $1\n",
    "do\n",
    "    echo $i\n",
    "    FILE=\"$(basename -- $i)\"\n",
    "    FILE=${FILE%.fa}\n",
    "    PREFIX=MEME_\n",
    "    OUTFILE=$PREFIX$FILE\n",
    "    $DIR/meme-chip \\\n",
    "        -meme-minw 4 \\\n",
    "        -meme-maxw 10 \\\n",
    "        -meme-nmotifs 3 \\\n",
    "        -meme-p 8 \\\n",
    "        -meme-mod anr \\\n",
    "        -norc \\\n",
    "        -rna \\\n",
    "        -noecho \\\n",
    "        -oc $OUTFILE $i &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now doing it for the diffBUM_HMM and deltaSHAPE sites. Fragment sizes = fixedlength:\n",
    "All the DRNs that were located within a 5 nucleotide window from each other were grouped into a single interval and length of this interval was normalized (so either extended or trimmed) to the fixed length (30 bp for these analyses). The reason for doing this is because otherwise the analyes would generate for each DRN a sequence containing the DRN in the middle and many of these sequences could overlap if the DRNs are in quite close proximity. This would result in a much higher enrichment of motifs from regions that have a high concentration of DRNs. So this is why they were first grouped together. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    " intervallengthcounter = defaultdict(lambda: defaultdict(int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### diffBUM_HMM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "572\n",
      "1164\n",
      "108\n",
      "142\n"
     ]
    }
   ],
   "source": [
    "sequence = \"\".join(data.nucleotide.values)\n",
    "\n",
    "outfile = open(\"diffBUM_HMM_ex_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "\n",
    "### How many DRNs are there in the data?:\n",
    "ex_vivo_pos = data[data.ex_vivo >= threshold].index\n",
    "\n",
    "### intervalsFromClosePositions groups the DRNs together in intervals.\n",
    "intervals = intervalsFromClosePositions(ex_vivo_pos,)\n",
    "\n",
    "### This prints the number of intervals that were detected and how many DRNs there were in the data.\n",
    "print(len(intervals))\n",
    "print(len(ex_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"diffBUM_HMM_ex_vivo\"][length] += 1\n",
    "    ### These intervals are then set to a fixed length here:\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">diffBUM_HMM_ex_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()\n",
    "        \n",
    "in_vivo_pos = data[data.in_vivo >= threshold].index\n",
    "outfile = open(\"diffBUM_HMM_in_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "intervals = intervalsFromClosePositions(in_vivo_pos)\n",
    "\n",
    "print(len(intervals))\n",
    "print(len(in_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"diffBUM_HMM_in_vivo\"][length] += 1\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">diffBUM_HMM_in_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### deltaSHAPE rep 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n",
      "460\n",
      "91\n",
      "337\n"
     ]
    }
   ],
   "source": [
    "sequence = \"\".join(data.nucleotide.values)\n",
    "\n",
    "ex_vivo_pos = data[data.deltaSHAPE_rep1 > 0].index\n",
    "outfile = open(\"deltaSHAPE_rep_1_ex_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "intervals = intervalsFromClosePositions(ex_vivo_pos)\n",
    "\n",
    "print(len(intervals))\n",
    "print(len(ex_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"deltaSHAPE_rep1_ex_vivo\"][length] += 1\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">deltaSHAPE_rep1_ex_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()\n",
    "        \n",
    "in_vivo_pos = data[data.deltaSHAPE_rep1 < 0].index\n",
    "outfile = open(\"deltaSHAPE_rep_1_in_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "intervals = intervalsFromClosePositions(in_vivo_pos)\n",
    "\n",
    "print(len(intervals))\n",
    "print(len(in_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"deltaSHAPE_rep1_in_vivo\"][length] += 1\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">deltaSHAPE_rep1_in_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### deltaSHAPE rep 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "152\n",
      "510\n",
      "123\n",
      "428\n"
     ]
    }
   ],
   "source": [
    "sequence = \"\".join(data.nucleotide.values)\n",
    "\n",
    "ex_vivo_pos = data[data.deltaSHAPE_rep2 > 0].index\n",
    "outfile = open(\"deltaSHAPE_rep_2_ex_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "intervals = intervalsFromClosePositions(ex_vivo_pos)\n",
    "\n",
    "print(len(intervals))\n",
    "print(len(ex_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"deltaSHAPE_rep2_ex_vivo\"][length] += 1\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">deltaSHAPE_rep2_ex_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()\n",
    "        \n",
    "in_vivo_pos = data[data.deltaSHAPE_rep2 < 0].index\n",
    "outfile = open(\"deltaSHAPE_rep_2_in_vivo_%s_mers.fa\" % fixedlength,\"w\")\n",
    "intervals = intervalsFromClosePositions(in_vivo_pos)\n",
    "\n",
    "print(len(intervals))\n",
    "print(len(in_vivo_pos))\n",
    "\n",
    "count = 0\n",
    "\n",
    "for (i,j) in intervals:\n",
    "    length = (j-i)+1\n",
    "    intervallengthcounter[\"deltaSHAPE_rep2_in_vivo\"][length] += 1\n",
    "    start,end = normalizeIntervalLength(i,j,chromosomelength,fixedlength)\n",
    "    sequence = data.loc[start:end+1,\"nucleotide\"].values\n",
    "    sequence = \"\".join(sequence)\n",
    "    outfile.write(\">deltaSHAPE_rep2_in_vivo_%s\\n%s\\n\" % (count,sequence))\n",
    "    count += 1\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmerfasta = \" \".join(glob.glob(\"*%s_mers.fa\" % fixedlength))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['diffBUM_HMM_ex_vivo_30_mers.fa', 'deltaSHAPE_rep_1_in_vivo_30_mers.fa', 'deltaSHAPE_rep_2_in_vivo_30_mers.fa', 'deltaSHAPE_rep_1_ex_vivo_30_mers.fa', 'diffBUM_HMM_in_vivo_30_mers.fa', 'deltaSHAPE_rep_2_ex_vivo_30_mers.fa']\n"
     ]
    }
   ],
   "source": [
    "print(kmerfasta.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$kmerfasta\"\n",
    "\n",
    "echo $1\n",
    "\n",
    "DIR=$HOME/meme/bin\n",
    "\n",
    "for i in $1\n",
    "do\n",
    "    FILE=\"$(basename -- $i)\"\n",
    "    FILE=${FILE%.fa}\n",
    "    PREFIX=MEME_V2\n",
    "    OUTFILE=$PREFIX$FILE\n",
    "    $DIR/meme-chip \\\n",
    "        -meme-minw 4 \\\n",
    "        -meme-maxw 10 \\\n",
    "        -meme-nmotifs 20 \\\n",
    "        -meme-p 8 \\\n",
    "        -meme-mod anr \\\n",
    "        -norc \\\n",
    "        -rna \\\n",
    "        -noecho \\\n",
    "        -o $OUTFILE $i &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at the lenghts of the intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['diffBUM_HMM_ex_vivo',\n",
       " 'diffBUM_HMM_in_vivo',\n",
       " 'deltaSHAPE_rep1_ex_vivo',\n",
       " 'deltaSHAPE_rep1_in_vivo',\n",
       " 'deltaSHAPE_rep2_ex_vivo',\n",
       " 'deltaSHAPE_rep2_in_vivo']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(intervallengthcounter.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(<function __main__.<lambda>()>,\n",
       "            {'diffBUM_HMM_ex_vivo': defaultdict(int,\n",
       "                         {1: 297,\n",
       "                          6: 23,\n",
       "                          8: 20,\n",
       "                          2: 76,\n",
       "                          5: 24,\n",
       "                          9: 12,\n",
       "                          4: 38,\n",
       "                          14: 6,\n",
       "                          11: 6,\n",
       "                          3: 43,\n",
       "                          12: 4,\n",
       "                          18: 1,\n",
       "                          23: 1,\n",
       "                          25: 1,\n",
       "                          13: 1,\n",
       "                          22: 1,\n",
       "                          7: 10,\n",
       "                          21: 1,\n",
       "                          15: 1,\n",
       "                          10: 4,\n",
       "                          19: 1,\n",
       "                          30: 1}),\n",
       "             'diffBUM_HMM_in_vivo': defaultdict(int,\n",
       "                         {2: 13, 1: 80, 5: 4, 7: 2, 4: 2, 6: 5, 3: 2}),\n",
       "             'deltaSHAPE_rep1_ex_vivo': defaultdict(int,\n",
       "                         {3: 70,\n",
       "                          5: 13,\n",
       "                          1: 12,\n",
       "                          2: 9,\n",
       "                          4: 8,\n",
       "                          6: 4,\n",
       "                          8: 1,\n",
       "                          27: 1,\n",
       "                          21: 1,\n",
       "                          19: 1,\n",
       "                          11: 3,\n",
       "                          7: 3,\n",
       "                          13: 1,\n",
       "                          9: 1,\n",
       "                          15: 1}),\n",
       "             'deltaSHAPE_rep1_in_vivo': defaultdict(int,\n",
       "                         {4: 8,\n",
       "                          3: 34,\n",
       "                          2: 13,\n",
       "                          7: 4,\n",
       "                          5: 8,\n",
       "                          1: 8,\n",
       "                          8: 2,\n",
       "                          6: 5,\n",
       "                          13: 1,\n",
       "                          10: 3,\n",
       "                          20: 1,\n",
       "                          9: 1,\n",
       "                          12: 1,\n",
       "                          15: 1,\n",
       "                          11: 1}),\n",
       "             'deltaSHAPE_rep2_ex_vivo': defaultdict(int,\n",
       "                         {1: 15,\n",
       "                          3: 83,\n",
       "                          8: 3,\n",
       "                          2: 17,\n",
       "                          4: 12,\n",
       "                          6: 5,\n",
       "                          11: 1,\n",
       "                          5: 5,\n",
       "                          9: 2,\n",
       "                          33: 1,\n",
       "                          10: 2,\n",
       "                          7: 2,\n",
       "                          16: 1,\n",
       "                          17: 1,\n",
       "                          15: 1,\n",
       "                          12: 1}),\n",
       "             'deltaSHAPE_rep2_in_vivo': defaultdict(int,\n",
       "                         {7: 5,\n",
       "                          3: 56,\n",
       "                          2: 14,\n",
       "                          1: 13,\n",
       "                          4: 9,\n",
       "                          5: 8,\n",
       "                          16: 1,\n",
       "                          10: 2,\n",
       "                          6: 5,\n",
       "                          9: 5,\n",
       "                          8: 1,\n",
       "                          15: 1,\n",
       "                          23: 1,\n",
       "                          13: 1,\n",
       "                          14: 1})})"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intervallengthcounter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This plot was generated to check on average how many DRNs were found in the 30bp sequences used for MEME analyses. \n",
    "As you can see, the majority of the DRN intervals still consist of single-nucleotide DRNs, which means that most DRNs are actually quite far apart in the data. This is not the case for deltaSHAPE that consitenly picks up stretches of 3 DRNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAALICAYAAAAOvn6rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xu8ndO97/HPd0VCRVzapC5tZMV2V1vK2sQlbnXZnLhUaXZVKlTzoqi6bzZtSrVayjkt6iQqqeA0W4oSJSEnISjpysUtKnWJoNUuTpFsZEv8zh9jTHkys+aalzHXvKz1e79e67Wey3jG85uX3xzjeeZ8xiMzwzlXmZZ6B+BcM/MEci6BJ5BzCTyBnEvgCeRcAk8g5xJ4AjmXoFcnkKRjJc2WdJmkb8RlYyQtlTRd0mBJz0paKGlPSavi9EJJz0h6UtLemfpM0sDO9pGp2yT9IK+MJL0s6dkSYp4k6bxOlpukgZJa4/TDBbYtq1yxeLqbpDZJU+sdRyG9OoFyzOx7ZnZLnP0GcLGZHQocALxpZsOAvwIfmNmw+Lcz8DNgUpm7WwqckLdsBLB+xQ9gbR8C20kaklsgqT+wd4Xl6sbM2s3s2HrHUUivS6DY2rwkaS7w5bhskqTzJF0L7A5cLuls4IfAv0iaVaC6zxASqxzPAMsl7ZVZdiJwa5n1dGUVMAX4embZMcDvKizXKUn9JF0rab6kp+LzuKGkTSX9VdLhsdzlkh6UVPD9Jul2Sedm5k+TNEXS/rEXsJGk9yRtlinzpKTD4rpbY7lnJP1U0jqlPIZUvSqBJB0FfAUYBuwFbJRdb2ZnA+3A+WZ2LfA9YI6ZHRCLfCrThXsV+F/AjysI5RZgdIxpfUIL9EAZ25+diWOhpIVd7SM6kc5by1LLdebfgZXAbma2C/AX4Eoz+1usZ7ykL8fp483s4y7qmgCMycyPicsAMLN3gbuIrbekHYDNgOnAz4G3gZ2BNmAXYK1ubnfoVQkEHATcaWbLzGwlcHOZ22e7cEOAw4E7JA2N6zv7YWEL4ZM+6zbgGEn9CK3gPYQ3YqmuzcQxLHYx12Bm84BVknaTNBgYYGZrHWOVWq6AkcBRwIKYxEcDO8Z6ZxBat98CXzezjiJ1zQbWi8c8OwKDgJl5ZW4iJCPAScDNMSkPA66zYAVwY1zW7WrSzDUYZabLedOuxcwekvQiodv3CvAWoVv3VqbYpoRPx+x2b0qaT3iRTwTOAbrjgH0y4RO7I06nlsvXBzjLzO4HkLQBsF6cFrAD8DdgODCnq4rMzCT9inAMugL4VVyWLTNH0jqSdgeOB/aMq1pY88OrBehbxuOoWG9rge4HjpO0ceyPjy62QVckbQu0Agsy9X8n19eXtAkhQX7fyea3AOcCG5XxiV+uW4HjgFHA7VUol286cEY8FmohdLlyXdqzgQ0IXapzJP1LCfVNAo6MsUwsUOYm4BfA02b2Wl4ckrQuMBZ4sIzHUbFelUBm9ntCt60deBJ4t8wqssdAC4GpwFgzWxzXn0X4BH5W0tPAI4RuzK87qetuQl+9nE/8spjZG8DzwJ/N7P+lluvE5cASwgfIIkLrfq6kLwIXAyfGur8L/B9JA4rE+yYwn5AcfylQ7NeEY9ibMsu+A3yWcILmGeAF4IoyHkfF5NcDOVe53ngM1NAkfR04v8Dq28zsqhrGsh2hBe3MC2Y2qsz6pgDbFVg9ysxeKKe+RuAtkHMJetUxkHPV1nRduIEDB1pra2u9w3A93Lx5894ys0HFynVLAknqQziluR3hS8STCGdoJhHO1z8LnG5mH0v6PvA/CN/JfNfM5nZVd2trK+3t7d0RtnOfiL80Kaq7WqAjAMxsb0n7A9cQEugSM5st6UbgqBjkfsAewGDCt9alfF/gXEPolmMgM7ub8GUWwBDCt9G7Abmfzt9P+FnNPsCM+BOMpcA6koo2m841im47iWBmKyX9mvCt8VTCGb/cKb9lhB9ybsiaX2bmlq9B0lhJ7ZLaOzqK/aTKudrp1pMIZnaipAsJ3/p/KrNqAPAO8F6czl+eX894YDxAW1ubDR26pMv9vvJKa0rYzpWsW1ogSaMlXRRn3wc+Btrj8RCEH1HOAR4DDpXUImlLoMXM3lqrQucaVHe1QHcCEyU9QvhV7HcJv7WaEH/C/zww1cxWSZoD/IGQzKd3UzzOdYtuSSAz+y/gq52s2q+TsuOAcd0Rh3PdzX+J4FwCTyDnEngCOZfAE8i5BJ5AziXwBHIugSeQcwk8gZxL4AnkXAJPIOcSeAI5l8ATyLkEnkDOJfAEci6BJ5BzCTyBnEvgCeRcAk8g5xJ4AjmXwBPIuQRVTyBJfSVNljRH0lxJR0raWtKjcdkvM7dA/H4s83i876VzTaWkUXkk9Qc2AT4iDNl7i5kVGnz7BOBtMxst6TOE2/8txMfFdj1QqS3QbYSxra8iJNH4LsreAVyamV9J4rjYPrSva1SlJtAmwD3A58zsSmDdQgXNbLmZLYs3lJ0KXELCuNixzvFm1mZmbYMG+djzrnGUmkD9CLdkny9pR8LtywuSNBiYBUw2s9sJQ/vmlDUutnONrNQEOpdwG/ErgAOA0woVlLQpMAO40MxujosX+LjYricqdWjfA83sgjh9vaQfA38sUPZiQpfvUkm5Y6GzgJ/7uNiup+nyLt2SvgmcAuwALIqLW4B+ZrZr94e3tra2Nnv77aldlvHbm7hUkuaZWVuxcsVaoFuBmYRW5Yq47GPg72nhOdczdHkMZGYrzGwJcCqwKeF2jUMJ39041+uVegw0lXAS4bU4b8Aj3RKRc02k1ATazMz26tZInGtCpZ7G/pOkLbo1EueaUKkt0AhgqaTc72jMzDyhXK9XUgKZ2TbdHYhzzajUX2NPJJw4+ISZndwtETnXRErtwv0m/hewK+DdN+covQs3PTP7gKQZ3RSPc02l1C7cIZnZzQlfqjrX65XahftaZvpDwI9/nKP0LtxJkr4A7AgsNrOF3RuWc82hpC9SJZ0JTAD2AsZLOq9bo3KuSZTahTseGGFmKyX1BR4Hru6+sJxrDqX+lEdmthLAzD4iDCziXK9Xagv0qKSphEux9yFcju1cr1e0BZI0FrgImEgYNedhMzu/uwNzrhl0mUCSxgGHAH3N7D7gFuDAzFgHzvVqxVqgw4DjzOx9gHh16ijgyGIVS9pD0uw47UP7uh6pWAItt7xRR+JJhGVdbSTpAuAmYL246BrC0L4jCL+nO0rSrqwe2vffgOvLD9+5+iqWQB9I2iq7IM4XHsoneAk4JjOfNLSvc42q2Fm4C4G7Jc0EXga2BA4FTuxqIzP7raTWzKJCQ/u+nSmTW77W4NfxRMZYgC233JIWvymLaxDFRuV5jnA16gKgPzAf2NvMFpS5n6ShfX1sbNeoin4PZGbvEs6+pVggaX8zm004MTELeBH4qaSrgc/jQ/u6JlTqF6mpzgUm+NC+rqfptgSKp7yHx+nFhDNu+WXGAeO6KwbnupsfjjuXwBPIuQS1Ogaqi6FDl3S53u/i4FJ5C+RcAk8g5xJ4AjmXwBPIuQSeQM4l8ARyLoEnkHMJPIGcS9Cjv0gtVbEvXMG/dHWd8wQqk/+6wWV5F865BJ5AziXwBHIugR8DdZNSj5X8BEZz8wTqgfxER+14F865BN4C9WLeUqWrewLFcbJvAHYBVgCnmNmL9Y2qMfkbvvHUPYGAo4H1zGxPScOBnwFH1Tkml1HOiY5qnTyp5MOgHh8wjZBA+wAPAJjZE5La6hyPayDdcZaymonWCAm0IfBuZn6VpHVyt5SENcfGBpbD0Bfy6hgIvLW6fMF9lVquO+qsqFw9993Ln58hBUtnmVld/wi3PvlqZv71Cupor2a57qizJ8XYW5+fzv4a4TT2Y8DhAPEY6Jn6huNc6RqhC3cXcLCkxwk33zqpzvE4V7K6J5CZfQycmljN+CqX6446e1KM9dx3PWNci2If0DlXgUY4BnKuaXkCOZfAE8i5BJ5AziXwBHIugSeQcwk8gZxL4AnkXAJPIOcSeAI5l8ATyLkEnkDOJej1CSTpWEmzJV0m6Rtx2RhJSyVNlzRY0rOSFkraU9KqOL1Q0jOSnpS0d6Y+kzSws31k6jZJP8grI0kvS3q2hJgnSTqvk+UmaaCk1jj9cIFtyypXJJZPnrdqknSTpIOqXW+11f1yhkZhZt/LzH4DuNjMbo1vjjfN7CBJrcAHZjYsV1DSV4FJwDZl7G4pcALw/cyyEcD6wPsVPYC1fQhsJ2mImb0KIKk/sHeF5TqV97xVjZmd0h31VluvbIHip+ZLkuYCX47LJkk6T9K1wO7A5ZLOBn4I/IukWQWq+wzw1zJDeAZYLmmvzLITgVvLrKcrq4ApwNczy44BfldhuU5lW0NJH0oaJ+lxSa9IOq3ItttKektSvzjfR9IbknaIvYJjJf1I0i8y2xwm6ck4fbSkBZKekvSopN1Libmael0CSToK+AowDNgL2Ci73szOBtqB883sWuB7wBwzOyAW+VSmC/cq8L+AH1cQyi3A6BjT+oQW6IEytj87E8dCSQu72kd0IqG1rLRcMesCb5nZXsCxwLWS1itU2MwWA88BR8ZFhwCvmNnzmWI3Af+WSzJgDDBB0vbAjcBXzGwXwuv0O0kbVhB3xXpdAgEHAXea2TILI//cXOb2H5jZsPg3hDCewx2Shsb1nV2h2EL4pM+6DTgmvjG+DNwDrMzfsAvXZuIYlu1W5pjZPMIoR7tJGgwMMLO1jrFKLVeiXMs1n5BQ/YuUv4mQFBAu55+QF9vLwNPAkZI2AQ4EfhP/z4zrMbP/C/wd2K3CuCvSW4+BsgMclfOmXYuZPSTpRUK37xXC8EifYc2hlzYF3s7b7k1J84HDCJ/45xCGV6q2yYTjrY44nVqumA8AzMwUxofqenAsuAO4RtIOwH6sTqasCYTj0k2Bu81suaQ+rP1h1QL0rTz08vXGFuh+4DhJG8dhhUcX26ArkrYFWoEFmfq/E+smfmqeCPy+k81vAc4FNkr4xC/mVuA4YBRwexXKVZWZfUhoUSYBvzWzzk6i3EVoWb7F6hZqJnCopK0AJB0IDAae7O6Ys3pdC2Rmv5e0M+E45x/AU8CgMqr4VN7xRgswNvbnAc4iDE/8rKSVhE/gW4Bfd1LX3YR+/H+U9yhKZ2ZvSHoeeNfM/l9quW4yATgD6PSkg5mtkDQFOMjM5sZliyR9G7hT0jqEs5dHmNm7ndXRXXxQEecS9LoWqBlI+jpwfoHVt5nZVTWMZTvCae7OvGBmo4psfwBwbYHVs+JZz6blLZBzCXrjSQTnqqbpunADBw601tbWeofherh58+a9ZWZFTy41XQK1trbS3t5e7zBcDxd/ZVKUd+GcS+AJ5FyCpuvCgd9s1zUOb4GcS1CzBIpXYs6Of0/Ea0eOidfl5JbvV6t4nKuGmnXhzGwS8RoTSdcTLiPYFbjAzH5bqzicq6aad+EUbmO/k5mNJ/zC9mRJcyT9LP4o0LmmUY9joIuB3IAaDwJnAvsCG1DgVo+Sxkpql9Te0dFRmyidK0FNE0jSxsD2ZpYbX+BmM3vZwg/yfgd8sbPtzGy8mbWZWdugQeVceeBc96p1C7Qv8BCEYZyApyV9Pq77EjCvxvE4l6TWxxzbAblr2E3SKYQLoj4AFpF3Pbxzja6mCZR/HYuZzQBm1DIG56rJv0h1LoEnkHMJPIGcS+AJ5FwCTyDnEngCOZfAE8i5BJ5AziXwBHIugSeQcwk8gZxL4AnkXAJPIOcSeAI5l8ATyLkEnkDOJag4gST1l/R5SZtKulTSkGoG5lwzSGmBbiMMS3UV8BEwvioROddEUhJoE+Ae4HNmdiWwbrENJC3IjEI6UdJwSU9KekzS9xNica4uUsZE6Ee4Rft8STsSxnUrSNJ6AGa2f2bZQuArhIFG7pO0q5nNT4jJuZpKaYHOBT4LXAEcQIFblGfsAqwvaYak/ytpX2BdM3spjgs3nTC0lXNNIyWBDjSzC8zsHTO7HjimSPn3gauBQwkjkE6My3KWARt1tqGPTOoaVdldOEnfBE4BdpB0eFzcQujSXdTFpouBF2Nrs1jSu8CnM+sHAO90tmEcR3s8QFtbm739drlRO9c9KjkGuhWYSRjj+oq47GPg70W2OxnYGfi2pC2A9YH/kvRPhGOgQ1k9ZrZzTaHsBDKzFcASSacCbcB6cdVQ4JEuNv0VMEnSo4AREupjwunwPsAMM3uy3Hicq6eUs3BTCScRXovzRhcJZGb/DRzfyarhCTE4V1cpCbSZme1VtUica0IpZ+H+FI9lnOu1UlqgEcBSSbnzymZmnlCuV6k4gcxsm2oG4lwzqjiBJE0knDj4hJmdnByRc00kpQv3m/hfhLtte/fN9TopXbjpmdkHJPmNslyvk9KFOyQzuzmwaXo4zjWXlC7c1zLTHxJ+WeBcr5LShTtJ0heAHYHFZrawemE51xxSxkQ4k3BX7b2A8ZLOq1pUzjWJlC7c8cAIM1spqS/wOOF6H+d6jZSf8sjMVgKY2UeEgUWc61VSWqBHJU0F5gD7AI9VJyTnmkdFLZCksYSrTycSLsN+2MzOr2ZgzjWDshNI0jjgEKCvmd0H3AIcKOnSKsfmXMOrpAU6DDjOzN4HMLMlwCjgyCrG5VxTqCSBlseBQT4RTyIsq05IzjWPShLoA0lbZRfEeStQPlemr6TJkuZImivpSEm7SnojM1rpqArica5uKjkLdyFwt6SZhNF0tiSMqHNike1OAN42s9GSPgMsAC4DrjGzn1UQh3N1V3YLZGbPEa5GXQD0B+YDe5vZgiKb3gFkTzSsJAxO/z8kPSLpV5IGdLahD6zoGpXyDme6f4chSe4h/AxoXeBpM5sn6T+ATcysy58EhYEVp3a5j1deaa1StK63kjTPzNqKlUv5IrVskgYDdwE3mNntkjY2s9xopHcBv6hlPM6lqtkd6iRtCswALjSzm+Pi6ZJ2j9NfAubVKh7nqqGWLdDFhHsKXZr50vUc4H9K+m/gTWBsDeNxLlnNEsjMzgLO6mSVD87ompbfZNi5BJ5AziXwBHIugSeQcwk8gZxL4AnkXAJPIOcS1PSnPLU2dOiSLtf7b+ZcKm+BnEvgCeRcAk8g5xJ4AjmXwBPIuQSeQM4l8ARyLkGP/h6oGRT7rgr8+6pG5gnUTfxL3N6h7l04SS2SbpT0hzi44tb1jsm5UjVCC3Q0sJ6Z7SlpOPAz4KhaBlBON6qeLYu3ao2n7i0Q4d5CDwCY2RNA0bG4nGsUjdACbQi8m5lfJWmd3N3v4JP7EeVG7FkOQ1/Iq2Mg8Nbq8gX3VWq57qizonLdVWdiue6os9FiHFLSHsysrn/ANcBXM/OvV1BHezXLdUedPSnG3vr8dPbXCF24x4DDAeIx0DP1Dce50jVCF+4u4GBJjwMCTqpzPM6VrO4JZGYfA6cmVjO+yuW6o86eFGM9913PGNdS87szONeTNMIxkHNNyxPIuQSeQM4l8ARyLoEnkHMJPIGcS+AJ5FwCTyDnEngCOZfAE8i5BJ5AziXokQkk6VhJs4uUWSKpLU7PkDSwhHqHS5ol6WlJz0q6X9JOmfWWX4+kMZKm5S07IpYdlbd8nKQOSQslLYj7eUjStp2sz/5dWfRJqSJJu0j6ay332RVJW8Rf89dc3X+N3SAOLlZA0rrANOAQM5sfl50A3C9pqJmtKmN/3wZuA84GpuStm2JmZ2T2eyZwO6svdV9jfS1JWgc4E7gQ2KAeMXTGzP4C7FWPffeYFkjSZZJekjQX+HJc1k/StZLmS3pK0iRJG+ZtNzFOzpI0WNJISY9Lape0VNLlcf36wMas+ca5DTgD6FNGnFsB+wPnANvEiwi7MhPYvtT6C+yzVdJrsaVdLGlzSXtJmhOfmz9KGhnLjpF0X2xdF0l6UNIWsapdgZ2Jz28Z+/+PuJ+Fku6OLcanJD0n6duxzDclPSNp/S7q+ZGkX2TmD5P0ZHx8yyX1iY9zt0yZKZJOk9RX0i/iY3pG0k2SBpTzODpV6aWsjfRHGMXnOWAAoVWdBswGvgdcxerLNn4E3BCnlwBtcdoI18ULmAVsE5dvAawEBsb5c4D3gZeBycDJwPqZOIxwRe3CzN9SYFqmzE+BqXH6ekKLkls3DrguM78OYZSiezPrO/LqXwgcWuT5aY2xjYjzmwAvAK2Zx/kasCUwBlgObBvXXZmLN6++5SW+Nt8AfgOsE+fHAr+P0zsTxiI4HHgT2K5IXVvFx98vzk8BTsnGA/wg9xzGx/k2sFFc/lugL6HhuBm4Mfm9V+83f5US6BfA5Zn5Y2ICzQUWZ95oi4DZhRIoTm8A/Bvw/fjCrwKGZOoeAIyMifBC/Nsov55M+TG5BALWjW+AI+L87sBHwOACCbKIkKibdZZgZTw/rXE/uTfx4cB7rJ3oI2O8d2a2/SdgWUIC/Wd8rnP7eQZ4KbP+TOBjYHSJ9c0Ejo3J0RFfr2wCDQH+DvQjdJUnxuVzgYMz9XwReDX1vddjunCE1iMnN6JPH+AsMxtmZsMIb9hjC1Yg9QcWELoq84HzCW88Sdpb0vlmtszMppnZBcBOhKQpegwVfZXwwl8naQnhzWWEN1HOlFy8ZrajmY02szdLrL8rK2z1SEd9gOcz+xkGDAemx/UrM9u1ED5EKtUH+ElmP23A3pn1OwF/i/svxQRCq3Y8cLeZLc+uNLNXCa/dSMLwADdl4shePdpCaI2S9JQEuh84TtLGklqA0XH5dOCMeCzUQnjyf9zJ9qsIT+Y2hGG2LjGzewnHKusSnvwO4BJJ+2S225zQPSh1IJTTgCvMbIiZtZpZK+Fy9m/F5K2VJwjHX/sCSBoG/Bn4XFz/JUm56VOBexP2NR04JXPseRmhVUXSMcCBwD8Dh0g6uoT67gJ2A75FeD07M4FwoqO/mT0Wlz0A5I6FWoDTgQcreDxr6BFn4czs95J2BtqBfwBPAYOAy4GrCa1KH0IX4txOqrgDeJjQOk0D/iRpBSExFgFbm9n0+AL/SNLngQ8J49mdZGb549StRdIuwDDgyLxVtwCXELpOpRiVl8QAS80sv96CzKxD0leAqyStR/ggHW1mSxQGl3sdmCxpc8LjH1u4tqJuIiTmE5KM0FUcI2kwcCOhO9sh6UTgbkntZvZ6F7GvkDQFOMjM5hYodg9wA/CTzLIfEt4LCwnv+7ms2fJXxMdEcGuQNAY41sxG1juWZtAjWiAHkq4FDiiw+mwzm9VN+x0AzCmwepmZjSizvro8jkp5C+Rcgp5yEsG5uvAEci5B0x0DDRw40FpbW+sdhuvh5s2b95aZDSpWrukSqLW1lfb29nqH4Xo4Sa+WUs67cM4l8ARyLkHTdeG6g99q3lXKWyDnEngCOZfAE8i5BJ5AziXwBHIugSeQcwk8gZxL4AnkXIKaJpCkiyT9QdK8OA7Y1pIejeOT/TJeq+5c06jZG1bS/oTRI/cG9gMGA9cQBvAYQRhV56haxeNcNdTyE/9QwiAddxFGeZlGGF3l4bj+fuCgGsbjXLJa/hZuIGHQu5HAUMLIKS22+pryZYQhotYiaSxxZJgtt9yy+yN1rkS1bIHeBqab2X/HYaA+ZM2EGQC809mGZjbezNrMrG3QoKLXODlXM7VMoEeBf1WwBdAfmBmPjQAOo/DoLs41pJp14cxsWhwJcy4hcU8HXgEmSOoHPA9MrVU8zlVDTa8HiuNJ59uvljE4V03+vYtzCTyBnEvgCeRcAk8g5xJ4AjmXwBPIuQSeQM4l8ARyLoEnkHMJKk4gSf0lfV7SppIulTSkmoE51wxSWqDbCNfzXEW4Ffz4qkTkXBNJSaBNCNf0fM7MriTcDt65XiUlgfoRbhk/X9KOwAbVCcm55pGSQOcCnwWuINxV+bSqRORcE0lJoAPN7AIze8fMrgeOqVZQzjWLsq8HkvRN4BRgB0mHx8UthC7dRVWMzbmGV8kFdbcCM4GLCd03gI+Bv5eysaTPAvOAg4GVwCTAgGeB083s4wpicq4uyu7CmdkKM1sCnApsShhpZyiwR7FtJfUF/jfwQVzk48K5ppZySfdUwkmE1+K8AY8U2eZq4EZWd/Xyx4U7hDBunHNNISWBNjOzvUotLGkM0GFm0yXlEkg+LpxrZiln4f4Uh6cq1cnAwZJmA8OAWwgtWI6PC+eaTkoLNAJYKqkjzpuZFUwoM9s3Nx2T6FTgKkn7m9lswrhwsxLica7mKk4gM9umCvs/Fx8XzjWxihNI0kTCiYNPmNnJpWxrZvtnZn1cONe0Urpwv4n/BewKlHM85FyPkNKFm56ZfUDSjCrE41xTSenCHZKZ3ZzwpapzvUpKF+5rmekPCaepnetVUrpwJ0n6ArAjsNjMFlYvLOeaQ8qYCGcCEwj3PR0v6byqReVck0jpwh0PjDCzlfFHoo8TfuvmXK+R8lMemdlKADP7iDCwiHO9SkoL9KikqYTbMu4DPFadkJxrHhW1QPHX0RcBEwm/oH7YzM6vZmDONYOyE0jSOMJ1O33N7D7Cr6oPlHRplWNzruFV0gIdBhxnZu8DxKtTRwFHVjEu55pCJQm0PHMRHPDJSYRl1QnJueZRSQJ9IGmr7II4bwXKO9djVXIW7kLgbkkzgZeBLYFDgROrGZhzzaCSUXmeI1yNugDoD8wH9jazBV1tJ6mvpMmS5kiaK+lISVtLejQu+6Ukv92KayoVfQ9kZu8Szr6V4wTgbTMbLekzhARcSBjWarakGwnDWvmoPK5p1PIT/w4ge6p7JWsPa3VQDeNxLlnNEsjMlpvZMkkDCGMfXEIZw1pJapfU3tHR0VkR5+qipscckgYTRt6ZbGa3E4YEzvFhrVzTqVkCSdoUmAFcaGY3x8ULJO0fpw8j/K7OuaaR8mPScl1MuKvdpZmf/ZwF/NyHtXLNqmYJZGZnERImnw9r5ZqWf+/iXAJPIOcSeAI5l6CWJxF6laFDl3S5/pVXWmsSh+teTZlA/uZ0jcK7cM4l8ARyLoEnkHMJmvIYqFR+rOS6W49OoO7gSemyvAvnXAJPIOcSeAI5l8ATyLkEfhKhzoqdlAA/MdHIPIGaSL3OAHqSF1b3LpykFkk3SvqDpNmStq53TM6VqhFaoKOB9cxsT0nDgZ8RxodzFSq1peofeTNBAAAgAElEQVSOFq3a+y6n9avH426EBNoHeADAzJ6Q1JZfIN6PaGycXQ5DX8grMhB4a3X5gvsqtVx31FlRuXruu5c/P0MKls4ys7r+ATcBh2XmlwLrlFlHezXLdUedPSnG3vr8dPZX92Mg4D3CmHA5LRbvvepco2uEBHoMOBwgHgM9U99wnCtdIxwD3QUcLOlxQMBJFdQxvsrluqPOnhRjPfddzxjXotgHdM5VoBG6cM41LU8g5xJ4AjmXwBPIuQSeQM4l8ARyLoEnkHMJPIGcS+AJ5FwCTyDnEngCOZfAE8i5BD06gSQdK2l2kTJLclfBSpohaWAJ9Q6XNEvS05KelXS/pJ0y6y2/HkljJE3LW3ZELDsqb/k4SR2SFkpaEPfzkKRtO1mf/buy6JNSRZJ2kfTXEsueKunfuyGGyyR9o9r1lqoRLmdoJAcXKyBpXWAacIiZzY/LTgDulzTUzFaVsb9vA7cBZwNT8tZNMbMzMvs9E7gdaOtsfS1JWgc4E7gQ2KCUbczsxu6Ixcy+1x31lqrHtUDxE+klSXOBL8dl/SRdK2m+pKckTZK0Yd52E+PkLEmDJY2U9LikdklLJV0e168PbMyab5zbgDOAPmXEuRWwP3AOsE28mLArM4HtS62/wD5bJb0WW9rFkjaXtJekOfG5+aOkkbHsGEn3xdZ1kaQHJW0Rq9oV2Jn4/Ja473GSrovTS+L8HEmvZp7bQttuJOk9SZtllj0p6bD4Wp4naaykezPrt5f0hqQ+kkZIeiK25O2S/rWMp61rlV4L3oh/hNF8niNcIr4OoaWYDXwPuIrV1z/9CLghTi8B2uK0EQaYEDAL2CYu3wJYCQyM8+cA7wMvA5OBk4H1M3EY4crahZm/pcC0TJmfAlPj9PWEFiW3bhxwXWZ+HcJoRfdm1nfk1b8QOLTI89MaYxsR5zcBXgBaM4/zNWBLYAywHNg2rrsyF29efctLfG0+eUzxOb86Tn8O+AAYWmT7XwPnxekdgFcJDcAk4Lz4mv8D2CyW+QnwQ+AzwN+APeLynQgDiHS5v1L/eloX7iDgTjNbBiDpZuA7wEhCq3GwwpAr/YC/F6rEzEzSEcBISccTXjAB/YG3zOwaSROA/YB9CV2ZCyXtbmbvxmoOMLPMSC8aAxwbp9clXHl7clz9a+AxSYPN7LW4bJSkfeJ0P2Ae8K1MmJV24VYCf4jTewKbA3dr9ZA1BvxznJ5hZovj9ARCklbL7wDM7A1Jfwc+DbzSRfmbgBuAqwnP3c1m9nEubjNbJulO4ARJ1wJfB0YAewAvmtmTsdxzkh4jtP5d7a8kPa4LR3ij5+QGJ+kDnGVmw8xsGLA78c3caQVSf2ABoasyHzgf+Cis0t6SzjezZWY2zcwuIHyqGSUcQ0VfJXz6XydpCfCfcfszM2Wm5OI1sx3NbLSZvVli/V1ZYasHbekDPJ/ZzzBgODA9rs8O7tIClHN8V8wHmWljzddtLWY2B1hH0u7A8cDNnRSbAHwD+FfC43qF8BjzL7tuAfpWGPdaFfUk9wPHSdpYUgswOi6fDpwRj4VaCE/0jzvZfhXhid0G2BC4xMzuJXxarUt4MTqASzKtA4RP8Y0ofUCU04ArzGyImbWaWStwKvCtmLy18gTh+GtfAEnDgD8TulUAX5KUmz4VuHftKmrqJuAXwNOZlvoTZvYEIRG/R3iNIbS228fEI54t3ZfQtU/Wo7pwZvZ7STsD7YT+8FPAIOByQtO/gJAEC4FzO6niDuBhQus0DfiTpBWExFgEbG1m0yUdDfxI0ueBD4F3gZPMLH/Ax7VI2gUYBhyZt+oW4BLCsUcpRuUlMcBSM8uvtyAz65D0FeAqSesRPlBHm9mS2DV6HZgsaXPC4x9buLaa+DXh+PVrXZSZAFwK3A1gZm9JOg74haT1gY8Jr9XiLuoomQ8q4jqVO2Yzs5H1jqWR9agWyEE8gD6gwOqzzWxWN+13ADCnwOplZjaiyPbbsfZ3YTkvmNmoAuvqylsg5xL0tJMIztWUJ5BzCZruGGjgwIHW2tpa7zBcDzdv3ry3zGxQsXJNl0Ctra20t7fXOwzXw0l6tZRy3oVzLoEnkHMJmq4L1x38LtSuUt4COZfAE8i5BJ5AziXwBHIugSeQcwk8gZxLULPT2JL6EC522o5w5edJhKsHJxEuuX0WON3MPq5VTM6lqmULdASAme1NuOT2mvh3SbxWRIRRdZxrGjVLIDO7m9WXBA8hDDW0G+ESagjjGRxUq3icq4aaHgOZ2UpJvyYMDDGVcEFf7oq+ZYSBOdYSB81rl9Te0dFRo2idK67mJxHM7ERgW8Lx0KcyqwYA7xTYZryZtZlZ26BBRX9h7lzN1CyBJI2WdFGcfZ8wOkq7pP3jssMofE29cw2plj8mvROYKOkRwthr3wWeByZI6henp9YwHueS1SyBzOy/CCNy5tuvVjE4V23+RapzCTyBnEvgCeRcAk8g5xJ4AjmXwBPIuQSeQM4l8ARyLoEnkHMJPIGcS+AJ5FyCihNIUn9Jn5e0qaRLJQ2pZmDONYOUFug2whWlVxFuAT++KhE510RSEmgT4B7gc2Z2JeE28M71KikJ1I9wq/j5knYENqhOSM41j5QEOhf4LHAF4a7Qp1UlIueaSEoCHWhmF5jZO2Z2PXBMV4Ul9ZU0WdIcSXMlHSlpa0mPxmW/lORnBV1TKfuKVEnfBE4BdpB0eFzcQujSXVRwQzgBeNvMRkv6DLAAWEgYF262pBsJ48LdVW5MztVLJZd03wrMBC4mdN8gDBDy9yLb3cGaYx6sZO1x4Q7BE8g1kbK7TGa2wsyWAKcCmxIGSRwK7FFku+VmtkzSAEIiXYKPC+eaXMqgIlMJJxFei/MGPNLVBpIGE1qYG8zsdkk/zazuclw44vdMbW1t1lkZ5+ohJYE2M7O9Si0saVNgBnCGmc2MixdI2t/MZhPGhZuVEI9zNZeSQH+StIWZ/aXE8hcTvny9VNKlcdlZwM99XDjXrFISaASwVFLuoMTMbItChc3sLELC5PNx4VzTqjiBzGybagbiXDOqOIEkTSScOPiEmZ2cHJFzTSSlC/eb+F/ArkDB7ptzPVVKF256ZvYBSTOqEI9zTSWlC3dIZnZzwpeqzvUqKV24r2WmPwT8+Mf1OilduJMkfQHYEVhsZgurF5ZzzSFlTIQzCbdp3AsYL+m8qkXlXJNI6cIdD4yINw7uCzwOXF2dsJxrDikXsMnMVgKY2UeEgUWc61VSWqBHJU0l3Bh4H+Cx6oTkXPOoqAWSNJZw9elEwjU8D5vZ+dUMzLlmUMkl3eOALwC3mtl9kp4DrpG0iZldXu0AUwwduqTL9a+80lqTOFzPVUkLdBhwnJm9DxCvTh0FHFnFuJxrCpUk0PLMZdjAJycRllUnJOeaRyUJ9IGkrbIL4rxfau16nUrOwl0I3C1pJvAysCVwKHBiKRtL2gP4iZntL2lrYBIh+Z4FTjezjyuIybm6qGRUnucIV6MuAPoD84G9zWxBsW0lXQDcBKwXF11DGBduBOGyiKPKjce5eqroeyAzexe4pYJNXyKMYDo5zvu4cK6p1XQoXTP7LWv+YsHHhXNNrd5jUWePd7ocF87M2sysbdCgQbWJzLkS1DuBFkjaP04fRvhZkHNNI+W3cNVwLjDBx4VzzarmCRR/uTA8Ti/Gx4VzTazeLVBF6vkbt1L37b/D6x3qfQzkXFPzBHIugSeQcwk8gZxL4AnkXAJPIOcSNOVp7J6k2Olu8FPejcxbIOcSeAvURPzL2cbjLZBzCTyBnEvgXThXlJ/oKMxbIOcSeAvUA/nJhtrxFsi5BHVvgSS1ADcAuwArgFPM7MX6RtU7dEdLVe3rpco5/qrHtVqN0AIdDaxnZnsC/w78rM7xOFeyRkigfYAHAMzsCaCtvuE4V7q6d+GADYF3M/OrJK2Tu/sdfHI/orFxdjkMfSGvjoHAW6vLF9xXqeW6o86KytVz3738+RlSsHSWmdX1jzC871cz869XUEd7Nct1R509Kcbe+vx09tcIXbjHgMMBJA0HnqlvOM6VrhG6cHcBB0t6nDDA/El1jse5ktU9gSzczuTUxGrGV7lcd9TZk2Ks577rGeNaFPuAzrkKNMIxkHNNyxPIuQSeQM4l8ARyLoEnkHMJPIGcS+AJ5FwCTyDnEngCOZfAE8i5BJ5AziXokQkk6VhJs4uUWSKpLU7PkDSwhHqHS5ol6WlJz0q6X9JOmfWWX4+kMZKm5S07IpYdlbd8nKQOSQslLYj7eUjStp2sz/5dWfRJqQJJgyTdGeNaJOmqOKZF3cXnYeNa77fuv8ZuEAcXKyBpXWAacIiZzY/LTgDulzTUzFaVsb9vA7cBZwNT8tZNMbMzMvs9E7id1Ze6r7G+xq4FFpnZMZLWA2YAY4Cb6xTPJ8xsWD322xCfHtUg6TJJL0maC3w5Lusn6VpJ8yU9JWmSpA3ztpsYJ2dJGixppKTHJbVLWirp8rh+fWBjYIPM5rcBZwB9yohzK2B/4Bxgm3gRYVdmAtuXWn+BfbZKei22tIslbS5pL0lz4nPzR0kjY9kxku6LresiSQ9K2iJWdRdwHYCZfQg8SwmXPkv6pqR5sVV9SNL2klokzZT001jmIEmvS9q0i3rGSro3M7+9pDck9cm1/vG1+0qmzE8k/SROXxof09OSpkrarPxnM0+ll7I20h9wFPAcMIDQqk4DZgPfA65i9WUbPwJuiNNLgLY4bYTr4gXMAraJy7cAVgID4/w5wPvAy8Bk4GRg/UwcRriidmHmbykwLVPmp8DUOH09oUXJrRsHXJeZX4cwStG9mfUdefUvBA4t8vy0xthGxPlNgBeA1szjfA3YktCiLAe2jeuuzMWbV+cXgXeAYUX2vR/wSO55Ag4Bno/TmwNvxtfvdWDfInUNAP4BbBbnfwL8MO81PCn3fBM+2N4AtonLHwf6Z57LB5Lfe/V+81cpgX4BXJ6ZPyYm0FxgceaNtgiYXSiB4vQGwL8B3wd+A6wChuS9iCNjIrwQ/zbKrydTfkzmBV03JsARcX534CNgcIEEWURI1M0y66+r4PlpjftZJ84fDrzH2ok+MsZ7Z2bbfwKW5dV3aIxzdAn7/inwl7x9vQl8Oq4/AvgYuLTEx/Ir4LyYHK8DQ/MSaH3gbWCz+DhnxfX/CXwrU88m8Tnpl/Le60nHQNkxV3Ij+vQBzjKz+wEkbQCsV7ACqT+wgNBVmUPo2x8dVmlvYC8zu4rQwk2TdDGhG3MwMLWEGL9KeOGuk/SLuMyAM4EL4nx3HeOssNUjHfUhtAJ75FbGbloH8HVWP38QuvmrMuXOIYzf9zUze6iE/fYBJpvZhXH7FkKL94+4fifgb8AenW++lgmEK0ifj4/hlexKM3tf0h3A8cCewE2ZOLJXj7YQWviuxx4qoqccA90PHCdp4/gCjY7LpwNnxGOhFsKT/+NOtl8F9CU09RsCl5jZvYRjlXUJT34HcImkfTLbbQ5sROkDoZwGXGFmQ8ys1cxaCZezfysmb608QTj+2hdA0jDgz8Dn4vovScpNnwrcG8udDpwODC8xeSC8Bl+TtHmmvpmxvt2BswgnSDaSdFaxyiyMHShC93xCgWITgBOBvYHfxmUPACdnnufvAI+Y2YoSH0enekQLZGa/l7Qz0E74ZHsKGARcDlxNaFX6ELoP53ZSxR3Aw8CxhNblT5JWEBJjEbC1mU2XdDTwI0mfBz4kjGd3kpnlj1O3Fkm7AMOAI/NW3QJcQug6lWJUXhIDLDWz/HoLMrOOeKB9VTyb1kLoji1RGBTtdWByfNMvAsZK6kc4HnoPuFOrB1m7w8yu6GJfM+JB/IOSPo7bH0PoKv8f4Ewze0PSGGCupEfMbEGRhzABuBS4u8A+50laRTh2+zAu/hUwOO6jBXiR0Nom8TER3BriG/lYMxtZ71iaQY9ogRxIuhY4oMDqs81sVjfuew7h5EpnRpjZsjLqOp/CLcNVZnZbufF1J2+BnEvQU04iOFcXnkDOJWi6Y6CBAwdaa2trvcNwPdy8efPeMrNBxco1XQK1trbS3t5e7zBcDyfp1VLKeRfOuQSeQM4laLouXHco50a2zmV5C+RcAk8g5xJ4AjmXwBPIuQSeQM4lqFkCxcEqZse/JyR9KOmYOBBIbvl+tYrHuWqo2WlsM5sETAKQdD3hculdgQvM7LeFt3SucdW8C6cwmOFOZjYe2I1wme0cST+T1GlCx+GM2iW1d3R01DRe57pSj2Ogi4EfxOkHCQNq7Eu4xLfT292b2XgzazOztkGDiv6+z7maqWkCKQy9un3m6sibzexlC1f1/Y4w1phzTaPWLdC+wEMQxokCno4DdAB8CZhX43icS1Lr38JtRxjVEzMzSacQRnj5gDD6S6FhipxrSDVNoDgoYXZ+BmGAcueakn+R6lwCTyDnEngCOZfAE8i5BJ5AziXwBHIugSeQcwk8gZxL4AnkXAJPIOcSeAI5l8ATyLkEFSeQpP6SPi9pU0mXShpSzcCcawYpLdBthEuyrwI+Itx63LleJSWBNgHuAT5nZlcSbgfvXK+SkkD9CLeMny9pR8KYBs71KikJdC7wWeAKwt2hTyu2gaQFmTHgJkoaLulJSY9J+n5CLM7VRcoVqQea2QVx+npJPwb+WKiwpPUAzGz/zLKFwFcIl3nfJ2lXM5ufEJNzNVV2Akn6JnAKsIOkw+PiFkKX7qIuNt0FWF/SjLjfccC6ZvZSrHc6YWARTyDXNCppgW4FZhLGd7siLvsY+HuR7d4HrgZuArYB7gfeyaxfBmzV2YaSxgJjAbbccssKQnaue5R9DGRmK8xsCWEQxE2BIcBQYI8imy4GbrVgMfAu8OnM+gGsmVDZffrAiq4hpRwDTSWcRHgtzhvwSBflTwZ2Br4taQtgfeC/JP0T4RjoUFaPWOpcU0hJoM3MbK8yyv8KmCTpUUKynUzo+t0G9AFmmNmTCfE4V3MpCfQnSVuY2V9KKWxm/w0c38mq4QkxOFdXKQk0AlgqKXe7BDOzLaoQk3NNo+IEMrNtqhmIc82o4gSSNJFwLPMJMzs5OSLnmkhKF+438b8Id5rz7pvrdVK6cNMzsw/EXxg416ukdOEOycxuTvhS1bleJaUL97XM9IeE73Wc61VSunAnSfoCsCOw2MwWVi8s55pDypgIZxLuKLcXMF7SeVWLyrkmkdKFOx4YYWYrJfUFHif82tq5XiPlilSZ2UoAM/uIMLCIc71KSgv0qKSpwBxgH+Cx6oTkXPOoqAWKF7hdBEwENgIeNrPzqxmYc82g7ASSNA44BOhrZvcBtwAHSrq0yrE51/AqaYEOA44zs/cB4tWpo4AjqxiXc02hkgRabmb5PyL9iDCmgXO9SiUJ9IGkNQb/iPNWoLxzPVYlZ+EuBO6WNJMwlsGWhPEMTuxqo/hd0c1AK2EY4B8CrwP3An+OxX5pZlMqiMm5uig7gczsOUkjgKMIlzDMBy4zs2JduBOAt81stKTPAAuAy4BrzOxn5cbhXCOo6HsgM3uXcPatHHcQRvLJWUm4u8N2ko4itELfLSERnWsYNbvBlpktN7NlkgYQEukSYC5wvpntS+gOdjo+tqSxktoltXd0dHRWxLm6qOkd6iQNBmYBk83sduAuM5sXV98FfLGz7XxgRdeoapZAkjYFZgAXmtnNcfF0SbvH6S8B8zrd2LkGlfJbuHJdTLgp16WZXy2cA/xPSf8NvEkc/9q5ZlGzBDKzs4CzOllVzuimzjUUv0u3cwk8gZxL4AnkXIJankSouaFDl3S5/pVXWmsSh+u5vAVyLkFTtkDesrhG4S2QcwmasgWqJ2/9XJa3QM4l8ARyLoEnkHMJPIGcS+AJ5FwCTyDnEvhp7G5S6unuYuWyZV3j8RbIuQTeAvVA/mVv7dQ9gSS1ADcAuwArgFPM7MX6RtWYqp0Y3s1MV/cEAo4G1jOzPSUNB35GGLTRNaFqJWUlyVvtfZeiERJoH+ABADN7QlJbfoF4P6LcgCPLYegLeUUGAm+tLl9wX6WW6446KypXz3338udnSMHSWWZW1z/gJuCwzPxSYJ0y62ivZrnuqLMnxdhbn5/O/hrhLNx7wIDMfIvFe6861+gaIYEeAw4HiMdAz9Q3HOdK1wjHQHcBB0t6HBBwUgV1jK9yue6osyfFWM991zPGtSj2AZ1zFWiELpxzTcsTyLkEnkDOJfAEci6BJ5BzCTyBnEvgCeRcAk8g5xJ4AjmXwBPIuQSeQM4l6NEJJOlYSbOLlFmSu4hP0gxJA0uod7ikWZKelvSspPsl7ZRZb/n1SBojaVresiNi2VF5y8dJ6pC0UNKCuJ+HJG3byfrs35VFn5QqkDRI0p0xrkWSroqX5ne1zZGSft4NsZwq6d+rXW+pGuHX2I3k4GIFJK0LTAMOMbP5cdkJwP2ShprZqjL2923gNuBsYEreuilmdkZmv2cCtwNtna2vsWuBRWZ2jKT1gBnAGODmQhuY2T3APdUOxMxurHad5ehxLZCkyyS9JGku8OW4rJ+kayXNl/SUpEmSNszbbmKcnCVpsKSRkh6X1C5pqaTL4/r1gY2BDTKb3wacAfQpI86tgP2Bc4Bt4rVQXZkJbF9q/QX22SrptdjSLpa0uaS9JM2Jz80fJY2MZcdIui+2roskPShpi1jVXcB1AGb2IfAsRS6BzrbAkmZL+rGkR2IP4KauWjBJfWLcu2WWTZF0WmyNr5N0iKRnMus3lvQPSZtI2inTY3hK0jcqfQ7XUumlrI34RxiM5DnCFa7rEFqK2cD3gKtYffnGj4Ab4vQSoC1OG+H6eAGzgG3i8i2AlcDAOH8O8D7wMjAZOBlYPxOHES4MXJj5WwpMy5T5KTA1Tl9PaFFy68YB12Xm1yEMtnJvZn1HXv0LgUOLPD+tMbYRcX4T4AWgNfM4XwO2JLQoy4Ft47orc/Hm1flF4B1gWJF9j8k9/via/CfhA3xD4A3ggCLb/yD3nMS43wY2yj1X8TV7JfNangbcGp+7l4BjMo/xdWDParzneloX7iDgTjNbBiDpZuA7wEhCq3GwwogR/YC/F6rEzEzSEcBISccDOxBeoP7AW2Z2jaQJwH7AvsCFwIWSdjezd2M1B5hZZqAKjQGOjdPrEi4cPDmu/jXwmKTBZvZaXDZK0j5xuh8wD/hWJsxKu3ArgT/E6T2BzYG7tXrEDQP+OU7PMLPFcXoCIUk/IelQwpv0TDNbY10J7jWzj4H3JL0IfLpI+ZuBP0o6B/gacI+ZvZuLO75mNxMStZ3w/J4PbEsY9enOWO4vkn4L/Curn4eK9bguHOGNnpMbW6EPcJaZDTOzYcDuxDdzpxVI/YEFwK7AfMIL8VFYpb0lnW9my8xsmpldAOxEeOMVPYaKvkr4FL1O0hLCp7EBZ2bKTMnFa2Y7mtloM3uzxPq7ssJWjznRB3g+s59hwHBgelyfHZuiBfjk+C6+kScDXzOzyRXE8UFm2ljzdVuLmb1KeC1GEpLjpk6K3QwcJ2kYsLGZPUx4jPlXjbYAfSuIeS09LYHuJzyBG8c+9ei4fDpwRjwWaiF8mv64k+1XEZ7YbQhdi0vM7F7Cscq6hBejA7gk0zpA+BTfiNLHczgNuMLMhphZq5m1AqcC34rJWytPEI6/9gWIb7w/A5+L678kKTd9KnBvLHc6cDow3MweqmG8EwitfX8zeyx/pZm9AcwF/jerE+xPwEeSjgGIx3FfAR6sRkA9qgtnZr+XtDOhCf8H8BQwCLgcuJrQqvQhdEXO7aSKO4CHCa3TNOBPklYQEmMRsLWZTZd0NPAjSZ8HPgTeBU4ys/zx6tYiaRdgGHBk3qpbgEsIXZBSjMpLYoClZpZfb0Fm1iHpK8BV8WxaCzDazJbErtHrwGRJmxMe/1hJ/QjHQ+8Bd2a6fneY2RWl7rtC9xBGsf1JF2UmAFOJz6+ZfRRfr59LGkd4z19mZrOqEZCPieA6lTtmM7OR9Y6lkfWoFsiBpGuBAwqsPrtan7wF9j2HNcf4yxqRO7lTYNsBwJwCq5eZ2YjU+LqDt0DOJehpJxGcqylPIOcSeAI5l6DpTiIMHDjQWltb6x2G6+HmzZv3lpkNKlau6RKotbWV9vb2eofhejhJr5ZSzrtwziXwBHIuQdN14bqD30TXVcpbIOcSeAI5l8ATyLkEnkDOJfAEci6BJ5BzCTyBnEvgCeRcAk8g5xJ4AjmXwBPIuQSeQM4l8ARyLoEnkHMJPIGcS+AJ5FwCTyDnEngCOZfAE8i5BJ5AziWoeQJJ+my8Yez2kraW9Gi8ye0vu7rRrHONqKZvWEl9CXcPy93e7xrCXeBGEG7xd1Qt43EuVa0/8a8GbgT+Eud3I9wRDsLtGQ/qbCNJY+Pt5ts7Ojq6P0rnSlSzBIp3POsws+nZxbb6BkXLCPcZXYuZjTezNjNrGzSo6HDFztVMLQdWPBkwSQcR7hF6C/DZzPoBwDs1jMe5ZBW3QP+/vfMPtqqq4vjnKz9kQhwdSpCRgAoVrLHE0lDhSYmiU5lNWlQoSIyMljY4OqA09mvGmaz8A7WBeGAY4yQJmg3KSBA8kGHE50xpgaWmZdaEaZCi/Fj9sfd7XB7nnHvf3fdx7/Wtz8ydt8+537PPOufuddY6+52zt6SBkk6SNETSfEkjivRmNsHMJppZC2GS32nAakktUTKF/Cn+HKchSUnhfkG4h/khsBdYWEUdc4DvSHoC6E+YXdlxmoaUFO54wrTj3zSzaZIuqnTDGIU6mJhgg+PUlZQI1J8QQZ6SNBY4pjYmOU7zkOJAcwidAD8gTKs+uyYWOU4TkeJAk8zsJjN73czuAi6rlVGO0yx0+x5I0tXATGCMpIvj6qMIKd3cGtrmOA1PNZ0I9wFrgXmE9A3gAPCvWhnlOM1Ct1M4M3vbzF4ErgGGACOAUcBZtTXNcRqflG7sFYROhJfjsgEbki1ynCYixYGGmtn4mlniOE1ISi/cnyQNq5kljngOXQAAAAegSURBVNOEpESg84CXJHW8X2Bm5g7l9CqqdiAzG11LQxynGanagSQtIXQcdGJmM5ItcpwmIiWFuz/+FXAG4Omb0+tISeFK3yx9VNKaGtjjOE1FSgo3uWTxRMI/VR2nV5GSwn25pLyH8Mq24/QqUlK46ZI+DIwFdpjZ07Uzy3Gag5QxEb4BLALGAwsl3VgzqxynSUhJ4aYC55nZvjhg4mbCuG+O02tIeZRHZrYPwMz2EgYWcZxeRUoEapO0gjAU1bnAptqY5DjNQ1UOJGkW4e3TycThec1sQZlt+gGtwEjgaOD7wLPAUsITDX8ArjWzA9XY5Dj1oNspnKTbCI7Tz8x+QxhhdJKk+WU2/SqwMw4kPwVYgA8u7zQ51dwDTQG+aGZvAsS3U68APltmuweAUifbR4WDyztOo1KNA+0uGRAe6OxE2FW0kZntNrNdkgYR3ma9lQoHl/fZGZxGpRoHekvSB0pXxGXL0ZfqhgPrgGVmtpwwGEkHuYPL++wMTqNSTSfCzcAqSWuB54H3AxcCVxZtJGkIsAa4zszWxtXtklrMbD0hNVxXhT2OUzeqGZXnGcLbqO3AQOAp4Bwzay+z6TzCeNrzJa2XtJ6Qxvng8k7TUlU3tpm9Qeh968421wPXZ3zlg8s7TYtP6us4CbgDOU4C7kCOk8CRnCP1iDNq1IuF37/wwsgjYofz7sUjkOMk0JQRyCOL0yh4BHKcBNyBHCcBdyDHScAdyHEScAdynASasheunngPoFOKO1AP4Y7WO/AUznEScAdynATcgRwnAb8HqjPl7pXg4P1Sre+r/D4tHY9AjpOAO5DjJOAO5DgJuAM5TgJ170SQdBRwN3A68DYw08z+XF+rnFK609HR26i7AwGXAgPM7JOSzgZ+hA8yn0Q9e9cq3XetdD1RZ3fOTyOkcOcCjwKY2RbgzPqa4ziV0wgR6FjgjZLl/ZL6dsx+B53zEc2Ki7th1PYudbwX+PdBfe6+KtX1RJ1V6eq5715+fkbkqksxs7p+CHMEXV6y/Lcq6niylrqeqPPdZGNvPT9Zn0ZI4TYBFwPEe6Df19ccx6mcRkjhVgIXSNpMmKVuep3tcZyKqbsDWZgT9ZrEahbWWNcTdb6bbKznvutp42Eo5oCO41RBI9wDOU7T4g7kOAk0vQNJOivOdlek6SdpmaSNkrZKypxRXFIfSa2SNknaIOmDZeo9QdLLkk4to2vvmJVP0pIC3VxJT0jaJunqHM1VJXVtkbRH0nEFx71c0uZ47Jl2Sjo66rZIWiNpdIam8zxL+pCktljnPfFxrCzd2KjbJOluSX3y6ixZNzXOWJi37zMk/b3kHFyRoztB0kPxd9zU8VtmtYVydhZSbf93I3yAmwjd3lvK6KYDd8byYOClHN2lQGsstwAPFdTZj9CDuAM4tUA3AGiv4FhagF8TLmrHALdVsM1dwKyC7z8H/DKWLwB+laO7DlgYy6cAjxWdZ+BhoCWWfwp8Pke3CpgQy0s7dHm/HfBRYG2XdV3rnAnMKdcW4v4uj+XzgUvy2kKRneU+zR6B/gJcVoHuAWB+yfK+LJGZreLgEw8jgH8W1HkHofG8UmbfpwPviVf238b/dWVxIaEBrCQ40iNFlUo6EzjNzIp6kHYAfWOEOBbYm6MbC6wGMLPtwJgu33c9z+OA38XyauDTObovmNkGSf2BoRx6Pg/RShoM3A7cUMG+L4mRZbGkQTm6c4CTJD0OfAVYH9dntYUiO4updxRJ/QAjKROBSrSDCDOBTy2juxf4LzA55/urgFtjeT3FEegjwNcJ/+M6Of7QfTN0iwizmPcnRIHtxF7SnHofBM4vcxzDCZNA7wBeA8bn6GYBi6ONZwP7gT555xl4pWT9JOC+vN+DcCF6DngSGJxVJ9CHEAXGZP2eXfY9HRgXy7cAd+To9gLTY/nbwHeL2kKRnUWfZo9AFSNpOOGELTOz5UVaM7uS0NgXSRqYIZlB+OfvekLa8XNJQ3Oq20FoYGZmO4CdwIkZup2E1OkdC1FgD/C+nGM5juC064qOA/hWrPNkQiS8V9KADF0r4YKxDvgMsM3M9hfUe6CkPAh4PU9oZn81s9GEaP3jHNk4YDRwD3A/MFbSnTnalWa2raMMfCxHt5OQakKI6J0PKWe1hQrtPIxe4UCShhCu7jebWWuB7muS5sbFNwkN5bCGZGYTzGyimbUATwPTzOzVnGpnEF7RQNIwQir1jwxdG3CRAsOAgYRGkMUE4PG84yjhPxx8UPc1wn1b1g3yx4G2eDwrgefL1NsuqSWWpwAbs0SSHi7pkNjFoY7XiZltNbPT4v6/BDxrZl1TuQ4ek/SJWP4UsC1H10Z8RIxwvp6JNh3WFiq1M4u6P4lwhJgHHA/Ml9SR/04xs7e66B4ElkjaQGhsN5jZnsR9LwaWSmoDDJhhJU+ad2Bmj0iaAGwlXNiuLYgCp1C+kQP8BGiVtJGQGs4zs/9l6J4DvifpRkI0yewBLGEOITr3B/4IrMjR3U449ncIF6SZFdhcjtnAgljnqxy8Z82y8WeSZhMuIlPj+qy2cEu1dvqTCI6TQK9I4Rynp3AHcpwE3IEcJwF3IMdJwB3IcRJwB3KcBNyBHCeB/wNVmgriwf5k6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 216x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "numberofplots = len(intervallengthcounter.keys())\n",
    "samples = list(intervallengthcounter.keys())\n",
    "fig,ax = plt.subplots(numberofplots,figsize=[3,10],sharex=True)\n",
    "for i in range(numberofplots):\n",
    "    sample = samples[i]\n",
    "    data = intervallengthcounter[sample]\n",
    "    x = list(data.keys())\n",
    "    y = list(data.values())\n",
    "    sns.barplot(x,y,ax=ax[i],color=\"blue\")\n",
    "    ax[i].set_title(sample)\n",
    "    ax[i].set_ylabel('Counts')\n",
    "plt.tight_layout()\n",
    "fig.savefig(\"Distribution_of_stretches_of_diff_mod_nucleotides_v2.pdf\",dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now running everything through fimo to get the coordinates for the motifs identified by deltaSHAPE and diffBUM-HMM. \n",
    "This would be a useful resource for people studying Xist and Xist RBPs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "directories = glob.glob(\"MEME_V2*\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['MEME_V2deltaSHAPE_rep_1_in_vivo_30_mers',\n",
       " 'MEME_V2deltaSHAPE_rep_1_ex_vivo_30_mers',\n",
       " 'MEME_V2diffBUM_HMM_ex_vivo_30_mers',\n",
       " 'MEME_V2diffBUM_HMM_in_vivo_30_mers',\n",
       " 'MEME_V2deltaSHAPE_rep_2_ex_vivo_30_mers',\n",
       " 'MEME_V2deltaSHAPE_rep_2_in_vivo_30_mers']"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "directories = \" \".join(directories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$directories\"\n",
    "\n",
    "DIR=$HOME/meme/bin\n",
    "\n",
    "for i in $1\n",
    "do\n",
    "    $DIR/fimo \\\n",
    "    --oc $i \\\n",
    "    --verbosity 1 \\\n",
    "    $i/meme_out/meme.txt \\\n",
    "    ../../Reference_sequences/Xist.fa &\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}