{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing various Bismark subsetting options and how they relate to _C.virginica_ genome-wide MBD CpG coverage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Variables to be set by user\n",
    "\n",
    "Re: two \"work_dir\" variables -\n",
    "Both are needed, as one is needed by Bash and the other by Python (respectively)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "env: work_dir=/home/sam/Downloads/20190418_cvir_mbd_coverage_analysis\n"
     ]
    }
   ],
   "source": [
    "%env work_dir = /home/sam/Downloads/20190418_cvir_mbd_coverage_analysis\n",
    "work_dir = \"/home/sam/Downloads/20190418_cvir_mbd_coverage_analysis\"\n",
    "output_plot = \"20190418_cvir_mbd_cov_comparison.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Import necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import os\n",
    "import numpy\n",
    "from IPython.display import display\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Make new working directory, download files and rename using wget ```--output-document``` argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 3.7G\n",
      "-rw-rw-r-- 1 sam sam 934M Apr 17 10:42 avg_reads_CpG_coverage.txt\n",
      "-rw-rw-r-- 1 sam sam 934M Apr 17 10:58 half_avg_reads_CpG_coverage.txt\n",
      "-rw-rw-r-- 1 sam sam 936M Apr 17 11:07 half_total_reads_CpG_coverage.txt\n",
      "-rw-rw-r-- 1 sam sam 937M Apr 17 11:21 total_reads_CpG_coverage.txt\n",
      "  28917406 avg_reads_CpG_coverage.txt\n",
      "  28917406 half_avg_reads_CpG_coverage.txt\n",
      "  28917406 half_total_reads_CpG_coverage.txt\n",
      "  28917406 total_reads_CpG_coverage.txt\n",
      " 115669624 total\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "mkdir ${work_dir}\n",
    "cd ${work_dir}\n",
    "wget \\\n",
    "--quiet \\\n",
    "--output-document avg_reads_CpG_coverage.txt \\\n",
    "http://gannet.fish.washington.edu/Atumefaciens/20190312_cvir_gonad_bismark/avg_reads_bismark/bismark_cytosine_coverage.txt.CpG_report.txt\n",
    "\n",
    "wget \\\n",
    "--quiet \\\n",
    "--output-document half_avg_reads_CpG_coverage.txt \\\n",
    "http://gannet.fish.washington.edu/Atumefaciens/20190312_cvir_gonad_bismark/half_avg_reads_bismark/bismark_cytosine_coverage.txt.CpG_report.txt\n",
    "\n",
    "wget \\\n",
    "--quiet \\\n",
    "--output-document half_total_reads_CpG_coverage.txt \\\n",
    "http://gannet.fish.washington.edu/Atumefaciens/20190312_cvir_gonad_bismark/half_total_reads_bismark/bismark_cytosine_coverage.txt.CpG_report.txt\n",
    "\n",
    "\n",
    "wget \\\n",
    "--quiet \\\n",
    "--output-document total_reads_CpG_coverage.txt \\\n",
    "http://gannet.fish.washington.edu/Atumefaciens/20190312_cvir_gonad_bismark/total_reads_bismark/bismark_cytosine_coverage.txt.CpG_report.txt\n",
    "\n",
    "\n",
    "ls -lh\n",
    "\n",
    "wc -l *.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/sam/Downloads/20190418_cvir_mbd_coverage_analysis\n"
     ]
    }
   ],
   "source": [
    "cd $work_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop through coverage files to calculate percent sequencing coverage for each Bismark subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------\n",
      "----------------------------------------------\n",
      "avg_reads_CpG_coverage.txt\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Total coverage: 4410786'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'5x coverage: 918437'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'10x coverage: 518052'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Mean coverage: 0.7'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'No coverage: 24506620'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent coverage: 15.3'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 5x coverage: 3.2'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 10x coverage: 1.8'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>chrom</th>\n",
       "      <th>pos</th>\n",
       "      <th>strand</th>\n",
       "      <th>meth</th>\n",
       "      <th>unmeth</th>\n",
       "      <th>C-context</th>\n",
       "      <th>trinucleotide</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>49</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>50</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGA</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>51</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGG</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>52</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>88</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGT</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         chrom  pos strand  meth  unmeth C-context trinucleotide  coverage\n",
       "0  NC_007175.2   49      +     0       5        CG           CGC         5\n",
       "1  NC_007175.2   50      -     0       0        CG           CGA         0\n",
       "2  NC_007175.2   51      +     0       5        CG           CGG         5\n",
       "3  NC_007175.2   52      -     0       0        CG           CGC         0\n",
       "4  NC_007175.2   88      +     0       5        CG           CGT         5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------\n",
      "----------------------------------------------\n",
      "total_reads_CpG_coverage.txt\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Total coverage: 13087284'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'5x coverage: 4263854'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'10x coverage: 2617681'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Mean coverage: 7.3'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'No coverage: 15830122'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent coverage: 45.3'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 5x coverage: 14.7'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 10x coverage: 9.1'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>chrom</th>\n",
       "      <th>pos</th>\n",
       "      <th>strand</th>\n",
       "      <th>meth</th>\n",
       "      <th>unmeth</th>\n",
       "      <th>C-context</th>\n",
       "      <th>trinucleotide</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>49</td>\n",
       "      <td>+</td>\n",
       "      <td>2</td>\n",
       "      <td>158</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>50</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGA</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>51</td>\n",
       "      <td>+</td>\n",
       "      <td>2</td>\n",
       "      <td>167</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGG</td>\n",
       "      <td>169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>52</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>18</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>88</td>\n",
       "      <td>+</td>\n",
       "      <td>5</td>\n",
       "      <td>483</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGT</td>\n",
       "      <td>488</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         chrom  pos strand  meth  unmeth C-context trinucleotide  coverage\n",
       "0  NC_007175.2   49      +     2     158        CG           CGC       160\n",
       "1  NC_007175.2   50      -     0      15        CG           CGA        15\n",
       "2  NC_007175.2   51      +     2     167        CG           CGG       169\n",
       "3  NC_007175.2   52      -     0      18        CG           CGC        18\n",
       "4  NC_007175.2   88      +     5     483        CG           CGT       488"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------\n",
      "----------------------------------------------\n",
      "half_total_reads_CpG_coverage.txt\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Total coverage: 9694519'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'5x coverage: 2649503'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'10x coverage: 1647876'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Mean coverage: 3.5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'No coverage: 19222887'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent coverage: 33.5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 5x coverage: 9.2'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 10x coverage: 5.7'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>chrom</th>\n",
       "      <th>pos</th>\n",
       "      <th>strand</th>\n",
       "      <th>meth</th>\n",
       "      <th>unmeth</th>\n",
       "      <th>C-context</th>\n",
       "      <th>trinucleotide</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>49</td>\n",
       "      <td>+</td>\n",
       "      <td>2</td>\n",
       "      <td>34</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>50</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGA</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>51</td>\n",
       "      <td>+</td>\n",
       "      <td>1</td>\n",
       "      <td>35</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGG</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>52</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>88</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>96</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGT</td>\n",
       "      <td>96</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         chrom  pos strand  meth  unmeth C-context trinucleotide  coverage\n",
       "0  NC_007175.2   49      +     2      34        CG           CGC        36\n",
       "1  NC_007175.2   50      -     0       4        CG           CGA         4\n",
       "2  NC_007175.2   51      +     1      35        CG           CGG        36\n",
       "3  NC_007175.2   52      -     0       5        CG           CGC         5\n",
       "4  NC_007175.2   88      +     0      96        CG           CGT        96"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------\n",
      "----------------------------------------------\n",
      "half_avg_reads_CpG_coverage.txt\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Total coverage: 3227641'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'5x coverage: 678974'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'10x coverage: 357365'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Mean coverage: 0.5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'No coverage: 25689765'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent coverage: 11.2'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 5x coverage: 2.3'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Percent 10x coverage: 1.2'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>chrom</th>\n",
       "      <th>pos</th>\n",
       "      <th>strand</th>\n",
       "      <th>meth</th>\n",
       "      <th>unmeth</th>\n",
       "      <th>C-context</th>\n",
       "      <th>trinucleotide</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>49</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>50</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGA</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>51</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGG</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>52</td>\n",
       "      <td>-</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGC</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NC_007175.2</td>\n",
       "      <td>88</td>\n",
       "      <td>+</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>CG</td>\n",
       "      <td>CGT</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         chrom  pos strand  meth  unmeth C-context trinucleotide  coverage\n",
       "0  NC_007175.2   49      +     0       4        CG           CGC         4\n",
       "1  NC_007175.2   50      -     0       0        CG           CGA         0\n",
       "2  NC_007175.2   51      +     0       4        CG           CGG         4\n",
       "3  NC_007175.2   52      -     0       0        CG           CGC         0\n",
       "4  NC_007175.2   88      +     0       5        CG           CGT         5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Variable declaration\n",
    "bismark_subset_list = []\n",
    "mean_seq_coverage = []\n",
    "percent_seq_coverage = []\n",
    "percent_5x_seq_coverage = []\n",
    "percent_10x_seq_coverage = []\n",
    "\n",
    "# Create list of coverage files in current directory\n",
    "cov_files = os.listdir()\n",
    "\n",
    "# Loop through coverage files\n",
    "for file in cov_files:\n",
    "    print(\"----------------------------------------------\")\n",
    "    print(\"----------------------------------------------\")\n",
    "    print (file)\n",
    "    subset_name = file[:-4] # Remove file suffix (.txt)\n",
    "    bismark_subset_list.append(subset_name)\n",
    "    #\n",
    "    #\n",
    "    # Create dataframe and add column names (taken from Bismark documentation)\n",
    "    dataframe = pandas.read_csv(\n",
    "        file,\n",
    "        sep='\\t',\n",
    "        header=None,\n",
    "        names=[\"chrom\", \"pos\", \"strand\", \"meth\", \"unmeth\", \"C-context\", \"trinucleotide\"])\n",
    "    \n",
    "    dataframe['coverage'] = dataframe['meth'] + dataframe['unmeth'] # Sum of methylated and unmethylated coverage for each position.\n",
    "    \n",
    "    total_CpG = len(dataframe) # Count of all CpGs in genome.\n",
    "    \n",
    "    \n",
    "    coverage = sum(dataframe['coverage']>0) # Count of all CpG positions with sequence coverage\n",
    "    coverage_5x = sum(dataframe['coverage']>=5)\n",
    "    coverage_10x = sum(dataframe['coverage']>=10)\n",
    "    mean_coverage = round(dataframe[\"coverage\"].mean(), 1)\n",
    "    \n",
    "    display(\"Total coverage: \" + str(coverage))\n",
    "    display(\"5x coverage: \" + str(coverage_5x))\n",
    "    display(\"10x coverage: \" + str(coverage_10x))\n",
    "    display(\"Mean coverage: \" + str(mean_coverage))\n",
    "    \n",
    "    no_coverage = sum(dataframe['coverage']==0) # Count of all CpG posiitions with no sequence coverage\n",
    "    percent_coverage = round((coverage / total_CpG * 100.0), 1) # Rounds to 1 decimal\n",
    "    percent_5x_coverage = round((coverage_5x / total_CpG * 100.0), 1) # Rounds to 1 decimal\n",
    "    percent_10x_coverage = round((coverage_10x / total_CpG * 100.0), 1) # Rounds to 1 decimal\n",
    "    \n",
    "    mean_seq_coverage.append(mean_coverage)\n",
    "    percent_seq_coverage.append(percent_coverage)\n",
    "    percent_5x_seq_coverage.append(percent_5x_coverage)\n",
    "    percent_10x_seq_coverage.append(percent_10x_coverage)\n",
    "    \n",
    "    display(\"No coverage: \" + str(no_coverage))\n",
    "    display(\"Percent coverage: \" + str(percent_coverage))\n",
    "    display(\"Percent 5x coverage: \" + str(percent_5x_coverage))\n",
    "    display(\"Percent 10x coverage: \" + str(percent_10x_coverage))\n",
    "    display(dataframe.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create new dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>Bismark Subset</th>\n",
       "      <th>Mean Coverage</th>\n",
       "      <th>Percent Coverage</th>\n",
       "      <th>Percent 5x Coverage</th>\n",
       "      <th>Percent 10x Coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>avg_reads_CpG_coverage</td>\n",
       "      <td>0.7</td>\n",
       "      <td>15.3</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>total_reads_CpG_coverage</td>\n",
       "      <td>7.3</td>\n",
       "      <td>45.3</td>\n",
       "      <td>14.7</td>\n",
       "      <td>9.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>half_total_reads_CpG_coverage</td>\n",
       "      <td>3.5</td>\n",
       "      <td>33.5</td>\n",
       "      <td>9.2</td>\n",
       "      <td>5.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>half_avg_reads_CpG_coverage</td>\n",
       "      <td>0.5</td>\n",
       "      <td>11.2</td>\n",
       "      <td>2.3</td>\n",
       "      <td>1.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  Bismark Subset  Mean Coverage  Percent Coverage  \\\n",
       "0         avg_reads_CpG_coverage            0.7              15.3   \n",
       "1       total_reads_CpG_coverage            7.3              45.3   \n",
       "2  half_total_reads_CpG_coverage            3.5              33.5   \n",
       "3    half_avg_reads_CpG_coverage            0.5              11.2   \n",
       "\n",
       "   Percent 5x Coverage  Percent 10x Coverage  \n",
       "0                  3.2                   1.8  \n",
       "1                 14.7                   9.1  \n",
       "2                  9.2                   5.7  \n",
       "3                  2.3                   1.2  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coverage_dataframe = pandas.DataFrame(\n",
    "    {\n",
    "        'Bismark Subset': bismark_subset_list,\n",
    "        'Mean Coverage': mean_seq_coverage,\n",
    "        'Percent Coverage': percent_seq_coverage,\n",
    "        'Percent 5x Coverage': percent_5x_seq_coverage,\n",
    "        'Percent 10x Coverage': percent_10x_seq_coverage,\n",
    "        \n",
    "    })\n",
    "\n",
    "coverage_dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create line plot overlayed on bar chart, showing percent sequencing coverage for each Bismark subset option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAGdCAYAAAACKIgZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XlcVPX+x/HXAYZVwA0VBUVTQ4UB991EQk1Tq6uiecvcqLym3X7XbLmamjeXylyyvFpmt8XlZrlkLllm11wyEhVRcckUBAUVBIaBGeb8/hgYRVERR2fh83w8esAczsx8zkBvv/Odcz5fRVVVhBBCOD4XWxcghBDCOiTQhRDCSUigCyGEk5BAF0IIJyGBLoQQTkICXQghnIQEuhBCOAkJdCGEcBIS6EII4STc7ueT1axZUw0JCbmfTymEEA4vPj4+U1XVgNvtd18DPSQkhN9+++1+PqUQQjg8RVH+LM9+MuUihBBOQgJdCCGchAS6EEI4ifs6h14Wg8FASkoKer3e1qUIB+Xp6UlQUBAajcbWpQhhUzYP9JSUFHx9fQkJCUFRFFuXIxyMqqpcvHiRlJQUGjZsaOtyhLApm0+56PV6atSoIWEuKkRRFGrUqCHv8ITADgIdkDAXd0X+fsSd2nhqIz2/6on2Uy09v+rJxlMbbV2SVdh8ykUIIe6njac2MnXXVPRF5nd1aXlpTN01FYC+jfrasLK7ZxcjdFtzdXUlMjKSsLAwBg0ahE6ns0kd8+bNu+lzv//++zRu3BhFUcjMzLRsV1WV8ePH07hxY7RaLb///vsdPWdycjJ9+vShcePGNGvWjMGDB3P+/Pm7Og4h7Nn83+dbwryEvkjP/N/n26gi65EROuDl5UVCQgIAw4YNY/Hixbz00kvlum9RURGurq5WqWPevHn89a9/xdvb+4afde7cmUcffZTu3buX2r5p0yaOHz/O8ePH2bt3L88//zx79+4t1/Pp9Xr69u3L3Llz6devHwDbt28nIyOD2rVr3/XxXM9oNOLmJn9ywnZyCnNIy0sr82fpeemlbod/Gl6h5zg0/FCF7mcNMkK/TteuXTlx4gQAn3/+Oe3atSMyMpJnn32WoqIiAKpUqcKUKVNo3749u3fvZt++fXTq1ImIiAjatWtHTk4ORUVFTJw4kbZt26LVavn3v/8NwE8//UT37t0ZOHAgoaGhDBs2DFVVWbBgAefOnSMqKoqoqKgb6mrZsiVl9cFZt24dTz/9NIqi0KFDB7KyskhLS+Obb77h4YcfRlVV0tLSaNq0Kenppf9gv/zySzp27GgJc4CoqCjCwsLQ6/WMGDGC8PBwWrZsyfbt2wFo3749hw8ftuzfvXt34uPjycvLY+TIkbRt25aWLVuybt06AJYvX86gQYPo168fPXv2JDc3l+joaFq1akV4eLhlP4A333yT0NBQYmJiGDp0KO+88w4AJ0+epHfv3rRu3ZquXbty9OjRO/69CrH9zHYeW/fYTX9ex6fOfazm3rCr4dK0DYdJOnfFqo/ZvK4fb/RrUa59jUYjmzZtonfv3hw5coRVq1bxyy+/oNFoGDt2LF988QVPP/00eXl5hIWFMX36dAoLCwkNDWXVqlW0bduWK1eu4OXlxccff4y/vz/79u2joKCAzp0707NnTwD279/P4cOHqVu3Lp07d+aXX35h/PjxzJ07l+3bt1OzZs1yH19qairBwcGW20FBQaSmpvL444+zZs0aFi1axObNm5k2bRp16pT+g01MTKR169ZlPu6iRYsAOHToEEePHqVnz54kJyczZMgQVq9ezbRp00hLS+PcuXO0bt2a1157jR49erBs2TKysrJo164dDz/8MAC7d+/m4MGDVK9eHaPRyDfffIOfnx+ZmZl06NCB/v37Ex8fz5o1a9i/fz9Go5FWrVpZaouLi2Px4sU0adKEvXv3MnbsWH788cdyv0aicsvMz2TWr7PYcnoLTao1YWCTgSxLXFZq2sXT1ZMJrSbYsErrsKtAt5X8/HwiIyMB8wh91KhRLFmyhPj4eNq2bWvZp1atWoB5zv0vf/kLAMeOHSMwMNCyn5+fHwBbt27l4MGDfPXVVwBkZ2dz/Phx3N3dadeuHUFBQQBERkZy+vRpunTpUqHaVVW9YVvJWR8LFy4kLCyMDh06MHTo0Dt63J07d/LCCy8AEBoaSoMGDUhOTmbw4MHExMQwbdo0Vq9ezaBBgyzHu379esuoWq/Xc+bMGQBiYmKoXr26pd7XXnuNn3/+GRcXF1JTUzl//jw7d+5kwIABeHl5AVjeNeTm5rJr1y7L8wAUFBTc0bGIyklVVdafXM/bv72NzqBjXOQ4RoaNROOqob5ffeb/Pp/0vHTq+NRhQqsJDv+BKNhZoJd3JG1t186hl1BVleHDhzNz5swb9vf09LTMm6uqWuZpc6qqsnDhQnr16lVq+08//YSHh4fltqurK0ajscK1BwUFcfbsWcvtlJQU6tatC5hH7y4uLpw/fx6TyYSLS+kZthYtWrBjx44yH7esfygA6tWrR40aNTh48CCrVq2yTCWpqsqaNWt48MEHS+2/d+9efHx8LLe/+OILMjIyiI+PR6PREBISgl6vv+nzmUwmqlatesPvR4hbSc1NZfru6ew6t4uWtVoyteNUGlVtZPl530Z9nSLArydz6DcRHR3NV199xYULFwC4dOkSf/55YwfL0NBQzp07x759+wDIycnBaDTSq1cvPvzwQwwGA2A+myQvL++Wz+nr60tOTs4d1dm/f3/+85//oKoqe/bswd/fn8DAQIxGIyNGjODLL7+kWbNmzJ0794b7Pvnkk+zatYuNG6+eg7t582YOHTpEt27d+OKLLyy1nzlzxhLWQ4YMYc6cOWRnZxMebv7gqFevXixcuNASzPv37y+z3uzsbGrVqoVGo2H79u2W17RLly5s2LABvV5Pbm6upSY/Pz8aNmzIf//7X8D8D8eBAwfu6DUSlUeRqYjPkj7j8XWPk3Ahgdfav8by3stLhbkzk0C/iebNmzNjxgx69uyJVqslJiaGtLQbPx13d3dn1apVvPDCC0RERBATE4Ner2f06NE0b96cVq1aERYWxrPPPnvbkXhcXByPPPJImR+KLliwgKCgIFJSUtBqtYwePRqAPn360KhRIxo3bsyYMWP44IMPAHjrrbfo2rUrXbt2Ze7cuXz00UccOXKk1GN6eXnx7bffsnDhQpo0aULz5s1Zvnw5tWrVYuzYsRQVFREeHk5sbCzLly+3vLMYOHAgK1euZPDgwZbHmjx5MgaDAa1WS1hYGJMnTy7zGIcNG8Zvv/1GmzZt+OKLLwgNDQWgbdu29O/fn4iICJ544gnatGmDv78/YB7Vf/zxx0RERNCiRYtSH6QKUeL45eM8velp5uybQ5vabVg7YC1DQ4fiolSemFNu9lb3XmjTpo16/QIXR44coVmzZvetBmG/cnNzqVKlCjqdjm7durFkyRJatWpVrvvK31HlVVhUyNJDS/no0Ef4anyZ1G4SfRr2uesriO3ptEVFUeJVVW1zu/3sag5dVG5xcXEkJSWh1+sZPnx4ucNcVF4JFxJ4Y9cbnMo+Rd9GfXm57ctU96xu67JsRgJd2I0vv/zS1iUIB6Ez6Jj/+3xWHF1BbZ/afBD9AV2Dutq6LJuTQBdCOJSdqTuZvns66XnpDAkdwoRWE/DR+Nz+jpWABLoQwiFc1l9mzr45fHvqWxr6N+Q/j/yHyFqRti7LrkigCyHsmqqqbPpjE7P3zeZKwRWe1T5LnDYOd1d3W5dmdyTQhRB2Kz0vnTf3vMnPKT8TXjOcpT2X0rRaU1uXZbcqzwmat+AI7XOfeeYZGjZsSGRkJJGRkXd85eSmTZto06YNzZo1IzQ0lH/84x/WKFmIe8Kkmlh5dCWPrXuMfen7mNhmIp898pmE+W1IoHP10v/ExETc3d1ZvHhxue9b0oHRGm4V6ABvv/02CQkJJCQkWHrPlEdiYiLjxo3j888/58iRIyQmJtKo0b27cu5uWhkIcSr7FCM2j+Bfe/9FeM1wvu7/NU+3eBpXF+u0qXZmEujXsdf2uTczfvx4pk+fDsCWLVvo1q0bJpOp1D5z5szh9ddft1yV6ebmxtixYwH4888/iY6ORqvVEh0dzZkzZ8jOziYkJMTyODqdjuDgYAwGw01b2T7zzDO89NJLREVFMWnSJH799Vc6depEy5Yt6dSpE8eOHbM81uDBg9FqtcTGxtK+fXtKLjbbunUrHTt2pFWrVgwaNIjc3Nw7/wUKh2UwGVhycAkD1w/kRNYJ3uz8JktilhDkG2Tr0hyGfc2hb3oF0q18lVWdcHhkVrl2tff2ua+//jrTp08nOjqaWbNm4eHhwaxZs2jbti1du3Zl/PjxfPfddzc04UpMTOT//u//ynzMcePG8fTTTzN8+HCWLVvG+PHjWbt2LREREezYsYOoqCg2bNhAr1690Gg0t2xlm5yczLZt23B1deXKlSv8/PPPuLm5sW3bNl577TXWrFnDBx98QLVq1Th48CCJiYmWdxqZmZnMmDGDbdu24ePjw+zZs5k7dy5Tpkwp1+9OOLbDmYeZsmsKyZeT6dmgJ6+2f5WaXuVvIy3M7CvQbcQR2ufOnDmTOnXqUFhYSFxcHLNnz2bKlCl4e3uzdOlSunXrxnvvvccDDzxwR8e+e/duvv76awCeeuopXn75ZQBiY2NZtWoVUVFRrFy5krFjx962le2gQYMsXSizs7MZPnw4x48fR1EUS5OynTt3MmGCue90WFgYWq0WgD179pCUlETnzp0BKCwspGPHjnd0LMLx5BvzWbR/EZ8d+YyanjWZHzWfHvV72Losh2VfgV7OkbS1OUL73MDAQAA8PDwYMWKEpe84mBehqFGjBufOnSvzvi1atCA+Pp6IiIjbPk/JsfTv359XX32VS5cuER8fT48ePcjLy7tlK9tr2+ROnjyZqKgovvnmG06fPm1ZOu9mvYNUVSUmJoYVK1bctkbhHPak7WHarmmk5KYwsOlAXmr9Er7uvrYuy6GVew5dURRXRVH2K4rybfHthoqi7FUU5biiKKsURXGqk0LtrX1uSadHVVVZu3YtYWFhgHkO/N1332X//v1s2rSpzPVEJ06cyFtvvUVycjJg7jFe0k63U6dOrFy5EjB3NSx5p1ClShXatWvHhAkTePTRR3F1db2jVrbZ2dnUq1cPMC9DV6JLly6sXr0agKSkJA4dMk+xdejQgV9++cXy+YVOp7PUK5xLdkE2U36ZwpitY3BRXFjWaxlvdHxDwtwK7uRD0QnAtf1XZwPvqaraBLgMjLJmYbZmb+1zhw0bRnh4OOHh4WRmZvLPf/4TVVUZNWoU77zzDnXr1uXjjz9m9OjR6PWlVzTXarXMmzePoUOH0qxZM8LCwizHsmDBAj755BO0Wi2fffYZ8+dfXfk8NjaWzz//nNjYWMu28rayffnll3n11Vfp3LlzqTOBxo4dS0ZGBlqtltmzZ6PVavH39ycgIIDly5czdOhQtFotHTp0kLVDndD3f37PY+seY/3J9YwMG8ma/mtoW6etrctyGuVqn6soShDwKfAv4CWgH5AB1FFV1agoSkdgqqqqvW7xMNI+V1BUVITBYMDT05OTJ08SHR1NcnIy7u539wZP/o7s2wXdBd7a+xY/nPmBZtWbMbXTVJrXaG7rsm7JmdvnzgNeBkreE9UAslRVLRlypgD17rhKUenodDqioqIwGAyoqsqHH35412Eu7Jeqqnx9/Gve/e1dCk2FvNjqRYa3GI6bi319fOcsbvuqKoryKHBBVdV4RVG6l2wuY9cyh/qKosQBcQD169evYJnCWfj6+nL9uzThnM5cOcO03dP4Nf1X2tRuw9ROU2ng18DWZTm18vwz2RnoryhKH8AT8MM8Yq+qKIpb8Sg9CCjzFAtVVZcAS8A85WKVqoUQdstoMvJZ0mcsSliExkXDlI5T+EuTv1SqpeBs5baBrqrqq8CrAMUj9H+oqjpMUZT/AgOBlcBwQBZ6FKKSO3rpKG/seoOki0lEBUfxevvXqe1T29ZlVRp3M5E1CVipKMoMYD/wsXVKEkI4moKiAhYfWMwniZ/g7+HPOw+9Q88GPe96XU9xZ+4o0FVV/Qn4qfj7U0A765ckhHAkv6X/xrTd0zh95TQDHhjAxLYT8ffwt3VZlZJMauEY7XMNBgOvvPIKTZo0ISwsjHbt2rFp06b7XKEQV+UW5vLm7jcZsWUEBpOBf8f8mxldZkiY25AEOo7RPnfy5MmkpaWRmJhIYmIiGzZsuOlVpdZgzeMSzuensz8xYN0Avjr+FU81f4qv+39Np7qdbF1WpSeBfh17bJ+r0+lYunQpCxcutPSBqV27NoMHDwZgxYoVhIeHExYWxqRJkwD48MMPLY22wHz5/QsvvHBHxzV9+nTatm1LWFgYcXFxlj4s+/btQ6vV0rFjRyZOnGhpQ3CzYxbO42L+RSbumMgLP76An7sfnz/yOS+3fRlvjbetSxPYWXOu2b/O5ugl617uHVo9lEntJpVrX3ttn3vixAnq169v6eR4rXPnzjFp0iTi4+OpVq0aPXv2ZO3atQwcOJCOHTsyZ84cAFatWsXrr79e7uMCc/uDkva1Tz31FN9++y39+vVjxIgRLFmyhE6dOvHKK69YarnZMTds2PDOf3HCrqiqyoZTG5izbw46g46/Rf6NUWGj0LhqbF2auIZdBbqtOEL73JvZt28f3bt3JyAgADD3fPn555957LHHaNSoEXv27KFJkyYcO3aMzp07s2jRonIdF8D27duZM2cOOp2OS5cu0aJFC7p27UpOTg6dOpnfXj/55JN8++23tzxmCXTHlpqbyvTd09l1bheRAZFM6zSNRlXv3YpXouLsKtDLO5K2Nntvn9u4cWPOnDlDTk4Ovr6lO9LdqhdPbGwsq1evJjQ0lMcffxxFUcp9XHq9nrFjx/Lbb78RHBzM1KlT0ev1t3y+mx2zcExFpiJWHF3Bgv0LUFB4td2rDAkdIhcI2TH5zdyEPbXP9fb2ZtSoUYwfP57CwkLA3E73888/p3379uzYsYPMzEyKiopYsWIFDz30EABPPPEEa9euZcWKFZaOieU9rpKOjTVr1iQ3N9cy6q5WrRq+vr7s2bMHwNJ6F6jQMQv7dOLyCZ7e9DSz982mde3WrB2wliebPSlhbufsaoRuT65tn2symdBoNCxatIgGDUr3ori2fW5+fj5eXl5s27aN0aNHc/r0aVq1aoWqqgQEBLB27dpbPmdJ+9zAwEC2b99e6mczZszgn//8J82bN8fT0xMfHx+mT59OYGAgM2fOJCoqClVV6dOnDwMGDADM4du8eXOSkpJo167dHR1X1apVGTNmDOHh4YSEhFimaMA8Vz5mzBh8fHzo3r07/v7m09QqcszCvhQWFfLRoY9YemgpVTRVmNl1Jn0b9pULhBxEudrnWou0z3UOubm5VKlSBYBZs2aRlpZWqo+6Lcjf0d1LuJDA1F1TOZl9kr6N+vJy25ep7lnd1mXZjDO3zxXCYuPGjcycOROj0UiDBg1KrUgkHI/OoGPB/gV8eeRLavvUZlH0IroFdbN1WaICJNDFHYuNjS21ipFwXL+k/sL03dNJy0sj9sFYXmz9Ij4an9vfUdglCXQhKqEsfRZz9s1hw6kNNPRvyKePfErLWi1tXZa4SxLoQlQiqqqy+fRmZv06iysFV3hW+yxjtGPwcPW4/Z2F3ZNAF6KSSM9LZ8aeGexI2UFYjTCW9lxK02pNbV2WsCIJdCGcnEk18d9j/+W939+jyFTExDYTGdZsGK4urrYuTViZ3QV6yCsbrfp4p2f1ve0+rq6uhIeHYzQaadasGZ9++ine3ve/2dC8efOIi4sr87nff/995s2bx8mTJ8nIyLD0e1FVlQkTJvDdd9/h7e3N8uXLadWqVbmfMzk5mRdffJHk5GQ0Gg3h4eEsXLiQ2rVllRln8Ef2H0zdNZXfL/xOh8AOTOk4hWDfYFuXJe4RuewLx2if27lzZ7Zt23bDBUCbNm3i+PHjHD9+nCVLlvD888+X+/n0ej19+/bl+eef58SJExw5coTnn3+ejIyMuzqOm7ldiwNhPQaTgaUHlzJw/UCOZx1neqfpLIlZImHu5CTQr2OP7XMBWrZsSUhIyA3b161bx9NPP42iKHTo0IGsrCzS0tL45ptvePjhh1FVlbS0NJo2bUp6enqp+3755Zd07NiRfv36WbZFRUURFhaGXq9nxIgRhIeH07JlS8uVq+3bt+fw4cOW/bt37058fDx5eXmMHDmStm3b0rJlS9atMy8xu3z5cgYNGkS/fv3o2bMnubm5REdH06pVK8LDwy37Abz55puEhoYSExPD0KFDeeeddwA4efIkvXv3pnXr1nTt2pWjR63bkdPZHM48zJBvh7Bg/wK6B3dn/WPrebzJ43K1ZyVgd1MutmSv7XNvJTU1leDgq6OuoKAgUlNTefzxx1mzZg2LFi1i8+bNTJs2jTp16pS6b2JiIq1bty7zcRctWgTAoUOHOHr0KD179iQ5OZkhQ4awevVqpk2bRlpaGufOnaN169a89tpr9OjRg2XLlpGVlUW7du14+OGHAdi9ezcHDx6kevXqGI1GvvnmG/z8/MjMzKRDhw7079+f+Ph41qxZw/79+zEajbRq1cpSW1xcHIsXL6ZJkybs3buXsWPH8uOPP5b/F1tJ5BvzWbR/EZ8d+YyanjWZHzWfHvV72LoscR9JoOPY7XPLat1QMhJbuHAhYWFhdOjQgaFDh97R4+7cudOyIEZoaCgNGjQgOTmZwYMHExMTw7Rp01i9ejWDBg2yHO/69esto2q9Xs+ZM2cAiImJoXr16pZ6X3vtNX7++WdcXFxITU3l/Pnz7Ny5kwEDBuDl5QVgedeQm5vLrl27LM8DUFBQcEfHUhnsTdvL1F1TSclNYWDTgfy99d/xc7+xf75wbhLo2H/73FsJCgri7NmzltspKSnUrVsXMI/eXVxcOH/+PCaTCReX0jNsLVq0YMeOHWU+7s16/NSrV48aNWpw8OBBVq1aZZlKUlWVNWvW8OCDD5baf+/evfj4XL3y8IsvviAjI4P4+Hg0Gg0hISG3bMtrMpmoWrXqDb8fYZZdkM3c+Ll8ffxr6vvWZ1mvZbSt0/b2dxROSebQb8Ke2ufeSv/+/fnPf/6Dqqrs2bMHf39/AgMDMRqNjBgxgi+//JJmzZoxd+7cG+775JNPsmvXLjZuvHpm0ebNmzl06BDdunXjiy++sNR+5swZS1gPGTKEOXPmkJ2dTXi4uYFRr169WLhwoSWY9+/fX2a92dnZ1KpVC41Gw/bt2y2vaZcuXdiwYQN6vZ7c3FxLTX5+fjRs2JD//ve/gPkfjgMHDtzRa+Ssvv/zex5b9xjrTqxjRNgI1vRfI2FeydndCL08pxneD/bWPnfBggXMmTOH9PR0tFotffr04aOPPqJPnz589913NG7cGG9vbz755BMA3nrrLbp27UrXrl2JjIykbdu29O3bt1RHQi8vL7799ltefPFFXnzxRTQaDVqtlvnz5zN27Fiee+45wsPDcXNzY/ny5ZZ3FgMHDmTChAlMnjzZ8liTJ0/mxRdfRKvVoqoqISEhlpWMrjVs2DD69etHmzZtiIyMJDQ0FIC2bdvSv39/IiIiaNCgAW3atLG05f3iiy94/vnnmTFjBgaDgSFDhhAREVHeX6VT2HhqI/N/n096XjoB3gHU9KxJ0qUkQquHsih6Ec1rNLd1icIOSPtcYTdK2vLqdDq6devGkiVLyn1OvTP/HW08tZGpu6aiL9KX2t67QW/e6vYWGhdZ1/NekPa5QtyFuLg4kpKS0Ov1DB8+/I4ukHJm83+ff0OYAxzIPCBhLkqRQBd248svv7R1CXYpPS+9XNvtaUQpbEM+FBXCjhUWFeLu6l7mz+r41Clzu6i8JNCFsFOFRYX8/ae/U1BUcMPUiqerJxNaTbBRZcJeyZSLEHaooKiAv2//O/9L/R+TO0zGR+NjOculjk8dJrSaQN9G9nFGmLAfEuhC2JmCogJe3P4iO1N3MqXjFAY1NV8lKwEubsf+An2qv5UfL/u2uzhC+9xnnnmGHTt2WM7NXr58uaVdQXls2rSJyZMnk5eXh6qqPProo5bL9IX9KCgqYML2CfyS+gtvdHyDgU0H2rok4UBkDh3HaJ8L8Pbbb5OQkEBCQsIdhXliYiLjxo3j888/58iRIyQmJtKoUSNrlFwmaZNbMQVFBUz4cQK7UncxrdM0CXNxxyTQr2Ov7XNvZvz48UyfPh2ALVu20K1bN0wmU6l95syZw+uvv265KtPNzY2xY8cC8OeffxIdHY1WqyU6OpozZ86QnZ1NSEiI5XF0Oh3BwcEYDIabtrJ95plneOmll4iKimLSpEn8+uuvdOrUiZYtW9KpUyeOHTtmeazBgwej1WqJjY2lffv2lFxstnXrVjp27EirVq0YNGgQubm5d/4LdFB6o57xP45n1zlzmD/R5AlblyQckAT6NUra54aHh5dqn5uQkICrq6ult0lJ+9y9e/fSrl07YmNjmT9/PgcOHGDbtm03tM/dt28fS5cu5Y8//gDMfU7mzZtHUlISp06dsrTPrVu3Ltu3b7/hsv8Sr7/+Olqtlr///e+WjoOzZs1i1apVbN++nfHjx/PJJ5/c0ITrVm1yx40bx9NPP83BgwcZNmwY48ePx9/fn4iICEvjrg0bNtCrVy80Gg1xcXEsXLiQ+Ph43nnnHcs/DGDu+bJt2zbeffddQkND+fnnn9m/fz/Tp0/ntddeA+CDDz6gWrVqHDx4kMmTJxMfHw9AZmYmM2bMYNu2bfz++++0adOmzP4zzijfmM8LP77A7nO7mdZpGo83edzWJQkHZX9z6DbgCO1zZ86cSZ06dSgsLCQuLo7Zs2czZcoUvL29Wbp0Kd26deO9997jgQceuKNj3717N19//TUATz31FC+//DIAsbGxrFq1iqioKFauXMnYsWNv28p20KBBli6U2dnZDB8+nOPHj6MoiqVJ2c6dO5kwwXy6XVhYGFqtFoA9e/aQlJRE586dASgsLKRjx453dCyOqCTMf037lTc7v8mAxgNsXZJwYBLoOEb73MDAQAB58DGOAAAgAElEQVQ8PDwYMWJEqQ80Dx06RI0aNTh37lyZ923RogXx8fHlamhVciz9+/fn1Vdf5dKlS8THx9OjRw/y8vJu2cr22ja5kydPJioqim+++YbTp0/TvXt34OZteVVVJSYmhhUrVty2RmeRb8znhR9e4Nf0X5nRZQb9H+hv65KEg5Mpl5uwt/a5aWlpgDn41q5dS1hYGGCeA3/33XfZv38/mzZtYu/evTfcd+LEibz11lskJycD5h7jJdMZnTp1YuXKlYC5q2HJO4UqVarQrl07JkyYwKOPPoqrq+sdtbLNzs6mXr16gPmMnBJdunRh9erVACQlJXHokPmy8w4dOvDLL79YPr/Q6XSWep2RzqBj3A/j+DX9V/7V5V8S5sIqbjtCVxTFE/gZ8Cje/ytVVd9QFKUhsBKoDvwOPKWqauFdV1SO0wzvB3trnzts2DAyMjJQVZXIyEgWL16MqqqMGjWKd955h7p16/Lxxx/zzDPPsG/fPjw9PS331Wq1zJs3j6FDh6LT6VAUhb59zec0L1iwgJEjR/L2228TEBBgab8L5mmXQYMG8dNPP1m2lbeV7csvv8zw4cOZO3cuPXpcXQZt7NixDB8+HK1WS8uWLdFqtfj7+xMQEMDy5csZOnSoZRpnxowZNG3a9Da/KcejM+gY9+M44s/H868u/6LfA/1ufychyuG27XMV83twH1VVcxVF0QA7gQnAS8DXqqquVBRlMXBAVdUPb/VY0j5XFBUVYTAY8PT05OTJk0RHR5OcnIy7e9n9SsrLUf6OdAYdf/vhb/x+4Xfe6vKWVS8WkuZc1mVPr6fV2ueq5sQvOX9MU/yfCvQAnize/ikwFbhloAuh0+mIiorCYDCgqioffvjhXYe5o9AZdIz9YSz7L+xnZpeZ9GnUx9YlCSdTrg9FFUVxBeKBxsAi4CSQpapqyad5KUC9e1KhcCq+vr5c/y6tMtAZdDy/7XkOZBxgdtfZ9G7Y29YlCSdUrg9FVVUtUlU1EggC2gFlvbctc+5GUZQ4RVF+UxTlt4yMjIpXKoSDyjPkWcJ8VrdZEubinrmjs1xUVc0CfgI6AFUVRSkZ4QcBZZ4zp6rqElVV26iq2iYgIOBuahXC4Vwb5rO7zaZ3iIS5uHduG+iKogQoilK1+Hsv4GHgCLAdKGk2MRxYd6+KFMIR5Rbm8tz3z3Ew4yBzus2hV0iv299JiLtQnjn0QODT4nl0F2C1qqrfKoqSBKxUFGUGsB/4+B7WKYRDyS3M5bltz3E48zBvP/Q2MQ1ibF2SqATKc5bLQaBlGdtPYZ5Pt6qKnip0M+U5hcgR2ucaDAYmT57MmjVr8PDwwNvbm2nTpvHII4/c9zrFreUU5vDctudIykzinYfeIbpBtK1LEpWEXCmKY7TPnTx5MmlpaSQmJpKYmMiGDRtuelWpNVjzuCqTK4VXePb7Z81h3l3CXNxfEujXscf2uTqdjqVLl7Jw4UJLH5jatWszePBgAFasWEF4eDhhYWFMmjQJgA8//NDSaAvMl9+/8MILd3Rc06dPp23btoSFhREXF2fpw7Jv3z60Wi0dO3Zk4sSJljYENzvmyuJK4RWe3fosRy4d4d3u7xJdX8Jc3F8S6New1/a5J06coH79+pZOjtc6d+4ckyZN4scffyQhIYF9+/axdu1aBg4caOmiCLBq1SpiY2PLfVxdunRh3Lhx7Nu3j8TERPLz8/n2228BGDFiBIsXL2b37t2WJmXALY/Z2WUXZBO3NY6jl48y96G59Kjf4/Z3EsLKpNsijtE+92b27dtH9+7dKTkldNiwYfz888889thjNGrUiD179tCkSROOHTtG586dWbRoUbmOC2D79u3MmTMHnU7HpUuXaNGiBV27diUnJ4dOnToB8OSTT1qC/mbH3LBhwwodm6PILsgm7vs4ki8n81739+ge3N3WJYlKSgId+2+f27hxY86cOUNOTg6+vr43PM/NxMbGsnr1akJDQ3n88cdRFKXcx6XX6xk7diy//fYbwcHBTJ06Fb1ef8vnu9kxO7PsgmzGbB3DiawTzOs+j4eCH7J1SaISkymXm7Cn9rne3t6MGjWK8ePHU1hobmiZlpbG559/Tvv27dmxYweZmZkUFRWxYsUKHnrIHCpPPPEEa9euZcWKFcTGxt7Rcen1egBq1qxJbm6uZdRdrVo1fH192bNnD4Cl9S5QoWN2ZKXCPErCXNie3Y3Q7aXzm721z50xYwb//Oc/ad68OZ6envj4+DB9+nQCAwOZOXMmUVFRqKpKnz59GDDAvOpNtWrVaN68OUlJSbRr1+6Ojqtq1aqMGTOG8PBwQkJCLFM0YJ4rHzNmDD4+PnTv3h1/f3+ACh2zo8rSZzHm+zGcyjrF/Kj5dA3qauuShLh9+1xrkva5ziE3N5cqVaoA5jVN09LSmD9/vk1rup9/R6XCvMd8utSr2Ocf1mZP7V6dgT29nlZrnyvE9TZu3MjMmTMxGo00aNCg1IpEzu6y/jJjto7hj+w/WNBjAZ3rdbZ1SUJYSKCLOxYbG2uZk69MLusvM3rraP688icLeyykU71Oti5JiFLs4kPR+zntI5zP/fj7uaS/xKitoyTMhV2zeaB7enpy8eJFCXVRIaqqcvHixVJrqFrbxfyLjNoyirNXzvJ+9Pt0rNvxnj2XEHfD5lMuQUFBpKSkIItfiIry9PS0XKhlbRfzLzJ662hSclJ4P/p92ge2vyfPI4Q12DzQNRqN019JKBxTZn4mo7eMJjU3lUXRi2gXaPXmokJYlc0DXQh7lJmfyagto0jLS+ODhz+gbZ22t7+TEDYmgS7EdTJ0GYzaOor0vHQWRS+SMBcOQwJdiGtk6DIYuWUk53Xn+SD6A9rUue21HELYDQl0IYpd0F1g1JZRnNed58OHP6R17da2LkmIOyKBLgRXw/yC7gKLH15Mq9qtbF2SEHdMAl1UeufzzjNq6ygydBksjllMy1o3LKErhEOQQBeVWnpeOqO2jOKi/iL/jvk3kbUibV2SEBUmgS4qrfS8dEZuGckl/SUWP7xYwlw4PAl0USml56UzYvMIsgqy+HfMv4kIiLB1SULcNQl0Uemk5aYxcstIS5hrA7S2LkkIq5BAF5XKudxzjNwykisFV1gSs4TwgIotYiCEPZJAF5VGam4qo7aMMod5zyWE1QyzdUlCWJUEuqgUUnNTGbl5JDmGHJb2XEqLmi1sXZIQVieBLpxeSk4KI7eMJM+QZw7zGhLmwjnZfIELIe6lszlnJcxFpSEjdOG0zl45y8itI8k35vNRz49oVqOZrUsS4p6SQBdO6eyVs4zYMgJ9kZ6Pen5EaPVQW5ckxD0ngS6czpkrZxixZQSFRYV83PNjHqz+oK1LEuK+kDl04VT+vPInI7aMwFBk4KOeH0mYi0pFRujCaZzOPs2oLaMwmAx81OsjmlZrauuShLivJNCFU/gj+w9GbRlFkVrEx70+pkm1JrYuSYj7TqZchMM7lX3qapj3lDAXlZeM0IVDKwlzk2piWa9lPFD1AVuXJITN3HaErihKsKIo2xVFOaIoymFFUSYUb6+uKMr3iqIcL/5a7d6XK8RVp7JOMXLzSFRVlTAXgvJNuRiB/1NVtRnQAfiboijNgVeAH1RVbQL8UHxbiPviZNZJRmwZASBhLkSx2wa6qqppqqr+Xvx9DnAEqAcMAD4t3u1T4LF7VaQQ1zpx+QQjt4zERXFhWe9lNKrayNYlCWEX7uhDUUVRQoCWwF6gtqqqaWAOfaDWTe4TpyjKb4qi/JaRkXF31YpK7/jl44zaOgpXxZVlvZbRyF/CXIgS5Q50RVGqAGuAF1VVvVLe+6mqukRV1TaqqrYJCAioSI1CAJB8OZnRW0fjprixrNcyGvo3tHVJQtiVcgW6oigazGH+haqqXxdvPq8oSmDxzwOBC/emRCHg2KVjjN5SHOa9lxHiH2LrkoSwO+U5y0UBPgaOqKo695ofrQeGF38/HFhn/fKEKA7zraPRuGpY1nsZDfwa2LokIexSec5D7ww8BRxSFCWheNtrwCxgtaIoo4AzwKB7U6KozErC3MPVg2W9llHfr76tSxLCbt020FVV3QkoN/lxtHXLEeKqo5eOMnrraLzcvFjWcxnBfsG2LkkIuyaX/gu7dOTiEUZtGYW3mzfLekmYC1EeEujC7iRdTGL01tH4aHzMYe4rYS5EeUigC7ty+OJhRm8dTRVNFZb1WkaQb5CtSxLCYUhzLmE3DmceZsz3Y/Bz9+PjXh9Tr0o9W5ckhEOREbqwC4mZiYzZag7zZb2WSZgLUQES6MLmDmUcMoe5hx+f9PqEulXq2rokIRySBLqwqYMZB4n7Po6qHlX5pNcnBFYJtHVJQjgsCXRhMwcyDhD3fRzVPKvxSW8JcyHulgS6sImECwk8+/2zVPeszrJey6jjU8fWJQnh8CTQxX1XEuY1PGvwSa9PJMyFsBIJdHFf7b+wn2e/f5YA7wCW9VpGbZ/ati5JCKchgS7um9/P/86z3z9LLe9aEuZC3AMS6OK+iD8fz3PbnqO2d22W9VpGLe8yF7gSQtwFuVJU3BMbT21k/u/zSc9Lp7pnda4UXCHIL4iPe35MgLesXCXEvSCBLqxu46mNTN01FX2RHoCL+osoKDwZ+qSEuRD3kEy5CKtSVZV3f3vXEuaW7agsS1xmo6qEqBxkhC7uis6g4/DFwyRcSOBAxgEOZhzkcsHlMvdNz0svdTv80/AKPeeh4YcqdD8hLA6uhh+mQ3YK+AdB9BTQDrZ1VXdNAl2Um6qqpOSmcCDjAAcuHOBAxgGSLydTpBYBEOIXwkPBD7H9zHayC7NvuL+cby7swsHVsGE8GPLNt7PPmm+Dw4e6BLq4Kb1RT9LFJBIyEiwBflF/EQBvN2/Ca4YzMmwkkbUiiQiIwN/DH4CNgaXn0AE8XT2Z0GqCTY5DCIvcDNj08tUwL2HIN4/YJdCFM1BVlfS8dPPoO+MACRcSOHrpKEbVCECwbzCd6nYiIiCCiFoRNK7aGDeXsv98+jbqC2A5y6WOTx0mtJpg2S7EfaOqMK1q+fbNPgtT/a/ebuh4C5JLoFdShUWFJF1MsgT4gYwDXNBdAMyj6RY1WzC8xXAiAiLQBmip4VXjjh6/b6O+EuDCNowFcPp/cGwzJG+2dTX3lQR6JXFBd8Ey8j6QcYCki0kYTAYA6lWpR+varYkIiCCyViRNqzVF46KxccVC3IHcDDi+xRzgJ7dDYS64ecEDUeaRdyUhge6EDCYDxy4ds3x4mZCRQFpeGgDuLu40r9GcYc2GmadPAiLk3HDheFQVLiTBsU3mEE/5DVDBty6ED4IHH4GG3UDjVXoaxclJoDuBi/kXzaPv4g8vky4mWT6QrOVdi8iASP7a7K9E1ooktHoo7q7uNq5YiAowFsDpneYAT94MWWfM2+u2hO6vQNPeEBgBimLbOm1IAt3BGE1Gjl8+XurDy5TcFADcXNxoVr0ZA5sOJKJWBJEBkXKqoHBseZlwfKt5JH7yx6tTKY26Q9f/gya9wE8WRikhgW7nsvRZHMw8aJn7PpR5iHyj+ZSrml41iQiIYPCDg4msFUmz6s3wdPO0ccVC3AVVhQtHIHkTJG+Bs79inkoJhPCB0LR4KsXd29aV2iUJdDtSZCriZPbJUhfunL5yGgBXxZWm1ZryWOPHLHPf9arUQ6nEby+FkzAWwp87r56VkvWneXtgBDw0CR7sDYGRlXoqpbwk0G3oSuEVDmUcssx9H8o8RK4hF4BqHtWICIhgQOMBRARE0KJGC7w1MioRTiLvonkqJXkTnPgRCnPAzdM8ldLl79C0F/jVtXWVDkcC/T4xqSZOZ58uNfd9KvsUKiouigtNqjahT8M+lrnvYN9gGX0L56GqkHH06lkpJVMpVepA2BPFZ6U8JFMpd0kC/R7JM+RxMOOgJcAPZhzkSuEVAPzc/dAGaHmk4SNE1IogvGY4PhofG1cshJUZC+HPX8wBfmzTdVMpLxeflRIJLtL01Vok0K1AVVXO5JwpdeHOiawTmFQTAI2rNiamQYzlsvkQvxBcFPkjFk5Id6n0WSkFV8xTKQ0fgi4vmkNcplLuGQn0CihpGXvth5clLWOraKoQXjOcaG00EQERhAeE4+fuZ+OKhbhHVBUyjl09N/zsXlBNUKU2tHjMfFZKo4fAXd6B3g8S6Ne4dtm0koZSfRr2ITU3tdTo+/qWsd2Culnmvhv5N8LVxdXGRyLEPVRkME+lHNts/lDz8mnz9jpa6DbR/IFmYEuZSrEBCfRi1y+blpaXxms7X+PN3W+SZ8wDwMvNq1TLWG1NLVU9y9nJTQhHprsEx78vPivlB/NUiquHefTdabx5KsW/nq2rrPQk0IvN/33+DcummVQTJky83v51ImtF3rJlrBBORVUhM7n4A83NcHaPeSrFpxY0H2A+K6VRd5lKsTOSTsWuXx6thN6oZ0joEMttWTZNOK0iA/y56+pZKZf/MG+vEw5d/2EehdeVqRR7dttAVxRlGfAocEFV1bDibdWBVUAIcBoYrKpq2QtJOog6PnUsHQmv3y6Es/IrKjIvyXasZCol2zyV0rAbdBpXPJUSZOsyRTmVZ4S+HHgf+M81214BflBVdZaiKK8U355k/fLunwmtJsiyaaJSCCk08FB+Pt11+UTqC+DMmOKplH7FZ6V0B48qti5TVMBtA11V1Z8VRQm5bvMAoHvx958CP+HggS7Lpgln5aaqtNQX0F2XTzddPiFG87KCR901LPP3I27QN1C3lUylOIGKzqHXVlU1DUBV1TRFUWpZsSabkWXThLPwKzLRJT+fh3T5dMnPx8+kUgj86uXJ5/6+7PD2It3N/L9/XFAb2xYrrOaefyiqKEocEAdQv77jLboqhKNoYDDwkM48ldJSX4AbcNHFhR+8vfnJ24vdXp7kyyjcqVU00M8rihJYPDoPBC7cbEdVVZcASwDatGmjVvD5hKjcDq6GH6ZDdor5Q8roKdDiCTiz23KV5rcXzR/qHyueStnh7cUhD3dUafJ2g7XGTrxtjOUcNajLRSa6reIxt122LuuuVTTQ1wPDgVnFX9dZrSIhRGkHV8OG8WAwL2xC9ln45jlYPx6M+eDqDiFd+RcX2eHlRZpGzka+lbXGTrxqHEM+HgCkEsCrxjEADh/qt33/pSjKCmA38KCiKCmKoozCHOQxiqIcB2KKbwsh7oVtb1wN8xJqkXnBh8Gfwcun4KmvWennK2FeDnOMQyxhXiIfD942xtqoIuspz1kuQ2/yo2gr1yKEgPKvUm/Qweqnrt5u6NyfUYW8srFC9/tUoyVJbcARU32S1Aaco0aZ+91suyORf86FsBdZZ+HIeltX4XSGG14BoB4ZNHM5w3m1Gjnc2LKgLhfvd2lWJ4EuhC1d+sMc4knrIDXe1tU4pRWaN2nmcoaqirnJ3vVz6ABeFDDRbZWtSrQaCXQh7rfME5C01hzi6QfN2+q2hIenQrP+sLCVLauzG6qqkpqVz5G0HJLOXanw43R0PVLqdskHn3KWixCiYi4cNQd40jq4cNi8Lagd9PwXNOsH1RrYtj4bKzAWcfx8LklpVziSdoWkc+avV/Tmq1qtfeblY267nCLAryeBLsS9oKpw/vDVEM88BihQvyP0nm0O8UraP/xSXmGp0E5Ku8KJC7kYTebLVLw0roQG+vJoRF2aB/rRLNCP0Dq+tHhji40rt38S6EJYi6pCWsLVEL90ChQXaNAZ2o0xh7hv5eneaTKpnL6YZ54yScsuDvAc0q9cbYBX28+D5oF+9AitRfO65vAOqeGDq4tcDFUREuhC3A1VNX+YWTInnnUGFFfzSj6dJ8CDfaFKgK2rvOd0hUaOpueUGnUfS89BV2heqtHVRaFJrSp0fKCGZdTdLNCXGlU8bvPI4k5IoAtxp0wm82LIR9ZD0nq4kgIuGnggCh6aBA/2Ae/qtq7ynlBVlfNXCiyhXRLgf1zMQy1u7OHr6UbzQD8GtwmmeV0/mgf60aR2FTzcZK3de00CXYjyMBWZV/NJWgdHNkBuunkhiMbRED3ZvBCEl3OtL2soMnEyI7fUqPtIWg6X8got+9Sv7k2zQF8GRNajWaAvzev6Ua+qF4r0j7EJCXQhbqbICKf/Zw7xo99CXga4eUGTGPO6mk17gYevrau0imyd4eoZJsVfj5/PpbDIBIC7mwuhdXzp2bw2zQL9aF7X/EGlr6fGxpWLa0mgC3EtYyH8saM4xDdC/iXQ+JjDu/kAc5g78MLIqqqgGqpRpA/EVBCISV+XzrN+JDXraq+YmlXcaRbox4guITQPNE+ZNKzpg5urtN61dxLoQhj0cGp7cYh/Z15X093XvLJ98wHmaRWNl62rvGOqyQ1TQZ1rwjuQooJAMHkW72HCxT2D1g9U468dGhSfZeJLLV/PWz6usF8S6KJyMuTDiW3mED+2GQpzwNMfmj1qvlrzgShws58zMG7XvztD9cOY24Si4lG3SR+IqTAAS0NVFz2uHmlo/H/HxSMNV880XDzSUVyMLBh6yDYHJaxOAl1UHgW5cOJ7c4gnbwVDHnhVh7DHzSPxkG7g5m7rKm9QVv/uicZnWV/UEaPiRpKpAZlUhbPm/RW3y7h6puHml4iLxzlcPdNQNJdRFFlfxtlJoAvnpr8CyVvM54mf2AZGPfgEQESsOcQbdAFX+/zfQFXhT7U204xP39C/24CGH9VWtOA03V0P0Fz5k7fr6s3h7Zp/k0cUzs4+/5KFuBv5l83TKEnr4OQPUFQIvoHQarg5xOt3AJd7f050Rft3T3D9CwnqAxwwPUAWNz+LRgE2erxuuf2ej3P3Qxe3J4EunEPeRTi20Xyhz6mfwGQAvyBoO8Yc4kFtwUEWSF5Y9DhNlRR6u+4jQjnJe8aBXKDaDfvVJdMG1Ql7JoEuHFduBhzdYB6J//E/87JsVRtAh+eh+WNQr5X12/TdB4c8RuGjFFhue1HgtP27hXVJoAuHEmA0wq9LzSH+5y+gmqD6A9DlRfNIvI7W5iF+IUfPgbPZFb7/tWEOzt2/W1iXBLqwe7WNRmLydMTk5RNZUABn/wEBodBtojnEazW3WYjrCo0cSsnmQEoWCWezOHA2u9RFOtbirP27hXVJoAu7VM9g5GGdjpg8HREF5t4hx9w1fFDVn3FPboVaofe9piKTyvELOSScyeJAShb7z2SRfD6H4jbeBFf3olWDaozoHELL+lX5y4e773uNonKTQBd2o77BYBmJtyg0h/hhd3fmVfPnex9vzmjMfUPG3YcwV1WV9Ct6Es5kkZCSRcKZLA6lZlvawfp7aYgIrkrPFnWIDPYnIqiqtIIVNieBLmyqYaGBnnk6HtbpCC00AHDAw513qldlm7c3qZr78yeaozdwKCWb/WezOHDWPH1yIcc8l+3u6kKzuuZ2sBHB/kQGVyOkhrd0FBR2RwJd3BsHV8MP0yE7BfyDIHoKaAebr5a5kARJ6/gmJY3GBgMmIMHDg9nVq7LNx5t0t3v7Z2koMnEsPYeEsyXz3lmcyMi19PNuVNOHzo1rEhlclYjgqjQL9JVe3sIhSKAL6zu4GjaMN/dLAcg+C+vHweFvIDMZLp4AxYVLHhr+5VeNH7y9yLhHIa6qKimX80uNvBNTsykwmtvC1vBxJyK4Kv0i6hIZXBVtkD9Vve3v8n8hykMCXVScsRAKc83/FeRCYZ65ydXmV66GuWXfAjj2HTTqDh3/BqGPMmpNtNVLytYZSEi5Gt4HzmZxsXhBBg83F8Lr+fPXDg2IDK5KZHBVgqrJYgzCeUigVxaqag7ZGwI4Fwpyrn5/y5/lFf+seFtR4e2f93qnfjL/9+3foeHdXapeYCziSFoOCWcucyAlm4SzWfyRmQeYz2JsHFCFHqG1iCgO7wfr+KKRnt7CiUmgX+tm8762YCq6LkhzrgvV3Nv/7NpgLsw1X4RTHi4a8Khi7gnu7lP8vQ9UqQ3uVYpvV7n6/eZX7u1rgfnfI7WwJkX6IN5Yl0hCSjZHzl2xrKhTy9eDyOCqDGoTRGRQVcKD/GU1HVHpSKCXKGved8N48/flCHU3VcXbZMLHpOKjmvAu/uplUvExmcxXN14/NXGrADboyl+7xtscuNeGrXdN82XwNwRzyX8+N//ZnbaQvQeBbjL6YMoPpkgfTFF+MEX5QWDyBuC/mSmE1/NnRJcQIoOqElm/KnX8PGXqRFR6Euglfph+47yvIR82vgRnf7VMP/w77YIlsL1Vc4B7m0zcNgK/+4f5q+JSenRbEsR+QaVve/he/f76fUv9zOe+dA68l1STG0X6etcEeBCqoUbxT024eKSj8TuEi9dZXD3PcmjM97i6SHgLcb3KG+hT/cu3X0EO7Ftquent4U6e4sIFjYLOxQWdopDn4oLORSFPMX+9uv3qtk1D/mcOX42XzXuN3A83W2HHpCqcVOuy39QYfVoERfnBmArqAOZ/lBS3y7h6ncW12p7iAE9FcTGUemwJcyHKVnkDvYKeqlunYnesEmDdQuxYWSvs/MP4HIuMA0ijBrmYp064osfV6yzuVXbg6nUWF6+zuLjl2rByIRybBPo1jKoLborpprcdUZFJpcBYhN5gKvW1wGBCbyiiwHj1a+nvS+9ruc91XwuMJvQFcyhAg151pwAN2fhgXn7hKiNunKYOQ1y3E+FykkjlJI810siyaEJYkQR6sbXGTuw0hfGi29fUVS5yTq3BPOMTdHFJvOsud6pKmQFZ+nZZAXuTIDaaKDAUWb6W3KfQeOPjGE33PjAfcUnFk0I8XAx4Usjyol5l7mfElTc1yy23FUVW2BHCmiTQi71tjCWVAL4q7F5q+/emNqRSkwJVYx6FptcA1Q3VpDF/VTVgcgNVg6q6gan4q6pBLd6OquHBf26ucG1uLgoebi54aFzxLP7qcc1XPy8NAb4eeJZsd3O55ntXPDXXbNNcu+3qVw83F2Le+7lC9X3oPr/0a1bUmlRunIhaAbAAABqkSURBVGKqy8UKPb4QonwcPtArum7jac/St89Ro8z9sqnC28YhAHhQiCG7CEUxgIsBRTGCiwEUI4piRHHNN29TDCguRrjm+wmtx94QyteGs+d1X6/93s3BLoaZ6LZKVtgRwgYcPtCtpS4XyxxVBpLJdo//wwMDigLhFby68W9Rc++2RIchK+wIYRt3FeiKovQG5mM+5+wjVVVnWaUqG7jZqHKS20o8FcMt7inKIivsCHH/VTjQFUVxBRYBMUAKsE9RlPWqqiZZq7j7SUaVQghHdzcj9HbACVVVTwEoirISGAA4ZKCDjCqFEI7tbj5tqwecveZ2SvE2IYQQNnA3I/Syrr++4aRnRVHigLjim7mKohy7i+e0mttcPF4TyCz7R4kVe75nnPdydXktrUteT+tyktezQXl2uptATwGCr7kdBJy7fidVVZcAS+7iee47RVF+U1W1ja3rcAbyWlqXvJ7W5Wyv591MuewDmiiK0lBRFHdgCLDeOmUJIYS4UxUeoauqalQUZRywBfNpi8tUVT1stcqEEELckbs6D11V1e+A76xUiz1xqCkiOyevpXXJ62ldTvV6Kqoq3e6EEMIZOFaTECGEEDclgS6EEE5CAl0IIZyEBLoQDkBRlAaKojxc/L2Xoii+tq7JUTnzaymBXkxRlKaKovygKEpi8W2toij/tHVdjkheS+tSFGUM8BXw7+JNQcBa21XkuJz9tZRAv2op8CpgAFBV9SDmi6XEnZPX0rr+BnQGrgCoqnocqGXTihyXU7+WEuhXeauq+ut124w2qcTxyWtpXQWqqhaW3FAUxY0y+iaJcnHq11IC/apMRVEeoPiXqyjKQCDNtiU5LHktrWuH8v/t3Xm4XVV9xvHvCw0ECUlAAVGUQSLIVKYwRAqiEK2IyCQWUChi6wit1cehWhFUphaLY6siAhIU1ArUYkCk5mEmkyQiiBqoiooiYIKE8e0fex3P4XJJzg07Z2Wt8/s8T56bve894b1vuCv7rD0s6YPAWpL2Ay4GLsucqVRVdxk3FiWSNqe5a2wacB+wCDjK9p05c5UoumyXpNWANwPTaR4eOJNmhbD44R2j2ruMAX0ESWsDq9lenDtL6aLLEAYrBvRE0rtH2f0AMMf2/EHnKVl02S5JC3jqPO8DwGzgY7bvHXyqMtXeZQzoiaQZwC5059P2p3lE8FbAxbZPz5WtNNFluySdDjwOzEi7OlcM/RHY0/YBWYIVqPYuY0BPJM0EDrG9JG1PoLle9SCaI8utc+YrSXTZLknX2n7paPskLbC9Xa5spam9y7jKpeuFwCM9248Cm9h+CHg4T6RiRZftmiBpt86GpF2BCWkzLgcdm6q7fEbPQ6/MDOAGSZek7QOAC9OJvVvzxSpSdNmu44Avp3c6opkeOC71eUrWZOWpusuYcukhaReau8gEXGN7duZIxYou2ydpEs3P7P25s5Su1i5jQB9B0gbA+M627f/LGKdo0WV7JO0PbMOT+zwpX6Jy1dxlzKEnkl4r6Q6am2B+kD5enjdVmaLLdkn6D+Bw4F0073gOAzbJGqpQtXcZA3rXycDuwE9sbwbsC1ybN1Kxost2TbP9JuA+2x8F9gBekDlTqaruMgb0rkfTTQWrSVrN9tXADrlDFSq6bNfS9PFPkp5Hc9XQZhnzlKzqLuMql67705nvWcAFku6hgsuYMoku23WZpMnAGcBcmjsdv5g3UrGq7jJOiibpsqWHaN61HAlMAi4o/VbgHKLL9qSHSe1u+7q0vSYw3vYDeZOVZxi6jAEdkLQ6MNP2vrmzlC66bJ+k623vkTtHDWrvMubQAduP08ypTcqdpXTR5UpxhaRDJCl3kApU3WUcoSeSLqK5MuNK4MHOftvHZwtVqOiyXZIWA2vTPFTqIZrL7Wx7YtZgBaq9yzgp2vWd9Cs8c9Fli2xXsyp9brV3GUfoPSStBbzQ9u25s5QuumxPmh44EtjM9smSXgBsNMq6rWE5au8y5tATSQcA84Hvpu0dJF2aN1WZosvWfY7mBpgj0vYS4LP54hSt6i5jQO86EdgVuB8graxTzQ0HA3Yi0WWbdrP9DtJNMbbvA9bIG6lYVXcZA3rXY6NcjxrzUSsmumzXo+lyUANIWh94Im+kYlXdZQzoXQslHQGsLmmKpE8D1+UOVajosl2fAv4L2EDSx4FrgE/kjVSsqruMk6KJpGcB/wxMT7tm0iwau/TpXxVGE122T9JWwCtoLrO7yvaPM0cqVs1dxoCeSNrR9rzcOWoQXbZL0lnA1zu3rIcVV3uXMeXSdaak2ySdLGmb3GEKF122ay7wIUk/lXRGWg0qrJiqu4wj9B6Sngu8nuYB+BNp/iX/WN5UZYou2ydpPeAQ4A001/hPyRypWLV2GUfoPWz/xvangLfSXEf9L5kjFSu6XCm2ALYCNgVuyxuleFV2GUfoiaSX0BxNHgb8Hvga8E3b92QNVqDosl2STgMOBn4GXAR8q7bFjQel9i7jWS5d5wAXAvvZvjt3mMJFl+1aBOxh+/e5g1Sg6i7jCL2HpDWAF6fN220/mjNPyaLLdkl6LbBX2vyB7cty5ilZzV3GgJ5I2hs4D7iT5vrUFwBH256VM1eJost2STqF5lEKF6RdfwPMtv2BfKnKVHuXMaAnkuYAR3SeDijpxcCFtnfOm6w80WW7JN0C7GD7ibS9OjDP9vZ5k5Wn9i7jKpeucb2PerX9E2Bcxjwliy7bN7nn97Ea1DNTbZdxUrRrtqSzgfPT9pHAnIx5ShZdtusUYJ6kq2mmsPYCqpgiyKDqLmPKJUkrgL8D2JPmL3oW8DnbD2cNVqDosn2SNgKm0vR5o+3fZI5UrJq7jAE9kbQ2sDQtctyZW1vT9p/yJitPdNkuSQcB3+88kljSZOBltr+dN1l5au8y5tC7rgLW6tleC/hepiyliy7b9ZHe58unG2E+kjFPyaruMgb0rvG2l3Q20u+flTFPyaLLdo32cxrnv1ZM1V3GgN71oKSdOhuSdgYeypinZNFlu2ZLOlPSiyRtLumTxEnmFVV1lzGHnkiaSvPMkc6t6hsBh9uu5i97UKLLdqVzEh8G9k27rgA+bvvBfKnKVHuXMaD3kDQO2JLm7PdtvberS9rP9pXZwhUmuhwcSZ+2/a7cOWpQepcxoPdJ0lzbOy3/K8PyRJftij7bU3qXMYfeP+UOUJHoMoSVIAb0/sVbmfZElyGsBDGgh1C+eMfTnqK7jAG9f3fmDlCRO3MHqMxZuQNUpOguh/6kqKSJwIa270jbh9G9y3Gm7d9mC1eY6LJdkvYENrd9Xtr+BrBe+vTHbH8/W7jCDEuXMaBLXwCus/2VtP1T4HKagegx22/NGK8o0WW7JF0FvMv2rWl7AXAMsDbwQduvyhivKMPSZTW3vD4DU4G/79le3LkOVdI1eSIVK7ps18TOAJTc0bk5K628E/o3FF3GHDr8hZ/8NuWNPb+fPPKLwzJFl+16Ume2D+7Z3HDAWUo3FF3GgA5PSHpuZ8P2QgBJzweeyJaqTNFlu26TtP/InZJeA9w+yteHpzcUXcaUC5wBXCbpn4B5ad9OwL+mz4X+RZft+kfgO5IOBeamfTsD04DXZEtVpqHocuhPigJIehXwQWCbtGshcKrty/OlKlN02a60+tORdPv8ETDD9tJ8qco0DF3GgB5CAdJU1m40U1c317Rs2qDV3GXMoSfp2ciXSvqdpHskXSJp89y5ShRdtkvSm4EbgYOAQ4EbJB2bN1WZau8yjtATSTcAnwUuTLveQHPd6m75UpUpumyXpNuBabbvTdvPprnef8u8ycpTe5dxhN4l2+fbfiz9+irxEKkVFV2265fA4p7txcAvMmUpXdVdxhF6IulU4H6alXYMHA6sSXOkie0/5EtXluiyXZLOA7YDLqHp80DgJuAnALbPzJeuLLV3GQN6ImnRMj5t2zEH3Kfosl2Slrkqve2PDipL6WrvMgb0EFZRksYD69j+3Yj9GwIP1HS53co2LF0O/Ry6pKMkvXGU/W+RdESOTKWKLlv3KeCvRtm/L/DJAWcp3VB0OfRH6JLmAXvZXjxi/0Tgats750lWnuiyXZJutb3103zuR7a3Ge1z4amGpcuhP0IHVh85AAHY/iMwLkOekkWX7VrW6jnxszs2Q9FlNd/IMzBO0tojd0paB1gjQ56SRZftukfSriN3SpoK/G6Urw9Pbyi6jIdzwdnANyS9zfadAJI2pbnE7ux8sYoUXbbrvcBFkr4CzEn7dgHeRHOzVujfUHQ59HPoAJLeCnwAmEBzbeqDNA+U+nzWYAWKLtslaQPgHcC2adePgM/YvidfqjINQ5cxoPeQNIGmk6fMA0s62va5GWIVKbocHEnftH1I7hw1KL3LmEPvYXvJaANQcsJAwxQuuhyouFGrPUV3GQN6/5Z1ljyMTXTZrnib3Z6iu4wBvX9F/0WvYqLLEFaCGND7F0eV7Yku2xV9tqfoLmNA79+1uQNUJLpcQZLWlbT9iN3vyxKmcDV2GQN6IukESRPVOFvSXEnTO5+3/c6c+UoSXbZL0v+mPtcDfgicI+nPj3m1fUW+dGWpvcsY0LuOTbeoTwfWB/4WODVvpGJFl+2alPo8GDgnPRNn38yZSlV1lzGgd3Xmzl5N8xf9QwqfT8soumzXX0jaCHg98N+5wxSu6i5jQO+aI+kKmkFoZnr+yBOZM5UqumzXScBM4Ke2b04Lbt+ROVOpqu4y7hRNJK0G7AD83Pb9afHY59u+JXO04kSXIeQx9A/nkrTTiF2bSzE7sCKiy3ZJ+jTLuGbf9vEDjFO0Yely6Ad04N/Sx/HAzsAtNPO92wM3AntmylWi6LJds9PHlwJbA19P24fRfWJg6M9QdBlTLomkrwEft70gbW8LvMf2MVmDFSi6bJekq4Hpth9N2+OAK2zvkzdZeWrvMk6Kdm3VGYAAbC+kmQcOYxddtut5wDo92xPSvjB2VXcZUy5dP5b0JeCrNHNtRwE/zhupWNFlu04F5qWjS4C9gRPzxSla1V3GlEsiaTzwNmCvtGsW8HnbS/OlKlN02T5JzwV2S5s32v5Nzjwlq7nLGNBDKICkdYEpNCecAbA9K1+ictXcZUy5JJKmAKfQnAHv/Ysu+oH3OUSX7ZJ0HM2iIBsD84HdgeuBl+fMVaLau4yTol3nAJ8HHgP2Ac4Dzs+aqFzRZbtOAKYCd6WrMXakopXqB6zqLmNA71rL9lU001B32T6RSv7VziC6bNfSzvkHSWvavg3YMnOmUlXdZUy5dC1Nt6zfIemdwK+ADTJnKlV02a5fSpoMfBu4UtJ9wN2ZM5Wq6i7jpGgiaSrNpXWTgZOBicAZtm/IGqxA0eXKI2lvYBLwXduP5M5Tshq7jAF9BElr234wd44aRJftkbQnMMX2OZLWBybYXpQ7V4lq7jLm0BNJe0i6lXQDjKS/lPS5zLGKFF22S9JHaJZG+0DaNY7mpq0wRrV3GQN6178DrwTuBUiLMuy1zFeEpxNdtusg4LXAgwC27+bJt6+H/lXdZQzoPWz/YsSux7MEqUB02apH3MyNGpqprMx5SlZ1lzGgd/1C0jTAktaQ9B7i+SMrKrps10WS/hOYLOktwPeAL2bOVKqqu4yToomk5wBn0SwYK+AK4ATb92YNVqDosn2S9qNZdFvATNtXZo5UrJq7jAEdkLQ6cLztT+bOUrrosl2pz5m2q1mZPpdh6DKmXADbjwMH5s5Rg+iyXanPP0malDtL6Yahy7hTtOtaSZ+hWZrqz9dO256bL1Kxost2LQUWSLqSJ/dZxTqYA1Z1lzHlkvQ88L6XbcczSMYoumyXpKNH22/73EFnKV3tXcaA3idJR9fyl55bdNkuSd+0fUjuHDUovcuYQ+/fCbkDVCS6bFc8Z749RXcZA3r/lDtARaLLdsXb7PYU3WUM6P0r+i96FRNdhrASxIDevziqbE902a7osz1FdxkDev+uzR2gItFlu96XO0BFiu4yrnJJJL17lN0PAHNszx90npJFl+2QtIDRp6dEcxno9gOOVKxh6TIG9ETSDGAX4LK0a3/gZmAr4GLbp+fKVprosh2SNlnW523fNagspRuWLmNATyTNBA6xvSRtTwC+QfP85Dm2t86ZryTRZQh5xBx61wuB3nUFHwU2sf0Q8HCeSMWKLlskaXdJN0taIukRSY9L+mPuXCWqvct4lkvXDOAGSZek7QOAC9MD8G/NF6tI0WW7PgO8AbiYZirrTcAWWROVq+ouY8qlh6SdgT1pTpRcY3t25kjFii7bI2m27V0k3dI5eSfpOtvTcmcrTe1dxhF6Iuks4Ou2z8qdpXTRZev+JGkNYL6k04FfA1UtnTZAVXcZc+hdc4EPSfqppDMk7ZI7UMGiy3a9keZn9Z00j3x9AXBw1kTlqrrLmHIZQdJ6wCE082wvtD0lc6RiRZftkHTCyHc7o+0Ly1d7l3GE/lRb0FwvvSlwW94oxYsu2zHaM7yPGXSISlTdZcyhJ5JOo3nr9TPgIuBk2/fnTVWmEV1+nehyhUj6G+AIYDNJl/Z8aiIQC26PwbB0GQN61yJgD9u/zx2kAouAaTTPll4T2F4StmfljVWc62hO2j0H+Lee/YuBW7IkKtdQdBlz6D0krQtMAcZ39sUgNHaS3gIcD2wMzAd2B66PJehWnKQNgalp8ybb9+TMU7Kau4w59ETSccAsYCbw0fTxxJyZCnY8zQ/MXbb3AXYEfpc3UrkkHQbcBBwGvB64UdKheVOVqfYuY8ql6wSaQegG2/tI2opmYA9jt9T2UklIWtP2bZK2zB2qYB8CpnaOJCWtD3yP5vk4YWyq7jIG9K4YhNrzS0mTgW8DV0q6D7g7c6aSrTZiWuBe4t31iqq6yxjQu2IQaontg9JvT5R0NTAJ+G7GSKX7bnqC5YVp+3Dg8ox5Snb5KF3+T8Y8rYqToqOQtDdpELL9SNq3ru378iYLw0rSwXSfjTPL9n9ljlSkdEntjfR0Cexuu+iVijpiQO+TpLm2d8qdIwwfSaeNHHBG2xeWb7Sf494HdZWumrmjASh68dhQtP1G2ffXA09RMElvS8vQbSnplp5fi6joOvSYQ+9fvJUJAyXpbcDbgc0l9Q466xALbY/VDJrzDqcA7+/Zv9j2H/JEal9MufQpplzCoEmaBKzLcgahOL8TOoZ+QJe0me1FfXzdPNs7DiJTCGMRBxuhI+bQ0w0Fkq5azte9YgBZQlgRcX4nADGHDrCapI8AL5b07pGftH1m+ljNPFuoznC/zQ5/FkfozeILS2n+cVtnlF8hhFCEoZ9D75D017bj7ruwyojzO2Gshn5AH22apVdnyiWEQZM0x/bOkq6y/bTncCStF1OCAWIOHWJaJay64vxOGJOhH9BtxyNyw6rqDcDr6J7fCWGZhn7KpUPSeODNwDY8ecWiY7OFCoE4vxP6N/RH6D3Op1mZ/pXAScCRwI+zJgpDrXeaRdJLRn4+zu+EkWJA79rC9mGSDrR9rqQZNMvQhZBLTLOEMYkBvevR9PF+SdsCvwE2zRcnDLs4vxPGKgb0ri9IWpdmzcFLgQnAh/NGCiHO74T+xZ2iXefTPGN6T+Bc4LPAhlkThdA4H3guzfmdHwAbA4uzJgqrpBjQuy4BDgQeA5akXw9mTRRCYwvbHwYetH0usD+wXeZMYRUUUy5dG9t+Ve4QIYwizu+EvsQRetd1kuKoJ6yKRp7fuRU4LW+ksCoa+huL0jqDpnm3MgX4OfAwzTOmXcvisaFcktYEDqE5Kh+Xdtv2SdlChVVSTLnAa3IHCGE5LgEeAObQHGyEMKqhP0IPYVUnaaHtbXPnCKu+mEMPYdUX53dCX+IIPYRVVJzfCWMVA3oIqyhJmyzr87bvGlSWUIYY0EMIoRIxhx5CCJWIAT2EECoRA3oYKEmPS5ov6YeS5kqalvY/T9I3BpjjZZL+ezlf8yxJF0haIGmhpGskTVjOa/5X0i7PMNtkSW9/Jn9GGE5xY1EYtIds7wAg6ZXAKcDetu8GDh1EAEn9/n9/AvBb29ul121J97kqK9Nk4O3A5wbw3woViSP0kNNE4D4ASZtKWph+v42km9KR/C2SpqTP3ybpS+lo+QJJ+0q6VtIdknZNr91V0nWS5qWPW6b9x0i6WNJlwBW9ISRNTV+/+Yh8GwG/6mzYvt32w71Z0+vfI+nEntcdlf7bC3ty7Z2+n/npv7VO2v9eSTen77OzoMWpwIvS157xTEsOwyOO0MOgrSVpPs1CDRsBLx/la94KnGX7AklrAKvTPJt+C+Aw4O+Am4EjaJ5f/1rgg8DraNaF3cv2Y5L2BT5B8xwUgD2A7W3/QdLLANKUz6eBA23/34gcXwaukHQocBVwru07+vge17Y9TdJe6c/YFngP8A7b16Zpm6WSptNcX74rzbXll6bXvB/YtvNOJoR+xYAeBq13ymUP4Lz0SNhe1wP/LGlj4Fu275AEsMj2gvTaHwFX2Xa6AWfT9NpJwLmSptDclDOu58+90vYferZfAnwBmJ6mfJ7E9vx01D4d2Be4OWV+aDnf44Xp9bMkTZQ0GbgWOFPSBel7+mUa0KcD89LrJtAM8CP/YQmhLzHlErKxfT3wHGD9Eftn0Bx1PwTMlNQ5iu99MNUTPdtP0D04ORm4Oj375AB6lmzjqQuW/BpYCuy4jIxLbH/L9tuBrwKvplkEpfdnZ/zIlz31j/GpwHHAWsANkraiOSo/xfYO6dcWts9+uiwhLE8M6CGbNKitDtw7Yv/mwM9tf4rm+d9jucV9Et1572OW87X306z+84nOFMyIHC9NzyEnTf1sDdwF/BbYQNKz06NtRz6x8/D0mj2BB2w/IOlFthfYPg2YDWwFzASO7Vw5I+n5kjagWV5unTF8zyEAMeUSBq8zhw7NEerRth9PUyodh9OcWHyUZnWek2hOoPbjdJopl3cD31/eF9v+raQDgMslHWv7xp5Pvwj4vJpwqwHfAb6ZpnlOAm4EFtHM2/e6T9J1KXNnIed/kLQP8DjNAhWXpxOsLwGuT9//EuAo2z9LJ3sXpq97b5/fexhycet/CCFUIqZcQgihEjGghxBCJWJADyGESsSAHkIIlYgBPYQQKhEDegghVCIG9BBCqEQM6CGEUIn/BwHkMdZRQePJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sort data\n",
    "coverage_dataframe = coverage_dataframe.sort_values(by=['Percent Coverage'])\n",
    "\n",
    "\n",
    "exclude = ['Mean Coverage']\n",
    "coverage_dataframe = coverage_dataframe.loc[:, coverage_dataframe.columns.difference(exclude)] # Don't plot Mean Coverage\n",
    "ax = coverage_dataframe.plot(x='Bismark Subset',marker='o') # Line plot with marker points.\n",
    "coverage_dataframe.plot(kind='bar', x='Bismark Subset', ax=ax) # Overlays bar chart using axes defined in line plot\n",
    "\n",
    "# Save figure to file.\n",
    "# \"bbox_inches\" argument needed to prevent x asis lables from getting cut off.\n",
    "plt.savefig(output_plot, bbox_inches = \"tight\")"
   ]
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}