mortonjt/ancom_example.ipynb

## ancom_example.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from skbio.stats.composition import ancom\n",
    "from skbio.stats.composition import ancom, closure, clr\n",
    "import biom\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will load the metadata information and OTU tables from memory.  `metadata.txt` in this case contains information about treatment and diseased cohorts, and the fecal samples corresponding to each individual.  `table.biom` contains information about the OTU abundances per fecal sample.\n",
    "\n",
    "There are some scenarios where reading the metadata directly from the tables may cause issues.  So we will take advantage of the additional parameters available in `pandas` to not infer data types, and read all of the variables in as strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "map_fp = 'metadata.txt'\n",
    "map_ = pd.read_csv(map_fp, \n",
    "                   sep='\\t',\n",
    "                   dtype=str, \n",
    "                   encoding = \"ISO-8859-1\",\n",
    "                   na_values=['NA', 'Not applicable', 'not applicable', \"Don't Know\"]\n",
    "                   )\n",
    "map_.set_index('#SampleID', inplace=True) \n",
    "\n",
    "otu_biom = biom.load_table('table.biom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we load the biom table in memory, we'll want to convert it from the biom format into a pandas dataframe.\n",
    "See my [blogpost](http://mortonjt.blogspot.com/2016_07_01_archive.html) about more information about the biom format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def exploding_panda(_bt):\n",
    "    \"\"\"BIOM->Pandas dataframe converter\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    _bt : biom.Table\n",
    "        BIOM table\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    pandas.DataFrame\n",
    "        The BIOM table converted into a DataFrame\n",
    "        object.\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    Based on this answer on SO:\n",
    "    http://stackoverflow.com/a/17819427/379593\n",
    "    \n",
    "    Function written by Yoshiki Vazquez Baeza \n",
    "    (from https://github.com/biocore/biom-format/issues/622)\n",
    "    \"\"\"\n",
    "    m = _bt.matrix_data\n",
    "    data = [pd.SparseSeries(m[i].toarray().ravel()) for i in np.arange(m.shape[0])]\n",
    "    out = pd.SparseDataFrame(data, index=_bt.ids('observation'),\n",
    "                             columns=_bt.ids('sample'))\n",
    "\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "otu_ = exploding_panda(otu_biom).to_dense().T\n",
    "otu_ = otu_.astype(np.int)  # convert the counts into integers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will only consider samples that are either case or control.  The reason why we want to do this is because\n",
    "some of the samples may be blank, so we'll want to exclude those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "map_ = map_.dropna(subset=['case_control'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will perform some filtering and remove the OTUs that show up in less than 10% of the samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# remove all OTUs in less than 10% of samples\n",
    "otu_ = otu_.loc[:, (otu_ > 0).sum(0) > 0.1*len(map_.index)].astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some samples that were fishy, so we'll remove those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Defines the samples to filter\n",
    "keep_samples = (set(map_.index).intersection(set(otu_.index)) - \n",
    "                {'10122.BLANK.1.1A', '10122.BLANK.3.3C', '10122.BLANK.4.12H', \n",
    "                 '10122.BLANK.4.4E'})\n",
    "map_ = map_.loc[keep_samples]\n",
    "otu_ = otu_.loc[map_.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll run ANCOM.  By default, this will run an ANOVA on all pairs of log ratios. For more information about the parameters in ANCOM check out [scikit-bio](http://scikit-bio.org/docs/0.4.2/generated/generated/skbio.stats.composition.ancom.html).  For more information about the caveats behind ANCOM check out my [blogpost](http://mortonjt.blogspot.com/2016_06_01_archive.html) on it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Don't forget to add pseudocounts.\n",
    "res, pcts = ancom(otu_ + 1, map_.case_control)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's take a look at the significant hits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Reject null hypothesis</th>\n",
       "      <th>W</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4465907</th>\n",
       "      <td>True</td>\n",
       "      <td>651</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Reject null hypothesis    W\n",
       "4465907                   True  651"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ids = res.loc[res['Reject null hypothesis']].index\n",
    "res.loc[res['Reject null hypothesis']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we'll visualize some diagnostic plots, namely volcano plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myu_pval = lambda x: -np.log(mannwhitneyu(x[map_.case_control=='Case'], x[map_.case_control=='Control']))[1]\n",
    "fchange = lambda x: x[map_.case_control=='Case'].mean() - x[map_.case_control=='Control'].mean()\n",
    "\n",
    "#pvals = otu_.apply(myu_pval, axis=0)\n",
    "pvals = res.W\n",
    "clrotu_ = otu_.apply(lambda x: clr(x+1), axis=1)\n",
    "fold_change = clrotu_.apply(fchange, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAJiCAYAAABKE490AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xt43HWd9//nJ03bNJNBzi20HOxBhHXd2+Jy0ChFBKRF\nDoVtgIKieLt7re7Swt60dHXpLbgIqy78AH/+lvjbmy5VBrEIaLBUpUBcES0rHkBIjZyaUlRWOk3o\niXzuP2bSJs3k2Elm5pvn47pykXznO9/5zIFrXv0c3p8QY0SSJEmVr6rUDZAkSVJxGOwkSZISwmAn\nSZKUEAY7SZKkhDDYSZIkJYTBTpIkKSHKKtiFEN4WQvivEMKT+f++HkL4+xDCfiGEh0IIz4YQVocQ\n3tLtPleHEFpCCM+EEE4rZfslSZJKKZRrHbsQQhXwMnA88GngjzHGG0MIS4D9YoxLQwjHACuBvwSm\nAd8HZsVyfVKSJEkjqKx67PbwQeC3McaXgLOBO/LH7wDOyf9+FnBXjHFnjPF5oAU4brQbKkmSVA7K\nOdg1AF/P/z45xrgJIMb4CnBw/vhU4KVu99mQPyZJkjTmlGWwCyGMJ9cb9838oT2HVh1qlSRJ2kN1\nqRvQhzOAdTHGP+T/3hRCmBxj3BRCmAK8mj++ATis2/2m5Y/1EkIwDEqSpIoRYwxDvU9Z9tgBFwLf\n6Pb3/cCl+d8/CtzX7fgFIYQJIYS3AjOBJ/q6aIzRnwr8ueaaa0reBn98/8bqj+9f5f743lX2z3CV\nXY9dCKGW3MKJT3Y7fANwdwjh48ALwAKAGOPTIYS7gaeBHcDfxr15NSRJkipY2QW7GGMHcNAex14j\nF/YKnX89cP0oNE2SJKmsletQrLTLnDlzSt0E7QXfv8rm+1e5fO/GprItUFxsIQRHaSVJUkUIIRCH\nsXii7IZiJUkqtSOPPJIXXnih1M3QGHDEEUfw/PPPF+169thJkrSHfG9JqZuhMaCvz9pwe+ycYydJ\nkpQQBjtJkqSEMNhJkiQlhMFOkiSVjeuvv55PfvKTA584CO94xzt49NFHd/39sY99jP33358TTjiB\n5uZmjj766KI8Tjkx2EmSVKHmzJnD/vvvz44dO3ocv/TSS6mqquJnP/vZrmO//e1vqarq+bW/evVq\nTjrpJPbZZx8mT57MySefzAMPPLDr9g0bNnDxxRdz4IEHkk6nOeGEE/jud7/b4xpVVVVMmTKFzs7O\nXcd27tzJwQcfzLhx44b8nK6++mr+7d/+bcj3K+RXv/oV73//+wFobm7mBz/4AW1tbTz++OPU19fz\nzDPPFOVx+vLII49w2GGHDXxiERnsJEkqks7OTr73ve+x+FOf4qrFi3n88cdHbHXtCy+8QHNzM1VV\nVdx///09bgshcMABB/CZz3ym1/Eu99xzDwsWLODSSy9lw4YNbNq0ic997nN85zvfAeC///u/qa+v\np6amhmeeeYY//OEPLFq0iIsuuohVq1b1uO5+++3Hgw8+uOvvBx98kP3337/YT3mvPP/88xx55JHU\n1NSM2mPGGHu85qP2oGPhJ/dUJUka2HC+M7LZbDzp3e+O/6OuLn4B4vKqqjg9lYoXz58fd+7cWfQ2\nfu5zn4v19fXxyiuvjGeeeWaP2y699NJ45ZVXxkMOOSQ++uijMcYY169fH6uqqnadc/jhh8cvfelL\nfV7/M5/5TPzzP//zXsdvuOGGeMQRR+z6O4QQP//5z8e/+qu/2nXs/PPPj//8z//c4/H29IUvfCFO\nnTo1ptPp+Pa3vz3+8Ic/jDHGuHz58njxxRfvOu+OO+6IRxxxRDzwwAPjtddeG4888sj4gx/8YNe5\nCxYsiB/5yEdiOp2O73jHO+K6det23bfr3K997WuxpqYmVldXx3Q6HZcvXx7Xrl0bp02btuvcl156\nKc6fPz8edNBB8cADD4x/93d/F2OM8be//W38wAc+EA844IB40EEHxYULF8bXX3+9x2N88YtfjO98\n5zvjvvvuGxsaGuK2bdtie3t7nDRpUhw3blysq6uL6XQ6bty4sdfr0NdnLX98yHnHHjtJkopg6aJF\nHP7LX7JuyxaWANd0dvKr9nZe/N73+Mqttxb98VasWMHFF1/MRRddxOrVq/n973/f4/ba2lqWLVvG\nsmXLet33N7/5DS+//DLnnXden9f//ve/X/D2BQsW8OKLL9LS0gLkegHPOeccHn30UTZv3syf/vQn\nmpubOfvss/u89nPPPcdtt93GunXr2Lx5M6tXr+bII4/cdXtXL9fTTz/Npz71Kb7xjW+wceNGXn/9\nddra2npc64EHHuCiiy7i9ddf58Mf/jCf+tSnej3exz/+cb761a9y4oknsnnzZq655poej9PZ2cmZ\nZ57JW9/6Vl588UU2bNjABRdcAOQ6wJYtW8Yrr7zCM888w8svv8zy5ct7XP+b3/wmDz30EL/73e/4\nxS9+wf/5P/+H2tpaHnzwQQ499FCy2SybN29mypQpfb4mxWKwkyRpL23dupU7V67kC9u29fhinQRc\n19HB//flLxf18Zqbm3nxxRdZsGABs2fPZubMmXz961/vdd4nP/lJXnzxRVavXt3j+GuvvQbAIYcc\n0udj/OEPfyh4e9exP/zhD0Au+NTU1HDWWWdx1113kclkOOuss5g4cWKf1x43bhzbt2/nV7/6FTt3\n7uTwww/nrW99a6/zvvWtb3HWWWdx4oknUl1dzec+97le59TX13P66acTQuCSSy7hF7/4RZ+P25ef\n/OQnbNy4kRtvvJGamhomTJjAe97zHgBmzJjBKaecQnV1NQcccACLFy/mkUce6XH/yy+/nMmTJ7Pv\nvvvy4Q9/mJ///OdDbkOxGOwkSdpL//3f/83EEDi0wG3vBF7YtKmoj7dixQpOO+009ttvPwAuvPBC\n7rjjjl7nTZgwgc9+9rN89rOf7XH8gAMOAGDjxo19PsaBBx5Y8PauYwcddFCP45dccgkrVqzgP/7j\nP/jIRz7Sb/tnzJjBTTfdxPLly5k8eTIXXXQRr7zySq/z2traeiw+mDRp0q62d+neC1ZbW8vWrVt7\nLOQYjJdffpkjjjii1+ISgFdffZULL7yQadOmse+++3LxxRfvCrVdJk+e3KMNW7ZsGdLjF5PBTpKk\nvXTAAQfwZlUVrQVu+wlw1BFHFO2xtm7dyt13380jjzzCIYccwiGHHMJNN93EU089xS9/+cte53/s\nYx/jT3/6U48FD0cddRSHHXYY3/rWt/p8nA9+8IO9FkkAZDIZDj/8cGbOnNnj+Pve9z42btzIq6++\nynvf+94Bn8cFF1zAY489tmtP3iVLlvQ655BDDuHll1/e9fcbb7zBH//4xwGvPVSHHXYYL774YsFA\nuGzZMqqqqvj1r3/Nn/70J+68885BL4gZ9YUTGOwkSdprEyZM4JN/8zdcPmkS27odfw1Ykkrx6aVL\ni/ZY9957L9XV1TzzzDM89dRTPPXUUzzzzDPU19ezYsWKXuePGzeO5cuXc8MNN/Q4/qUvfYlrr72W\nO+64g2w2S4yR5uZm/vqv/xqAxYsX8/rrr3PZZZexadMmtm3bxje+8Q2uv/56vvjFLxZs23e+8x3u\nu+++XX/3FYCee+45Hn74YbZv386ECROYNGlSwd6y888/nwceeIDHH3+cHTt29JrbVshgQ1d3xx13\nHIcccghLly6lo6ODbdu28Z//+Z8AZLNZ6urqSKfTbNiwgX/5l38Z9HUnT57MH//4RzZv3jzkNg2X\nwU6SpCJY/s//TOqDH2RGbS2Xjx/PJydO5G01NZz2iU/w0UsvLdrjrFixgo9//ONMnTqVgw8+eNfP\npz/9aVauXFmw1+nCCy/kkEMO6dGDdN5555HJZPja177G1KlTmTJlCv/0T//EOeecA8D+++9Pc3Mz\nb7zxBscccwwHHnggN910E3feeSfnn3/+rut0v+bRRx/do+hvXz1W27ZtY+nSpRx00EEceuih/P73\nv+f666/vdd4xxxzDLbfcQkNDA4ceeij77LMPBx98cL/z97o/5mB7zKqqqnjggQdoaWnh8MMP57DD\nDuPuu+8G4JprrmHdunW75s/tuaCkv8c46qijuPDCC5k+fTr7779/weHmYgvDSbaVKIQQx8pzlSTt\nnRDCsOvPPfXUUzz00EOMHz+es88+u+CiAA1Pe3s7++67L+vXr+eIIg5vl1Jfn7X88SGP5RrsJEna\nw94EOxXXd77zHU455RQ6Ozu58sor+elPf8q6detK3ayiKXawcyhWkiSVrfvuu49DDz2UadOm8dvf\n/pa77rqr1E0qa/bYSZK0B3vsNFrssZMkSVJBBjtJkqSEMNhJkiQlRHWpGyBJUrk54ogjSrJrgMae\nYpdtcfGEJElSmXHxhJTX1tbGRz/6Ud5z/Il89KMfpa2trdRNkiRpVBjslChf+cpXmDF1GhtXrGTe\nE4+zccVKZkydxle+8pVSN02SpBHnUKwSo62tjRlTp3E/kVO7HV8DnEXgdxvbmDJlSqmaJ0nSoDkU\nqzHv6quv5n1U9Qh1AKcC76OKpUuXlqJZkiSNGoOdEqPlN89xEm8WvO39vEnLs8+NcoskSRpdBjsl\nxqy3v41HGFfwtkcZx6yj3jbKLZIkaXQ5x06J4Rw7SVJSOMdOY96hhx7Kl267lbMInMY4rgNOYxxn\nEfjSbbca6iRJiWePnRLnlVdeYenSpbQ8+xyzjnobX/jCFwx1kqSKMtweO4OdJElSmXEoVpIkaYwz\n2EmSJCWEwU6SJCkhDHaSJEkJYbCTJElKCIOdJElSQhjsJEmSEsJgJ0mSlBAGO0mSpIQw2EmSJCWE\nwU6SJCkhDHaSJEkJYbCTJElKCIOdJElSQhjsJEmSEsJgJ0mSlBAGO0mSpIQw2EmSJCWEwU6SJCkh\nDHaSJEkJYbCTJElKCIOdJElSQhjsJEmSEsJgJ0mSlBAGO0mSpIQou2AXQnhLCOGbIYRnQgi/DiEc\nH0LYL4TwUAjh2RDC6hDCW7qdf3UIoSV//mmlbLskSVIplV2wA24GmmKMRwN/AfwGWAp8P8Z4FPBD\n4GqAEMIxwALgaOAM4CshhFCSVkuSJJVYWQW7EMI+wPtijP8OEGPcGWN8HTgbuCN/2h3AOfnfzwLu\nyp/3PNACHDe6rZYkSSoPZRXsgLcCfwgh/HsI4ckQwr+FEGqByTHGTQAxxleAg/PnTwVe6nb/Dflj\nkiRJY065BbtqYDZwW4xxNtBObhg27nHenn9LkiSNedWlbsAeXgZeijH+LP/3t8gFu00hhMkxxk0h\nhCnAq/nbNwCHdbv/tPyxgpYvX77r9zlz5jBnzpzitVySJGmY1q5dy9q1a/f6OiHG8ur8CiE8AvzP\nGONzIYRrgNr8Ta/FGG8IISwB9osxLs0vnlgJHE9uCHYNMCsWeFIhhEKHJUmSyk4IgRjjkBeElluP\nHcDfAytDCOOBVuBjwDjg7hDCx4EXyK2EJcb4dAjhbuBpYAfwt6Y3SZI0VpVdj91IscdOkiRViuH2\n2JXb4glJkiQNk8FOkiQpIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhLCYCdJkpQQBjtJkqSEMNhJ\nkiQlhMFOkiQpIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhLCYCdJkpQQBjtJkqSEMNhJkiQlRHWp\nGyB1l81myWQytLS0MmvWdBoaGkin06VuliRJFSHEGEvdhlERQohj5blWqubmZubOnU9nZz3t7ceS\nSq2jqqqZpqZV1NfXl7p5kiSNmhACMcYw5PuNlbBjsCtv2WyWqVNnkM2uBE7tdssa0umFtLW1UldX\nV6rmSZI0qoYb7Jxjp7KQyWTo7KynZ6gDOJXOznoymUwpmiVJUkUx2KkstLS00t5+bMHb2ttns359\n6yi3SJKkymOwU1mYNWs6qdS6grelUk8yc+b0UW6RJEmVxzl2KgvOsZMkabfhzrGz3InKQjqdpqlp\nVbdVsbNJpZ7ctSrWUCdJ0sDssVNZ2bJlC5lMhvXrW5k5M1fHzlAnSRprLHcyAIOdJEmqFJY7kSRJ\nGuMMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIkJYTBTpIkKSEMdpIk\nSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIkJYTBTpIkKSEMdpIkSQlRXeoG\nSGNJNpslk8nQ0tLKrFnTaWhoIJ1Ol7pZkqSECDHGUrdhVIQQ4lh5ripPzc3NzJ07n87OetrbjyWV\nWkdVVTNNTauor68vdfMkSWUkhECMMQz5fmMl7BjsVErZbJapU2eQza4ETu12yxrS6YW0tbVSV1dX\nquZJksrMcIOdc+ykUZDJZOjsrKdnqAM4lc7OejKZTCmaJUlKGIOdNApaWlppbz+24G3t7bNZv751\nlFskSUoig500CmbNmk4qta7gbanUk8ycOX2UWyRJSiLn2EmjwDl2kqShGO4cO8udSKMgnU7T1LSq\n26rY2aRST+5aFWuokyQVgz120ijasmULmUyG9etbmTkzV8fOUCdJ2pPlTgZgsJMkSZXCcieSJElj\nnMFOkiQpIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhKi7IJdCOH5EMJTIYT/CiE8kT+2XwjhoRDC\nsyGE1SGEt3Q7/+oQQksI4ZkQwmmla7kkSVJplV2wAzqBOTHGd8UYj8sfWwp8P8Z4FPBD4GqAEMIx\nwALgaOAM4CshhCHXfJEkSUqCcgx2gd7tOhu4I//7HcA5+d/PAu6KMe6MMT4PtADHIUmSNAaVY7CL\nwJoQwk9DCJ/IH5scY9wEEGN8BTg4f3wq8FK3+27IH5MkSRpzqkvdgALeG2PcGEI4CHgohPAsubDX\n3bD2Blu+fPmu3+fMmcOcOXOG20ZJkqSiWbt2LWvXrt3r65T1XrEhhGuALcAnyM272xRCmAI8HGM8\nOoSwFIgxxhvy538PuCbG+JMC13KvWEmSVBESsVdsCKE2hFCX/z0FnAb8ErgfuDR/2keB+/K/3w9c\nEEKYEEJ4KzATeGJUGy1JklQmym0odjJwbwghkmvbyhjjQyGEnwF3hxA+DrxAbiUsMcanQwh3A08D\nO4C/tVtOkiSNVWU9FFtMDsVKkqRKkYihWEmSJA2fwU6SJCkhDHaSJEkJYbCTJElKCIOdJElSQhjs\nJEmSEsJgJ0mSlBAGO0mSpIQw2EmSJCWEwU6SJCkhDHaSJEkJYbCTJElKCIOdJElSQhjsJEmSEsJg\nJ0mSlBAGO0mSpIQw2EmSJCWEwU6SJCkhDHaSJEkJYbCTJElKCIOdJElSQhjsJEmSEsJgJ0mSlBAG\nO0mSpIQw2EmSJCWEwU6SJCkhDHaSJEkJYbCTJElKiOpSN0BSecpms2QyGVpaWpk1azoNDQ2k0+lS\nN0uS1I8QYyx1G0ZFCCGOlecq7a3m5mbmzp1PZ2c97e3Hkkqto6qqmaamVdTX15e6eZKUeCEEYoxh\nyPcbK2HHYCcNTjabZerUGWSzK4FTu92yhnR6IW1trdTV1ZWqeZI0Jgw32DnHTlIPmUyGzs56eoY6\ngFPp7Kwnk8mUolmSpEEw2EnqoaWllfb2Ywve1t4+m/XrW0e5RZKkwTLYSeph1qzppFLrCt6WSj3J\nzJnTR7lFkqTBco6dpB6cYydJpTfcOXaWO5HUQzqdpqlpVbdVsbNJpZ7ctSrWUCdJ5cseO0kFbdmy\nhUwmw/r1rcycmatjZ6iTpNFhuZMBGOwkSVKlsNyJJEnSGGewkyRJSgiDnSRJUkIY7CRJkhLCYCdJ\nkpQQBjtJkqSEMNhJkiQlhMFOkiQpIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhLCYCdJkpQQBjtJ\nkqSEMNhJkiQlhMFOkiQpIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhLCYCdJkpQQBjtJkqSEqC51\nAwoJIVQBPwNejjGeFULYD8gARwDPAwtijK/nz70a+DiwE7g8xvhQaVotDU02myWTydDS0sqsWdNp\naGggnU6XulmSpAoWYoylbkMvIYTFwLHAPvlgdwPwxxjjjSGEJcB+McalIYRjgJXAXwLTgO8Ds2KB\nJxVCKHRYKonm5mbmzp1PZ2c97e3Hkkqto6qqmaamVdTX15e6eZKkEgshEGMMQ75fuYWdEMI04N+B\nzwNX5IPdb4CTYoybQghTgLUxxreHEJYCMcZ4Q/6+DwLLY4w/KXBdg53KQjabZerUGWSzK4FTu92y\nhnR6IW1trdTV1ZWqeZKkMjDcYFeOc+z+FfhfQPcUNjnGuAkgxvgKcHD++FTgpW7nbcgfk8pWJpOh\ns7OenqEO4FQ6O+vJZDKlaJYkKQHKKtiFEOYBm2KMPwf6S6l2valitbS00t5+bMHb2ttns3596yi3\nSJKUFOW2eOK9wFkhhLnAJCAdQvgP4JUQwuRuQ7Gv5s/fABzW7f7T8scKWr58+a7f58yZw5w5c4rb\nemkQZs2aTirVRHt779tSqSeZOXPe6DdKklRSa9euZe3atXt9nbKbY9clhHAScGV+jt2N5BZP3NDH\n4onjyQ3BrsHFEypzzrGTJA1kuHPsyq3Hri9fAO4OIXwceAFYABBjfDqEcDfwNLAD+FvTm8pdOp2m\nqWlVt1Wxs0mlnty1KtZQJ0karrLtsSs2e+xUbrZs2UImk2H9+lZmzszVsTPUSZIgQeVORorBTpIk\nVYoklTuRJEnSMBjsJEmSEsJgJ0mSlBAGO0mSpIQw2EmSJCWEwU6SJCkhDHaSJEkJYbCTJElKCIOd\nJElSQlTKXrGSRkA2myWTydDS0sqsWbltzdLpdKmbJUkaJrcUk8ao5uZm5s6dT2dnPe3tx5JKraOq\nqpmmplXU19eXunmSNKa5V+wADHbSbtlslqlTZ5DNrgRO7XbLGtLphbS1tVJXV1eq5knSmOdesZIG\nLZPJ0NlZT89QB3AqnZ31ZDKZUjRLkrSXDHbSGNTS0kp7+7EFb2tvn8369a2j3CJJUjEY7KQxaNas\n6aRS6wrelko9ycyZ00e5RZKkYnCOnTQGOcdOksrbcOfYWe5EGoPS6TRNTau6rYqdTSr15K5VsYY6\nSapM9thJY9iWLVvIZDKsX9/KzJm5OnaGOkkqPcudDMBgJ0mSKoXlTiRJksY4g50kSVJCGOwkSZIS\nwmAnSZKUEAY7SZKkhDDYSZIkJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKU\nEAY7SZKkhDDYSZIkJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKk\nhDDYSZIkJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIk\nJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIkJYTBTpIk\nKSEMdpIkSQkxYLALIZwcQpgwGo2RJEnS8A2mx+4HwJ9CCD8IIXwmhPCeEMK4kWhMCGFiCOEnIYT/\nCiH8MoRwTf74fiGEh0IIz4YQVocQ3tLtPleHEFpCCM+EEE4biXZJkiRVghBj7P+EEGYAHwDm5H8O\nAdqBZuCHwMPAujjQhQbboBBqY4wd+fD4I+DvgfOAP8YYbwwhLAH2izEuDSEcA6wE/hKYBnwfmFWo\nLSGEYjVRkiRpRIUQiDGGod5vwB67GONvY4y3xxgXxhinAscAVwGvA1cCPwH+ONQH7ufxOvK/TgSq\ngQicDdyRP34HcE7+97OAu2KMO2OMzwMtwHHFaoskSVIlqR7qHWKMvwkhvAa8Ri7cXQDUFatBIYQq\nYB0wA7gtxvjTEMLkGOOm/OO/EkI4OH/6VODH3e6+IX9MkiRpzBlUsAshHEBuGPZkcsOy08mFr0eA\nBeSGZYsixtgJvCuEsA9wbwjhz8j12vU4rViPJ0mSlBQDBrsQwi+AWcDPyAW5y4EfdRsyHRExxs0h\nhLXAh4BNXb12IYQpwKv50zYAh3W727T8sYKWL1++6/c5c+YwZ86cIrdakiRp6NauXcvatWv3+jqD\nWTzRAfyJ3MKEtcDaGGPrXj9y4cc6ENgRY3w9hDAJWA18ATgJeC3GeEMfiyeOJzcEuwYXT0iSpAo3\n3MUTgxmKfQu54HQycDFwWwjhVXK9d2spbtA7BLgjP8+uCsjEGJtCCI8Dd4cQPg68QG74lxjj0yGE\nu4GngR3A35reJEnSWDVgj12vO4QwETiR3XPujgM2xRiPLHbjiskeO0mSVClGrNxJAZ3dfiIQ6DnP\nTZIkSSUwmMUT1eR65U7O/5wITCI3JPow8LX8fyVJklRCg1k8sYVckNtILsA9DPwwXxC4YjgUK0mS\nKsVILp64Ang4xtgy9GZJkiRptAx58USlssdOkiRVipHssZMEZLNZMpkMLS2tzJo1nYaGBtLpdKmb\nJUnSLvbYSYPQ3NzM3Lnz6eysp739WFKpdVRVNdPUtIr6+vpSN0+SlDDD7bEz2EkDyGazTJ06g2x2\nJXBqt1vWkE4vpK2tlbq6ulI1T5KUQKNZx04aUzKZDJ2d9fQMdQCn0tlZTyaTKUWzJEnqxWAnDaCl\npZX29mML3tbePpv160dk62RJkobMYCcNYNas6aRS6wrelko9ycyZ00e5RZIkFeYcO2kAzrGTJI02\ny51IIySdTtPUtKrbqtjZpFJP7loVa6iTJJULe+ykQdqyZQuZTIb161uZOTNXx65SQ501+SSpvFnu\nZAAGOynHmnySVP4MdgMw2EnOF5RU/hxRyLGOnaQBWZNPUjlrbm5m6tQZLFrUxI03pli0qImpU2fQ\n3Nxc6qZVDBdPSGOINfkklatsNsvcufN7jCi0twOsYe7c+Y4oDJI9dtIYYk0+SeXKEYXiMNhJY0hD\nQwNVVc3Amj1uWUNVVTMNDQ2laJYkOaJQJAY7aQzpqsmXTi8klZoPXEcqNZ90eqE1+SSVlCMKxeGq\nWGkMSlJNPknJ4Kr9nix3MgCDnSRJ5a1nnc2eu/yMtTqbBrsBGOwkSSp/jijkGOwGYLCTJEmVwgLF\nkiRJY5zDn8UNAAAgAElEQVTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7\nSZKkhDDYSZIkJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDY\nSZIkJYTBTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIkJYTB\nTpIkKSEMdpIkSQlhsJMkSUoIg50kSVJCGOwkSZISwmAnSZKUEAY7SZKkhDDYSZIkJUR1qRsgSVKl\nymazZDIZWlpamTVrOg0NDaTT6VI3S2NYiDGWug2jIoQQx8pzlSSNvObmZubOnU9nZz3t7ceSSq2j\nqqqZpqZV1NfXl7p5qnAhBGKMYcj3Gythx2AnSSqWbDbL1KkzyGZXAqd2u2UN6fRC2tpaqaurK1Xz\nlADDDXZlNccuhDAthPDDEMKvQwi/DCH8ff74fiGEh0IIz4YQVocQ3tLtPleHEFpCCM+EEE4rXesl\nSWNFJpOhs7OenqEO4FQ6O+vJZDKlaJZUdnPsdgJXxBh/HkKoA9aFEB4CPgZ8P8Z4YwhhCXA1sDSE\ncAywADgamAZ8P4Qwy645SdJIamlppb392IK3tbfP5umnf0NjY6Nz7zTqyqrHLsb4Sozx5/nftwDP\nkAtsZwN35E+7Azgn//tZwF0xxp0xxueBFuC4UW20JGnMmTVrOqnUuoK31dT8lNtu+yqLFjVx440p\nFi1qYurUGTQ3N49yKzUWle0cuxDCkcBa4B3ASzHG/brd9lqMcf8Qwi3Aj2OMX88fbwSaYoyrClzP\njjxJUlH0N8cu1xdxF7m+h93HnXunoUjEHLsu+WHYe4DL8z13eyYyE5okqWTS6TRNTatIpxeSSs0H\nriOVmk9NzQJqat5Nz1AHzr3TaCm3OXaEEKrJhbr/iDHelz+8KYQwOca4KYQwBXg1f3wDcFi3u0/L\nHyto+fLlu36fM2cOc+bMKWLLJUljSX19PW1trWQyGdavb2XmzHn86lfTuemm/Qqe394+m/XrW0e5\nlaoUa9euZe3atXt9nbIbig0hrAD+EGO8otuxG4DXYow35BdP7Bdj7Fo8sRI4HphKrg+84OIJh2Il\nSSOtsbGRRYuaaG/vNSOIVGo+N988j8suu6wELVOlSUQduxDCe4FHgV+SG26NwDLgCeBucr1zLwAL\nYox/yt/nauAyYAe5oduH+ri2wU6SNKKsb6diSUSwG0kGO0nSaOi5I8VsUqkn3ZFCQ2awG4DBTpI0\nWrZs2dJt7l2ujp09dRoKg90ADHaSJKlSJKrciSRJkobOYCdJkpQQBjtJkqSEMNhJkiQlhMFOkiQp\nIQx2kiRJCWGwkyRJSgiDnSRJUkIY7CRJkhKiutQNkCRJg5PNZslkMrS0tDJrVm6rsnQ6XepmqYy4\npZgkScNUKGgBIxK+mpubmTt3Pp2d9bS3H0sqtY6qqmaamlZRX1+/19cvBoNn8bhX7AAMdpKkYioU\ntDo7HyaEKkI4qajhK5vNMnXqDLLZlcCp3W5ZQzq9kLa2Vurq6vb6Oe2NSgielcRgNwCDnSSpWAoH\nrSxwJHAXxQ5fjY2NLFrURHv7ql63pVLzufnmeVx22WXDunYxVELwrDTDDXYunpAkaYgymQydnfX0\nDDEZ4KQ9jgGcSmdnPZlMZtiP19LSSnv7sQVva2+fzfr1rcO+djEUfj2gGM9dQ2OwkyRpiAoHrVZg\nZMLXrFnTSaXWFbwtlXqSmTOnD/vaxVDuwXMsMdhJkjREhYPWdGBkwldDQwNVVc3Amj1uWUNVVfOu\nRRulUu7Bcyxxjp0kSUM02nPsYM/FCbNJpZ4sm8UJzrErPhdPDMBgJ0kqpkJBq+eq2OKHry1btpDJ\nZFi/vpWZM3PlRMolMJVz8KxEBrsBGOwkScVWKGgBZRu+Rlo5B89KY7AbgMFOkiRVCsudSJIkjXHu\nFSup4rhtkSQV5lCspIritkWSxgLn2A3AYCdVvnIqqWCvoaSR5Bw7SYlXLtsWNTc3M3XqDBYtauLG\nG1MsWtTE1KkzaG5uHpXHl6S+OMdOUsUoh22Lstksc+fO79Fr2N4OsIa5c+dbiFVSSdljJ6lilMO2\nReXSayhJhRjsJFWMctgvsxx6DZUM2WyWxsZGlixZRmNjI9lsttRNUgIY7CRVjHQ6TVPTKtLphaRS\n84HrSKXmk04vpKlp1agMgZZDr6Eqn/M0NVJcFSup4pRy26JyWpmryuRnSIMx3FWxLp6QNCTlUOaj\nrq6Oyy67bFQfs0tXr2Ffm537hayBDGaeZqk+36p8BjtJg9a7OHATV1yxbMwVB66vr6etrbVbr+E8\nGhpWGOo0KM7T1Egy2EkaFMt89FTKXkNVttw8zab8/z895eZpzhv9RikxXDwhaVAs8yEVRzms7lZy\nGewkDSibzXLPPatob38NaAR6lmVw+EgavHJY3a3kcihWUr+65tVt3348ud66JmAZsArIzatz+Ega\nGudpaqRY7kRSn/orywALgVbgx5ZokKQis9yJpKLrb14dHMeECe9n4sSXHT6SpDJhsJPUp/7KMsAJ\nnHLK49x996OJDHXlUK9PkobKxROS+jTQ9lnnnXduIkOd2z1JqlTOsZPUp7G49dFYfM6Sys9w59jZ\nYyepT2OxLIP1+iRVMufYSerXWCvL4HZPkiqZwU7SgMbS9llu9ySpkjnHTpK6cY6dpHJgHTtJZaHS\ny4R0zSucO3c+nZ31tLfPJpV6kqqq5sTOK5SUHPbYSSqaru3HcoHoWFKpdbsCUX19fambNyRbtmzp\nNq8wF1ANdZJGy3B77Ax2koqiXIcwK70HUdLYZLkTSSVVjmVCLDQsaaxxjp2koii3MiHZbJa5c+f3\n6EHMrXRdw9y5810EocSxd1pgj52kIhlo+7GZM6cPeI1sNktjYyNLliyjsbGRbDY77PaUYw+iNFLs\nnVYXg52komhoaKCqqhlYs8cta6iqaqahoaHf+xf7i6ncehClkdK9d7q9fRXwj7S3ryKbXcncufPZ\nsmVLqZuoUeRQrKSi2JsyISMxbGqhYRVTOQ9zDqZ3eqwUGJc9dpKKqGv7sZtvnsfSpW9w883zaGtr\nHbDUyUgMm+5tD6LUpdyHOe2dVnf22EkqquFsPzYSX0wWGlYxVMIiHHun1Z09dpJKrhgLLwoZbg+i\n1KUSFuHYO63u7LGTVHINDQ1cccUycl9MPYsb576YVgz72sPpQZS6VMIwp73T6s5gJ6nk/GJSuaqU\nYc6u3und2+DNo6Fhhf/vjEFuKSapbLg/q8pNuW6Vp+RLzF6xIYSvAWcCm2KM78wf2w/IAEcAzwML\nYoyv52+7Gvg4sBO4PMb4UB/XNdhJJVTO5SKk/jQ3N/fZm+x8TY2UJAW7emALsKJbsLsB+GOM8cYQ\nwhJgvxjj0hDCMcBK4C+BacD3gVmFEpzBTiqdnl+Mx5JKrfOLURXF3mSNtsQEO4AQwhHAA92C3W+A\nk2KMm0IIU4C1Mca3hxCWAjHGeEP+vAeB5THGnxS4psFOKoGkD2XZE6mh8POiwRpusKuUcicHxxg3\nAcQYXwEOzh+fCrzU7bwN+WOSykQllIsYqq49bRcuvISDDz6cyy//blkWrlV5KfdCx0qGSl0VO6yu\nt+XLl+/6fc6cOcyZM6dIzZHUl0ooF9HdQD0qXcPKb755PB0djwHfpFwL16p8ZLNZzjjjHLZs+QhQ\nA0ymvf0O4HE/LwJg7dq1rF27dq+vUynBblMIYXK3odhX88c3AId1O29a/lhB3YOdpNFRKeUioNBc\nwCauuGLZrrmAPXcheAEYj/tzajCuu+46tmzZRm7937FAE7AMWOXnRUDvDqf//b//97CuU65DsSH/\n0+V+4NL87x8F7ut2/IIQwoQQwluBmcATo9VISQOrlKr43UNbe/sq4B9pb19FNruSuXPn75o8v3tY\nuZXcF3Rv5dgTqdLJZrPcdNNXgVX5n3/M/3clMJ/29j/z86KiKbtgF0L4OvCfwNtCCC+GED4GfAE4\nNYTwLHBK/m9ijE8DdwNPk/vnz9+6QkIqL13Fh9PphaRS84HrSKXmk04vLKviw4OZC9hzWHk6UPxt\n0JQ8mUyGEN5Poc8W1DNhwnf9vKhoym4oNsZ4UR83fbCP868Hrh+5FknaW5VQFX8wcwF7Dis3kBtK\nK/42aEqWlpZWtm07oY9bZxPjQ2XTc63KV3bBTlIylfuerYOZC7hgwYI99rRdBcwHTgT+Ml+f70dl\n1ROp0uvvswWPccUVn/bzoqIpyzp2I8E6dpL6M9h6e3vuQlBb+wRvvvkw559/LieffJKFa9VLf5+t\nurqFbNzoilj1lqgCxSPBYCdpIIPdOspdCCrXUAsEF6ug8GA+WxYvVncGuwEY7CQNhqEtuYa6tV2x\nt8Lr77Pltnvak8FuAAY7SeXAXpnSGOrWdqO5FV7St93T8CR9SzGprHVtMbVkyTIaGxvJZrOlbpLK\nkFtKlc5Qt7Ybza3wkrjtnkrHVbHSXhpotwJVpmL3rPXctcItyEbbULe2G82t8Cpt2z2VN3vspL0w\nmN0KVHlGomfNXpnSypUcGXxB6aGeP5ptk/pjsJP2gl/WyTNSYd1emdIa6tZ2o7kVXqVsu6fK4FCs\ntBf8sq5shYZbBxPWh1NoeTAFkDVyura266vkyJ7D4EM9fzTbJvXHYCftBb+sK1dfcyM//OEPjUhY\nb2ho2GPXii5uQVYMg5kTOdSt7UZzK7xK2HZPlcFyJ9JesExBZervfZs4cT7V1R+kvf3eXvdLpeZz\n883zhr012mALIGtorAGnJLKO3QAMdhopfllXnsbGRhYtasrPoeuptvbDdHY2s3Xr3ewZ+mpqFvA3\nf3MZf/Znbx/2KlkLIBfXaP3jyvqDGm0GuwEY7DSS/LKuLEuWLOPGG1PAPxa49ToWLnyO++//3q6w\nXlPzU7ZuXUNNzbvZuvV0e4TKSH8hfW97WLvYI1gaYz1MDzfYOcdOKoK6urq9/vLQ6BlobuTJJ8/j\nq1/9CplMhqef/g233fZD4C62bj0LsP5cORnpBUzWHywN64MOn+VOJI05gykv0RXWjz76KKqrTwXO\n2uNcS9qUg5GuAWdJo9FnfdC9Y4+dpMQrNKQz2PISlrQpbz1XG58AZIBWYCshPLrXq42T/P6X61Dn\nSJUcGivssZOUaH3tIgHQ1tbKzTfPY+nSN7j55nm0tbX2GuZxV4Dy1lUDbtKkBcChwLeBFPBrOjsj\nP//5z/fq+kl9/8t53+Ikh+nR4OIJSRVlKL0MxVgxaUmb8pfNZjn00Bls2VL89yiJ73+5P6fRWBBT\nCYa7eMIeO0kVY6i9DMWYH9XVI5ROLySVmg9cRyo1n3R6obsClIlMJkOMIzMPLonvf7nPG3SLtb3j\nHDtJFWE4qxOLNaQz0rsClOtcp0ox0kN3SdsVotyHOt1ibe8Y7CRVhOFMqC7mlm8jVdLGsg57bzS2\n9ktSSaNK2AoxaWF6NDnHTlJFGKio8NKlb3D99Z/vcbTc5xKVe/sqRbFex7HSc+rnrjJYoFhSog2n\nl6HUQzoDBQXLOhTHYN7ngd6L5uZmzjhjPjt2HM+2bScwceJ9LF68jAcfTF7Paan/v9DIssdOUkXY\nm16G/rZ8G6lemsFsQzWcXkj1ra/3ec/3orb2Z3R2rmX+/DM5+eSTmDt3LrNm/TkdHXex52ertvYC\nNm16IZFhx60Qy5t7xQ7AYCdVvp5f0D17GYbTqzJSe4AONoRa1mHk9fdewF9RW/t+3nzzMbZtmwU8\nUeAKp3PrrWfxqU99alTaK3Wx3ImkxOuaUD1QUeHBGMltiwZbTsKyDiOvv/cCPkBHx+ls2/Yx4FfA\nrUB2j/Pex3e/++AotFQqDufYSaooxVqdOJLz2wZbTqJrrtMZZ5zD9u2Hs337PkyYsJkJE16kqenb\nDouxe6j817/+Da+99nv22+9A3vGOowc9ZN7fewEHAMuAU8gNh68GPgesArr+sfA44GiPKofBTtKY\nNJK1vAov9MgCGaqr76St7Xja2tpoamri4YcfZceObcBkoJ4QHgdeHvZjJ0U2m+W6667jppu+Sozv\nY8eOE4GngfuoqXlnnyVh9pwzuX37G8A6es9jzAJ3kQtxew7Rzie33+yPgUeYN+8LfbZxLKyiVWVx\njp2kMam/+W21tedy/vn7MGXK1GF9Yfee19VMLiycCBxHTc1P2Lr1+0yc+G62bTud3NyuH7O7p2hs\nl53IrVA9hy1btlE4eC0EGkmnP0FbWysxRjKZDA8//Cjf+ta9jBs3h46OE6itfZyOju8DteRCXPfr\nXAH8mlwv3Z7OBNqAF5g0qZNXX32p1/swUvMzpS4unhiAwU5Sd/1Pqj+b2toP0tFxfK8VlIMNeV1f\n/G++eTwdHY8B3yzwOAvJ9QzV9fp7rC6e2P2+XAz8Dri3wFnzgXmkUt/l059+G1/5yv/Pm2++l46O\nd5MLyT8CLiM3b24DcCNwMfBe4Fjgp8APgKvpa0VyVdVt1Nbu5MEH7y3YK2gdOI00F09I0hAU2gO0\ntvZc4Gzg83R03A/8Ix0d97J16918/esPcPnl9/e7N213XQs9zj9/f8aPfy+FJ+/XA5mCf5fD1k6l\nsHvu4x+Ad/dx1myglfb22Xz5y7eSza6ko+NeciHtf5GbE/dz4D3ANHKh7k5yPXHPAT8EDiP3Wi8D\nGum+aGL8+Me55JLT2bjxdwV738p9r1WNbc6xk1SRijG/ac9ti9ra9uGeez5IR8fiPc7sWkE5D/i7\nPvem3VNdXR07d0Z27HhPH2fkAkqhv1OpJ5k27RQaGxvH1Byu3NzHdwD/Cnygj7PWAWcyceK3ifFt\n7A5YWXK9eb3r0e3uDb0MeAg4l1wPXgpoIhfwVgFvUFPzBLfe2vf7O9J7rTp3T3vDYCep4hRzf9Xu\nq2yXLFlGR8fb+jizK3RdVnDVbKEv42w2y91330Pvnp0uTwLzCvy9hhgfYenSx+jsPIGOjmqqq1fx\n6U8v5v77v8Vpp502pOdYWd4E7gPmkJt3uIbeIW0t8DFifJTt2y/vdluGXK9nf72jC4ALgW8XuO4Z\n1NXV0NR0b7+hvb9dUMaPf5y2tgPIZrPDCmPuHay9FmMcEz+5pyppT5s3b4633357vOqqq+Ptt98e\nN2/eXOom9Wvz5s0xnT4owkMRYrefh2I6fVDMZrPDvvbtt98eU6lz97hu18+5ERrzv18bly5dtut+\njz32WEynD8rf97qYSp0ba2v3j+PHT4pwaoTC7c0dz3b7+y2xtvbDsa7uoDhp0r4Rvpw/J3ddOCvC\npLh69epivJRlZ/PmzbGu7oAItfnn+1i3539t/r9viePGHRbT6YPiVVddtcf7dXX+foXev2sjnBTH\njfvzCKf1cc5pcfHixYNqZ1+fwa73MJ0+KD722GNDfv4j9dlW5cnnlqHnneHcqRJ/DHZSb4UCyXC+\nkEZTf+ErlTo3NjY2Dvva/X9h7w5h3R+n//v0FVByAW3cuBMjXBtTqXNiTc2+ceHCS2JjY2O85ZZb\nYm3tmX0Gwqqq2njKKafGW265peyD+FDsfm8vygfimH/NGyMsy//3lPie97wnZrPZAq/97fnXuFBo\nOzVWVU2OcFS/4W/8+NSua/f3D56u/3dqa8/p9p4elH+vN0dYHMePTw3pPRrJz7Yqz3CDnYsnpDFq\nJHdeGEkjOb+p+4KK3EKK68gtplhIbv5VbvVq910h+t/Z4ChyBW7ryQ3jzgPeAM5i0qQ5XHzx2/I7\naJzJ73//EnfeuYLLLruMl15qo6Ojmr6GFTs7T+YHP3iFv/u7q5k8+fBBLeaoBLvn160GfkZueLSO\n3Ly4zwOHAz/lhBPeQ11dXYEFMC/k79N7Jw9YR2fneuAfyL0nhTxJCG/juuuuY+rUGSxa1MSNN6ZY\ntKip16KZ3Ytj9qG6+k7gLHbPl5wBPM+OHVfzD/+wetALbkZ67p7GBufYSWPUSO68MJL6m9+USj3J\nzJnzet8wBN0XVDz88CPcc88PGTfuZDo61pJKfXlXrbKuOVj972wwD/h/2D1PrOv1XEN19c/6nKA/\na9Z0qqtXsXPnJX1c9wTgL4A5vPHGBZxxxrls3Pi7siyxMZSFALNmTWfChK+wffv7ydWZm08u3M4m\nN/9wDdXVh3DMMW/fdZ89F8DE+Cluu20hMeb2E5448Sds27YWeJBcSGwArqTw3L1mtm//GP/6r/8v\n27fvnhuZ+6yt6bVoJsbIpk2/Z+fOKeRW4m7Ot3l3GZRt22Dbtt737ev5j+RnW2PEcLr5KvEHh2Kl\nHq66qv/5SN3nkJWTvZmHNJz5hNlsNjY2NsalS5fFxsbGXtcfeF7ekj2GYc+INTX79jvcvXnz5jhx\nYl1+eK+v6zbu+n3ChHeV5TDdUIf6N2/eHMePT3X7XO45DPuZQQ2Vdn/PPvShuRE+u8frd1WEuthz\n7l5uGHXChHfFiRPPHHA4tOu55c69Ln+Nvt+zwQylOsdO3eEcO4OdNBTlOp9nMOGrZ2C4dlBzA0dq\nPmH/c+z2z4eTroByYZw4sS5u3LhxwOuuXr06wqQ+rtt90cW1Ed4WFy++cq+eR7H1fF02x9z8t6sj\nLI51dQf0CCnd3/PTTz899re4YcmSJXH16tWxpmbfOH78vAjXxUmT5sWJE+viRRdd0uszU/hzvjnC\ngREWdwuN2QgP7REsC/+Dp+/3/IIB7zuQ4Xy2lUwGO4OdNCTl2DswlPA1UE9adyP9XAt9GdfW7h8n\nTdq31xf06tWrB9Vr+Nhjj8UJE1IRaiJ8KHb19uXC3m3dnsO8CO8YsBdwtO0OVF0LR7p6ts6MUBeX\nLFkSY+z9nucWjRQOtHV1B8Vvf/vbe9zedf2z8p+Zc3q8zosWXZnv/fz2Htf7coRJMZU6p8f703ul\nbe9/8PT9j6Lb8+/R3v1jaSifbSXXcIOdW4pJY1jPmlmzSaWeLNl+lyO5TVN/+8IOdeuuvuaMbdmy\nZdc8r5kzp/dYXNF17LDDDuP88y8ecH/RbDbLoYe+lS1bIvDvwO/JTcyfDhwEfILdm9Qv3PV7OW1n\ntWTJMm68sRq4jUIFg6urz+PFF5/lqKP+osB7/q/AP5JKnU57+7FMmvQjdux4mBkzZtLS8iydnR8g\nt8AiS26hQn/3fzc1NU+wdesaamrezdatp+36nN9zz5289NJLPd6zGGOfn8O6uoVs3NjKtdf+Mzfe\nmKL3dmRZcjta9N4+rpzeG1WG4W4pVvKetNH6wR47qaBy6R0YyaHhYs0n3Jvh3KH0Gt5+++1xwoR3\nxb5Ld8yL8K64u7xGcV6nYso9h3fG/oZVL7rooj7f89raD8ePfvSj8dRTT8/30J2Rfw8/1K3XslB5\nk82xrzIxEyfuGz/4wdPihz40t98yJI899lisrd0/3/Zr8693XZw0Kdcr2t9ntabmfbGmpndPbTn1\npqoyMMweO1fFSmNc950XSmkkSz0UY7Vh9/IwA62WLGQoq5BbWlrZvn0fchvWF3ICuZIdj5Jb6Um+\nPeVTEqOhoYG//utPk9vpoZD38eMf30F7+6UFb+3oOI66uldZs+ZucjtRdH/d7gP+itzK4Gnkesq6\nVtr2vfvEtm3H8eijv2f79vN47LEfsmzZ5wr2Tv/FX/wFVVVVwJ+RK09zLnAXb7zxY+bOnc+zzz7F\nFVcso9DK2vHjf8Nzzz3Dd7/73XxP4DwaGlbYU6dRY7CTVBZGstRDQ0NDn1/EuZp0K3rdZ88h161b\nt/Lmm+9luOVhhhJcc2U/vsf27ev6uFozuWDTMyyUU0mMdDrNu9/9Lp544sd9nPE427Ztpbb2cTo6\nuh/PAhmqq+/ku9/dBpxEz9e8GfifwCnkAtyPyA3HrmJ3vcC+AvH72L79DXI1G6GvUJ7JZIjxfcCX\n97h/7r1uamqiqWlVn9MYpkyZ0uuz4P6vGi0WKJZUFhoaGqiqaqZQcdnuBYGHo3ch2+tIpeaTTi/s\nUZOuS3Nzc68CtYsXL6Gj490Frz+YnrJccC0c1HKBbPquvxsaGpgw4UXgEQoX232U3Fy7nserqpqZ\nO3cujY2NLFmyjMbGRrLZbL/tGkmXXLKQ3HO4D7gVmEuutt+ngbW8+uo+dHSsJhfcbgUWAYcAt7Jz\n5wKef76KXFjrkmV3nbgHyc1xa8r/PR/YQm4e4hN9tOjJ/O1ddofy7gYTwuvr63n22ac477w0xx//\nIOedl+bZZ58qODe10OdpsEWLhyubzZbN50CjbDjjt5X4g3PspLI30qUeBjOfsO+5cIv7nC+WSp3T\nY4uxQqte+5tjV1Ozb2xra+tx36uuuirW1LwlwsQIsyKcHeEDsbZ2/3jbbbcVfJ2+9KUvxYkT62J1\n9VERzom1tWeWdH7X5s2b44QJdREmxFyNt3n5eXKn5o/tH7tWs+aOpSJ8Iu6uK3dV7FkXrr8tw86I\ncE6EuTG3knjP1/nb+TZcmb/O5j7nWA5mvudg51uWYvV5JW4VqN5wVWz/XBUrVYZCq0sHOz+pGMNd\nfa+gzQKHkhvy6zmcO3HiAr74xWv58Y9/wqpV36Gqag4dHe+mtvZndHauZf78Mzn55JM4/PDDOf/8\ni9mx40S2bv1L4D+B/6Sm5p1UVf2aGDupqjqZ9vZjqalZzdatP2P8+JPZseM9jBv3I6qrf8z992c4\n7bTTer1OtbW1XHTRZfm2HQesIzdseTXp9PWjsiJzz9f/wAMP5NxzLwRq6bkyNgscSaHVst1X+eZ+\nj8DX8+ctAwqtRoXc9m/XAovJ9RL+BjiZ3K4V9wO/zP/9XnLblf0IWEUq9eVeq6Kz2SxTphxJR0fv\n9tXWXsD69b/qYzVv79WvxVyRPRgjubpco8tVsfbYSWPa3q5YveWWW+Lpp58RjzxyZr7np6tHp3tx\n3foI+8SeOxakY3V1bb7+2lsK9sxAKo4f/85YV3dA/Pa3v52vq3Zh3F0Yt6tHqTbfo/QvEfYteK1U\nai7OqV8AACAASURBVP944YUXxuOOe0/8yEc+Ejds2BDvvffePnqpcsWMa2s/POKrZbte/9zrcE6+\n13BchD/Pv57d29Vfz1vPXTVgUcytgp0R4X/ke/b6ut/ifE/fZyL8Q/46F8S+Cz3vH+vqDuzVa7Z5\n8+Y4adK+Mdej2P29ztUmvOWWWwa9gnu0d3gpp8Ljw9npRbvhqlhJY9XerFhtbm7mtNM+zBtv7CQ3\nUf9Scj1dM4Dl+Z96chPy00AA3kZuteQpwEPs3Pltdu58ARhPocUVMIcdO9rYsSOyYMHFVFefTK4X\nCnK9V0vI1XubTq5m3a3kepf2vFYL7e0dfOMbfwLm8sQTP2LFipnAVuCDfTz28XR0PMu3vrWKBQsW\njMiE/d2v/9XA9cD/be/e46ws672Pf34MzAzMKCqi4jlSUVLzFKAgkKaooCiax5JO7p60187T46HH\nyr21du7cWqKWgbX1ycp2eMBA0VJSbEumpuYp3WImoI+HFATFAX7PH7/rdu5Zs9aaAzCHm+/79Vov\nZt3H677vWbN+XNfvuq4xaZ7bG4EaoDQ3sVoHh33SeoBBwE+I5zKGqOG8jxin7qzcPjHPa5zvJWA2\ncAbRI/dM4j6Uuzf7MnnyFjQ2NraobXz11UWYjQf+L9HL9kUiN/BG+vQ5ldmz72T58gPKlr5cR5iu\nnP91ffYu74iWY2TuS0PDHM4+++vdMkbmBqcz0WBvfKEaO5HC6mwtxdKlS72xcVCqmSmXk9XWlF75\nvLvqNTMxddXdHnle2Vyk2awJh6Z9j07nPLjMsRZVKU+9x5y0lc49zmtr11++3fTp072u7ghvPX7c\nNI8cwdL5U9tTY7c03aty19vf4ch0bYd4y9k4LvGo+ZzmkW83pOpzMav1M8880xsbB31Y2xu1jZX2\nuch3332PNJ1Z279vXZ1j1xNq7HrirDa9EZ2ssev2gKurXgrsRIqrs81d1QcCnu4xGG65Y050GJ6C\nlktTEHJqWl4tWMn2HezVBtKNZtgjS5af6pWmq4pyblulrFevty/WpUuX+vjxBzns6dHJI98xYZrD\nQWWuc6mXD6ZLg+ZKza4THbb36GzxEY9AbpO0334eAeEx6dns4ZUHSc46XByaypMN9jzdWwajWXP8\nZxwGeL9+h3mlZvdy97cr53/tCUFVTwgui6CzgZ2aYkWk1+tsc1f1gYBfpOVQG3mjgDeAy4FfAT8A\nRhLNga3HymtuJsz2nUsMz1F+XDwYD/y25FiPAZWGfBmdti137nnAAmIw30NYvnxfrrvuOs4555wP\nt2pvp5Nly5Zxww03cNttd/Dqq4upr+/PE0/8haYmS2UYRww/8nXgZ8SQJAcB/5KudwzR1PoosBI4\nlmhC3gd4gGhqPZF+/U6kqeke4JsVrndU+veXwPeBp4E7gLOBvxCdJbL78M9ER41y9+ZPxHNuJIZk\nORz4ctr+wbRN/1T2kamMt9HUdAjxTKcA+wOfSNPDPVh2+JwxY8awePGLuc4u62/Q4mxon0pj7HVF\nx4me0hy8wepMNNgbX6jGTqSwOltL0XaNXaUassM9auBKm0fzE9JfkmqDWk77FbVrRzsM9epNt/un\nGqm9Hf63Q98q5Tksre/vUZN4kUdzb6PDZ1Pt16BUo3WRg3m/fnVu1sche+3kMMn79h3mdXWNPnfu\n3Bb397zzzvO+fTdKNWID0zm/4dWbS/f05lrMZemefT39e4TDf3jURB7g8EmHvj527PjUueSMCs/F\nvXUHixnpPDXeutk3ey5ZE/gl6T7mn0v23LJpy6KjRNQEZtdXrvk4rqlv32E+derUHtXE2J1TBarG\nbt1ATbEK7EQ2ZB1p7sp665155jleW9vgsGmZ4KRajt3A9KU+3Vs3v2YBzFCPpsRlJfsOzu1bLghx\nj6bBvVOQMcmjyffAKuWpd9jFo8fuiLRdYwqesty9RofjPYLFLKfv0BS8bOQR3O3ozU2T/f3WW2/1\nuXPner9+AzzGnTsuHec2bw5+J3lzr+F8M+xRHk3AlZqby41tt5E3NGTlq9ZUnd3D/M+XOGzh1XLj\nIvgdmK4xey7VzjPAa2qy5uC2m/vVCzT0hObgIuhsYKdx7ESkMNozBl5pb736+j/y/vt3ET1axxPN\nfA8QzXSXEL1iR6TljxK9MocRc7W2Na7apcBOwOZE0+2rwG1Ek+Qyoudt6/HG4ETgbzRPGXYSsDuw\nKXBuKucBRHPhPKI5sx8xfttKYA2t51e9BziK6OV5VMnyTwN7E022g4G3iNkh/j0dqy+wA7AHMbvD\n48R4ftcRTaAHpm2eARYBM4nZH24AJhHNsqOI5syszDXArWXKeGRatgewGrgeGEs01z5E3P+TgTeJ\n5tBbiGbmkenco4nm4FJTiCbWa4kex3em5TPS9q3HmYNDied7a9XtGhqm8NWv7sK11/4k1wv0kQ+b\nPzvaC7QI04+1/Jy1bA5Wr9j26ew4dgrsRGSDUW3w1gh2NgY+QeRaTSUCq3eJYGE7YjL4I4DdgP8i\ngq/bieCm1BhikOBxRODzADEV2OXA6Wmb+UTgM4YIfBYQQc+dtMzvm0EEFHPS+Y5Lyz9J5O5tld5n\n6w6rUKbDiCCudEDcscQAvmOJIHYOEcQOA56j9cC+FwLfTvdmU+AfadvJRNB3DxEg7gE0pevejhjg\n+XXgBSL3LguuMvNTGceke5YNsvwFIhC8I5Xz08AEIufwz+m8+6X78Cta5thBy4GPjyMbnLg9gx73\n7fszVq16lmqBeGNjDKT87rs/b7Wuo4MCtx4mpPMBYndbm8HGRQMUqylWRNpULfcnmv8a2mj6i2X1\n9Qd6ff0mXlt7WGr+LNeMWzpg8FKPXLdah8tzTZZHOkz1yD07wiNvrbRsS1MT4m2pLGd55bzAj1Rt\nMoyBeysde7rD2R5Ns8eksp6VK2t2PzZJ15dv1s2mAXvgw2bMOM55Dr/wltN5HVumjG01vZ7h0SRd\n69FsfKhH82r2zLL9r0j/ZgMLH+XRNH1rOvfoVLYBHs3dE71S7uKAAUenfL/yOZS1tZO8tnZjnzBh\nQhqYufUxOpJTpiZMyaOTTbF91nWEKSLSEy1btoxf//oWli9/i6gBK50U/UBqatYQTXaTiGbUQ4na\noVtobha9h379nmXhwme44opJ1NR8kPaZkPaZRAyMO47m2ptswOOXiJ6ev0vvrySaF68masCOAR4u\nU/qNiKbY44gevPXp/Qyixim7nqyH5+8r3IX7gcUly24GdgY+Q9QmziKaX99NZX0plTWbsH4UUQs3\ngGiinEvUdt1C1GZlPUXHpW2vBD5P1MJtStQG3kM0SZeWYwzlewnvRTTBfoOo0atN13I3UbN3CNHs\nO5IYuDgbUPg9oiZ271SGG4les2PTfds6HSfrAZt3DzU1DzJr1kw22ugUGhqmAPMYMGAU/frdQ9++\n/wYs4oMPvsa99/ZlxYoHcvdoMVHjO5rlyxfy2GOP0x4333wza9aUvwdr1ozh5ptvbtdxZMOm4U5E\npPCy5q333x9BfGlmQ3LcQtbk2dDwKJdd9j0AZs+eAzzErrvuzvTpj+F+BcuXz2uRJ9TQ0EBdXR0n\nnHAcM2fewZo1r9LU9DP69HmbNWvWEAEHRMA1hfK5dJOJACcLGgcTTbHlhub4K9GcOQZYClxFzICx\nMXAXcAHRxLoYWFjhGP9N5PzNIIKfocB/Ak+lbT9OBGLlct+mpH1uTvvtRPkgbEza5oB07AFEs3W5\na59Fc75ftdkoxhJBWtZceg9wNHGPNyUCqguIJmKI+5k1Ny8DziECv9MqlOWo9BoPjGbAgEdYs2Ye\nRx45iZdffpnnnnucOXPm8MILL7LtthO44IKHePfd5ufZ1JS/RxcAF6VjHQHM55prZjB8+K6cfvrp\nVKNhQmRdKESOnZkdRgxm1Ae43t0vK7ONF+FaRSRkCeZPPfUsb731Optuujm7775bq0Tz6nl1zRPO\nNzaewpIlrXOhsjyhp556hjfffJ1BgwbTr18N11wzgw8+2Jampvcxewv3N9lss80ZNmwnttxyS267\nbSURQFZLzp9ATc0SVq8+ngEDHmH16nk0NvbjzTc/IHLQsjHf7icClduIP3NLiJyzcUQw9AhRS9cE\nrCBqp24hcscgat9WAh8BniWm6tokrcvG69s+HbeO1vl5y9K5hgAOvE3UiFXqNPIeMe7e20TeW7lr\nn0zUuB1M1OzdQgS2c8tsOymV+Wqi9hIiD3ILohPFH4kax78CvynZ9+q07DgqP4fDiKniXqFPnxeo\nqTH69fsUK1aMbJXjNmPGDM48cw7Ll5d/npFLWa7jymSWLHmRrbbaqsx+odqxGxqm8IMfTOSLXyzN\nj2y/InTK2JBssDl2xF+5F4guW/2ITNpdy2y3Fi3dItKTZEOb1NcflXK1DncY6HV1B3hdXaOffPJn\nPxxuonpe3cSUZ9Xo559/fqvhKhYtWuTTpk3zESP295qa+pRTlw0TUpvyu7IhRbLZCEamdVmOXVtT\nje2cctbMY9iVIR55crunfK7tveXQJYekY19RJhetwWEbh2Epr2ykN485l+3b4PAlj3Hy6tOy7B5m\n6/LHzfLKJua26+uVZ+U4xiMvb6DDqDLXns3iMM5jzL09PfILqw2NMtAjFzE/9lyWL9jokfO2NN3H\n/P4PeOTSXdqO5/D13Pk289JharIct+qznOxV5b4c7lOnTq36e70+c+xaDgd06Xqd/ULWDTqZY9ft\ngdnavoj/6t2Ze38BcH6Z7dbuDotIj1Dtyy+CmSPSF9fRvtFGg/3kkz/bxhf6RIeL/JRTPtvii6++\nPhs3LusccFQusGhrSqzbPYKmeoePeeWpxg7xGFS33iNg+pJHcJUFWxO8dSeM/HmWlSyf6BEc1ntz\nR4ty+w6qsq7R2x7jrdoYfw0Om3uMu3eqt+zkkQWJ2XRfk9JxrihZn80Dm7/npdd9qEcAmX++c9Px\njvLoWNHozR1N2jM/baX3zZ0gqv1HwWxzr/a7NmrU6DZ/v9fH9GPqlNE7bciB3bHAj3PvPwNcVWa7\ntbvDItIjVK+BO6rkC/lur6tr9Nra8j0Wsy/w2tpJJb0f2+qhOa0dQcIxHjMrfNIrB2ebevT0vLvC\nOasNYtw6+IhgaBuPILFaILO3R1BVbl0WMHkbxzjQo4Ys6316WLrOUzwC28EOi3PXVO2eZgF5dpy+\nHgMuz/DWwWtWI9jg0SN5d28ZOM/1COi2TPchO2+1QLc0SM7X4DUvywYhrhQk9e1b59VmK2mrxi6z\nrmeN0EwQvVNnAzt1nhCRXqVagnmMQZdPMD+EmprxrFw5j8pzuH6O1avnUVf3ydz6aj00xxDjrx1Q\noQz7pDLsQ+R/3Uv0Kp2S9s/y5uYTY8ZZWj6D6NWZP+eLaZtq58l7gMiRG031zggb0zzfaqkDge8S\n+Wrv0voeZA5N5d01naue6H36OjG+3HeInLxbiGvfOl1LuXs6jugpPIjIrhlF5L2VyyfbJ5Xvq8QA\nxo+kfbPneyiRK3gwcR82ypVhV6KX8wE0D/D8CC17PUM8n5bzC2dzDlebi/XGG3/JMcecTPnftXl8\n97vt6/zQ2Ni4Vrl0pdQpY8NShMBuEZH1m9k2LWvl4osv/vDn8ePHM378+PVZLhFZD3beeSgNDXNY\nvrzc2tZfyCtWjGTkyDdYsGAK0SHhE8Rk838A9gROYb/9dmfBgpG5vaoFRfsQvUsfqbA+K8PsXFm2\nIzpq7EDzUBw3Esn2R+TOWRpsDaX8LAoQw6KUziDxp3TMB4lAptK+S4lhViqV/zJihoaXiV6tlbab\nSHPwdSTRQ/XPRG/dfyUCrmzGiWyGi3JGE4HkLkSQtZLK9/d+4j42Ule3APdH+PSnj+XXv/40K1fu\nn471KPA0zQM3j6G5R++z1NT8inHjmnjwwYdYufImWg4GfQ/REeXGFsv69JnPCSfEsjFjxrB48Yu5\nwXcncsIJN9LY2Mg111zOGWdMpnl2kD8A87jmmsurdpxYn6p9ZrKAVbrfvHnzmDdv3tofqDPVfD3p\nRXQNyzpP1BJ/VXYrs93a1oqKSA9QPceudd5ZQ8MxfvXVV3tj4+apCe9cjwGBz3U4yxsbN/dp06aV\nNFW1lY91tbedY5cvy/DU1Fh6rFO9Odk+m3c1v76t5svDUrNh1iniAYdFaV1nc+yy8jemV7XrXFZm\nv0Fpv9s9mlK/nv49I5Wz3D2d5NH8elK6t4OqnHcTr61t8LPPPrdFM+XixYvTvL8npfPlm4FbHiPL\nKyuXzzZgwGbev/8ma5XjtmTJEp86daqPGjXap06d6kuWLFlfH4d2UY5d70Qnm2K7PTBbFy+izv45\n4HngggrbrN0dFpEeo2Wv2Es88prySfjt+xLPvrBbf/FVC6gGOhzkLXvFZrlhAz1mNihN+B/g5Tsb\n/MKb8+8qnfMKb+7EcUn6t95hhMdk9jumoCjfI/Qab90rNh8AlvaKzZd/v9w2WQ/fBo88tuwY9Wm7\n7N4P9Mi5y2Z+6F9y7AnePNtDpXs6IrfPl7x5Zov8eQd4//6bVAyyHnjgAW9sHOx1dZMcLvF+/Q5w\n6P/h70m5IK1cPtu6znHrCdZHpwxZvzob2BViHLv20Dh2IsWSjS/39NPP8uabr/PeeyuZNesuamrG\npzk2W086Xm3uytJJy+vr7+b99/9E5JwdSDYBfU3NGs4992uceeaZzJw5k9mz57By5QcsX76Md95Z\nxgsvLKS29uA0BtqjmM3njDO+wEsv/Y2ZM2djdiBNTfszYMAjrFp1L01NK3DvS8zHug3wC6IJb2w6\n5++BVcRsEJsDbwD9ify3QcRAwVmzZV+iKfJAYqDjeYATI0HVEePROZFP1kTkqK0i8vzi72NjYyMH\nHXQQQ4Zsw5NPPs777zexatVKFi58iXffXUGfPrD99juw994f53/+ZyGPP/4k0dT88XT8B4mGlJOJ\nZtV/EM25Q+jbdw2rVi2mX79P0tS0P336PEifPvPZa6/hDB68OZ/61MHU1dXxyiuLGTx4EA899BAP\nP/wITU1NDBu2C1OmHMPUqVOrzjda+ownTpzI7NmzNV8pmru1t+nsOHYK7ESkMNb2i6tcUBDB252A\nM3HiER0OLPJlKLcOYNq0aVx99TUsX/4eQ4fuwJAhQ1iw4E+8995yVq9eDRibbLIJo0fvz447DmX4\n8F1bBCzbbrs1K1eu5K675rJo0SLWrHEiUHPM+rDtttsxefKRHHvssR8Go2BMnHh4m9eztvcsH6wp\n0BJpPwV2bVBgJyIiIr1FZwO7PuujMCIiIiLS9RTYiYiIiBSEAjsRERGRglBgJyIiIlIQCuxERERE\nCkKBnYiIiEhBKLATERERKQgFdiIiIiIFocBOREREpCAU2ImIiIgUhAI7ERERkYJQYCciIiJSEArs\nRERERApCgZ2IiIhIQSiwExERESkIBXYiIiIiBaHATkRERKQgFNiJiIiIFIQCOxEREZGCUGAnIiIi\nUhAK7EREREQKQoGdiIiISEEosBMREREpCAV2IiIiIgWhwE5ERESkIBTYiYiIiBSEAjsRERGRglBg\nJyIiIlIQCuxERERECkKBnYiIiEhBKLATERERKQgFdiIiIiIFocBOREREpCAU2ImIiIgUhAI7ERER\nkYJQYCciIiJSEArsRERERApCgZ2IiIhIQSiwExERESkIBXYiIiIiBaHATkRERKQgFNiJiIiIFIQC\nOxEREZGCUGAnIiIiUhAK7EREREQKQoGdiIiISEEosBMREREpCAV2IiIiIgWhwE5ERESkIBTYiYiI\niBSEAjsRERGRglBgJyIiIlIQCuxERERECkKBnYiIiEhBKLATERERKQgFdiIiIiIFocBOREREpCAU\n2ImIiIgUhAI7ERERkYLoMYGdmR1nZn8xs9Vmtk/JugvN7Hkze8bMDs0t38fMnjCzv5rZ97u+1CIi\nIiI9R48J7IAngWOA3+cXmtluwPHAbsDhwLVmZmn1D4EvuvsuwC5mNqELyytdZN68ed1dBFkLen69\nm55f76Vnt2HqMYGduz/n7s8DVrJqMvBLd1/l7i8BzwMjzGwrYCN3fzhtdyNwdJcVWLqM/jj1bnp+\nvZueX++lZ7dh6jGBXRXbAH/PvV+Ulm0DvJJb/kpaJiIiIrJB6tuVJzOze4At84sAB/6Pu9/RlWUR\nERERKRpz9+4uQwtmdh9wjrs/mt5fALi7X5be3wV8C/gbcJ+775aWnwiMc/evVDhuz7pQERERkSrc\nvTQ9rU1dWmPXAfkLmQXcZGZXEk2tOwF/dHc3s3fMbATwMHAqcFWlA3bm5oiIiIj0Jj0mx87Mjjaz\nvwOjgN+Y2Z0A7v408CvgaWAOcLo3VzOeAVwP/BV43t3v6vqSi4iIiPQMPa4pVkREREQ6p8fU2K1r\n1QY8LtnuJTN73MweM7M/dmUZpbwOPLvDzOzZNED1+V1ZRqnMzDY1s7vN7Dkzm2tmAytsp89eD9Ge\nz5KZXZUGiv+zme3V1WWUytp6fmY2zszeNrNH0+ui7iintGZm15vZa2b2RJVtOvTZK2xgR4UBj8tY\nA4x3973dfcT6L5a0Q5vPzsz6AFcDE4CPASeZ2a5dUzxpwwXAb919GHAvcGGF7fTZ6wHa81kys8OB\nj7r7zsCXgR91eUGlrA78Lbzf3fdJr0u7tJBSzU+JZ1dWZz57hQ3sqgx4XMoo8H3ojdr57EYQeZV/\nc/cm4JfEYNbS/SYDN6Sfb6DywOH67PUM7fksTSYGgcfdFwADzWxLpCdo799CdSDsgdx9PvCPKpt0\n+LOnP6oxjt49ZvawmZ3W3YWRdisduFoDVPccW7j7awDu/iqwRYXt9NnrGdrzWao0ULx0v/b+Ldw/\nNeXNNrPhXVM0WQc6/NnrqcOdtMs6GvB4tLsvMbPBxJfMMymClvVIg1X3blWeX7ncnUo9tPTZE+ka\njwDbu/uK1LR3G7BLN5dJ1pNeHdi5+yHr4BhL0r+vm9mtRLW2vlzWs3Xw7BYB2+feb5uWSReo9vxS\nIvCW7v5amtP5/1U4hj57PUN7PkuLgO3a2Ea6R5vPz93fzf18p5lda2abuftbXVRG6bwOf/Y2lKbY\nsrkFZjbAzBrTzw3AocBfurJg0qZKeSEPAzuZ2Q5mVgucSAxmLd1vFvC59PNU4PbSDfTZ61Ha81ma\nRQwCj5mNAt7Omtul27X5/PI5WWlQf1NQ16MYlb/rOvzZ69U1dtWY2dHANGBzYsDjP7v74WY2BJju\n7pOIpqRb03RjfYGb3P3u7iu1QPuenbuvNrOvAncT/0G53t2f6cZiS7PLgF+Z2ReIqf+OB9Bnr2eq\n9Fkysy/Hav+xu88xsyPM7AVgOfD57iyzNGvP8wOOM7OvAE3Ae8AJ3VdiyTOznwPjgUFm9jIxZWot\na/HZ0wDFIiIiIgWxoTTFioiIiBSeAjsRERGRglBgJyIiIlIQCuxERERECkKBnYiIiEhBKLATERER\nKQgFdiLSLcxsnJmtMbPNcu9XZ+/X0zmPM7M1bWyz0MzOXl9lyJ1nmJn9wczeM7MX27nPt8zsyTa2\nmWZm961FuX5qZhrsW6SXKuwAxSLSK+QH0nwQGLKeR8R3Ks9d29UuJQYc3QVY0YH92lP+nnKNItLF\nFNiJSFlmZsQg5lVruNYVd19FhXllC2on4DZ3/3t3F0REikNNsSIFYGb3mdk1ZvZtM3vdzF4zs++V\nbLOJmd1gZm+Z2Qozu8fMhufWTzWzZWZ2eGruWwnsmprm7jCz88xsiZm9bWbfsXBxOtcSMzuv5Hxn\nmdnjZvaumb1iZtPNbGCVayhtml2Y3q9JTbTZz9un9Rub2Y/T+Zeme7BvyTFPNbOXUhlmEVOZdfTe\nbmdmt6ZzLDWzmWa2Tck2F5rZq2b2jpldb2bfMLOFVY65BtgT+Fa6tm+m5Xuk57LCzN5M937jKsfp\nY2aXp2f6ppldCdS045qGmdnt6VkuM7MHzexjJdv8c3pub5nZT8ysPrdugpndnzvvXWa2a279DulZ\nTTGzu81suZk9ZWafKjnHRDN7NjVH32tmx+efcdrmADObl47xisUE9hu1dY0iGyoFdiLFcTIxF+T+\nwBnAmWaWnxPyBuATwJHp3xXAXWZWl9umHrgI+CdgOPByWj4W2BEYB3wZOB+YA/QDRgMXA981s71z\nx1oNfC0d56R0zqvauIZ8E+J+wFbpNQT4DfA0kE2APSetOwLYC7gf+J2lCc/NbCTwU+BHaf0dwL+2\ncf4WUq3lLGBwuvbxwNbArbltTgS+CVwI7As8D5xN9ebQrYC/Apena7vczAYAdwFL07UfDRwAXF/l\nOOcCXwROI557DXBKG9c0BJhPPJ+DgY8TzyUfEI4FPpbWHw8cQzzLTANwZSrnOOBt4A4zK20FuhT4\nPhHEPgz8Il0nZrYdMJN4LnsCVwP/Tu6+mdkewFzgNmCPVI6Pt3FPRDZs7q6XXnr18hdwH/BgybK7\ngR+nn3cG1gCjc+s3Jr6Qv5DeTyW+7PcqOc5Pgb+R5pZOyx4GHivZbiFwdpUyTgDey70fl863Wbn3\nJfueTzTT7pjeH0QEQHUl2z0GnJt+vgmYW7J+OrC6jXv54XUAhxDB8na59R9J5Twovf8DcE3JMeYC\nL7ZxnieBb+benwb8AxhQco/WAEPT+28BT+TWLwIuyL034Dng3irn/Xa6xpoK68s97x8Dd1c5ZgOw\nCjggvd8hlftLuW22Tsuybf4NeKrkOBeme7t9en8DML1km73ScTbv7s+dXnr1xJdq7ESK44mS94uB\nLdLPuxJfmA9lK919KRFcDM/tswp4vMyxn3b3fA3Ua8BfSrZ5LXc+zOyg1Az3dzNbCtwC1JrZVu2/\nJDCzI4mAZoq7v5QW70MEE2+kpsRlZraMqGUamrbZDfjvksOVvm/LrsBiz+XBuftC4t4Oz23zcMl+\nCzp4nuw4T7h7viPFH4ggZnjpxqmJdggtn6m349x7AfPdfXWVbUqfd/53CTMbamY/N7MXzOwdApJn\nVQAAA6lJREFU4FUiqNy+5Dgf9uB198Xpx+w4w2j7vu0LfKbkGc8navU+WqX8IhssdZ4QKY6mkvdO\n+9It8l/gK0u+0Ksdu+L5Uo7Ub4DrgG8AbxJf0j8HattRJtJxdgd+Bpzu7vNzq/oQwcQYIqDIW9re\n46+lrux52tW9XNv6XZpNNNP/E1FruAp4htbPtvQ40LEUoD7ADOAKWj/nRR04jsgGQzV2IhuGZ4jP\n+/7ZglTjswfw1Ho4335E/t3Z7r7A3V8AtmljnxbMbHMiv+06d//PktWPEh0h3N1fLHm9kbZ5BhhV\nst/+dMwzwNYlyfxDiWbF7L49S+QP5o3s4Hmyc+1hZg25ZaOJgOaZ0o1TjesSWl/jiDbO8xgwpkw+\nXLukzi3DgO+4+73u/hwwkI5XFDxL/J7kld63R4GPufvCMs95ZWfKL1J0CuxENgApsJoFXGdmY1JS\n+s+Ad4BfrIdTPk/8fTnLzHY0s5NomXyfKa2Fyb+fCbwCXGlmW+Ze5u6/Jca9u93MDkvn2N+il+7o\ntP9VwKfM7AIz28nMTiM6JLRbOs+TwE1mtq+Z7Ufctz+5+7y02Q+Az5nZ59N5ziOCq47Wst1EdGi5\n0cx2N7OxRMePme5eaQDjHwDnmdmxZraLmX2faJ6t5lqgEfgvM9vPzD5qZiea2Z7tLOc/gDeA09K+\n44AfUr52rpofAR81s++lsk8hagCh+d5dBowwsx+a2V7pfJPM7EcdPJfIBkOBnUgxtCeI+BzwR+B2\nIi+rDjhsHdZ8fFgGd3+SCOTOImq2vgCcU22fMu8PJGqsXiFyvJakf7dL648A7iUS+58FfkkM9rs4\nlWEB0WP0fxF5g0cTuXrtvo7kKOD1dK7fpeMfk7vWm4FLiM4AjxL5cD8C3u/Iedz9PaKDycZErtmt\nRPD6xSrH+A+is8N04pkaEXhWPmnkuo0lalTvTWX+KtGc2qbUVH880ZP1SWAa0ZO69Peo3O9k/nfk\nZeBYopf2n4nfl39Jq99P2zyZyroDMC9t922iGV5EyrDy6TQiItJZZnYL0et0cneXpTcxs68BF7v7\npt1dFpHeSp0nRETWgpn1B75CjEG3mqiFOgqY0p3l6g3M7HSiZ+zrRP7jRUQNpIh0kgI7EZG148Dh\nxBhs/Yn8wlPcfVa3lqp32An4OrAZ0eR+LdGsLSKdpKZYERERkYJQ5wkRERGRglBgJyIiIlIQCuxE\nRERECkKBnYiIiEhBKLATERERKQgFdiIiIiIF8f8BqIq04xW4JyYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103f5a1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "plt.scatter(fold_change.values, pvals.values, s=50)\n",
    "plt.scatter(fold_change.loc[ids].values, pvals.loc[ids].values, c='r', s=50, label='ANCOM significant')\n",
    "#plt.ylim([0, 20])\n",
    "plt.ylabel('W', fontsize=14)\n",
    "plt.xlabel('normalized log fold change', fontsize=14)\n",
    "plt.legend(scatterpoints=1)\n",
    "plt.savefig('../results/volcano_plot.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}