{
 "metadata": {
  "name": "",
  "signature": "sha256:92294c5c0c897aecf6b97a46362eb0334053fbb5f5a82903db02b17af81288bc"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Example PE-GBS (w/ merged reads) notebook\n",
      "\n",
      "Paired-end GBS (or similarly, paired Ez-RAD) data may commonly contain first and second reads that overlap to some degree and therefore should be merged. This notebook will outline how to analyze such data in _pyRAD_ (v.3.03 or newer). I will demonstrate this example using an empirical data set of mine. To begin, we will create an empty directory within our current directory in which to perform the assembly."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "mkdir -p analysis_pyrad"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we create an empty params file by calling pyRAD with the -n option."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "\tnew params.txt file created\n"
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Enter params file arguments"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then we enter some arguments into the params file, which can be done using any text editor. I use the script below to automate it here. Some important arguments for dealing with PE-GBS data in particular include the following:  \n",
      "\n",
      "+ Set the correct datatype. In the case of this analysis, we will use 'pairgbs' for de-multiplexing in step 1, but then switch it to 'merged' for steps 2-7 after we merge the paired reads.  \n",
      "\n",
      "+ The filter setting (param 21) is not needed in this case because we will be using an external read merging program that applies this filtering. \n",
      "\n",
      "+ We will, however, still probably want to trim the overhanging edges of reads (param 29), which helps to remove barcodes or adapters that were missed. \n",
      "\n",
      "+ It can also be helpful with this data type to use a low mindepth setting in conjunction with majority rule consensus base calls (param 31). This will tell pyrad to make statistical base calls for all sites with depth >mindepth, but to make majority rule base calls at sites with depth lower than mindepth, but higher than the majority rule minimum depth. I first run the analysis with the option turned off, and then at the end of the tutorial show the difference with it turned on. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "sed -i '/## 1. /c\\analysis_pyrad/      ## 1. working directory   ' params.txt\n",
      "sed -i '/## 2. /c\\/home/deren/Longiflorae/*.gz   ## 2. working directory     ' params.txt\n",
      "sed -i '/## 3. /c\\/home/deren/barcodes/cyatho_longi.barcodes  ## 3. barcodes ' params.txt\n",
      "sed -i '/## 7. /c\\12                   ## 7. N processors   ' params.txt\n",
      "sed -i '/## 9. /c\\6                    ## 9. NQual     ' params.txt\n",
      "sed -i '/## 10. /c\\.85                 ## 10. Wclust     ' params.txt\n",
      "sed -i '/## 11. /c\\pairgbs             ## 11. datatype   ' params.txt\n",
      "sed -i '/## 14. /c\\c85d6m4p3           ## 14. output prefix  ' params.txt\n",
      "sed -i '/## 29./c\\1,1                  ## 29. trim edges ' params.txt\n",
      "sed -i '/## 30./c\\p                    ## 30. output formats ' params.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 83
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Let's view the params file"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cat params.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "==** parameter inputs for pyRAD version 3.0.3  **======================== affected step ==\r\n",
        "analysis_pyrad/      ## 1. working directory   \r\n",
        "/home/deren/Longiflorae/*.gz   ## 2. working directory     \r\n",
        "/home/deren/barcodes/cyatho_longi.barcodes  ## 3. barcodes \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",
        "TGCAG                       ## 6. Restriction overhang (e.g., C|TGCAG -> TGCAG)     (s1,s2)\r\n",
        "12                   ## 7. N processors   \r\n",
        "6                           ## 8. Mindepth: min coverage for a cluster              (s4,s5)\r\n",
        "6                    ## 9. NQual     \r\n",
        ".85                 ## 10. Wclust     \r\n",
        "pairgbs             ## 11. datatype   \r\n",
        "4                           ## 12. MinCov: min samples in a final locus             (s7)\r\n",
        "3                           ## 13. MaxSH: max inds with shared hetero site          (s7)\r\n",
        "c85d6m4p3           ## 14. output prefix  \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",
        "                       ## 18.opt.: loc. of de-multiplexed data                      (s2)\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",
        "1,1                  ## 29. trim edges \r\n",
        "p                    ## 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"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step 1: De-multiplex reads"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "\n",
        "\n",
        "     ------------------------------------------------------------\n",
        "      pyRAD : RADseq for phylogenetics & introgression analyses\n",
        "     ------------------------------------------------------------\n",
        "\n",
        "\n",
        "\tstep 1: sorting reads by barcode\n",
        "\t ..........................................................."
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Let's look at the first 100 lines of the stats output for demultiplexing\n",
      "It looks like most reads from each file were matched to a barcode."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "head -n 100 analysis_pyrad/stats/s1.sorting.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "file    \tNreads\tcut_found\tbar_matched\n",
        "lane2_NoIndex_L002_049.fastq.gz\t4000000\t3721343\t3563126\n",
        "lane2_NoIndex_L002_039.fastq.gz\t4000000\t3775744\t3746367\n",
        "lane2_NoIndex_L002_038.fastq.gz\t4000000\t3749702\t3712199\n",
        "lane2_NoIndex_L002_008.fastq.gz\t4000000\t3788251\t3763415\n",
        "lane2_NoIndex_L002_030.fastq.gz\t4000000\t3756116\t3721905\n",
        "lane2_NoIndex_L002_042.fastq.gz\t4000000\t3821669\t3761365\n",
        "lane2_NoIndex_L002_027.fastq.gz\t4000000\t3809029\t3786687\n",
        "lane2_NoIndex_L002_046.fastq.gz\t4000000\t3800963\t3717126\n",
        "lane2_NoIndex_L002_037.fastq.gz\t4000000\t3746059\t3721532\n",
        "lane2_NoIndex_L002_045.fastq.gz\t4000000\t3698917\t3593579\n",
        "lane2_NoIndex_L002_019.fastq.gz\t4000000\t3778683\t3750877\n",
        "lane2_NoIndex_L002_044.fastq.gz\t4000000\t3751756\t3661419\n",
        "lane2_NoIndex_L002_058.fastq.gz\t4000000\t3755612\t3729931\n",
        "lane2_NoIndex_L002_032.fastq.gz\t4000000\t3823466\t3804376\n",
        "lane2_NoIndex_L002_025.fastq.gz\t4000000\t3821273\t3800999\n",
        "lane2_NoIndex_L002_055.fastq.gz\t4000000\t3789040\t3767778\n",
        "lane2_NoIndex_L002_023.fastq.gz\t4000000\t3832294\t3814024\n",
        "lane2_NoIndex_L002_014.fastq.gz\t4000000\t3832280\t3812920\n",
        "lane2_NoIndex_L002_013.fastq.gz\t4000000\t3835097\t3816100\n",
        "lane2_NoIndex_L002_031.fastq.gz\t4000000\t3827911\t3808997\n",
        "lane2_NoIndex_L002_011.fastq.gz\t4000000\t3845888\t3828585\n",
        "lane2_NoIndex_L002_022.fastq.gz\t4000000\t3837573\t3819811\n",
        "lane2_NoIndex_L002_024.fastq.gz\t4000000\t3829959\t3810628\n",
        "lane2_NoIndex_L002_034.fastq.gz\t4000000\t3611235\t3590798\n",
        "lane2_NoIndex_L002_056.fastq.gz\t4000000\t3802234\t3779975\n",
        "lane2_NoIndex_L002_026.fastq.gz\t4000000\t3819201\t3799143\n",
        "lane2_NoIndex_L002_021.fastq.gz\t4000000\t3835070\t3816075\n",
        "lane2_NoIndex_L002_052.fastq.gz\t4000000\t3828882\t3810331\n",
        "lane2_NoIndex_L002_041.fastq.gz\t4000000\t3835218\t3791068\n",
        "lane2_NoIndex_L002_051.fastq.gz\t4000000\t3833143\t3814474\n",
        "lane2_NoIndex_L002_047.fastq.gz\t4000000\t3759517\t3607421\n",
        "lane2_NoIndex_L002_059.fastq.gz\t4000000\t3757946\t3719770\n",
        "lane2_NoIndex_L002_054.fastq.gz\t4000000\t3819129\t3798453\n",
        "lane2_NoIndex_L002_028.fastq.gz\t4000000\t3795559\t3771229\n",
        "lane2_NoIndex_L002_048.fastq.gz\t4000000\t3748548\t3623691\n",
        "lane2_NoIndex_L002_035.fastq.gz\t4000000\t3266305\t3243997\n",
        "lane2_NoIndex_L002_043.fastq.gz\t4000000\t3817002\t3744305\n",
        "lane2_NoIndex_L002_057.fastq.gz\t4000000\t3743088\t3719312\n",
        "lane2_NoIndex_L002_018.fastq.gz\t4000000\t3794176\t3769286\n",
        "lane2_NoIndex_L002_017.fastq.gz\t4000000\t3807601\t3784147\n",
        "lane2_NoIndex_L002_050.fastq.gz\t4000000\t3745707\t3651533\n",
        "lane2_NoIndex_L002_033.fastq.gz\t4000000\t3821688\t3802696\n",
        "lane2_NoIndex_L002_040.fastq.gz\t4000000\t3761738\t3722996\n",
        "lane2_NoIndex_L002_053.fastq.gz\t4000000\t3822961\t3803422\n",
        "lane2_NoIndex_L002_029.fastq.gz\t4000000\t3781432\t3755276\n",
        "lane2_NoIndex_L002_016.fastq.gz\t4000000\t3822901\t3802934\n",
        "lane2_NoIndex_L002_036.fastq.gz\t4000000\t3805280\t3783717\n",
        "lane2_NoIndex_L002_020.fastq.gz\t4000000\t3732223\t3681151\n",
        "lane2_NoIndex_L002_005.fastq.gz\t4000000\t3814967\t3793880\n",
        "lane2_NoIndex_L002_006.fastq.gz\t4000000\t3812538\t3792209\n",
        "lane2_NoIndex_L002_001.fastq.gz\t4000000\t3838930\t3821090\n",
        "lane2_NoIndex_L002_007.fastq.gz\t4000000\t3799427\t3776091\n",
        "lane2_NoIndex_L002_015.fastq.gz\t4000000\t3822598\t3802203\n",
        "lane2_NoIndex_L002_004.fastq.gz\t4000000\t3823816\t3803572\n",
        "lane2_NoIndex_L002_010.fastq.gz\t4000000\t3761201\t3730496\n",
        "lane2_NoIndex_L002_009.fastq.gz\t4000000\t3777890\t3750566\n",
        "lane2_NoIndex_L002_002.fastq.gz\t4000000\t3826291\t3807448\n",
        "lane2_NoIndex_L002_003.fastq.gz\t4000000\t3824097\t3804683\n",
        "lane2_NoIndex_L002_012.fastq.gz\t4000000\t3838714\t3820626\n",
        "\n",
        "\n",
        "sample\ttrue_bar\tobs_bars\tN_obs\n",
        "d31733    \tAAAAGTT    \tAAAAGTT\t5114026    \n",
        "d31733    \tAAAAGTT    \tAANAGTT\t813193     \n",
        "d31733    \tAAAAGTT    \tAAGAGTT\t90500      \n",
        "d31733    \tAAAAGTT    \tAACAGTT\t57040      \n",
        "d31733    \tAAAAGTT    \tAAAAGGT\t19502      \n",
        "d31733    \tAAAAGTT    \tAAAAGAT\t19214      \n",
        "d31733    \tAAAAGTT    \tAATAGTT\t13373      \n",
        "d31733    \tAAAAGTT    \tAAAGGTT\t9539       \n",
        "d31733    \tAAAAGTT    \tAAAAGCT\t8946       \n",
        "d31733    \tAAAAGTT    \tGAAAGTT\t8772       \n",
        "d31733    \tAAAAGTT    \tNAAAGTT\t8297       \n",
        "d31733    \tAAAAGTT    \tAAAAGTC\t7090       \n",
        "d31733    \tAAAAGTT    \tAAAAATT\t6140       \n",
        "d31733    \tAAAAGTT    \tAGAAGTT\t6042       \n",
        "d31733    \tAAAAGTT    \tAAAATTT\t4801       \n",
        "d31733    \tAAAAGTT    \tCAAAGTT\t4089       \n",
        "d31733    \tAAAAGTT    \tAAAAGTG\t3887       \n",
        "d31733    \tAAAAGTT    \tAAACGTT\t3255       \n",
        "d31733    \tAAAAGTT    \tACAAGTT\t2946       \n",
        "d31733    \tAAAAGTT    \tTAAAGTT\t1997       \n",
        "d31733    \tAAAAGTT    \tAAAAGTA\t1857       \n",
        "d31733    \tAAAAGTT    \tATAAGTT\t1145       \n",
        "d31733    \tAAAAGTT    \tAAATGTT\t994        \n",
        "d31733    \tAAAAGTT    \tAAAACTT\t896        \n",
        "d31733    \tAAAAGTT    \tANAAGTT\t512        \n",
        "d33291    \tAACGCCT    \tAACGCCT\t949865     \n",
        "d33291    \tAACGCCT    \tAANGCCT\t150664     \n",
        "d33291    \tAACGCCT    \tAAGGCCT\t18149      \n",
        "d33291    \tAACGCCT    \tAAAGCCT\t14152      \n",
        "d33291    \tAACGCCT    \tAATGCCT\t3499       \n",
        "d33291    \tAACGCCT    \tAACGACT\t2636       \n",
        "d33291    \tAACGCCT    \tAACGCTT\t2091       \n",
        "d33291    \tAACGCCT    \tAACGTCT\t2053       \n",
        "d33291    \tAACGCCT    \tGACGCCT\t1775       \n",
        "d33291    \tAACGCCT    \tAACGCCC\t1714       \n",
        "d33291    \tAACGCCT    \tAGCGCCT\t1543       \n",
        "d33291    \tAACGCCT    \tNACGCCT\t1474       \n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Merge paired-reads\n",
      "For this we will use the program [__PEAR__](https://github.com/xflouris/PEAR). This step is very important. I can pretty much assure you that you will have many paired end reads in your library that overlap, and if you do not merge these reads they can really mess up your analysis. In this case the vast majority of paired reads are overlapping (~80%), and in fact, the merged reads are the only portion of the data that we will be using for analysis.  \n",
      "\n",
      "The first thing we have to do is unzip the reads if they are gzipped. If you demultiplexed your data in pyrad they will be in the fastq/ directory. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "## For this example I subselect only the samples from \n",
      "## this library that start with the letter 'd'\n",
      "gunzip analysis_pyrad/fastq/d*.gz"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I prefer using a stringent filter setting, including the `-q` option which trims parts of the read that are low quality. This often results in removing all reads that are not merged. You can assemble merged and non-merged reads separately, as described in the 'non-merged PE-GBS tutorial', but __here we will focus on just analyzing the merged reads.__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "for gfile in analysis_pyrad/fastq/*_R1.fq;\n",
      "    do pear -f $gfile \\\n",
      "            -r ${gfile/_R1.fq/_R2.fq} \\\n",
      "            -o ${gfile/_R1.fq/} \\\n",
      "            -n 33 \\\n",
      "            -t 33 \\\n",
      "            -q 10 \\\n",
      "            -j 20  >> pear.log 2>&1;\n",
      "done"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Set the data location to the de-multiplexed & merged data "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "## set location of demultiplexed data that are 'pear' filtered\n",
      "sed -i '/## 18./c\\analysis_pyrad/fastq/*.assembled.fastq ## 18. demulti data loc ' ./params.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Set the datatype to \"merged\""
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "sed -i '/## 11./c\\merged            ## 11. data type ' ./params.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step 2: Quality filtering"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "\n",
        "\n",
        "     ------------------------------------------------------------\n",
        "      pyRAD : RADseq for phylogenetics & introgression analyses\n",
        "     ------------------------------------------------------------\n",
        "\n",
        "\tsorted .fastq from analysis_pyrad/fastq/*.assembled.fastq being used\n",
        "\tstep 2: editing raw reads \n",
        "\t........................"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Let's take a look at the results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "cat analysis_pyrad/stats/s2.rawedit.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "sample \tNreads\tpassed\tpassed.w.trim\tpassed.total\n",
        "d19long1.assembled\t2973258\t2865913\t0\t2865913\n",
        "d30181.assembled\t4178962\t4048613\t0\t4048613\n",
        "d30695.assembled\t4435144\t4249855\t0\t4249855\n",
        "d31733.assembled\t4757784\t4598646\t0\t4598646\n",
        "d33291.assembled\t847503\t819072\t0\t819072\n",
        "d34041.assembled\t988252\t969387\t0\t969387\n",
        "d35178.assembled\t2479633\t2429555\t0\t2429555\n",
        "d35320.assembled\t1308443\t1269046\t0\t1269046\n",
        "d35371.assembled\t572512\t559953\t0\t559953\n",
        "d35422.assembled\t3069039\t3016381\t0\t3016381\n",
        "d39103.assembled\t2629656\t2523785\t0\t2523785\n",
        "d39104.assembled\t3415956\t3353178\t0\t3353178\n",
        "d39114.assembled\t2215501\t2133972\t0\t2133972\n",
        "d39187.assembled\t907321\t875889\t0\t875889\n",
        "d39253.assembled\t5665909\t5444296\t0\t5444296\n",
        "d39404.assembled\t2567688\t2447216\t0\t2447216\n",
        "d39531.assembled\t3254639\t3128372\t0\t3128372\n",
        "d39968.assembled\t3561187\t3441543\t0\t3441543\n",
        "d40328.assembled\t1651265\t1601896\t0\t1601896\n",
        "d41058.assembled\t2444960\t2395622\t0\t2395622\n",
        "d41237.assembled\t1847033\t1815168\t0\t1815168\n",
        "d41389.assembled\t2491305\t2433467\t0\t2433467\n",
        "d41732.assembled\t2344042\t2291187\t0\t2291187\n",
        "decor21.assembled\t5056601\t4926349\t0\t4926349\n",
        "\n",
        "    Nreads = total number of reads for a sample\n",
        "    passed = retained reads that passed quality filtering at full length\n",
        "    passed.w.trim= retained reads that were trimmed due to detection of adapters\n",
        "    passed.total  = total kept reads of sufficient length\n",
        "    note: you can set the option in params file to include trimmed reads of xx length. \n"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step 3: Clustering"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Let's take a look at the results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "cat analysis_pyrad/stats/s3 "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### and let's look at an example cluster file"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "gzip -cd analysis_pyrad/clust.85/d35422.assembled.clustS.gz | head -n 100 | cut -b 1-80"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ">d35422.assembled_2120_r1;size=51;\n",
        "TGCAGCGCTACCTGCCCGATCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_274175_r1;size=15;+\n",
        "TGCAGCGCTACCTGCCCGATCTGTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_2265219_r1;size=2;+\n",
        "TGCAGCGCTACCTGCCGGATCTGTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_216079_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGATCTGCTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_2495995_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGGTCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_590416_r1;size=1;+\n",
        "TGCAGCGCTATCTGCCCGATCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_1947641_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGATCTTTTCGAGAAGCGGACGAGCTGAGAATGAACCTGCA\n",
        ">d35422.assembled_1601389_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCGGATCTGTTCGAGAAGCGGACGAACTGAGAATGACACTGCA\n",
        ">d35422.assembled_190367_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGATCCTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        ">d35422.assembled_2631365_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGATCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCG\n",
        ">d35422.assembled_126737_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCCGATCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTNCA\n",
        ">d35422.assembled_288692_r1;size=1;+\n",
        "TGCAGCGCTACCTGCCGGATCNGTTCGAGAAGCGGACGAACTGAGAATGACACTGCA\n",
        ">d35422.assembled_2194406_r1;size=1;+\n",
        "TGCAGCGCTANCTNCCCGATCTTTTCGAGAAGCGGANNANCTGAGAATGAACCTGCA\n",
        ">d35422.assembled_831862_r1;size=1;+\n",
        "TGCAGCGCTACCTGNCCGATCTTTTCGAGAAGCGGACGAACTGAGAATGAACCTGCA\n",
        "//\n",
        "//\n",
        ">d35422.assembled_61639_r1;size=24;\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCCCGCGACGCCCTGGTGCCGGCC\n",
        ">d35422.assembled_825340_r1;size=8;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGTCAGTTCATCCAGTCGCATCCGCGCGACGCCCTGGTGCCGTCC\n",
        ">d35422.assembled_88436_r1;size=2;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCCCGCGACGCCCTGGTGCCGGCC\n",
        ">d35422.assembled_354389_r1;size=2;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCCCGCGACGCCCTGGTGCCGGCC\n",
        ">d35422.assembled_2229078_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTACGTCAGTTCATCCAGTCGCATCCGCGCGACTCCCTGGTGCCGTCC\n",
        ">d35422.assembled_1405977_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATNCAGTCGCANCCCCGCGACGCCCTGGTGCCGNCC\n",
        ">d35422.assembled_2738651_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCCCGCGACGCCCTGGTGCCGGCC\n",
        ">d35422.assembled_2565401_r1;size=1;+\n",
        "TGCAGCGTCAGTACGAGCTGGCGGAGATGGGCCTGCGTCAGTTCATCCAGTCCCATCCGCGTGATGCACTGGTACCGGCC\n",
        ">d35422.assembled_2060896_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCGCGCGACGCCCTCGTGCCGGCC\n",
        ">d35422.assembled_27774_r1;size=1;+\n",
        "TGCAGCGTCAGTACGAGCTGGCGGAGATGGGCCTGCGTCAGTTCATCCAGTCCCATCCGCGTGATGCACTGGTACCGGCC\n",
        ">d35422.assembled_2776613_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGAGATGGGCCTGCGCCAGTTCATCCAGTCGCATCCCCGCGACGCCCTGGTGCCGGCC\n",
        ">d35422.assembled_2566934_r1;size=1;+\n",
        "TGCAGCGCCAGTACGAGCAGGCCGGGATGGGCCTGCGTCAGTTCATCCAGTCGCATCCGCGCGACGCCCTGGTGCCGTCC\n",
        "//\n",
        "//\n",
        ">d35422.assembled_781655_r1;size=2;+\n",
        "TGCAGGACCAGGGCTTCGCCGTGACCGCGAACATGGTGACCAGCCGCATCGTCGTGGACGCCGTGGGCACCTTCGACGCC\n",
        ">d35422.assembled_1218064_r1;size=2;\n",
        "TGCAGGACCAGGGCTTCGCCGTGACCGCGAACATGGTGACCAGTCGTGCTGTGGTCGAAGCCGTCGGCGCCTTCGACGCC\n",
        ">d35422.assembled_81417_r1;size=1;+\n",
        "TGCAGGACCAGGGCTTCGCCGTGNCCGCGAACATGGTNACCACCCGCGCCGTTGTCGACGCCGTCGGCACCTTCGACGCC\n",
        ">d35422.assembled_1447392_r1;size=1;+\n",
        "TGCAGGACCAGGGCTTCGCNGTGACCNCGAACATGGTGACCAGCCGCATCGTCGTGGACGCCGTGGGCACCTTCGACGCC\n",
        ">d35422.assembled_1916059_r1;size=1;+\n",
        "TGCAGGACCAGGGCTTCGCCGTGACCGCGAACATGGTGACCAGCCGCGTCGTCGTGGACGCCGTGGGCACCTTCGACGCC\n",
        "//\n",
        "//\n",
        ">d35422.assembled_1086507_r1;size=1;\n",
        "TGCAGNNANC-AAAAAAAATATATG-TTTTTTTTTTTTGCTGCA\n",
        ">d35422.assembled_1935566_r1;size=1;+\n",
        "TGCAGNTAACNAAAAAAAATATATGTTTTTTTTTTTTTGCTGCA\n",
        "//\n",
        "//\n",
        ">d35422.assembled_105061_r1;size=2;\n",
        "TGCAGGTCTCCTCCGACGTCTACTACTACACGGTCGGCTCGAAGCTCTTCTCGAGGCCGAACAAGCCCCTGCA\n",
        ">d35422.assembled_405948_r1;size=2;+\n",
        "TGCAGGTCTCCTCCGACGTCTACTACTACACGGTCGGCGCGAACCTCTTCGAGAGGAAGAACAAGCCGCTGCA\n",
        ">d35422.assembled_1112236_r1;size=1;+\n",
        "TGCAGGTNNCCTCNGACGTCTACTACTACACGGTCGGCGCGAACCTCTTCTCCAAGCCGAACAAGCCGCTGCA\n",
        ">d35422.assembled_129605_r1;size=1;+\n",
        "TGCAGGTGTCCTCCGACGTCTACTACTACACCGTCGGGTCGAAGCTCTTCTCGAGGCCGAACAAGCCCCTGCA\n",
        ">d35422.assembled_108168_r1;size=1;+\n",
        "TGCAGGTGTCCTCCGACGTCTACTACTACACGGTCGGCGCGAACCTCTTCTCGAAGCCGAACAAGCCCCTGCA\n",
        ">d35422.assembled_1796775_r1;size=1;+\n",
        "TGCAGGTGTCCTCCGACGTCTACTACTACACGGTCGGCGCGAACCTCTTCTCCAAGCCGAACAAGCCGCTGCA\n",
        "//\n",
        "//\n",
        ">d35422.assembled_2015519_r1;size=1;+\n",
        "TGCAGGGGGNCGCGGGGAGCAACCCCAGCTCGGTGCACGTGAGACAGGTGTCGATCACGGGCGGGGACCTCCTGGGGCCT\n",
        ">d35422.assembled_1326334_r1;size=1;\n",
        "TGCAGGGGGTCGCGGGGAGCAACCCCAGCTCGGTGCACGTGAGACAGGTGTCGATCACGGGCGGGGACCTCCTGGGGCCT\n",
        "//\n",
        "//\n",
        ">d35422.assembled_735556_r1;size=2;\n",
        "TGCAGCGCCGCCAGCGGCTGCAGAACTGCCCCCGAGTTTCCAGCTTCTGCA\n",
        ">d35422.assembled_2203593_r1;size=1;+\n",
        "TGCAGCGCCTGGACGCCCTGCAGAACTGCCCCCGAGTTTCCAGCTTCTGCA\n",
        "//\n",
        "//\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "\n",
        "gzip: stdout: Broken pipe\n"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step 4: Parameter estimation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "cat analysis_pyrad/stats/Pi_E_estimate.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "taxa\tH\tE\n",
        "d35371.assembled\t0.0182707\t0.00219242\t\n",
        "d39187.assembled\t0.02781985\t0.00272968\t\n",
        "d39104.assembled\t0.02073342\t0.0015139\t\n",
        "d35320.assembled\t0.02621882\t0.00175474\t\n",
        "d34041.assembled\t0.02897603\t0.00381241\t\n",
        "d33291.assembled\t0.02675046\t0.00311296\t\n",
        "d40328.assembled\t0.02060546\t0.00267209\t\n",
        "d35422.assembled\t0.04092042\t0.00224104\t\n",
        "d19long1.assembled\t0.02280653\t0.00095164\t\n",
        "d41389.assembled\t0.02985094\t0.00193657\t\n",
        "d41237.assembled\t0.05196614\t0.00204698\t\n",
        "d39114.assembled\t0.02263828\t0.00168706\t\n",
        "d39404.assembled\t0.04004115\t0.00165295\t\n",
        "d41058.assembled\t0.03910496\t0.00179648\t\n",
        "d39531.assembled\t0.03286906\t0.0014618\t\n",
        "d41732.assembled\t0.0367377\t0.00251955\t\n",
        "d35178.assembled\t0.0248532\t0.00261068\t\n",
        "d31733.assembled\t0.03005955\t0.00205113\t\n",
        "d30181.assembled\t0.02449057\t0.00222622\t\n",
        "d39968.assembled\t0.03075244\t0.00329168\t\n",
        "d39253.assembled\t0.03181217\t0.00343958\t\n",
        "decor21.assembled\t0.03917313\t0.00253141\t\n",
        "d30695.assembled\t0.02899983\t0.00255832\t\n",
        "d39103.assembled\t0.03842947\t0.00396182\t\n"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step 5: Consenus base calling\n",
      "\n",
      "You may want to adjust:\n",
      " + max cluster size (param 33)\n",
      " + max number of Ns\n",
      " + max number of Hs"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "cat analysis_pyrad/stats/s5.consens.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "taxon          \tnloci\tf1loci\tf2loci\tnsites\tnpoly\tpoly\n",
        "d35371.assembled\t60885\t7839\t6432\t489780\t2313\t0.0047225\n",
        "d41237.assembled\t293759\t10550\t5791\t433633\t2696\t0.0062172\n",
        "d35320.assembled\t156419\t12616\t9039\t777103\t3193\t0.0041089\n",
        "d41058.assembled\t329719\t13087\t8587\t620904\t3583\t0.0057706\n",
        "d35422.assembled\t169511\t13976\t8440\t588597\t3865\t0.0065665\n",
        "d19long1.assembled\t161936\t11666\t8611\t763559\t2339\t0.0030633\n",
        "d33291.assembled\t169929\t14920\t10442\t913599\t3976\t0.004352\n",
        "d34041.assembled\t197867\t19992\t14191\t931867\t6604\t0.0070868\n",
        "d39187.assembled\t96755\t12685\t8535\t752923\t3913\t0.0051971\n",
        "d39104.assembled\t101266\t8902\t6827\t546599\t1931\t0.0035328\n",
        "d41732.assembled\t452240\t23593\t15770\t1129545\t7020\t0.0062149\n",
        "d41389.assembled\t211387\t21245\t14809\t1053031\t6331\t0.0060122\n",
        "d39404.assembled\t243761\t23389\t13769\t1106220\t7083\t0.0064029\n",
        "d40328.assembled\t136690\t24752\t18911\t1544641\t7199\t0.0046606\n",
        "d39114.assembled\t218883\t25257\t18878\t1562653\t7444\t0.0047637\n",
        "d39531.assembled\t332792\t26717\t16869\t1393571\t7705\t0.005529\n",
        "d35178.assembled\t450842\t47869\t38116\t2467261\t21384\t0.0086671\n",
        "decor21.assembled\t577428\t46634\t28116\t2021202\t14304\t0.007077\n",
        "d31733.assembled\t439816\t40309\t27103\t2178840\t11820\t0.0054249\n",
        "d30181.assembled\t516149\t49126\t36637\t2879253\t13306\t0.0046213\n",
        "d39103.assembled\t642379\t48620\t27395\t2247944\t13442\t0.0059797\n",
        "d30695.assembled\t485553\t48503\t32034\t2797840\t15100\t0.005397\n",
        "d39253.assembled\t431026\t51267\t31529\t2654672\t14942\t0.0056286\n",
        "d39968.assembled\t456174\t61773\t39746\t3309319\t18475\t0.0055827\n",
        "\n",
        "    ## nloci = number of loci\n",
        "    ## f1loci = number of loci with >N depth coverage\n",
        "    ## f2loci = number of loci with >N depth and passed paralog filter\n",
        "    ## nsites = number of sites across f loci\n",
        "    ## npoly = number of polymorphic sites in nsites\n",
        "    ## poly = frequency of polymorphic sites\n"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Cluster across samples"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 6"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Assemble final data set"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "pyrad -p params.txt -s 7"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Examine results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "cat analysis_pyrad/stats/c88d6m4p3.stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "17401       ## loci with > minsp containing data\n",
        "16418       ## loci with > minsp containing data & paralogs removed\n",
        "16418       ## loci with > minsp containing data & paralogs removed & final filtering\n",
        "\n",
        "## number of loci recovered in final data set for each taxon.\n",
        "taxon\tnloci\n",
        "d19long1.assembled\t2232\n",
        "d30181.assembled  \t5279\n",
        "d30695.assembled  \t8028\n",
        "d31733.assembled  \t5717\n",
        "d33291.assembled  \t3529\n",
        "d34041.assembled  \t4291\n",
        "d35178.assembled  \t4252\n",
        "d35320.assembled  \t3451\n",
        "d35371.assembled  \t2535\n",
        "d35422.assembled  \t2776\n",
        "d39103.assembled  \t6100\n",
        "d39104.assembled  \t2311\n",
        "d39114.assembled  \t4695\n",
        "d39187.assembled  \t3249\n",
        "d39253.assembled  \t7324\n",
        "d39404.assembled  \t5203\n",
        "d39531.assembled  \t6176\n",
        "d39968.assembled  \t6724\n",
        "d40328.assembled  \t4921\n",
        "d41058.assembled  \t1872\n",
        "d41237.assembled  \t1308\n",
        "d41389.assembled  \t4676\n",
        "d41732.assembled  \t4007\n",
        "decor21.assembled \t7367\n",
        "\n",
        "\n",
        "## nloci = number of loci with data for exactly ntaxa\n",
        "## ntotal = number of loci for which at least ntaxa have data\n",
        "ntaxa\tnloci\tsaved\tntotal\n",
        "1\t-\n",
        "2\t-\t\t-\n",
        "3\t-\t\t-\n",
        "4\t5937\t*\t16418\n",
        "5\t3301\t*\t10481\n",
        "6\t2050\t*\t7180\n",
        "7\t1386\t*\t5130\n",
        "8\t910\t*\t3744\n",
        "9\t588\t*\t2834\n",
        "10\t412\t*\t2246\n",
        "11\t329\t*\t1834\n",
        "12\t241\t*\t1505\n",
        "13\t167\t*\t1264\n",
        "14\t140\t*\t1097\n",
        "15\t121\t*\t957\n",
        "16\t97\t*\t836\n",
        "17\t106\t*\t739\n",
        "18\t91\t*\t633\n",
        "19\t76\t*\t542\n",
        "20\t93\t*\t466\n",
        "21\t107\t*\t373\n",
        "22\t105\t*\t266\n",
        "23\t98\t*\t161\n",
        "24\t63\t*\t63\n",
        "\n",
        "\n",
        "## nvar = number of loci containing n variable sites (pis+autapomorphies).\n",
        "## sumvar = sum of variable sites (SNPs).\n",
        "## pis = number of loci containing n parsimony informative sites.\n",
        "## sumpis = sum of parsimony informative sites.\n",
        "\tnvar\tsumvar\tPIS\tsumPIS\n",
        "0\t2451\t0\t6616\t0\n",
        "1\t2053\t2053\t3490\t3490\n",
        "2\t1932\t5917\t1952\t7394\n",
        "3\t1655\t10882\t1254\t11156\n",
        "4\t1376\t16386\t865\t14616\n",
        "5\t1193\t22351\t611\t17671\n",
        "6\t959\t28105\t383\t19969\n",
        "7\t806\t33747\t266\t21831\n",
        "8\t660\t39027\t210\t23511\n",
        "9\t499\t43518\t184\t25167\n",
        "10\t429\t47808\t118\t26347\n",
        "11\t324\t51372\t110\t27557\n",
        "12\t266\t54564\t77\t28481\n",
        "13\t250\t57814\t70\t29391\n",
        "14\t214\t60810\t48\t30063\n",
        "15\t172\t63390\t39\t30648\n",
        "16\t171\t66126\t31\t31144\n",
        "17\t139\t68489\t24\t31552\n",
        "18\t104\t70361\t16\t31840\n",
        "19\t107\t72394\t10\t32030\n",
        "20\t79\t73974\t16\t32350\n",
        "21\t87\t75801\t10\t32560\n",
        "22\t59\t77099\t4\t32648\n",
        "23\t58\t78433\t2\t32694\n",
        "24\t43\t79465\t5\t32814\n",
        "25\t35\t80340\t2\t32864\n",
        "26\t36\t81276\t2\t32916\n",
        "27\t27\t82005\t0\t32916\n",
        "28\t33\t82929\t2\t32972\n",
        "29\t36\t83973\t1\t33001\n",
        "30\t23\t84663\t0\t33001\n",
        "31\t23\t85376\t0\t33001\n",
        "32\t13\t85792\t0\t33001\n",
        "33\t20\t86452\t0\t33001\n",
        "34\t12\t86860\t0\t33001\n",
        "35\t14\t87350\t0\t33001\n",
        "36\t10\t87710\t0\t33001\n",
        "37\t7\t87969\t0\t33001\n",
        "38\t8\t88273\t0\t33001\n",
        "39\t7\t88546\t0\t33001\n",
        "40\t3\t88666\t0\t33001\n",
        "41\t2\t88748\t0\t33001\n",
        "42\t2\t88832\t0\t33001\n",
        "43\t2\t88918\t0\t33001\n",
        "44\t4\t89094\t0\t33001\n",
        "45\t2\t89184\t0\t33001\n",
        "46\t4\t89368\t0\t33001\n",
        "47\t0\t89368\t0\t33001\n",
        "48\t2\t89464\t0\t33001\n",
        "49\t3\t89611\t0\t33001\n",
        "50\t2\t89711\t0\t33001\n",
        "51\t0\t89711\t0\t33001\n",
        "52\t1\t89763\t0\t33001\n",
        "53\t1\t89816\t0\t33001\n",
        "total var= 89816\n",
        "total pis= 33001\n"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "head -n 100 analysis_pyrad/outfiles/c88d6m4p3.loci | cut -b 1-88"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ">d31733.assembled      GTACAGCACCGGGTAGCGGCGGRGGCTAGCGGCGTAGCCGGGCGGCAGGTACACCCACACCCGGC\n",
        ">d39114.assembled      GTACAGCACCGGGTAGCGGCGGRGGCTAGCGGCGTAGCCGGGCGGCAGGTACACCCACACCCGGC\n",
        ">d39968.assembled      GTACAGCACCGGGTAGCGGCGGGGGCTAGCNGCGTAGCCGGGCGGCAGGTACACCCACACCCGGC\n",
        ">decor21.assembled     GTACAGCACCGGGTAGCGGCGGRGGCTAGCGGCGTAGCCGGGCGGCAGGTACACCCACACCCGGC\n",
        "//                                           *                                          \n",
        ">d30695.assembled      AATTGTGAGTACCGTCGTGACCGA-GGCACCGAGGC\n",
        ">d31733.assembled      AATTGTGAGTACCGTCGTGACCGANGGCACCGAGGC\n",
        ">d39103.assembled      AATTGTGAGTACCGTCGTGACCGA-GGCACCGAGGC\n",
        ">d39531.assembled      AATTGTGAGTACCGTCGTGACCGA-GGCACCGAGGC\n",
        ">d41389.assembled      AATTGTGAGTACCGTCGTRACCGA-GGCACCGAGGC\n",
        ">d41732.assembled      AATTGTGAGTACCGTCGTGACCGA-GGCACCGAGGC\n",
        ">decor21.assembled     AATTGTGAGTACCGTCGTGACCGA-GGCACCGAGGC\n",
        "//                                       -                 |\n",
        ">d30695.assembled      CCRTGGTGCGGCTGGCGCAAGGTGCGSGASCTSGCGGTGGCGCAGGGC\n",
        ">d31733.assembled      CCATGGTGCGGCTGGCGCAAGGTGCGSGACCTSGCGGTGGCGCAGGGC\n",
        ">d33291.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGAGCTCGCGGTGGCGCAGGGC\n",
        ">d34041.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGARCTCGCGGTGGCGCAGGGY\n",
        ">d35320.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGAGCTCGCGGTGNCGCAGGGC\n",
        ">d35422.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGAGCTCGCGGTGGCGCAGGGC\n",
        ">d39253.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGAGCTCGCGGTGGCGCAGGGY\n",
        ">d39531.assembled      CCGTGGTGCGGCTGGCGCAAGGTGCGCGAGCTCGCGGTGGCGCAGGGC\n",
        "//                       *                       *  *  *              *|\n",
        ">d30181.assembled      GCCCGGGAAGCACAAGGCTCCACGTTCCTTGGCCAGGGCATCGCCATTCCCCATGGTACGCCGCA\n",
        ">d30695.assembled      GCCCGGGAAGCACAAGGCTCCACGTTCCTTGGCCAGGGCATCGCCATTCCCCATGGTACGCCGCA\n",
        ">d31733.assembled      GCCCGGGAAGCACAAGGCTCCACGTTCCTTGGCCAGGGCATCGCCATTCCCCATGGTACGCCGCA\n",
        ">d39114.assembled      GCCCGGGAAGCACAAGGCTCCACGTTCCTTGGCCAGGGCATCNCCATTCCCCATGGTACGCCGCA\n",
        ">d39404.assembled      GNCCGGGAAGCACAAGGCTCCACGTTCCTTGGCCAGGGCATCGCCATTCCCCATGGTACGCCGCA\n",
        "//                                                                                      \n",
        ">d30181.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGGCAGCAGGGGCAGCAGCTTGACCTTGCAGGTTCACGACG\n",
        ">d30695.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGGCAGCAGGGGCAGCAGCTTGACCTTGCAGGTTCACGACG\n",
        ">d31733.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGGCAGCAGGGGCAGCAGCTTGACCTTGCAGGTTCACGACG\n",
        ">d39114.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGGCAGNAGGGGCAGNAGCTTGACCTTNCAGGTTCACGACG\n",
        ">d39253.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGGCAGCAGGGGCAGCAGCTTGACCTTGCAGGTTCACGACG\n",
        ">d39404.assembled      CCAGGGCAGCGAAGGCCTGAGCGTCGNCAGCAGGGGCAGCAGCTTGACCTTGCAGGNTCACGACG\n",
        "//                                                                                      \n",
        ">d35320.assembled      GCCCGAGTCGGCCCCGACCAGCCAGCCCAGGACGTCGGAGCGGTAGGTGGCGAAGTCGTAGAGCT\n",
        ">d39114.assembled      ACCGGAGTCCGCGCCGACCAGCCAGCCCAGGACGTCGGAGCGGTAGGTGGCGAAGTCGTAGAGCT\n",
        ">d39531.assembled      RCCGGAGNCCGCRCCGACCAGCCAGCCCAGGACGTCGGAGCGGTAGGTGGCGAAGTCGTAGAGCT\n",
        ">decor21.assembled     GCCGGAGNCCGCACCGACCAGCCAGCCCAGGACGTCGGAGCGGTAGGTGGCGAAGTCGTAGAGCT\n",
        "//                     *  -     -  *                                                    \n",
        ">d19long1.assembled    ATATTTTTGGAGGATGAAGAGTTTCTGTCTTCTGCATTCTCTTCTTCTTCTTCTTCCTCTTAGGT\n",
        ">d35422.assembled      ATATTTTTGGAGGATGACGAGTTTCTGTCTTCTGCATTC---TCTTCTTCTTCTTCCTCTTAGGT\n",
        ">d39104.assembled      ATATTTTTGGAGGATGACGAGTTTCTGTCTTCTGCATTC---TCTTCTTCTTCCTCCTCTTAGGT\n",
        ">d41058.assembled      ATATTTTTGGAGGATGAAGATTTTCTGTCTTCTGCATTC---TCTTCTTCTTCCTCCTCTTAGGT\n",
        "//                                      *  -                                *           \n",
        ">d35320.assembled      CGCCGGGCGTTCGAACCGAGCGCGGGCGTCGCGCCGCACGCGGCGCGGCGCGNC\n",
        ">d39103.assembled      CGCCGGGCGTTCGAACCGAGCGCGGGCGTCGCGCCGCNCGCGGCGCGGCGCGCC\n",
        ">d39404.assembled      CGCCGGGCGTTCGARCCGAGCGCGGGCRTCGCGCCGCACGCGGCGCGSCGNGCC\n",
        ">d39531.assembled      CGCCGSGCGTTCGARCCGAGCGCGGGCRTCGCGCCGCACGCGGCGCGGCGCGCC\n",
        ">d39968.assembled      CGCCGGGCGTTCGAACCGAGCGCGGGCGTCGCGCCNCACGCGNCGCGGCGCGCC\n",
        "//                          -        *            *                   -      |\n",
        ">d30695.assembled      GGCGTACGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">d31733.assembled      GGCGTAYGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">d34041.assembled      GGCGTACGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">d39114.assembled      NGCGTACGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">d39253.assembled      GGCGTACGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">d41732.assembled      GGCGTAYNACTCGTCGTTCAGCCAGGCCAGCGT\n",
        ">decor21.assembled     GGCGTACGACTCGTCGTTCAGCCAGGCCAGCGT\n",
        "//                           *                          |\n",
        ">d30695.assembled      ACGCGAGGCCCATTCCAAA-TTCGAGAGGACTTACACAGATATCGCCGCAGATAACTTTTACGAT\n",
        ">d39253.assembled      ACGCGAGGCCCATTCCAAAGTTYGAGAGGACTTACACAGATATCGCCGCAGATAACTTTTACGAT\n",
        ">d39968.assembled      ACGCGAGGCCCATTCCAAAGTTCGAGAGGACTTACACAGATATCGCCGCAGATAACTTTTACGAT\n",
        ">d40328.assembled      ACGCGAGGCCCATTCCAAAGTTCGAGAGGACTTACACAGATATCGCCGCAGAYAAYTNTTACGAT\n",
        "//                                           -                             -  -         \n",
        ">d30181.assembled      TCGGTGCTTCCGGAACTCCGCAACGNG\n",
        ">d30695.assembled      TCGGTGCTTCCGGAACTCCGCAACGCG\n",
        ">d39103.assembled      TCGGTGCTTCCGGAACTGCGCAACGCG\n",
        ">d39968.assembled      TCGGTGCTTCCGGAACTCCGCAACGCG\n",
        "//                                      -         |\n",
        ">d33291.assembled      GTTGCGCGCGATCTCCACCATCTGCTGCTGGCCGATGCCCAGCGTGCCGACCCGCGTGCCCGGGT\n",
        ">d35178.assembled      GTTGCGCGCGATCTCCACCATCTGCTGCTGGCCGATGCCSAGCGTGCCGACCMGCGTGCCCGGGT\n",
        ">d35320.assembled      GTTGCGSGCGATCTCCACCATCTGCTGCTGGCCGATGCCSAGCGTGCCGACCCGCGTGCCCGGGT\n",
        ">d39253.assembled      GTTGCGCGCGATCTCCACCATCTGCTGCTGGCCGATGCCCAGCGTGCCGACCNGCGTGCCCGGGT\n",
        ">d40328.assembled      GTTGCGCGCGATCTCSACCATCTGCTGCTGSCCGATGCCGAGCGTGCCGACCMGCGTGCCSGGGT\n",
        "//                           -        -              -        *            *       -    \n",
        ">d30695.assembled      CGACGACACGTCGGAGACCGGCTCCTCGACGAGAACCACTCCGGCAGCGCGACACGCACTCAGTA\n",
        ">d39531.assembled      CGACGACACGTCGGAGACCRGCTCCTCGACGAGAACCACTCCGGCAGCGCGACACGCANTCAGTA\n",
        ">d41389.assembled      CGACGACACGTCGGAGACCGGCTCCTCGACGAGAACCACTCCGGCAGCGCGACACGCACTCAGTA\n",
        ">decor21.assembled     CGACGACACGTCGGAGACCRNCTCCTCGACGAGAACCACTCCGGCAGCGCGACACGCANTCAGTA\n",
        "//                                        *                                             \n",
        ">d31733.assembled      CGGCAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d35178.assembled      CNGCAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d35320.assembled      CNGCAGCTCCAGCTGGGGCCGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d35422.assembled      CGGCAGCTCNAGCTGGGGCCGCGCCAGCCGCTGCCGCNGNTC\n",
        ">d39103.assembled      CGGCAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d39187.assembled      CGGCAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d39253.assembled      CGGCAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d39404.assembled      CGGYNGCNCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        ">d39968.assembled      CGGCAGCTCCAGCCGGGGCTGCNCCAGCCGCTGNCGCAGCTC\n",
        ">decor21.assembled     CGGNAGCTCCAGCCGGGGCTGCGCCAGCCGCTGCCGCAGCTC\n",
        "//                        -         *     *                      |\n",
        ">d30181.assembled      AGCCTGCTSGGCTACGTNGTCCGCTGGGTCGACGCCGGGGTGGG\n",
        ">d39187.assembled      AGCCTGCTSGGCTACGTCGTCCGCTGGGTCGACGCCGGNGTGGG\n",
        ">d39253.assembled      AGCCTGCTGGGCTACGTCGTCCGCTGGGTCGACGCCGGNGTGGG\n",
        ">d39968.assembled      AGCCTGCTGGGCTACGTCGTCCGCTGGGTCGACGCCGGCGTGGG\n",
        "//                             *                             -     |\n",
        ">d30181.assembled      CAGCGAGGTGTCCACGCGGCTGCCCTGGCCGGTGCGCAGCCGCTTCACGTAGGCGGTGACCACGC\n",
        ">d33291.assembled      CAGCGAGGTGTCCACGCGGCTGCCCTGGCCGGTGCGCAGCCGCTTCACGTAGGCGGTGACCACGC\n",
        ">d39103.assembled      CAGCGAGGTGTCCACGCGGCTGCCCTGGCCGGTGCGCAGCCGCTTCACGTAGGCGGTGACCACGC\n",
        ">d39187.assembled      CAGCGAGGTGTCCACGCGGCTGCCCTGGCCGGTGCGCAGCCGCTTCACGTAGGCGGTGACCACGC\n"
       ]
      }
     ],
     "prompt_number": 81
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Infer a tree"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "raxmlHPC-PTHREADS-AVX -s analysis_pyrad/outfiles/c88d6m4p3.phy \\\n",
      "                      -m GTRGAMMA \\\n",
      "                      -f a \\\n",
      "                      -x 12345 \\\n",
      "                      -p 12345 \\\n",
      "                      -N 100 \\\n",
      "                      -w /home/deren/Documents/PEGBS/analysis_raxml/ \\\n",
      "                      -n c88d6m4p3 \\\n",
      "                      -o \"d35320.assembled\" \\\n",
      "                      -T 20 &> /dev/null"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext rmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "library(ape)\n",
      "tre <- read.tree(\"analysis_raxml/RAxML_bipartitions.c88d6m4p3\")\n",
      "plot(tre)\n",
      "nodelabels(tre$node.label)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nOzdd1xT5/4H8AfCJkCC\nhiF7iFtZAiJaoVVx1YKgFRTEnxUR6y1S1KoIXkctakVUZKogoChYx1VrRcGBIho2ogwZCmGFoYwI\nJOf3x2lzueQEUUYifN8v/0i+efLkSdJ+PZ6c8zkiGIYhAAAAwkdU0AsAAABADBo0AAAIKWjQAAAg\npKBBAwCAkIIGDQAAQgoaNAAACClo0AAAIKSgQQMAgJCCBg0AAEIKGjQAAAgpaNAAACCkoEEDAICQ\nggYNAABCCho0AAAIKWjQAAAgpKBBAwCAkIIGDQAAQgoaNAAACClo0AAAIKSgQQMAgJCCBg0AAEIK\nGjQAAAgpaNAAACCkoEEDAICQggYNAABCCho0AAAIKWjQAAAgpKBBAwCAkIIGDQAAQgoaNAAACClo\n0AAAIKSgQQMAgJCCBg0AAEIKGjQAAAgpaNAAACCkoEEDAICQggYNAABCCho0AAAIKWjQAAAgpKBB\nAwCAkIIGDQAAQgoaNAAACClo0AAAIKSgQQMAgJCCBg0AAEIKGjQAAAgpaNAAACCkoEEDAICQggYN\nAABCCho0AAAIKWjQAAAgpKBBAwCAkIIGDQAAQgoaNAAACClo0AAAIKSgQQMAgJCCBg0AAEIKGjQA\nAAgpaNAAACCkoEEDAICQggYNAABCCho0AAAIKWjQAAAgpKBBAwCAkIIGDQAAQgoaNAAACClo0EC4\nXL58+dtvv0UIlZeXL1++nEajUSgUT09PDoeDEKqvrxcTExP5x5EjR/iNbGtr+/HHH1VVVSkUiru7\nO5vNxuc/evSovr6+srLy6dOnBfcukYeHR2hoaI/isWPHNm7c2Jencz8lMLxBgwbC5fHjx9OnT8cw\nzMnJyczMrLCw8NGjRxEREW/fvkUIJSYmrly5svQfGzZs4DfS29u7tLT06dOnycnJiYmJN27cQAgF\nBwcnJCTcunULb4Wtra2Cepvp6ekzZ87sUXz27Nn06dP78nT8UxqEdQEhgwEgaElJSebm5lQq1c/P\nb8aMGbdu3Wppably5QqHw2GxWLGxscbGxhwOB8Mwa2vru3fvdn8u4cimpiYymVxdXY2PWbBgwenT\npzs7O7W0tEpKSjAM43A4FAqluLi4+1SPHj0yMzOTkZExMjK6dOkSXjx+/LiWlpaMjMzq1as7OzsJ\nK4WFhQsXLlRQUDA2Ns7OzsbnV1VVPX78uK6uroaGxv37911dXUeNGjVnzhw2m93e3j5q1Kjdu3dr\naWnp6uo+ePAAfy0DA4Pc3FzCCQk/pcH6PoDQgAYNBOzJkyfKysr37t1jMpmzZ88WFRWtq6vDHwoJ\nCUEIUSiU2tpaDMMYDIa4uLi1tbWMjMz06dNfvnzJnaTHSBaL9fbtW/yh/Px8KpX69u3bq1evzpw5\nEy92dnZKSEiwWCzuDO3t7TQa7erVq83NzTt27MBH3rlzx8DAoKys7PXr1/Ly8hkZGbwVBoOhpKR0\n+vTp5ubmLVu2mJubYxhWWVmJEDp06BCLxbK3t5eTk6PT6XV1dWQyOTc3Ny0tDSG0b9+++vr6ffv2\n6enpYRjW2NhIJpO7uroIJ+zlUwLDmJhgt98H1bZt2xQVFQW9CsCXiIjI5s2b/fz8/P39ra2tEUIr\nV6588+bN6NGj8QHu7u62trbu7u5BQUF79+5NT0+3tbU9cOAAlUr19PTcuXNnQkIC4UhJSUk1NbWu\nrq6oqCg/P7/o6Gg1NbWgoKD58+fj40tLS1VUVCQlJbmLERUVffbsmZaWVmtra0dHx7hx4xBCubm5\no0ePFhER0dHRKS4uplAoJ06c6FHx8/P7/vvv3dzcEELr1q2Ljo5GCGVlZU2cOHHLli2ioqJKSkr/\n+te/jI2NMQyTkJCg0WjJycl2dnY7d+5ECK1fv37Xrl0dHR10Ot3IyIhEIgUFBfFO2MunBIax4dyg\n7969e/DgQUGvAvAlIiIiJiZ2//79kydP4pXW1tbp06d/+PAhIyNjxowZCCEtLS1bW9tXr14hhL79\n9lvuL2Ourq7BwcH8RiKEUlNTPTw8pk2blpaWpq6ujhAqKSmxt7fHH719+3aPXcAlJSW+vr4VFRXq\n6uoFBQU//fQTQsjFxaWoqGjq1Km6urr79+9fsGABb+XPP/8sLCyMiopCCGEYNnbsWIRQZmamra2t\nqKgofvvo0aMIoTdv3khLSysrKz979szR0RF/3fb2dj09PQkJCe4OaN4JORwO76c0KF8JEDYC3oIf\nTHPmzBH0EsBHVFdXi4qKtrW1YRjG4XDMzc0PHTpUWloqJyfHZrPxMU5OTuHh4Uwm8/jx49zi7t27\nfX19CUdiGBYTEzNhwoTHjx93fy1TU9MbN25gGMZms83Nzbl7fjEMY7FY6urqd+7cYbPZHA5n6tSp\ndDq9ubn5yZMnGIZ1dnaePHly/PjxvBUMw+Tl5UtLS/F5SktL6+vrMQxzcHCIjY3FR5LJ5NbWVgzD\n/vjjjyVLlmAYNn78+Js3b+JPCQgIWL9+PYZh9vb258+fJ5yQ8FMauC8BCC9o0ECQOjo6yGRySEhI\nXV2dt7e3mJhYSkoKi8WSkpI6efIkk8kMDQ3V19dvb29vbGwkkUhhYWGNjY3nz58fN24ck8kkHMlg\nMGg0WmpqKvdgD7y1LVu2zNnZubKy8pdffpk9ezb+q2NycnJWVlZubq6cnFx9fX1FRQW+56G9vf3x\n48dSUlKZmZmNjY0eHh729va8FQzDrKysvLy8ampqEhISFBUV8d6qp6dXUFCAYVhOTs6kSZPwN+vr\n6+vv79/U1CQtLe3s7NzU1BQfH6+kpMRgMDAM09DQwH+05J2Q8FMSzBcGhhY0aCBgMTExNBptypQp\noaGhZDL53bt3GIadO3duzJgxVCp16dKllZWV+MgzZ85oampSKJTFixdXVFTgRd6RgYGBPf6Z+OLF\nCwzDCgoKjIyMKBSKo6NjU1MT/nR9ff2wsLCOjo7FixeTyWRjY+PLly+PHTt2z549HR0dbm5uZDKZ\nSqUuX768urqat4JhWF5enqGhoZSUlIWFRWpqKoZhzc3N3O36s2fPuri44K+1aNGia9eu3b17d/Xq\n1c7OzrKyspaWlunp6RiGVVdXU6lU/O8M3gn5fUpg2BPBMEwQe1aGgrW1dXJysqBXAQAAnwlOVAEA\nACEFDRoAAIQUNGgAABBS0KABAEBIQYMGX4beU+5aW1vd3NyUlJQ0NDRiY2O7P7GsrMzQ0LB7pZdA\nu6NHj9rZ2Q3yW+kNpNyB7qBBgy9DLyl3GIYtW7ZMRUUlPz9/69at69ate//+PUKIwWAcOXJk4cKF\nBgYG3Hl6CbS7fv26t7e3kZGRAN7ePyDlDnQHDRoIr7t371pYWCgqKvr7++Otp62tbevWrd7e3jIy\nMjk5OZMnT9bQ0EhLS6urqztw4ACNRvP09JSQkKipqUEIFRQUlJaWstlsExMTfMKurq6AgIBz586N\nHTt2xYoV0tLS1dXV+ENZWVl+fn7jx4/nbdCpqanm5uaysrLGxsbc9I8TJ05oa2vLysq6uLh0dXUR\nVoqKihYtWkShUExMTHJychBCGIaNGTPmxIkTenp6mpqaDx48WLNmzejRo62trfFAvvLy8kuXLmlr\na+vp6T18+BB/LW6D5p2Q8FMazO8EDC3BHoY9qOBElS9a31PuAgMDXV1d8YcqKirExcU7Ojq485ib\nmyclJeG3+QXaVVZWTp06taCgQFZW9s2bN92XASl3QICGc1gSGHClpaUNDQ2ED4WGhnZ0dPQoiouL\njxo16lNfRVRU1Nvbu+8pdzQa7c6dO3l5eZKSkmvXrtXX1xcXF8dHdnZ2ZmdnGxsb43dTU1N5A+3a\n2tocHBwCAwNFRERkZGTU1NR6LAZS7oCgQIMGn8Bm4eKxU437Pj7rzt3oyHAxsU/+z0xaWrrvKXcO\nDg7Xrl2ztLQ0MzPT0dHR0NDgzpOfn6+mpkalUvG7vIF2GIatXr3azc3N2tr6woULRkZGIiIi3VcC\nKXdAkAS7AT+oYBfHgJtkZpn4sqrvfwyt5nQPxe+7vqfcYRhWVFTU1dWFYRibzZ4+fTo3KA7DsIiI\niOXLl3Pv8gbalZSUIISUlZWVlZXJZLKUlNSCBQu44yHlDggW/EgIPtPTO7d+XDDL1WxC0LbNH9rb\nuj9U+brYyUi/P5MrKirKyMhER0fX19f7+PjQ6fTp06erqqp2dnaGhIQ0NDSEhYWlp6evWrUKIWRq\nahoVFcVkMn18fJSVlefNm8ed5/nz56ampty7WlpacXFxVVVVu3btkpSUtLKyUldXZzAYWVlZ+M6H\n0NBQ/Ci9lJSU7OzsoqKi5uZmIyOjyspKX1/fnJyciRMn5ufnW1tbZ2VltbS05OXlEVYQQlOnTg0K\nCqqtrU1MTDQxMcEPLMnMzMT3txQUFODXzUIIZWRkmJiYNDc3l5eXx8bGNjc3X7x48fDhw3v27EHd\nfiHknZDwU+rPxw6EjqD/hhhEsAU94Lhb0EE3H8jKy/87OuHMkzzzuQsWuazjbjVfyn8zztBElETq\nzxY09ikpd2fPnlVWVlZUVPTw8OgR82Zqatr9GoaEgXa4jo4OaWnpmpoa/C6k3AFhAGl24BNMNp/5\n7+gEhNDV06feFBVu+vUoQqiqtGTHym/PpuXjY/4TFf4q83nanZuX8t/sXeeUlvRn9ytLAQD6DnZx\ngM/R1dFBEiPht8XExd83NeJ7Oaorym5fiF7lvUOgqwNgmIAGDT7HVMvZabdvFOdmvW9qPHfkAELo\nQ3s7xuGc8vVZs81Pmiwn6AUCMBxAgwafY+xUI6ctv/y2ae1Pi+Zoj58oSiKRFShJCXGKSsomc74R\n9OoAGCagQYPP0d7aMt16Xvj9jMjU7MnmMzX0DERJpNy0R8/u/bXGYtKP82dy2Ow1FpNa378bpAX0\nSAXqkYhEmJ3U98ihvuDNYBp6kKw07EGDBp+DWc3YZGtVQE+vZ1TGHNk/f6ULQuiH3b8eu3H/yJWk\nf5+7LCIqeuRKkowseZAWwA2d4E1EwvhkJ/U9cqh3hBlMAgHJSsMeNGjwN39//+Ufwx2srjd29c87\nA378v1++XzLFfOa8710QQnIU6igV1VEqqlQlZRERkVEqqiKiA/kfGGEqEG8iEr/sJG7nysnJmTdv\nnoKCgrGxcX5+PkIIwzAVFZXQ0NBJkybJy8vjWR8IoadPn1pZWdFoND8/PzMzs7KyMsJX7A6SlcBA\nEuxRfoMKjoP+JLNnz45KL+j9z5CdScir91Sg7olIhNlJ3MihoqIiGo0WGxvb2Ni4cePGuXPnYhhW\nXl6OENq+fTuLxQoKCjIwMMAwrKCgQFFR8Y8//igvL7e0tKTRaPjxyLyvyAXJSmBgQRYH+JuoqChZ\nXqH3MS1NjSG7t/Z9zsrXxT2iLT5bL6lAPRKRCLOTuJFDgYGBmzdvdnJyQgi5urquXLkS/ZNhtH//\nflFR0SlTpuAzHzhwwN/f/7vvvkMIOTo6Xrt2jfteerwiFyQrgYEFDRrwtX/9qjXb/dV0/z5puyQv\nm81hP7pxZZLZDLv1P4pLSCCE3pYUJp461lBbgxcvBP42WVdTQkICf4rWQtvdu3f3fyWGhoa9pAL1\nSEQizE7i/sP/1q1bf/75Jz6ysbFRT08PIZSZmbl48WJuhhG+7+LWrVt4mBFCqKGhoXtOdI9X5IJk\nJTCwoEEDArlPHj26eSXjwb012/3xCpvddWjzDw4ePxnPtgnavjk//fEy981sdtehH9d1L8rKK/zw\nww/y8vIDux4Oh9PZ2YkHgWIYdunSJQcHB+6jdDq9++7gioqK2NhYEonE4XAsLCy4iRaOjo5dXV2l\npaUUCgUfefHiRVtbW4RQVlYWvk2NEMrIyPjmm28aGxubmprwN4JhGH6xFX6viPvw4cO8efPOnDlj\nY2MjIiJiaGg4ffr0d+/eFRUVBQcHBwUFhYWFbdmyZebMmT0qCxYsKCkpycvL09bWRgiVlZXJycnh\nq8Kvv9XV1ZWfnz9t2jR8efiW+7Nnz/DNf4RQfHz8119/zX2bCCHeCevq6nr5DIGQEuT+lUEG+6A/\nyZw5c7j7jl22+i5wdpOUkg66+QCv+J+5qKarj9/+d3SCqpYOYdFi7kImkznga+vo6CCTySEhIXV1\ndd7e3mJiYikpKdxHN2zYEBAQwL2roKAQGRlZX1+/ZcuWxYsX40F3GhoaxcXFGIZpaWkdPHiQyWQG\nBgZqamo2NzfjRfxRDMMmTpyYm5vb1dUlLy+fkJDw7t27n3/+GSGUm5vL7xWTk5OzsrJyc3Pl5OTq\n6+srKirwPQ/t7e2PHz+WkpLKzMxsbGz08PCwt7fnrWAYZmVl5eXlVVNTk5CQoKioiKfW6enpFRQU\nYBiWk5MzadIk/LV8fX39/f2bmpqkpaWdnZ2bmpri4+OVlJQYDEb3t8k7Ye+fIRBOcBQHILB0rcc6\n3/2yCv/dJV1b+UbLYAJ+W9NgQl3VWwzDeIuDtB5xcfGQkBBfX18bGxsDAwMpKanu+3+fP3/efXv2\n2LFjO3bsMDAwaG9vj4uLI5FINTU1LS0turq6CKGoqKjTp0/r6+vfu3cvNTVVXl4e31jW0dFBCLW2\ntlZUVIwfP55EIkVERHh6ehoaGlKpVCqVOmHCBH6v+MMPP6Snp48bN+6rr77S1tb+7rvvTExMxo4d\nGxAQYGpqunLlylmzZunq6jKZzODgYN4KQigkJCQ5OVlLS+vw4cPXr1/X1tZ+9+5dbW0tfiQfHneH\nvxa+BU2n0/HtXzU1tWPHjv3nP/9RUVHp/jZ5J+z9MwTCCcKSRqIFCxa0tLRISUl1L9Lp9Ign+d0r\nP3xljG8jI4T+CD9ZVVbiuf93hBC7q3P5ZK2Y56/+PB/dozjdet6fCecVFRWH8N0MChaLJSoqiu9M\nj4qKunr16uXLlwW9KDDiwD7okUhLS8vLyws/xoAL/3GfH7KCAqvt79DntpYWEklMSkaWt0gikQZp\nzUNs27ZtjY2NYWFhb9682bt3b0xMjKBXBEYi2MUB+kRZXetN0Sv89tuSIpqauoioKG8RDdBBdQLn\n6emZn5+vpKS0bNmygIAACwsLQa8IjESwBQ36ZJL5jNZ3TU9u/8dolvXVyOCvvl1GWCx/VSDolQ4M\nAwMDOp0u6FWAkQ62oEGfkEhi24PPJpwKdLeeTlag2q3/kV8RADBQYAt6WGlsbIyIiOBwOL0Py8rK\nCg8Pp9Fo3Yvl5eUt75q7V45eT0YIcYvKmtp7ohPx2x9Y7R9Y7bzFrs7OgXgfAACEoEEPM5mZmVeS\nH82Yv6j3YVYr1yKE6v63qKQ/fv/6VQghdmdn5/smc3Nz3ic2NzeXlpb2ErM5SWsMmTwwCXaXL18+\ne/bstWvXuJWjR48+ePDgjz/+4FbKysq+++67rKys7mNOnjz5/v37X3/9de3atQih+vp6FRUVNpuN\nDzh8+LC3t3djY6OPj8+1a9c+fPiwcuXKU6dODdQp6Z/Kw8PD0NDQ3d29e/HYsWOvXr3Cj8DrHe+n\nBIYTaNDDjaa+gaXtks94IvdZjXU1d4MD4uLieMfk5ORERkYeO3asX0vsmx5Za/i5fP7+/vhdBoMR\nFxcXGRk5efJk7pjg4OCEhIRbt27R6fQ1a9asWLFCVlY2MTFx5cqVe/fuxcfg/25Ys2YNvpe5tbXV\nzMxs8+bN+KW4h156erqnp2eP4rNnz/AzAz8KEumGN9gHPcztX7+q8nUx925JXvbP9vNczMaf+MWr\ng8XqpSgQ/MIws7Ky/Pz8xo8fzw3E4M387OrqCggIOHfu3NixY1esWCEtLV1dXY0Qio+Pd3Nz0/6H\nrKwsQigtLW3WrFmqqqolJSU0Gg0/JZoLIkOBsBDwmYyDaQSe6n337t3vf/yZe3L2N45OCCHu6doX\n8ytoY9Q99h4Ov58xZYaVk9d2wmLEw8yVK1cSzp+dnb158+ZBWjy/MMzKysqpU6cWFBTIysq+efOm\n+1O6Z35evXoVj/fEMKyzs1NCQoLFYjEYDHFxcWtraxkZmenTp798+RIfgJ+9LSIiIi8vX1lZ2X1O\niAwFwgN2cQxbrwtyxSUkJaWkuZUX6WkSUlJ413b0+OmUr88y9828RRv7FfzmZLFYT5482b59+4Cv\nVl9f/9KlS7xhmG1tbQ4ODoGBgSIiIjIyMnjWD65H5mdqaur8+fPx26WlpSoqKpKSkunp6ba2tgcO\nHKBSqZ6enjt37kxISAgPD3/48GFxcfHo0aMXLlwYFhbG3XOCIDIUCBNo0MPW0rUeCKGnSbe4lT7m\naWD8z/6fOnXq/v37B+NMbllZ2U2bNvUIw8QwbPXq1W5ubtbW1hcuXDAyMur+U16PzM+SkhJ7e3v8\n9u3bt/FrQX377bfcy+65uroGBwdzOJxdu3alpKTgQaNOTk4pKSndVwKRoUB4DOcGXVhYaGpqKuhV\nDKn3798bz+d7GdB3jQ1SsrL4bRkyuauzk9XaQlD85+xtXlJSUnPnzh3YNeNqamp4wzBLS0svX76c\nmprq6+vb2tra1dW1cOHCmzdv4k/pkflZXl6O54hyOJyYmJhDhw41NDTExcVt3LgR749ZWVkzZsxo\naGhgMpl4ohBC6NWrV90PSoHIUCBUhnODNjAwGGlhSffu3Qu/covfo33M05CSluY3w+BRVFSUkZGJ\njo5etmzZwYMH6XT6b7/9pq6uzmAw8AFLly719PRcsuS/B6g8f/68+1/AWlpacXFxhoaGJ06ckJSU\ntLKyam5u/umnnyQlJR0dHf/888/4+PjHjx9TKBQqlXro0KGNGzf++eef169fz8jIQAilpKRQqVQS\nidTc3GxkZFRZWRkaGpqTkzNx4sTMzEwbG5snT55oa2vn5eVNnDgxPz+/RwUhNHXq1KCgoO3btz98\n+HD9+vV0On3UqFGZmZn4ASQFBQVaWloyMjLon3S65ubm8vLy2NhYS0vL27dvHz58ODs7GyH07Nmz\ngIAAwgnV1NR4P6Wh+5LAkIOjOEaQPuZpDOyVXvuIMAxTQkJCRUVFRUVl1KhRubm5tra23S9i0iPz\nc9++fS9evJg0aVJxcTF+eSoKhRIREbFv3z4dHZ3Y2Ng7d+4oKiqKiopeuHAhLi5OW1s7Ojo6OTkZ\nnxMiQ4EQgrjRYQXfgnb03MKtdI8MZbO7NlhPX7tzr9Es68CfPXUnTlm+yZu3OHfFKn7HQQMAhhJs\nQY8gkKcBwJdlOO+DBgih8PsZ3e/qTZ525EpSjzGERQCAwEGDHlbk5eWTr1x6du+v/kzC4bCXLVow\nUEsCAHw2aNDDiqmpaXVF2UeHbdiwgfeKKgAAYQP7oIFwuXz5Mn5qSXl5+fLly2k0GoVC8fT0xDNU\n6+vrxcTERP5x5MiR7s89evQoftAxQqitre3HH39UVVWlUCju7u54mh1hUSA8PDxCQ0N7FI8dO7Zx\n48a+PJ37KYHhDRo0EC54+g+GYU5OTmZmZoWFhY8ePYqIiHj79i1CCI+mK/3Hhg0buE/E4+64aUre\n3t6lpaVPnz5NTk5OTEy8ceMGv6JApKen4+c6dsc9h/CjICNphIAGDQSPN56tra1t69at3t7eMjIy\nOTk5kydP1tDQQHyi6RBP3F1zc3NMTExkZKSmpqaRkZGZmRmTySQsdl8GhNgBoSOgkKahMALT7PrI\n3d2dm+smcL3Es4WEhCCEKBRKbW0thmH8oul44+5YLNbbt2/xR/Pz86lU6tu3bwmL3GVAiB0QQvAj\n4UhUVVV15MiRwcg8+lRff/314cOH+cWzubu729rauru7BwUF7d27lzCajjDuTlJSUk1NraurKyoq\nys/PLzo6Gq8TFnEQYgeEEDTokejw4cPv3r0T1EWeulNRUeGNZ/vw4UNGRsaMGTMQQlpaWra2tq9e\nvUJE0XQY/7i71NRUDw+PadOmpaWlqaur91LEQYgdEEaC3YAfVLCLQ/hVV1eLioq2tbVhGMbhcMzN\nzQ8dOlRaWionJ8dms/ExTk5O4eHhTCbz+PHj3OLu3bt9fX1LSkoQQsrKysrKymQyWUpKasGCBRiG\nxcTETJgw4fHjx91fi7CIY7FY6urqd+7cYbPZHA5n6tSpdDq9ubn5yZMnGIZ1dnaePHly/PjxvBUM\nw+Tl5UtLS/F5SktL6+vrMQxzcHCIjY3FR5LJ5NbWVgzD/vjjjyVLlmAYNn78+Js3b+JPCQgIWL9+\nPYZh9vb258+fJ5yQ8FMamC8ACDdo0ECQOjo6yGRySEhIXV2dt7e3mJhYSkoKi8WSkpI6efIkk8kM\nDQ3V19dvb29vbGwkkUhhYWGNjY3nz58fN24ck8n88OED4x9mZmZRUVENDQ0MBoNGo6WmpnIP9mhr\nayMsYhiWnJyclZWVm5srJydXX19fUVGB73lob29//PixlJRUZmZmY2Ojh4eHvb09bwXDMCsrKy8v\nr5qamoSEBEVFRby36unpFRQUYBiWk5MzadIk/M36+vr6+/s3NTVJS0s7Ozs3NTXFx8crKSkxGAwM\nwzQ0NIqLiwknJPyUBPSNgSEFDRoIWExMDI1GmzJlSmhoKJlMfvfuHYZh586dGzNmDJVKXbp0KfeS\nVGfOnNHU1KRQKIsXL66oqOg+SUdHh7S0dE1NDYZhgYGBPf6Z+OLFC8IihmH6+vphYWEdHR2LFy8m\nk8nGxsaXL18eO3bsnj17Ojo63NzcyGQylUpdvnx5dXU1bwXDsLy8PENDQykpKQsLi9TUVAzDmpub\nuf8COHv2rIuLC77IRYsWXbt27e7du6tXr3Z2dpaVlbW0tExPT8cwrLq6mkqlcjgcwgn5fUpg2IM0\nOwAAEFJwHDQAAAgpaNAAACCkoEEDAICQggYNAABCCho0AP3Se/xea2urm5ubkpKShoZGbGws/hTC\nUD3CoD7COQUC4vcEAho0AP3SS/wehmHLli1TUVHJz8/funXrunXr3r9/j/iE6vEG9RHOKai3CfF7\ngiHYo/wGFRwHDQZJUlKSubk5lUr18/ObMW1hL8AAACAASURBVGPGrVu3Wlparly5gifVxcbGGhsb\nczicx48f4zcwDGOz2fLy8kVFRU1NTWQyGT+GGsOwBQsWnD59GsMwa2vru3fvdn8Vwjm7D3j06JGZ\nmZmMjIyRkdGlS5fw4vHjx7W0tGRkZFavXt3Z2UlYKSwsXLhwoYKCgrGxcXZ2NoZhHA5HVVX1+PHj\nurq6Ghoa9+/fd3V1HTVq1Jw5c9hsdnt7+6hRo3bv3q2lpaWrq/vgwQP8tQwMDHJzcwknJPyUBuv7\nGL6gQQPwafoevxcYGOjq6oo/VFFRIS4u3tHRQRiqxy+oj3dOLojfGwkgLAl8RFNTEx550buAgABp\naeleBtBoNBKJNHDrEgBRUVEvL69eguV6xO/RaLQ7d+7k5eVJSkquXbtWX19fXFwcEYXqXbt2jTeo\nj3DO7ouB+L1hDxo0+Ihff/01vbCMSqN9ZJwCrZP/g+WvCuzHjrW1tR3QpQmAvLx83+P3HBwcrl27\nZmlpaWZmpqOjg19zABGF6vEG9fGbkwvi90YCaNDgIzgczsJVbnqTp/VnkrsJ5zU1aSYmJgO1KkGp\nqanp7OzEg6QxDLt06ZKDgwODwZg/f35TUxPe4J49e4ZvOVZUVMTGxpJIJA6HY2FhsWfPHoRQbGzs\n/v37IyMj8eaLEGpoaIiLi9u4cSP+9KysrBkzZvCbE/fhw4d58+adOXPGxsZGRETE0NBw+vTp7969\nKyoqCg4ODgoKCgsL27Jly8yZM3tUFixYUFJSkpeXp62tjRAqKyuTk5PDXxS/nGNXV1d+fv60adMQ\nQhkZGcbGxvirr1y5En/p+Pj4r7/+Gi/iXZt3wrq6Ot5PafC/nGEIjuIAn2z/+lWVr4u5d0vysn+2\nn+diNv7EL14dLFYvxWFAUVFRRkYmOjq6vr7ex8eHTqdPnz5dVVW1s7MzJCSkoaEhLCwsPT191apV\nCCFTU9OoqCgmk+nj46OsrDxv3rzq6movL6+IiAhVVdWysrKysrL29nZRUdGffvopMjKyqanpwoUL\n8fHxP/30E785U1JSsrOzi4qKmpubjYyMKisrfX19c3JyJk6cmJ+fb21tnZWV1dLSkpeXR1hBCE2d\nOjUoKKi2tjYxMdHExAQ/sCQzMxPvxQUFBfgvigihjIwMExOT5ubm8vLy2NjY5ubmixcvHj58GP+b\nhrsFzTsh4ackuC/tSybQPeCDC34kHBA///xzQMKtxJdViS+r/M9c/MbRCSEUdPMBXrmYX0Ebo+6x\n93D4/YwpM6ycvLYTFjfuOxIdHS3otzIw+h6/d/bsWWVlZUVFRQ8PD3wYv1A9wqA+wjkhfm9EgTQ7\n8BE+Pj5KFl/juziunj5VV/n2XuKFQ5dvq+nqI4RynzwK37sj6OYDhFB++uNTvj4nbqfyFu1+2GSh\nSVu9erVg3wsAXxbYBw0+wdK1Hgihp0m3uJXayjdaBhPw25oGE+qq3mIYxltEw3c7AIDBAw0a/Nfm\nzZvxnY/dPXjwwMHia35PedfYICUri9+WIZO7OjtZrS28xY6OD4O0ZgCGMWjQ4L/S09O5h0ZxVVVV\n9fIUsoICq60Nv93W0kIiiUnJyPIWxSUkB2PBAAxv0KDBf0lLS/MeCaesrNzLU5TVtd4UReK335YU\n0dTURURFCYpCcAVxAL44cJgd6JdJ5jNa3zU9uf0fVlvr1cjgr75dxq842D4jVQ5XVlZmaGjYvXL0\n6FF9fX1lZeXTp0/jFUiVAwIBDRr0C4kktj34bMKpQHfr6WQFqt36H/kVB9tnpMoxGIwjR44sXLjQ\nwMCAO09wcHBCQsKtW7fwrtfa2ko45xC8I0KQKjeiQIMGnyz8fgZ+jB1Ob/K0I1eSop6+2PTrUXEJ\niV6KA+7u3bsWFhaKior+/v5462lra9u6dau3t7eMjExOTs7kyZM1NDTS0tLq6uoOHDhAo9E8PT0l\nJCRqamoQQgUFBaWlpWw2m7tjp6urKyAg4Ny5c2PHjl2xYoW0tHR1dTXhnN2XkZqaam5uLisra2xs\nzM3QOHHihLa2tqysrIuLS1dXF2GlqKho0aJFFArFxMQkJycHIYRh2JgxY06cOKGnp6epqfngwYM1\na9aMHj3a2toaj7UrLy+/dOmStra2np7ew4cP8dfiNmjeCQk/pUH6OsCAg33Q4CNIJFJiaJA8dVR/\nJqkqLfnq558Gakm4tLQ0Z2fn8+fPT5s2zc7O7unTp6amprKyskuXLg0NDd2wYQOFQiksLBQREUlP\nT58yZQq+H7yysrK9vV1LSwshZGNjY2NjY2FhYWpqis958+ZNdXV1XV1dhBCbzW5ra1NXV5eUlOSd\nk7sMFotlZ2cXERExZ86c3377LTAw0MHBISkp6fjx4/fv3+dwOIaGhrm5uUwms0dFVVXVysrq4MGD\n58+f37Nnz/r169PS0hgMBoPBYLFYL168cHJyWrx4cUpKyuHDh3V0dF68eNHa2spkMiUkJOh0ekhI\niJubW3FxcVNTU1VV1YQJE6qrq3knJPyUBvaLAIMHGvSIUF9fX15e/tFh79+/p9PpPYqLFi2ytLTE\nb+vq6vIeh4ebO3fu3Llz+c2sS51Kp9MzMjL6vOSP8PDw6H+qHEKos7MzOzsbP8sZIZSamjp//nz8\ndmlpqYqKiqSkJOGc3JVAqhwYPNCgR4Tff//98YsiRaXejsdACGkaW+wPPcPv0YqiVz+uWeXq6kr4\naHx8/FAeqkGj0fqfKocQys/PV1NTo1Kp+N2SkhJ7e3v89u3bt2fOnAmpckCAoEGPCBiG2Tq5Gkzr\nV5hcypWLvRy9MMT/cO5/qhyOTqd3P7KwvLycQqEghDgcTkxMzKFDhyBVDggQ/Eg4Qn3piXT9TJXj\nzvP8+fPuf7VoaWnFxcVVVVXt2rVLUlLSysoKUuWAIAkmo2lIQJod1/bt23+Nv97PRLofDwbiV88T\nEv1JleMyNTXtfiXAgoICIyMjCoXi6OjY1NSEFyFVDggKpNmNCL/88ouC0Ux8F8dnJ9I5ePzLSEke\n/w0KADAEYB/0iAOJdAB8KaBBf6k6OjpWrVqFH7H7Uffv319q1PP0My5IpANAOEGD/lK1tLQ0NDRs\n27atL4PxE+f4gUQ6AITTMGzQXV1datq6Klra7S0tgl7L4FJQUOjjZVhVVFR6eRQS6QAQTsPwMDsO\nh6Oqret3Ol5CSlrQa/kyCE8iXQ+9B9RxHT16FD92mHu3RxYdv5GIKMpu6EFAHeBnGDZo8KmEJ5Gu\nh14C6rhjrl+/7u3tbWRkhN/lzaLjN5Iwyk4gIKAO8AMNeoQSnkS6HvoYUIcPzsrK8vPzGz9+PN52\nCbPoCEcioii77iCgDgiDYbgPGvAikUh/hJ9UUOxXSk5V2WuLf3kO1JII9T2gDiFUVVXl6uoaHx9v\namqKt13CLDrCkYgoyo4LAuqAkIAGPSJs3759bW0tv0cZDAaTyUQI7dy5c//+/YRj2Gz2kSNHbt++\n/ddff/GbR1lZ+bN/SBQREfHz8+t7QF1bW5uDg0NgYKCIiIiMjAweN0GYRUc4Etcjyo4LAuqAkIAG\nPSKQyWQymczv0c2bN5OUNEhipLEzvjr7H779V2mi4Xv+L5GblvrvX7Z+9v5cERERCQmJPgbUYRi2\nevVqNzc3a2vrCxcuGBkZ4X8x8GbR8RuJ6xFlxwUBdUBIQIMGiMPhrPrXVgkpqf5Mcva3PePGjePd\nGu27vgfUlZaWXr58OTU11dfXt7W1taura+HChTdv3uTNouM3En/FHlF2OAioA8IDfiQEBASSddf3\ngDp1dXUGg5GVlYXvPQgNDcWvA8ubRcdvJK5HlB0E1AGhI6CQpkH04cOHaTO/SnxZNcnMUtBrGURM\nJtPe3n5AplqwYMH5rNf9zLpb4uZOp9P7uZK+B9ThOjo6pKWla2pq8LuEWXSEI3E9ouwgoA4Im2GY\nZtfR0WFmM2935PndLg55T1MFvZzB0tDQ8MMPPyQmJvZ/qoULF7r8ehLfxfHZWXemNvP8N63vzy4O\nAEAPsA8a/I/Pz7oDAAy04dmgmTWMOxdj2lt6OehAqBUXF5eVlfU+5v3797W1tUlJSfwGkEikr776\nCv9trT/6mHXX1dnRzxcCAPQwDBu0hITE8UO/dXZ2pkmQBL2Wz7Rp0yYpdT1R0Y8cU0ybMC344hV+\nj+Y8fnjlQuzkyZP7uZg+Zt2RxMT7+UIAgB6GYYNGCH333XcIoeDgYEEv5DNxOBxnr239bHmnWb69\nXOO17yDrDgBBgcPsvgwCvMbrEGfdDXg2W1tb248//qiqqkqhUNzd3dlsNvchiLIDQg4atLDLffLo\nlO/PGQ/ucStsdtehzT/YrlwTeD2lvrryelQYv+KAGOKsuwGP/vH29i4tLX369GlycnJiYuKNGzcQ\nRNmBLwQ0aGH3uiBXXEJSslu29Yv0NAkpqW8cnRSVVRw9fkq+HM+v+NmGOOuOMJuNMMUtLS3tq6++\nkpOTmzlzJn6ZmJycnHnz5ikoKBgbG+fn5yOEMAxTUlKKiopSU1PLyMiIiYmJjIzU1NQ0MjIyMzPD\nU0cgyg58EYbnPujhZAiOeyORSOF7d5BI/fqPoTCLLun9OdvRhNlshCluWVlZixcvjoyMnDFjxrp1\n60JDQ52cnL755pvAwMCLFy/u3LnTy8vrr7/+Ki0tbWxsvHLlSmJi4sSJE1++fKmsrIwQevHiRVpa\nWnh4OIIoO/CFgAYtRJYtW4bHkhUUFPQyrM/HvXX28XXPnj3b3Nzc+5iWlpYlS5ZYWlryGzDT1Cgo\nKKiPr4jT0dHZvn07YTbbjh07eFPcDhw4sHPnzqVLlyKEDh48WFNTExgYuHnzZicnJ4SQq6srHmqR\nlZWloaFx4cIFSUlJhJCamlpXV1dUVJSfn190dDQ3zQ6i7IDwgwYtRN68ebNjxw6EUHZ2di/D+nzc\nW1+/3FGjRo0aNeqjw2JjY6WlB/IqYnJycvyy2XhT3DAMu3v37u+//46PnDhx4sSJE9etW/fnn3/i\nlcbGRj09PYRQZmbmsmXL8O6MEEpNTfXw8Jg2bVpaWhoeD42DKDsg/KBBCxFZWVl8l6i8vHwvwwR1\n3JuVldXAToj4JNghohS39+/fNzQ0KCoqIoQ6Ozv37NmzZcuW0tJSPL4OIXTx4kVbW1uEUFZWlouL\nC16MjY3dv38/vlekx0tDlB0QfvAj4ZdHaK/x+hn4ZbPxprjJycmpqqoGBQVVVVV5e3unp6crKipq\namqePn26oaHh2LFjSUlJ69atQ93S46qrq728vCIiIlRVVcvKysrKytrb27kvDVF24AsgmIymITFn\nzhxBL+HTcBc8d+7ci3nleJgc/kdRWYUbL5f4siog4Zb2+IlkBYq13YoLOWWExUWr/y87O1uw76gv\nCLPZCFPc/vzzTy0tLSqV6ujoWFtbi2FYSkqKgYEBlUr99ttv37x5g2FYbW2tgoICnh4XGBjY47/2\nFy9ecF8XouyA8BuGaXZc1tbWycnJgl7FJ+AueN68eT8cjejvmYT7fQ/8vHnq1KkDtDoAwFCDfdDC\niEQi/b5lYz9zjkoL8iR2/DxQSwIADD1o0MIoJibmo8e9NTU17dixo5e8ERKJpKWlNdBLAwAMHWjQ\nwqgvx701NDTIysrq6uoOzZIAAEMPjuIAAAAhBQ0afGG4+W3l5eXLly+n0WgUCsXT0xPPVm1tbXVz\nc1NSUtLQ0Oh+fVjc0aNH8aOSuXf19fWVlZVPnz7d+8ihByl3AEGDBl8cPB4IwzAnJyczM7PCwsJH\njx5FRES8ffsWw7Bly5apqKjk5+dv3bp13bp1+MHIuOvXr3t7exsZGeF3g4ODExISbt26hXe91tZW\nfiMFAlLuAEJwHLQw+aQFD+BVvYVfUlKSubk5lUr18/ObMWPGrVu3Wlparly5ggfCxcbGGhsbczic\nx48f4zcwDGOz2fLy8kVFRfgMmZmZRkZGEyZMuHbtGoZhnZ2dWlpaJSUlGIZxOBwKhVJcXEw4srtH\njx6ZmZnJyMgYGRldunQJLx4/fhw/A2X16tWdnZ2ElcLCwoULF+KRe/jB6RwOR1VV9fjx47q6uhoa\nGvfv33d1dR01atScOXPYbHZ7e/uoUaN2796tpaWlq6v74MED/LUMDAxyc3MJJyT8lAbvGwFDAxq0\nEIEGTejJkyfKysr37t1jMpmzZ88WFRWtq6vDHwoJCUEIUSgU/LyVwMBAV1dX/KGKigpxcfGOjg4M\nwyorK6dOnVpQUCArK4ufz3L16tWZM2fiIzs7OyUkJFgsFuFIrvb2dhqNdvXq1ebm5h07duBPv3Pn\njoGBQVlZ2evXr+Xl5TMyMngrDAZDSUnp9OnTzc3NW7ZsMTc3x18IIXTo0CEWi2Vvby8nJ0en0+vq\n6shkcm5ublpaGkJo37599fX1+/bt09PTwzCssbGRTCZ3dXURTtjLpwS+XHAUB/h8bDb7wYMH3a9R\nQigxMfH169fi4nzPu5GTk+MXP2JoaHjlyhV++W3u7u62trbu7u5BQUF79+6l0Wh37tzJy8uTlJRc\nu3atvr6+uLh4W1ubg4NDYGCgiIiIjIwMHmSRmpo6f/58fJLS0lIVFRVJSUnCkVyQcgeGHjRo8Pmy\ns7M3bdthPPvrj4yTHaU0he9Rg12dHXX5Gfv37yd8VFpa2svLq0d+24cPHzIyMvD8Iy0tLVtb21ev\nXiGEHBwcrl27ZmlpaWZmpqOjo6GhgWHY6tWr3dzcrK2tL1y4YGRkhGdIlZSU2Nvb43Pevn0b3xwm\nHMkFKXdg6EGDBp+Pw+GMnWJk94Nnfyb50N4Wu+tf/A7oJoy7YzAY8+fPb2pqwhvcs2fP8C3HioqK\n2NhYEonE4XAsLCz27NlTWlp6+fLl1NRUX1/f1tbWrq6uhQsX3rx5s7y8HI/B43A4MTExhw4d4jfy\n70VCyh0QBDiKAwykAb+4LWF+m6qqamdnZ0hISENDQ1hYWHp6+qpVqxBCpqamUVFRTCbTx8dHWVl5\n3rx56urqDAYjKysL36UQGhqKH3unpaUVFxdXVVW1a9cuSUlJKysrfiMh5Q4IkgD3fw82+JFwsD17\n9mzpWg88S8//zMVvHJ0QQtzUvYv5FbQx6h57D4ffz5gyw8rJazthMS6zeNGiRb28CmF+27lz58aM\nGUOlUpcuXVpZWYmPPHv2rLKysqKiooeHR4+Yt46ODmlp6ZqaGvxuQUGBkZERhUJxdHRsamrqZSSk\n3AEBggYtRL7oBu2y1XeBs5uklDS3Qfufuaimq4/f/nd0gqqWDmHxow0agBEL9kGDgfHZF7fFhm/g\nLQD9BA36C5abm+vu7i7ABdTV1SGqCr9H+3hxW1ZbK78ZABjhoEF/qRQVFS9duoQfEiAoubm5Z67d\n5vdoHy9uKyUtMxRrBeALBA36C4YfmyVADQ0NvTza14vb9u+6BAAMY/D/Bhgsg3Fx28+OsisrKzM0\nNOwxG2ERQZQdEBrQoMFgIZHEtgefTTgV6G49naxAtVv/I79i331GlB2DwThy5MjChQsNDAy48xAW\ncRBlB4QHNGgwkMLvZ6jp6nPv6k2eduRKUtTTF5t+PSouIdFLsRd37961sLBQVFT09/fHW09bW9vW\nrVu9vb1lZGRycnImT56soaGRlpZWV1d34MABGo3m6ekpISFRU1ODECooKCgtLWWz2SYmJtw5CYsI\noaysLD8/v/Hjx/M26NTUVHNzc1lZWWNj44SEBLx44sQJbW1tWVlZFxeXrq4uwkpRUdGiRYsoFIqJ\niUlOTg5CCMOwMWPGnDhxQk9PT1NT88GDB2vWrBk9erS1tTWez1deXn7p0iVtbW09Pb2HDx/ir8Vt\n0LwTEn5Kff7SgPCCfdDg84mJiT1L/quu6m1/JmGzu9TkpPk9mpaW5uzsfP78+WnTptnZ2T19+tTU\n1FRWVnbp0qWhoaEbNmygUCiFhYUiIiLp6elTpkzBAzQqKyvb29vxSzLa2NjY2NhYWFiYmppypyUs\nVlVVubq6xsfHm5qa9mjQLBbLzs4uIiJizpw5v/32W2BgoIODQ1JS0vHjx+/fv8/hcAwNDXNzc5lM\nZo+KqqqqlZXVwYMHz58/v2fPnvXr16elpTEYDAaDwWKxXrx44eTktHjx4pSUlMOHD+vo6Lx48aK1\ntZXJZEpISNDp9JCQEDc3t+Li4qampqqqqgkTJlRXV/NOSPgp9edLAUICGjT4fNOmTcujP/vosHXr\n1nV2dioqKvIbICYmRni8oJ2d3dGjR/sTZYeP7OzszM7Oxk+t5upRhCg7IISgQQ+FjIyM58+ff3RY\nVVVVWFgYv0dFRUVXrFgh2OPqehAREaFSqR8dFhQU1NbWRiKRPnV+KpXKG9LW9yg77jz5+flqamo9\nltq9iEGUHRBK0KCHwpEjR+THG8rKK/Q+bMlG78J2vo+m/XVdR0cHDzb7svTYGu27fkbZceeh0+k9\n9jX3KEKUHRBO0KCHiMmcuYpKyv2ZoXtK3AjBDWlbtmzZwYMH6XT6b7/9xo2y+/777xMSEtLT0yMj\nIxFCpqamv//++9KlSw8cOIBH2XHnef78Oe8+2e5FPMoOv7106VJPT88lS5YghFJSUqhUKolE4kbZ\nhYaG4lF2mZmZNjY2T5480dbW5kbZ9aigf5Lntm/f/vDhw/Xr19Pp9FGjRmVmZu7duxf1GmVnaWl5\n+/btw4cPZ2dnI4SePXsWEBBAOKGamhrvpzTI3wwYKgLKABkKwhOW5OTkFP4gE08Iwv8Yz7bhhgol\nvqwKSLilM3GyrLy8td2K81mvCYvfb/ZJSkoS9FsZagMSZWdqanr37t0eMxMWIcoOCBURbPhG1Vhb\nWycnJwt6FQgh5OzsbL3BB9+Czn3y6NHNK0mX4oJuPsCPSGOzuzznWjp4/GQ82yZo++YpFlbL3Dfz\nFjls9rpv53+JuzgAAJ8HjoMeaq8LcsUlJCWl/ntg2Yv0NAkpqW8cnRSVVRw9fkq+HM+vCAAYUWAf\n9FD77FjOoV8qAECwYAt60Pn5+d25c6eXAX2M5ezs6Bj0tQIAhAk06EGHYVjvsXN9jOUU++e0CwDA\nCAENWvCU1bXeFL3Cb/9vLOf/Fv/31Akg/HrP3sP1SNSztrYW6UZeXp7D4RA+vZc5hxhk7w0eaNCC\nNxixnEAY9JK9h/gk6p07d670H46Ojv7+/iIiIoTRfYRzCgRk7w0eaNCCNxixnEBQ+pi9h/gk6qmr\nq2tra2traz979kxRUdHLy4vw6fzm5ILsvWFCkAdhDzIhOVHF19f3m2++6XGiymf8GZknqnxZnjx5\noqysfO/ePSaTOXv2bFFR0bq6OvyhkJAQhBCFQqmtre3+FHNzc96v9dWrV1ZWViwWi1shfDq/Odvb\n22k02tWrV5ubm3fs2DFz5kwMw+7cuWNgYFBWVvb69Wt5efmMjAzeCoPBUFJSOn36dHNz85YtW8zN\nzTEMq6ysRAgdOnSIxWLZ29vLycnR6fS6ujoymZybm5uWloYQ2rdvX319/b59+/T09DAMa2xsJJPJ\nXV1dhBP28imBHuAwu6EgKSm57wcnMbF+/crXWF/rvXr5QC2JV2lpaUlJSe9jMAzz9/fvJSBJRERE\nRUVFdERexWrNmjX+/v59zN7Di4Qxe21tbatWrYqIiJCUlOQWCZ9OWESQvTeMQIMeCtu2bZs1a9ZH\nhwn21EcfHx+SqraEpFTvw3RnWPfyaNajlDVr1gj2ciSCQqPR+p69hyOM2du9e/fXX389depUhBDh\n03ufE0H23jACDRr8jcPhfOuyjvyxyL3etbe8V1VV1dXVHahVfUE+KXsPxxuzV1lZmZCQwN1XS/j0\n3ueE7L3hZCT+UxT00f71q7pH6JXkZf9sP8/FbPyJX7w6WKxeiiMTN3uvvr7ex8eHTqdPnz6dm73X\n0NAQFhaWnp6+atUq7lN4Y/b279+/adMmeXl5/C7h0/nNmZKSkp2dXVRUxM3e8/X1xbP38vPzra2t\ns7KyWlpauNl7PSron6i82traxMREExMT/KKOmZmZeC/uJXuvubn54sWLhw8fxlNeuVvQvBMSfkpD\n8wV9kQS5A3yQCc+PhA8ePOjLSMEu2M7OLiq9AP9B0v/MxW8cnRBC3Mi9i/kVtDHqHnsPh9/PmDLD\nyslrO2FxmftmbsTaCNT37D1cj0S90tJSBQUFJpPZfQzh0wmLkL03/ECa3aDbvXv33LlzhX8ftL29\n/Xe/HMB3cVw9faqu8u29xAuHLt/GI/dynzwK37sj6OYDhFB++uNTvj4nbqfyFi1tl2xxWWFpaSmo\ndwHAcAL7oAEBSHQCQBhAgx5Brl+/HhMTQ6FQCB/NzMz8jv9zIdEJgKEHDXoEqaysNDMzw3+R51VW\nVtbLcyHRCYChBw16BBEXF1dSUuJ3DJzsP9vChJTVtd4UReK3/zfR6X+LkOgEwMCBw+xAn0CiE/pY\nOl1ra6ubm5uSkpKGhkZsbGz3J/aIrGtsbFy3bp2SkpKCgsKGDRvwH+oJc+yG9v39DdLphAdsQQ83\nZ8+ePXHiBOFD9fX1IiIix44d61EPCgr66HEXeHhT8C7vkN1bp9vM757o1L146eTvA/IuhFP3dDo7\nO7vQ0NDKykoTE5Nt27ZpaGgsW7bMyMgoPz//woUL69at+/bbb+Xk5BgMRlxcXGRk5OTJk7nzrFmz\nxsDAgE6nt7a2mpmZbd68eeLEiefOncPjihBCW7dutbCwENQZ8+np6Z6enj2Kz5496+P1MCH8aADB\nFvRwk5+fv+KXfTtir/P++f32kyN/Pu5RNLdz4hfBEX4/Az/GDqc3edqRK0lRT19s+vWouIREL8Xh\npI/pdGlpaXV1dQcOHKDRaJ6enhISEjU1NYhPZF1aWtqsWbNUVVVLSkpoNBp+oh1vjl33ZUA63cgE\nW9Dgb2JiYid/+UlMvF9NtqLopdR614FaksClpaU5OzufP39+2rRpdnZ2T58+NTU1lZWVXbp0aWho\n6IYNGygUSmFhoYiISHp6+pQpU/BdJ57cngAAIABJREFU8JWVle3t7VpaWgghGxsbGxsbCwuL7mcM\nuri4LF26VERERE5OrqCgAD83D1dYWBgUFJSUlNR9bz6LxbKzs4uIiJgzZ85vv/0WGBjo4OCQlJR0\n/Pjx+/fvczgcQ0PD3NxcJpPZo6KqqmplZXXw4MHz58/v2bNn/fr1aWlpDAaDwWCwWKwXL144OTkt\nXrw4JSXl8OHDOjo6L168aG1tZTKZEhISdDo9JCTEzc2tuLi4qampqqpqwoQJ1dXVvBMSfkpD+C0N\nZ9Cgh7Ond27F/H7gHbPexHquu/9BSen/NoLK18U+y2zjMv97Jnd0dHR7e3vvE3Z2dlpYWMyePZvf\nAAOVmby7L79ECgoKAQEBveSu9UiSo9Fod+7cycvLk5SUXLt2rb6+vvg/B7T0iKwLDw9/+PBhcXHx\n6NGjFy5cGBYW5u/vjz9EmGOHIJ1uBIMGPWxVvi4+udNr24nTGmPHh+z2iT16cO2Of+MPcdjskzu8\nOjs+dB8vJSUlJfWRKDuE0F9//TUS0kTFxcUJc9f4Jck5ODhcu3bN0tLSzMxMR0ene3x+98g6Doez\na9eulJQUPT09hJCTk1NKSgp3ZPccu+4gnW7EggY9bD1PuWP29YJJZpYIoVVbduxY+S23Qd+MOT1K\nZUxRbtZnTKuvr//xQcPCJ6XTVVRUxMbGkkgkDodjYWGBZwbhukfWNTQ0MJlM7pGOr1694h7d0SPH\njgvS6Uay4b8pNPywWKyioqLXfDQ3N+PDujo6SGIk/LaYuPj7psYP7W0IoeqKstsXold57xDYG/hC\nfFI6nampaVRUFJPJ9PHxUVZWnjdvHnee7pF1ioqKVCr10KFDDQ0NcXFx169f5x671iPHDkE6HUCQ\nZjf4BjzNztfX19LS0pGPcePGBSTcSnxZdfDiDbK8wm+Xbp5Ny7dc8C1C6Mzj3IQXbyebz9wREn3m\nSZ4oiZT4smrjviPR0dH9e4vDVt/T6c6ePausrKyoqOjh4dEjnq1HZF1SUtKECRPk5OTmz59fWlqK\nFwlz7CCdDkCa3aAb8DS7kJAQCoXy/fffEz7q4+OjZPG13uRpCKHbF6ITTgVyutgLXf7vQtCh+Jyy\nu4nnXzxL+9ehE+8aG/7Patql/Dd3E85baNJWr179qe8LADDYYB/0sNXe2jLdet78710QQq+y6Bp6\nV0VJpNy0Rxn3762xmIRxOBw2e43FpHkrViNNmqAXCwAgAA162GJWM7Y6LPCNiKONUYs5sn/+SheE\n0A+7f8Wve9LS3ORtN/fIlSR6SpKgVwoAIAYNethS1xu7+uedAT/+n5i4+FxH53nfuyCE5Ch/X59U\nXFJSRERklIoqSQz+GwBASMH/nEIEw7B///vfbf8EePKTmZkpJSWVlUV8kNyDBw9mKGkhERGEkIGh\nya7wv1N7Xr/I7THy4MUbJfk5dVVvYRcHAMIJGrQQYbPZl27+5bLVt/dh35hZ9/LoZDFy/O/7CePE\nXr58SSKR8LMVuKgIwXFRI8Hly5fPnj177dq18vJyHx+f5OTkzs5OZ2fn48ePi4qKEhbb2tq2bduW\nkJDQ3t6+YsWK4OBgEonU2tq6adOmGzduSEpKHjx40NnZmfsSZWVl3333Hb9Nh6Hh4eFhaGjo7u7e\nvXjs2LFXr14FBwd/9OncT2nQFvhpoEELF1l5eb1JPU8k+yR6k6a+Tr1LeL51ZGSkhIQEHLAxMvUe\nxcdb1NTU9Pb2fvPmzdOnT5lM5ty5cxctWrRkyZK+h/YJxDCL4oMTVYTa/vWrKl//Ny6jJC/7Z/t5\nLmbjT/zihf/Wx68IAOpzFB9hsbm5OSYmJjIyUlNT08jIyMzMjMlkflJoHxdE8X0+wR6GPai+uBNV\nrKyszOcuSHxZlfiyyv/MxW8cnRBCQTcf4JWL+RW0Meoeew+H38+YMsPKyWs7vyK/Nx4REQHnpIwc\nT548UVZWvnfvHpPJnD17tqioaF1dHf5QSEgIQohCodTW1nLH9yiyWKy3b9/iD+Xn51Op1Ldv3wYG\nBrq6uuLFiooKcXHxjo4O7gzm5uZJSUk9ltHe3k6j0a5evdrc3Lxjx46ZM2diGHbnzh0DA4OysrLX\nr1/Ly8tnZGTwVhgMhpKS0unTp5ubm7ds2WJubo5hWGVlJULo0KFDLBbL3t5eTk6OTqfX1dWRyeTc\n3Ny0tDSE0L59++rr6/ft26enp4dhWGNjI5lM7urqIpywl09JGMAuDiH1uiBXXEJSUkqaW3mRniYh\nJYV3bUePn075+ixz30xY5Dcng8FITk5+9OjREKwfCNbu3bv7HsVHWJSUlFRTU+vq6oqKivLz84uO\njlZTU+t7aB8XRPH1BzRoIbV0rQdC6GnSLW6ltvKNlsEE/LamwYS6qrcYhhEW+c3p5eW1YsUKEok0\nmAsHQkFJSamPUXz88vkQQqmpqR4eHtOmTUtLS1NXV0d9Du3rDqL4+gMatBCpra2lasnxe/RdY4PU\nP9d1lSGTuzo7Wa0thEV+M8jKyvY4hAMMV32P4uOXzxcbG7t///7IyEi8d+P6GNrHBVF8/QQ/EgoR\nfA8aP2QFBdY/h0i3tbSQSGJSMrKExUFfKBB6fY/iIyxWV1d7eXlFRESoqqqWlZWVlZXhF3PoY2gf\ngii+AQIN+ouhrK71pujvf3u+LSmiqamLiIoSFgW3RiAsxMXFQ0JCfH19bWxsDAwMpKSkjI2NJSUl\nw8PD9+/fr6+vf/Pmzfv370tJSREW4+Pj6+rqZs6cqfOPsrIyhNCxY8d27NhhYGDQ3t4eFxfXfXfZ\n8+fPu29B//DDD+np6ePGjfvqq6+0tbW/++47ExOTsWPHBgQEmJqarly5ctasWbq6ukwmMzg4mLeC\nEAoJCUlOTtbS0jp8+PD169e1tbXfvXtXW1trYGCA/mnK+GvhW9B0Oh3f/lVTUzt27Nh//vMfFRWV\nmpqalpYWPICbd0LCT2novqQ+gDS7Qdf3NLvGxsYFK5y3Ho/kVn74ytj/zEX8yq1sdtcG6+lrd+41\nmmUd+LOn7sQpyzd5ExaPb3AShjcOAOgn2Nr6YpBIYtuDzyacCnS3nk5WoNqt/5FfEQAwPMCPhEIt\n/H5G97t6k6cdudIzfI6wCAAYBqBBCxEpKSnx9venNrn0Puzt27fi4uLKysr8Bujo6Az00gAAAgAN\nul8qKirwc1J70djYWFVV9fr1a34DFBUVKRQKQkhaWpp7cmover+iCgBg2IAG/fkKCwu/mjtfd+KU\nj47MKC5HZ2MJH8IwjCLKvnnz5kCvDvTUe5wbYUgbv+S2Cxcu/P777wUFBY6OjuHh4SQSqbGx0cfH\n59q1ax8+fFi5cuWpU6dEREQE8jaHWZzbCAcN+vN1dHQYzbJe57u/P5Nw2OyQzWsGaEWgN73HufGG\ntJHJZMLktsTExBMnTpw7d45EIllaWt6/f9/GxmbNmjUGBgZ0Or21tdXMzGzz5s34wbxDb5jFuY1w\ncBTHQILwOWHTxzg3wpA2wiKbzd65c+eFCxfGjRunr6+fmZlpZmaGEEpLS5s1a5aqqmpJSQmNRsNP\nV+OCODfwmQQU0jQUBjvNLjc3d4GzWz/D5y7lv/n6668/6XVPnTp1/vz5QXpTw0nf49wIQ9oIizdv\n3rS2tnZyclJQUDA1NS0sLMQH/PzzzwghEREReXn5ysrK7suAODfw2WAXx8D47PA5u3U9/zX6UVVV\nVYmJiXAqSu8WL1584sSJPsa5EYa0ERZv3LhRUFCwY8eOI0eObNiwwdfX98KFC+Hh4Q8fPiwuLh49\nevTChQvDwsL8/f25K4E4N/DZoEEPjMEIn+Nn69atDAYDQul6p6io6Ojo2Jc4N8QnpI2wWFhYuHfv\n3m+++QYhtGrVqlOnTnE4nF27dqWkpOjp6SGEnJycUlJSuq8E4tzAZ4MGTeyXX37BD33rBX4tCX76\nHD7X+qlrI5PJEEr3UX2Pc0N8QtoIi2VlZdwPPzMz09LSsqGhgclk4mkPCKFXr14ZGhpylwFxbqA/\n4EdCYjdu3ND9GDwhl58+h8/JDPqbGZH6HueG+IS0ERZ1dHSCg4Orq6uvXLkSHx+/efNmRUVFKpV6\n6NChhoaGuLi469evb9y4EUGcGxgQAtz/Pdj68yNhX57b/UdC/I+isgr3R0K/0/Ea+uPw2/tir6ho\nahMWP+NHQtBHMTExNBptypQpoaGhZDL53bt3GIadO3duzJgxVCp16dKl3F/zzp49q6ysrKio6OHh\ngQ/jV8zPzzc1NZWTk5s7d25RURFeTEpKmjBhgpyc3Pz580tLS/Givr5+WFhYR0fH4sWLyWSysbHx\n5cuXx44du2fPno6ODjc3NzKZTKVSly9fXl1dzVvBMCwvL8/Q0FBKSsrCwiI1NRXDsObmZjk5OTab\njS/PxcUFf61FixZdu3bt7t27q1evdnZ2lpWVtbS0TE9PxzCsurqaSqVyOBzCCfl9SkBIQJrd/7d3\n53FNXHv/wE/YArIFagwUilhQ4bYKKAIVtYJV3BVErWAVvFakqI+IuK/X5VHRK7ggCKiobIq4VS0/\nsOKCBGgQy6YgBhEMEIJAC6aEZH5/zO08uUkIUZYM8H2/7h/k5MyZyXD9NpyZ85lP37awsHDToX+L\n3wf9CeFznv7rI9b5pKdDmAYAQBJMcfQUCJ8DAHQRXCTsTp8QPicSCnvjyAAAfRAU6E9HoVBePvvt\n8rEuLfXGMExZoQ0AAJKDAv3prKysIsP+3Wm32NjYMWPGjBrVYaYSpIMCAGRT9lXKHtTTd3EoaOfO\nnY8ePequ0YBSXLt2bc6cORiGVVRULFy4cPDgwfr6+j/99BN+QwWXyxVfN3T06FFiQzabbWNjQ7yc\nPHmy+L8+/JaMhoaGf/7zn3Q6XU9Pz8/PD7/jQilWr14dEREh0RgaGurv76/I5sRZAt0FLhIC0Dnx\nJDwHB4fS0tInT55ER0dXVVUhhK5du7ZkyRL231avXo0Q4nA4x44dmzlzJv6QU9ylS5eIbgsXLtyz\nZ4+KioqPj4+BgQGLxcrOzo6Pjy8pKVHWx8zJyXF2dpZoJG6j7hQELXU7KNAAyKZgEh5CKCkpydfX\n1/xv2traCKGSkhI2my0UCsWfdW1qaor3yc3NNTQ0DAwMRJCEB+RQ9lf4HgRTHOCTKZ6Ex+Fw1NXV\nXVxcBg0aNG7cuBcvXoiP4+jomJ6eLjH4y5cvJ0yYwOfz8ZeQhAc6AhcJQe8pKCh48eKF/D5CoXDf\nvn36+voddaBQKMbGxniYRk/Q0NCIiYmRk/EmkYSXk5Mzffr0gwcPGhgYBAQEbN++nfieKxAInj9/\njq/MJrS2ti5dujQ6OppKpSKEIAkPyAEFGvSerVu36lr8Q01dQ363MTPc5byb9+j+5mXLvv766249\ntP+joqKirq4unfHWURLe3Llz586di/dcvny5+GOlioqKTExMDAwMxMfftWvXlClTRo8ejRCCJDwg\nHxRo0Kvm+vpRtbqUD/XH+4bPP/+cSI/rCYon4eEBST/99BPemJ+fj1dwHIvFEp+ARghVV1cnJycT\nE7iQhAfkg4uEQJnI+ZAwxZPwVFRU1q9fHxMT09jYmJiYmJSUhH/Jxf3222/29vb/9XkPHFizZo2e\nnh6xI0jCA/IodQa8Z8FFQrKZNWtW/LNXXXxI2LwV/rm5uT19qIon4Z0/f97MzIxGo82ePbuyslJ8\nEHt7+/v37xMv2Wy2vr4+j8cT7wNJeEAOSLPr/m0l7Nq1a+rUqRMnTuyW0fq02bNne+8Pw6c4bp47\nw62u+vVaYkhKKp7/V5D1JGrfthN3HyGEinKentkZfCo1U7rRYcr0Hf4rJL6ZAtAvwRw0UI5PfkhY\n7x8qAMoCBboPa21tnT59upGRkbIPRFF5eXneHb+r4EPC2gVtPX+kAJACFOg+rKWlRU9P79ChQ8o+\nEEWtXLlSzrsKPiRMVU29xw8UAHKAAt23DRo0qEdvOOteg+Q+gJFhOvRtWQz+c1V5Gd3ElKKiIqMR\n0lnBgAG32QGy+Mrxm5bmxqzUn/mtLTdjwr+du6CjxgErJSUFXxTz5s2bRYsW0el0Go0WEBAgEokQ\nQvX19WpqapS/HTt2rKOeLi4uFDF6enp4O0KooqJC/F5spfD394+MjJRoDAsLw+9B7BRxlvoBKNCA\nLOAhYZ362FC9jnrKDNWTGb+nFBCqR4ACDZQp6mEefo8dDn8eWGx28Zr/Pa6uoSGnceDoSqheRz1l\nhurJjN8jQKiecij1Luye1e8XqtTV1S1cuLAnRu4h7u7uXzs6jx4/sSv/G2LyRVFRkbI/Si/pllA9\niZ4EiVA9nMz4PQjVUxa4SEhS+fn5ZWVl8vs0NzdXVVVdvXq1ow5qampz584Vf9iHcsXFxfE7W6vd\n2tr6zTffSCTAiftyyGfieW/9lbGxcVhYWLeE6kn0xBslQvVwMuP3EITqKQ8UaJIKCgoa/q1bp3cs\n2Mxw//VFRUfvPv75xsiRI/FwBjLQ0tLS0tKS38fAwCA3N1dj4E1lSFBVVZUZF6d4qF5HPXHioXoE\nmfF7CEL1lAcKNHlNXehN6Vrq8ZuXJVgfXMrPYDCUfQik0MVQPZk98ZElQvUI0vF7CEL1lAouEvYZ\n5Ax+Az2ni6F6MnviI0uE6hEk4vcgVE/5lDkB3sP69EVCV1fX5OKqLga/zfD2LSws7K4PAnpfF0P1\nZPaUGaqHk4jfg1A9pYM0u+7fVsKnpdlNmTLlp1MX8SmOTw5+s50wOWRr0FdffdUtHwQA0MtgDroP\ngOA3AAYmKNCk4OrqKnFjUGFhoZz+Cge/CXrogAEAvQAKNCkIBAKJ8IH58+fL6a9w8Bv8fgHow+Af\nMCmoqalJ3HyqJre2QvAbAAMB3GbXJ/Wb4Df58WwIoePHj1taWjIYjHPnzhFbJSYmOjg46Orqrlix\nQigUIoRaWlp8fX2HDBnyxRdfxMXF4d1kprspBcSzgU8DBbpP6jfBb/Lj2cLDw5OTk+/du4fXspaW\nFoTQtWvXTp06denSpWfPnv38888PHz7EMGzBggVGRkZFRUWbNm1auXIlfr+tdLqbsj4mxLOBTwMF\nus/oN8FvCsaztbe3Hzly5NKlS8OHD1+8eLGWllZNTY1QKNy+fXtiYuLIkSMtLS2fPXvm4ODAZDK5\nXO7BgwfpdHpAQICGhkZtbS2Sle4mfhgQzwb6AKXehd2z+tBCFendTZ8+fdQ3E7oe/Pb69evu+iDd\nQvF4tps3b+KpaRiGCQQCDQ0NPp9/9+5dFxcXLy8vfX19e3v70tJSDMNCQ0OXL1+O96ysrFRXV29r\na5OT7oZBPBvoI+AiIUndu3ev0z5cLjcgIODKlSud9kxNTW1ubpbfJzc395dfftHR0emog6amZldC\nwlRUVE6cOKF4PFtmZqabmxv+FpvNNjIyolKpd+7cKSkp2bZt27Fjx1avXr1z587ExEQ6nZ6WllZY\nWEilUlesWGFpaamuri4n3Q1BPBvoI6BA93/t7e0rA9bOXv5jJ/106c4Lf5Dz/v9LvHTt2rWuHIme\nnp7i8Wzl5eUeHh54z9TUVHwOt7S0dN++fd999x1CaOnSpWfOnEEIeXp63rp1a/z48Q4ODsOGDcMz\n6aXT3cSPBOLZQJ8ABXpAMGQYTV20tIuDPP75hnQQ5UdRPJ4NIfTmzRsajYYQEolEly9fDgkJQQhV\nVFTgNREh9OzZs/HjxyOEKisr4+Li8HxOJyenvXv3ykx3Iw4D4tlAXwEXCQcopWTjKR7PhhAaOnRo\nfHz8u3fvduzYQaVSJ0yYgBAaNmxYeHh4TU3NjRs3kpKS1q1bhxCyt7ePjY3l8XjBwcEMBmPatGky\n090QxLOBPkepM+A9iwwXCf/1r3+Zm5vb29t/J5eBgYGcd52dnYuLi6UHV/CRVwKBYPT4iXgMXley\n8b5yGN/1E6J4PFtJSYmdnR2NRlu4cGFjYyPeWFRUZG9vr6urO3Xq1LKyMrzxwoULDAbD0NDQ39+f\nSDKTme4G8Wygb4E0u+7fVtz3338/bd022mf0rgySdOqY/4LZkydPlmhX8CJhe3v72G9dd59Lwl9+\ncjaeIcO4MDuzKx8EAPBRYA56wPnkbDxDhnHvHy0AAxkU6N52YNVSny17iCUn5YXPz+wKrquqdJgy\nY9Xu/9XQ1JRu7GgoCoXCYrEWLVokf4+d/pGkYDaeSChU8DMCALoFFOjeU5D15MndG3mPfvXZsgdv\nEQrbQ9b96Om/fswk1xNb1t2OPbvAb510Y0cDDh48OC8vj4it6Eh7e/t38+XlciiYjadCmqeDAzBA\nDLgCXVhYWFJS0mk3Lpd79epVOR3c3NykH+km3+uSAnUNKlXz/x5rXZzD1NDUxC/ZLfRff2Zn8AK/\nddKNE+d4dDSmvr5+p/vFFyjLoWA2Xqc7AgB0rwFXoENCQkSDTXVoNPndJi5e/uuLio7effb4gZ6e\nHrHOTUGkfTAKEYNnN9FFOhuPaCzMyerpI5EvJSXlwoULt27dIlqOHz/+6NGj69evEy0VFRXz58/P\nz8/HX7q4uGRkZBDv6urqNjY2vn37Njg4+MGDBwKBwNvb++TJkyoqKq2trZs3b05OTv7w4cPixYvD\nw8NVlfQXg7+/v62trZ+fn3hjWFjYy5cvJZbbyCR9lkDfNeAKNEJo4hz3wcYmXRnhj8b33XIkCk7+\nCtraumV3HcFj8MJ3BEXs2jTO1U08G0+8UekFWiIz6Pbt20FBQXv27MFfcjic+Pj4mJiYr7/+muhz\n6dIl4g+ITZs2OTk5USgULy8vd3f3yMjI6urqsWPHbt682czMLCgo6O3bt9nZ2Tweb+rUqbNmzVJW\nwmdOTk5AQIBEY25uLr7CpVOQrNSfwEIV5SzZwCk4+aumrt5deyT0lWy8jkLd8vPzd+/ebWVlZWdn\nh7eUlJSw2WyhUDh27Fhic1NTUzzKLjc319DQMDAwUGZyXlNT0+XLl2NiYszMzOzs7BwcHHg8nvhh\nQPQdUIqB+A2a8MlX7brrSSW982AUCoVS+fLF3hWLu3i0H/78o4sjfCwmk+nt7Z2QkGBjY+Pu7p6d\nnW1vb48Qevfu3fLly5OSkuzt7YkC7erq6urq6uTkhPcRV1paeuLEifT0dAqFoq2tPW/evMjIyNWr\nV9NotNLSUgqFoqmp+eLFCwaDgRAqLi5mMplRUVHE5nw+393dPTo6evLkyYcPHw4NDfX09ExPTz95\n8uTDhw9FIpGtrW1BQQGPx5NoMTY2njBhwqFDhxISEvbu3btq1Somk8nhcDgcDp/PLy4u9vLymj17\ndkZGxtGjR4cNG1ZcXNzS0sLj8TQ0NFgsVkREhK+v76tXrxobG9+9e2dtbV1TUyM9YEdnCfQDA7pA\nf/JVO9cF33fLASg4+dvFpUSqqqrv6+s67Xb8+PGEhAQ8FEI2E6NOb+nrRu7u7hcuXJAOdWttbfX0\n9AwNDaVQKIMGDcIzK3ACgeD58+f4qmtCa2vr0qVLo6OjqVQq0SiRnEelUk1MTNrb22NjY3fv3n3x\n4kXxYSH6DijLgC7QSr9qp+Dkb8rZk921RznWrVvn4+PTCztSEJVK9fX1lQh1wzDshx9+8PX1dXFx\nSUxMtLOzE//zoqioyMTERCLRadeuXVOmTBk9ejRCqKPkPIRQZmamv7+/jY0Nk8k0Nf2vW1Yg+g4o\nS38u0G/fvp06dapEY1FRkbNvh8+CUviq3V+KH8ablyU8fQ7xckv4Bf6H1vKi/8weIgrlpwP/xn+s\nLHshs/F9Xa3iu/tkqqqqXQyr614yo+/YbHZKSkpmZubOnTtbWlra29tnzpx59+5dfBMWiyU+AY0Q\nqq6uTk5OJuZqO0rOi4uLO3DgQExMjHjoHQ6i74AS9ecC/fz58zap+x/kP6ZT4at2il4x+/7775lM\nZqflPCEhgfhnKW3sl6ajRo1ScI/9BhHqtmDBgkOHDrFYrMOHD5uamnI4//mv3bx58wICAubMmUNs\n8ttvv0lMvx44cGDNmjXEHetEct7333+fnJyck5MTExNTU1MTGBh448YNY2PjiooKhBCDwdDS0srI\nyDAwMFBVVSWi7yIjI/Hou2fPnrm6umZlZZmbmxPRdxIt6O+kui1btjx+/HjVqlUsFuuzzz579uzZ\nvn37kNzou/Hjx6emph49evT58+cIodzc3CNHjsgc0MTERPos9fwvB/QW5WQ0Kc+yZcsiH+QS0W7X\nXrwzZBgRuW67zyV9YTkS/3l/3A0jM3OZjd4btv7yyy/deFTd+ISt/kROqFtbW5uWllZtba14f3t7\n+/v37xMv2Wy2vr4+j8cT7yOdnBcaGirxjwLPDoToO6B0/TnNTqbly5c7+64Vvw/6x2/H7Dl/Bb/n\nTChsX+0ybsX2fXYTXUI3Bnz5j1GL1gRJN6ppaPww7duPXagiR3eF5wEA+hO4D/q/4Bfoks+E+rmM\n09E3EL9qJ9EIAAA9rT/PQSso6mGe+Et8dYZEH5mNAADQowZcgdbS0grdGKAudkvsJ6irert20bzu\nOiQAAJBpwBXoiIgIRbrBpDAAQOlgDhqATqSkpODBSW/evFm0aBGdTqfRaAEBAUQS9/Hjxy0tLRkM\nxrlz58Q3rKiosLW1JV6+f/9+5cqVQ4YM0dfXX716tfj1eYmeSuHv7x8ZGSnRGBYWJv/OVAJxlkA3\nggINQCfw+CEMw7y8vBwcHEpLS588eRIdHV1VVYUQCg8PT05OvnfvHl7LWlpaEEIcDufYsWMzZ84c\nMWIEMY6Pj4+BgQGLxcrOzo6Pj8dzyWX2VIqcnBxnZ2eJRmIRY6cgpKknQIEGQAbpfDiZMXjt7e1H\njhy5dOnS8OHDFy9erKWlVVNTgzqI1mMymRMnTjQ2Ni4vL6fT6fiCQJk9CZCiN9Ap9zZs0urllSOw\nUIVUsrKyGAzGr7/+yuPxJk2apKKiwuVy8bfwaxg0Gq2urg7DsJs3bzo7O+NvCQQCDQ0NPp9PjOPo\n6Jienk683LhxI0KIQqHo6elsShjqAAAbIElEQVTha2Q66on78OEDnU6/efNmU1PTtm3b8B2lpaWN\nGDGioqLi9evXenp6eXl50i0cDmfIkCHnzp1ramrasGGDo6MjhmHV1dUIoZCQED6f7+Hhoaury2Kx\nuFyujo5OQUEBk8lECO3fv7++vn7//v0WFhYYhr1//15HR6e9vV3mgHLOEuguA+4iITlVVFT0ZlAc\n6AiFQjl8+LCcfDiJGLzMzExivRKbzTYyMiIy8ySi9aKioh4/fvzq1avBgwfPnDnz7NmzxKMGZIbw\nIUjRAwPwLg5yev78uRCemU0Ourq60vlwHcXglZeXe3j854mRqamp4nO44tF6IpFox44dGRkZFhYW\nCCEvLy/xB3HJDOFDkKIHEExxdADmHAasmpoa/BGFGIaJRCJHR8eQkBA2m01kaGAY5uXlFRUVhWGY\nvb39nTt3MAwTCoWOjo6PHj0ixomOjl60aBH+M5fLVVVVJWY/1q5du3//fpk9CXw+39TUNC0tTSgU\nikSi0aNHs1ispqamrKwsDMMEAsHp06etrKykWzAM09PTY7PZ+DhsNru+vh7DME9Pz7i4OLynjo5O\nS0sLhmHXr1+fM2cOhmFWVlZ3797FNzly5MiqVaswDPPw8EhISJA5oMyz1NVTD6TARUIA/guRoldf\nXx8cHMxiscaNG0fE4DU0NJw9ezYnJ2fp0qUIoaFDh8bHx797927Hjh1UKnXChAnEOOLReoaGhgYG\nBiEhIQ0NDfHx8bdv3xa/d00ihC8jI+P58+dlZWVEit7OnTvxFL2ioiIXF5f8/Pw///yTSNGTaEF/\nh97V1dVdu3Zt7Nixf/zxB0Lo2bNn+CyKnBS9pqamK1euHD16dO/evUjsCqH0gDLPUs//cgYeZf8X\ngqTgG/RAJjMfTjoGD8OwkpISOzs7Go22cOHCxsZG8UEkovXS09Otra11dXXd3NyIb6Mye0KKHiAM\nuDQ7BcFKQgCA0sEUBwAAkBQUaAAAICko0AAAQFJQoAEAgKSgQINPJD/jraWlxdfXd8iQIV988UVc\nXBy+SX19vZqaGuVvx44d62hzHGS8gQEOCjT4RHIy3jAMW7BggZGRUVFR0aZNm1auXInfinvt2rUl\nS5aw/4ZHbsqMiIOMNwAQFGjwURTMeGMymVwu9+DBg3Q6PSAgQENDo7a2FiGUlJTk6+tr/jdtbW2Z\nmyPIeAMAp9S7sMkLFqpIUzzjLTQ0dPny5fhblZWV6urqbW1tHA5HXV3dxcVl0KBB48aNe/HiBTGy\nxOYEyHgDAxyEJfUZf/zxxy+//NJpt5iYmIaGBjW1Dn+zurq6+vr6H7VrFRWV2NhYxTPe6HR6Wlpa\nYWEhlUpdsWKFpaWlurp6Tk7O9OnTDx48aGBgEBAQsH37duLLr8TmeCNkvAEABbrPyMnJCYmIHjv5\nO/ndhjlPGdbxuy3NzR8qXkhf9ZKPQqGoq6srnvHm6el569at8ePHOzg4DBs2DJ+1mDt3LnG5bPny\n5eHh4R1tjoOMNwCgQPclFl/bTF20tCsjNNTVPogIka56naqtrRUIBCYmJgghDMOuXr3q6enJ4XDc\n3NwaGxvxApebm4t/c6ysrIyLi1NVVRWJRE5OTnv37sVDgn766Se8Z35+/jfffNPR5jgWiyU9Af3X\nX39Nmzbt/Pnzrq6uFArF1tZ23Lhxzc3NZWVl4eHhJ06cOHv27IYNG5ydnSVaZsyYUV5eXlhYiD/H\npKKiQldXFz8Sd3d3hFB7e3tRUZGNjQ1CKC8vD//mnpubu2TJEnzXSUlJU6ZMwRvxqi09IJfLlT5L\nH3uqASDARcI+7MCqpdWvXxEvywufb/SYtszB6tTWwDY+X07jJ/iojDd7e/vY2FgejxccHMxgMKZN\nm6aiorJ+/fqYmJjGxsbExMSkpKT169d3tDkOMt4AgIuEspHwImF6evr364KvvXh37cW7PeevfLfQ\nCyF04u4jvOVKUSX9c1P/fUejHuaN+maCV+AWmY1Rj555eXl92gEonvF24cIFBoNhaGjo7+9PhJyd\nP3/ezMyMRqPNnj27srISb5S5OQ4y3gCANDvZSJhmd//+/ehbqQt/CkQI3Tx3hltd9eu1xJCUVJMv\nLRFCBVlPovZtO3H3EUKoKOfpmZ3Bp1IzpRv/dSnlQUQIsXIEAEBmMAfdJ81b4Y8Qyk6/R7TUVb8d\nOsIa/9lshDX3XRWGYTIbe/9oAQCfBgo0qe3atSsvLw+fGK2rqzO2ceioZ/P7Bk1tbfznQTo67QIB\nv+VPGY2trb1w2ACAbgEFmtTa2trWrl3r4OCAEHr48GHSg6cd9dTR1yeKb+uff6qqqmkO0pbRqKXV\nC4cNAOgWcBcHqamqqurq6hoYGBgYGOC3hXWEYTr0bdl/biKuKi+jm5hSVFRkNvb4QQMAugn8c+0n\nvnL8pqW5MSv1Z35ry82Y8G/nLuiosf+RTow7fvw4fnczQSIYr7W1de3atcbGxjQazc/PTygUog5y\n9VxcXChi9PT0xPP2ehPk6g1AMMWhHFVVVe7u7jQaTX638vLytLQ0fGV2Q0PDiAlTOuqpqqq2JfxC\n+I6giF2bxrm6ua9aK7Pxj8b33ftByEAikOj27dtBQUF79uzBX3I4nPj4+JiYmK+//proExQU9Pbt\n2+zsbB6PN3Xq1FmzZs2ZM8fLy8vd3T0yMrK6unrs2LGbN282MzO7dOkSHrSEENq0aZOTk5OKkv4K\nycnJCQgIkGjMzc3Fl890CmKb+iIo0MrB4/EY/7Bdse1fim9SkPXkxbNc8Zaoh3niLy2+tjl2I11i\nK5mN/cD9+/e3b99eWlq6bt26p0+f7tq1C2/Pz8/fvXu3lZWVnZ0d3iIdjNfU1HT58uVXr14xGAwz\nMzMHBwcej4fn6s2dO7etrU08V8/U1BTf6urVq4aGhoGBgeKHkZmZuWHDhsLCwpEjR27btg1fN3jq\n1KmjR49yudwFCxacO3dOTU1NuqWsrGz9+vWZmZkWFhbnz58fPXo0hmEmJibbtm07fvy4QCC4fPny\nuXPnfv7551GjRt2/f7+trQ3P1Zs9e7aqquqFCxcmTpyIEMrNzd2yZQtCSHpAOWcJ9BVQoPsMDS2t\njBtXi39jdmUQQVubi4OMAM++hclkent7JyQk2NjYuLu7Z2dn42sO3717t3z58qSkJHt7e6JAu7q6\nurq6Ojk5EesSNTU1X7x4wWAwEELFxcVMJjMqKkpbW3vevHmRkZGrV6+m0WilpaUUCoXYY2lp6YkT\nJ9LT08Ub+Xy+u7t7dHT05MmTDx8+HBoa6unpmZ6efvLkyYcPH4pEIltb24KCAh6PJ9FibGw8YcKE\nQ4cOJSQk7N27d9WqVUwmk8PhcDgcPp9fXFzs5eU1e/bsjIyMo0ePDhs2rLi4uKWlhcfjaWhosFis\niIgIX1/fV69eNTY2vnv3ztrauqamRnrAjs4S6EOgQJPFgVVLfbbswVedIITKC5+f2RVcV1XpMGXG\nqt3/q6GpOdJ27IZ/nxFvVFFVjQn6MTU1VXwcNze31tZWTU1N2bsZRC0tLZ06dWpPf5weoqamlpyc\nLDMxrrW11dPTMzQ0lEKhDBo0CA/EwEkE41GpVBMTk/b2djyi7+LFi0Rnmbl6ra2tS5cujY6OplKp\n4gcDuXqgp0GBVr6CrCdP7t7Ie/Srz5Y9eItQ2B6y7kdP//VjJrme2LLuduzZBX7rpBvx5SoSbt26\n1dqvb3bW0tKSTozDMOyHH37w9fV1cXFJTEy0s7MT/6orHYyXmZnp7+9vY2PDZDJNTU3l5+rt2rVr\nypQp+KSBOMjVAz1OeavMSa2nszjy8/NnLVuJx2gs27RzhrcvVVOLCNbYc/6KyZeW+M//uphsPHSY\nzMakgjfTpk3r0eMkoZqaGhUVldbWVgzDRCKRo6NjSEhIeXk5QojBYDAYDB0dHU1NzRkzZhCbREdH\nL1q0iHh5+fJla2vrp0+fEi1sNpsI5cAwzMvLKyoqCv+5qqpq6NChTU1NEofB5/NNTU3T0tKEQqFI\nJBo9ejSLxWpqasrKysIwTCAQnD592srKSroFwzA9PT02m03sur6+HsMwT0/PuLg4vKeOjk5LSwuG\nYdevX58zZw6GYVZWVnfv3sU3OXLkyKpVqzAM8/DwSEhIkDmgzLPUxTMPeh/cZtf9uFxuTk4OS67i\n4mKi/7wV/it3HtAWC9GHddtyyEyMMzU15XA4+fn5+ERBZGSkeN6IeDBeTU1NYGBgdHS0sbFxRUVF\nRUXFhw8f5OTqHThwYM2aNXp6esRokKsHeg1McXS/w4cPV1ZWfvnll3L61NXVIaTa0buKrttuaenG\nw+4r1NXVIyIiAgMDT58+vWbNGk1NzTFjxmhoaBgZGSGEBAJBQUHB9OnTxSc0fvvtN2J+ICkpicvl\nij8Htri42NraOioqavPmzTt27Jg0adLDhw/xSfyKior4+PjXr1+LH8CPP/64adMmHx+fb7/91tzc\nfMSIETt27Bg+fPiRI0e2bt26ZMmSiRMnqqurT506NTw83NDQUKIFIRQREbF06dIzZ87Y2trevn3b\n3Ny8ubm5rq4Of0IuXpTxfeXl5fn5+bFYLE9PT5FIZGJiYmNj8/PPPxsZGdXW1v7555/4/82kB8Qb\nJc5SD/5WQM+ANDvZupJmd/ToUSsrq9mzZ8vp8/z58+3/PiV+m92P347BJzEQQmlXLv+e9SToeARC\n6I/G9/90tkkqqEhPjpdojMsrOxfsJ3GREADQb8AUBxnBum0AAIICTU4Ded02AIAABZqM8CXayWdC\n/VzG6egbiK/blmgEAPRjMActm8w56NbW1qCgoE6zcn7//XcajWZmZtZRB2traxcXF7/grdOX+HTl\nIIXC9pyrsTAHDUB/BXdxfAQej5f3qsJn82753WwXdjJO4r7NPj4+7lO+RXVs+T1TU1OtrKyGDh3a\nUYf/+Z//6WRnvSUlJeXChQu3bt0iWo4fP/7o0aPr168TLRUVFfPnz8/Pz8dfvnnzJjg4+MGDBwKB\nwNvb++TJkyoqKu/fvw8ODr5169Zff/21ZMmSM2fOUCiUlpaWNWvW3Llzh0qlHjp0yNvbu7c/3t/8\n/f1tbW39/PzEG8PCwl6+fInfoSGf9FkCQA4o0B9HS1uH8UWH5VJBKioqNBpt8+bNnfZsbm6eNWvW\n+PHju7jHXvCxkXIYhslMj/Px8RkxYgSLxWppaXFwcFi3bp21tfWCBQvs7OyKiooSExNXrlw5d+5c\n+enYPQci5UBvgjnoLjmwamn161fEy/LC5xs9pi1zsDq1NbCNz5fT2D/cv3/fycnJ0NBwz5494qVH\nkUg5PD0uKCho0KBB4ulxTCZz4sSJxsbG5eXldDrd3NycyWRyudyDBw/S6fSAgAANDY3a2lrxw8jM\nzHR0dNTW1h4zZkxycjLeeOrUKXNzc21t7WXLluF5odItZWVls2bNotFoY8eO/f333xFCGIZ9/vnn\np06dsrCwMDMze/TokY+Pz+DBg11cXEQiEZ/PxyPlzM3NLSwsHj9+jO+LWDAiPaCcswRA55S5jJHE\nZC71rqys/HaeJ7Ea+7uFXgghYn32laJK+uem/vuORj3MG/XNBK/ALR01Kr6OfNu2bZmZmd36ybpN\nVlYWg8H49ddfeTzepEmTVFRUuFwuhmHV1dWjR48uKSnR1tZ++/at+CaOjo7p6eniLREREQghGo1W\nV1eHt2zcuBEhhEfjV1dXYxgWGhq6fPly/N3Kykp1dfW2tjZihA8fPtDp9Js3bzY1NW3bts3Z2RnD\nsLS0tBEjRlRUVLx+/VpPTy8vL0+6hcPhDBky5Ny5c01NTRs2bHB0dMQPHiEUEhLC5/M9PDx0dXVZ\nLBaXy9XR0SkoKGAymQih/fv319fX79+/38LCAsOw9+/f6+jotLe3yxywo7MEgCJgiuMTvS4pUNeg\nUjX/7xF/xTlMDU1NvGov9F9/ZmfwAr91MhsV30tFRUVgYKD4OmOSmDdv3u3bt7sSKYeTSI+Liop6\n/Pjxq1evBg8ePHPmzLNnz+7Zs4dOp6elpRUWFlKp1BUrVlhaWqqrqxMjQKQc6MegQH8iPEkuO/0e\n0dITARrnz59vIeV6bg0NjY0bN35ypJzM9DiRSLRjx46MjAwLCwuEkJeXV0ZGBkLI09Pz1q1b48eP\nd3BwGDZsGD4TQoBIOdCPQYGWDX8SkkQjn89X/cyoo00UDtD4U/HD0NDQ0NDQ+Mhj7w21tbUCgQD/\ngoxh2NWrVz09PdlsdkpKSmZm5s6dO1taWtrb22fOnHn37l18ExaLRUxAczgcNze3xsZGvBTm5ua6\nuLg0NDTweDwiw+Tly5f4UwQrKyvj4uJUVVVFIpGTkxOeE4T766+/pk2bdv78eVdXVwqFYmtrO27c\nuObm5rKysvDw8BMnTpw9e3bDhg3Ozs4SLTNmzCgvLy8sLMRjKyoqKvCrjvn5+fjDDNvb24uKimxs\nbBBCeXl5+Bf/3NzcJUuW4LtOSkrCLwwSVVt6QC6XK32WevT3AvoZuEgo2++//54mJT4+Xs4mOvr6\n/L+DmFv//FNVVU1zkLbMxh4/+p7XxUg5melxhoaGBgYGISEhDQ0N8fHxt2/fxh+Ham9vHxsby+Px\ngoODGQzGtGnTEETKgYEBCnS3GVABGnik3M6dO11dXUeMGCEeKWdkZPTZZ5/JjJQjvkFTqdSoqKgD\nBw5YWlrevXsXT49TUVFJTEyMj483Nze/ePHigwcP8M3DwsK2bds2YsSIDx8+xMfHq6qqIoR+/PHH\nnJyckSNH4pFy8+fPHzt2LB4pZ29vjwfIffnllzweLzw8XLoFIRQREfHgwYOhQ4cePXq000i5MWPG\n4JFyCCETE5OwsDCZkXISA8o8S736ewJ9HKwk/Ahv3779Ye2GdYdPEC3iEXRCYftql3Ertu+zm+gS\nujHgy3+MWrQmSGbjydVenxyVBwAYOPrDtzmSgAANAED3gouEXRL1ME/8pcXXNsdupEv0kdkIAACd\nggL9EfT19akfms+sWSa/25s3b7S1teXc7tonlm4DAJQO5qC7nyJPVAEAgE7BHDRQppSUlLlz54q3\nHD9+HL8TmXhpaWnJYDDOnTuHt7x582bRokV0Op1GowUEBODpr62trWvXrjU2NqbRaH5+fkKhECH0\n/v37lStXDhkyRF9ff/Xq1Ur8LuLv7x8ZGSnRGBYWht9H2CnpswQGCCjQQJlkZuAREUvh4eHJycn3\n7t3Da1lLSwuGYV5eXg4ODqWlpU+ePImOjq6qqkIIBQUFsdns7OzsBw8eXLt27c6dOwghHx8fAwMD\nFouVnZ0dHx9fUlKilM+IEMrJyRF/TC2OuIG6UxCxNGBBgQa9TcEMvPb29iNHjly6dGn48OGLFy/W\n0tKqqamRmYHX1NR0+fLlmJgYMzMzOzs7BwcHHo+HZAXjiR8GZOCBPkApEU39W0hIyO3bt5V9FCSl\neAbezZs38Wg6DMMEAoGGhgafz8dfSmTg8fn8qqoq/K2ioiIDAwP8pXQwHgEy8ECfAHdxdL+3b98m\nJCSEhYUp+0DIhUajXblyRWa6m8wMvMzMTDc3N3xbNpttZGREpVLxlxIZeFQq1cTEpL29PTY2dvfu\n3RcvXjQxMZEZjEccDGTggT4BCnT3g9LcEZnpblgHGXjl5eUeHh54z9TUVGdnZ5kZeHiHzMxMf39/\nGxsbJpNpamraUTAeATLwQN+g3C/wYECpqalRUVFpbW3FMEwkEjk6OoaEhJSXlyOEGAwGg8HQ0dHR\n1NScMWMGhmH29vZ37tzBMEwoFDo6Oj569IjNZuvq6gqFQnw0Ly+vqKgoDMMuX75sbW399OlTYkdc\nLldVVZWYElm7du3+/fuJd/l8vqmpaVpamlAoFIlEo0ePZrFYTU1NWVlZGIYJBILTp09bWVlJt2AY\npqenx2az8XHYbHZ9fT2GYZ6ennFxcXhPHR0d/GLm9evX58yZg2GYlZXV3bt38U2OHDmyatUqDMM8\nPDwSEhJkDijzLHX3rwL0DVCgQe9pa2vT0dGJiIjgcrlBQUFqamoZGRl//fUX528ODg6xsbENDQ0Y\nhi1YsMDb27u6unrr1q2TJk3Cr7ZpamqePn2ax+NFRkZaWlp++PCBw+HQ6fTMzEz231pbW4VC4eDB\ng/ft28fj8eLi4szNzfExHzx4kJ+fX1BQoKurW19fX1lZic88fPjw4enTp5qams+ePXv//r2/v7+H\nh4d0C4ZhEyZMCAwMrK2tTU5ONjQ0xGurhYVFSUkJhmG///77V199hX/YnTt37tmzp7GxUUtLy9vb\nu7GxMSkpaciQIRwOB8OwL7744tWrVzIHlHmWlPMLA8oGBRr0qsuXL9Pp9FGjRkVGRuro6DQ3NxNv\ntbW1aWlp1dbW4i9LSkrs7OxoNNrChQsbGxvxxkuXLn3++ecGBgbz5s0jHogl8UdhcXExhmHp6enW\n1ta6urpubm7EV1RLS8uzZ8+2tbXNnj1bR0dnzJgxKSkpw4cP37t3b1tbm6+vr46OjoGBwaJFi2pq\naqRbMAwrLCy0tbXV1NR0cnLCn0bW1NREfK+/cOHCsmXL8H3NmjXr1q1b9+/f/+GHH7y9vbW1tceP\nH5+Tk4NhWE1NjYGBgUgkkjmg/LMEBhRYSQgAACQF90EDAABJQYEGAACSggINAAAkBQUaAABICgo0\nAACQFBRoAAAgKSjQAABAUlCgAQCApKBAAwAASUGBBgAAkoICDQAAJAUFGgAASAoKNAAAkBQUaAAA\nICko0AAAQFJQoAEAgKSgQAMAAElBgQYAAJKCAg0AACQFBRoAAEgKCjQAAJAUFGgAACApKNAAAEBS\nUKABAICkoEADAABJQYEGAACSggINAAAkBQUaAABICgo0AACQFBRoAAAgKSjQAABAUlCgAQCApKBA\nAwAASUGBBgAAkoICDQAAJAUFGgAASAoKNAAAkBQUaAAAICko0AAAQFJQoAEAgKSgQAMAAElBgQYA\nAJKCAg0AACQFBRoAAEgKCjQAAJAUFGgAACApKNAAAEBSUKABAICkoEADAABJQYEGAACSggINAAAk\nBQUaAABICgo0AACQFBRoAAAgKSjQAABAUlCgAQCApKBAAwAASUGBBgAAkoICDQAAJAUFGgAASAoK\nNAAAkBQUaAAAICko0AAAQFJQoAEAgKSgQAMAAElBgQYAAJKCAg0AACQFBRoAAEgKCjQAAJAUFGgA\nACApKNAAAEBSUKABAICkoEADAABJQYEGAACSggINAAAkBQUaAABICgo0AACQFBRoAAAgKSjQAABA\nUlCgAQCApKBAAwAASf1/bJDiBTRNTUoAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Other ways to improve assembly"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you look at the stats output for step 5 you can see that we discarded most of our data for being too low of coverage. One way to salvage these data is to set the minimum depth for majority rule consensus base calls to 2 (param 31). In this way, statistical base calls are still made for all stacks in which the depth is >mindepth (param 8), but for everything between (2-mindepth) majority rule is used to make base calls. This will probably increase the size of this data set ~10X, as you can see in the example S5 stats file below.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%bash\n",
      "## This is the step 5 results after I moved the consens.gz files out of \n",
      "## the clust.85/ directory and replaced them with new consens calls \n",
      "## generated by running step 5 with the majority rule consensus \n",
      "## option turned on. \n",
      "\n",
      "tail -n 32 analysis_pyrad/stats/s5.consens.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "taxon          \tnloci\tf1loci\tf2loci\tnsites\tnpoly\tpoly\n",
        "d35371.assembled\t60885\t24564\t23062\t1780882\t2350\t0.0013196\n",
        "d41058.assembled\t329719\t90853\t84689\t6244501\t3656\t0.0005855\n",
        "d41237.assembled\t293759\t90845\t84120\t5741117\t2794\t0.0004867\n",
        "d33291.assembled\t169929\t67772\t62376\t5502056\t4098\t0.0007448\n",
        "d34041.assembled\t197867\t89535\t83219\t5919976\t6740\t0.0011385\n",
        "d35422.assembled\t169511\t56563\t50274\t3627878\t3974\t0.0010954\n",
        "d35320.assembled\t156419\t52912\t48780\t4266120\t3268\t0.000766\n",
        "d39187.assembled\t96755\t45614\t41058\t3775082\t4079\t0.0010805\n",
        "d39104.assembled\t101266\t34065\t31273\t2410273\t1958\t0.0008124\n",
        "d19long1.assembled\t161936\t39462\t35591\t3167415\t2402\t0.0007583\n",
        "d41389.assembled\t211387\t77211\t70021\t5475971\t6483\t0.0011839\n",
        "d40328.assembled\t136690\t65435\t59329\t5145197\t7386\t0.0014355\n",
        "d41732.assembled\t452240\t181118\t170999\t12626514\t7201\t0.0005703\n",
        "d39114.assembled\t218883\t92037\t85134\t7595069\t7721\t0.0010166\n",
        "d39404.assembled\t243761\t89814\t78903\t6913918\t7306\t0.0010567\n",
        "d39531.assembled\t332792\t112016\t101057\t8879145\t7975\t0.0008982\n",
        "d35178.assembled\t450842\t201307\t190564\t13878623\t21667\t0.0015612\n",
        "d31733.assembled\t439816\t175627\t160987\t13694620\t12159\t0.0008879\n",
        "decor21.assembled\t577428\t227921\t207326\t15697724\t14727\t0.0009382\n",
        "d30181.assembled\t516149\t203114\t189290\t15726272\t13709\t0.0008717\n",
        "d39103.assembled\t642379\t235667\t211348\t18531562\t13982\t0.0007545\n",
        "d30695.assembled\t485553\t174196\t155764\t14331377\t15685\t0.0010945\n",
        "d39253.assembled\t431026\t172635\t151210\t13507456\t15403\t0.0011403\n",
        "d39968.assembled\t456174\t186053\t162430\t14528464\t19145\t0.0013178\n",
        "\n",
        "    ## nloci = number of loci\n",
        "    ## f1loci = number of loci with >N depth coverage\n",
        "    ## f2loci = number of loci with >N depth and passed paralog filter\n",
        "    ## nsites = number of sites across f loci\n",
        "    ## npoly = number of polymorphic sites in nsites\n",
        "    ## poly = frequency of polymorphic sites\n"
       ]
      }
     ],
     "prompt_number": 70
    }
   ],
   "metadata": {}
  }
 ]
}