{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "### Intro to downstream ipyrad data visualization\n", "All of the modules below are included in the ipyrad conda installation. " ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ipyrad 0.4.1\n", "numpy 1.11.1\n", "h5py 2.6.0\n", "toyplot 0.13.0\n" ] } ], "source": [ "## import libraries\n", "import ipyrad as ip ## ipyrad\n", "import numpy as np ## array operations\n", "import h5py ## access hdf5 database file\n", "import toyplot ## my fav new plotting library\n", "import toyplot.html ## toypot sublib for saving html plots\n", "\n", "## print versions for posterity\n", "print 'ipyrad', ip.__version__\n", "print 'numpy', np.__version__\n", "print 'h5py', h5py.__version__\n", "print 'toyplot', toyplot.__version__" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Import the ipyrad Assembly class object\n", "I assembled the data set under three minimum cluster depth settings (6, 10, 20), \n", "and import their Assembly objects as data1, data2, and data3. These are ipyrad Assembly class objects which have many features and functions available to them. To see these type the object name (e.g., data1) followed by a period (.) and then press tab to see list of all the available options. \n" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " loading Assembly: pedictrim5\n", " from saved path: ~/Downloads/pedicularis/pedictrim5.json\n" ] } ], "source": [ "## import Assembly objects\n", "data = ip.load_json(\"/home/deren/Downloads/pedicularis/pedictrim5.json\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Open database files\n", "\n", "We also open the database file for each data set. This is the file with the suffix \".hdf5\" that should have a file name like: \"[project_dir]/[name]_consens/[name].hdf5\". The assembly class objects save the database file path under the attribute \".database\", which can be used to access it more easily. Below we open a view to the hdf5 database file for each of the three assemblies. " ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "## load the hdf5 database \n", "io5 = h5py.File(data.database, 'r')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### What's in the database file?\n", "The hdf5 data base is compressed and sometimes quite large. If you moved your JSON file from a remote machine (e.g., HPC cluster) to a local machine you will have to update the `data.database` path to the location of the database file on your local machine. " ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "location of my database file:\n", " /home/deren/Downloads/pedicularis/pedictrim5_outfiles/pedictrim5.hdf5\n", "keys in the hdf5 database\n", " [u'edges', u'filters', u'snps']\n" ] } ], "source": [ "print 'location of my database file:\\n ', data.database\n", "print 'keys in the hdf5 database\\n ', io5.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Access arrays\n", "The hdf5 data base contains the following five arrays with the following dimensions." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "\n" ] } ], "source": [ "## This doesn't actually load them into memory, they can be very large.\n", "## It just makes a reference for calling keys more easily\n", "\n", "#hcatg = io5[\"catgs\"] ## depth information (not edge filtered)\n", "#hseqs = io5[\"seqs\"] ## sequence data (not edge filtered)\n", "hsnps = io5[\"snps\"] ## snp locations (edge filtered)\n", "hfilt = io5[\"filters\"] ## locus filters\n", "hedge = io5[\"edges\"] ## edge filters\n", "\n", "## arrays shapes and dtypes\n", "#print hcatg\n", "#print hseqs\n", "print hsnps\n", "print hfilt\n", "print hedge" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A function to filter the SNPs array" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def filter_snps(data):\n", " ## get h5 database\n", " io5 = h5py.File(data.database, 'r')\n", " hsnps = io5[\"snps\"] ## snp locations\n", " hfilt = io5[\"filters\"] ## locus filters\n", " hedge = io5[\"edges\"] ## edge filters\n", "\n", " ## read in a local copy of the full snps and edge arrays\n", " snps = hsnps[:]\n", " edge = hedge[:]\n", "\n", " ## print status\n", " print \"prefilter {}\\nshape {} = (nloci, maxlen, [var,pis])\"\\\n", " .format(data.name, snps.shape)\n", " print \"total vars = {}\".format(snps[:,:,0].sum())\n", "\n", " ## apply edge filters to all loci in the snps array\n", " for loc in xrange(snps.shape[0]):\n", " a, b = edge[loc, :2]\n", " mask = np.invert([i in range(a, b) for i in np.arange(snps.shape[1])])\n", " snps[loc, mask, :] = 0\n", "\n", " ## get locus filter by summing across all filters\n", " locfilter = hfilt[:].sum(axis=1).astype(np.bool)\n", "\n", " ## apply filter to snps array\n", " fsnps = snps[~locfilter, ...]\n", "\n", " ## print new shape and sum\n", " print \"postfilter {}\\nshape {} = (nloci, maxlen, [var,pis])\"\\\n", " .format(data.name, fsnps.shape)\n", " print \"total vars = {}\".format(fsnps[:,:,0].sum())\n", " \n", " ## clean up big objects\n", " del snps\n", " del edge\n", " \n", " ## return what we want\n", " return fsnps" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def filter_snps(data):\n", " ## get h5 database\n", " io5 = h5py.File(data.database, 'r')\n", " hsnps = io5[\"snps\"] ## snp locations\n", " hfilt = io5[\"filters\"] ## locus filters\n", " #hedge = io5[\"edges\"] ## edge filters\n", "\n", " ## read in a local copy of the full snps and edge arrays\n", " snps = hsnps[:]\n", " #edge = hedge[:]\n", "\n", " ## get locus filter by summing across all filters\n", " locfilter = hfilt[:].sum(axis=1).astype(np.bool)\n", "\n", " ## apply filter to snps array\n", " fsnps = snps[~locfilter, ...]\n", "\n", " ## clean up big objects\n", " del snps\n", " \n", " ## return what we want\n", " return fsnps" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [1258, 940],\n", " [1379, 1061],\n", " [1456, 1008],\n", " [1474, 999],\n", " [1551, 1055],\n", " [1596, 1089],\n", " [1532, 1103],\n", " [1530, 1180],\n", " [1594, 1116],\n", " [1532, 1128],\n", " [1620, 1148],\n", " [1565, 1205],\n", " [1642, 1216],\n", " [1676, 1231],\n", " [1667, 1123],\n", " [1625, 1115],\n", " [1620, 1217],\n", " [1614, 1185],\n", " [1594, 1247],\n", " [1548, 1215],\n", " [1641, 1230],\n", " [1616, 1145],\n", " [1676, 1202],\n", " [1643, 1227],\n", " [1636, 1167],\n", " [1641, 1237],\n", " [1661, 1304],\n", " [1640, 1165],\n", " [1696, 1272],\n", " [1680, 1173],\n", " [1683, 1191],\n", " [1703, 1214],\n", " [1619, 1253],\n", " [1692, 1236],\n", " [1704, 1269],\n", " [1773, 1164],\n", " [1641, 1236],\n", " [1710, 1308],\n", " [1646, 1255],\n", " [1671, 1232],\n", " [1682, 1279],\n", " [1755, 1212],\n", " [1670, 1285],\n", " [1715, 1279],\n", " [1721, 1229],\n", " [1808, 1251],\n", " [1757, 1294],\n", " [1725, 1259],\n", " [1731, 1227],\n", " [1765, 1279],\n", " [1764, 1266],\n", " [1784, 1283],\n", " [1821, 1396],\n", " [1864, 1340],\n", " [1796, 1397],\n", " [1921, 1362],\n", " [1902, 1363],\n", " [1936, 1432],\n", " [1990, 1511],\n", " [2060, 1471],\n", " [2111, 1514],\n", " [2048, 1644],\n", " [2244, 1688],\n", " [2310, 1818],\n", " [2480, 1951],\n", " [2554, 2053],\n", " [2634, 2105],\n", " [2776, 2237],\n", " [2971, 2473],\n", " [ 747, 552],\n", " [ 482, 329],\n", " [ 319, 236],\n", " [ 205, 131],\n", " [ 106, 96],\n", " [ 61, 54],\n", " [ 32, 23],\n", " [ 9, 7],\n", " [ 2, 1],\n", " [ 1, 1],\n", " [ 1, 0],\n", " [ 1, 0],\n", " [ 0, 1],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0],\n", " [ 0, 0]])" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fsnps.sum(axis=0)\n" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [], "source": [ "## apply filter to each data set\n", "fsnps = filter_snps(data)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,\n", " 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,\n", " 44, 45, 46, 47, 48, 49])" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.arange(0, 40)#.reshape(10,4)\n", "b = np.arange(10, 50)#.reshape(10, 4)\n", "#np.concatenate((a,b), axis=1)\n", "np.array([a,b]).max(axis=0)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The snps array" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3 2 7 5 1]\n", "[[1606 1095]\n", " [1539 1106]\n", " [1543 1186]\n", " [1604 1119]\n", " [1539 1132]]\n" ] } ], "source": [ "## the last dimension has two columns (var, pis)\n", "\n", "## how many snps per locus\n", "varlocs = fsnps[:, :, 0].sum(axis=1)\n", "pislocs = fsnps[:, :, 1].sum(axis=1)\n", "print varlocs[:5]\n", "\n", "## how many snps per site (across loci)\n", "persite = fsnps[:, :, :].sum(axis=0)\n", "print persite[10:15]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Get a color map" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "colormap = toyplot.color.Palette()\n", "colormap" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/html": [ "
02468101214n variable (or pis) sites025005000750010000n nloci w/ n var (or pis) sites
    \n", "
  • \n", "
  • \n", " Save as .csv\n", "
  • \n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## deconstruct array into bins\n", "vbars, vbins = np.histogram(varlocs, bins=range(0, varlocs.max()+2))\n", "pbars, pbins = np.histogram(pislocs, bins=range(0, varlocs.max()+2))\n", "\n", "## setup canvas and axes\n", "canvas = toyplot.Canvas(width=350, height=300)\n", "axes = canvas.cartesian(xlabel=\"n variable (or pis) sites\",\n", " ylabel=\"n nloci w/ n var (or pis) sites\",\n", " gutter=50)\n", "\n", "## set up x axis\n", "axes.x.domain.max = 16\n", "axes.x.spine.show = False\n", "axes.x.ticks.labels.style = {\"baseline-shift\":\"10px\"}\n", "axes.x.ticks.locator = toyplot.locator.Explicit(\n", " range(0, 16, 2), \n", " map(str, range(0, 16, 2)))\n", "\n", "## set up y axis\n", "axes.y.ticks.show=True\n", "axes.y.label.style = {\"baseline-shift\":\"35px\"}\n", "axes.y.ticks.labels.style = {\"baseline-shift\":\"5px\"}\n", "axes.y.ticks.below = 0\n", "axes.y.ticks.above = 5\n", "axes.y.domain.min = 0\n", "axes.y.domain.max = 10000\n", "axes.y.ticks.locator = toyplot.locator.Explicit(\n", " range(0, 11000, 2500), \n", " map(str, range(0, 11000, 2500)))\n", "\n", "## add bars\n", "axes.bars(vbars, color=colormap[1], opacity=0.5)\n", "axes.bars(pbars, color=colormap[0], opacity=0.5)\n", "\n", "## or as a filled/smoothed plot\n", "#x = np.arange(0, len(pbars))\n", "#fill = axes.fill(x, vbars, color=colormap[0], opacity=0.5)\n", "#fill = axes.fill(x, pbars, color=colormap[1], opacity=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variation along the length of reads\n", "This data includes 75bp reads sequenced on an Illumina GAIIx. We know that the error rate increases along the length of reads, and that the error rate was a bit higher in this older type of data than it is in more recent sequencing technology. " ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0 0 0 0 0 1266 1384 1462 1485 1560 1606 1539 1543 1604 1539\n", " 1631 1575 1655 1685 1678 1633 1624 1626 1602 1559 1648 1630 1687 1656 1641\n", " 1651 1673 1646 1708 1694 1689 1712 1624 1696 1712 1778 1638 1719 1653 1675\n", " 1685 1757 1670 1724 1730 1809 1760 1732 1736 1764 1769 1794 1827 1865 1784\n", " 1910 1889 1912 1954 2032 2042 528 334 255 166 83 47 16 7 1\n", " 1 0]\n", "[ 0 0 0 0 0 941 1068 1008 1001 1060 1095 1106 1186 1119 1132\n", " 1151 1209 1221 1235 1128 1118 1225 1192 1252 1226 1241 1153 1206 1232 1172\n", " 1241 1306 1169 1277 1179 1192 1222 1262 1248 1277 1170 1243 1315 1265 1232\n", " 1279 1223 1289 1287 1232 1257 1297 1260 1232 1285 1269 1289 1397 1341 1407\n", " 1372 1358 1425 1492 1443 1485 433 251 179 87 66 38 21 8 2\n", " 0 1]\n", "maxlen = 76\n" ] } ], "source": [ "## the snps array is longer than the actual seq length (it's a bit padded)\n", "## and so we want to lop the extra on the end off. Let's get the max values w/ data.\n", "maxend = np.where(fsnps[:, :, :].sum(axis=0).sum(axis=1) != 0)[0].max()\n", "\n", "## all variables (including autapomorphies)\n", "distvar = np.sum(fsnps[:, :maxend+1, 0].astype(np.int), axis=0)\n", "print(distvar)\n", "\n", "## synapomorphies (just pis)\n", "distpis = fsnps[:, :maxend+1, 1].sum(axis=0)\n", "print(distpis)\n", "\n", "## how long is the longest seq\n", "print 'maxlen = ', maxend" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot distribution of variation along length of RAD loci" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def SNP_position_plot(distvar, distpis):\n", " ## set color theme\n", " colormap = toyplot.color.Palette()\n", "\n", " ## make a canvas\n", " canvas = toyplot.Canvas(width=800, height=300)\n", "\n", " ## make axes\n", " axes = canvas.cartesian(xlabel=\"Position along RAD loci\",\n", " ylabel=\"N variables sites\",\n", " gutter=65)\n", " ## x-axis\n", " axes.x.ticks.show = True\n", " axes.x.label.style = {\"baseline-shift\":\"-40px\", \"font-size\":\"16px\"}\n", " axes.x.ticks.labels.style = {\"baseline-shift\":\"-2.5px\", \"font-size\":\"12px\"}\n", " axes.x.ticks.below = 5\n", " axes.x.ticks.above = 0\n", " axes.x.domain.max = maxend\n", " axes.x.ticks.locator = toyplot.locator.Explicit(\n", " range(0, maxend, 5), \n", " map(str, range(0, maxend, 5)))\n", " \n", " ## y-axis\n", " axes.y.ticks.show=True\n", " axes.y.label.style = {\"baseline-shift\":\"40px\", \"font-size\":\"16px\"}\n", " axes.y.ticks.labels.style = {\"baseline-shift\":\"5px\", \"font-size\":\"12px\"}\n", " axes.y.ticks.below = 0\n", " axes.y.ticks.above = 5\n", "\n", " ## add fill plots\n", " x = np.arange(0, maxend+1)\n", " f1 = axes.fill(x, distvar, color=colormap[0], opacity=0.5, title=\"total variable sites\")\n", " f2 = axes.fill(x, distpis, color=colormap[1], opacity=0.5, title=\"parsimony informative sites\")\n", "\n", " ## add a horizontal dashed line at the median Nsnps per site\n", " axes.hlines(np.median(distvar), opacity=0.9, style={\"stroke-dasharray\":\"4, 4\"})\n", " axes.hlines(np.median(distpis), opacity=0.9, style={\"stroke-dasharray\":\"4, 4\"})\n", " \n", " return canvas, axes" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAEsCAYAAAA7Ldc6AAAJNmlDQ1BkZWZhdWx0X3JnYi5pY2MA\nAHiclZFnUJSHFobP933bCwvssnRYepMqZQHpvUmvogJL7yxLEbEhYgQiiog0RZCggAGjUiRWRLEQ\nFBSxoFkkCCgxGEVUUPLDOxPn3vHHfX49884755yZA0ARBQBARQFSUgV8Pxd7TkhoGAe+IZKXmW7n\n4+MJ3+X9KCAAAPdWfb/zXSjRMZk8AFgGgHxeOl8AgOQCgGaOIF0AgBwFAFZUUroAADkLACx+SGgY\nAHIDAFhxX30cAFhRX30eAFj8AD8HABQHQKLFfeNR3/h/9gIAKNvxBQmxMbkc/7RYQU4kP4aT6edi\nz3FzcOD48NNiE5Jjvjn4/yp/B0FMrgAAwCEtfRM/IS5ewPmfoUYGhobw7y/e+gICAAh78L//AwDf\n9NIaAbgLANi+f7OoaoDuXQBSj//NVI8CMAoBuu7wsvjZXzMcAAAeKMAAFkiDAqiAJuiCEZiBJdiC\nE7iDNwRAKGwAHsRDCvAhB/JhBxRBCeyDg1AD9dAELdAOp6EbzsMVuA634S6MwhMQwhS8gnl4D0sI\nghAROsJEpBFFRA3RQYwQLmKNOCGeiB8SikQgcUgqkoXkIzuREqQcqUEakBbkF+QccgW5iQwjj5AJ\nZBb5G/mEYigNZaHyqDqqj3JRO9QDDUDXo3FoBpqHFqJ70Sq0ET2JdqFX0NvoKCpEX6ELGGBUjI0p\nYboYF3PAvLEwLBbjY1uxYqwSa8TasV5sALuHCbE57COOgGPiODhdnCXOFReI4+EycFtxpbga3Alc\nF64fdw83gZvHfcHT8XJ4HbwF3g0fgo/D5+CL8JX4Znwn/hp+FD+Ff08gENgEDYIZwZUQSkgkbCaU\nEg4TOgiXCcOEScICkUiUJuoQrYjexEiigFhErCaeJF4ijhCniB9IVJIiyYjkTAojpZIKSJWkVtJF\n0ghpmrREFiWrkS3I3uRo8iZyGbmJ3Eu+Q54iL1HEKBoUK0oAJZGyg1JFaadco4xT3lKpVGWqOdWX\nmkDdTq2inqLeoE5QP9LEado0B1o4LYu2l3acdpn2iPaWTqer023pYXQBfS+9hX6V/oz+QYQpoifi\nJhItsk2kVqRLZETkNYPMUGPYMTYw8hiVjDOMO4w5UbKouqiDaKToVtFa0XOiY6ILYkwxQzFvsRSx\nUrFWsZtiM+JEcXVxJ/Fo8ULxY+JXxSeZGFOF6cDkMXcym5jXmFMsAkuD5cZKZJWwfmYNseYlxCWM\nJYIkciVqJS5ICNkYW53txk5ml7FPsx+wP0nKS9pJxkjukWyXHJFclJKVspWKkSqW6pAalfokzZF2\nkk6S3i/dLf1UBiejLeMrkyNzROaazJwsS9ZSlidbLHta9rEcKqct5ye3We6Y3KDcgryCvIt8uny1\n/FX5OQW2gq1CokKFwkWFWUWmorVigmKF4iXFlxwJjh0nmVPF6efMK8kpuSplKTUoDSktKWsoByoX\nKHcoP1WhqHBVYlUqVPpU5lUVVb1U81XbVB+rkdW4avFqh9QG1BbVNdSD1Xerd6vPaEhpuGnkabRp\njGvSNW00MzQbNe9rEbS4Wklah7XuaqPaJtrx2rXad3RQHVOdBJ3DOsOr8KvMV6Wualw1pkvTtdPN\n1m3TndBj63nqFeh1673WV9UP09+vP6D/xcDEINmgyeCJobihu2GBYa/h30baRjyjWqP7q+mrnVdv\nW92z+o2xjnGM8RHjhyZMEy+T3SZ9Jp9NzUz5pu2ms2aqZhFmdWZjXBbXh1vKvWGON7c332Z+3vyj\nhamFwOK0xV+WupZJlq2WM2s01sSsaVozaaVsFWnVYCW05lhHWB+1Ftoo2UTaNNo8t1WxjbZttp22\n07JLtDtp99rewJ5v32m/6GDhsMXhsiPm6OJY7DjkJO4U6FTj9MxZ2TnOuc153sXEZbPLZVe8q4fr\nftcxN3k3nluL27y7mfsW934Pmoe/R43Hc09tT75nrxfq5e51wGt8rdra1LXd3uDt5n3A+6mPhk+G\nz6++BF8f31rfF36Gfvl+A/5M/43+rf7vA+wDygKeBGoGZgX2BTGCwoNaghaDHYPLg4Uh+iFbQm6H\nyoQmhPaEEcOCwprDFtY5rTu4bircJLwo/MF6jfW5629ukNmQvOHCRsbGyI1nIvARwRGtEcuR3pGN\nkQtRblF1UfM8B94h3qto2+iK6NkYq5jymOlYq9jy2Jk4q7gDcbPxNvGV8XMJDgk1CW8SXRPrExeT\nvJOOJ60kByd3pJBSIlLOpYqnJqX2pymk5aYNp+ukF6ULMywyDmbM8z34zZlI5vrMHgFLkC4YzNLM\n2pU1kW2dXZv9ISco50yuWG5q7uAm7U17Nk3nOef9tBm3mbe5L18pf0f+xBa7LQ1bka1RW/u2qWwr\n3Da13WX7iR2UHUk7fiswKCgveLczeGdvoXzh9sLJXS672opEivhFY7std9f/gPsh4YehPav3VO/5\nUhxdfKvEoKSyZLmUV3rrR8Mfq35c2Ru7d6jMtOzIPsK+1H0P9tvsP1EuVp5XPnnA60BXBaeiuOLd\nwY0Hb1YaV9YfohzKOiSs8qzqqVat3le9XBNfM1prX9tRJ1e3p27xcPThkSO2R9rr5etL6j8dTTj6\nsMGloatRvbHyGOFY9rEXTUFNAz9xf2pplmkuaf58PPW48ITfif4Ws5aWVrnWsja0Latt9mT4ybs/\nO/7c067b3tDB7ig5BaeyTr38JeKXB6c9Tved4Z5pP6t2tq6T2VnchXRt6prvju8W9oT2DJ9zP9fX\na9nb+aver8fPK52vvSBxoewi5WLhxZVLeZcWLqdfnrsSd2Wyb2Pfk6shV+/3+/YPXfO4duO68/Wr\nA3YDl25Y3Th/0+LmuVvcW923TW93DZoMdv5m8lvnkOlQ1x2zOz13ze/2Dq8ZvjhiM3LlnuO96/fd\n7t8eXTs6/CDwwcOx8DHhw+iHM4+SH715nP146cn2cfx48VPRp5XP5J41/q71e4fQVHhhwnFi8Ln/\n8yeTvMlXf2T+sTxV+IL+onJacbplxmjm/Kzz7N2X615OvUp/tTRX9KfYn3WvNV+f/cv2r8H5kPmp\nN/w3K3+XvpV+e/yd8bu+BZ+FZ+9T3i8tFn+Q/nDiI/fjwKfgT9NLOcvE5arPWp97v3h8GV9JWVn5\nBy6ikLxSF1/9AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAHXRFWHRT\nb2Z0d2FyZQBHUEwgR2hvc3RzY3JpcHQgOS4xMJremEEAACAASURBVHic7d1/cNv3nd/5lwiSAn8I\nkPhDsiVBqrNyNhYZ79448UbMXa+TbSpp27ttmjWl67Yb2/rR63Q265XXM3ft2omT3c52FWn3dub+\nqK2kzXYydeSmaTvXSMole531iHQca38ohOjEtFQS4g/xlwAIIED84v1BfWGABECCBL7fD4jnY8Zj\nE18RePNjivy+8Pl83p9tS0tLSwIAAAAAGzQ4XQAAAACA+kEAAQAAAGAbAggAAAAA2xBAAAAAANiG\nAAIAAADANgQQAAAAALYhgAAAAACwDQEEAAAAgG0IIAAAAABsQwABAAAAYBsCCAAAAADbEEAAAAAA\n2IYAAgAAAMA2BBAAAAAAtiGAAAAAALANAQQAAACAbQggAAAAAGxDAAEAAABgGwIIAAAAANsQQAAA\nAADYhgACAAAAwDYEEAAAAAC2IYAAAAAAsA0BBAAAAIBtCCAAAAAAbEMAAQAAAGAbAggAAAAA2xBA\nAAAAANiGAAIAAADANgQQAAAAALYhgAAAAACwDQEEAAAAgG0IIAAAAABsQwABAAAAYBsCCAAAAADb\nEEAAAAAA2IYAAgAAAMA2BBAAAAAAtiGAAAAAALANAQQAAACAbQggAAAAAGxDAAEAAABgGwIIAAAA\nANsQQAAAAADYhgACAAAAwDYEEAAAAAC2IYAAAAAAsA0BBAAAAIBtCCAAAAAAbEMAAQAAAGAbAggA\nAAAA2xBAAAAAANiGAAIAAADANgQQAAAAALYhgAAAAACwTaPTBZRy7eoVXXr9Nfn9Q3mPHz12XKfP\nnFVPT69DlQEAAADYCGMDyODAgL70yst69Stf1ZG+Pnk8XklSIBDQ4OB1nXruWX3/Bz/MPg4AAADA\nfOYGkMHr6j9xQkePHc973Ofzyec7qVtDfvmH/DrS1+dQhQAAAADKZewekP0+nwKBQMFr4XBIg4PX\ntd/ns7kqAAAAAJth7AxIf/9JXbt6RZ96+hPq6e3J7vfw+4fkH/Lr9Jmz8hFAAAAAgJqybWlpacnp\nIkrx+4c0ODCgcDgkSerp6dXhnl7CBwAAAFCDjJ0BkeiCBQAAAGw1xgYQumABAAAAW4+5AYQuWAAA\nAMCWY2wA2e/zaXBgoOA1qwvWqTNni37+xQvnJUnnXnypKvUBAACYZDYeUTSV0EIqIUlaSCWUzKSV\nzKQVXIxJkv6nRw9pZ3OLk2UC5gYQumABAACsbSGV0M35cc3EImv+2VQmbUNFQGnGBhBJ+vo3vrmq\nC1Z//0kdfpUuWAAAABMLIf1kflyJNMECtcP4NrwbdfHCeV16/bVVj996730HqgEAAKis4eCURkIz\nZX3OU90HtLeVBj5wltEzIJttw3v6zFn2gAAAgC1lIZXQjdmx7L6OcoQSMQIIHGdsAKENLwAAQL7R\nyLzeC06x5Ao1zdwAQhteAAAASVIyk9ZwcEqjD+adLgXYtAanCyhmv8+nQCBQ8JrVhnc/G9EBAMAW\nlsykNRqZ159PjlQkfCTpggUDGDsDQhteAABQr2bjEU0shCo+42GdEQI4ydgAItGGFwAA1A9rtmP0\nwTxBAVua0QFEerjn48QJNpsDAIAtaWIhpNl4hP0dqBtGB5Avv/KyLl9+Q9LqlrqHP/Y4Z3oAAADH\n3Jwf1+iDeXW3tMvb3KLWxmZ5m1u0s7ml6OcEEzGFEjEtpBKajUc21Ep3M9ZzWjpQbcYGkMuX35Df\nP6S333lXknTquWd1Uec51wMAADhuJDyTnbGYiUVW3dh3t7SrtbFZTQ0uScvnb4QSMdrnAjK4C9bd\nQCB7/ofH49XFP/4TXbt6JbsXBAAAwAkTCyEN358q+WdmYstLqkZCMxoJzWgmFiF8AA8ZG0A8Hm9e\nG16fz6ejx47r4te+5mBVAACg1iUzaQUTG1v6FEzE9JP58QpXZC82uMNpxgaQ/hMndDcQ0DOf/1x2\nH8i5F1+S3z+kU89/weHqAABALbE6TN2YHdPVwC29NTmi4WDpWYxCz3FjZqzmZzIIIHCasXtAPB6v\n3vzOdzU4MJD3+Jvf+a4uX34jey4IAABAMVaHqcmF0KrgMBKaUSgR01NdB7J7NYpJZtJ6e/oON+9A\nBRgbQCxH+vpWPdbff9KBSgAAQC0IJmKaXAhpLDK/5mzFTCyiP58c0VPdB0p2rxoOTtnesQrYqowP\nIAAAAGuZWAgplIhpIhoqe5ZiIZXQW5MjemLXIzrk6V51fTg4taXO6AgmYupytztdBuoYAQQAgC1q\nODil1sZmHWzvcLqUiktm0pqJR4our9qI4ftTCiVierJjX3ZJ1mhkuZPVVpLM1PYeFtQ+AggAAFtM\nMpPWjdmx7NkUI6EZ7W3z6pCne829DiazDvGbjUc0Ea1OW35rBuXjHfskSe+VuVEdwNoIIAAAbCGz\n8YhuzOZ3alpIJTQSmtFYZF4H2jt0sL1DrY3NDla5PrPxSF7osKv7VHAxph9N35Gkmu94BZiIAAIA\ncMRIeEYfhGf0S7sfK7n5F+s3HJwquVwokU5nD8Y7uKNDhzzdjgaR2Xhk1WMz8YhCidiqk8XttpWD\nR2iDZ6AAlUIAAQDYyjrIzeoodGNmTH/z0UM1vTQo10IqYftN/UIqoRuzY2V1aRp9MK/RB/M6uKND\n3uYWdbvbq1L3bDyS3atBFykAEgEEAGCTZCatkfDMqnfoF1IJ3Zwf11NdBxyqrDKSmbRuzo9rIhrS\nwR0devLhHoJqm1gI6Sfz4xt+xz63u1NrY7P2tnmzgWQjodBqgWvCLAYAMxFAAABVN7EQ0vD9qaLt\nUSeiIXmbZwq2QK0FK0PA6IN5hRIxfWr3Y1Wd2VlryVW5rL0ilp3bW9Tlbld3TsvWYCJWsIvSQiph\n6z4NbBxLsOA0AggAoGoWUgkNB6fW1bFo+P6UutztG9oP4sSyJyl/1mOl4GJsXQfcFTIbj2hiYfV5\nFnbPKAQXYwouxrZcG9p6R0iE0wggAGCAhVRCC6lEVQ8HszoKSaraTIP1dVjvkq/nJOpc5e4HSWbS\nGg5OaXIhpAPtHXpi5yMbLb1s61n6tJBK6EfTd/SxnY+s6yyO2XhEI+EZli4B2NIIIADgsOHg1Kob\n9e6W5fX3rY3NampwaWdzixof/ns9rBAQevjPyhvalYetlSuZSWffoV9IJZYPhavATXM5+0FWbmYf\nCc1oNh7REzsfqWqQK7aXpZhEOq2bc+NaSCWKBiSCB+yWzKS3TOMH1B4CCOCAZCa9ag2ut7mFXwbr\nNBuPKJpKVK1rj11K7YtY741os8sl74pQEkrE1px1sA5be6rrQNljuNlNz2s+/zr2g4xG5vVecGpV\nDcHFmAbv3dEhb3dFD90LJmJKZdJKZNIl97KUMhKaUTKT1hM7H8nWRfCAU0KJWFWDOlAKAQSwSfDh\nQVqz8ci6bja6W5Z/MTza6l3X0o16UOhmzdok+2ir15GzJKxzDMoJkNa7/JW46UykNz7zEFyM6a2p\nkXWfw1HuO/+bUWo/yM358bzOTYVYh+59vGOf9rZ613w9q1Ws9OEG3WgysaGgUYq1Of3nPN0ai8wT\nPADUJQIIUCXJTDrb+34mFin7Rsa6MZmJRfRecKqmTi+utIVUQiPhmYI3nbmbZJtdLj3autxCdG+r\nt+IzStaSpoVUouCyJisMFWthaucN/Hol0mm9NTmiJzv3lQy6s/GI/vrhMiK7rNwPUu5ZF4l0Wjdm\nxjTW0q4nO/Zl/+5Ye2Gs/492n00RXIzpxsyYra8JACYxOoBcu3pFl15/TX7/UN7jR48d1+kzZ9XT\n01uR1zn327+V9/HFP/q/uM71DV//Vxcuajg4lb1Z/tbvX8i7/uv/4sW8j8u9/s//4Pe1t9WbnTo3\n7etf6/or/+oP1PZwX8Nan5/MpPVPfvOfaiH54U3vesfv5ty4dm5v0aVX/0BNrgZJkruhSb/3ta/l\nvateqD5rqU0oEdOF331FyXRGi+lUWa8vLQeSb/7eH8rd0KTUUkbhZFz/2/95bt2fb+f1b0lqbWrW\nN/7v1/Kuv/DCFxVJxRVJJGyvbyGV0Ol/9k/U6W5TLJ1UKBHb8Pg1NGxT0zaXfu3/+K2C16tRP9e5\nbvL1aCqhf2nY74d6u17PjA0ggwMD+tIrL+vVr3xVR/r65PEsT6EHAgENDl7Xqeee1fd/8MPs45vx\n67/+j7nO9YpcjyQX9fb0nbx3VPt+9XjJzy/3unV6cWtjsw55u3XiH/5DNW4r/k6/SeMzEp7RW5Mj\nkj7cu/A//uqvyLWtQe1NbrU3Lb9Dbc0UjEXm9Yu/8ssln7/U+AUXY3rq73027zHr9Vsbm9XW1Kxf\n/JVfVjydUjydzLtuefL45l5/M5/vxPW3p+/oqa4DampwaTYe0eFjfys7Nk7U99Tf+6zam7Yrklx0\n5PW5zvWten0hlTDq90M9Xq9n25aWlpacLqKQixfOS5LOvfhSwetffuVlHT12XEf6+jb0+UAhEwsh\nTS6EdLC9o+zNecU2xdpl5/YP9yDkbkruLvJ1RFP569u9zS3a2dxSlSVeyUxaN2bH1rXevdm1/DXQ\np9451lIyk5aKAaisQ95uW9tWA7mMnQHZ7/NpcGCg4LVwOLQ8C3LmrM1VYSsbjczr5ty4pOUuPN0r\n1o0XY51DsNam2GrLnXXJvdEv9yaytbFZO7e3yNvcsuFD4fLqWtEqdS0ED+dZ+2oAAKgGYwNIf/9J\nXbt6RZ96+hPq6e3J7vfw+4fkH/Lr9Jmz8vl8DleJUhZSCY1Glm/KK9kOsxBryc5GN2kPB6dW3ajP\nxCL64fhPdXBHhw55ugs+b7k317XAOtch92Tn7pbljdXljm8wEdOPpu8QKgDAMCtbwQN2MnYJlsXv\nH9LgwIDC4eWboZ6eXh3u6V0zfLAEyzmz8YgmFkJ5MwLNLte622GWayQ8ow/CM9mb3HKmlcuZvTi4\noyOvf7/TS66c0t3Svq7WwLkzSgAAs3S3tOtTux9zugzUKWNnQCw9Pb0V63aF6hqNzGtyIVRwnX+x\ndpibUawt6EhoRhPRkH6hc1/JfRzJTHrVhvFSRh8sf30HHt541+v6+JlYZM3WwIVmlAAAACTDA4hd\nbXixcdbSJ+tU5bXMxCJ6a2pEP+fpLnnKcSnrOcRtIZXQ4L072tvm1ZMd+1Yt/9ro0qlEOs2N9UPW\nWIyEZrKzIntbvbo5P563fAsAYJ5o0r4zfYCVjA0gdrbhxcasXPq0Xol0WsP3pzS5ENLHO/ate5Pz\nRg5xm4iGNBuP5AUe9iVUnjUrwpIrAKgNdh4qCqxkbgAZvK7+Eyd09Fh+D2ufzyef76RuDfnlH/IX\nbcOL6klm0hV5lzu4GNNbkyM65O1e1Sp2ZYtYSRqLzG8oNOQGni53+4afBwAAAJtnbAChDa+ZZuMR\n3Zgdq+gNvLWMp9poLQoAAOA8YwMIbXjNspHlTwAAwFzBRGzTZz0BG7Gl2/D++J139Mmnn857nLa8\n5duKZ10AAFDvjux5rGS3SKBajJ0BseS24Q2HQwoEAsx82GgkPKPh+1NOlwEAQF1zx6OKu9ucLgOo\nCGMDyOXLb0haXoolSV9+5eXsYx6PV69+5aurNqiv9Mmnn2bGo0zJTFqhREwz8Yhm4xFmPQAAcFBz\nclG7J2+rJRpU3N2u8K49Cu3cvaHn8dy/p/CuPUo0bZckJTI0ZIEzjA0gdwOB7H9fvvyG/P4hvf3O\nu/J4vPL7h3TuhS/mtefFxgQTMc3GI1pIJRRKsEkbAABTdE2PaVdOe3N3PCL3ZEQds+Oa79qnyI4O\npV2lb+W8wWntCM+pJRqUJC20ebMBJJSIaW8r91Gwn7EBJNfdQECnz5zNho2enl7t9/low7tByUxa\nw8EpjT6Yd7oUAACwQms0pD2Tt9WYjBe83piMa/fkB+q8N6rwrj2637k3L4hYsx2e+/fkyqTsKhtY\nt5oIINasR+6Sq3AoLI/X42BVtSmYiOnGzBgHEAEAYBhXOqVdcxN5sx4l/3wmpV1z4/Lcv6eop1Px\nlva82Y5CWqMhLbQx6wFnGRtA9vt8unb1ig5/7HF5PF6Fw6Hsfo5Tz39BHq8nuzkd6zMcnKKNLgAA\nBtrxYE7dE7c3NGPhyqTkCd6TJ3ivrM/jzUg4xdgA0t9/MrsBPRwOyT/kz147cuTT6j9xwqnSak4y\nk9aN2THNxCJOlwIAAHLkbjK3gyv9YcBJsgkdDjE2gOTyeLx5ez1OcwL6ulXj5HIAAOpBc3Ixu2G7\nGjrmJrRzdtzWfRpNyUXbXgsoxtgAEggE8jphFcIG9NJYcgUAqEebPTMjdxN43N2uqKdTEU9nxcKI\nOx5V9+RtuePOrkxgBgROMTaA3PIP6eKF8wqUCCG33nvfxopqRzKT1tvTd2ipCwCoO97gtHZPfqBU\nk1vzXfvKOjPDlU6pa3osby+FOx6ROx5R5/TopsNIuZvMqyF3qRf3CXCKsQHk6LHjOtzTq6Of/QxB\nowyEDwBAvbLCh/Rhq9qO2XHde/Qja3Z+8gan1XlvtORyqEJhJO5u06K7bc3zONZqrQvUE2MDiCT5\nfL41TzvHhwgfAIB6lRs+cjUm49o3dkuxtp2a79y7KohsdBO4FUZyxd3tWnI1Zpd/Wa/lDU6rPTxb\n1vNXU7X3tgBrMTqASNIf/fGfOF1CzRgOThE+AAB1p1j4yNUSDWpfNKhY207N7j6guLtt1Unjm2UF\nEivMOLnUqpTGRDwbQBZSCbU2NjtcEeqN8QEE63NzfpyTzQEAdWc94SNXSzQo352gUk1ulkOJAAJn\nNDhdADZvJDxD+AAA1J1yw0eueg4frdGQ0yWgzhFAatxoZF7D96ecLgMAHNOcXFTH3ETeAWvY+jYT\nPgA4iyVYNWxiIaSbhq4vBeqFNzitpkRcs7sPOF1KXco9yG3n7LjCu/bofufeNTsSobYRPjYnN6zP\nxCPqcrc7WA3qET+ha1QwEdNP5gkfgFNWds5JNrvLOm8Am1Ooc5Erk9KuuXHtmhtXeOce3e/aZ2Sn\nH1c6JW9wes2D8lLN7qrWv+PBnFojIXmC9xR3tyu8a48iOzqMDm/ueFSt0ZA6p0edLqWmcRo6nGbu\nTxkUFUzE9KPpO0qkOcEUsFuxg8Q6741q0d224dOX3fGokk3bjb75K4crnVLrQkgLrd6Kf007Hsyp\ne+J2yfMaPMF78gTvVTSIWMFhMydil3sWxP3OfRWd0WlOLspz/552hOfyanDHI3JPRrR78gNFPF1a\naPNWNIy0RkPqmJtQSzSY16o22exWsml7wXM0WqMhueNRNSXiakoult0mF4C5tsZvujpzY2aM8AE4\noNTNoyuTUvfkbU0ceKLsmzZ3PKq9o7eUbHZv6POrwZVObejraA/PqSUayrYjjbvbK/Y1udIp7Z66\nU9Z5ClYQiXi6FPF06MGOzrJftzm5qF2z49nTsTunR8sONhupXVpu4+q5f29TS8tc6ZTaH8zLc//e\nqnMrCmkPz6o9PJsNIxsdNyk/eFhWtqrNFWvbqaZEvK43iNshd+yTGe4nYD/nf8uhbAuphNMloA65\n41FlXI1GLmmptubkojqnx9a8eXTHI9o9dUeT+x5f93O70intHb0lVyYlVzyigyN/qYmDhzc8k7IZ\nhd4dL3aw2kKbV83JRbVEQ2qNhtQSCRackXDHI9o7NrzpELKeWY9SrJvq7obG5dOrW9rXfIe/0M2z\npZwZltx9KhthLS0rJ4hYS5W2x6ObOgAvd9xi7Tu16G7TQpt3ze/PUmNXCrMc9uOeAk4ggNSY2XW8\newVU0sp3f1NN7uUbkJb2spYcNScX1ZhYvqldeRKxybqmx+S5f2/dN4/t4Vl1uNs037l3zT/rSqe0\nd2w477ldmeVAYlcIWevd8UocrLaRYJZb30ZmDoo+XyaVDQ+5y41ibd5siPAGp9UxO76ud+E9wXtq\nC88VDAbueFRd02MVu6leGURyGx+0PgyC7nhUzbHIhsNOqde2wkinpHROkMsdu40GD9iP09DhJAKI\npBvvvpv38VOf+ISx14OJmLRvZ971O0O38j5+rPcw1w26PnbzJ2pKxNWQSatpMa4nDh1QU3JR29Ip\nzTz6EQ2PjJb8fKfqt/Y6jN94VwuZtKzbzic+sl+eYDwbSIb++6QWW9qVaG5Rstmtn3vio3LHo2pI\np+SOR/Wzn36g5kRMqZzPv9+5L3vzVOn6Dz3xUbUuhLQ9FtViS5tujt7b0PO741F1T97WnVvvKbfR\n9RMf2Z/354dv3837+ImP7Ffn9KjiD98pLjW+e8eGdefWe6s+PzeEVOv7ozUa0o7wnMb/4oYWJc2V\n8fVt7PpdRUfG5P2ff3nd9bdGQ5p/678pkP7wHdpK1xf4q7/K/neiuUW9B/bk3byv//n/u9INLu17\n6hMK79ojz/172jU3XsHx+/C6FUR2hOf0k9EpueMRhSWFK/T867luBbnxv7ghSUq5mpVxufT4/s51\nfT7Xnb/emIjrpz/9QPebW9S0a/kngMn3P1vxej0jgEj61rf+Xd7HK79BTLqeSKd09MX/Pe/6wH++\nkvfxyhsQrlf3+tv/6f9RQyatbUtL2pZO65O72+WORyUtv2v8X771X/L+/C89+79m/3vv6C392//4\nZ3nvmla6vh//h+9q++KCUo3NSrsadeiJj+a9XqHPz30H+Bv/7cd511f+gru24vrfzPn6JOkHP7i+\n6vN3zY2rKbmo6Uceq8j4u1JJbctk1JhO6u/8o+N51//031/V4vbW7Ne8nufvmh7LvtP/vT//i1X1\n5yp2/ZG7P9Poof+h4PNb4cMdjxT9fCuE/Ol/+L6STe6iX38549eUWtQnd7fnLbF6fYNf30avn/no\nY7r36EdK1p+70f+N/+9tW+v7hRXfv+V+/ot/49G8WaJK15d3A5mM68++/+dVe/6NXN/s+HHd3usD\n//mKGhq26a9aPJLMvv/Zitfr2balpaUlp4uohosXzkuSzr34ksOVVNZwcEojoRmny6hrXdNjFV3m\nkG5orNpym2K98gstO5HsXz4Rd7fr3v6Plr0MwHr3vjUaWvdm1fDO5SUrpdbOW7Me69moux5xd7sC\nj30877Hc8LFe04/+3IZa/Fr7AFqjIaOWxIR37skLIbnc8agevfszNiEDVZY7E/2/HPz4Gn8aqCxm\nQIB1KmddeDmqtea/1EFd1lpuafkmOerp3PRm1Y1wxyPaf/vmur/2zQSk3LX6Kw8NLNZad7Pc8Yi6\npseyr7eR8CEp+/9xrRDSGg3JlUmpNRIqK5zZzVq+tzKE5M48AQC2LgIIjGfdHGZcjeva2Ftpy4de\njVXsXfFCKh1Cyjkl2B2PVPVrW4srk5Lvzs2S7/JXamYmd+38fNc+hXbuLvtchnLtmhvXYkubFlq9\nGwoflt2TH6gpEddiS5u2x5aX+OUu9as1uSGk0jNPANaWexhhMpNWU4PLwWpQbwggNSaUiDldgq1W\n3pi0RkOa3X3Alu5Ale5gs5ZKhZBywodJrBvslZ19qrEkrDEZ1+7JD9Z9LsJmdU/cVrLZvenX2mqz\nA57gPTUlF6vStQlAaa70h3/nQomYutztDlaDekMAgbG8wWl13hvNuzFpiQbluxPMW7taaSvbztpp\nsyGkVsOHxdqcHvF0yHt/uurhz6533K0zPrBaLc7eAFvBtjShH84xOoBcu3pFl15/TX7/UN7jR48d\n1+kzZ9XT0+tQZagmVzq1fPZCiQBgLaO59+hHNnWmhCud0vZ4NK9trNM3RBsNIbUePiy5+1MAANXB\nkkc4ydgAMjgwoC+98rJe/cpXdaSvTx7P8k1mIBDQ4OB1nXruWX3/Bz/MPl4vZmJb+wdGOWvBG5Nx\n7Ru7ta7uRlbQcGVS2h5bDhxNibixm3TLDSFbJXwAAOwXZAkWbGZuABm8rv4TJ3T0WH5Pf5/PJ5/v\npG4N+eUf8utIX59DFaLSCi25Wg+ru9HcnoPKuFx5G3S3pVM1+y6PtTk73dCoRMvyL4bcMGLN/Ljj\nUXVOjxZ8DgAAinHHo4q725TMpJ0uBXXG2ACy3+fT4MBAwWvhcGh5FuTMWZurQjWsZ8nVms+RSW3Z\nGQBXJpVdFpa7PGyrbUgGANirgX0gcIixAaS//6SuXb2iTz39CfX09mT3e/j9Q/IP+XX6zFn5fD6H\nq7TXbA2+k28tfWpKLqopsbzcyWodSucbAACcw+9gOMXYACJJX//GN+X3D2lwYEDhcEjScjA5/Gpv\n3YWPWuJKp9T+YF6t0RCbiQEAMNT2WFQPdnTWXYt/OM/oACJJPT29dLuqAYQOAAAArIfRAYQ2vPlm\n1rkEa8/kbTVk0prbfUCJpu1VrckbnCZ0AABQg3JPQwfsZGwAoQ3vxniD09nN3O3hWd3v3Kf7nXtL\ntqjd6Ot0zI4b28YWAACUZp2GHk0mHK4E9cbcAEIb3rK1RkOrOkHtmhuX5/49Bbv2ab5zb0Veo2Nu\nwvHD+gAAwOZYp6EvpAggsFeD0wUUs9/nUyAQKHjNasO7v842opfaJOaOR/XI3Z8VvObKpNQ5ParH\nRv5SrdHQhl67ObmoPZO3tW/sFuEDAIAtoFbPyULtM3YGhDa86+dKp9Q9eXvNdnrWyeGxtp2afvQj\n69of4kqntGtugjMnAAAAUBHblpaWlpwuopSVbXh7enp1uGftNrwXL5zXle/9Vx04eDDv8a9/45tV\nq7Xa3p6+o5nY6ncr9o8Nb3hWItXkVrLZnf0496RtSfLcv0efcAAAtqjAY08q7m7TkT2Pqcvd7nQ5\nqBPGzoBIm++C9eQv/KJ+7df6q1mirQqFjz2Ttze1JKoxGc/bSM7yKgAA6genocMJxgaQSnTB2rdv\n35bepJ7b8QoAAKBcrHKAE8wNIHTBKqlQxysAAIByWKehR1MJdTldDOoGXbBqxGxOp4pSHa8AAADK\nRSte2MnYGRC6YBW3no5XAAAAa+E0dDjB2AAiLXesWtkFq7//pA6/unYXrK3KlU7RtxsAAFSEi03o\ncIDRAURabru7VrerejDzMHRsj0cdrgQAqrPB/AAAE0VJREFUAGwVnIYOJxi7BwQAAADVZa2qSGbS\nDleCemLsDEggENDdIpvQLfXYAas1GnK6BAAAAGDDjA0gt/xDunjhfNFOWJJ06733bazIWaFEzOkS\nAADAFuSOR5Xc3uJ0GagjxgaQo8eO63BPr45+9jN1FTTWQrcKAABQSQ3plIKLvNEJ+xi9B8Tn8606\niLDe0a0CAABUEq39YTejA4gk/dEf/4nTJRhhJra8SawpEXe4EgAAsJVsj9FhE/YyPoAgX2OSAAIA\nACovyH5T2IQAUkOa2f8BAAAqzNpfmqIVL2xCAKkBsw97dDey/AoAAFQY+0thNwJIDWGTGAAAqLRt\nBBDYjABSQ9gkBgAAKs06DX3m4b+BaiOA1AB+IAAAAGCrIIDUEHecGRAAAFB53GPATgQQAACAOteQ\nTilJFyzYhABSA0IP+3I3x1iKBQAAKq8puaiFVMLpMlAnCCA1hC5YAACgGppo9Q8bEUBqBIcQAgAA\nYCsggNSAmViEQwgBAEDVuONRzbDUGzYhgNQIll8BAABgKyCA1AgOIQQAANXCHhDYiQBiuFkOIQQA\nAFXWmFwOILTihR0IIDWCA4IAAEC1Wa3/gWoigAAAAECt0ZDTJaBOEEAMN/NwCRaHEAIAAGArIIDU\nCLpgAQCAampKLirIEizYoNHpAkq5dvWKLr3+mvz+obzHjx47rtNnzqqnp9ehyuzFIYQAAKDamhJx\nNqHDFsYGkMGBAX3plZf16le+qiN9ffJ4vJKkQCCgwcHrOvXcs/r+D36YfXyrCiViHEIIAACALcPc\nADJ4Xf0nTujoseN5j/t8Pvl8J3VryC//kF9H+vocqtA+TcyAAACAKqPjJuxi7B6Q/T6fAoFAwWvh\ncEiDg9e13+ezuSr7JTNpDgcCAAC2oA0v7GDsDEh//0ldu3pFn3r6E+rp7cnu9/D7h+Qf8uv0mbPy\n1UEACS7G1OV0EQAAYMtrSsTFmgvYwdgAIklf/8Y35fcPaXBgQOHwcm/q/v6TOvxqb12EDwtTogAA\noNqs09CBajM6gEhST09v3XS7WolOFAAAwE4swYIdjA4g9d6G1/oh0BINOlwJAACoB43hec3GI+py\ntztdCrYwYzehW214T585q7ffeVe33ntft957X9f+3z/Tkb4+nXru2eyyLAAAAFTGSHjG6RKwxRk7\nA1KJNryXXn9Nl15/Le+xW++9X5V6qyGRSbP/AwAA2KYpuaiZWIRZEFSVsTMglWjDe/rM2ezMifVP\nLQklYmpIp5wuAwAA1Amr9f9fz407XAm2MmNnQGjDu4xDCAEAgN0WUgmNRuZ1sL3D6VKwBRkbQKTV\nbXgXFxd15Min9eVXf68uwockDiEEAAC2yV36PRKa0d5Wr5oaXA5WhK3I2CVYly+/ocuX31BPT69O\nnzmrcCisf/en39TFC+f1zOc/p2tXrzhdYtXRCg8AADhlIZVgQzqqwtgZkLs5+z8uX35Dfv+Q3n7n\nXXk8Xvn9Qzr3whd1pK9PHo/XwSqrj03oAADALitXXoxF5nXI080sCCrK2BmQXHcDAZ0+czYbNnp6\nerXf55N/yO9wZdXFQYQAAMBOjcm4WqMfHnOQSKeZBUHF1UQAsWY9coVDYXm8HocqskdwMcYhhAAA\nwFZ7Jm/nfTwSmtFCKuFQNdiKjA0g+30++f1DOvyxx1ed53Hq+S/I4/Vs+ZPQAQAA7NaYjKtjbiLv\nMWZBUEnG7gHp7z+p/v6TkpbP/chdbnXkyKfVf+KEU6XZIskhhAAAwCE7Z8cV2rlbadfyreLog3nt\nbfVyOCEqwtgZkFwejzfvxPPc/SBbFYcQAgAAp7gyKXVNj+U9xiwIKqUmAki94hBCAADgFE/wXt5q\njJlYRLPxiIMVYasggBgqkUlzCCEAAHDUylmQv54bd6gSbCUEEEOFEjG5WIIFAAAc1BINyhuczn68\nkEpoNDLvYEXYCgggBmMJFgAAcFrH7Hjem6LvBac0sRAq8RlAaQQQg21jBgQAADisMRnXrpy2vIl0\nWjdmxjQcnHKwKtQyAoihQomY3Gz0AgAABtg1N67mFSszRkIzemtqhEMKUTYCCAAAANbUuWJDuiQF\nF2N6a2qEJVkoCwHEUK7wnNMlAAAAZLWHZ9UaXR00WJKFchFADBWmBS8AADBMoVkQC0uysF4EEEPl\nHvwDAABgAnc8oo6cDekrsSQL60EAMVAyk1YDHbAAAICBds6OlwwhLMnCWgggBuIQQgAAYCpXJqXO\n6VHtHxte1Rkr10hoRjfnx5XMpG2sDrWAAGIoDiEEAAAma4kGtf/2zZKzIaMP5vX29B1CCPIQQAyU\nyKQ5hBAAABhvPbMhwcWY/mzipwomYjZXB1MRQAzEIYQAAKCWWLMhOx4UPkYgkU7rR9N32JwOSQQQ\nAAAAVIArk9Ijd3+mR8ffL7iX1dqcPhKecaA6mIQAYqCmB/NOlwAAALAh7eFZHRz5y6JHCgzfn9LN\n+XGbq4JJCCAGesAhhAAAoIa5MintHb1VdEnW6IN5vTU1wub0OkUAMVBTjP0fAACgtllLsrzB6YLX\n2ZxevwggBsqkEk6XAAAAUBG7Jz9Q1/RYwWuJdFpvTY6wL6TOEEAMlFgIO10CAABAxeyaG9eeydtF\nD1q29oWwJKs+EEAMxCnoAABgq/EE72nv2HDR+xzr0EKWZG19BBDDJDNpNbEJHQAAbEHueER7x4aL\ndsgKLsY4L6QOEEAME0rE1JgkgAAAgK3JHY9kO2SVOi9kODjlQHWwQ6PTBWCFVNLpCgAAAKrK6pAl\nSakmt+It7Uo2bddCm1epZrcSTds1EprRRDSkvW1ePdrq1c7mFoerRqUQQAyTjtx3ugQAAADbNCbj\nan+4+mPX3PIBhemGRiVa2hV3t+mup1Mj7ja1NjYTRrYIAohhHrD8CgAA1DlXJqWWaFAt0aB2zY0r\n3dCoWPtOzbvbdLfNq4b2XTq4o0N7W71qbWx2ulyUiQBimKbwvOiBBQAA8CFXJqX28Kzaw7Pq1PKy\nrfk2ryZa2tXQtlOdXfvU7W5Xl7vd6VKxDkYHkGtXr+jS66/J7x/Ke/zoseM6feasenp6HaoMAAAA\nTmlMxuUJxuUJ3lt+YOSG7rrbdXd7q1rbd2rHjg51eXeradcjzhaKgowNIIMDA/rSKy/r1a98VUf6\n+uTxeCVJgUBAg4PXdeq5Z/X9H/ww+/hWEecQQgAAgLK54xEpHlEmNK2QpJAkt6tJ7javtrd61eBu\n01Jjk7w7urTN3Sptb5PcrU6XXZfMDSCD19V/4oSOHjue97jP55PPd1K3hvzyD/l1pK+v6HP8+J13\ndPHC+bzHzr34UlXqrZSGdFKcAQoAALB58XRS8fCsFJ7NPpbb3LfVtV3p7W6lmt1qbWzWkqtJ6Wa3\nJMnd2KTtrkZlPF0lX2Npe4uWtq8OMq2Nzcv7U+IL0mLOuSfe7k19TVuBsQFkv8+nwYGBgtfC4dDy\nLMiZs0U//9yLL60KH7Xg8U/9qtMlAAAAoFLcrcy0rLBtaWlpyekiijn1/BfkH/Krp7cnu9/D7x+S\nf8iv02fO6nSJAAIAAADAPEYHEGk5cAwODCgcDkmSenp6dbinVz6fz+HKAAAAAJTL+AACAAAAYOto\ncLoAAAAAAPWj7gPIxQvna2qzeq3VOzgwoFPPf8HpMtaNeqvv8Mced7qEslBvdZ16/gtFG46YiHqr\nq9bqrbXfydRbXbX4O9kpdR9AAAAAANiHAAIAAADANgQQAAAAALap+y5YFy+c14/feUeffPppp0tZ\nlx+/844k1Uy94+PjuvnXf6Xjv/J3nS5lXai3+i69/lpNneFDvdV15Xv/VU/+wi9q3759TpeyLtRb\nXbVWb639Tqbe6lrv7+RzL75kU0XmqvsAIqmmNjgBAACgdhFACCAAAAAAbMQeEAAAAAC2aXS6AKdd\nvvyG7gYCOnLk0zrS1+d0OSVZtVr2+3zq7z/pYEWFXbt6RYd7euXz+fIeN3msV9Zs6liHwyFd/va3\nFQ6HVtVk4vgWq7cWxnflOJo+vsXqtZgyxpaLF86vWoZg4hhbVtZr4viuVZNp41uqXhPHV1o+52Fw\n8HrBekwbX6l4vaaNr1XnSs/0n1z1e9mU8S1V8+DgdaPG10R1PQOSe+DRb7/wRV27esXhikozvT5J\nCgQCunjhfN5fPMnssS5Us0n15Tr13LPy+4ckSW9++9v68isvLz9u6PgWq9eU+nKFwyH9nb/9ywqH\nQ5LyD8AycXxL1SuZOcaWixfO69Lrr+U9ZuIYWwrVa1J9llI1mTi+pWowob6VLl44r0uX/rWk5Z9n\nuX/fTBzfUvWaUN9acv/OmTi+hVg1m1qfUZbq1NDQT5Z+6ZNPZT++euV7S7/2D/6+gxWt7YmfP+R0\nCSW98Fu/ufTEzx9aeuLnDy0NXL+efdzksS5Ws4ljPXD9et64hULBpSd+/tDS228PGjm+xeq1/m2a\nb3/73y89/9xvZD++euV7S88/9xvGfv8Wq9di4hgvLS1/Xzz/3G/k1WfqGC8tFa53acnM8S1Wk6nj\nW2oMTRvfsbGxpV/65FNLoVAw+/ELv/WbS0tLZo5vqXqXlswb35UufO0Ply587Q+XlpbMHN9Ccms2\nfXxNULdLsAYHBtTT25P9+Ehfn377hS86WFFp4XBIHo9Xl15/reDyGxP80R//iSTp8Mcez3vc5LEu\nVLOpY+3xevTMiROrHn/3x+8YOb7F6g2HHxg5vkeOfFo9Pb3ZjwOBgDwer7Hfv8Xqlcz9Hg6HQ/ry\nl35Xl77xTR397Geyj5s6xsXqNXF8S9Vk4viWqtfE8R0cvK6e3h6FQmH5h/za7/Nlf3+YOL6l6jVx\nfHMFAgENDgzoze98V5KZ47tSbs2mj68p6nYJVjgcyvvlbf3iNpV/yK9wOJStO3c5i+kY68ro6enN\n/hALBAI69dyzOn3mrBYXF40c32L1BsbGjBxfn8+nnp5e+f1DOvrZz+jS66/p3IsvGfv9W6xeydzv\n4S+98rLOvfjSqv1hpo5xsXpNHN9SNZk4vqXqNXF87wYCCofCOvfCF3Xp0r/WM5//nC5ffkOSmeNb\nql4TxzfXl7/0u3n7rUwc35VyazZ9fE1RtzMgteZIX59uvfd+3sefevoTOvc7v2PkX8ZaZvpYX7xw\nXteuXtHRY8d17sWXjD/HZmW9kowe356eXr35ne8u39C/8EUjNjuWsrLeN7/zXSO/hy+9/pq8Hq+O\nHjvuyOuXq1S9Jo5vqZpMVKpeE8dXWr4RfvM735XH41UgENDRz35Gxwz+fi5Wr6njK0l+/5DCobDx\nP3dzrazZ5PE1Sd3OgFjvHFoCKzZNm876Jq6Fuhnrynnm859TOBTWm9/5bvZm3uTxLVTvSqaM7+DA\nQHaTo8fj1bkXX5LfP2Ts+BartxATxnhw8LouX35Dhz/2eHbJ4+GPPb68vMLAMS5V70omjO9KuTWZ\nOL4rlRpDE8Z3v8+n/T5fthZrVsw/5DdyfEvVu5IJ42t589vfXhX6TRzfXIVqzmXS+JqkbgPI4Z7e\n7DSZtPzLxuR35i69/ppOPf+F7Md+/5A8Hm/etKSpGOvKuHz5DXm8Hn35K1/NexfF1PEtVq+p4+v3\nD2WXKFgfezxeY8e3WL2SmWP89W98U7feez/7j7Q8E3akr8/IMS5Vr4njW6omE8e3VL0mju+RI5/O\nG0OrJlO/f0vVa+L4Wq4+nC3PZeL45lpZs8nja5K6XYLl8/l0+sxZPfP5z+nIkU/r6tUr+vq/+bdO\nl1VU/4kTGhy8rmc+/zl5vB7dDQT06le+6nRZ68JYV8bdh5vcVm7yv/Xe+0aOb7F6337nXSPHt//E\nCZ167kpeXed+53eM/f4tVq91zcQxLsbUMS7GxPEtVZOJ41uqXhPH1+fzqf/ECT3z+c9pv88n/5A/\n+/fNxPEtVa+J4ystz+p6vZ5Ve65MHF9LoZpNHV/TbFtaWlpyuggnBQIB3Q0E1NPbUxNr86y1hrVS\nby7GuroY38qwltisrMvU8S1Wr2TuGBdj6hgXY+L4lqrJxPEtVa+J41tqDE0c31I1mTi+pZg4vqXU\n2vjare4DCAAAAAD71O0eEAAAAAD2I4AAAAAAsA0BBAAAAIBtCCAAAAAAbEMAAQAAAGAbAggAAAAA\n2xBAAAAAANiGAAIAAADANgQQAAAAALYhgAAAAACwDQEEAAAAgG0IIAAAAABsQwABAAAAYBsCCAAA\nAADbEEAAAAAA2IYAAgAAAMA2BBAAAAAAtiGAAAAAALANAQQAAACAbQggAAAAAGxDAAEAAABgGwII\nAAAAANsQQAAAAADYhgACAAAAwDYEEAAAAAC2IYAAAAAAsA0BBAAAAIBtCCAAAAAAbEMAAQAAAGAb\nAggAAAAA2xBAAAAAANiGAAIAAADANgQQAAAAALYhgAAAAACwDQEEAAAAgG0IIAAAAABsQwABAAAA\nYBsCCAAAAADbEEAAAAAA2IYAAgAAAMA2jU4XAACVdvnyG7obCKx63OPx6uix4/L5fFV53YsXzkuS\nnuk/WfQ11vNnnGBqXcX+X0rL/z/7T5yQx+MteP3a1Svy+4ckFf66ij33fp9Px44dL/q8lmqOman/\nPwCgEv5/DiWFLoeHPtMAAAAASUVORK5CYII=\n", "text/html": [ "
total variable sitesparsimony informative sites051015202530354045505560657075Position along RAD loci0500100015002000N variables sites
    \n", "
  • \n", "
  • \n", " Save as .csv\n", "
  • \n", "
" ], "text/plain": [ "" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "canvas, axes = SNP_position_plot(distvar, distpis)\n", "\n", "## save fig\n", "toyplot.html.render(canvas, 'snp_positions.html')\n", "\n", "## show fig\n", "canvas\n", "#axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Possible reasons for the uptick in SNPs towards the end of reads: \n", "+ Sequencing error rate increases along the length of reads\n", "+ Poor muscle alignment at the end of reads\n", "\n", "I think the more likely explanation is that poor alignment towards the end of reads is causing the excess SNPs for two reasons. First, if it was sequencing errors than we would expect an excess of autapomorphies at the end of reads equally or more so than we observe synapomorphies, but that is not the case. And second, increasing the minimum depth does fix the problem. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Apply a more stringent edge filter\n", "Positive values of the edge filter " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "## make a new branch of cyatho-d6-min4\n", "data5 = data1.branch('cyatho-d6-min4-trim')\n", "\n", "## edge filter order is (R1left, R1right, R2left, R2right)\n", "## we set the R1right to -10, which trims 10 bases from the right side\n", "data5.set_params('edge_filter', (\"4, -10, 4, 4\"))\n", "\n", "## run step7 to fill new data base w/ filters\n", "data5.run('7', force=True)" ] }, { "cell_type": "code", "execution_count": 169, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " loading Assembly: cyatho-d12-min4\n", " from saved path: ~/Downloads/pedicularis/cyatho-d12-min4.json\n" ] } ], "source": [ "data5 = ip.load_json(\"~/Downloads/pedicularis/cyatho-d12-min4\")" ] }, { "cell_type": "code", "execution_count": 170, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "prefilter cyatho-d12-min4\n", "shape (54652, 150, 2) = (nloci, locus_len, [var,pis])\n", "total vars = 158199\n", "postfilter cyatho-d12-min4\n", "shape (29859, 150, 2) = (nloci, locus_len, [var,pis])\n", "total vars = 86868\n", "[ 0 0 0 0 0 922 1050 1022 982 1066 1108 1018 1090 1089 1104\n", " 1134 1082 1085 1182 1131 1133 1134 1127 1119 1157 1168 1133 1150 1174 1082\n", " 1162 1161 1084 1178 1186 1146 1128 1151 1149 1177 1218 1142 1196 1146 1172\n", " 1226 1211 1167 1215 1181 1239 1244 1188 1180 1217 1242 1222 1297 1321 1245\n", " 1313 1337 1322 1404 1436 1507 1519 1603 1624 1727 1871 1878 2041 2199 377\n", " 182 132 78 41 24 16 4]\n", "[ 0 0 0 0 0 428 503 479 484 544 535 522 608 554 523\n", " 556 595 588 628 535 535 584 579 576 576 541 551 629 555 597\n", " 584 602 543 596 585 540 584 591 616 601 578 588 597 597 553\n", " 591 599 606 631 580 617 623 625 594 646 593 606 672 599 673\n", " 657 659 675 717 721 686 800 776 844 957 954 1024 1120 1182 198\n", " 100 65 33 20 11 6 1]\n", "maxlen = 81\n" ] } ], "source": [ "## filter the snps\n", "fsnps = filter_snps(data5)\n", "\n", "## trim non-data \n", "maxend = np.where(fsnps[:, :, :].sum(axis=0).sum(axis=1) != 0)[0].max()\n", "\n", "## all variables (including autapomorphies)\n", "distvar = np.sum(fsnps[:, :maxend+1, 0].astype(np.int), axis=0)\n", "print(distvar)\n", "\n", "## synapomorphies (just pis)\n", "distpis = fsnps[:, :maxend+1, 1].sum(axis=0)\n", "print(distpis)\n", "\n", "## how long is the longest seq\n", "print 'maxlen = ', maxend" ] }, { "cell_type": "code", "execution_count": 171, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAEsCAYAAAA7Ldc6AAAJNmlDQ1BkZWZhdWx0X3JnYi5pY2MA\nAHiclZFnUJSHFobP933bCwvssnRYepMqZQHpvUmvogJL7yxLEbEhYgQiiog0RZCggAGjUiRWRLEQ\nFBSxoFkkCCgxGEVUUPLDOxPn3vHHfX49884755yZA0ARBQBARQFSUgV8Pxd7TkhoGAe+IZKXmW7n\n4+MJ3+X9KCAAAPdWfb/zXSjRMZk8AFgGgHxeOl8AgOQCgGaOIF0AgBwFAFZUUroAADkLACx+SGgY\nAHIDAFhxX30cAFhRX30eAFj8AD8HABQHQKLFfeNR3/h/9gIAKNvxBQmxMbkc/7RYQU4kP4aT6edi\nz3FzcOD48NNiE5Jjvjn4/yp/B0FMrgAAwCEtfRM/IS5ewPmfoUYGhobw7y/e+gICAAh78L//AwDf\n9NIaAbgLANi+f7OoaoDuXQBSj//NVI8CMAoBuu7wsvjZXzMcAAAeKMAAFkiDAqiAJuiCEZiBJdiC\nE7iDNwRAKGwAHsRDCvAhB/JhBxRBCeyDg1AD9dAELdAOp6EbzsMVuA634S6MwhMQwhS8gnl4D0sI\nghAROsJEpBFFRA3RQYwQLmKNOCGeiB8SikQgcUgqkoXkIzuREqQcqUEakBbkF+QccgW5iQwjj5AJ\nZBb5G/mEYigNZaHyqDqqj3JRO9QDDUDXo3FoBpqHFqJ70Sq0ET2JdqFX0NvoKCpEX6ELGGBUjI0p\nYboYF3PAvLEwLBbjY1uxYqwSa8TasV5sALuHCbE57COOgGPiODhdnCXOFReI4+EycFtxpbga3Alc\nF64fdw83gZvHfcHT8XJ4HbwF3g0fgo/D5+CL8JX4Znwn/hp+FD+Ff08gENgEDYIZwZUQSkgkbCaU\nEg4TOgiXCcOEScICkUiUJuoQrYjexEiigFhErCaeJF4ijhCniB9IVJIiyYjkTAojpZIKSJWkVtJF\n0ghpmrREFiWrkS3I3uRo8iZyGbmJ3Eu+Q54iL1HEKBoUK0oAJZGyg1JFaadco4xT3lKpVGWqOdWX\nmkDdTq2inqLeoE5QP9LEado0B1o4LYu2l3acdpn2iPaWTqer023pYXQBfS+9hX6V/oz+QYQpoifi\nJhItsk2kVqRLZETkNYPMUGPYMTYw8hiVjDOMO4w5UbKouqiDaKToVtFa0XOiY6ILYkwxQzFvsRSx\nUrFWsZtiM+JEcXVxJ/Fo8ULxY+JXxSeZGFOF6cDkMXcym5jXmFMsAkuD5cZKZJWwfmYNseYlxCWM\nJYIkciVqJS5ICNkYW53txk5ml7FPsx+wP0nKS9pJxkjukWyXHJFclJKVspWKkSqW6pAalfokzZF2\nkk6S3i/dLf1UBiejLeMrkyNzROaazJwsS9ZSlidbLHta9rEcKqct5ye3We6Y3KDcgryCvIt8uny1\n/FX5OQW2gq1CokKFwkWFWUWmorVigmKF4iXFlxwJjh0nmVPF6efMK8kpuSplKTUoDSktKWsoByoX\nKHcoP1WhqHBVYlUqVPpU5lUVVb1U81XbVB+rkdW4avFqh9QG1BbVNdSD1Xerd6vPaEhpuGnkabRp\njGvSNW00MzQbNe9rEbS4Wklah7XuaqPaJtrx2rXad3RQHVOdBJ3DOsOr8KvMV6Wualw1pkvTtdPN\n1m3TndBj63nqFeh1673WV9UP09+vP6D/xcDEINmgyeCJobihu2GBYa/h30baRjyjWqP7q+mrnVdv\nW92z+o2xjnGM8RHjhyZMEy+T3SZ9Jp9NzUz5pu2ms2aqZhFmdWZjXBbXh1vKvWGON7c332Z+3vyj\nhamFwOK0xV+WupZJlq2WM2s01sSsaVozaaVsFWnVYCW05lhHWB+1Ftoo2UTaNNo8t1WxjbZttp22\n07JLtDtp99rewJ5v32m/6GDhsMXhsiPm6OJY7DjkJO4U6FTj9MxZ2TnOuc153sXEZbPLZVe8q4fr\nftcxN3k3nluL27y7mfsW934Pmoe/R43Hc09tT75nrxfq5e51wGt8rdra1LXd3uDt5n3A+6mPhk+G\nz6++BF8f31rfF36Gfvl+A/5M/43+rf7vA+wDygKeBGoGZgX2BTGCwoNaghaDHYPLg4Uh+iFbQm6H\nyoQmhPaEEcOCwprDFtY5rTu4bircJLwo/MF6jfW5629ukNmQvOHCRsbGyI1nIvARwRGtEcuR3pGN\nkQtRblF1UfM8B94h3qto2+iK6NkYq5jymOlYq9jy2Jk4q7gDcbPxNvGV8XMJDgk1CW8SXRPrExeT\nvJOOJ60kByd3pJBSIlLOpYqnJqX2pymk5aYNp+ukF6ULMywyDmbM8z34zZlI5vrMHgFLkC4YzNLM\n2pU1kW2dXZv9ISco50yuWG5q7uAm7U17Nk3nOef9tBm3mbe5L18pf0f+xBa7LQ1bka1RW/u2qWwr\n3Da13WX7iR2UHUk7fiswKCgveLczeGdvoXzh9sLJXS672opEivhFY7std9f/gPsh4YehPav3VO/5\nUhxdfKvEoKSyZLmUV3rrR8Mfq35c2Ru7d6jMtOzIPsK+1H0P9tvsP1EuVp5XPnnA60BXBaeiuOLd\nwY0Hb1YaV9YfohzKOiSs8qzqqVat3le9XBNfM1prX9tRJ1e3p27xcPThkSO2R9rr5etL6j8dTTj6\nsMGloatRvbHyGOFY9rEXTUFNAz9xf2pplmkuaf58PPW48ITfif4Ws5aWVrnWsja0Latt9mT4ybs/\nO/7c067b3tDB7ig5BaeyTr38JeKXB6c9Tved4Z5pP6t2tq6T2VnchXRt6prvju8W9oT2DJ9zP9fX\na9nb+aver8fPK52vvSBxoewi5WLhxZVLeZcWLqdfnrsSd2Wyb2Pfk6shV+/3+/YPXfO4duO68/Wr\nA3YDl25Y3Th/0+LmuVvcW923TW93DZoMdv5m8lvnkOlQ1x2zOz13ze/2Dq8ZvjhiM3LlnuO96/fd\n7t8eXTs6/CDwwcOx8DHhw+iHM4+SH715nP146cn2cfx48VPRp5XP5J41/q71e4fQVHhhwnFi8Ln/\n8yeTvMlXf2T+sTxV+IL+onJacbplxmjm/Kzz7N2X615OvUp/tTRX9KfYn3WvNV+f/cv2r8H5kPmp\nN/w3K3+XvpV+e/yd8bu+BZ+FZ+9T3i8tFn+Q/nDiI/fjwKfgT9NLOcvE5arPWp97v3h8GV9JWVn5\nBy6ikLxSF1/9AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAHXRFWHRT\nb2Z0d2FyZQBHUEwgR2hvc3RzY3JpcHQgOS4xMJremEEAACAASURBVHic7d1rcFt3fub5h1eBFwGy\nbrZlQd5O3IltMNndaccx6Zrdmk48JGc2teN1TCmVnbhtUZzdpOPYUiubmVi2LGcmU5HF7nhmUhVd\n3O2t2hoKSpdfpCYk02lX7fQKcJz2i7EJyYllu8gj6kaKAiCQBEGC3BfgOQJIAARIXA6I76dKtoAD\nAj/+AYrnOf9bzdLS0pIAAAAAoARqy10AAAAAgOpBAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAA\nACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQ\nAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQ\nMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEA\nAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVD\nAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAA\nACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQ\nAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQ\nMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMgQQAAAAACVDAAEAAABQMps2gPSfOqn+UyfL\nXQYAAEDV+3hyrNwlwEY2bQABAABA+V2bCenadKjcZcBGCCAAAAAoivnFuD6dGi93GbCZ+nIXkM3w\n0KDOnjmtQGAk5f7Orm71HuqTx9NWpsoAAACwlk+mxhWLxyUlwkhDbV2ZK4Id2LYHxO/z6Y3Xj6n3\nUJ8+/OinuvTZ57r02eca/tEHau/o0MEXv6VwmO48AAAAO1o59CoUmy1jNbAT2/aA+P0X1bN/vzq7\nulPud7vdcrsP6NJIQIGRgNo7OspUIQAAANJh6BWysW0PyF63W4ZhpD0WDofk91/UXre7xFUBAABg\nLclDr4CVbNsD0tNzQMNDg3rqySfkafNY8z0CgREFRgLqPdQnNwEEAADAVkYjU6x6haxsG0Ak6dy7\n7ykQGJHf57Pme/T0HNDjb7YRPgAAAGxmZiGmz4I30h6biEa009Fa4opgR7YOIKyCBQAAUDkYeoVc\n2DaAmKtgvXniLbV3dMjpdEmSDMOQ339RB1/8lv76b35s3Q8AAIDyuRKe0MRspNxloALYN4CwChYA\nAEBFmFmI6YvwRLnLQIVgFSwAAABsSC5Dr2YWYiWqBnZn2x4QVsECAACwv1yHXs0vMjcECbYNIBKr\nYAEAANgZQ6+wHrYOIJLk8bTltNrV449+fdV9vYf6ilESAAAAxKpXWB9bB5B8luG99NnnKY/pP3Wy\nJDUCAABUo3xXvZqeZw4IEmwbQFiGFwAAwJ7mF+N5D71iEjpM9g0gLMMLAABgSwy9wkawDC8AAABy\ndm0mpGvToXKXgQpm2x4QluEFAACwl/nFuD6dGt/Q1zfU1hWwIlQi2wYQiWV4AQAA7GSjQ69CsVnt\ndLQWsCJUog0FEMMwdNUw5GnzFG0yuNvtlnv/fiabAwAAlBFDr1Aoec0B8XoH9NSTT8jv88nv86nz\nmW/q4Esv6Kknn5DXO1Dw4o6/fkxPPfmEnnryiVXL6qbb9wMAAACFt9GhV0CynAOI3+fT8dePWUOh\nVgaO/rffto4Vgtc7oEBgRB9+9FN9+NFP5ff52NsDAACghIKxWU1GI6x6hYLKeQjW8NCgpMTu4u0d\nHTr40guJ+3/0gS54BxIbBhZwWdyrhpGy/0f/995R70svqPdQH8OxAAAACmQyGtGV5T09QrHZogaN\niWiEOSDIvQfEuDomSTp85KgVRjye4k0GdzpdKcvwut1udXZ1q//tt4vyegAAANVmfjGu/3Z7XBOz\nEU3MRujlQEnkHEDMXgdzVSpJRd0EsGf/fl01DD3/3LPWcK/DR44qEBixel8AAACwfp9MjbNDOUou\n5yFY7R0dGh4a1PPPPWvd19nVLa93QN7z5yVJnjZPwQpzOl268MP3rbBjuvDD9+X1Dlj7ggAAACB/\nV8ITrGqFssg5gJgbAyb3fng8bRoeGlQ4HCra3Ix0vSw9PQcK/joAAADVIhib1eU7N0r+uvS2QMpz\nHxBzY0BJVg9Ee/vT8nja1NnVXfjqAAAAUFDzi3F9PDFWttcG8toHRLoXPPw+nwzD0N7lyeEAAACw\nP+Z9oNzy6gHx+3w6/sZr1upUvYf6rPvPff8HLI8LAABgY8z7gB3kvgyvYejVV15OWRrXFAiM6NVX\nXi5oYQAAACiccs37SDY9T88L8gggF7wD1mTzS599bt1/+MhRdXZ1y+/zWfNDAAAAYB/lnPeRjKFf\nkPIIIGa4OHzk6KpjXctzQMKhcIHKAgAAQKFcDt7g5B+2sa6NCFei5wMAAMCeRiNTGr07Ve4yAEvO\nk9C7urqtjQjNVa/MXcn9Pp+cTldRd0YHAABAfkYjU/rk9ni5y0gxvxhXQ21ductAGeXcA9LZ1W2t\nejU8NCgpsfqVGT7ePPFWcSoEAABA3i4Hb9gufEhSKDZb7hJQZnktw2tOOB8eGlQ4FJbT5ZTT6VJn\nV7fcbnexagQAAEAePpkaZ9gVbCvnAOL1Dmh4aFDn3n3P2ozQlBiK9ZoOHzm66hgAAABKh/ABu1sz\ngHi9A7pqGAoERuT3+dR/6uSqx5jHwr2sggUAAFAuhA9UgjUDyPDQoPw+n3X77JnTGR/rdDkLUxUA\nAADyUinhYyIa0U5Ha7nLQBmtGUDMIVVXDUOGYWRc6aqzq5vhVwAAAGVQKeEDkHIIIObGg8lzQAAA\nAFA+84txhWKzCsZmFYrN6tp0qNwlATnLGkAMw9BVw9Bet1vt7U/LvXdfynCsldgHBAAAoLDMsDER\njWhmIabg3Cy7mqOiZQ0gF7wDOnvmtLX/R7b5H5J06bPPC1cZAABAlbscvKEroYlyl1FQhCfktQ8I\nAAAAim9+Ma6PJ8c0MRspdykFN78YL3cJKLOapaWlpXIXUQzmcsHmHBYAAIBKEIzN6uOJsU3bU7Cr\nqVVP7f5auctAGdVu5IsNw0js/xFm4hMAAMBGjUam9Le3vtq04QOQ8gwgXu+AnnryCfl9Pvl9PnU+\n800dfOkFPfXkE/J6B4pVIwAAwKb3ydS4Prk9rlh8cw9Rmp4nXFW7nAOI3+fT8dePWb0dKwNH/9tv\n0xMCAACQp/nFuH5y40rV7ONB7w5yDiDDQ4OSpN5DfWrv6LBuD//oA/Ue6lM4HFJgJFCcKgEAADaJ\nmYWYJqMRXZsJ6XLwhj649vcKzs2WuyygZHJeBcu4OiYpManbDB8eT5vcbndxKgMAAKhwV8ITmowm\nVrIKxWY3/fAqIBc5BxCn0yVJCgRGrM0I2XgQAABgtWBsVp9OjdOzkcH8YlwNtXXlLgNlknMAMYdd\nPf/cs9Z9nV3d8noH5D1/XpLkafMUvkIAAIAKshk3Dyy0UGxWOx2t5S4DZZLzHJCengMpPR7tHR3y\neNp01TAUDofUe6jP6iUBAACoNsHYrH5y4wrhA1hDXjuhn3v3PQUCI5IS8z8kqb39aXk8bers6i58\ndQAAADY3vxjXlfAEwQPIUV4BRLoXPEzMAwEAANVqMhrRf7s9ztKyQB7yDiAAAADVbDIa0UQ0omvT\nIYLHOk1EI8wBqWIEEAAAgCzmF+O6NhNSKDar6zMhltIFNogAAgAAkEYwNqvPgjc0MRspdynApkIA\nAQAAWOFKeEKX79wodxnAppTzMrzpGIYhv8+ncDhUqHoAAADKZn4xro8nxwgfRcbcmeqWVwDxegf0\n1JNPyO/zye/zqfOZb+rgSy/oqSefkNc7UKwaAQAAii4Ym9WHt77StWkurBbb/CLzaKpZzgHE7/Pp\n+OvHrN6OlYGj/+236QkBAAAVaTQypb+99ZWCc7PlLgXY9HIOIMNDg5Kk3kN9au/osG4P/+gD9R7q\nUzgcUmAkUJwqAQAA8jAZjehKeEKT0YiCseyh4pOpcX1ye5zVrYASyXkSunF1TJJ0+MhRK3x4PG1y\nu93FqQwAAGxawdistjU2FeS55hfjmohGFIrNKhSbzbhqVXN9o1oaGtVc36iG2jq5Gpv0RXiCXo8y\nmJ5nDkg1yzmAOJ0uSVIgMCK/zyeJXdABAEBuzM37zJAQi8e1q6lV39i5Tw21dXk/38xCTFfCE5qY\njeQ8oXlmIcbkZ5vgfahuOQcQc9jV8889a93X2dUtr3dA3vPnJUmeNk/hKwQAABVnNDKlmYVYYghU\nhh6GidmIPrj29/qF7Q9pT7Mr5+e+HLyhK6GJQpUKoMRyngPS03MgpcejvaNDHk+brhqGwuGQeg/1\nWb0kAACgOpnL2H5ye1xXQmsPb4rF4/p4YkyXgzfWXBlpMhrRT25cIXwAFS6vjQjPvfueAoERSYn5\nH5LU3v60PJ42dXZ1F746AABQMa7NhPTp1Pomc18JJSaM/8L2h1bNDZlfjOtKeILgscnML8bXNfwO\nlS/vjQjN4OH3+WQYhva63YQPAACq2PxiXJeDN/TxxNiGVpIKzs3qJ9ev6Er4XtC4NhPSf71Or8dm\nFFpjdTJsXnn1gPh9Ph1/4zUZhiEpsSSvef+57/+AIVgAAKxTMDarheUhSK7Gpoq5MhyMzerjibGC\nTiq+fOeGJqMRNdTWsSkgsAnlvgyvYejVV15Ou9lgIDCiV195Wefefa+gxQEAsJlcmwnltOxrY12d\nHt32gB5u3V6SujItY7ttS5NcjU1qrm/UtsamVcGomJPBMy2lC6Dy5RxALngHrMnmh48c1eOPfl1S\nYl8QwzA0PDSoQGDEGqIFANg8ZhZiaq5vLHcZGZlzBMwhHeaJ8rblk+dy125uipfrSXUsHrcmcT/i\n2lXwIDKzEEsJHJkCUXBu9bHm+kZt29KkmYUY+2cAWJecA4g5+fzwkaOrjnV1dWt4aFDhULhwlQEA\nym40MqXrMyFNzEa0bUtT2gnC5WQGj7HIVMrcg3Qn+tu2JMLIg82uvJZ8NY1GpqzX2dPi0i5Hq3Y6\nWrN+zcxCTJeDN9Y9jGhmIVawIDKzENO1mZCuz4Q2FBzYSwOFMhGNrPkzhM1pXRsRruzlMMMJAKDy\nBWOzGlsOHskn9cG5Wf3tra/0s85desS5K6/nNE9+H27dXpC5DZmCRzbm1fxr0yFdrm/UnhaXHm7d\nnrV3ZH4xrtHIlEbvTqWcdF8JJVZkaqyr04PNLu10tKaEGnOTvNG7U+v/JpOYQeT6TEiPOHflfNI2\nvxjXtZmQJqMR5lIAsI2cA4jZy/H8c89aq14FAiM6+NIL8vt8cjpd7IwO6N5EUq7qVIf5xXjKSi4T\n0XtX3mcWYtYyk+u96l4q5on2WlfHY/G4Lt+5oVBsVr+4/aE1w8TKE/HLd27o4a3b9Yhz17qGRa0n\neGSsazlEPLx1u/Ysh4jk48k9HpnE4nGN3k0ElE+Xw0hDbd2G68tkYjaiidmIGuvq5ErqiUr+u9mu\nhA7Yyda7t3V3645ylwGbqFlaWlrK9cH9p07q7JnTq+53Ol1688RbtlqOt//USUnph4wBxTIZjejj\nycQylLuaWvO6UlkpklfqKcbY+uSx6TMLMT267QFbDflJNr8Y14e3vsp5OEtjXZ32tW7Xg80uW3xP\nk9GIJqKRrDtVZ9NYV6dv7NyX9jOeSw/Aw1u3a1/r9jXbwvxMzCzEinZiLyU+z3taXIkwVqCeC6Da\nNc7Paff1L9U0HdTnj7WnHHvEtUuPbXugTJWhnPJahvfwkaPqTJrv4XQ55XS61NnVLbfbXfDihocG\ndfbM6VVDvDq7utV7qK9gE94Pv/p7Kbf7v/unHOd43sfNCZm/+YdHJN27UvkX//5PtbVhixrr6m1d\nf/LxYGxW12dCml+M69QfviFJmosvSJL1/ZnO//F35aivV1Ndoxpq6/J+/W+//DuKLs4rvrikucUF\n/ca/Ppxy/J3Xf1vbG5vX1X6Xgzf05u//gXW7vrZW3z5xLOXxf37838nZsEV1y1fyc33+mYWYPp4c\n0396/Y9Sjq9sn//n355adfxKaELbtjTpwWaX/uPrb6m2piZt/Wt9f/kej8UX9DsnjqWscpSuvrXq\nv/d8cX375W+rtbFRrfUO1dbU6I9OnkwJHtm+fvTulP7dH/yhttTVWz8j/d/905QQ+if/5pgWFhfX\nVR/HOc7x8h2viy9o4M0/1pa5Geu+X/+jf6RYwxbr9p/8m2O6b0uzddsOv/9KebyaZQ0ghmHo6vKe\nH8na259OuX11+XGFHILl9/n0xuvH9OaJt9Te0WHNQTEMQ37/RR188Vv667/5cUH2HvnN3/yXHOf4\nho7fmA3ri3D6pSi/8b88I0lqbdii+5ucZakvl+PR+LwmoxH95MaVlKvhZv2ZPPlrndbfHXUNuhy8\nkXKF33z+YGxW0fi8LgcTw3em5xMTWX+u83/K+vxP/VqXJOlnnbv0QI7tN78Y18eTY5qYjajjf03t\nmV05Ofl//Oe/qvraWj3cuiPn5w/GEnMhYvH4qudfKdNxcz7CzzzztHY6WuWoq1dL/RYFY7MpPQKZ\nvj/zBP1/+Ge/osjCnHWS/pejn6Y89r/71dR/l1cumbre+lceb23Yopb6Lfrx+N+v+/lbGxJfnzzX\nIvnztZH6OM5xjpfu+Pbb17RtclzP/vKjKffXx6IpAeSfPPtr8ty3J+Pz2+X3Y7GOV7OsQ7AyDbnK\n5NJnnxekKPO1pcxDqI6/fkydXd0ZQw9DsJDMPGGbXB6fv97x5+l8MjWe13ANc5hHOuacgXyeK3kM\nvquxSY1Jt9ca/mVOUB2LTBV8Oc3m+ka1NDQqFJst2JCZXLrrk4fB5WtPi2vNeQ3XZkL6dGq8aMOA\nTOaKTc31jXI1JpY8NZdLZQUiAHbUPB3SjltjckTTLzc9vu9xzST9/tvV1Kqndn+tVOXBRvIaglVK\ne91u+X2+tMfC4VCiF2R5J3YgnWBsVpPLoWPlVe/Ru1Mbmghryjd8SPcmv5aaGQgkWSfYxZygWoyl\nOq+EJjS/GNdj2x5IGxI2uinatenEBOz/fsdDacPbaGRKn9weX/fz5yPd/gsAYFf3X/9SzuDNrI9p\nng6lBBBUr7wmoZfawZdeUGAkIE+bx5rvEQiMKDASUO+hPvVmCSD0gJSWubb8zzp3lWVyrXmyO718\nlXjl8qHZrCeIzC/G9cnUOCvMlMm2LU16avfXrBCSPOSqUFb2thRzx2cAqFR18QXtGbucsdcj2Z0d\nD2ly9z7r9rYtTfrHDzxSzPJgU+vuATEMoygTz5Ode/c9BQIj8vt8CocTJ3o9PQf0+JttRX9trM1c\npvLadMi60j0ZjegXtj+0ruVGVy5nmkkwNms9dn4xvuGrxOYSmrkGkXxXPkLhBedm9V+vX9E3du3T\nwnL4KPSQqCuhCU1GI/rGzn0F3c8BADaLfMKHJDXMz6Xc5vdo9corgBiGoXNnTsvrHbDuc7vd6uzq\nLlpPg8fTVrDVrlAYZm9Huqv/sXhcH0+Maea+B/LaqCwYm018XRnHtptBZFfT6qE3hQg6KKyZhZg1\nEbxYgnOzqyZVAwAkR3RaD179B9XPR3P+mrrl1RSBnAOIYRh6/rlnrZ6I5PvPnjmtcCis4yfeKmhx\npVqGF7m5HLyR8xr8l+/c0MxCTL+4/aGcntdOQ1sKOYwHxVXsieAAgNUc0WntGb2kusX8AkVDLPew\ngs2tNtcHXvAOKBwOyeNp03e/944+/OinGv7RB9Y8DK93YFVQ2AhzGd7eQ3368KOf6tJnn+vSZ59r\n+EcfqL2jQwdf/NaqMFTN8l09KR/zi3H95MYVXQlN5HXCN3p3Sh/e+ipjXcHYrPW8AADA/lzBW+sK\nH5Ly6i3B5pZzD4gZLnoP9Vk7njudLh0+clThUFhe74DCoXDBCvP7L6pn//5Vu6u73W653Qd0aSSg\nwEigoHuPVKrkPQl2NbXK1dikXY7WguzAvZElTaVEb8KHt77SL2x/KGVyut16PQAAQHau4C3tvv7F\nhp6jLr6geN2908/JaKQg5yuoLDn3gOQy3MnpWr2J13rtdbtlpNkEUbq3DO9eJqJrfjGesifBxGxE\nV0IT8t/8Sn85+qk+vPWVroQnFMxhcvdKV8KJ59noMJfgXCIgTUYj9HoAAFCBdt4a23D4kKQt0ekC\nVINKl/NO6B5Pm5xOl954/ZhC4ZDcexPLqPn9F+X1DhR8snhPzwENDw3qqSefyLgMLythSR9PjmWd\nHD0xe28PDHMDvIdbt2dd6Wl+Ma7LwRsFXfUnFo/Lf/Orgj0fAAAovsb5Oe2+/qWapoPlLgWbiG13\nQjetXIbX42nT457Vy/AefOmFlNtjo6Pq/mf/fFPvA7KeTfBM27Y0aV/rdu1pdqVs6BaMzerTqXFW\nfAIAoMrtvDUm552b65rvkcnKvUDa7/8aQ7CqkG13QpfyWwWrt/dfpTzmL/7CW5Iay2U0MrWhHorE\nLsvj+uT2uB7eut364U8ezgUAAKqPIzqtXde/zHl/j40IxmYJIFXItjuh+30+vfrKy3rzxFtq7+iQ\n05nY2M4wDPn9F9X/9tv667/5sXX/Spt5J/TJaIThTAAAoOB23hrTfbfHi/b84W336+aDP2PdfsS1\nS49te6Borwd7KlgPSKF3RmcVrPSCsVl9PDlW7jIAAMAmUqpej5W7oaM65RVA+k+dzLjXh9/nK+gc\nkL1ut/w+X9pj5ipYB5f3IKkWK1e8AgAAyKYuvqA9Y5czBouoo1VLdfUlm2Rew27oUB4BxJyPUSqs\ngrXaWiteAQAAmNYKH5JKMs8j2+sVaxNl2FvOAcTsjeg91CfDMDQ8NKjjJ97SVcPQ2TOnrR3RC+nc\nu++tWgWrp+eAHn9z9SpYm93l4A1rOV0AAIC1PDj+eckDRr5mFmLlLgFlkHMAMa4m5h0cPnJUgcCI\nhocG5d67Tz09ByRJ3vPn1XuoL+Ok8PUq9P4ilSq0jo0EAQBAdbrfxnt3NE+HNNNS2PNFVJacd0I3\n+X0+KxD4/Ret+8PhkAIjgcJVhhT0fgAAgFzcf/1LOYM3y10GkFHOPSDt7U9bS+Oe+/4P5PG0yXv+\nvDUno9CSd2HPWFOVrYAFAACQTSWED0d02uoBYQ5Idco5gPTs36/hoUEFAiMKh8LqPdSnV1952Zob\n4na7CxoILgVG1H/qpIwsIaQYO6/bUZDhVwAAYA2VED4kqTZpJSwW16lOOQcQp9OlCz98PzEEq80j\np9Olc+++J7//opxOl3r27y9oYZ1d3Xrc06bOZ75ZNUEjkwWuDgAAgCy2375WEeFDSqzOheqWNYCY\nw6D2Lq84ZQ6JSh5y1d7+tHVfoYdEud3uVRsRViN6QAAAQCbbb1/Tjluj5S4jZ2xGiKwB5IJ3IGWJ\n3bX2ASlGT8V3v/dOwZ+z0jA+EgAArNQ4P6fdNl7tCsgkr53QUR4EEAAAkGznrTHdd3u83GWsy8rA\nNBmNaKejtUzVoByyBpDDR47q8JGjKbdRemzSAwAApMQeGjtujdl+g0Egm5z3ATn40gt6/NGvq//U\nyWLWgzToAQEAYHO6//qXcn/1qXbeGpMreEuNGeZH1MUXtPPWmB4au7QpwocjOl3uElBGOQ/B6uzq\ntpbcRWmxRB0AAJtP8spVyaEiXluv2dZtmm/YopkWl+oWF7Tz5pjq56PlKrXgalkJq6rltRFhZ1e3\nvOfPa6/bLffefasfw8aABUfvBwAAm48reCvjylV1iwtqDU9KUsXO81hL8kpYwdgsc0CqTM4B5IJ3\nQMNDg5Kk468fS/uYat+voxhCLMELAMCm4gre0u7rX5S7jLJqiN3rzeFia/XJeQ4IAAAANobwAeTR\nA7JyRSyUxsQmmGgGAAASE68JHwlMQq9uBdsHxDAMuZd3TAcAAKh0W+/e1pbZeyfKDfNzqkszeXqm\nxaWIc4diDVsyPpcjOq09o5eKUmelYwhW9ckrgPSfOqlAYCTtMb/PxxyQImAOCAAApdM4PyfnnZva\nGr6d86pTTdNB7bg1qqijVdPOHavCiBk+6hZZ+cmUvBkh+51Vn5wDyPDQoM6eOV3MWgAAAMrCFbyl\n5umQtfrUejiiETmiEe24NaqFBofuOndorqlFu659SfgAkuQcQMw9QHoP9ckwDA0PDer4ibd01TB0\n9sxp9R7qK1qR1YweEAAAiqN5OqTm6ZCcd24WPCDUz0c37RK6hdI4P5d12Bo2r5wDiHF1TFJiMnog\nMKLhoUG59+5TT88BSZL3/Hn1HuqT0+kqTqVVKhZnXCQAoLLsvDWmreHbmtr5kELbduf99Y3zc7r/\n6j/IEY1oocGh+UaHoo4WSYn5Fot19dbttTii02qYj2rL7LTq4gtqmJ9T42yEHgkbqI9FFWvYwhyQ\nKpT3JHS/z2dtOOj3X7T+Hg6HFBgJsBlhATEmEgBQSZqnQ9pxa8za1Xv39S/UPB3SrQe+pnhdbqcc\n229f07bJcSsg1M9HVT8fteYMJPcqLDQ4NtXu4NXGfI+Dc4z2qDZ57YTu9/n06isv69z3fyCPp03e\n8+cVCIwoMBIoZo1ViwACACikuviC7rt9Tc47NzXt3KE7Ox8qyBCYuviCdt4akzN4c9Wx1vCkmiJB\nTez5Gd3duiPrczw4/nnK5OS1ED4q25bZ6ayfCWxeOQeQnv37NTw0qEBgROFQWL2H+vTqKy9bc0Pc\nbje9HwU2TQCpOObqKRHnjpyHB1QrR3RareHburNjT85XRoFyqIsvaEt0Ws3TITmi05pv2KLQffcX\n/GfcfJ21OKLTmm/ckteJW3LwMK86O4M33KnUwgAAHXlJREFU5QzeVHjb/RsKIq7gLe24OZp1SFPd\n4oIeuPoPanXuTNsbsvXubSZqA1Uk59/6TqdLF374vvw+nzxtHjmdLp179z35/RfldLrUs39/Meus\nSvSAVIa6+IJa705pa/i2deXOeefmmlf7qlnyEIut4du6+eDPaKaF+WOFtvLEOeposcbOLzQ6mPyZ\nQeP8nLZEI9oyO62m6ZA1nMjUpMTJ+2zLNk3t2JPzZ7cuvmCttCRpw/MQHpAUce7UnKMl4x4U6YLH\nSusNIo7otHbeGsurx8LsDbmx9+c00+JSXXxBu298taGVp1C52IywetUsLS0t5fLAQGBETqerYjYb\n7D91UpIqevf2y8EbuhKaKHcZJbPz1picd24qfN/9FXFVvHk6pK3h22oJ3874i/3Wgz+7rgmYm1W2\nk407Ox7K+31vnJ9TbXwh7yvRjfNzqo9FrRPB8H33r/tk3NyUzA6f18b5OTVNh9QQi6Y9cU4n6mjV\n0nIomW90aLbFVZRgYrb1WooRRJNf2xGdVm3SRnIrN5bL52TalC2ImBcoNrq8ay6Sl32dc7Tqvsnx\nrP8+ZZIcRMyflbrFBWtDPkd0WjXxhZw+X2u9TvN0iGFUVWy2ZZuu7ntMktR+/9e009Fa5opQKnnv\nA+LxtKmzq1udXd0VE0bW8j//49ShY//vT3y2OG4uwftHv3Ew5fhr//lcyu1KP/5vf+OgauJx1S4l\nVsE4d+K3U4LIm//7vyp6fTWSahbjql1a0rk3/w/r2EKDQ33/5nvW7aWaGvX/hxMpG1QdfP3PUp7v\n3Inftv6++/oX+oO+31e8rq6g9ZsnBYt19XrtxZfz/vpyHK9dWlTtYlzn3vw/U46ntF9NrY7/33+W\nEijSPf/Wu7fVGp5Sa3jy3tfX1GhJNfpPb/9Bymo5rx18RVpaUs2SJC3p3aT313r9mhot1tRqqbZO\nf5jj92deWf7Ot1+TzOs4NTX6j2//a0myajj8u29oqaZGSzU1Wdtn5fMnH6/RkmqWX6P/nTdTjh/5\n9rHES2tJWlpK+fxZ31+Slcd/5/f/ZNXxeG29Yk2tijpadPjlN7RUUyvzSlXy978lOr3cvtLyf/Qf\n+v/QOqFvnI2o77V3sr7+yvr+/I8P665zhxUKc/181SwtqWZpUX928v9KOTle6/sv6PGaWv3x6ZOK\nbN1uhY7fe+VE6V6/QMfNHpFiv74zeNOW3z/HS3f829/5Yy3UNUiSGmvr9ZP/zx7nX6U6Xs1y7gHp\nP3Vy1UaEZhjp2b/fdsvvboYekJ/cuLLpV4ZwRKe16/qXGa+kxWvri9ojUoiNp3IR3na/Jnfv2/D3\n0Dwd0vbb19a8Smte1c5kvmFLTrXMNbUUZBjZzltjea2Hf2fHQ5rcvS/lPnN+TTHW6zfFa+t1+/6H\ns/Za5TKkJZ2oo1Xh++7Pq0ds693bct25ta6r8oVmLoVayloizp0KbdudsVekcX5Orcu9kBu9Gg+g\nPD5/rF0SPSDVJucAIkmGYcjvvyi/zye/z6dw+F63dmdXt777vXeyfHVpbYYA8pejn5a7hKLKZ9Jh\nIYNILkOniiHqaNW1fY+tq35X8FbKHJNSitfWa9q5QzOtrrzDSOP8nHZf/3JddUcdrbq59+fUtPx+\nlfJ7X2hwrJqXst7gsZL5Wc407MsMWsm9bNVuocGhUFJ4WznnCkDlGn3kHynWsIUAUmXyCiArDQ8N\nqv/USRmGIUm69NnnBStsowgghVcXXyhYL0S+V8RN5slwtKlVka3bc65n693bao6Eyj7eOOpo1cSD\nP5PTnAVz7Pj2yXFbnYhGnDs10+JKaX9zXH3yOPG6+ELJQ16hzbZsU+i+3doyO12UnpfkK/xmME63\njCkAbFbj+x7XTItLj933gB5x7ip3OSiRvM8mzV3Q/T6fAoER6367DcGqdJM2Gk7QOD+n+ybHU06M\nzOEYUupwnrmmFsVr6zXnaEkbDja64knd4oI1Pnn39S8UdbRqtsW1atnb5ImfTZGgbU6CHdGI9oxe\nUvi++9d8bDGHGm1Ea3hSreFJ7b7+RblLKbqm6WBRr7KbbQkA1Y7d0KtLzgHE6x3QuTOnrd4OKRE6\nurq61d7Roc6u7qIUiPLJNuTE3JlWSixJmYkZVMyQkuvKPLlyRCNyRCO67/a4FhocmmlxqWF+ztZD\nM+oWF9bV+wMAwGbTPB1iGfYqlHMAuWoYVvjoXA4dXV3d9HwUSTBW3snnuWwslQszqGQLKYVSPx+V\nM2ifoUoAAABYLecA4vG06fiJtwgdJZJPV2TzdEj3X/9S840ORR0t1hrw61nL33wuO805AAAAm1PD\n/JwkhmBVm5wDCEOsSivXXdCTJ3PXz0dThh4tNDgUbWrVnKPF2gE5eQMuKfGD3xBLhA1HdNrWQ5cA\nAMDmYu4blOt5DzaH8m/di7TWuhKw1v4ZUiKQtM5HmeQKAABsybwIiupCALGp6fnMVwK2376mbZPj\ntlwhCQAAIFcM+a5OBBCbStcVWRdf0IPjnzNMCgAAbBp18QXmgFQZAogNpfshzGfXcAAAgEqxJTqt\nYIE2OkZlyPpu+30++f0Xc36ySt513E5CaZbg3XlzjPABAACAipc9gPgv6uyZ0zk/GQGkMGIrekAa\n5+cYIwkAADYlNiOsPlkDyF63W+0dHRmPB0YCCodDBS+q2q3sAWmapo0BAACwOWQNID09B9TTc2DV\n/YZhqP/USSt8tHd00PtRRI7ZzEvtAgAAVDJzL5DJaEQ7Ha1lrgalkNeMn3A4pLNnTst7/rzC4ZDc\nbrcOHznKJoUFtrIHpJkeEAAAsEmZu6FPEECqRs4BZHhoUP2nTsowDDmdLvUe6lPvoT45nYzZKybm\nfwAAgM2sZrkH5EpoQg+3bldzfWOZK0KxrRlAAoER9Z86Kb/PJykxLOvgoT653e6iF1etJpKGXDH/\nAwAAbGaO6L3znivhCf3i9ofKWA1KIWsA8XoHdPz1Y9bt9o4OOV1OXfAOpH0880AKj/kfAACgWoze\nndKeZhdDsTa5rAHkqmGk3Pb7fFZPSDoEkI0LMv8DAABUmeSleK+EJwggm1xtuQtAqoWkPUCY/wEA\nAKrNxGxE12a4ALuZZe0BOXzkKL0aJTa9ELP+zvwPAABQDRzR6ZTNCC/fuaE9zSx0tFnRA2IzM0kB\nhPkfAACgGtQur4RlmlmIaTQyVaZqUGwEEJuZTxqCxfwPAABQDepWBBBJ+ix4I+W8CJsHAcRmzB4Q\n5n8AAIBqYW5GmCwWj+tKeKIM1aDYCCA2YyZ95n8AAIBqNxaZohdkEyKA2ExwLrEML/M/AABAtWia\nDqa9PxaP63LwRomrQbERQGyE+R8AAACpRu9OpSzSg8pHALGR0PImhMz/AAAA1SbbxddPpsZLWAmK\njQBiQ8z/AAAA1eaBq/8gR3Q67bGJ2YgmowxP3ywIIDYysfyDxfwPAABQbeoWF7Rn9FLGnhD/za/0\nydS4gssjRlC5CCA2xPwPAABQjeoWF/TQ2CW5grfSHh+9O6WfXL+iD299RY9IBasvdwG4Z2YhxvwP\nAABQ9XZf/0KSFNq2O+3xidmIJmYj2tXUqn2t27Wn2VXK8rBBBBAbmV+MM/8DAABAiRDSEItqcve+\njI8xg8jl+kY94tqlh1u3l7BCrBdDsGwkFJtl/gcAAMCy+26P6/7rX675uJmFmD65Pa4fj/+9RiNT\nJagMG0EAsZFYPM78DwAAgCTO4M2cQohEEKkUDMGyCeZ/AAAApOcM3lRjdFrTzh2KOHco1rAl6+PN\nIHIlNMHQLBsigNjEzEKM+R8AAAAZOKIROaIR7bg1qqijNacwQhCxJwKITcQW48z/AAAAyEG6MDLT\n4lLU0ZL28WYQ+Sx4Qw82u7Svdbu2NTaVuGqYCCA2EYrNMv8DAAAgT1YYkbTQ4NBd5w7NNbXo7tYd\nqx4bi8c1endKo3en1FzfqD0tLj3Y7CKMlBgBxCbq5mZVw/wPAACAdaufj+q+2+OSpF219Zpt3aaZ\nFpciW7crXpd62juzENOV0ISuhCasMLLL0aqdjtZylF5VbB1AhocGdfbMaQUCIyn3d3Z1q/dQnzye\ntjJVVnhzU9fVUO4iAAAANom6xQW1hifVGp7U7utfKOpo1WxLYsPCmZX/TwojktRc36iWhka5GpvU\nUFunbY1N1t+xcbYNIH6fT2+8fkxvnnhL7R0dcjoTHxDDMOT3X9TBF7+lv/6bH1v3V7qG2bvlLgEA\nAGDTModqSbJ6SUxRR6uW6upT5pBEm1o0XZs4Vf4HR4vVg7Jty70g0lzfmBJKdi33njTXN6q5vrF4\n30yFs28A8V9Uz/796uzqTrnf7XbL7T6gSyMBBUYCau/okCT1nzqZ8ri/++gj/dKTT5as3o1auHPT\nvm8GAADAJmYGk6bp4JqPXWhwaL7RIUmKL/+RpKijRSu/Ot7sVMtyqNm6/DWPuB8vRMkVzbbnvHvd\nbvl9vrTHwuFQohfkUF/Gr/+lJ5/U4SNHi1VewbX9ym+VuwQAAACg6GqWlpaWyl1EJgdfekGBkYA8\nbR5rvkcgMKLASEC9h/rUmyWAAAAAALAfWwcQKRE4/D6fwuHEErUeT5se97TJ7XaXuTIAAAAA+bJ9\nAAEAAACwedSWuwAAAAAA1YMAAgAAAKBkCCAAAAAASoYAAgAAAKBkCCAAAAAASsa2GxGW0spd1AEA\nAIBiqKSNsoul6ntA+k+d1N999FG5y8hqfHxcg3/1X8pdxpoG/+q/aHx8vNxlZPV3H31k+/dbks6e\nOV3uEtZUCe93pfzs8H4XBu934VTCv5WV8n5Xws8O73fhVML7bQf0gEj6pSeftHUa9ft8Cgbv2LpG\nKbFp5K//eo/aOzrKXUpGZm+X3dvy7JnTtq+xEt7vSvnZ4f0uDN7vwqmEfysr5f2uhJ8d3u/CqYT3\n2w6qvgcEAAAAQOlUfQ/IZ5cv6+7du/L7fLZMq17vgP72ww81Njqq/lMntdftVk/PgXKXZRkeGtTj\nnja53e6U+73eAV01DLW3P132dk1Xo1mfqZztGg6H5D1/XuFwaFUddmnHXGo0lastk2tM1152bMuV\ntdilLZP1nzq56oqjXdoyWf+pk2pvf1qSvdoxXS3pjpezLbPVaKe2lBJXwf3+i9rrdsu9d591vx3a\n0ZRco93+rTRrWykcDsvpdEqyR1umq3N8fFxz0TlbtKOpUn7v2FFV94AcfOkFXbuWGKf36isva3ho\nsMwVrWbHmkyGYaj/1MmUfwikRLv6fT5J5W/XTDXaqV0PvvgtBQIjkqQL58/r+OvHEvfbqB0z1SjZ\noy3D4ZD+6a/+isLhkKTEyWjy4hJ2acu16rRDWybrP3Vy1XwFu7RlspV12qEmU7Za7NKW2V7XTm3Z\nf+qkzp79c0mJf4f+4i+8kuzTjulqtPPPtyn5Z8dObbmSOf/DTjUl/2606+8d21qqUiMjny798i99\nw7o9NPhXS7/+v/2LMlaU3mM//0i5S0jrld/73aXHfv6Rpcd+/pEl38WL1v12atdMNS4t2addfRcv\nprRPKBRceuznH1n68EO/bdoxU42hUHBpackebXn+/H9eeunF37JuDw3+lXXbTp/JbHUuLdmjLU2+\nixeXXnrxt1JqslNbmtLVaad2zFSLndoyW3vZpS3HxsaWfvmXvmH9uzM2Nrb0yu/9rq3aMVONJru0\nZbJTb//J0qm3/2Rpaclen8mVkuu0Szv6Ll5c+qe/+k9Sbtvx945dVe0QLL/PJ0+bx7rd3tGhV195\nuYwVrRYOh+R0unT2zOm0Q1/K6bvfe0eS9PijX0+5307tmqlGO7Wr0+XU8/v3r7r/p3/3kW3aMVON\nkn3asr39aXk8bdZtwzDkdLok2eszma1Ou7SlWcvxN17T2XffU+cz37Tut1NbSunrtFs7ZqrFLm2Z\nrUY7taXff1GeNo9CobACIwHtdbv13e+9o7NnTtuiHbPVKNmrLU2GYcjv8+nCD9+XZJ/P5ErJddqp\nHfe63TIMQ4HAiDyeNgUCI7b8vWNXVTsEKxwOpZwImB8aOwmMBBQOh6xaVw59sSPaNT8eT5v1j6dh\nGDr44rfUe6hPc3NztmnHTDU6nS7btKXb7bZ+AXQ+882UVYbs9JnMVqdd2lKS3nj9mA4fObpqbped\n2lJKX6ed2jFbLXZpy2w12qktrxqGwqGwDr/yss6e/XM9/9yz8noHbNOO2WqU7NWWpuNvvJYyv8tO\nbZksuU47taPb7VbvoT49/9yzevzRr6fMl7NrW9pJ1faAVIL2jg5d+uzzlNtPPfmEDn/nO3yYN8CO\n7dp/6qSGhwbV2dWtw0eO2nJzzJU1SvZrS4+nTRd++H7ixP6Vl60re3aTrk67tOXZM6flcrrU2dVd\nstdcj0x12qUd16rFLrLVaKe2lBIndRd++L6cTpcMw1DnM9/Uv/ytF7Rly5aS15JJuhq7urpt15aB\nwIjCobDtJ0avrNNO7ej3+TQ8NKjhH30gt9tt+987dlO1PSDmVUiTsWKSsh2ZP1x2rpV2zd/zzz2r\ncCisCz983zqxt1s7pqsxnXK1pd/nsyb7OZ0uHT5y1Go/O7VltjpXKltb+i/K6x3Q449+3Rq++Pij\nX08MKbBTW2apM1m5f76TJddip7ZMlq29ytmWe91u7XW7rRrMXq/W1lbbtGOmGgMjgVWPLffn8sL5\n86vCux0/k+nqTFbOdvT7L6qzq9t6n3sP9SkQGEkZkmWyQ1vaTdUGkMc9bVZXnnTvg2QnZ8+c1sGX\nXrBum+MLk7v17IZ2zY/XOyCny6njJ95KuXpjp3bMVKNkn7YMBEasoQ7JdUj2astsddqlLc+9+54u\nffa59UeSLn32udo7OmzVlpnqDARGbNGOUvb31C5tma1Gu3wmpcT8qeT2Mmv5F88+Z4t2zFZje0eH\nrdpSkoaWe7OT2eUzmWxlnXZqx71ud0rIMP/udrtt2ZZ2U7VDsJLH7rW3P62hoUGd+/4Pyl1Wip79\n++X3X9Tzzz0rp8upq4ahN0+8Ve6ysqJd83N1eXLdyonylz773DbtmK1Gu7Rlz/79OvjiYEod5jAX\nO30ms9Vpl7bMxk5tmYmd2jFbLXZpy2w12qkt3W63evbv1/PPPZs48RsJ6PB3vmObdsxWo2SvtvT7\nfHK5nKvmeNmpLaX0ddqpHXt6DiSGJT/zTev9Ni/UOZ0uW7WlHdUsLS0tlbuIcjIMQ1cNQ542j23n\nVZhjIO1c40q0a2FUQjtK9mlLc/hNujrs1JbZ6rRLW2Zjp7bMxE7tmK0Wu7Rlthrt1JaZ2ssu7bhW\nLXZqy0zs1JaZ2KkdK+Hn246qPoAAAAAAKJ2qnQMCAAAAoPQIIAAAAABKhgACAAAAoGQIIAAAAABK\nhgACAAAAoGQIIAAAAABKhgACAAAAoGQIIAAAAABKhgACAAAAoGQIIAAAAABKhgACAAAAoGQIIAAA\nAABKhgACAAAAoGQIIAAAAABKhgACAAAAoGQIIAAAAABKhgACAAAAoGQIIAAAAABKpr7cBQDAeni9\nA7pqGKvudzpd6uzqltvtLsrr9p86KUl6vudAxtfI5THlYNe6Mr2XUuL97Nm/X06na9Wx4aFBBQIj\nkjJ/T5mee6/bra6u7rTPm8yubQYAlaxmaWlpqdxFAEC+Dr70gvw+X9pjTqdL3/3eO2rv6Cj46z7+\n6NclSefefU/tHR3qP3VSgcCIOru61dNzIO1j7MKudWV7L6XE+/nXf/PjVWGh85lvylgOF72H+nT4\nyNG8ntvpdOnc938gj6ct42vbtc0AoJLRAwKgorV3dKScQHrPn1c4HFL/qZO60PF+wV+v91CfpMQV\ndEkKBEbk9/lSalj5GOTG7Xars6s75b6zZ04rHA7Je/681a5Sot2NpJ6N4aHBtAHElOlzcviVlzX8\now8K+F0AANZCAAFQ0TyetpQTz71ut46/fswammMyDEMXvAMKBEbkdLrk8bSlHdrj9Q7I7/MpHA7J\nvXefOru6M175Th7eEwiMyOsdsHpB0smlhuQhP8NDg/L7LyaGIfUcyHoF3jxJ9/svSpLce/epvaNj\n1Ql9sWvy+3zyegckyRri5PdfVHv702v2IOx1u1eFCL/Pp0BgROFwKOX+4aFBSYlg4ff5ZBiGAoGR\njL0ZKz8nnV3dev65Z2UYxprv20q5fpb8Pp+GhwZlXB1L+1lKN7zL/EytbK98PpcAYHcEEACbiivD\nXIFXX3l51X3DQ4M69/0fWCeOK4fr+JU4mT5+4i3rBPXsmdOSpPb2pxMnl8sBxPy6np4DKY8xTyxz\nrcH8WvPEO/mx2YYBHXzxWymPN2vPNDSpGDWtfL7hoUF5PG3W16znhNkMHquC4vnzkhJBwul0aXho\nUBfOn5fnRObhVMk8njYrvFwaCUg9udWTa5t5vQM6/vox6zHp3o9MnxPzs2S2Vy6fSwCoJKyCBaCi\nhUNh+X0+6495Umee0IXDIb2xfCLY03NAH370U5179z05nS4FAiPW482hVJJ04Yfv69Jnn1tDfvrf\nfjvtayeffPce6tO5d99LX2OONaR7/nPvvpdycpqOGQycTpfOvfteSu0re4KKWVPy85179z0dPnI0\n4+unc9Uw1H/qpPXn4EsvyDAMq4fBNDw0aAWTrqSegKEM7ZOJ2VtiXB3L6fG5tlk4HLI+M72H+nTp\ns89TQsfK3pxs1vO5BAC7I4AAqGhe74AOvvSC9cc84TVP+AIjAeuE7/B3viOn06X2jg4d/s53JKU/\nqb9w/rw1p+Dcu+/pu997Z0M1rqeG3kN9au/oSBlGlelE2dPmSQSD7/9A7R0dq+ZHlKImc3iQ+Xzt\nHR3W43NlGIbOnjlt/fH7fHI6XTp85GhKD4h5Qm72fnQt1xIOhzKGtELItc3MtjBrlxJtZ/4JhcLr\nev1Cfy4BoFwYggWgornd7lWTvXt6DlgnyOaciPaOjpSTWPfefZJknagnD8nxegeseQzJq1utV641\nJFtrediVjzUn3ufa41DomjI9n8fTlnWFq2QrJ6GbV/9ffeVlffjRTyUlQobZ0+F2u63nNod6+X2+\nNee9mNYKaSvl2mbme+Bp86R8fbZJ8pkU83MJAOVCAAFQ0Tq7urOe2JkniuEcrjqfe/c9BQIj1jh8\n8+/DQ4Ma/tEH694HIp8a1sM8SXc6XYleh/anZVwdS5mDUKqaNvJ8Kyehm8PqwuGQ/D6f2js6NJQ0\n/MrsKUk2NDSo4yfeyun1ksNLLor9PmZSrM8lAJQLQ7AAbGrmEKBAYGTVBGrp3smn3+ezehAOHzmq\nCz98P2V51kwb5RWyhvUyr8x72jw6fOSo2js61qy30DW1tz9tPZ/ZE1CMIVFmaHC73dZwMPNPrq8Z\nDod0/PVjVpB5PseehFzbzGwLc3Uu8zWfevIJPf7o17P2CK0MN8X8XAJAudADAmBTSx7CcvDFb6ln\n/36FQ2FrKIs5oTccDunsmdNyOl1pT+pWDqdZKdsysLnWsF7mlfnASMBa2jXTJPJi1dTe0SG32y3D\nMNT70gvq7OpOOQEvhORwcfBQ36ohSM8/92zGYVjpekukxPeZaw9Crm1m7jkSCIzo8CsvW/NywuGQ\nFZyke0Pnjr/xmjq7ulcFG/N73sjnEgDsiB4QAJueuSu6eTKXfMJonqh2dnWr91Cf9Rjzj9Pp0vET\nb2Wc/2B+/fDQoHXyv94a1qtn/3653e57z520aV+2EFDomvq/946cTpc1mdwwjA19b8lzLfz+iymr\nXHWleV7ztZKHaWXi8bTp+Im38p6XkWubHT/xltxut7U6lt/nk9vtVn/SxHFz8rrZXuFQeFWAXe/n\nEgDsrGZpaWmp3EUAQCkEAiPWEBdPmyftyVs4HFJgJGDdzvS4dM/rdDnXHLqUSw3rZc1pyPN5C1lT\ncvt52jzWCXO2/UgqUa5tljzcKt2KYIZhWD0ba200me/nEgDsigACANiwQGBEzz/3rKREL4E5pOjg\ni99KDDNi0zwAwDICCACgIPpPnUw7z6Kn50DOK1MBADY/AggAoGAMw7B2Knc6Xers6maZWABACgII\nAAAAgJJhFSwAAAAAJUMAAQAAAFAyBBAAAAAAJUMAAQAAAFAyBBAAAAAAJUMAAQAAAFAy/z+Iw9pi\n9/M9rQAAAABJRU5ErkJggg==\n", "text/html": [ "
05101520253035404550556065707580Position along RAD locus010002000N variables sites
  • Save as .csv
" ], "text/plain": [ "" ] }, "execution_count": 171, "metadata": {}, "output_type": "execute_result" } ], "source": [ "SNP_position_plot(distvar, distpis)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why the right side but not the left?\n", "I found early on that leaving the cut site attached to the left side of reads improved assemblies by acting as an achor, which then allowed gap openings to arise in the early parts of the read but not to be treated differently, i.e., as terminal gap openings. For some reason it didn't occur to me to similarly anchor the right side of reads. Let's see what happens if I had an invariant anchor to the right side of reads. \n", "\n", "Another point to note, though I don't show the results here, this increase in variation along the length of reads is not observed in simulated data, suggesting it is inherent to real data. \n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " loading Assembly: cli\n", " from saved path: ~/Documents/ipyrad/tests/cli/cli.json\n" ] }, { "data": { "text/plain": [ "(,\n", " )" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/html": [ "
total variable sitesparsimony informative sites051015202530354045505560657075808590Position along RAD loci0102030N variables sites
  • Save as .csv
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "simdata = ip.load_json(\"~/Documents/ipyrad/tests/cli/cli.json\")\n", "\n", "## filter the snps\n", "fsnps = filter_snps(simdata)\n", "\n", "## trim non-data \n", "maxend = np.where(fsnps[:, :, :].sum(axis=0).sum(axis=1) != 0)[0].max()\n", "\n", "## all variables (including autapomorphies)\n", "distvar = np.sum(fsnps[:, :maxend+1, 0].astype(np.int), axis=0)\n", "\n", "## synapomorphies (just pis)\n", "distpis = fsnps[:, :maxend+1, 1].sum(axis=0)\n", "\n", "SNP_position_plot(distvar, distpis)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results when I add a right-side anchor\n", "I now add a five base anchor just before muscle alignment, and then strip it off again after aligning. I have an added check to make sure that the anchor is definitely not left on the read in case muscle did try to stick an opening into the anchor. " ] }, { "cell_type": "code", "execution_count": 172, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.2.1\n" ] } ], "source": [ "## import the new version of ipyrad w/ this update.\n", "import ipyrad as ip\n", "print ip.__version__" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }