{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook 8: \n",
    "\n",
    "This is an IPython notebook. Most of the code is composed of bash scripts, indicated by %%bash at the top of the cell, otherwise it is IPython code. This notebook includes code to download, assemble and analyze a published RADseq data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "### Notebook 8\n",
    "### Data set 8: Barnacles\n",
    "### Authors: Herrera et al. 2015\n",
    "### Data Location: SRP051026"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the sequence data\n",
    "Sequence data for this study are archived on the NCBI sequence read archive (SRA). Below I read in SraRunTable.txt for this project which contains all of the information we need to download the data. \n",
    "\n",
    "+ Project SRA: SRP051026 \n",
    "+ BioProject ID: PRJNA269631\n",
    "+ SRA link: http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP051026"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## make a new directory for this analysis\n",
    "mkdir -p empirical_8/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For each ERS (individuals) get all of the ERR (sequence file accessions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    BioSample_s  MBases_l  MBytes_l                   Organism_s  \\\n",
      "0  SAMN03256053        39        32    Ashinkailepas seepiophila   \n",
      "1  SAMN03256054        87        72  Ashinkailepas kermadecensis   \n",
      "2  SAMN03256055        17        14             Leucolepas longa   \n",
      "3  SAMN03256056        78        65        Neobrachylepas relica   \n",
      "4  SAMN03256057         7         6         Neolepas sp. SH-2014   \n",
      "\n",
      "  ReleaseDate_s       Run_s SRA_Sample_s Sample_Name_s    collected_by_s  \\\n",
      "0  Dec 10, 2014  SRR1702480    SRS785645         AsOK3  HOV Shinkai 2000   \n",
      "1  Dec 11, 2014  SRR1702481    SRS785647        AsNiN2    ROV Quest 4000   \n",
      "2  Dec 11, 2014  SRR1702482    SRS785649         VuTO2       HROV Nereus   \n",
      "3  Dec 11, 2014  SRR1702484    SRS785650         bar06       ROV Jason 2   \n",
      "4  Dec 11, 2014  SRR1702485    SRS785651         NeIn2  HOV Shinkai 6500   \n",
      "\n",
      "   collection_date_s    ...     Platform_s SRA_Study_s         breed_s  \\\n",
      "0               1997    ...       ILLUMINA   SRP051026  not applicable   \n",
      "1               2012    ...       ILLUMINA   SRP051026  not applicable   \n",
      "2               2009    ...       ILLUMINA   SRP051026  not applicable   \n",
      "3               2009    ...       ILLUMINA   SRP051026  not applicable   \n",
      "4               2009    ...       ILLUMINA   SRP051026  not applicable   \n",
      "\n",
      "  g1k_analysis_group_s  g1k_pop_code_s identified_by_s       isolate_s  \\\n",
      "0       <not provided>  <not provided>      S. Herrera  not applicable   \n",
      "1       <not provided>  <not provided>      S. Herrera  not applicable   \n",
      "2       <not provided>  <not provided>      S. Herrera  not applicable   \n",
      "3       <not provided>  <not provided>      S. Herrera  not applicable   \n",
      "4       <not provided>  <not provided>      S. Herrera  not applicable   \n",
      "\n",
      "  isolation_source_s        source_s  tissue_s  \n",
      "0  Hydrothermal vent  <not provided>    muscle  \n",
      "1  Hydrothermal vent  <not provided>    muscle  \n",
      "2  Hydrothermal vent  <not provided>    muscle  \n",
      "3  Hydrothermal vent  <not provided>    muscle  \n",
      "4  Hydrothermal vent  <not provided>    muscle  \n",
      "\n",
      "[5 rows x 33 columns]\n"
     ]
    }
   ],
   "source": [
    "## IPython code\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import urllib2\n",
    "import os\n",
    "\n",
    "## open the SRA run table from github url\n",
    "url = \"https://raw.githubusercontent.com/\"+\\\n",
    "      \"dereneaton/RADmissing/master/empirical_8_SraRunTable.txt\"\n",
    "intable = urllib2.urlopen(url)\n",
    "indata = pd.read_table(intable, sep=\"\\t\")\n",
    "\n",
    "## print first few rows\n",
    "print indata.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def wget_download(SRR, outdir, outname):\n",
    "    \"\"\" Python function to get sra data from ncbi and write to\n",
    "    outdir with a new name using bash call wget \"\"\"\n",
    "    \n",
    "    ## get output name\n",
    "    output = os.path.join(outdir, outname+\".sra\")\n",
    "    \n",
    "    ## create a call string \n",
    "    call = \"wget -q -r -nH --cut-dirs=9 -O \"+output+\" \"+\\\n",
    "           \"ftp://ftp-trace.ncbi.nlm.nih.gov/\"+\\\n",
    "           \"sra/sra-instant/reads/ByRun/sra/SRR/\"+\\\n",
    "           \"{}/{}/{}.sra;\".format(SRR[:6], SRR, SRR)\n",
    "        \n",
    "    ## call bash script\n",
    "    ! $call "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we pass the SRR number and the sample name to the `wget_download` function so that the files are saved with their sample names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for ID, SRR in zip(indata.Sample_Name_s, indata.Run_s):\n",
    "    wget_download(SRR, \"empirical_8/fastq/\", ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 1012150 spots for empirical_8/fastq/166_40.sra\n",
      "Written 1012150 spots for empirical_8/fastq/166_40.sra\n",
      "Read 1581225 spots for empirical_8/fastq/72638_22.sra\n",
      "Written 1581225 spots for empirical_8/fastq/72638_22.sra\n",
      "Read 2024006 spots for empirical_8/fastq/82121_15.sra\n",
      "Written 2024006 spots for empirical_8/fastq/82121_15.sra\n",
      "Read 1020530 spots for empirical_8/fastq/AsNiN2.sra\n",
      "Written 1020530 spots for empirical_8/fastq/AsNiN2.sra\n",
      "Read 458811 spots for empirical_8/fastq/AsOK3.sra\n",
      "Written 458811 spots for empirical_8/fastq/AsOK3.sra\n",
      "Read 909109 spots for empirical_8/fastq/bar06.sra\n",
      "Written 909109 spots for empirical_8/fastq/bar06.sra\n",
      "Read 570988 spots for empirical_8/fastq/bar22.sra\n",
      "Written 570988 spots for empirical_8/fastq/bar22.sra\n",
      "Read 898436 spots for empirical_8/fastq/barJC6731B1_1.sra\n",
      "Written 898436 spots for empirical_8/fastq/barJC6731B1_1.sra\n",
      "Read 84629 spots for empirical_8/fastq/NeIn2.sra\n",
      "Written 84629 spots for empirical_8/fastq/NeIn2.sra\n",
      "Read 77842 spots for empirical_8/fastq/NeOg1.sra\n",
      "Written 77842 spots for empirical_8/fastq/NeOg1.sra\n",
      "Read 139189 spots for empirical_8/fastq/SEPR_3.sra\n",
      "Written 139189 spots for empirical_8/fastq/SEPR_3.sra\n",
      "Read 845743 spots for empirical_8/fastq/VuMaU1.sra\n",
      "Written 845743 spots for empirical_8/fastq/VuMaU1.sra\n",
      "Read 208597 spots for empirical_8/fastq/VuTO2.sra\n",
      "Written 208597 spots for empirical_8/fastq/VuTO2.sra\n",
      "Read 9831255 spots total\n",
      "Written 9831255 spots total\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## convert sra files to fastq using fastq-dump tool\n",
    "## output as gzipped into the fastq directory\n",
    "fastq-dump --gzip -O empirical_8/fastq/ empirical_8/fastq/*.sra\n",
    "\n",
    "## remove .sra files\n",
    "rm empirical_8/fastq/*.sra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 724M\n",
      "-rw-rw-r-- 1 deren deren  75M Nov 24 12:49 166_40.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 117M Nov 24 12:49 72638_22.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 151M Nov 24 12:49 82121_15.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  75M Nov 24 12:49 AsNiN2.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  34M Nov 24 12:49 AsOK3.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  67M Nov 24 12:49 bar06.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  42M Nov 24 12:49 bar22.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  66M Nov 24 12:49 barJC6731B1_1.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 6.3M Nov 24 12:49 NeIn2.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 5.8M Nov 24 12:49 NeOg1.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  11M Nov 24 12:49 SEPR_3.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  63M Nov 24 12:49 VuMaU1.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  15M Nov 24 12:49 VuTO2.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "ls -lh empirical_8/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a params file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyRAD 3.0.63\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "pyrad --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\tnew params.txt file created\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## remove old params file if it exists\n",
    "rm params.txt \n",
    "\n",
    "## create a new default params file\n",
    "pyrad -n "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note: \n",
    "The data here are from Illumina Casava <1.8, so the phred scores are offset by 64 instead of 33, so we use that in the params file below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## substitute new parameters into file\n",
    "sed -i '/## 1. /c\\empirical_8/           ## 1. working directory ' params.txt\n",
    "sed -i '/## 6. /c\\TGCAGG                 ## 6. cutters ' params.txt\n",
    "sed -i '/## 7. /c\\20                     ## 7. N processors      ' params.txt\n",
    "sed -i '/## 9. /c\\6                      ## 9. NQual             ' params.txt\n",
    "sed -i '/## 10./c\\.85                    ## 10. clust threshold  ' params.txt\n",
    "sed -i '/## 12./c\\4                      ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 13./c\\10                     ## 13. maxSH            ' params.txt\n",
    "sed -i '/## 14./c\\empirical_8_m4          ## 14. output name      ' params.txt\n",
    "sed -i '/## 18./c\\empirical_8/fastq/*.gz ## 18. data location    ' params.txt\n",
    "sed -i '/## 29./c\\2,2                    ## 29. trim overhang    ' params.txt\n",
    "sed -i '/## 30./c\\p,n,s                  ## 30. output formats   ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==** parameter inputs for pyRAD version 3.0.63  **======================== affected step ==\r\n",
      "empirical_8/           ## 1. working directory \r\n",
      "./*.fastq.gz              ## 2. Loc. of non-demultiplexed files (if not line 18)  (s1)\r\n",
      "./*.barcodes              ## 3. Loc. of barcode file (if not line 18)             (s1)\r\n",
      "vsearch                   ## 4. command (or path) to call vsearch (or usearch)    (s3,s6)\r\n",
      "muscle                    ## 5. command (or path) to call muscle                  (s3,s7)\r\n",
      "TGCAGG                 ## 6. cutters \r\n",
      "20                     ## 7. N processors      \r\n",
      "6                         ## 8. Mindepth: min coverage for a cluster              (s4,s5)\r\n",
      "6                      ## 9. NQual             \r\n",
      ".85                    ## 10. clust threshold  \r\n",
      "rad                       ## 11. Datatype: rad,gbs,pairgbs,pairddrad,(others:see docs)(all)\r\n",
      "4                      ## 12. MinCov           \r\n",
      "10                     ## 13. maxSH            \r\n",
      "empirical_8_m4          ## 14. output name      \r\n",
      "==== optional params below this line ===================================  affected step ==\r\n",
      "                       ## 15.opt.: select subset (prefix* only selector)            (s2-s7)\r\n",
      "                       ## 16.opt.: add-on (outgroup) taxa (list or prefix*)         (s6,s7)\r\n",
      "                       ## 17.opt.: exclude taxa (list or prefix*)                   (s7)\r\n",
      "empirical_8/fastq/*.gz ## 18. data location    \r\n",
      "                       ## 19.opt.: maxM: N mismatches in barcodes (def= 1)          (s1)\r\n",
      "                       ## 20.opt.: phred Qscore offset (def= 33)                    (s2)\r\n",
      "                       ## 21.opt.: filter: def=0=NQual 1=NQual+adapters. 2=strict   (s2)\r\n",
      "                       ## 22.opt.: a priori E,H (def= 0.001,0.01, if not estimated) (s5)\r\n",
      "                       ## 23.opt.: maxN: max Ns in a cons seq (def=5)               (s5)\r\n",
      "                       ## 24.opt.: maxH: max heterozyg. sites in cons seq (def=5)   (s5)\r\n",
      "                       ## 25.opt.: ploidy: max alleles in cons seq (def=2;see docs) (s4,s5)\r\n",
      "                       ## 26.opt.: maxSNPs: (def=100). Paired (def=100,100)         (s7)\r\n",
      "                       ## 27.opt.: maxIndels: within-clust,across-clust (def. 3,99) (s3,s7)\r\n",
      "                       ## 28.opt.: random number seed (def. 112233)              (s3,s6,s7)\r\n",
      "2,2                    ## 29. trim overhang    \r\n",
      "p,n,s                  ## 30. output formats   \r\n",
      "                       ## 31.opt.: maj. base call at depth>x<mindepth (def.x=mindepth) (s5)\r\n",
      "                       ## 32.opt.: keep trimmed reads (def=0). Enter min length.    (s2)\r\n",
      "                       ## 33.opt.: max stack size (int), def= max(500,mean+2*SD)    (s3)\r\n",
      "                       ## 34.opt.: minDerep: exclude dereps with <= N copies, def=1 (s3)\r\n",
      "                       ## 35.opt.: use hierarchical clustering (def.=0, 1=yes)      (s6)\r\n",
      "                       ## 36.opt.: repeat masking (def.=1='dust' method, 0=no)      (s3,s6)\r\n",
      "                       ## 37.opt.: vsearch max threads per job (def.=6; see docs)   (s3,s6)\r\n",
      "==== optional: list group/clade assignments below this line (see docs) ==================\r\n"
     ]
    }
   ],
   "source": [
    "cat params.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Assemble in _pyrad_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 234567 >> log.txt 2>&1  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "sed -i '/## 12./c\\2                    ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 14./c\\empirical_8_m2       ## 14. output name      ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 7 >> log.txt 2>&1  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "We are interested in the relationship between the amount of input (raw) data between any two samples, the average coverage they recover when clustered together, and the phylogenetic distances separating samples.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Raw data amounts\n",
    "The average number of raw reads per sample is 1.36M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count         13.000000\n",
      "mean      702926.615385\n",
      "std       544381.441228\n",
      "min        72418.000000\n",
      "25%       197210.000000\n",
      "50%       775823.000000\n",
      "75%       945017.000000\n",
      "max      1854154.000000\n",
      "Name: passed.total, dtype: float64\n",
      "\n",
      "most raw data in sample:\n",
      "2    82121_15\n",
      "Name: sample , dtype: object\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "## read in the data\n",
    "s2dat = pd.read_table(\"empirical_8/stats/s2.rawedit.txt\", header=0, nrows=42)\n",
    "\n",
    "## print summary stats\n",
    "print s2dat[\"passed.total\"].describe()\n",
    "\n",
    "## find which sample has the most raw data\n",
    "maxraw = s2dat[\"passed.total\"].max()\n",
    "print \"\\nmost raw data in sample:\"\n",
    "print s2dat['sample '][s2dat['passed.total']==maxraw]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at distributions of coverage\n",
    "pyrad v.3.0.63 outputs depth information for each sample which I read in here and plot. First let's ask which sample has the highest depth of coverage. The std of coverages is pretty low in this data set compared to several others. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "summary of means\n",
      "==================\n",
      "count    13.000000\n",
      "mean     33.165923\n",
      "std      17.043971\n",
      "min      12.552000\n",
      "25%      18.726000\n",
      "50%      35.771000\n",
      "75%      40.778000\n",
      "max      67.946000\n",
      "Name: dpt.me, dtype: float64\n",
      "\n",
      "summary of std\n",
      "==================\n",
      "count     13.000000\n",
      "mean     109.360538\n",
      "std       76.550954\n",
      "min       19.386000\n",
      "25%       57.572000\n",
      "50%       94.316000\n",
      "75%      154.275000\n",
      "max      256.698000\n",
      "Name: dpt.sd, dtype: float64\n",
      "\n",
      "summary of proportion lowdepth\n",
      "==================\n",
      "count    13.000000\n",
      "mean      0.241166\n",
      "std       0.100536\n",
      "min       0.147948\n",
      "25%       0.164743\n",
      "50%       0.182743\n",
      "75%       0.310341\n",
      "max       0.422289\n",
      "dtype: float64\n",
      "\n",
      "highest coverage in sample:\n",
      "2    82121_15\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "## read in the s3 results\n",
    "s8dat = pd.read_table(\"empirical_8/stats/s3.clusters.txt\", header=0, nrows=14)\n",
    "\n",
    "## print summary stats\n",
    "print \"summary of means\\n==================\"\n",
    "print s8dat['dpt.me'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of std\\n==================\"\n",
    "print s8dat['dpt.sd'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of proportion lowdepth\\n==================\"\n",
    "print pd.Series(1-s8dat['d>5.tot']/s8dat[\"total\"]).describe()\n",
    "\n",
    "## find which sample has the greatest depth of retained loci\n",
    "max_hiprop = (s8dat[\"d>5.tot\"]/s8dat[\"total\"]).max()\n",
    "print \"\\nhighest coverage in sample:\"\n",
    "print s8dat['taxa'][s8dat['d>5.tot']/s8dat[\"total\"]==max_hiprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "highest prop coverage in sample:\n",
      "2    82121_15\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "maxprop =(s8dat['d>5.tot']/s8dat['total']).max()\n",
    "print \"\\nhighest prop coverage in sample:\"\n",
    "print s8dat['taxa'][s8dat['d>5.tot']/s8dat['total']==maxprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Means for sample 82121_15\n",
      "67.9461337257 256.697997144\n",
      "79.3792210284 276.499000416\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "## print mean and std of coverage for the highest coverage sample\n",
    "with open(\"empirical_8/clust.85/82121_15.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "print \"Means for sample 82121_15\"\n",
    "print depths.mean(), depths.std()\n",
    "print depths[depths>5].mean(), depths[depths>5].std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the coverage for the sample with highest mean coverage\n",
    "Green shows the loci that were discarded and orange the loci that were retained. The majority of data were discarded for being too low of coverage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<toyplot.mark.BarMagnitudes at 0x7fe4346dd410>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<div align=\"center\" class=\"toyplot\" id=\"t9d3142543d40498bbbd98da7cf240c60\"><svg height=\"300.0px\" id=\"t3bbd5c03d5944190988b98dbd2283e08\" preserveAspectRatio=\"xMidyMid meet\" style=\"background-color:transparent;fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:12px;opacity:1.0;stroke:rgb(16.1%,15.3%,14.1%);stroke-opacity:1.0;stroke-width:1.0\" viewBox=\"0 0 350.0 300.0\" width=\"350.0px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:toyplot=\"http://www.sandia.gov/toyplot\"><g class=\"toyplot-axes-Cartesian\" id=\"t1fc4948f3a4f4062a8401066ebe07c68\"><clipPath id=\"t66fdc3eb2e9244e7afd483243dc22144\"><rect height=\"200.0\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect></clipPath><g class=\"toyplot-coordinate-events\" clip-path=\"url(#t66fdc3eb2e9244e7afd483243dc22144)\" style=\"cursor:crosshair\"><rect height=\"200.0\" style=\"pointer-events:all;visibility:hidden\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect><g class=\"toyplot-mark-BarMagnitudes\" id=\"t4469e76256d04277b67ffc595b9f375a\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"180.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"60.0\"></rect><rect class=\"toyplot-Datum\" height=\"80.161943319838059\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"159.83805668016194\"></rect><rect class=\"toyplot-Datum\" height=\"48.201272411798726\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"191.79872758820127\"></rect><rect class=\"toyplot-Datum\" height=\"35.396182764603822\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"204.60381723539618\"></rect><rect class=\"toyplot-Datum\" height=\"33.105841526894153\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"206.89415847310585\"></rect><rect class=\"toyplot-Datum\" height=\"26.963562753036427\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"213.03643724696357\"></rect><rect class=\"toyplot-Datum\" height=\"25.089647194910356\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"214.91035280508964\"></rect><rect class=\"toyplot-Datum\" height=\"21.654135338345839\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"218.34586466165416\"></rect><rect class=\"toyplot-Datum\" height=\"24.152689415847334\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"215.84731058415267\"></rect><rect class=\"toyplot-Datum\" height=\"22.695199537304802\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"217.3048004626952\"></rect><rect class=\"toyplot-Datum\" height=\"21.23770965876227\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"218.76229034123773\"></rect><rect class=\"toyplot-Datum\" height=\"25.089647194910356\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"214.91035280508964\"></rect><rect class=\"toyplot-Datum\" height=\"24.048582995951449\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"215.95141700404855\"></rect><rect class=\"toyplot-Datum\" height=\"24.777327935222644\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"215.22267206477736\"></rect><rect class=\"toyplot-Datum\" height=\"20.925390399074615\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"219.07460960092538\"></rect><rect class=\"toyplot-Datum\" height=\"28.525159051474873\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"211.47484094852513\"></rect><rect class=\"toyplot-Datum\" height=\"23.840370156159622\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"216.15962984384038\"></rect><rect class=\"toyplot-Datum\" height=\"20.717177559282817\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"219.28282244071718\"></rect><rect class=\"toyplot-Datum\" height=\"25.714285714285694\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"214.28571428571431\"></rect><rect class=\"toyplot-Datum\" height=\"23.840370156159622\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"216.15962984384038\"></rect><rect class=\"toyplot-Datum\" height=\"26.234817813765204\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"213.7651821862348\"></rect><rect class=\"toyplot-Datum\" height=\"24.048582995951449\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"215.95141700404855\"></rect><rect class=\"toyplot-Datum\" height=\"25.610179294389809\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"214.38982070561019\"></rect><rect class=\"toyplot-Datum\" height=\"22.486986697513004\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"217.513013302487\"></rect><rect class=\"toyplot-Datum\" height=\"21.862348178137665\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"218.13765182186233\"></rect><rect class=\"toyplot-Datum\" height=\"20.613071139386932\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"219.38692886061307\"></rect><rect class=\"toyplot-Datum\" height=\"22.278773857721234\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"217.72122614227877\"></rect><rect class=\"toyplot-Datum\" height=\"44.141122035858871\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"195.85887796414113\"></rect></g></g><g class=\"toyplot-mark-BarMagnitudes\" id=\"tc286bd492dd84271842cdf6ddcc183e3\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"26.963562753036427\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"213.03643724696357\"></rect><rect class=\"toyplot-Datum\" height=\"25.089647194910356\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"214.91035280508964\"></rect><rect class=\"toyplot-Datum\" height=\"21.654135338345839\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"218.34586466165416\"></rect><rect class=\"toyplot-Datum\" height=\"24.152689415847334\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"215.84731058415267\"></rect><rect class=\"toyplot-Datum\" height=\"22.695199537304802\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"217.3048004626952\"></rect><rect class=\"toyplot-Datum\" height=\"21.23770965876227\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"218.76229034123773\"></rect><rect class=\"toyplot-Datum\" height=\"25.089647194910356\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"214.91035280508964\"></rect><rect class=\"toyplot-Datum\" height=\"24.048582995951449\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"215.95141700404855\"></rect><rect class=\"toyplot-Datum\" height=\"24.777327935222644\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"215.22267206477736\"></rect><rect class=\"toyplot-Datum\" height=\"20.925390399074615\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"219.07460960092538\"></rect><rect class=\"toyplot-Datum\" height=\"28.525159051474873\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"211.47484094852513\"></rect><rect class=\"toyplot-Datum\" height=\"23.840370156159622\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"216.15962984384038\"></rect><rect class=\"toyplot-Datum\" height=\"20.717177559282817\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"219.28282244071718\"></rect><rect class=\"toyplot-Datum\" height=\"25.714285714285694\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"214.28571428571431\"></rect><rect class=\"toyplot-Datum\" height=\"23.840370156159622\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"216.15962984384038\"></rect><rect class=\"toyplot-Datum\" height=\"26.234817813765204\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"213.7651821862348\"></rect><rect class=\"toyplot-Datum\" height=\"24.048582995951449\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"215.95141700404855\"></rect><rect class=\"toyplot-Datum\" height=\"25.610179294389809\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"214.38982070561019\"></rect><rect class=\"toyplot-Datum\" height=\"22.486986697513004\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"217.513013302487\"></rect><rect class=\"toyplot-Datum\" height=\"21.862348178137665\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"218.13765182186233\"></rect><rect class=\"toyplot-Datum\" height=\"20.613071139386932\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"219.38692886061307\"></rect><rect class=\"toyplot-Datum\" height=\"22.278773857721234\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"217.72122614227877\"></rect><rect class=\"toyplot-Datum\" height=\"44.141122035858871\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"195.85887796414113\"></rect></g></g></g><g class=\"toyplot-coordinates\" style=\"visibility:hidden\"><rect height=\"14.0\" style=\"fill:rgb(100%,100%,100%);fill-opacity:1.0;opacity:0.75;stroke:none\" width=\"90.0\" x=\"200.0\" y=\"60.0\"></rect><text style=\"alignment-baseline:middle;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"245.0\" y=\"67.0\"></text></g><line style=\"\" x1=\"60.0\" x2=\"282.3333333333333\" y1=\"250.0\" y2=\"250.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"60.0\" y=\"250.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"136.66666666666666\" y=\"250.0\">10</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"213.33333333333331\" y=\"250.0\">20</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"290.0\" y=\"250.0\">30</text></g><text style=\"alignment-baseline:middle;baseline-shift:-200%;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"250.0\">Depth of coverage (N reads)</text><line style=\"\" x1=\"50.0\" x2=\"50.0\" y1=\"60.0\" y2=\"240.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 240.0)\" x=\"50.0\" y=\"240.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 187.94679005205322)\" x=\"50.0\" y=\"187.94679005205322\">500</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 135.89358010410643)\" x=\"50.0\" y=\"135.89358010410643\">1000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 83.84037015615962)\" x=\"50.0\" y=\"83.84037015615962\">1500</text></g><text style=\"alignment-baseline:middle;baseline-shift:200%;font-weight:bold;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 150.0)\" x=\"50.0\" y=\"150.0\">N loci</text><text style=\"alignment-baseline:middle;baseline-shift:100%;font-size:14px;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"50.0\">dataset8/sample=82121_15</text></g></svg><div class=\"toyplot-controls\"><ul class=\"toyplot-mark-popup\" onmouseleave=\"this.style.visibility='hidden'\" style=\"background:rgba(0%,0%,0%,0.75);border:0;border-radius:6px;color:white;cursor:default;list-style:none;margin:0;padding:5px;position:fixed;visibility:hidden\"><li class=\"toyplot-mark-popup-title\" style=\"color:lightgray;cursor:default;padding:5px;list-style:none;margin:0;\"></li><li class=\"toyplot-mark-popup-save-csv\" onmouseout=\"this.style.color='white';this.style.background='steelblue'\" onmouseover=\"this.style.color='steelblue';this.style.background='white'\" style=\"border-radius:3px;padding:5px;list-style:none;margin:0;\">Save as .csv</li></ul><script>\n",
       "(function()\n",
       "{\n",
       "  if(window.CSS !== undefined && window.CSS.supports !== undefined)\n",
       "  {\n",
       "    if(!window.CSS.supports(\"alignment-baseline\", \"middle\"))\n",
       "    {\n",
       "      var re = /\\s*alignment-baseline\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t9d3142543d40498bbbd98da7cf240c60 text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          if(match[1] == \"middle\")\n",
       "          {\n",
       "            var style = getComputedStyle(text[i]);\n",
       "            var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "            var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "            dy += 0.4 * font_size;\n",
       "            text[i].setAttribute(\"dy\", dy);\n",
       "          }\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "    if(!window.CSS.supports(\"baseline-shift\", \"0\"))\n",
       "    {\n",
       "      var re = /\\s*baseline-shift\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t9d3142543d40498bbbd98da7cf240c60 text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          var style = getComputedStyle(text[i]);\n",
       "          var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "          var percent = 0.01 * match[1].substr(0, match[1].length-1);\n",
       "          var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "          dy -= percent * font_size\n",
       "          text[i].setAttribute(\"dy\", dy);\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var data_tables = [{\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1729.0, 770.0, 463.0, 340.0, 318.0, 259.0, 241.0, 208.0, 232.0, 218.0, 204.0, 241.0, 231.0, 238.0, 201.0, 274.0, 229.0, 199.0, 247.0, 229.0, 252.0, 231.0, 246.0, 216.0, 210.0, 198.0, 214.0, 424.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"t4469e76256d04277b67ffc595b9f375a\", \"filename\": \"toyplot\"}, {\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 259.0, 241.0, 208.0, 232.0, 218.0, 204.0, 241.0, 231.0, 238.0, 201.0, 274.0, 229.0, 199.0, 247.0, 229.0, 252.0, 231.0, 246.0, 216.0, 210.0, 198.0, 214.0, 424.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"tc286bd492dd84271842cdf6ddcc183e3\", \"filename\": \"toyplot\"}];\n",
       "\n",
       "  function save_csv(data_table)\n",
       "  {\n",
       "    var uri = \"data:text/csv;charset=utf-8,\";\n",
       "    uri += data_table.names.join(\",\") + \"\\n\";\n",
       "    for(var i = 0; i != data_table.data[0].length; ++i)\n",
       "    {\n",
       "      for(var j = 0; j != data_table.data.length; ++j)\n",
       "      {\n",
       "        if(j)\n",
       "          uri += \",\";\n",
       "        uri += data_table.data[j][i];\n",
       "      }\n",
       "      uri += \"\\n\";\n",
       "    }\n",
       "    uri = encodeURI(uri);\n",
       "\n",
       "    var link = document.createElement(\"a\");\n",
       "    if(typeof link.download != \"undefined\")\n",
       "    {\n",
       "      link.href = uri;\n",
       "      link.style = \"visibility:hidden\";\n",
       "      link.download = data_table.filename + \".csv\";\n",
       "\n",
       "      document.body.appendChild(link);\n",
       "      link.click();\n",
       "      document.body.removeChild(link);\n",
       "    }\n",
       "    else\n",
       "    {\n",
       "      window.open(uri);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function open_popup(data_table)\n",
       "  {\n",
       "    return function(e)\n",
       "    {\n",
       "      var popup = document.querySelector(\"#t9d3142543d40498bbbd98da7cf240c60 .toyplot-mark-popup\");\n",
       "      popup.querySelector(\".toyplot-mark-popup-title\").innerHTML = data_table.title;\n",
       "      popup.querySelector(\".toyplot-mark-popup-save-csv\").onclick = function() { popup.style.visibility = \"hidden\"; save_csv(data_table); }\n",
       "      popup.style.left = (e.clientX - 50) + \"px\";\n",
       "      popup.style.top = (e.clientY - 20) + \"px\";\n",
       "      popup.style.visibility = \"visible\";\n",
       "      e.stopPropagation();\n",
       "      e.preventDefault();\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  for(var i = 0; i != data_tables.length; ++i)\n",
       "  {\n",
       "    var data_table = data_tables[i];\n",
       "    var event_target = document.querySelector(\"#\" + data_table.id);\n",
       "    event_target.oncontextmenu = open_popup(data_table);\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var axes = {\"t1fc4948f3a4f4062a8401066ebe07c68\": {\"x\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 30.0, \"min\": 0}, \"range\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 290.0, \"min\": 60.0}, \"scale\": \"linear\"}], \"y\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 1729.0, \"min\": 0.0}, \"range\": {\"bounds\": {\"max\": -Infinity, \"min\": Infinity}, \"max\": 60.0, \"min\": 240.0}, \"scale\": \"linear\"}]}};\n",
       "\n",
       "  function sign(x)\n",
       "  {\n",
       "    return x < 0 ? -1 : x > 0 ? 1 : 0;\n",
       "  }\n",
       "\n",
       "  function _mix(a, b, amount)\n",
       "  {\n",
       "    return ((1.0 - amount) * a) + (amount * b);\n",
       "  }\n",
       "\n",
       "  function _log(x, base)\n",
       "  {\n",
       "    return Math.log(Math.abs(x)) / Math.log(base);\n",
       "  }\n",
       "\n",
       "  function _in_range(a, x, b)\n",
       "  {\n",
       "    var left = Math.min(a, b);\n",
       "    var right = Math.max(a, b);\n",
       "    return left <= x && x <= right;\n",
       "  }\n",
       "\n",
       "  function to_domain(projection, range)\n",
       "  {\n",
       "    for(var i = 0; i != projection.length; ++i)\n",
       "    {\n",
       "      var segment = projection[i];\n",
       "      if(_in_range(segment.range.bounds.min, range, segment.range.bounds.max))\n",
       "      {\n",
       "        if(segment.scale == \"linear\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          return _mix(segment.domain.min, segment.domain.max, amount)\n",
       "        }\n",
       "        else if(segment.scale[0] == \"log\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          var base = segment.scale[1];\n",
       "          return sign(segment.domain.min) * Math.pow(base, _mix(_log(segment.domain.min, base), _log(segment.domain.max, base), amount));\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "\n",
       "  // Compute mouse coordinates relative to a DOM object, with thanks to d3js.org, where this code originated.\n",
       "  function d3_mousePoint(container, e)\n",
       "  {\n",
       "    if (e.changedTouches) e = e.changedTouches[0];\n",
       "    var svg = container.ownerSVGElement || container;\n",
       "    if (svg.createSVGPoint) {\n",
       "      var point = svg.createSVGPoint();\n",
       "      point.x = e.clientX, point.y = e.clientY;\n",
       "      point = point.matrixTransform(container.getScreenCTM().inverse());\n",
       "      return [point.x, point.y];\n",
       "    }\n",
       "    var rect = container.getBoundingClientRect();\n",
       "    return [e.clientX - rect.left - container.clientLeft, e.clientY - rect.top - container.clientTop];\n",
       "  };\n",
       "\n",
       "  function display_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var data = axes[dom_axes.id];\n",
       "\n",
       "    point = d3_mousePoint(e.target, e);\n",
       "    var x = Number(to_domain(data[\"x\"], point[0])).toFixed(2);\n",
       "    var y = Number(to_domain(data[\"y\"], point[1])).toFixed(2);\n",
       "\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "    {\n",
       "      coordinates[i].style.visibility = \"visible\";\n",
       "      coordinates[i].querySelector(\"text\").textContent = \"x=\" + x + \" y=\" + y;\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function clear_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "      coordinates[i].style.visibility = \"hidden\";\n",
       "  }\n",
       "\n",
       "  for(var axes_id in axes)\n",
       "  {\n",
       "    var event_target = document.querySelector(\"#\" + axes_id + \" .toyplot-coordinate-events\");\n",
       "    event_target.onmousemove = display_coordinates;\n",
       "    event_target.onmouseout = clear_coordinates;\n",
       "  }\n",
       "})();\n",
       "</script></div></div>"
      ],
      "text/plain": [
       "<toyplot.canvas.Canvas at 0x7fe4346ab310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import toyplot\n",
    "import toyplot.svg\n",
    "import numpy as np\n",
    "\n",
    "## read in the depth information for this sample\n",
    "with open(\"empirical_8/clust.85/82121_15.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "## make a barplot in Toyplot\n",
    "canvas = toyplot.Canvas(width=350, height=300)\n",
    "axes = canvas.axes(xlabel=\"Depth of coverage (N reads)\", \n",
    "                   ylabel=\"N loci\", \n",
    "                   label=\"dataset8/sample=82121_15\")\n",
    "\n",
    "## select the loci with depth > 5 (kept)\n",
    "keeps = depths[depths>5]\n",
    "\n",
    "## plot kept and discarded loci\n",
    "edat = np.histogram(depths, range(30)) # density=True)\n",
    "kdat = np.histogram(keeps, range(30)) #, density=True)\n",
    "axes.bars(edat)\n",
    "axes.bars(kdat)\n",
    "\n",
    "#toyplot.svg.render(canvas, \"empirical_8_depthplot.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print final stats table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\r\n",
      "15502       ## loci with > minsp containing data\r\n",
      "15502       ## loci with > minsp containing data & paralogs removed\r\n",
      "15502       ## loci with > minsp containing data & paralogs removed & final filtering\r\n",
      "\r\n",
      "## number of loci recovered in final data set for each taxon.\r\n",
      "taxon\tnloci\r\n",
      "166_40       \t12509\r\n",
      "72638_22     \t12551\r\n",
      "82121_15     \t12957\r\n",
      "AsNiN2       \t8260\r\n",
      "AsOK3        \t5942\r\n",
      "NeIn2        \t2368\r\n",
      "NeOg1        \t636\r\n",
      "SEPR_3       \t2401\r\n",
      "VuMaU1       \t10889\r\n",
      "VuTO2        \t4745\r\n",
      "bar06        \t2308\r\n",
      "bar22        \t9924\r\n",
      "barJC6731B1_1\t12071\r\n",
      "\r\n",
      "\r\n",
      "## nloci = number of loci with data for exactly ntaxa\r\n",
      "## ntotal = number of loci for which at least ntaxa have data\r\n",
      "ntaxa\tnloci\tsaved\tntotal\r\n",
      "1\t-\r\n",
      "2\t-\t\t-\r\n",
      "3\t-\t\t-\r\n",
      "4\t2793\t*\t15502\r\n",
      "5\t2943\t*\t12709\r\n",
      "6\t3110\t*\t9766\r\n",
      "7\t2763\t*\t6656\r\n",
      "8\t2107\t*\t3893\r\n",
      "9\t1187\t*\t1786\r\n",
      "10\t469\t*\t599\r\n",
      "11\t116\t*\t130\r\n",
      "12\t14\t*\t14\r\n",
      "\r\n",
      "\r\n",
      "## nvar = number of loci containing n variable sites (pis+autapomorphies).\r\n",
      "## sumvar = sum of variable sites (SNPs).\r\n",
      "## pis = number of loci containing n parsimony informative sites.\r\n",
      "## sumpis = sum of parsimony informative sites.\r\n",
      "\tnvar\tsumvar\tPIS\tsumPIS\r\n",
      "0\t356\t0\t3376\t0\r\n",
      "1\t692\t692\t3028\t3028\r\n",
      "2\t1043\t2778\t2719\t8466\r\n",
      "3\t1249\t6525\t1967\t14367\r\n",
      "4\t1419\t12201\t1495\t20347\r\n",
      "5\t1502\t19711\t1012\t25407\r\n",
      "6\t1441\t28357\t714\t29691\r\n",
      "7\t1356\t37849\t426\t32673\r\n",
      "8\t1225\t47649\t319\t35225\r\n",
      "9\t1065\t57234\t169\t36746\r\n",
      "10\t843\t65664\t120\t37946\r\n",
      "11\t722\t73606\t76\t38782\r\n",
      "12\t642\t81310\t47\t39346\r\n",
      "13\t511\t87953\t17\t39567\r\n",
      "14\t396\t93497\t7\t39665\r\n",
      "15\t302\t98027\t3\t39710\r\n",
      "16\t205\t101307\t7\t39822\r\n",
      "17\t152\t103891\t0\t39822\r\n",
      "18\t148\t106555\t0\t39822\r\n",
      "19\t74\t107961\t0\t39822\r\n",
      "20\t56\t109081\t0\t39822\r\n",
      "21\t31\t109732\t0\t39822\r\n",
      "22\t33\t110458\t0\t39822\r\n",
      "23\t17\t110849\t0\t39822\r\n",
      "24\t12\t111137\t0\t39822\r\n",
      "25\t5\t111262\t0\t39822\r\n",
      "26\t1\t111288\t0\t39822\r\n",
      "27\t2\t111342\t0\t39822\r\n",
      "28\t2\t111398\t0\t39822\r\n",
      "total var= 111398\r\n",
      "total pis= 39822\r\n",
      "sampled unlinked SNPs= 15146\r\n",
      "sampled unlinked bi-allelic SNPs= 14830\r\n"
     ]
    }
   ],
   "source": [
    "cat empirical_8/stats/empirical_8_m4.stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\r\n",
      "26325       ## loci with > minsp containing data\r\n",
      "26325       ## loci with > minsp containing data & paralogs removed\r\n",
      "26325       ## loci with > minsp containing data & paralogs removed & final filtering\r\n",
      "\r\n",
      "## number of loci recovered in final data set for each taxon.\r\n",
      "taxon\tnloci\r\n",
      "166_40       \t14910\r\n",
      "72638_22     \t16286\r\n",
      "82121_15     \t16879\r\n",
      "AsNiN2       \t12110\r\n",
      "AsOK3        \t9262\r\n",
      "NeIn2        \t2676\r\n",
      "NeOg1        \t832\r\n",
      "SEPR_3       \t2729\r\n",
      "VuMaU1       \t12796\r\n",
      "VuTO2        \t5385\r\n",
      "bar06        \t2969\r\n",
      "bar22        \t11523\r\n",
      "barJC6731B1_1\t14297\r\n",
      "\r\n",
      "\r\n",
      "## nloci = number of loci with data for exactly ntaxa\r\n",
      "## ntotal = number of loci for which at least ntaxa have data\r\n",
      "ntaxa\tnloci\tsaved\tntotal\r\n",
      "1\t-\r\n",
      "2\t7376\t*\t26325\r\n",
      "3\t3447\t*\t18949\r\n",
      "4\t2793\t*\t15502\r\n",
      "5\t2943\t*\t12709\r\n",
      "6\t3110\t*\t9766\r\n",
      "7\t2763\t*\t6656\r\n",
      "8\t2107\t*\t3893\r\n",
      "9\t1187\t*\t1786\r\n",
      "10\t469\t*\t599\r\n",
      "11\t116\t*\t130\r\n",
      "12\t14\t*\t14\r\n",
      "\r\n",
      "\r\n",
      "## nvar = number of loci containing n variable sites (pis+autapomorphies).\r\n",
      "## sumvar = sum of variable sites (SNPs).\r\n",
      "## pis = number of loci containing n parsimony informative sites.\r\n",
      "## sumpis = sum of parsimony informative sites.\r\n",
      "\tnvar\tsumvar\tPIS\tsumPIS\r\n",
      "0\t3730\t0\t13784\t0\r\n",
      "1\t2679\t2679\t3314\t3314\r\n",
      "2\t2320\t7319\t2799\t8912\r\n",
      "3\t2157\t13790\t1998\t14906\r\n",
      "4\t2056\t22014\t1510\t20946\r\n",
      "5\t2029\t32159\t1015\t26021\r\n",
      "6\t1857\t43301\t714\t30305\r\n",
      "7\t1653\t54872\t426\t33287\r\n",
      "8\t1542\t67208\t319\t35839\r\n",
      "9\t1320\t79088\t169\t37360\r\n",
      "10\t1061\t89698\t120\t38560\r\n",
      "11\t921\t99829\t76\t39396\r\n",
      "12\t805\t109489\t47\t39960\r\n",
      "13\t632\t117705\t17\t40181\r\n",
      "14\t457\t124103\t7\t40279\r\n",
      "15\t320\t128903\t3\t40324\r\n",
      "16\t221\t132439\t7\t40436\r\n",
      "17\t167\t135278\t0\t40436\r\n",
      "18\t154\t138050\t0\t40436\r\n",
      "19\t80\t139570\t0\t40436\r\n",
      "20\t58\t140730\t0\t40436\r\n",
      "21\t33\t141423\t0\t40436\r\n",
      "22\t34\t142171\t0\t40436\r\n",
      "23\t17\t142562\t0\t40436\r\n",
      "24\t12\t142850\t0\t40436\r\n",
      "25\t5\t142975\t0\t40436\r\n",
      "26\t1\t143001\t0\t40436\r\n",
      "27\t2\t143055\t0\t40436\r\n",
      "28\t2\t143111\t0\t40436\r\n",
      "total var= 143111\r\n",
      "total pis= 40436\r\n",
      "sampled unlinked SNPs= 22595\r\n",
      "sampled unlinked bi-allelic SNPs= 22162\r\n"
     ]
    }
   ],
   "source": [
    "cat empirical_8/stats/empirical_8_m2.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Infer ML phylogeny in _raxml_ as an unrooted tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_8/ \\\n",
    "                      -n empirical_8_m4 -s empirical_8/outfiles/empirical_8_m4.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_8/ \\\n",
    "                      -n empirical_8_m2 -s empirical_8/outfiles/empirical_8_m2.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "This is RAxML version 8.0.16 released by Alexandros Stamatakis on March 21 2014.\n",
      "\n",
      "With greatly appreciated code contributions by:\n",
      "Andre Aberer      (HITS)\n",
      "Simon Berger      (HITS)\n",
      "Alexey Kozlov     (HITS)\n",
      "Nick Pattengale   (Sandia)\n",
      "Wayne Pfeiffer    (SDSC)\n",
      "Akifumi S. Tanabe (NRIFS)\n",
      "David Dao         (KIT)\n",
      "Charlie Taylor    (UF)\n",
      "\n",
      "\n",
      "Alignment has 63494 distinct alignment patterns\n",
      "\n",
      "Proportion of gaps and completely undetermined characters in this alignment: 51.86%\n",
      "\n",
      "RAxML rapid bootstrapping and subsequent ML search\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "head -n 20 empirical_8/RAxML_info.empirical_8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the tree in R using `ape`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAMgCAIAAABUEpE/AAAgAElEQVR4nOzdeTyV6f8/8Ou2C0WW\nkCwpbSoNsqRiLGmSNJX2QkWloaimVdqmtO8lVOqDVmYqaaVJhRxqUMnIUkKLJWWJc879++P+fs7P\nR2i7OdLr+eiPc973dd/3dZ+Zc7zOfV/nviiapgkAAAAAsEdE2B0AAAAAaG8QsAAAAABYhoAFAAAA\nwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUA\nAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaA\nBQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABY\nhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAA\nAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAA\nAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQ\nsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAA\nyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAA\nAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIW\nAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAn1ZdXe3j46OlpSUr\nKzts2DAOh8PUf/rpJ1FRUbF6srKyCCF2dnYiIiJiYmIiIiIqKirBwcFM+8GDBzPtRUVF5eXlt23b\nJrRDAgCAloSABfBpixcvrqurS09Pz87OHjNmzMSJEwkhNTU16enpz58/59ajp6dH0zSHw0lPT+dy\nuTU1NevXr1+6dCkhpLa29p9//snPz2daRkRErFu3rsGOKioqZs+e3aVLF3V19Y+XAgDA9wIBC+AT\naJoOCwtbu3Ztx44dVVVVPTw8hg8fTghJS0tTUlJSV1dv0D43N5fL5fbu3ZsQIiEhYW9vLy0tTQjJ\nyMhQUlLS0NAghFAU1a1bNxUVlQbrTps2jcfjZWRk3L59e8OGDXw+vzWOEAAA2IaABfBpKioq8+bN\nu3fvHk3TsrKyR44cIYRwOJyffvrp48YcDoe5FEjTdH5+/vz586dPny6oM20qKip27tw5derU+ity\nuVw3N7f9+/fzeLzo6GhjY2MREbxDAQC+S2LC7gD7xo8fjz9LwBZfX18TE5OYmBg/Pz9ra2tFRcXg\n4GAbGxtCCIfDuXTpEkVRgsYPHz7s27cvh8O5ceMGU1dRUZkyZcratWuZ9n/++aegff/+/YOCgurv\nS0xMzMnJKTU11dDQkBBy+fLlVjtMAABgVzsMWMXFxcePHxd2L6A9oCiKuaKnq6sbFhZWXV39xx9/\nuLq6Pnv2jKIoDodz/vz50aNHN1iLw+GcOXNm/PjxH9cvX748YsQI5vGwYcP4fP7HXwYGDRpUVlZ2\n9OjR2bNnP3/+vMUODgAAWhBF07Sw+8AyKyuruLg4YfcC2gk+ny8hIfH+/XspKSlCyOvXrwcOHPji\nxYvq6mo5ObkXL16oqqo2aK+goJCamqqrq1u/XlNTw7QXjLuiKKq6uprZLGPq1KnBwcHMgK2cnBwb\nG5ucnJwWP0IAAGgBuJQG0BwRERFNTc0DBw68ffu2oKBg/vz5Xl5eFEU9ePBAXV29QboihGRnZxNC\ndHR0GtTT0tJUVVU/HtVe39WrVyMiIng83osXLzw9Pb28vNg9FgAAaDUIWACfcOLEiRMnTqirq9vY\n2Jibmy9ZsoQQwuFwCgoK6t8By97enqkbGBh8fOGv0RHx9cdvEUIOHjy4fv16ZWXlUaNGTZgwwdvb\nuyUPCwAAWhAuEQIAAACwDGewAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAA\nAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhY\nAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBl\nCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAA\ngGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsA\nAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELIBPWLFixcaNG4XdCwAA\n+J4gYAF8AofD6dev31esmJCQYG5u3rlz5+7du58/f76ZIgAAtDMIWADNoWn66wJWUVGRg4ODn5/f\nmzdvVq9e7eHh0VQRAADaHzFhd4B9dXV1OTk5wu4FtBMfPnyorKwMDw8PDAwUFxcPDAy0t7cnhKSn\np3t7e6elpcnKyk6dOpW5hrhgwQJFRcULFy5Mnz69a9euLi4uTGMDAwNRUVFCSHx8/MdFAABofyia\npoXdB5bp6ura2NgIuxfQTnTr1s3f3//o0aOTJk0KCQnZsGFDQUEBIcTKymrq1KkuLi6PHj0yNjau\nrq4WERExMjKSkZHZsWOHvr6+pKQks4Xs7OwZM2b89ttvkydPFmy20SIAALQb7TBgWVlZxcXFCbsX\n0E4sXbq0tLQ0ODiYEFJeXq6mplZZWSkiIlJRUfHgwYO8vLzY2NjMzMzExMSqqqrOnTs/ffq0a9eu\nzLqlpaV+fn6pqamHDh0aMGBAM0UAAGhnMAYLoDkcDsfCwoJ5nJSUZGRkJCIiEhISYm1tHRcXp6ys\nrKamNmTIEEJIamqqgYGBIF1FR0c7OTlZW1vfuXNHEKQaLQIAQPvTDsdgAbCFz+ffv3+/a9eukyZN\nysrK8vLyOnDgACFk6dKlsbGx/fv353A4p06dCggIIIQkJiaam5szK8bFxfn7+8fExEhJSVVWVlIU\nJSMj02hRmIcHAAAtBmewAJqUnZ1taGgoJyfXrVu3WbNm7dy509ramhDi4uLi4ODQt2/fv/76q3fv\n3seOHSOEJCYmmpmZMStu2bKFw+EoKyvLycnJyckxo9obLQIAQLuEMVgAAAAALMMZLAAAAACWIWAB\nAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYh\nYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAA\nliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAA\nAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQs\nAAAAAJYhYAEAAACwTEzYHWBTVVXVhw8f+Hy+sDsCAAAAP7R2FbC0evRU19EtfVks7I4AAADAD61d\nXSLs0k17zZFTnRSVhd0RAAAA+KG1q4AFAAAA0BYgYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAA\nYBkCFsBnqa2tNTc3z8zMFFTKy8snTpyopqZmamp65MiRZooN3Lx5U11dncfjMU/v3bs3aNAgZWXl\n48ePt/RRAABA60DAAvgELpe7fPlyExOTzMxMPT09psjn8x0dHUeMGPHixYtt27ZdvHixqWIDt27d\nmjhx4sCBA0VFRQkh7969Gz9+/NGjR8PCwhYvXtyaxwUAAC2nXd1oFKAl1NXVDRkyhMvlZmRkiIj8\n33eS2NhYSUlJNzc3QoiFhYWFhUVTxfru3LmzdetWJyenbt26MZXTp0+PHj3awMCgqqqqtLSUz+cL\ndgEAAN8vBCz4XlVVVdnb26uqqrboXnx9fU1MTBwcHG7dumVmZiaoh4aGamtr6+vrFxcXb9261dXV\ntamiQGJi4tatW0+ePPnzzz9PmjSJKZ4+fXrBggWEkOLiYk1NTaQrAID2AQELvlc1NTWysrKbN29u\nuV1QFKWhocE8TkhIWL16tWBRamrqwIED4+Lirl27Nm/evJkzZ4qIiDRaZNrfu3dv69at4eHhFEWl\np6cPHjyYqf/777/9+/cnhKSlpfXs2bPljgUAAFpTuwpYJcWFUUH7K0pLhN0RaCXS0tLdu3dvhR3V\n1tampKQIUhGfz8/JyYmLi1NWVrazs5OUlKQoqtGiYAsrV668fv16ZGQk81RWVrasrKxjx47Pnz/v\n0qULIeTy5cvOzs6tcCwAANAK2tX1iJCD+x1NBirKSgu7I9DePHjwQEdHR15ennlKUZS0tHRSUlJt\nba2fn5+rqytFUY0WBVu4du0aTdM0Te/atcvd3Z2maXl5eYqiOnbseOfOndTU1IsXLyJgAQC0G+0q\nYP3yyy8TJkwQ/BUEYMvdu3frD8CiKGr//v3u7u66urqysrLr1q1rqvixe/fuCc6EURS1ffv2SZMm\nubm5RUZGysnJtcKxAABAK6BomhZ2H1hmZWUVFxcn7F5AiystLZ0zZ865c+eE3REAAICG2tUZLAAA\nAIC2AAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFg\nAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACW\nIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAA\nAJYhYAEAAACwDAELAACgBXG5XGlp6UWLFgkqkZGRdnZ2zazi6Oh44MCBlu8atCAxYXcAoBE7duxI\nTExsvk1tbS2Hw3F2dm6qAUVRe/fuVVFRYbt3AABf4NGjRzwe7/jx45s2bZKSkiKEcDgcIyOjZlbh\ncDirVq1qamlCQoKvr29mZqa8vPyuXbscHR3Z7zR8MwQsaIsuX7482X+7qLh4881+bXbpf7ZtLC4u\nRsACAOHicDiOjo5paWmRkZFTpkxhKh4eHk21LywsfP369YABAxpdWlRU5ODgEBYWZmdnFxoa6uHh\ngYDVNiFgQRsl26mTqNgnAlbzJCQl2eoMAMBXS0lJMTY2NjU1PXz48JQpU2iaTklJYc5gxcfHr1q1\n6tGjR7KysmvXrp0xYwbTvn///lJSUg4ODurq6tXV1enp6SNGjAgICGBWcXFxsbe3J4QYGBiIiooK\n9+igKRiDBQAA0II4HI6hoaGLi0tiYuKTJ0/y8vIoitLU1ExISFiyZMmBAwdevnx57NgxDw8PPp9P\n/nsBkabppKSk/v37nzhxIiYmZufOncxSZ2fn7du3E0Kys7M9PT23bt0q5MODJiBgwXegrrZ2nrVJ\n+ZtXgkp5yWt/V+eZg/tsnu9S9a6imSIAgBDV1dWlpaUZGhoqKSn9+uuvQUFBzOkriqLWrFmTlJSk\nr68vKipqaWkpIyMjIiJCCGEa5OTkSEpKenp6EkJKS0v79OnDLGWeLliwYMaMGYcOHZo8ebIwDw+a\nhoAFbV308eDV08a+evG8fjF081otvT7Bt/+RlO5w9uDuZooAAEL08OFDdXV1BQUFQoiHh8exY8fu\n3LljaGhICElKSnr58iVN0zRNV1VVFRQUEEJommbOYCUlJVlaWjKhKiEhwczMjNlgdHS0k5OTtbX1\nnTt3mhqnBW0BAha0dVq9+oyfv7B+habpezeumNqNEpeQMLEdmRx3takiAIBwCYZbEUKGDRumpKQU\nGBjIVLS1tYOCgj58+JCcnNyzZ8/U1FRCSGFhYVlZmb6+fmJioqmpKbOi4HFcXJy/v39kZKStrW1l\nZWVlZaWQDgs+DYPcv1hpaenAgQN79+4t7I60Z6mpqXP++1jfZEiDpTVVlTVVlV279yCEaOj2LHv1\nsqkiAIBwMQOwmMcURc2dO3fRokV3794dO3bssWPHZs+eHRAQ0L1793379pmbmzPtBw4cKCEhkZSU\nxIx5J4QkJiYuXryYELJlyxYOh6OsrMzULSwsNm3atGnTpnv37klJSS1evNjb27uFDuTZs2daWlqC\nMfV6enqbN2/GDxibQ7c7lpaWLbr9N2/ejBs3rkV3Aba2tqcz8s9lFgr+EUJCbj9gHv8nJYsQcizx\n4bnMwh3nb0hKSTdaHDV91j///CPsQwGA9ubDhw9mZmaPHz8WVMrKypydnVVVVU1MTEJCQpop0jRt\na2sbFRXFPI6Li1NTU+NyuczTpKQkAwMDJSWl0NDQpvYeExPTt29f5vF//vMfPT29q1evVlVVJScn\nd+rU6e+//25q73v27GkQAB4/fvz27dtZs2apqKioqamtXbu2mcM8d+6csbEx8/jVq1eBgYGKiorN\nvCbNv2I/AlwihO+PVAcZqQ4yhblPCSFFeTmKqmpNFQEAWMTlcpcvX25iYpKZmamnp8cU+Xy+o6Pj\niBEjXrx4sW3btosXLzZVJP8dYtWvXz9CyK1btyZOnDhw4EDmtNC7d+/Gjx9/9OjRsLAw5nxVo5KT\nk42NjQkhr1+/XrVq1dWrV21tbaWlpY2MjObMmcPhcJrau7u7+7v/CggImDt3bq9evaZNm8bj8TIy\nMm7fvr1hwwbmh4qNHmb9U3HKysojRoxgBpY12rj5V+wHgUuE8P2hKMr4Z7uUv2/oGRjej48ztrZv\ntMit/SDsngJAu1JXVzdkyBAul5uRkSH4TV9sbKykpKSbmxshxMLCwsLCoqkiISQ3N7eysjI8PHzv\n3r1VVVWWlpbMlcH09HRnZ+fXr187OTk5OzuXlpby+XwvLy9FRcULFy5Mnz5dMNNOcnLyiBEjCCFb\ntmyZMWOGlpaWoHuCWzY0undJSUlJSUlCyNWrV+Pi4s6fP8/j8dzc3Ozs7CoqKqKjo42NjZmDavQw\nORzOhAkTCCG1tbVZWVl//PHHwYMHm2rc/Cv2g/ixjhbaDZdla57c57hbGlaUlkyY591MEQCALdLS\n0g4ODhRFCX7TRwgJDQ3V1tbW19dXUlI6evRoM0VCCIfD4fF4fD7f3Nx806ZNN27cGDJkCCHEy8tL\nVFQ0IiLi/Pnzu3bt0tTUFBERSUxMvHnzZlBQ0Pz585nVaZpOTk5mxsjfvHnz118bn8+iqb0TQh4/\nfuzn5xceHi4uLi4mJubk5JSZmammpubl5eXn59fUYTIn3tzd3SmKkpSU7N+/v5OTk42NTVOvSfOv\n2A8CAQvaqPdv376v+P//Qu89FpOQFDwVk5Bcsjd454U4z027eHx+o8XaDziDBQDsS0hIEPy+jxCS\nmppaWVkZFxe3Z8+ehQsXMlfZGi0SQjgczqhRozIyMk6ePDlx4sTa2lomLS1ZsuTly5c5OTk7duzQ\n0dHp2bNnVVVVRkZGeHg4c2GO2cKLFy/evHkzcOBAPp//8OHDHj16CLpRW1t78uTJ6urqZvbODMw6\nfvw4c3WPMWjQoLKysh07dsyePbupw8zNza2qqvrw4QNN01wuNyIiwtPTk6bppl6T5l+xHwQuEUJb\nZG1tHbNrffNtmMmembPrjZInRPBbGwAAVtTW1qakpAwePJh5yufzc3Jy4uLilJWV7ezsJCUlKYpq\ntMi053A4L1++fPToUWRkJFPp2LHjnj17/vjjj9LSUl1d3devX0tJSTk7O6emphoYGHTt2pWmaTU1\ntYcPH6qpqUVHR//0009SUlI0TYuLixcUFPTq1YvZzqFDh/7888+JEyc2s/e1a9dOnz5dMBZq6tSp\nwcHB0tLS8vLyY8aM2bt3b1OHyeFwfvrpJwkJCUKIqKiotbW1YJsfN27+FftxIGBBW/T7779/sk1p\naemcOXNOnz7dCv0BAGA8ePBAR0dHXl6eeUpRlLS0dFJS0ogRI/z8/FxdXZnk0WiRz+ffv3/fwcEh\nJSUlKyvLxsZm8ODBFy9eVFRUNDExSUhIeP78+YkTJ6qrq52dnQMDA5kvkDRNV1ZWZmdnV1RUbN++\nnRn/TlHUmDFjvL29Q0NDO3bsePbs2YCAgOvXrzez98zMzCtXrqSlpQmO5erVqxERETNnziwuLvb0\n9PTy8mrqMDkcjiAhlZeX+/r6Tp8+XZCxGjRu/hX7ceASIQAAwOe6e/du/eFEFEXt37/f3d1dV1dX\nVlZ23bp1TRUJIdnZ2YaGhnJyct26dZs1a1a/fv3Gjh1LCHFxcYmNjZWSklq6dGlFRUXv3r3l5OQS\nExOZHYmIiGzcuNHBwcHKymr69Olz5vzfXQJ37dqlqKjYp08fLS2tc+fOXbt2rU+fPs3snblLlri4\nuKDzBw8eXL9+vbKy8qhRoyZMmFD/HloNDpPD4ezZs0dMTExMTKx3797KysobN25sqnHzr9iPg6p/\nDbV9sLKyiouLa7ntl5SUeHh4nD17tuV2AZ+DOYN17ty5lt4Rl8uVk5ObO3fuzp07mUpkZOShQ4eu\nXm3yZvGOjo729vaCcakA0GatWLFCRkZm5cqVwu1GS//lgtaHM1gAn/Do0SMej3f8+PGamhqmwswU\n1swqzTdISEgwNzfv3Llz9+7dz58/z3J3AX5g1dXVPj4+WlpasrKyw4YNY24KRQj56aefREVFxerJ\nysoihNjZ2W3atMnPz09ERERFRSU4OJhpP3jwYKa9qKiovLz8tm3bPt5Xo2/kr3t319XV1T+xBO0D\nAhbAJ3A4HEdHR0VFRcGg1Po33PtYYWHh69evm5qEtaioyMHBwc/P782bN6tXr/bw8GiRTgP8kBYv\nXlxXV5eenp6dnT1mzJiJEycSQmpqatLT058/f86tR09Pj7n1QMeOHTMzM2tqatavX7906VJCSG1t\n7T///JOfn8+0jIiIEFxlE8jKyrKysnr8+LGUlNTAgQOZN/JXv7sLCgo0NDRYfSVA+DDIHX4ghYWF\ngrNQn0ldXT0lJcXY2NjU1PTw4cNTpkyhaVowe2t8fPyqVasePXokKyu7du1aZuKwlJSU/v37S0lJ\nOTg4qKurV1dXp6enjxgxIiAggFnFxcXF3t6eEGJgYCCY2AsAvhFN02FhYTk5OR07duzYsaOHh8fD\nhw8JIWlpaUpKSurq6g3a5+bm1tbW1tXVhYeHBwYGUhTF3AkzIyOjU6dOM2bMSEtLk5WVtbOzU1FR\nIYQsWLBAcNvPsLCwHj16xMXFvXv3rlevXl26dCHf8O7Oy8vT0tJ69eoVsyNoHxCw4AcyevTogQMH\nftGp+JkzZ3I4nLFjxxoYGKxaterJkycSEhIURWlqaiYkJCxZsiQkJKRPnz7x8fH29vbTpk0TERFh\nrg/SNJ2UlOTn5/fbb78VFRVpaWlt2rRJRETE2dnZ2dmZEJKdne3p6Sm48zIAfDsVFZV58+b5+voa\nGxvLysoeOXKE/PcWAx835nA4Ojo6jx8/1tHRiY+Pd3Jyys/PZ+o0TU+ZMuXq1avJyckWFhbMCK3E\nxEQZGZmgoKDevXvr6OgwN0A/duyYtLQ080b+ond3bW3tnDlzQkNDCSHMfufPn4/Rve2KMCZAbFmY\n7PkHUVJS8uuvv37RKmPGjCkvL/+iVWpra6WkpEpLS2manjx5sq+v75kzZ0aMGEHTtK2tbf23kmDe\n01GjRgUGBmZnZ3ft2pXH49E0nZGRMWDAgPo99/T0NDMzw1zUAOzKzs6eMmWKrKyslpbWtWvXmKKr\nq2uDP3wPHz6kaXrJkiWCCpPMpKSkeDye4Gd6jA4dOvB4vMrKSklJyYKCAsG+BGPSDx48KCh+0bvb\nyckpLy+Ppunly5f37dv35cuXLL8cIFQYgwXQnIcPH6qrqzN3Pfbw8Dh27NidO3eYAVhJSUmCD8Sq\nqqqCggLy3wkljIyMkpKSLC0tmSsOCQkJgl8pR0dHOzk5WVtb37lzp6lxWgDwdXR1dcPCwl69ejV9\n+nRXV1eapgkhHA7n/Pnz9f/y9e3bl6nr6+sfPXqUpumXL1+OGTPGyMhIRETkypUrenp6/v7+ly5d\ncnFxqa2t5fP5gtt+MjuKjo728/M7ceLE9u3bBTcs+NJ398yZM48fP04ICQ8PnzFjBq4PtjMIWADN\nEQy3IoQMGzZMSUkpMDCQqWhrawcFBX348CE5Oblnz56pqamEkMLCwrKyMn19/cTERMHUEILHcXFx\n/v7+kZGRtra2lZWVlZWVQjosgPaGz+eLiYkxgyylpaW9vLx4PB4hpKqq6uHDh8bGxh+3T0lJefbs\n2Y0bN2pqatLS0ry8vPz9/Wtqap49e3bo0KHVq1crKir+/fffzFD3xMREwbwRNjY2fn5+kZGRTk5O\ndnZ2YmJi5Kve3Q4ODjExMe/fv1dRUfmcuyuz6NSpUw1+6TxnzhzmFqYfc3FxYX56SVGU4MeYzFfK\nhISE8ePHq6ioaGpq+vj4VFVVtUbvvxMIWADNqf+DQYqi5s6dW11dzVSOHTsWGRmprKw8Z86cffv2\nMR++HA5n4MCBEhISSUlJHwesLVu2cDgcZWVlOTk5OTk5ZjAsAHw7ERERTU3NAwcOvH37tqCgYP78\n+V5eXhRFPXjwQF1dXVVVtUH77OxsPp9vZGQkuO3nzp07ra2t09LS5OTkZsyY0bdv37/++qt3795M\ne8FtPwkh8fHxqampzBu5f//+zInqRt/dPB7P1dW1oqKi0T6LiYmZmZkFBAQwI7daU48ePXJzcwVP\ns7Ozo6Kili9f3mjjY8eOcbncyspKCQkJwVOdG/YAACAASURBVI8xNTQ0wsLCXF1d3dzccnNzr1+/\nnpCQsGLFCmYV3IyGEIzB+nIYg9VGtM4YLAD4Xty+fdvAwKBDhw69evXasWMHl8ulaXr37t2EENF6\nmDGUYWFhw4YN+3gj+/fvd3R0rF8hhNTU1NSvnDlzRltbW0FBYeDAgSEhIXw+X7AoIyOjQePr16+b\nmpoyY60+9uDBAzU1tcLCwq896K9UVlZGCHn79i3zdMqUKRs2bGh+FQ6Ho6amJnj67NkzTU3NoqIi\nQSU1NVVVVZWm6cLCws6dO8fExPB4vCNHjjDFHxAC1hdDwGojELAAoK1Zv369paWljY3NmjVr4uLi\nqqqqaJpOS0szMjJKTU39uD2fz9fW1m71btI0TSsqKjIj8dPT09XV1d+/f5+bm9ujRw9m6b1790xM\nTOq3P3To0OjRowVP58yZs23btvoN3r9/Ly4uzuPxTp065ePjwxRTU1OZ+ap/QLhECAAAwI5Vq1bF\nxcXFxMSMHj06KSlp8uTJ5ubm+/btYy6lXb58uUH74uJiwViCVqarq5uXl0cI8fPzW7lypYyMTP1x\nZnfv3m3QsQY3WL59+/bIkSPrN8jNze3ZsydzM5rt27eTH/5mNLgPFgAAAJvExMQMDQ2ZOFJVVXX3\n7t1bt27JyspOmDDh119/DQwMlJKSYlrm5+draWkJpZPMMCwOh5OWlnby5ElCSIOAxUxELcDhcMaM\nGcM85vP5WVlZ2tra9RtcunRJsHppaamfn19qauqhQ4d+2J9L4wwWAABAS+nQoYONjc26devi4+Of\nPXtmbGxsaWl55MgRLpdLhBqwmDNYq1atWr9+vYSEBKkXsMrLy69du1b/DFZ1dXVGRobgDJaIiIiy\nsvK///4raFBUVLR7925mkDtuRsNAwAIAAGgNCgoKCxYsiI+P53K5Q4YMOXz4cE5OjhDPYEVFRb18\n+ZKZsZEQkp6eLiYmxuVyly5dKiEhUb9jaWlpKioqampqgsqkSZPmz5//9OnT9+/fR0dHDx061N/f\nX0dHBzejEUDAAgAAaD3i4uLu7u5XrlzJy8vbvn07c/aoFdA0HRUVxQxUJ4T06NEjPz+fmcKLaeDp\n6Tl8+HBLS8tOnTqZmppSFCVY9+MZ7jdu3GhsbDxs2DBdXd3g4OCgoCDmDvi4GY0AxbzQ7YmVlZVg\nBoOWUFJS4uHhgRmjhK60tHTOnDnnzp37/FWcnJxCQ0M7derUcr0CAPh8ISEhZWVlTd3hky00TZ89\ne3bXrl2Ojo4LFiyQkZFp0d0BA2ewAAAAhGPq1Kl//vlni+4iNjb2559/TklJOX/+/O+//4501Wrw\nK0IAAADhkJKS6t+//7179wYPHsz6xjkczurVqzU1NU+cOKGhocH69qF5CFgAAABC4+rqevToUXYD\nVlpamp+fn6qq6sGDBxvcTAFaDQIWAACA0AwePNjX17eqqqpDhw7fvrWcnJxly5aJiIhs3LixX79+\n375B+GoIWAAAAMI0ZsyYv/76a/Lkyd+ykeLi4jVr1uTl5a1du1ZYd4eH+hCwoJ1g5i5tXl1dXXl5\nOZ/Pb6pBx44dRUVFWe0XAMAnzJgxw9XV9asDVnl5+ebNmxMSElavXm1jY8Nu3+CrIWBBe/Dnn3/+\ntnhp5y6qn2zpOGV6U4uqKysnjLLfuHEjq10DAPgEFRUVCQmJZ8+eaWpqftGK1dXVe/bs+fPPP318\nfDZt2lT/zlUgdAhY0B5UVVU5zJxj6zztWzaS+yjjRXwMW10CAPh8M2fOPH78+KpVqz6z/YcPH3bt\n2hUVFeXr63vnzh3BzUKh7cB/EgAAACFzcHC4dOnS59z6m8fjHT58eNiwYQoKCrdu3ZowYQLSVduE\n/yoAAABCJiYmZmpqGh8f30wbmqbPnDkzbNiwsrKy2NhYd3f3VptmB74CAha0Q3W1tfOsTcrfvBJU\nykte+7s6zxzcZ/N8l6p3Fc0UAQCEws3N7ejRo00tFdyQ/cKFC7gh+3cBAQvam+jjwaunjX314nn9\nYujmtVp6fYJv/yMp3eHswd3NFAEAhEJfXz8nJ+fdu3cN6hwOx8HB4a+//vrPf/6zefPmzp07C6V7\n8KUQsOB7RdP0s2fPDh8+fPjw4djYWEFdq1ef8fMXNmh578YVU7tR4hISJrYjk+OuNlUEABAiZ2fn\ns2fPCp6mpaXZ29sHBQXt27dv9+7dXbt2FWLf4EvhV4TwverUqdMvv/zycV3fZEiDSk1VZU1VZdfu\nPQghGro9y169bKoIACBEkydPnjRpkqurq+CG7Nu3b8cN2b9TOIMF3ysxMbG1a9e6u7u7u7v//PPP\nn2zP3CGGpmk+j9d8EQBAKDp37iwmJubm5jZr1ix3d/eTJ0+2aLo6deqUkZFR/cqcOXMWL17caOPa\n2lpJSUkHB4f6RRcXF4qiXr5s8guqlpZWcnKy4On+/fvHjBnzbb3+biBgQfsn1UFGqoNMYe5TQkhR\nXo6iqlpTRQAA4erbt29dXV1sbGwr3JO9R48eubm5gqfZ2dlRUVHLly9vtHFGRoaIiEh+fr6g8ujR\no7Nnz3br1q1Lly6NrlJSUlJcXDxgwABB5d69e8bGxoKntbW15ubmmZmZ33okbRICFrR/FEUZ/2yX\n8vcNmqbvx8cZW9s3VQQAEK5ff/1VU1Ozde7JrqurW1paWlHxf7+hXrNmzaJFixQVFRttzOFw7Ozs\nnj17JrhZFzMzT4NzYPXdv39/wIABkpKSgkpycjITsLhc7vLly01MTDIzM/X09Fg7pLYEAQt+CC7L\n1jy5z3G3NKwoLZkwz7uZIgCAEGlqaubl5bXOvuTl5RUVFZndZWRk3Lx5c+HChXl5eT179mQaJCcn\nC+aNTklJsbCwkJaWLi8vJ4RwOJznz59raGgwASs9Pf3nn39WUlLS1tZeuXIls8r9+/frn6969+7d\nkydPmPZ1dXVDhgyxsbExMTFprzdKxSB3aCcexN98//at4Ok035VxUWfqNzCwsDSwsCSEXI443mix\n7FWxnlLH1uktAECj1NTUioqKWm13urq6eXl5AwYM8PPzW7lypYyMzIULF8zNzZmld+/eFQQsDocz\nadIkLS2t/Px8BQWFFStWbNy40c/Pz9HRkRDi5eU1derUq1evPnr0yNjYeP369SIiIvfv37e1tRXs\nKyUlRVtbmzlDJi0t7eDgcOvWLTMzs1Y72FaGgAXtgZ2dnbi4+Cebbd261dPTs0OHDk0sH2hoaMhu\nxwAAvoioqGheXt67d+/k5ORaYXfMMCwOh5OWlnby5ElCSGJiYv2ANXbsWELIhw8fHj58+NNPPzEB\nq7y8nMvlWlpaPnjwgPnY/Ouvvx48eBAeHh4bGzto0CDmpNSDBw9WrFgh2Jfg+qBAQkLC6tWrW+Ew\nhQIBC9oDJSWlCRMmfLJZWFiYk5NTp06dWqFLAABfITY2VkVFxc7OztnZed68eVJSUi26O+YMVkxM\nzPr165mJdxITE2fNmkUIKS8vv3btWkBAACEkPT1dS0urU6dOTMA6efLk9u3bHz16pKampqioGBIS\ncujQIQcHh8GDBzMVQgifz8/OztbV1RXs6/LlyzNmzBA8ra2tTUlJGTx4cIseoBC1zwufAAAA352S\nkhJ/f/8bN27cuXNHQ0PDysoqICCgpqam5fbYo0ePqKioly9fTpw4kamkp6eLiYlxudylS5dKSEho\naWkRQlJSUpixU1paWocPH1ZWVjYzM+NwOExx6dKlwcHBq1evVlRUPHXqFHNVUURERENDY/v27a9f\nvy4rK1u9enVpaamzs7Ng1w8ePNDR0ZGXl2+5oxMunMECAABoE+bNm7dp0yZmnsEJEyY4ODgcOHDA\nysrK0dGxwcW1b6epqamnp9ejR4/8/PxDhw4JRpp7enoOHz5cT0/PzMzM1NSU+T0jh8NhLgVqa2s/\nevQoPDycKTIBy8XFxcHBQUZGZty4cb179z527BhzSeH06dO//fbbpk2blJSUHB0dr127Ji0tLejA\n3bt32/EALEIIJfi9ZbthZWUVFxfXctsvKSnx8PCoP5sBfC+cnJxCQ0NxiRAA2qCIiIjU1NStW7c2\nqL97927kyJGVlZW2trYs/uBuwIABU6ZMYWtr8DGcwQIAABCy58+fHzx48Nq1ax8vkpOTu3379rlz\n544ePRoWFoaviN8LBCwAAABhoml67ty5e/bsqX9PzgbGjRvXuXNnBweHiIgIDQ2N1uwefB0McgcA\nABCmAwcOWFhYGBgYNN/Mysrq4MGDEyZMSEtLa52OwbdAwAIAABCaJ0+eXLx48ffff/+cxvr6+ufO\nnfP29r5161ZLdwy+ES4RAgAACEddXd3cuXOPHDny+aPX1dXVo6KiJk6cWFhYOGnSpBbtHnwLnMEC\nAAAQjs2bN0+aNElHR+eL1pKXlz9//vz58+eZu4BC24SABQAAIATJyckpKSkeHh5fsa6kpGRYWFhp\naam3tzefz2e9b/DtELAAAABaW3V1tY+Pz8GDB796CxRFBQQEdO/efeLEidXV1Sz2DViBgAUAANDa\nVqxY4ePjo6am9o3b8fb2njRp0qhRo968ecNKx4AtGOQOAADQqmJjY0tKSsaOHcvK1phbZDk6OkZE\nRDBTB0JbgDNYAJ+wYsWKjRs3fvXqb968ERUVfffuHfP09u3bo0aNUlZW7tat2+7du1nqYyOePXtG\nUZTYf/Xt2/f8+fMttzsA+EzMzMd79+5lcZtWVlZBQUHOzs73799ncbPwLRCwAD6Bw+H069fvW1bv\n3bu3nJwcISQsLGzWrFkLFy589uxZVFTUmjVrBDezKS8vnzhxopqamqmp6ZEjRwghe/fupf5XZmZm\nRUXF7Nmzu3Tpoq6uvm7duvo7qq2tNTc3z8zMFOzX2NiYy+VyudyioqKFCxe6ubk11biB5pcCwLfw\n8fFZt24d6zPe9OvXLyoq6rfffouJiflk47dv33p7e2tpacnLy69Zs0ZQT0hIMDc379y5c/fu3et/\nJfv4M+HzW9bXYClN0/r6+oKPOFdX1y896jaNbncsLS1bdPtv3rwZN25ci+4CWsiYMWPKy8u/aBU+\nn6+goJCVlfXVO123bt3MmTNpmn716pW2tnZeXp5g0eLFi7dv307TNI/HGzp0aEhICI/Hi4+PHzt2\nLE3TNTU17/4rICBg7ty5fD5/9OjRLi4ur169evr0qbi4OI/Ho2m6rq5u2bJlBgYGCgoKTIWm6eXL\nl8+dO1ewr7y8vB49ejTVWKD5pQDwjU6fPv3bb7+13PbLysqsra2DgoKabzZq1KgNGza8fv368ePH\ngk+SwsLCzp07x8TE8Hi8I0eOqKqq0k18Jnx+S4FGl8bHxw8bNkzwQVdTU8PuqyFcGIMF0Jzc3NzK\nysrw8PDAwEBxcfHAwEB7e3tCSHp6ure3d1pamqys7NSpU5lriAsWLFBUVLxw4cL06dMXLVrEbCE5\nOXnEiBGEkC1btsyYMaP+CImtW7cyD2JjYyUlJZkzTBYWFhYWFoQQSUlJZmKyq1evxsXFnT9/nsfj\nubm52dnZVVRUREdHGxsbMzcnrKurGzJkCJfLzcjIENyukMPhTJgwgRBSW1ublZX1xx9/ML9XarSx\nQPNLAeBbFBcX79+///Llyy23C3l5+ejoaFdX14KCAn9//6aa3blzZ/bs2aKiomfOnDEyMmLe7PHx\n8S4uLsxHnIGBgaioKGniM+HzWwo0ujQwMLCysrJr166dO3c+cODAyJEjWX9BhAgBC34gBQUFixYt\nEhcX//xVtLS0eDxe9+7d8/PzQ0JCZs+eXVBQQAjx8vKaOnXq1atXHz16ZGxsvH79ehERkcTERBkZ\nmaCgIH19fWZ1mqaTk5NXrlxJCLl582ZwcHCjewkNDdXW1tbX1y8uLt66dWv98+SPHz/28/OLiYlh\nuu3k5JSammpoaEgIEXxMS0tLOzg43Lp1y8zMTLBfDodz7do1d3d3pnLq1CkbG5tGG9fX/FIA+Go0\nTc+ZM2fHjh1SUlItuiNJSckTJ054eno6OzufOHGi0Qmk3dzcBEPsb9y4wTxwdnZ2dnYmhGRnZ3t6\nejLfABv9TPj8lgIfL62qqsrOzl65cuXw4cN37949e/bsFy9esPg6CJ+Qz6C1AFwihKZkZWU9/UIL\nFy6cNWsWs3pZWZmUlBRzcvvt27d///13aGjozJkzTUxMaJqurKyUlJQsKCigabqmpoZp9vz5czEx\nserqah6PJy0t/f79e0FnPnz4EBERUVVVRdN03759J0+e/OrVq7CwsI4dOwrOn5eWlurr6z958qT+\nUfD5/LKysh07dmhoaNSvW1hYXLlyhXn89OlTSUnJDx8+0DTN5XIjIiKUlJT4fH6jjT/W/FIA+ApB\nQUFr1qxpzT2amJhoa2uHhoZyudz69f379//888+ZmZmvX7/euXNn7969BYtKSko8PT3NzMz++eef\n+qt8/Jnw+S0/Z2lpaam0tHQ7G5aAgPXFELB+KFZWVkePHmUeX7582cLCgqbp4OBgIyMjf3//S5cu\nLVu2zMfHh6bp+Ph4Jmkxw7YKCwtpmj506NDgwYOZYseOHTMzMwVb3r17t5WVFZ/P5/F4UlJSL1++\npGn69evXysrKgiTk7e0dEBAgWGXKlClMIKNp+unTpzo6OoJFHz58kJaWLisrY56eOnXKzMxMsPTV\nq1f1N9ugcQPNLwWAr5CVlfXzzz/X1dW15k7r6uqGDh26Zs2awYMHBwYGCmJWnz59srOzmce5ubla\nWlrM44sXLw4dOjQyMrL+lzG6sc+Ez2/ZzNLIyMinT58yj+/duzdkyJBvPN62BgHriyFg/Th4PJ68\nvPy0adOqq6v/+ecfPT2969ev0zTduXPnBw8e8Hi8pKQkHR2d06dP0zS9devWRYsWMWtJSEjcunUr\nMzOzZ8+egYGBzNamT58+YsSI4uLiqqqq48ePq6urP3r0iP5vIDt//vyHDx/mzZu3dOlSpv3jx497\n9+5dW1sr6I+SklJISAiXyy0oKLC3t9+5c6dgUVJSUt++fQVPlyxZ4u3tzTwuKyubPn06kwIbbdxA\n80sB4EtxuVwbG5vHjx+3/q6zs7OHDRtWWlq6efNmc3Nz5sNKRUXl6NGjPB4vNzfX1tZ23759NE3H\nxsYaGRm9fv2aGW9e/3R7g8+Ez2/ZQIOlzs7OU6dOff/+fX5+vqmpaWxsLLvHLnQIWF8MAevH8eTJ\nE2tr63nz5ikpKRkZGUVHRzN1Hx8fDQ2NXr16rVixYuTIkb/88gtN0+PGjWM+vGia3rp1a8eOHdXU\n1NatWyf4hldSUjJlyhQFBQVlZeUxY8Y8fPhQsKPw8HBVVVUNDY0lS5YIfkczatSogwcP1u/PmTNn\ntLW1FRQUBg4cGBISUv+7486dOwWXMmmatrKyoihKVFRUVFS0S5cuPj4+1dXVTTVuoPmlAPClAgIC\ndu/eLay979u3j7k0+eLFCw8Pj19++eX48eM9evTo3LmzoaFhREQE80nCjFgXYM7WMxp8Jnx+ywYa\nLH3y5ImJiUmnTp0MDQ0vX77M6kG3CRRN06076KvFWVlZxcXFtdz2S0pKPDw8zp4923K7AACA9uGf\nf/5Zs2ZNVFQURVFC6QBN005OTitWrDAxMSGEFBQUKCoqSktLC6UzPxT8DBsAAKBF1NTULFiwYP/+\n/cJKV4QQiqIOHjzo6+tbVVVFCNHQ0EC6ah0IWAAAAC3C39/fw8Oja9euwu2Gurq6r6/v8uXLhduN\nHw0CFgAAAPv+/vvv/Pz8adOmCbEPghm3fv311z179jAzbpEmJrr5eMKuRufmamqaHYFGN97U1Drt\nmbAHgbEPg9wBAEC4KioqzM3N37x5I9xu1J9xa9OmTR4eHnw+v9GJbhqdsKvRubkanWZHoNGNN1ps\n93AndwAAAJYtXrx4zZo1ioqKwu1G/Rm3/v777/Pnz1MU1ehENx9P2MXlchudm6vRaXYEGt14o8X2\nT9gJj304gwUAAEJ04cIFDw8PYffi/3v06JGJiUlpaWmD+r///mtmZhYeHk7T9LRp02bPnt2vXz9F\nRcUjR44I2qSkpDBpQXAnBR8fH0GEuHHjRlM7rb/x5ovtFQLWF0PAAgCAphQXFw8ZMqT+7TeFq9EZ\ntz6e6KapCbsazM3VzDQ7zWy8qWL7hkHuAAAArGEmP5aRkRF2R/7P2rVrp0+frqenJ6hER0c7OTlZ\nW1vfuXNnwIABhBA+n5+Tk7Nr1y5lZWU7OztJSUmKoqZOnVpdXU1RlLy8/JgxY5j55vft23f48OFe\nvXopKSk5OTlVV1c32N3HG2+q2O5hDBYAAAA7Tpw4oaenZ2ZmJuyO/J/MzMwrV66kpaUJKnFxcf7+\n/jExMVJSUpWVlRRFycjIUBQlLS2dlJQ0YsQIPz8/V1dXiqKuXr0aERExc+bM4uJiT09PLy8vQkhJ\nSUl8fLyOjs6zZ8/c3d2XLFlSf3eNbrzRYmu/EEIh7FNo7MMlQgAAaH15eXnW1tb15w8Vuo9n3Gpq\nopuPJ+xqdG6uv/766+NpdprfeDNT67RvmCrni2GqHAAAaIDP5zs6Om7ZsqVv377C7gu0CRiDBQAA\n8K327t07fPhwpCsQwBgsAACAb5KRkRETE3Pp0iVhdwTaEAQsAACAr1dbW7tgwYJjx441uOUm/ODw\nfwMAAMDX27hx44wZM7S1tYXdEWhbELAAAAC+UlJSUmZmJjPDDEB9uEQIAADwNd6/f+/j4/Pnn38K\nuyPQFuEMFgAAwNdYvnz50qVLlZWVhd0RaIsQsAAAAL5YTEzM+/fvx4wZI+yOQBuFS4QAAABfprS0\ndMOGDTExMcLuCLRdOIMFAADwZby9vQMCAjp27CjsjkDbhYAFAADwBU6dOqWiomJhYSHsjkCbhkuE\nAAAAn+vcuXO+vr5ZWVnC7gi0dQhYAAAAn5aUlOTv79+hQwcHB4cOHToIuzvQ1iFgAQAANIfD4axe\nvVpTUzMoKOju3buvX78Wdo/gO4CABQAA0Li0tDQ/Pz9VVdWDBw8yk+Hk5+f36dNH2P2C7wACFgAA\nQEM5OTnLli0TERHZuHFjv379BPX8/Hx7e3shdgy+FwhY/6O2traysrL5NuXl5bW1tWVlZU01EBUV\nxW93AQC+U0VFRf7+/vn5+f7+/qampg2W5ufna2lpCaVj8H1BwPof48ePz3tTRlHUJ1taOjR5997C\n3Jy7f8f17NmT1a4BAEDLKi8v37x5c0JCwurVq21sbBptU1FRga/Q8DkQsP5HdXX16uBwUdFvellC\nNqyqqalhq0sAANDS3r17t3379ps3b/r6+m7atOlzvmYDNA83GgUAgB/Xhw8fAgICbG1t+/XrFxsb\nO3r06GbSVXl5uYKCQmt2D75fCFgAAPAj4vF4hw8fHjZsmIKCwq1btyZMmCAi8om/iXl5eRiABZ8J\nAas5dbW186xNyt+8ElTKS177uzrPHNxn83yXqncVzRQBAKBtomn6zJkzQ4cOLSsri42NdXd3l5CQ\n+JwVMcIdPh8CVpOijwevnjb21Yvn9Yuhm9dq6fUJvv2PpHSHswd3N1MEAIA26Pr167a2tg8fPrx4\n8eLvv/8uIyPz+esiYMHna4eD3MvKyg4fPvx167548ULwWKtXny6aWpvmzhRUaJq+d+PKqqAwcQkJ\nE9uREbu3zFi6+uOiwZDh33oMAADAtlu3bq1bt87IyCgiIkJZWfkrtpCfn29ubs56x6BdaocBy8XF\nhZXt6JsMaVCpqaqsqars2r0HIURDt2fZq5dNFQEAoO148ODBsmXLtLS0jhw5oqmp+dXbwRks+Hzt\nMGAtXLjwq9c9c+bMJ9swPzChaZrP4zVfBAAA4crOzl65ciVFUTt37vz2KW7evHnzdae+4AfUDgNW\ny5HqICPVQaYw92mvQUZFeTmKqmpNFQEAQLgKCwvXrl37/Plzf3//wYMHs7JNmqZZ2Q78CDDI/QtQ\nFGX8s13K3zdomr4fH2dsbd9UEQAAhOXVq1fe3t5ubm7Tp0+/dOkSW+mqsrKyQ4cOrGwKfgQIWF/G\nZdmaJ/c57paGFaUlE+Z5N1MEAGhnaJo2NzeXkpLi8/nNt0xISBg/fryKioqmpqaPj09VVRUhhM/n\nd+rU6dGjR0ybf//9t1evXgcPHmSxhxUVFcuWLRs3btzo0aMvX75sYWHB4safP3+OAVjw+Sic8KzP\n1ta2ywBjivqm3Jkce+VM6JH+/fuz1SsAgLYgMjJy27ZtKSkpubm56urqTTULCwtbv379jh07hg8f\n/uLFi5kzZ5qYmOzatSsrK2vQoEEVFRWioqKJiYmTJ0/ev3//L7/8UlFR4ePjc+HCBVFR0blz5/r5\n+X1F32pqanbv3v3nn3/6+PiMGzfuk7cM/QpXrlxJTU1dvnw561uGdgkB638kJyfn5eU13+bdu3eH\nDx/29fVtqoG4uPjIkSMlJSVZ7hwAgPDU1dUNHDjwxIkTkyZNOnr0aFMnh54/f25hYZGUlKSqqspU\n7t+//8svvxQVFYWHhx84cOD27dt//fWXr6/vmTNnBg0aRAhxdHRUVFTcsmXLu3fvevfuXVNT80Xx\niMvlHjlyJCQkZNasWa6uruLi4t9+sI06fPiwrKzslClTWmj70N7Q8IXevHkzbtw4YfcCWs/JkycN\nDQ3rV2bPnu3r69to45kzZ4qKioqKihJCREREmMfPnz+nafru3bvjxo1TVlbu1q3bokWLKisrW6P3\nAN+Mz+ebmZmJiYm5ubnRNG1tbR0aGsosiouL69+/v7Ky8rJly6SkpBYuXDhnzpxt27bRNH3u3Dlb\nW1uapt+/fy8uLs7j8RYuXCgqKmplZaWurv7s2TOmjY2NTVRUVGVlZVFR0Z49e8zNzT+/YzweLzQ0\ndOjQobt27WqFN9SKFStu377d0nuBdgNjsAA+oUePHrm5uYKn2dnZUVFRTV0mOHbsGJfLrayslJCQ\neP78OZfL5XK5GhoaYWFhrq6ubm5uOIS0dAAAIABJREFUubm5169fT0hIWLFiBbNKQkKCubl5586d\nu3fvfv78+dY4JIAvERUVxePxeDyel5cXIURbWzsnJ4cQQtP05MmT9+3b9/jx4+joaB6Pd/z48fj4\n+JEjRxJCOByOkZERIeT06dNcLpfP58fHx/P5/Nu3b5eVlfF4PKaNsbGxk5NTZmammpqal5dXWlra\nZ74RmLluioqKLl265O3t3QrDz3ETLPgywk543x+cwfrRlJWVEULevn3LPJ0yZcqGDRuaX4XD4aip\nqQmePnv2TFNTs6ioSFBJTU1VVVWlabqwsLBz584xMTE8Hu/IkSNMEaDtqK2t7dOnj4eHR/0/HDNm\nzKBpms/nGxkZ2djYHD58eO/evePGjevZs6eIiAhzMsnW1vbs2bM0Ta9evdrJyYnH40lKSvbs2bNn\nz55Dhgzx9vau3+bFixcKCgru7u4aGhqffCPExcVZWVn9/vvvr1+/bo2X4L+srKy4XG5r7hG+a7gP\nVjtRWFhYU1Mj7F60Q5qamvLy8oqKinl5eQMGDMjIyLh58+bhw4fz8vJsbW3//fdfQkhycvJvv/2W\nmJgoWEvw3Z2xfv16Ly8vwZAUQoienl5JSQnznd7FxcXe3p4QYmBgwFxbBGg7QkJCBg0aFBMTY2lp\n6erqOmPGjNDQ0O3btw8YMKC4uHj27NlWVlZhYWHnzp1buXKlqanpypUr//333wEDBqSkpBgZGRUV\nFe3YsWPZsmVZWVl8Pt/NzU1MTCw8PDw4ONjPzy8lJeXEiRO//PLL7du3XV1dPT09r1271tQboaSk\n5MiRI1FRUTo6OoGBgT179mzll4LH4+EdCp8PAaudMDc3Z+teLyBAUdSiRYtMTU11dXWZgOXn57dy\n5UoZGZkLFy4IpiS7e/euqalp/RU5HI6hoaHg6e3btxtMMJCbm8t813d2dnZ2diaEZGdne3p6bt26\nteUPC+BzvX//fteuXebm5q6urgUFBcyVQTU1tYcPH8bFxTk6OkZFRfn7+5eUlJw8edLIyMjAwGDZ\nsmWurq47d+4khKSnpy9cuFBJSWnUqFEpKSlSUlJMm1WrVvXq1Wvz5s0URd25cyciImLmzJlDhgyZ\nPXv2lClT6r8RiouLb926FR8ff//+fYqicnJyNm/ePH369NZ/Kerq6lpu+Dy0SwhY7YSOjs7p06eF\n3Yt2ixmGxeFw0tLSTp48SQhJTEysH7DGjh1bvz2HwxkzZgzzmM/nZ2VlaWtr129w6dIlweqlpaV+\nfn6pqamHDh0aMGBASx8LwOfbsWPHkydPnjx5wjydMWMGIaRr1658Pn/t2rUjR47MyMjo0qVLz549\naZo2NDRUUFAYN25cdnb22LFjKysrQ0JC9u3bN3HiRH19/WPHjlVXVzNtfv3117dv3wYFBRkbG7u7\nuy9ZsmTx4sVdu3ZVUVGJjY1l3gjLli27ffu2srLy8OHDXV1dzczMgoOD09LSFBUVhfJSFBQUdOvW\nTSi7hu+VsK9Rfn/a5hgsS0tLYXehPVu9evXChQtHjBgRHh7OVExMTNLS0miaLisrU1BQyM3NFTSu\nqqoSExMrLCwUVFRVVR88eCB4WlhYqK6unpOTQ9P0xYsXhw4dGhkZyefzW+dYAD7Ty5cvNTU1Kyoq\nmKfHjh2zsLBgHnO53KtXr86cOVNTU5PP59+/f7979+7Mops3byoqKi5cuHDFihU0TcfHxw8dOpSm\n6Wba0I29EQoLCwWPAwICRo8e/f79+9Y47CbExcWtWbNGiB2A7w5+RQjwaT169IiKinr58uXEiROZ\nSnp6upiYGJfLXbp0qYSERP3fFqWlpamoqKip/f9ZKSdNmjR//vynT5++f/8+Ojp66NCh/v7+Ojo6\ncXFx/v7+kZGRtra2lZWVlZWVrX1gAE1bt26dq6urnJwc85T5VkAIUVBQuH79urW1tb29vbi4OEVR\nzHArptmwYcOUlJQCAwOZSmJioomJCSGkmTaNvhHU1NQoiuLz+d7e3rm5uVFRUTIyMq3+Gvx/+Akh\nfCkELIBP69GjR35+/qZNmwT3P/T09Bw+fLilpWWnTp1MTU0pihI0bjAAixCyceNGY2PjYcOG6erq\nBgcHBwUFzZkzhxCyZcsWDoejrKwsJycnJyfHDHUHaAuys7PDw8OZ+zIwunbtWlhYWF1dvXbt2t9/\n/11RUXHnzp0RERHkf/+fpyhq7ty5/6+9Ow+oqs7/P/65l01ABAVUXLi4oaa5AaKoJZpopqaNkmWj\ngFtpX9SZsknDtWZSK2sqM5cQE0xLzG1EJ8VCBBPQzBRHVBRRXBBcQAQu5/fHmR9Dxu6Bc5fn4697\n3/dzPud9zPDFOZ97jnw1UJQJWJWMqeR/hGnTptna2q5atUr11eUELNQUd3Kvsezs7OnTp3/33Xdq\nN/I7/v7+sbGxancBAEqKj4/v16+f2l0IIcTkyZPnzZvXrl07tRuB0eAMFgDAQBlIuhJCZGRksMgd\nNULAAgCgCoWFhdbW1mp3AWNCwAIAoDIlJSVl11kC1UHAAgCgMtevXy/7vWCgOghYAABUhq8QohYI\nWAAAVCY9PZ2AhZoiYAEAUBnOYKEWCFgAAFSGgIVaIGABAFCZS5cuubu7q90FjAwBCwCAyuTl5TVs\n2FDtLmBkCFgAAFSGZ8qhFghYAABU6Pbt2y4uLmp3AeNDwAIAoEKscEftELAAAKgQAQu1Q8ACAKBC\nBCzUDgELAIAKEbBQOwQsAAAqRMBC7RCwAACo0O3bt5s0aaJ2FzA+BCwAAACFEbAAACjf/fv37e3t\n1e4CRomABQBA+RISEvLz89XuAkaJgAUAQPn27t1bUFCgdhcwSgQsAADKIUlSYmIilwhROwQsAADK\ncfjw4d69excXF6vdCIwSAQsAgHKEh4dPnjxZo9FIkqR2LzA+BCwAAB517969tLS0J598slmzZtev\nX1e7HRgfAhYAAI/atm3buHHjhBA6ne7SpUtqtwPjQ8ACAOBRUVFREyZMEAQs1BYBCwCA30lLS3Ny\ncpKfkEPAQu0QsAAA+J2IiIigoCD5NQELtUPAAgDgf/R6/f79+4cOHSq/JWChdghYAAD8z8GDB/39\n/S0sLOS3DRs2vHfvnrotwRgRsAAA+J8NGzaEhISUrWg0GrWagfEiYAEA8F85OTlZWVmenp5li40b\nN87JyVGrJRgpAhYAAP/1zTffvPTSS48UWYaFWiBgAQDwX1u3bg0MDHykSMBCLRCwAAAQQoiTJ0+2\nbt26UaNGj9QJWKgFAhYAAEIIsXHjxuDg4D/WCVioBQIWAACiqKjoyJEjAwcO/ONHBCzUAgELAACx\nd+/eoUOHlntHBmdn59u3b9d/SzBqBCwAAH73eJw/kiSpHnuBKSBgAVXYsmWLt7d32crUqVPfeOON\ncgcXFhba2NiMGDGibDEoKEij0Vy/fr2iXeh0umPHjpW+/fzzz59//vnH6xpADWRlZd2/f1+n01U0\nwM7OLi8vrz5bgrEjYAFVaN++/cWLF0vfpqWlbd++/e233y538KlTp7RabdnlGqdPn/7uu+9at27d\nrFmzcjfJzs7Oysrq1q1baeXnn3/28fEpfVtYWOjn55eamvq4RwKgAps3b/7zn/9cyQB3d/fLly/X\nWz8wAQQsoArt2rW7ffv23bt35bcLFy6cM2eOs7NzuYOTkpICAgIuX75cekEhLCzsmWeeeeQcWFnH\njx/v1q2bjY1NaeXYsWNywCouLn777bd9fX1TU1MfubU0AAVFR0f/6U9/qmQA69xRUwQsoApOTk7O\nzs7p6elCiFOnTh06dGj27Nnp6ekdOnSQBxw7dqxPnz7y6+Tk5P79+9va2ubm5gohkpKSMjIyWrVq\nJQesX3/9ddCgQS4uLh4eHvPnz5c3OX78eNnzVffu3Tt79qw8vqioqF+/fs8884yvr69Wy/+tQJ1I\nTk5+4oknbG1tKxnj4eFBwEKNWKrdAFBLDx48GD16tIeHR53uJSQkxNfXt127dunp6d26dVuwYMH8\n+fPt7e137drl5+cnjzly5EhpwEpKSho/frz8y27jxo3nzZv33nvvLViwYNSoUUKI0NDQCRMm7N+/\n//Tp0z4+PkuXLtVqtcePHx8yZEjpHpOTkz08POQzZLa2tiNGjPjpp5/69u1bp4cJmLPw8PByb39V\nlk6n27VrV/30A9NAwIKxevDggVarfeutt+p0L/LCKXkZVlJS0smTJ7/55hshRGJiYtmANWbMGCHE\nw4cPf/vtt169eskBKzc3t7i4eODAgSdOnPDy8hJC7Nix48SJE1FRUQcPHuzZs6d8UurEiRPz5s0r\n3WPp9cFSCQkJYWFhdXqYgNkqKCg4ceLEZ599VvkwLhGipghYMGJ2dnZt27athx3JZ7D27t27dOlS\na2trIURiYuLkyZOFELm5uf/+97+XLVsmhPj11191Op2jo6P8s/ibb7758MMPT58+7ebm5uzsvH79\n+tWrV48YMaJ3795yRQhRUlKSlpbWrl270n3FxMRMnDix9G1hYWFycnLv3r3r4TABM7Rjxw75BHPl\nmjdvnpWVVQ/9wGSwqgOoWvv27bdv3379+vUXX3xRrvz666+WlpbFxcVz5861traWv92dnJwsr53S\n6XRr1qxxdXXt27dvUlKSXJw7d+66devCwsKcnZ23bNkiX1XUarWtWrX68MMPb968mZOTExYWdvv2\n7bLPmj1x4kSbNm2cnJxUOGzADGzatGnSpElVDtNqtSUlJfXQD0wGAQuoWvv27S9duvSPf/yjdKX5\nzJkzn3766YEDBzo6Ovbp00e++3NSUpJ8KdDDw+P06dNLly6Vi3LACgoKGjFixBNPPLFjx45OnTpt\n2LBBnmrr1q179uzx8PDo0aOHfD6s7GLbI0eOsAALqCMZGRkajaaiW6g8wsrKqqioqK5bgsnQcHfa\nmsrOzp4+ffp3332ndiO/4+/vHxsbq3YX9er27dtTp07dtm2b2o0AMFZ///vfO3bsWPkNGkqFhISE\nhYW1adOmrruCaeAMFgDAHEmStHv37pEjR1ZzPOvcUSMELACAOYqPj/fx8ZG/tlIdBCzUCAELAGCO\nqnP7q7IIWKgRAhYAwOzcv3//P//5T48ePaq/CQELNULAAgCYnW3bto0bN65Gm7Ru3TojI6OO+oHp\n4UajAACzExUVFRUVVaNNuE0DaoQzWAAA83L+/HlHR0f5aQo1otFouN0oqomABQAwLxEREUFBQbXY\n0M3NjQfmoJoIWAAAM6LX6/ft2zd06NBabMs6d1QfAQsAYEYOHjzo7+9vYWFRi20JWKg+AhYAwIxs\n2LAhJCSkdtsSsFB9BCwAgLnIycm5du2ap6dn7TYnYKH6CFgAAHOxZcuWl156qdab63S6y5cvK9gP\nTBj3wQIAmIutW7d+//33td7czs4uLy9PwX5gwjiDBQAwC7/++mvLli0bNWr0OJNIkqRUPzBtBCwA\ngFnYuHFjjZ7uXC4XF5fs7GxF+oFpI2ABAExfUVFRfHy8v7//Y87DOndUEwELAGD6YmJiAgICNBrN\nY85DwEI1EbAAAKYvIiLi8a8PCgIWqo2ABQAwcTdu3MjLy9PpdI8/FQEL1UTAAgCYuMjIyAkTJigy\nFQEL1cR9sAAAJi46Onr//v2KTNW4cePc3FxFpoJp4wwWAMCUJScnd+7c2dbWVqkJuRUWqoOABQAw\nZeHh4Yosby9lb29/7949BSeESSJgAQBMVkFBwYkTJ/r27avgnDqdLiMjQ8EJYZIIWAAAk7Vz586R\nI0cqOyfr3FEdLHKHIfrss89++umnyscUFhYmJSUFBgZWNMDKyuqf//yns7Oz0t0BMBqbNm1as2aN\nsnPqdLr09HRl54TpIWDBEO3cuXPUXxdaWFhUPmxopZ9u/XxlZmYmAQswW/KFvObNmys7rU6nO378\nuLJzwvQQsGCgmrVqbWFp9Tgz2Ds4KNUMAGO0adOmiRMnKj6th4cHlwhRJQKWoTt27Fh1flW6evVq\n5afBx4wZ4+rqqlxfAGDQJEnavXt3bGys4jM3a9bs+vXrik8LE0PAMnRr1qyRmns4uVSRjUbO+Ot/\nHlT4afKhH9zc3BRf6VlvigoLQ58d8I8tu5xcmsqV3OybH78x8+Jvv3b29g1d9k87h0Z/LKraMgCV\nHTlyxNvb29raWvGZNRoNt8JClQhYRqDXU4Oau3s8zgy3rmUq1IsK9mxcF7d7+43M330pOuL9xTrP\nzvO/3PTZ27O/++KTiXPD/lhUqV8ABiE8PPz111+vo8mtra0fPnxoY2NTR/PDBHCbBhg6XcfOY2fM\nLluRJOnnA/v6BDxnZW3tO+TZY7H7KyoCME/3798/e/Zsjx496mj+1q1bX7lypY4mh2ngDJaRUepi\n2bhx41JTUw32169z585N/f+vu/r2e+TTgvy8gvy8lm3bCyFateuQc+N6RUUA5mnbtm3jxo2ru/nl\nW2G1a9eu7nYBY0fAMia1vljm9Ifl7S1atAgLC+vWrVv9dV8TAQEBVY7RaDRCCEmSSvT6yosAzE1U\nVFRUVFTdzc+9RlElLhEaEy6WyRrY2Tews7968bwQ4lr6BefmbhUVAZih8+fPN2rUqE7vgUfAQpU4\ng1VjkiQlJCQMGTKkfnZ35syZnoEh8msulsk0Go3PoIDkHw949vA6HhfrM3hYucXiwodqdwpABRER\nEUFBQXW6CwIWqkTAqjEXF5fMzPr7Ut7UqVOrHGOGF8uC/rZw5V9nThvo1b5r96C3FpZb3PzJcnWb\nBFD/9Hr9vn37FixYUKd7admyJYvcUTkClnErvS7Wsaf3Hy+WlS0anft37lhY/e9O7hE/nxFC3L97\nR35raW3z5qfr5Nf6khK5/kix8CFnsACzExsbO3DgQEvLuv3XzcrKqri4uE53AWNHwDJu1bxYpnab\nNdavX7+ohX+pfExRUdGZM2cqWadvZWXFzesBc7Nhw4a6Pn0l02q1er2+ykemwmwRsIxedS6W/fvb\nSHWbrKmFCxdWOeb27dtTp07dtm1bPfQDwCjk5uZevXrV09OzHvbVokWLa9eutWrVqh72BWNEwDI+\n21Kvln3r5NJ0ccS3j4wptwgApm3Lli3jx4+vn33J69wJWKgIAcvQOTk5fTZvjtXjPU7rzq1bEwOe\nVqolc3Pnzp0FCxZ8//33d+7cmTVr1uLFi+V6QkLCX//619TUVCcnp48//njUqFFyvbCwcODAgV99\n9VWnTp1qOrKsRz6VJOnJJ5/87bff5E+DgoLCw8Pr9MABo7N169bt27fXz77kgNWv36Nf7gZkBCxD\nt2LFihXVGObv718XD42HEGLChAl9+/ZNTk6+detWt27dFi5cqNVqr127NmLEiMjIyICAgIiIiOnT\np48aNaq4uDgsLCwmJubSpUulFymqP7JUuZ/Gx8c7Ozvfu3dPfmtV5hsAAIQQp0+fbtGiRaNGjepn\ndzqdLjk5uX72BWNEwAKqEB8fP2XKFAsLi2+//dbb21ur1Qoh4uLigoKChg0bJoTo0aOHvNC1qKio\nX79+xcXFp06dkofVaGSpcj/98ssv8/LyWrZs2aRJk1WrVj377LP19QcAGIfw8PDg4OB6251Op4uO\njq633cHoELBgRnJycnJycmq0ibu7e0hIyJgxY+S3Bw4ckF8EBgYGBgYKIdLS0mbOnCmfZ7S1tR0x\nYsRPP/3Ut2/f0hmqP7LUHz/Nz89PS0ubP3/+008//cknn0yZMqU+b8YGGL6ioqL4+Pjly+vv7nfu\n7u6XL1+ut93B6BCwYEaGDBni7u5e/RvkaDQad3f3EydOpKamOjs7b9q0aebMmWfOnJE/vX379oIF\nC1JSUlavXl32bhEJCQlhYWFl56n+yLLKfmpnZ5eQkCC/nj179ooVK0pKSv546gswWzExMQEBAfIN\nluuHra3tgwcP6m13MDoELJiRVq1ahYeHOzo6Vn+TJ554YteuXe3atRNCjB49+uOPP5bre/bsWbZs\n2Zw5cz799NOyP9MLCwuTk5N79+5dWqn+yLIe+XT79u3du3dv27atECItLa1Xr16kK6CsiIiIDz74\nQO0ugP8hYAGVyc7OjouLa9OmzeXLl6dNm/bmm28KIWJjYxctWrR3794GDRrk5eVpNBp7e3t5/IkT\nJ9q0aePk5CS/rf7IRzzy6TfffLNt27Yvv/wyOzs7NDT073//e90eNmBUbt68ef/+fQ8Pj3rer4uL\ny61bt1xcXOp5vzAK/BIMVGbt2rXvvfeeq6vr2LFjQ0JCZsyYIYRYvnx5UlKSq6urg4ODg4ODvIBd\nduTIkbLLqqo/8hGPfLp06dK0tLSWLVu+8MILixYt8vf3V/g4AWMWGRk5YcKE+t+vTqdLT0+v//3C\nKHAGC6jMqFGjSm9bVWrv3r0VjZ89e3btRlb+qaenZ2JiYmWNAmZs27Zt+/fvr//9yrfC8vb2rv9d\nw/BxBgsAYMRSUlI6depka2tb/7v28PC4dOlS/e8XRoGABQAwYvV8+6uy5DNYquwaho+ABQAwVg8f\nPjx+/Lifn58qeydgoRIELACAsdqxY8fIkSPV2rujo+OdO3fU2jsMHAELAGCsNm3aNGnSJLW7AMpB\nwAIAGKWMjAwhRPPmzVXswcHB4e7duyo2AINFwAIAGKVNmzZNnDhR3R50Oh1PJES5CFgAAOMjSdKu\nXbv+eJu6esY6d1SEgAUAMD5Hjhzx9va2trZWtw0CFipCwAIAGB8Vb39VFgELFSFgAQCMzP3798+e\nPduzZ0+1GyFgoUIELACAkYmOjh47dqzaXQghRNOmTW/cuKF2FzBEPOwZAGBkoqKiIiMj1e5CCCE0\nGo0kSWp3AUPEGSwAgDE5f/58w4YNnZ2d1W7kv2xsbAoKCtTuAgaHM1gwBd9///3/vTG3SbOq7zf4\n1PAKn6rxIC9v3HPD3nvvPUVbA6CwjRs3BgUFqd3F/7i7u2dkZHTo0EHtRmBYCFgwBfn5+SMmTR0S\n+MrjTHLx9KnMuL1KtQSgLuj1+piYmLCwMLUb+R95nTsBC4/gEiEAwGjExsYOHDjQ0tKAzg7wRUKU\ni4AFADAaGzZsmDx5stpd/A4BC+UiYMEEFRUWvjbYN/fW/747nZt9c1Fw4KTend+fEZR/724lRQAG\nKzc39+rVq56enmo38jsELJSLgAVTs2fjurBXxtzIzChbjHh/sc6z87rDv9jY2n33xSeVFAEYrC1b\ntowfP17tLh7VsmXLq1evqt0FDA4BC8ZKkqSrV69+++2333777dGjR0vruo6dx86Y/cjInw/s6xPw\nnJW1te+QZ4/F7q+oCMCQGWbAsrCwKC4uVrsLGBwDWicI1IiDg8OAAQNycnKEEHl5eVYu/6139e33\nyMiC/LyC/LyWbdsLIVq165Bz43pFRQAG6/Tp0y1atGjUqJHajZTDwsJCr9dbWFio3QgMCAELxsra\n2nr58uXy64YNG8adz6x8vEajEUJIklSi11deBGCANmzYYAhPdy5Xy5YtMzMz3d3d1W4EBoRLhDB9\nDezsG9jZX714XghxLf2Cc3O3iooADFNRUdHhw4cHDRqkdiPlY507/oiABdOn0Wh8BgUk/3hAkqTj\ncbE+g4dVVARgmPbt2zdkyBD5lLMBImDhjwhYMAtBf1t49njStIFed29nj3ttViVFAAYoIiLCYK8P\nCgIWysMaLJiI00lHy759dcnyYwd/98XA/s89L784/K8d5RZvXbvauoGB/n4MmLObN2/evXvXw8ND\n7UYqpNPptm7dqnYXMCwELJgCf39/KyurKoetWLFi5syZdnZ25X/cyaNLly4KdwbgsUVGRr7yymM9\nabSuyc97VrsLGBYCFkyBm5vbuHHjqhwWGRk5evRoR0fHemgJgFK2bdu2f79B36nOxsamoKBA7S5g\nWFiDBQAwXCkpKZ06dbK1tVW7kapJkqR2CzAgBCwAgOEKDw835OXtpZo2bXrz5k21u4ABIWABAAzU\nw4cPU1JS/Pz81G6kanyREI8gYAEADNTOnTtHjhypdhfVQsDCIwhYAAAD9fXXX0+cOFHtLqpFp9Ol\np6er3QUMCAELAGCIMjIyJElq0aKF2o1UC2ew8AgCFgDAEEVGRhrL6SshhIeHBwELZXEfLACAwZEk\naefOnYcOHVK7kepycHC4d++e2l3AgHAGCwBgcBISEry8vKytrdVuBKglAhYAwOAYy+2vynJ0dLxz\n547aXcBQELAAAIbl/v37qampvXr1UruRmmGdO8oiYAEADEt0dPTYsWPV7qLGCFgoi0XuAADDEhUV\nFRkZqXYXNUbAQlmcwQIAGJDz5883bNjQ2dlZ7UZqjICFsghYQBU+/fRTze+lpqYKIRISEvz8/Jo0\nadK2bdudO3fKg3Nzc1988UU3N7c+ffp89dVXQoi7d+9OmTKlWbNmLVq0WLJkiTzszp07s2bN0ul0\nTk5OCxcu/ONOy5283CJgYjZu3BgUFKR2F7VBwMLvSDAJAwcOrNH40NDQX375pY6aMVjPP/98bm5u\nTbcqKCi49/8tW7bs1VdfLSkpuXr1apMmTfbu3avX67/66qvmzZtLkqTX6wcMGLB+/Xq9Xh8XFzdm\nzBhJkkaOHBkUFHTjxo3z589bWVnp9XpJkp577rl333335s2bZ86cKS2WKnfycouAidHr9X369Ckq\nKlK7kVry9/dXuwUYCtZgAVWwsbGxsbERQuzfvz82Nnbnzp0ajSYuLi4oKGjYsGFCiB49elhYWAgh\nDh48aGNjExISIoTo379///79i4uLQ0JCAgIC7t69u2fPHh8fH61WK4SIj4+fMmWKhYXFt99+6+3t\nLRdLlTt5uUXAxMTGxg4YMMDS0lj/bSopKVG7BRgKY/1LDNTClStXZs+eXaNbF4aEhPj6+gohzpw5\ns2DBgr1791pZWQkhAgMDAwMDhRBpaWkzZ85csWKFECIiIsLDw6Nr165ZWVkrVqwIDg62tLQcPXp0\nSkqKl5eXECImJqZ02jFjxshqCB+5AAAgAElEQVSvDxw48MhOy5283CJgYjZs2PDOO++o3UXt2dra\n5ufn29nZqd0I1EfAghnZvHlzTU/8NGvWTAiRk5MTGBi4bdu2xo0bl350+/btBQsWpKSkrF69ulu3\nbkKIlJSU7t27x8bG/vvf/37ttdcmTZokn5rq2bNnTk5OeHj4lClTMjIyVq1adeLEidTUVGdn502b\nNs2cOfPMmTOP7PePk1dUBExGbm7ulStXOnbsqHYjtafT6TIyMoz6EKAYta9RQhmswapTs2bNWrZs\nWdnK7t27BwwYEB0dXVJSIlf0en2DBg2uX78uSdLNmzddXV1LSkpefvnl/Px8ecD58+fbtGkjSVLn\nzp3T0tLk4sWLF3U63SO7++PkFRUBU/Lll1+uXr1a7S4ey9///vd9+/ap3QUMAmewgCqkpqbu27fv\n5MmTpZXY2NhFixbt3bu3QYMGeXl5Go3G3t5eo9HY2toePXp06NChCxYsCA4O1mg0+/fv37x586RJ\nk7KysmbOnBkaGiqEyM7OjouLa9OmzeXLl6dNm/bmm2+W3V25k5dbrO8/CKCObdmyJTo6Wu0uHgtf\nJEQpbtMAVOGNN96YNWuWvPRKtnz58qSkJFdXVwcHBwcHB3nhuUaj+fzzz6dNm9auXbuGDRvKd2T4\n4osvli5d6urq+txzz40bN27WrFlCiLVr17733nuurq5jx44NCQmZMWNG2d2VO3m5RcCUnDlzxs3N\nzdHRUe1GHgsBC6U0kiSp3QMU4O/vHxsbW/3xs2bNmjx5Mut4ABiIt956KyAgYPDgwWo38liuXLny\n9ttvf/3112o3AvVxBgsAoLKioqK4uLhBgwap3cjjcnNzu3btmtpdwCAQsAAAKtu3b9+QIUM0Go3a\njTwuCwsLvV6vdhcwCAQsAIDKIiIigoOD1e5CGZaWlsXFxWp3AfURsAAAarp58+bdu3c9PDzUbkQZ\nLVu2zMzMVLsLqI+ABQBQU1RU1IQJE9TuQjE6nS49PV3tLqA+7oMFAFDT1q1b9+/fr3YXiuFODZBx\nBgsAoJrjx4937tzZlG6c6+HhQcCCIGABAFQUHh5uMsvbZZzBgoyABQBQx8OHD5OTk/v166d2I0pq\n3bp1RkaG2l1AfQQsAIA6du3aNWLECLW7UJi1tXVhYaHaXUB9BCwAgDq+/vrriRMnqt2F8jQaHkMH\nAhYAQA0ZGRklJSUtW7ZUuxHlNWvW7Pr162p3AZURsAAAKoiMjDTJ01eCde4QQhCwAAD1T5KknTt3\nPv/882o3UicIWBAELABA/UtISPDy8rK2tla7kTpBwIIgYAEA6p/p3f6qLAIWBAELAFDP8vPzz549\n26tXL7UbqSsELAgCFgCgnm3btu2FF15Qu4s61LBhw3v37qndBVTGw54BAPUqKipq06ZNandRtzQa\njdotQGUELBOUk5NT5ZiHDx/evXu3kpF2dnY2NjaK9gUA4vz583Z2ds7Ozmo3UrcaN26ck5PTuHFj\ntRuBaghYpubYsWMvvDTBvUOnKke+9d7yij4qLips69p48+bNirYGAOLrr78OCgpSu4s6Jy/DImCZ\nMwKWqSksLPQdMvzl2W89ziT37+TuWPaOUi0BgKykpGTfvn3vvGP6P17kgNWjRw+1G4FqWOQOAKgn\nhw4dGjBggKWl6f9uzxcJQcACANST8PDwyZMnq91FfSBggYBl4ooKC18b7Jt760ZpJTf75qLgwEm9\nO78/Iyj/3t1KigCgoNzc3CtXrnTs2FHtRuoDAQsELFO2Z+O6sFfG3MjMKFuMeH+xzrPzusO/2Nja\nfffFJ5UUAUBBW7duHT9+vNpd1BNnZ+fbt2+r3QXURMAyETdu3FizZs2aNWt27txZWtR17Dx2xuyy\nwyRJ+vnAvj4Bz1lZW/sOefZY7P6KigCgrC1btphPwBJCSJKkdgtQk+mvNDQTU6ZMkb8P7ODgIO4+\nlItdffs9MqwgP68gP69l2/ZCiFbtOuTcuF5REQAUdObMmebNmzs6OqrdSP2xs7PLy8uzt7dXuxGo\ng4BlIubMmSO/iI+PP7Hp28oHy7cYliSpRK+vvAgAitiwYYMJP925XO7u7pcvX+7cubPajUAdXCI0\nLw3s7BvY2V+9eF4IcS39gnNzt4qKAKCU4uLiQ4cODRo0SO1G6hXr3M0cAcu8aDQan0EByT8ekCTp\neFysz+BhFRUBQCn79u0bOnSoVmte/+J4eHgQsMyZef11hxAi6G8Lzx5PmjbQ6+7t7HGvzaqkCACK\niIiIMLfrg4IzWGZPw9ccTEx8fPyrf32rR/+BjzPJw4IH+mvp27ZtU6gpAObr5s2bEyZM2L/f7L6e\nfPXq1TfffDMyMlLtRqAOFrmbmp49e658d1GVw1atWjVs2LC2bdtWNKBly5YKdgXASM2bN8/e3n7+\n/Pm1nmHz5s0TJkxQsCVj0bx586ysLLW7gGoIWKbGzs7umWeeqXLYrl27+vTp061bt3poCYDxSkpK\nmjFjRi02/Pnnn6dPn37lyhUXF5eff/5ZLubm5k6fPv2nn37S6XTTpk0LCQlRtFnDotVqS0pK1O4C\nqmENFgCgfJIkJSUldenSpaYb3rt3b+zYseHh4ZGRkVlZWQ4ODkKIkpKSUaNGDR06NDMz84MPPti9\ne3cdtGxYrKysioqK1O4C6iBgAQDKd/Hixby8vKioqBYtWuh0upiYGLn+66+/Dho0yMXFxcPDo/Tq\n4euvv75w4cJevXqtXLly69atI0eO7NGjR//+/e/duyefyDl48KCNjU1ISIhWq+3fv390dLRqB1Zf\nWrVqdeXKFbW7gDoIWACA8iUlJen1+rZt2166dOntt9+eMmWKXA8NDX355ZezsrJ27tz5wQcfyPkp\nMTHx0KFDa9eunTFjxtatW4cNGyaEyMrKcnd3l2/QEBER4eHh0bVrVxcXl/DwcBWPq97wRUJzRsAC\nAJQvKSkpKCjoz3/+s5WV1fjx47Ozs+UstWPHDk9Pz6ioqI8++qhnz55arTY/P//UqVNRUVFeXl42\nNjbnzp178sknhRAnT57s0KGDPFtKSkpeXl5sbOw///nP2bNnm8P6JAKWOSNgAQDKl5SU1L9/f/n1\n0aNHvb29tVrt+vXrBw8eHBsb6+rq6ubm1q9fPyFESkpKjx495G8fl5SUZGRkNGvWTAgRExMTGBgo\nFy9cuPDxxx+7uroGBATY2NjIj+cybTqdLj09Xe0uoA4CFlCF4uJiW1vb0qc9CiGio6MDAgIq3yoh\nIWHs2LFNmzZ1d3f/y1/+kp+fX8dtAgorKSk5fvz4gQMHCgoKTp48GRoaumjRIiHE3Llz161bFxYW\n5uzsvGXLlj59+gghEhMT/fz85A01Gk2jRo3i4+NTUlJ2794tByyNRmNra3v06NHCwsIFCxYEBweb\nScDiDJbZImABVTh9+rRer9+4cWNBQYFcSUpK8vb2rmSTyMjI4ODgkJCQixcv/vDDDwkJCfPmzat8\nL4WFhX5+fqmpqYr1DTyetLQ0Ly8vBweH1q1bT548eeXKlYMHDxZCBAUFjRgx4oknntixY0enTp02\nbNgghEhMTOzbt6+8oUaj+fDDD8ePHx8SEhIdHS1/hVCj0Xz++efTpk1r165dw4YNlyxZot6R1Z/W\nrVuzyN1scSd3MzVr1qzJkydzH6zq+Oqrr/71r3+dPHly0aJFL7/8shAiICBg+vTpf/rTn8odn5GR\n0b9//6NHjzZv3lyuHD9+fPjw4deuXSt3fHFxcVhYWExMzKVLl27dumVuz2sDTJu/v39sbKzaXUAF\n/CgHqpCcnOzj4zNt2rQ1a9YIISRJSk5Ols9gxcXFPf30066urm3atNm4caM8funSpaGhoaXpSgjh\n6ekprw5OT08fN25c69atO3XqtGzZMvn+jUVFRf369XvmmWd8fX1JV4CJ0Wg05rCcH3/ET3OgCklJ\nSV5eXkFBQYmJiWfPnk1PT9doNO7u7gkJCW+++eaqVauuX7++YcOG6dOnyz9GDx8+/Oyzz5ad4eLF\nix06dLh06VKPHj08PDxSU1PDw8P/9re/+fj4CCFsbW1HjBih0WhKr7AAMBlubm48MMc88agcmJG8\nvLzCwsIabWJra3vy5EkvL6/GjRu/8MILa9eu7dOnj7e3t0ajWbhw4dGjR7t27SqPdHZ2lp+M8Z//\n/MfDw6PsJP/617/8/Pzmzp3bunXrZcuWabXavn37Ojk5yQFLlpCQEBYW9tiHCMCwyOvcW7RooXYj\nqG8ELJiRp556SqfTWVrW4K/90KFDW7Ro0bhxYyGEvO5Kr9d7eXkJIY4ePXr9+vWmTZsKIR48eCB/\nJUqr1bq6up47d6579+7yDNeuXfvkk08OHz7s6+v74osvyhcBU1NTCwsLO3XqJI8pLCxMTk7u3bu3\noocLQH1ywOL8tBkiYMGMtG7dOjw83NHRsfqbrF+/vvQLg0899ZSLi8uXX34ZGRkphPDw8Fi7du0b\nb7xx8uTJMWPGbN26Vf6a+vjx42fMmLFx48ZmzZr9+OOPs2bNWrRoUZs2bezs7I4dO5afn//gwYNp\n06b16tWrNOqdOHGiTZs2Tk5OSh8xAJXpdLpffvlF7S6gAtZgAZWRF2DJrzUazauvvvrgwQO5smHD\nhujoaFdX16lTp3722WelNwF67733fHx8nnrqqXbt2q1bt27t2rVTp04VQmzatKmgoKBNmzYvvfSS\nk5NT2d9ojxw5wi+4gEniVlhmi9s0mCnzvE3D6NGjIyIianQGSyl6vX7Xrl3PP/+8RqORJKlr166f\nffaZv79//XcCoD7l5+cHBgbu3r1b7UZQ3ziDBdQHjUYzYcKEzZs337lz591337WxsRkwYIDaTQGo\nc3Z2dnl5eWp3ARUQsID6oNVq//nPf86ZM8fDw+P48ePbt2+v0Vp7AMaLK0XmiYAF1JPJkydfv349\nJycnOjpap9Op3Q6A2qvRI0pdXFyys7MFjyg1MwQsAABqpkaPKHV3d8/IyOARpeaGgAUAQM0kJSWN\nGjXK2dk5Ojq6tFL6jeNHvP322w4ODvPmzTt06NDw4cPt7e09PT1XrVq1ZcuWiuYvLi5+++23fX19\nU1NTPT096+QYUMcIWAAA1EyNHlHq6uq6bNkyHlFqbvjPBgBAzfCIUlSJgAUAQA0UFRXJjyh1cXGR\nH1Eqn74q+4hSCwuLgQMH2tvbV/MRpfb29uU+orRPnz71fXhQCF8UBwCgBn777TceUYoqcQYLAIAa\nKF1uJco8olSuyI8offjw4bFjxzp06JCSkiIPkx9Rev78+fv37+/Zs2fAgAGPPKI0OzubR5SaGAIW\nAAA1wCNKUR08i9BM8SxCAFAXjyg1bZzBAgBABTyi1LQRsAAAUAGPKDVt/LcEAEAdkydPnjx5stpd\noE4QsGAKvv/++/97Y26TZs2rHPnU8JEVffQgLy9wxLPvvvuuoq0BAMwRAQumID8/f8SkqUMCX3mc\nSS6ePpUZt1eplgAA5ow1WAAAAAojYAEAACiMgAUTVFRY+Npg39xbN0orudk3FwUHTurd+f0ZQfn3\n7lZSBADg8RGwYGr2bFwX9sqYG5kZZYsR7y/WeXZed/gXG1u77774pJIiAACPj0XuMBTHjh377bff\n7O3ta7Ht0aNHhUsr+bWuY+dm7rp/vDqp9FNJkn4+sO+dtZFW1ta+Q57d/MnyiXPD/lgcMGKMMkcC\nADB7BCwYiujo6Hv37tXu6T15eXlWLv993dW33yOfFuTnFeTntWzbXgjRql2HnBvXKyoCAKAIAhYM\nhaura9++fUeNGlWLbRs2bBh3PrPyMRqNRgghSVKJXl95EQCAx8QaLJi+Bnb2Dezsr148L4S4ln7B\nublbRUUAABRBwILp02g0PoMCkn88IEnS8bhYn8HDKioCAKAIAhbMQtDfFp49njRtoNfd29njXptV\nSREAgMfHGiyYiNNJR8u+fXXJ8mMH95et9H/uefnF4X/tKLd469rV1g00ddwmAMAsELBgCvz9/a2s\nrKoctmLFipkzZ9rZ2ZX/cSePLl26KNwZAMAsEbBgCtzc3MaNG1flsMjIyNGjRzs6OtZDSwAAc8Ya\nLAAAAIURsAAAABRGwAIAAFAYAQsAAEBhBCwAAACFEbAAAAAURsACAABQGAELAABAYQQsAAAAhRGw\nAAAAFEbAAgAAUBgBCwAAQGEELAAAAIURsIAqzJs377333lO7CwCAMSFgAVVISkrq0qVLLTb8+eef\ne/bs6erqunHjxtJibm7uiy++6Obm1qdPn6+++kq5NgEABoSABVRGkqTaBax79+6NHTs2PDw8MjLy\njTfekIslJSWjRo0aOnRoZmbmBx98sHv3bqX7BQAYBAIWUJmLFy/m5eVFRUW1aNFCp9PFxMTI9V9/\n/XXQoEEuLi4eHh7z58+Xi6+//vrChQt79eq1cuXKrVu3jhw5skePHv379799+3ZJSYkQ4uDBgzY2\nNiEhIVqttn///tHR0aodGACgLlmq3QBQfzIzM2fPnm1tbV39Tdq1a6fX69u2bXvp0qX169dPmTLl\nypUrQojQ0NAJEybs37//9OnTPj4+S5cu1Wq1iYmJ9vb2a9eu7dq166hRo15//XUhRFZWlru7u1ar\nFUJERER4eHh07do1KytrxYoVwcHBdXSkAAB1EbBgRqKioiwsLGq0yaeffhoUFPTnP/9ZCDF+/Pg5\nc+aUlJRotdodO3acOHEiKirq4MGDPXv21Gq1+fn5p06dOn/+fMuWLYUQ586de/LJJ4UQJ0+e7NCh\ngzxbSkpK9+7dY2Nj//3vf7/22muTJk2SgxcAwMQQsGBGSoNO9f3yyy8TJ06UXx89etTb21ur1a5f\nv3716tUjRozo3bu3m5ubs7OzECIlJaVHjx5yuiopKcnIyGjWrJkQIiYmJjAwUC5euHAhNjbW1dU1\nICDAxsZGo9EoeXgAAIPBb89AhUpKSo4fP37gwIGCgoKTJ0+GhoYuWrRICDF37tx169aFhYU5Oztv\n2bKlT58+QojExEQ/Pz95Q41G06hRo/j4+JSUlN27d8sBS6PR2NraHj16tLCwcMGCBcHBwQQsADBV\nnMECKpSWlubl5eXg4NC6dWsPD4+VK1cOHjxYCBEUFDRixAh7e/s//elPnTp12rBhw7hx4xITE198\n8UV5Q41G8+GHH44fP75Vq1bR0dEODg5y8fPPP582bZqlpeVLL720ZMkSNY8NAFCXNJIkqd0DVDBr\n1qzJkyd369ZN7Ub+56OPPmrfvv2oUaPUbgQAgMfFJUIAAACFEbAAAAAURsACAABQGAELAABAYQQs\nAAAAhRGwAAAAFEbAAgAAUBgBCwAAQGEELAAAAIURsAAAABRGwAIAAFAYAQsAAEBhBCwAAACFWard\nAEzfw4cPP/vss+Li4sqHxcXF/fLLL2fOnKloQJs2bQIDA5XuDgAA5RGwUOdu3boV8d33L0z7v8qH\nef+pjRDiQcUDPv30UwIWAMAoELBQHxq7Nu3mN+AxJ/lx4xeKNAMAQF1jDRYAAIDCCFiob0WFha8N\n9s29daO0kpt9c1Fw4KTend+fEZR/724lRQAAjAIBC/Vqz8Z1Ya+MuZGZUbYY8f5inWfndYd/sbG1\n++6LTyopAgBgFAhYqFe6jp3HzphdtiJJ0s8H9vUJeM7K2tp3yLPHYvdXVAQAwFiwyN1MZWZmzpgx\nw9bWth729fDhQ42Tq/y6q2+/Rz4tyM8ryM9r2ba9EKJVuw45N65XVAQAwFgQsMzUd999V2/7yszM\nfPm10MrHaDQaIYQkSSV6feVFAAAMH5cIobIGdvYN7OyvXjwvhLiWfsG5uVtFRQAAjAUBCyrTaDQ+\ngwKSfzwgSdLxuFifwcMqKgIAYCwIWFBf0N8Wnj2eNG2g193b2eNem1VJEQAAo8AaLKhgW+rVsm+d\nXJoujvj2kTHlFgEAMAoELNS5Bg0aXD6Xujjkxcecp2PL5or0AwBAXdNIkqR2D4AQQnz00Uft27cf\nNWqU2o0AAPC4WIMFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEA\nACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYQLVIkuTn59egQYOSkpKK\nxhQXF9va2s6ZM6e0Eh0dHRAQUNMxAABjR8ACqmX79u1CCEmSsrKyKhpz+vRpvV6/cePGgoICuZKU\nlOTt7S2EOHTokK2tbXFxcSVjhBAJCQl+fn5NmjRp27btzp076/aQAAB1hoAFVK2oqOidd9759NNP\n3d3dL1y4UNGwpKSkUaNGOTs7R0dHl1a8vLyEEJ07d46Li7O0tKxkzLVr10aMGLFgwYJbt26FhYVN\nnz697o8MAFAnCFgwI3l5eTk1pNfrhRDr16/38/Pz8vLS6XSlAevQoUPdunVr2rTp/Pnz5UpycrKP\nj8+0adPWrFkjhJAkKTk5WT47tXTp0qSkpMrHxMXFBQUFDRs2TKvV9ujRw8LCQo0/JACAAizVbgCo\nP0899ZROp7O0rMFf+5kzZ3p5eX388cc//vijEMLDw0MOWJIkvfTSS1u2bOnSpctzzz23ePFi+ezU\nmDFjevTo8c4775w9e9ba2lqj0bi7uwshEhMTJ0+eLISoZExgYGBgYKAQIi0tbebMmStWrKiLPwQA\nQD0gYMGMtG7dOjw83NHRsUZbLVmy5OzZs82bN5ffTpw4UX7RqlWrpUuXBgYGxsTEWFpaFhUVnTx5\n0svLq3Hjxi+88MLatWv79Onj7e2t0Wjy8/PPnTvXtWvXSsbIc96+fXvBggUpKSmrV6/u1q2bgscO\nAKhPXCIEKnPjxo3169ffvXtXkiRJkjZs2CCfwdJoNImJiXPnzo2Pj+/evbskSb/99luLFi0aN24s\nhJg+ffqGDRvi4+PlxVUpKSndu3e3srKqZIwQYs+ePaNHjx48eHB8fDzpCgCMGgELqMySJUuCg4Md\nHBzkty1atJADVuPGjX/44YfBgwcPGzbMyspKo9GULqUSQjz11FMuLi5ffvmlXElMTPT19RVCVDIm\nNjZ20aJF0dHRQ4YMycvLy8vLq/+DBQAohYAFVCgtLS0qKio0NLS00rJly6tXrz548GDx4sVvvfWW\ns7PzypUrN2/eLMp8GVAIodFoXn311QcPHsiV0oBVyZjly5cnJSW5uro6ODg4ODgMGzasng8WAKAg\njSRJavcACCHERx991L59+1GjRtXdLkaPHh0REVHTNVgAANQUZ7AAAAAURsACAABQGAELAABAYQQs\nAAAAhRGwAAAAFEbAAgAAUBgBCwAAQGEELAAAAIURsAAAABRGwAIAAFAYAQsAAEBhBCwAAACFEbAA\nAAAUZql2A4ACdu3aNWPOX5s0a17lyKeGj6zoowd5eeNHPbdkyRJFWwMAmCMCFkzBvXv3RkyaOiTw\nlceZ5OLpU5lxe5VqCQBgzrhECAAAoDACFgAAgMIIWDBBRYWFrw32zb11o7SSm31zUXDgpN6d358R\nlH/vbiVFAAAeHwELpmbPxnVhr4y5kZlRthjx/mKdZ+d1h3+xsbX77otPKikCAPD4WOQOQ3Hnzp34\n+PiHDx/WYtujR48Kl1bya13Hzs3cdf94dVLpp5Ik/Xxg3ztrI62srX2HPLv5k+UT54b9sThgxBhl\njgQAYPYIWDAUvr6+p06dysnJqcW2eXl5Vi7/fd3Vt98jnxbk5xXk57Vs214I0apdh5wb1ysqAgCg\nCAIWDMXw4cOHDx9eu20bNmwYdz6z8jEajUYIIUlSiV5feREAgMfEGiyYvgZ29g3s7K9ePC+EuJZ+\nwbm5W0VFAAAUQcCC6dNoND6DApJ/PCBJ0vG4WJ/BwyoqAgCgCAIWzELQ3xaePZ40baDX3dvZ416b\nVUkRAIDHp5EkSe0egMcVFRX15ZboJ7x9H2eSW9eutm6g+eijj5TqCgBgtljkDlPg7+9vZWVV5bAV\nK1bMnDnTzs6u/I87eXTp0kXhzgAAZokzWDAjo0ePjoiIcHR0VLsRAICJYw0WAACAwghYAAAACiNg\nAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAF\nAACgMAIWAACAwghYQLVIkuTn59egQYOSkpLKRyYkJIwdO7Zp06bu7u5/+ctf8vPzhRAlJSWOjo6n\nT5+Wx5w7d65jx45ffPFFnfcNAFADAQuolu3btwshJEnKysqqZFhkZGRwcHBISMjFixd/+OGHhISE\nefPmCSHS0tKKi4s7duwohEhMTAwICFi5cuVrr7129+7dKVOmNGvWrEWLFkuWLKmfYwEA1DUCFlC1\noqKid95559NPP3V3d79w4UJFwzIyMubNm3fo0KHhw4fb29t7enquWrVqy5YtQoikpKSePXtaWFjs\n2LHjlVdeiY6OHj58uBDilVde0ev1p06dOnz48Lvvvlvl6TEAgFEgYAFVW79+vZ+fn5eXl06nKw1Y\nhw4d6tatW9OmTefPny9Xli5dGhoa2rx589INPT09s7OzS0pKkpOTvby8Vq1aNWPGjNjY2J49ewoh\niouLQ0JCPv/8c71ev2fPHh8fH62W/yUBwBRYqt0AUH8yMzNnz55tbW1d/U1effXVDh06fPzxxz/+\n+KMQwsPDQw5YkiS99NJLW7Zs6dKly3PPPbd48WJLS8vDhw/Pnj277OYXL17s0KGDVqtNSkrKzs6+\ncOFCTk6OXq+XP7W0tBw9enRKSoqXl5cQIiYmRrFDBQCoSiNJkto9APXkt99+KygoqP54jUbTsWPH\nDz/8cOHChaXFiRMnRkRESJLUu3dvJyenwMDAcePGOTk5lZSUWFtb3717187OrnTw8uXLz5079+WX\nXzo6Ok6bNm3FihUTJ050cXH5+OOPS8dIknTnzp3w8PCPPvooIyNDkSMFAKiLgAVU5saNGz4+PqdO\nnXJwcBBCRERErFu3Li3hXrYAAAOdSURBVC4uTgih1+sPHjwYGRkZGxubnp6u0Wjc3NxiYmK6d+8u\nb3vt2jVvb+/Dhw8/fPiwd+/eubm5Wq326NGjgwcPzsjIaNy48YQJE9atW2drayuEuHDhwjPPPFPJ\nAi8AgBFhwQdQmSVLlgQHB8vpSgjRokULOQM1btz4hx9+GDx48LBhw6ysrDQajRBi/PjxM2bMOH/+\n/P379/fs2TNgwIBFixa1adMmOTm5V69e8voqX1/fLl26rF69Wgixf//+zZs36/X6zMzMmTNnhoaG\nqnegAAAlcQYLqFBaWlrv3r3T0tKaNGkiV06fPt2lS5f8/Py1a9d+9dVXly5d8vT0/Oyzz3x8fIQQ\n+fn58+bN+/bbb4uLi/38/EJDQ/39/YUQc+bM0Wq1H374oTxJZGTkG2+8kZ6evmvXrjfffPPOnTvu\n7u6hoaHBwcFyUAMAGDsCFgAAgMK4RAgAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACA\nwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAK\nI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiM\ngAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDAC\nFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghY\nAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2AB\nAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUA\nAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAA\ngMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAA\nCiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAo\njIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoLD/B2pEkaRaznKOAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -h 800 -w 800\n",
    "library(ape)\n",
    "tre <- read.tree(\"empirical_8/RAxML_bipartitions.empirical_8\")\n",
    "ltre <- ladderize(tre)\n",
    "\n",
    "par(mfrow=c(1,2))\n",
    "plot(ltre, use.edge.length=F)\n",
    "nodelabels(ltre$node.label)\n",
    "\n",
    "plot(ltre, type='u')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get phylo distances (GTRgamma dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 0.05226094\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "mean(cophenetic.phylo(ltre))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}