eating-responsive-bacteria.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-22T06:39:04.157456Z",
     "start_time": "2021-04-22T06:39:02.962351Z"
    }
   },
   "outputs": [],
   "source": [
    "import dill\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import calour as ca\n",
    "from bisect import bisect_right\n",
    "from itertools import combinations\n",
    "\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-22T06:39:04.162505Z",
     "start_time": "2021-04-22T06:39:04.159427Z"
    }
   },
   "outputs": [],
   "source": [
    "pd.set_option('max_colwidth', 200)\n",
    "pd.set_option('max_rows', 1000)\n",
    "pd.set_option('max_columns', 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# import relative abundance ASV table and events windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:31.134051Z",
     "start_time": "2021-04-02T01:18:31.101054Z"
    }
   },
   "outputs": [],
   "source": [
    "with open('../data/ASV-table-and-eat-windows-for-MA.pkl', 'rb') as f:\n",
    "    saliva_ma, ma_eat_windows = dill.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:32.768784Z",
     "start_time": "2021-04-02T01:18:32.762666Z"
    }
   },
   "source": [
    "**saliva_ma**: calour.AmpliconExperiment  \n",
    "relative abundance table of MA samples.  \n",
    "**ma_eat_windows**: tuple of tuples\n",
    "((eat_start_timepoint1, eat_end_timepoint1), (eat_start_timepoint2, eat_end_timepoint2)...)  \n",
    "each tuple is an eating window.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:32.971189Z",
     "start_time": "2021-04-02T01:18:32.965049Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AmpliconExperiment with 207 samples, 531 features"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# sample and feature counts of each subject\n",
    "saliva_ma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BarcodeSequence</th>\n",
       "      <th>Primer_Plate</th>\n",
       "      <th>Well_ID</th>\n",
       "      <th>Plating</th>\n",
       "      <th>Extraction_robot</th>\n",
       "      <th>TM1000_8_tool</th>\n",
       "      <th>Processing_robot</th>\n",
       "      <th>TM300_8_tool</th>\n",
       "      <th>TM50_8_tool</th>\n",
       "      <th>Orig_name</th>\n",
       "      <th>Description</th>\n",
       "      <th>drink</th>\n",
       "      <th>Host_ID</th>\n",
       "      <th>Barcode</th>\n",
       "      <th>eat</th>\n",
       "      <th>other</th>\n",
       "      <th>sleep</th>\n",
       "      <th>events</th>\n",
       "      <th>body site</th>\n",
       "      <th>Day</th>\n",
       "      <th>Hour</th>\n",
       "      <th>Min</th>\n",
       "      <th>dayhourmin</th>\n",
       "      <th>DATETIME</th>\n",
       "      <th>_sample_id</th>\n",
       "      <th>_calour_original_abundance</th>\n",
       "      <th>observed_otus</th>\n",
       "      <th>shannon</th>\n",
       "      <th>faith_pd</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>#SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>19122</th>\n",
       "      <td>ATCGTAGTGGTC</td>\n",
       "      <td>2</td>\n",
       "      <td>F6</td>\n",
       "      <td>GG</td>\n",
       "      <td>NEWE</td>\n",
       "      <td>311172B</td>\n",
       "      <td>ROBE</td>\n",
       "      <td>311313B</td>\n",
       "      <td>111075B</td>\n",
       "      <td>19122</td>\n",
       "      <td>Amnon_Med_2_19122_F6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>MA</td>\n",
       "      <td>19122</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>yes</td>\n",
       "      <td>sleep</td>\n",
       "      <td>saliva</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>1:07:00</td>\n",
       "      <td>2014-09-15 07:00:00</td>\n",
       "      <td>19122</td>\n",
       "      <td>13584.0</td>\n",
       "      <td>75</td>\n",
       "      <td>3.633408</td>\n",
       "      <td>6.924340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19123</th>\n",
       "      <td>GCTGCCCACCTA</td>\n",
       "      <td>2</td>\n",
       "      <td>C5</td>\n",
       "      <td>GG</td>\n",
       "      <td>NEWE</td>\n",
       "      <td>311172B</td>\n",
       "      <td>ROBE</td>\n",
       "      <td>311313B</td>\n",
       "      <td>111075B</td>\n",
       "      <td>19123</td>\n",
       "      <td>Amnon_Med_2_19123_C5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>MA</td>\n",
       "      <td>19123</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>saliva</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>30</td>\n",
       "      <td>1:07:30</td>\n",
       "      <td>2014-09-15 07:30:00</td>\n",
       "      <td>19123</td>\n",
       "      <td>5040.0</td>\n",
       "      <td>71</td>\n",
       "      <td>3.998352</td>\n",
       "      <td>7.172192</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          BarcodeSequence  Primer_Plate Well_ID Plating Extraction_robot  \\\n",
       "#SampleID                                                                  \n",
       "19122        ATCGTAGTGGTC             2      F6      GG             NEWE   \n",
       "19123        GCTGCCCACCTA             2      C5      GG             NEWE   \n",
       "\n",
       "          TM1000_8_tool Processing_robot TM300_8_tool TM50_8_tool  Orig_name  \\\n",
       "#SampleID                                                                      \n",
       "19122           311172B             ROBE      311313B     111075B      19122   \n",
       "19123           311172B             ROBE      311313B     111075B      19123   \n",
       "\n",
       "                    Description drink Host_ID  Barcode  eat other sleep  \\\n",
       "#SampleID                                                                 \n",
       "19122      Amnon_Med_2_19122_F6   NaN      MA    19122  NaN   NaN   yes   \n",
       "19123      Amnon_Med_2_19123_C5   NaN      MA    19123  NaN   NaN   NaN   \n",
       "\n",
       "          events body site  Day  Hour  Min dayhourmin            DATETIME  \\\n",
       "#SampleID                                                                   \n",
       "19122      sleep    saliva    1     7    0    1:07:00 2014-09-15 07:00:00   \n",
       "19123        NaN    saliva    1     7   30    1:07:30 2014-09-15 07:30:00   \n",
       "\n",
       "          _sample_id  _calour_original_abundance  observed_otus   shannon  \\\n",
       "#SampleID                                                                   \n",
       "19122          19122                     13584.0             75  3.633408   \n",
       "19123          19123                      5040.0             71  3.998352   \n",
       "\n",
       "           faith_pd  \n",
       "#SampleID            \n",
       "19122      6.924340  \n",
       "19123      7.172192  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# display sample metadata\n",
    "saliva_ma.sample_metadata.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>_feature_id</th>\n",
       "      <th>taxonomy</th>\n",
       "      <th>taxa</th>\n",
       "      <th>kingdom</th>\n",
       "      <th>phylum</th>\n",
       "      <th>class</th>\n",
       "      <th>order</th>\n",
       "      <th>family</th>\n",
       "      <th>genus</th>\n",
       "      <th>species</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG</th>\n",
       "      <td>AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG</td>\n",
       "      <td>k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Peptococcaceae;g__Peptococcus;s__sp._HMT_167</td>\n",
       "      <td>s__sp._HMT_167</td>\n",
       "      <td>k__Bacteria</td>\n",
       "      <td>p__Firmicutes</td>\n",
       "      <td>c__Clostridia</td>\n",
       "      <td>o__Clostridiales</td>\n",
       "      <td>f__Peptococcaceae</td>\n",
       "      <td>g__Peptococcus</td>\n",
       "      <td>s__sp._HMT_167</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG</th>\n",
       "      <td>AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG</td>\n",
       "      <td>k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Pasteurellales;f__Pasteurellaceae;g__Haemophilus;s__parainfluenzae</td>\n",
       "      <td>s__parainfluenzae</td>\n",
       "      <td>k__Bacteria</td>\n",
       "      <td>p__Proteobacteria</td>\n",
       "      <td>c__Gammaproteobacteria</td>\n",
       "      <td>o__Pasteurellales</td>\n",
       "      <td>f__Pasteurellaceae</td>\n",
       "      <td>g__Haemophilus</td>\n",
       "      <td>s__parainfluenzae</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                                                                                                                                                                                                                                                                   _feature_id  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG  AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG   \n",
       "\n",
       "                                                                                                                                                                                                                                                                          taxonomy  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG                    k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Peptococcaceae;g__Peptococcus;s__sp._HMT_167   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Pasteurellales;f__Pasteurellaceae;g__Haemophilus;s__parainfluenzae   \n",
       "\n",
       "                                                                                                                                                                     taxa  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG     s__sp._HMT_167   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  s__parainfluenzae   \n",
       "\n",
       "                                                                                                                                                            kingdom  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG  k__Bacteria   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  k__Bacteria   \n",
       "\n",
       "                                                                                                                                                                   phylum  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG      p__Firmicutes   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  p__Proteobacteria   \n",
       "\n",
       "                                                                                                                                                                         class  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG           c__Clostridia   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  c__Gammaproteobacteria   \n",
       "\n",
       "                                                                                                                                                                    order  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG   o__Clostridiales   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  o__Pasteurellales   \n",
       "\n",
       "                                                                                                                                                                    family  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG   f__Peptococcaceae   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  f__Pasteurellaceae   \n",
       "\n",
       "                                                                                                                                                                 genus  \\\n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG  g__Peptococcus   \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  g__Haemophilus   \n",
       "\n",
       "                                                                                                                                                                  species  \n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTAAATTAAGTCAGGTGTGAAAGTTCGGGGCTCAACCCCGTGATTGCACCTGATACTGATAAACTAGAGTGTTGGAGAGGTAAGTGGAATTCCTAG     s__sp._HMT_167  \n",
       "AGCGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTGACTTAAGTGAGGTGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGGTCGCTAGAGTACTTTAGGGAGGGGTAGAATTCCACG  s__parainfluenzae  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# display feature metadata\n",
    "saliva_ma.feature_metadata.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# filter ASVs\n",
    "Before eat-responding bacteria analysis, we removed ASVs with low prevalence and low unique counts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## filter ASVs with low prevalence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:40.590034Z",
     "start_time": "2021-04-02T01:18:40.585971Z"
    }
   },
   "outputs": [],
   "source": [
    "prev = 0.1\n",
    "cutoff = 0.0001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:40.816500Z",
     "start_time": "2021-04-02T01:18:40.804684Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(AmpliconExperiment with 207 samples, 531 features,\n",
       " AmpliconExperiment with 207 samples, 176 features)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# keep features present in more than 10% of MA samples\n",
    "saliva_ma_f = saliva_ma.filter_prevalence(prev, cutoff=cutoff) \n",
    "saliva_ma, saliva_ma_f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## filter ASVs with low unique counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:44.594068Z",
     "start_time": "2021-04-02T01:18:44.539142Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AmpliconExperiment with 207 samples, 133 features"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# count unique_counts for each ASV, and compare it with cutoff, return a boolen list \n",
    "b = [len(set(saliva_ma_f[:, i])) > n + 1 for i in saliva_ma_f.feature_metadata.index]\n",
    "# filter ASVs by the boolean list\n",
    "asvs = saliva_ma_f.feature_metadata.index[b]\n",
    "# keep filtered ASVs\n",
    "saliva_ma_f.filter_ids(asvs, axis='f', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ensure that samples are sorted in time order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "samples are already sorted in time order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-02T01:18:34.633652Z",
     "start_time": "2021-04-02T01:18:34.625810Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "saliva_ma.sample_metadata['DATETIME'].is_monotonic_increasing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# identify eat-responding bacteria based on a z-score method\n",
    "We developed a z-score method to measure how much the abundance change induced by eating deviates from normal fluctuations for a given ASV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we calulated the z-score value of one ASV in MA as a demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "asv = 'TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## calculate the distribution of normal fluctuations \n",
    "For each ASV in each subject, we computed a vector (D) consisting of the abundance differences between every 2 samples that were collected within 1 hr. The vector D represents the distribution of normal fluctuations in abundance for this ASV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "tps = saliva_ma.sample_metadata['DATETIME'].tolist() # sampling datetime\n",
    "abund = saliva_ma[:, asv] # relative abundance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-06T16:12:44.089862Z",
     "start_time": "2020-12-06T16:12:44.085935Z"
    }
   },
   "outputs": [],
   "source": [
    "window = '1H' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def monotonic(lst):\n",
    "    '''test if the list is monotonic'''\n",
    "    diff = np.diff(lst)\n",
    "    return np.all(diff <= 0) or np.all(diff >= 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_abund_diff(abund_lst, tps_lst):\n",
    "    '''calculate the monotonic abundance difference.\n",
    "    abund_list : list\n",
    "    tps_list : list\n",
    "        abund_list and tps_list should be paired\n",
    "    \n",
    "    returns\n",
    "    -------\n",
    "    max_value, t_start, t_end\n",
    "    max_value is the monotonic abundance difference within this time window.\n",
    "    t_start and t_end are timestamps.\n",
    "    max_value = the relative abundance at t_end - the relative abundance at t_start\n",
    "    '''\n",
    "    from itertools import combinations\n",
    "    v_abs_max = 0\n",
    "    v_a = 0\n",
    "    v_b = 0\n",
    "    # calculate the diff between every two data\n",
    "    for a, b in combinations(range(len(abund_lst)), 2): \n",
    "        if monotonic(abund_lst[a:b]): # make sure the list is monotonic\n",
    "            v = abund_lst[b] - abund_lst[a] # b is always greater than a\n",
    "            v_abs_max = max([v_abs_max, v], key=abs) # update the max value\n",
    "            if v_abs_max == v: # update the index to get tps\n",
    "                v_a = a\n",
    "                v_b = b\n",
    "    return v_abs_max, tps_lst[v_a], tps_lst[v_b]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the distribution of normal fluctuations\n",
    "v_ta_tb = [] \n",
    "for i in range(len(abund)):\n",
    "    # iterate from the first time point, slice time series to 1 hr window\n",
    "    tps_i = tps[i]\n",
    "    tps_end = tps[i] + pd.to_timedelta(window)\n",
    "    tps_end_index = bisect_right(tps, tps_end)\n",
    "    tps_within_window = tps[i:tps_end_index]\n",
    "    abund_within_window = abund[i:tps_end_index]\n",
    "    # calculate the monotomic abundance difference within the window\n",
    "    v_abs_max = 0\n",
    "    v_a = 0\n",
    "    v_b = 0\n",
    "    if len(abund_within_window) > 1:\n",
    "        v_abs_max, t_a, t_b = calc_abund_diff(abund_within_window, \n",
    "                                              tps_within_window)\n",
    "        v_ta_tb.append([v_abs_max, t_a, t_b]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.03916960438699569,\n",
       " Timestamp('2014-09-15 07:00:00'),\n",
       " Timestamp('2014-09-15 07:30:00')]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_ta_tb[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get vector D (distribution of normal fluctuations)\n",
    "D = [v_ta_tb[i][0] for i in range(len(v_ta_tb))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAM9klEQVR4nO3db4hl913H8ffHTdqKrTQxs3HJHyeBWBoFmzrUQtQHptFt0nZTNJCCZcHIPjGQgqJb86TSJ1vF4hOxrDa4aGgJtDVLgrTr2liEmjqb5i+buGmNNWbZ3aZIKko09uuDe9ZOJjM7d2bn/vkm7xcM95zfuTP3k5Mznz17/m2qCklSPz8w6wCSpK2xwCWpKQtckpqywCWpKQtckpq6YJofdskll9Ti4uI0P1KS2jt27Ni3q2ph9fhUC3xxcZHl5eVpfqQktZfkX9Ya9xCKJDVlgUtSUxa4JDVlgUtSUxa4JDVlgUtSUxa4JDVlgUtSUxa4JDU11Tsxpe22uP+BNcefPXDzlJNI0+ceuCQ1ZYFLUlMWuCQ1ZYFLUlMWuCQ1ZYFLUlMWuCQ1ZYFLUlMWuCQ1ZYFLUlNjF3iSHUm+nuT+Yf7iJEeSnBheL5pcTEnSapvZA78TOL5ifj9wtKquAY4O85KkKRmrwJNcDtwM/NmK4T3AoWH6EHDL9kaTJJ3LuHvgfwT8NvC9FWOXVtVJgOF151rfmGRfkuUky2fOnDmvsJKk79uwwJO8DzhdVce28gFVdbCqlqpqaWFhYSs/QpK0hnGeB3498IEkNwFvAn44yV8Cp5LsqqqTSXYBpycZVJL0ShvugVfVR6vq8qpaBG4D/raqfhU4DOwd3rYXuG9iKSVJr3I+14EfAG5McgK4cZiXJE3Jpv5Jtap6EHhwmH4BuGH7I0mSxuGdmJLUlAUuSU1Z4JLUlAUuSU1Z4JLUlAUuSU1Z4JLUlAUuSU1Z4JLUlAUuSU1Z4JLUlAUuSU1Z4JLU1KaeRih1sbj/gTXHnz1w85STSJPjHrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNbVhgSd5U5KvJXk0yZNJfm8YvzjJkSQnhteLJh9XknTWOHvgLwG/UFU/BbwD2J3k3cB+4GhVXQMcHeYlSVOyYYHXyH8MsxcOXwXsAQ4N44eAWyaSUJK0prGOgSfZkeQR4DRwpKoeAi6tqpMAw+vOdb53X5LlJMtnzpzZrtyS9Lo3VoFX1f9W1TuAy4F3JfnJcT+gqg5W1VJVLS0sLGw1pyRplU1dhVJV/w48COwGTiXZBTC8nt72dJKkdY1zFcpCkrcO0z8IvAd4CjgM7B3ethe4b1IhJUmvdsEY79kFHEqyg1Hh31tV9yf5KnBvktuBbwG3TjCnJGmVDQu8qh4Drltj/AXghkmEkiRtzDsxJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJakpC1ySmrLAJampDQs8yRVJvpzkeJInk9w5jF+c5EiSE8PrRZOPK0k6a5w98JeB36yqtwPvBn4jybXAfuBoVV0DHB3mJUlTsmGBV9XJqnp4mP4ucBy4DNgDHBredgi4ZVIhJUmvtqlj4EkWgeuAh4BLq+okjEoe2LnO9+xLspxk+cyZM+eXVpL0/8Yu8CRvBj4HfKSqXhz3+6rqYFUtVdXSwsLCVjJKktYwVoEnuZBRed9TVZ8fhk8l2TUs3wWcnkxESdJaxrkKJcCngeNV9ckViw4De4fpvcB92x9PkrSeC8Z4z/XAh4HHkzwyjP0ucAC4N8ntwLeAWycTUZK0lg0LvKr+Hsg6i2/Y3jiSpHF5J6YkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNWWBS1JTFrgkNXXBrANI41jc/8CsI0hzxz1wSWrKApekpixwSWrKApekpixwSWrKApekpixwSWrKApekpixwSWrKApekpixwSWrKApekpixwSWrKApekpnycrOaKj42VxrfhHniSu5OcTvLEirGLkxxJcmJ4vWiyMSVJq41zCOXPgd2rxvYDR6vqGuDoMC9JmqINC7yqvgJ8Z9XwHuDQMH0IuGWbc0mSNrDVk5iXVtVJgOF15/ZFkiSNY+InMZPsA/YBXHnllZP+OM2Z9U5KPnvg5iknkV57troHfirJLoDh9fR6b6yqg1W1VFVLCwsLW/w4SdJqWy3ww8DeYXovcN/2xJEkjWucywg/A3wVeFuS55LcDhwAbkxyArhxmJckTdGGx8Cr6kPrLLphm7NIkjbBW+klqSkLXJKassAlqSkLXJKassAlqSkLXJKassAlqSkLXJKassAlqSn/STVtik8XlOaHe+CS1JQFLklNWeCS1JQFLklNeRJTwpOz6sk9cElqygKXpKYscElqygKXpKZSVVP7sKWlpVpeXp7a573ebNeJuPV+jl5pVic4PeH6+pPkWFUtrR53D1ySmrLAJakpC1ySmrLAJakp78SU5tR2nUz2pOdrl3vgktSUBS5JTVngktSUx8AlvYLHzLfuXOctJrH+3AOXpKYscElqygKXpKYscElqypOYDW32Bg9PSr0+zOopktM+cafvcw9ckpqywCWpKQtckpqywCWpqTYnMV8LJ+JeC/8N2thm/z/P48nHSX+G2/z2cA9ckpo6rwJPsjvJ00meSbJ/u0JJkja25QJPsgP4Y+C9wLXAh5Jcu13BJEnndj574O8Cnqmqb1bVfwOfBfZsTyxJ0kZSVVv7xuRXgN1V9evD/IeBn6mqO1a9bx+wb5h9G/D0isWXAN/eUoDJMtfmmGtzzDW+ecwE08/1Y1W1sHrwfK5CyRpjr/rToKoOAgfX/AHJclUtnUeGiTDX5phrc8w1vnnMBPOT63wOoTwHXLFi/nLg+fOLI0ka1/kU+D8C1yS5KskbgNuAw9sTS5K0kS0fQqmql5PcAXwR2AHcXVVPbvLHrHloZQ6Ya3PMtTnmGt88ZoI5ybXlk5iSpNnyTkxJasoCl6SmZlLgSf4gyVNJHkvyhSRvXbHso8Ot+U8n+aUp57o1yZNJvpdkacX4YpL/SvLI8PWpecg1LJvZ+lqV42NJ/m3FOrpphlnm8hEPSZ5N8viwfpZnmOPuJKeTPLFi7OIkR5KcGF4vmpNcM9+uklyR5MtJjg+/h3cO4zNfZ1TV1L+AXwQuGKY/AXximL4WeBR4I3AV8A1gxxRzvZ3RzUYPAksrxheBJ2axrjbINdP1tSrjx4DfmtU6WpFjx7AergbeMKyfa2eda8j2LHDJHOT4eeCdK7dp4PeB/cP0/rO/k3OQa+bbFbALeOcw/Rbgn4bfvZmvs5nsgVfVl6rq5WH2HxhdQw6jW/E/W1UvVdU/A88wumV/WrmOV9XTG79zus6Ra6bra075iIcNVNVXgO+sGt4DHBqmDwG3TDUU6+aauao6WVUPD9PfBY4DlzEH62wejoH/GvDXw/RlwL+uWPbcMDYPrkry9SR/l+TnZh1mMG/r647hsNjdM/nr5Mi8rZOVCvhSkmPDIybmyaVVdRJGhQXsnHGeleZhuwJGh1OB64CHmIN1NrF/0CHJ3wA/usaiu6rqvuE9dwEvA/ec/bY13r+t1zmOk2sNJ4Erq+qFJD8N/FWSn6iqF2eca+Lr6xUfdo6MwJ8AHx8+/+PAHzL6w3naprpONun6qno+yU7gSJKnhr1OrW9etiuSvBn4HPCRqnoxWWtTm66JFXhVvedcy5PsBd4H3FDDQSSmcHv+RrnW+Z6XgJeG6WNJvgH8OLBtJ6K2kospP85g3IxJ/hS4f1I5NjC3j3ioqueH19NJvsDocM+8FPipJLuq6mSSXcDpWQcCqKpTZ6dnuV0luZBRed9TVZ8fhme+zmZ1Fcpu4HeAD1TVf65YdBi4Lckbk1wFXAN8bRYZV0qyMDz/nCRXM8r1zdmmAuZofQ0b8FkfBJ5Y770TNpePeEjyQ0necnaa0Yn8Wa2jtRwG9g7Te4H1/tY3VfOwXWW0q/1p4HhVfXLFotmvsxmd1X2G0XHKR4avT61YdhejqwieBt475VwfZLQH9xJwCvjiMP7LwJOMrmh4GHj/POSa9fpalfEvgMeBxxht2LtmmOUmRlcKfIPRIaiZ5FiV6eph+3l02JZmlgv4DKPDgv8zbFe3Az8CHAVODK8Xz0mumW9XwM8yOoTz2IrOumke1pm30ktSU/NwFYokaQsscElqygKXpKYscElqygKXpKYscElqygKXpKb+D14Fdbl+NDIGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(D, bins=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## abundance difference induced by eating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculated the abundance change induced by eating.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "v_abs_max_eat = []\n",
    "for t1, t2 in ma_eat_windows:\n",
    "    v = []\n",
    "    for v_abs_max, t_a, t_b in v_ta_tb:\n",
    "        # check if [t_a, t_b] in eat window [t1, t2]\n",
    "        if (t_a >= t1) and (t_b <=t2 + pd.to_timedelta(window)):\n",
    "            # if v_abs_max is negative value, and t_a is after t2, \n",
    "            # basically the ASV increased due to eating first, and decreased more after that. \n",
    "            # in this situation, we are more interested in how much it increased.\n",
    "            # so we delete the negative value even if t_a and t_b both locate within eat window.\n",
    "            if not ((v_abs_max < 0) and (t_a >= t2)):\n",
    "                v.append(v_abs_max)\n",
    "    \n",
    "    if len(v) > 0:\n",
    "        v_abs_max_eat.append(max(v, key=abs)) # the absolute max value within this window was taken as the abundance difference induced by eating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## calculate zscore"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFkAAAA9CAYAAADGbsXjAAAEfUlEQVR4Ae2aMWvjMBTH/TU8Z83WL5ApS+euHrp06B7IdtsNhk43HRj6BTQeHAVvXQ68laN4KyV4uaGEDKUU8z9k+9myLSmq07iVq0KJI1tW9NNf7z09yYP7OzoB7+gtuAbgII8gAgfZQR6BwAhNOCU7yCMQGKEJp2QHeQQCIzThlNyD/IqMXcDzvD3/MwTssVdbVuAgy6i8c9mEIIsKNFfZO/OUvm5CkHn/npFGZ/C8QyDvkIQLhalYYhUxxOlWClNVODHIwGsSYnYQZI7qH+LVSXuw8gwJCxH43FYvsY43yFVUO+UOcgdI+ZXU3J8R+YbhnIP2L8E2L9La3UK7IbfUJUYDfTjdjuu/qyE3KvexjO6N1Gwx5Cck4Sk8b44gusOOU9vdIQrm7Wmup6m4q4MsONiAIVO8QSy2FnKeRljyWHYZIa2NY45tvIbfssmPYMFM4chE9V+AZa8VGwe5AFE6OA+zMAGh4Tfex/GZQfZXMUziDEuVTIr9CMgUeXwBm1ybi5ZdJPjHc3xNdLFGvK3tlGiCe9eWKhlAfo9o6Qsxa45dyrBa8DIfi/UNMjMGPShNBCEO1hZpHFXv/0Jxcp7d4qqIJrgDmyMIf+Nv/B2zKrnTtdcSmpIisseiU6yuFytE7BeSzCw+ppfbq2TqgQWfDvLRB+nZHQk4OmOIkDOGQJuoNg9Zjv/D7WqhMRcVZD+4xv2O3DKFRB68xRWSulzTyb2DJXEofHBnIRJxVaFpwrZbHcinCJOnpg+7PwiLkGiOc/ZglAxpKrsrItBAppL6kxIwHvxzhg2Ju77vLkwJKCDn2CVXWPBp/Ia8qWmjQ57bv7GpMENaP/O+dVT9kkM+xEwcySZPDLIzEypFDi3vKPnzmYmhHftM9dqQlWaiTO8NywV8pu5+zG8RIKvMxAuy+FvhBB3kYYNUQ67zsxpv7CAfCHlYdUtr0S73zDzxTj0txMhTnnFabt7SDc1nrWTNM9O6VexoLxBc3VZJfU3+2JOfGCrz2CfGGwNfC3L+AHY+l6xgad9O2AkhtatODBVBwolRuuELQaZk1xmi9LkzO0nNAuTqiWZPr3ti6AUbdgnfYEVsJ+Rdipj9rPbb+ktjuYOu1No6p0Gs1ZCb/T5JqncbY+VLyum11ad1kPPsBmueGVyscJ1kyPMN4vUSfvADf3R7bwWQ/hGCkoMOsubEEG3mtnbMO4QB23ZGyHa2U7IUfuoOm5TP9M3BQZBRnU7akwu3S8mVGnsJflKUrw7J9CeLzJTcH8Sq3qQgU4av1yk67yaeZ2tP2+GQafbIbO8UIasUq1K4wHmouWiiC9ksmaK5QBU28RNCK4aU7zlWjo8fbtFukWkdH6lVtNkGJ4Zo0Kfl+Lg0xc5X4ZvRMlcVwpE97oeCPILRnhgqBk5mRoQpZF900f7xb/umW4y87U3l01NfjAxhwusol9UDXuiW1WpoZXJHTBCpn1XdcQkiFRmxnJI/LtUpUrH72q4Vn6WsHeQRBs5BdpBHIDBCE07JDvIIBEZowil5BMj/AQbrEvPT1LhZAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The z-score was computed as ![image.png](attachment:image.png) for each eating time window, where d is the abundance change induced by eating, D̅ is the mean of D, σ (D) is the standard deviation of D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_zscore(observed_value, null_values):\n",
    "    '''calculate z-score'''\n",
    "    from scipy.stats import norm\n",
    "    null_mean = np.mean(null_values)\n",
    "    null_std = np.std(null_values)\n",
    "    z_score = (observed_value - null_mean) / null_std\n",
    "    return z_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.29955509213708276,\n",
       " 0.44852896337395154,\n",
       " 1.2976800294241035,\n",
       " 1.7148068688873364,\n",
       " 1.9159215950571091,\n",
       " 1.7148068688873364,\n",
       " 2.5490605478138018,\n",
       " 1.3125774165477904,\n",
       " 1.6626660139544323,\n",
       " 1.707358175325493,\n",
       " 1.387064352166225,\n",
       " 1.4392052070991288,\n",
       " 0.19527338227127455,\n",
       " 1.4987947555938763,\n",
       " 4.105837502239081,\n",
       " 2.3107023538348117,\n",
       " 1.997857224237387]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "asv_zscores = [calc_zscore(d, D) for d in v_abs_max_eat]\n",
    "asv_zscores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the z-scores of all the eating windows was the final measure on how much this ASV responded to eating."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6210409616970718"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(asv_zscores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# temporal order of eat-increasing bacteria"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we intergated all the eat windows in subject MA together to reveal the temperoral response of the eat-increasing bacteria."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## eat-increasing ASVs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Top 8 ASVs that exist in at least two individuals with highest average z-score were used to visualize the temperoal order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4 of the ASVs were enriched in tongue samples\n",
    "asvs_tongue = ['TACGGAGGGTGCAAGCGTTACTCGGAATCACTGGGCGTAAAGGACGCGTAGGCGGATTATCAAGTCTCTTGTGAAATCCTATGGCTTAACCATAGAACTGCTTGGGAAACTGATAATCTAGAGTGAGGGAGAGGCAGATGGAATTGGTGG',\n",
    " 'TACGTATGTCGCGAGCGTTATCCGGAATTATTGGGCATAAAGGGCATCTAGGCGGATATACAAGTCAGGGGTGAAAACTTAGGGCTCAACTCAAAGCTTGCCTTTGAAACTGTATATCTAGAGTGCTGGAGAGGTGGACGGAACTACACG',\n",
    " 'TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGAATGGCAAGTCTGAAGTGAAATACCCGGGCTCAACCTGGGAACTGCTTTGGAAACTGTTGTTCTAGAGTGTTGGAGAGGTAAGTGGAATTCCTGG',\n",
    " 'TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGGGCTTGTAGGCGGTTGGTCGCGTCTGCCGTGAAATTCTCTGGCTTAACTGGGGGCGTGCGGTGGGTACGGGCTGACTTGAGTGCGGTAGGGGAGACTGGAACTCCTGG']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the other 4 ASVs were enriched in dental plaque samples\n",
    "asvs_plaque = ['TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA',\n",
    " 'TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCTGGGAATTGCGTTCTGAACTGGGTGGCTAGAGTGTGTCAGAGGGGGGTAGAATTCCACG',\n",
    " 'TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTTGTAGGCGGTTGGTCGCGTCTGCTGTGAAAGGCCGGGGCTTAACTCCGGTTTTGCAGTGGGTACGGGCTAACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGG',\n",
    " 'TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAACTCCGCA']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA</th>\n",
       "      <th>TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCTGGGAATTGCGTTCTGAACTGGGTGGCTAGAGTGTGTCAGAGGGGGGTAGAATTCCACG</th>\n",
       "      <th>TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTTGTAGGCGGTTGGTCGCGTCTGCTGTGAAAGGCCGGGGCTTAACTCCGGTTTTGCAGTGGGTACGGGCTAACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGG</th>\n",
       "      <th>TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAACTCCGCA</th>\n",
       "      <th>TACGGAGGGTGCAAGCGTTACTCGGAATCACTGGGCGTAAAGGACGCGTAGGCGGATTATCAAGTCTCTTGTGAAATCCTATGGCTTAACCATAGAACTGCTTGGGAAACTGATAATCTAGAGTGAGGGAGAGGCAGATGGAATTGGTGG</th>\n",
       "      <th>TACGTATGTCGCGAGCGTTATCCGGAATTATTGGGCATAAAGGGCATCTAGGCGGATATACAAGTCAGGGGTGAAAACTTAGGGCTCAACTCAAAGCTTGCCTTTGAAACTGTATATCTAGAGTGCTGGAGAGGTGGACGGAACTACACG</th>\n",
       "      <th>TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGAATGGCAAGTCTGAAGTGAAATACCCGGGCTCAACCTGGGAACTGCTTTGGAAACTGTTGTTCTAGAGTGTTGGAGAGGTAAGTGGAATTCCTGG</th>\n",
       "      <th>TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGGGCTTGTAGGCGGTTGGTCGCGTCTGCCGTGAAATTCTCTGGCTTAACTGGGGGCGTGCGGTGGGTACGGGCTGACTTGAGTGCGGTAGGGGAGACTGGAACTCCTGG</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DATETIME</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2014-09-15 07:00:00</th>\n",
       "      <td>0.039170</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.078339</td>\n",
       "      <td>0.03917</td>\n",
       "      <td>0.039170</td>\n",
       "      <td>0.078339</td>\n",
       "      <td>0.117509</td>\n",
       "      <td>0.039170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-09-15 07:30:00</th>\n",
       "      <td>0.078339</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.274187</td>\n",
       "      <td>0.03917</td>\n",
       "      <td>0.078339</td>\n",
       "      <td>0.195848</td>\n",
       "      <td>0.430866</td>\n",
       "      <td>0.156678</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                0.039170   \n",
       "2014-09-15 07:30:00                                                                                                                                                0.078339   \n",
       "\n",
       "                     TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCTGGGAATTGCGTTCTGAACTGGGTGGCTAGAGTGTGTCAGAGGGGGGTAGAATTCCACG  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                     0.0   \n",
       "2014-09-15 07:30:00                                                                                                                                                     0.0   \n",
       "\n",
       "                     TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTTGTAGGCGGTTGGTCGCGTCTGCTGTGAAAGGCCGGGGCTTAACTCCGGTTTTGCAGTGGGTACGGGCTAACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGG  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                0.078339   \n",
       "2014-09-15 07:30:00                                                                                                                                                0.274187   \n",
       "\n",
       "                     TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAACTCCGCA  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                 0.03917   \n",
       "2014-09-15 07:30:00                                                                                                                                                 0.03917   \n",
       "\n",
       "                     TACGGAGGGTGCAAGCGTTACTCGGAATCACTGGGCGTAAAGGACGCGTAGGCGGATTATCAAGTCTCTTGTGAAATCCTATGGCTTAACCATAGAACTGCTTGGGAAACTGATAATCTAGAGTGAGGGAGAGGCAGATGGAATTGGTGG  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                0.039170   \n",
       "2014-09-15 07:30:00                                                                                                                                                0.078339   \n",
       "\n",
       "                     TACGTATGTCGCGAGCGTTATCCGGAATTATTGGGCATAAAGGGCATCTAGGCGGATATACAAGTCAGGGGTGAAAACTTAGGGCTCAACTCAAAGCTTGCCTTTGAAACTGTATATCTAGAGTGCTGGAGAGGTGGACGGAACTACACG  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                0.078339   \n",
       "2014-09-15 07:30:00                                                                                                                                                0.195848   \n",
       "\n",
       "                     TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGAATGGCAAGTCTGAAGTGAAATACCCGGGCTCAACCTGGGAACTGCTTTGGAAACTGTTGTTCTAGAGTGTTGGAGAGGTAAGTGGAATTCCTGG  \\\n",
       "DATETIME                                                                                                                                                                      \n",
       "2014-09-15 07:00:00                                                                                                                                                0.117509   \n",
       "2014-09-15 07:30:00                                                                                                                                                0.430866   \n",
       "\n",
       "                     TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGGGCTTGTAGGCGGTTGGTCGCGTCTGCCGTGAAATTCTCTGGCTTAACTGGGGGCGTGCGGTGGGTACGGGCTGACTTGAGTGCGGTAGGGGAGACTGGAACTCCTGG  \n",
       "DATETIME                                                                                                                                                                     \n",
       "2014-09-15 07:00:00                                                                                                                                                0.039170  \n",
       "2014-09-15 07:30:00                                                                                                                                                0.156678  "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table = saliva_ma.filter_ids(asvs_plaque + asvs_tongue, axis='f').to_pandas(sample_field='DATETIME')\n",
    "table.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## intergate the eating windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((Timestamp('2014-09-15 10:00:00'), Timestamp('2014-09-15 10:30:00')),\n",
       " (Timestamp('2014-09-15 13:00:00'), Timestamp('2014-09-15 13:30:00')),\n",
       " (Timestamp('2014-09-15 16:00:00'), Timestamp('2014-09-15 16:35:00')),\n",
       " (Timestamp('2014-09-15 19:00:00'), Timestamp('2014-09-15 19:45:00')),\n",
       " (Timestamp('2014-09-16 10:00:00'), Timestamp('2014-09-16 10:30:00')),\n",
       " (Timestamp('2014-09-16 14:30:00'), Timestamp('2014-09-16 15:00:00')),\n",
       " (Timestamp('2014-09-16 19:10:00'), Timestamp('2014-09-16 20:00:00')),\n",
       " (Timestamp('2014-09-17 10:00:00'), Timestamp('2014-09-17 10:30:00')),\n",
       " (Timestamp('2014-09-17 14:00:00'), Timestamp('2014-09-17 14:40:00')),\n",
       " (Timestamp('2014-09-17 19:35:00'), Timestamp('2014-09-17 20:05:00')),\n",
       " (Timestamp('2014-09-18 09:15:00'), Timestamp('2014-09-18 10:00:00')),\n",
       " (Timestamp('2014-09-18 15:00:00'), Timestamp('2014-09-18 15:40:00')),\n",
       " (Timestamp('2014-09-18 19:35:00'), Timestamp('2014-09-18 20:10:00')),\n",
       " (Timestamp('2014-09-19 10:15:00'), Timestamp('2014-09-19 10:30:00')),\n",
       " (Timestamp('2014-09-20 10:00:00'), Timestamp('2014-09-20 10:30:00')),\n",
       " (Timestamp('2014-09-19 14:45:00'), Timestamp('2014-09-19 15:20:00')),\n",
       " (Timestamp('2014-09-19 19:35:00'), Timestamp('2014-09-19 20:10:00')))"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_eat_windows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each tuple within the ma_eat_windows was (eat_start_timepoint, eat_end_timepoint)  \n",
    "\n",
    "Within each eating time window, we set the eat_start_timepoint to be minute 0, then infered the abundance from 0 to 120 minutes using linear transformation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_trans(asv, start, table_reference, window='2H'):\n",
    "    '''estimate the relative abundance of the ASV within \"window\"-long eating window at 1-min interval.\n",
    "    start : timestamp\n",
    "    window : default '2H'\n",
    "    table_reference : inter the relative abundance based this table\n",
    "        index: timestamp, columns: ASV sequence \n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    pd.DataFrame\n",
    "        transformed relative abundance at every elapsed mintute.\n",
    "    '''\n",
    "    index = pd.date_range(start=start, end=start + pd.to_timedelta(window), freq='1min')\n",
    "    df = pd.DataFrame(index=index)\n",
    "    df['elapsed_minute'] = [int((i - df.index[0]).total_seconds() / 60) for i in df.index]\n",
    "    for i in df.index:\n",
    "        if i in table_reference.index:\n",
    "            df.loc[i, 'abund'] = table_reference.loc[i, asv]\n",
    "        else:\n",
    "            before = table_reference.index[table_reference.index < i][-1]\n",
    "            after = table_reference.index[table_reference.index > i][0]\n",
    "            abund = (table_reference.loc[after, asv] - table_reference.loc[before, asv]) * (i - before) / (after - before) + table_reference.loc[before, asv]\n",
    "            df.loc[i, 'abund'] = abund\n",
    "    df.set_index('elapsed_minute', inplace=True)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_asv_median = pd.DataFrame(index=range(120),\n",
    "                             columns=asvs_tongue + asvs_plaque)\n",
    "for asv in asvs_tongue + asvs_plaque:\n",
    "    df_asv = pd.DataFrame()\n",
    "    for start, _ in ma_eat_windows:\n",
    "        df = linear_trans(asv, start, table)\n",
    "        df_asv = pd.concat([df_asv, df], axis=1, sort=False)\n",
    "    df_asv_median[asv] = df_asv.median(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(120, 8)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_asv_median.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TACGGAGGGTGCAAGCGTTACTCGGAATCACTGGGCGTAAAGGACGCGTAGGCGGATTATCAAGTCTCTTGTGAAATCCTATGGCTTAACCATAGAACTGCTTGGGAAACTGATAATCTAGAGTGAGGGAGAGGCAGATGGAATTGGTGG</th>\n",
       "      <th>TACGTATGTCGCGAGCGTTATCCGGAATTATTGGGCATAAAGGGCATCTAGGCGGATATACAAGTCAGGGGTGAAAACTTAGGGCTCAACTCAAAGCTTGCCTTTGAAACTGTATATCTAGAGTGCTGGAGAGGTGGACGGAACTACACG</th>\n",
       "      <th>TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGAATGGCAAGTCTGAAGTGAAATACCCGGGCTCAACCTGGGAACTGCTTTGGAAACTGTTGTTCTAGAGTGTTGGAGAGGTAAGTGGAATTCCTGG</th>\n",
       "      <th>TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGGGCTTGTAGGCGGTTGGTCGCGTCTGCCGTGAAATTCTCTGGCTTAACTGGGGGCGTGCGGTGGGTACGGGCTGACTTGAGTGCGGTAGGGGAGACTGGAACTCCTGG</th>\n",
       "      <th>TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA</th>\n",
       "      <th>TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCTGGGAATTGCGTTCTGAACTGGGTGGCTAGAGTGTGTCAGAGGGGGGTAGAATTCCACG</th>\n",
       "      <th>TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTTGTAGGCGGTTGGTCGCGTCTGCTGTGAAAGGCCGGGGCTTAACTCCGGTTTTGCAGTGGGTACGGGCTAACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGG</th>\n",
       "      <th>TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAACTCCGCA</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.117509</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.235018</td>\n",
       "      <td>0.117509</td>\n",
       "      <td>0.430866</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.509205</td>\n",
       "      <td>0.078339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.123602</td>\n",
       "      <td>0.006715</td>\n",
       "      <td>0.274187</td>\n",
       "      <td>0.132306</td>\n",
       "      <td>0.420420</td>\n",
       "      <td>0.005223</td>\n",
       "      <td>0.500065</td>\n",
       "      <td>0.088784</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   TACGGAGGGTGCAAGCGTTACTCGGAATCACTGGGCGTAAAGGACGCGTAGGCGGATTATCAAGTCTCTTGTGAAATCCTATGGCTTAACCATAGAACTGCTTGGGAAACTGATAATCTAGAGTGAGGGAGAGGCAGATGGAATTGGTGG  \\\n",
       "0                                                                                                                                                0.117509   \n",
       "1                                                                                                                                                0.123602   \n",
       "\n",
       "   TACGTATGTCGCGAGCGTTATCCGGAATTATTGGGCATAAAGGGCATCTAGGCGGATATACAAGTCAGGGGTGAAAACTTAGGGCTCAACTCAAAGCTTGCCTTTGAAACTGTATATCTAGAGTGCTGGAGAGGTGGACGGAACTACACG  \\\n",
       "0                                                                                                                                                0.000000   \n",
       "1                                                                                                                                                0.006715   \n",
       "\n",
       "   TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGAATGGCAAGTCTGAAGTGAAATACCCGGGCTCAACCTGGGAACTGCTTTGGAAACTGTTGTTCTAGAGTGTTGGAGAGGTAAGTGGAATTCCTGG  \\\n",
       "0                                                                                                                                                0.235018   \n",
       "1                                                                                                                                                0.274187   \n",
       "\n",
       "   TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGGGCTTGTAGGCGGTTGGTCGCGTCTGCCGTGAAATTCTCTGGCTTAACTGGGGGCGTGCGGTGGGTACGGGCTGACTTGAGTGCGGTAGGGGAGACTGGAACTCCTGG  \\\n",
       "0                                                                                                                                                0.117509   \n",
       "1                                                                                                                                                0.132306   \n",
       "\n",
       "   TACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAATTCCGCA  \\\n",
       "0                                                                                                                                                0.430866   \n",
       "1                                                                                                                                                0.420420   \n",
       "\n",
       "   TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGGGCGCAGACGGTTACTTAAGCAGGATGTGAAATCCCCGGGCTCAACCTGGGAATTGCGTTCTGAACTGGGTGGCTAGAGTGTGTCAGAGGGGGGTAGAATTCCACG  \\\n",
       "0                                                                                                                                                0.000000   \n",
       "1                                                                                                                                                0.005223   \n",
       "\n",
       "   TACGTAGGGCGCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTTGTAGGCGGTTGGTCGCGTCTGCTGTGAAAGGCCGGGGCTTAACTCCGGTTTTGCAGTGGGTACGGGCTAACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGG  \\\n",
       "0                                                                                                                                                0.509205   \n",
       "1                                                                                                                                                0.500065   \n",
       "\n",
       "   TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGAGTGCGCAGGCGGTTTTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTGGTGACTGCAAGGCTAGAGTGTGTCAGAGGGAGGTGGAACTCCGCA  \n",
       "0                                                                                                                                                0.078339  \n",
       "1                                                                                                                                                0.088784  "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_asv_median.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moving_average(series, half=10, shift=10):\n",
    "    '''moving average.\n",
    "    series : pd.Series\n",
    "    half : half length of moving window\n",
    "    shift : shift length\n",
    "    \n",
    "    returns\n",
    "    -------\n",
    "    pd.Series\n",
    "    '''\n",
    "    index = range(int(series.index[0]), int(series.index[-1]), shift)\n",
    "    s_m = pd.Series(index=index)\n",
    "    for i in s_m.index:\n",
    "        index_inc = [sid for sid in series.index if (sid >= i - half) & (sid < i + half)]\n",
    "        s_m.loc[i] = series.reindex(index_inc).mean()\n",
    "    return s_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Normalized relative abundance')"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAE+CAYAAACA8heHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeVxU5f4H8M+ZhWHfURZZBQYRQVkGGJZETXFDEwoSUqTrkhq/1JvrLbebpRUqqdey615hUZpiZW7XBM0tU3EZwUAUAdn3gVme3x/jECLLjIBsz/v16lVz5sw5z0y8zvc85zzP5zCEEFAURVF9F6urG0BRFEV1LVoIKIqi+jhaCCiKovo4WggoiqL6OFoIKIqi+jhaCCiKovo4Tlc3oL3++OOPyRwOZykhxLyr20L1SYQQsn3o0KHru7ohFPW8enwhYLPZnw4cOLBWW1u7iGGYrm4O1cdIpVL2zZs3ZwOghYDqsXp8IQDA1tHREXd1I6i+icPhyEAvsVI9HP0D7kRpaWlaQqHQ2cfHhx8UFOR08+ZNXmfuLzw83C49Pb3ZfSQmJpokJCSYPs92ExMTTWQyWfsa94IlJiaanD17Vrur20FRPUFv6BF0S2KxmHn77bdtDx06lGljYyPNycnhlJaWsru6Xc9j//79pm+99VYxm91682UyGdpa50WJj48v7uo2UFRPQXsEKqqtrWVGjBjhGBQU5DRhwgSHxMREk9bWP3nypE5gYGCFjY2NFABsbGykHh4edcnJyfo+Pj58Nze3QVu2bDEBgIULF1pGRUXZCoVC56lTp9q8++67Fh4eHi7//Oc/LZTvh4eH2wmFQudXX33VDgBGjhw5sLCwkA0AcXFx1o3PfouLi9khISGO3t7e/NjYWGvl8hMnTugHBgY6jRw5cqBYLGays7O5vr6+zl5eXvyYmBgbQHEwj4yMtPXx8eEHBwc7nT59Wvv27dtaQqGQv3XrVuP09HReQECAk4+PD3/x4sUWgKInMm3aNJvg4GCnnJwczpIlS8wBYP78+Vaenp4uvr6+ztnZ2dy2fuOcnBxOcHCwk0Ag4M+bN88KAFasWGHu6enp4ufn55yRkaEBAHw+3/WVV16x4/P5rufOndMCgPXr15t5eHi4+Pr6Ol+7do23cOFCy0OHDukdP35cx93d3cXPz89506ZNJikpKXrx8fGWgKLXkJiYaJKens4bNmyYi6+vr/PSpUvpoAOqz+k1PQLPtae8OmI7f7w34kpzy/fv32/o5+dXtW7duvzo6GibtraTm5vLtbCwkDRdPmbMmKqIiAiRRCKBn58ff/78+cUAMGjQoNqkpKT7AQEBTqGhoRUff/xxnpub26BPPvkkDwAcHR3F33//fXZ0dLTNyZMndcLDw0v37dtn9H//939FN2/e1AoKCqrZtGkTAGDTpk2mERERJfPmzSuJjIy0PXXqlA4A6Ovry3766ae/VqxYYb53717DadOmlaWmpt7lcrmYNGmS/Y0bN3hXr17VMjMzkx44cECkPMMfNGhQbVpamojL5WL8+PEOe/bsyXZ0dJRMnDjR/t69e1wACAgIqNq7d28OAKxfvz4fAC5cuKB76dKlO2w2G3K5vM3ffuXKlRbvvPNOwZQpUypkMhlycnI4Z86c0fvjjz/uHDt2THfVqlXmX331VU5RURH3q6++up+amqr93//+19TW1jbv4MGDRleuXLnD4XDQ+DLWkSNHDNatW5c7YcKESrlcjp9++kmv6X5//fVXvRkzZhTGx8cXq9JOiuptaI9ARVlZWTwPD49aABg6dGhNW+tbWVlJHj16pNF0eVpamrZQKHQWCoX8jIwMLeVyd3d3MQCYm5tLhg0bVgsA2tracqlUCgDw9PRs2LdIJOJFR0eXHT161OCXX37RFQqFVY33ce/ePZ6vr28NAHh5edWIRCJe43Z7enrWZGZmahYUFHDGjh07UCAQ8C9fvqybk5PDFYlEmgEBAVUAmr3Mc+/ePc2pU6c6CAQCfkZGhtb9+/c1AMDf37+66bqLFi3KDw8Pt4+Li7Ouqqpq82/t3r17miNHjmzYd0ZGBm/w4MG1ABAUFFSdlZXFAwBbW9s6bW1tYmtrKykvL2eLRCKeu7t7DYejOK9p3O6FCxcWJiUlGU2ePNn+zJkz2gzDNMTtKpN3Y2NjS65fv64VFhZmn5ycrN9WOymqt+k1PYKWzuQ7ir29fd3169e1IiMjy69fv67t4+PzzIGvsZEjR1a/9957A3Jycjg2NjbS3NxcTlFREXvDhg3mu3btyrazs5M4ODi4KddvfIBqPAxWebD6888/tSIjI8uvXbumHRsbW2xgYCDX09OTb968uf/69etzG+/bwcGh7vfff9fx9vYWX7lyRXvWrFlFN2/e1Lx+/bo2AFy9elVr4MCBdTt37jQOCwsri4+PLw4LC7MnhDAuLi7i8+fP67z++uvlyh4Bh8ORS6VShsvlkoEDB4q3bdv2wNbWViKVSsEwDLZu3WrGYimO8xKJBPn5+Rxra2vphAkTKqKiosqXLl1qfuDAAYOZM2eWtvabOTo6ik+fPq0zefLkSplMBicnp7r09HQtADh79qyOnZ1dXdPfihACFxeXuhs3bmgr29u4R2Bqairdv39/TnZ2NnfatGl2H3300cP8/HwuANy4cUPL3d29VkNDg3z55ZcPxWIx4+3t7fLaa69VtNZOiuptek0h6GwxMTFl48ePHxgYGOiko6Mj53K5rT7IQVNTk3z22Wf3o6KiHCQSCcPj8eRffPHF/bCwsNKwsDDHwYMH1+jr66s8FOfOnTua/v7+zpaWlvWjRo2qBoCoqKiS1atXWw0bNuyp4bPvvPNOUXh4uP2uXbtMXV1da0eOHFl98+ZNzZKSEnZAQIATj8cjR48evffHH39ozpgxw/7w4cOGys9OnTq17MiRIwbe3t58HR0d2ZkzZzJDQ0PLx4wZ4xgXF1e4fv363DfeeMOuvr6e4XA45OjRo/ca71skEvH+/e9/myclJd0PDQ11rK2tZQHA3Llzn1qvOatWrcqbOnWq/bp16ywFAkHVli1bcoODgyuHDRvmwuVyyf79+7Oa+5ylpaV00qRJpZ6eni6ampryL7744r7yvYSEBLPDhw8b1dTUsBYsWJDv6+tbm5+fr/HSSy85GhkZyQDUfvPNN4bbt2/vJxaLWa+99hq9yUz1OUxPfzDNtWvXsj08PIpexL4kEgm4XC6io6NtZsyYUaw8IHe2hQsXWgYHB1dOnjy5svHyAwcOGKSnp2uuXbu24EW0QxW7d+82NDY2loWFhVW2vXbvcO3aNVMPDw+7rm4HRT0v2iNQQ0hIiFNNTQ3bzs5OvGbNGsvly5c3XPfW1NSUi8Xihtc7duzI9vDwqOustuzevdvws88+63/06NHMztrH84iNjS1rbvnRo0d1V69ebaV8zTAMIYQ0XAOzsbGpS05Ozn4BTaQoqgnaI6CodqI9Aqqno6OGKIqi+jhaCCiKovo4WggoiqL6OFoIOolIJNKYNGmS/fN8dt++fYYFBQUthvakpKTo3bp165nJagAwffp06+aWN45WUDp37pzWxo0bnyuIrjMp4yGUr728vPiAIspixIgRjsrl06ZNsxEIBPy1a9f2EwgEfCsrqyF8Pt9VIBDwjx49qtt4m6WlpSx/f39nb29vfkhIiGNpaSkLAEaMGOHo5eXF9/f3d1bOkl6yZIl5v3793Jv+XhTVW9FC0A0dOnTIsLCwsMURXadOndK7e/fuMymjMpkMe/bseaDqfoRCYe2CBQt61I32srIyTmlpKUsulyMvL48LAO+9997jixcviiIjI4vXr1//4OLFi6Lx48c/Nduax+ORpKSkrMuXL4smTJhQtm3bNlMA2LZtW86VK1dEixcvzv/oo4/6A8D8+fOL/vvf/zY7Z4GieiNaCFSkbuhcc06dOqUjEAj4np6eLps3bzYBAIFAwJ81a9YADw8Pl08++cQ0IyND48yZMwZTp061f++99/pnZGRo+Pn5OXt6erqsWLHCXCwWMwcOHDBZsmSJ9cyZMwckJiaajB8/3mH48OGOFy5c0FKePf/66686ysC3HTt2GAHAzZs3tUaMGOHo6enpUl5ezmrcSxg3bpyDj48PPyAgwKmkpIQF/N272Lt3r+GQIUMG+fn5OR84cMCgre/ZNOhNJBJp+Pr6Oo8cOXKgh4eHy507d5rtzaji5ZdfLktOTjY4ffq0jre3t8rzOJSRFADA4XAIm80mAODi4lLfaBkAwNraWkofckT1Jb1mHsH8g7c7JHRuyyuDOiR0rjnvv/++5S+//JJhaGgoDwwMdJ49e3YJAERERJRu27btoY+Pj8v8+fOLX3rppfLVq1fnubm51b3xxhs2q1atehQaGloVGBjoNHPmzOLIyMhi5QSzxMREEwMDA9nRo0f/aryvZcuWDTh69GimhYWFVCaT4eeff9YDgFOnTmUuWbLE/MiRI/qGhoYNM5sPHDiQraenJ09ISDDdtWuX8aJFi4qUvYvvv//eKDk5+R6fz69XJZStadBbRkaGRmFhITc1NfXu2bNnddasWWP+9ddf57S2jSVLllivW7dOBgCFhYUNyaVhYWHln376af/+/ftLZs2aVXzy5Mk2C1Nj5eXlrF27dpmdOHEiQ7lMKpXiww8/tNixY8f91j5LUb0V7RGoSN3QuebcuXNHOzQ01MnPz49fWFjIzcvL4wCAv79/LYfDwYABA+pyc3OfKs7Z2dk8oVBYAwDu7u41d+/efeZs2tPTs9kzYwsLCynwdwibq6urGFAE4jV+NoJUKsXcuXMHeHt787/44ot+jx49eioyetWqVXnvv/++RXh4uJ0qD9dpGvQGAC4uLrVcLhf+/v412dnZmm1tQ3mJ5+LFiyIzM7OGFFcDAwN5fX09SyQSaTWN1miLXC7H66+/brdmzZpcU1PThiI4a9Ys66lTpxYPHjy40yYAUlR31mt6BC2dyXcUdUPnmuPq6lqTkpJyT19fX15XV8fweDwCABcuXNAaPnx49cOHD3lWVlZSLpdLlMFptra2dWlpadpjx46tun79uva77777+OTJk3KZTNZw7UIZ+NYYwzDIz89nm5uby5Tbai55EwDOnz+vXVNTw7p8+bLo008/Nc3NzdUAgKysLK69vb3Eycmp7sCBA/ePHz+u89FHH/X/6quvWj2bbxr0tmPHjvsikUhLKpXi999/11KGxz2vmJiY4sePH6v9t7tgwQJLf3//qsbxFxs3bjRlGAbKOHCK6ot6TSHobOqGzgHAuXPn9IRCoTMAxMbGFq5evTp3zJgxjnK5HIaGhrJjx47dA4CkpCSjhQsXWsfExBRramqSMWPGVMyZM8f2lVdeKX3//ffzY2Ji7JcvX86EhoaW2dvbS0aOHFm5YsWKAefPn9exsbGpb27f69atezh27FgnDQ0N+T/+8Y9CZe+gOe7u7uLs7GxeUFCQk5WVVb2lpaUEACIjIx0uXrwoevfddy0vX76sW1NTw1q/fv3Dtr5306A3ADAxMZGOHj3asbi4mPP1119nAYp7EOrc3FaKiooqV/cz2dnZ3G3btpkPGzas+siRI0bh4eElS5YsKVyyZInNkCFDqgUCAT8gIKBy48aNjzZu3Gj65Zdfmj25Mc3Zt29fq4WPono6GjGhhs4InRMIBHzlQ1+6k/z8fPb69ev7b9y48VF7tyUSiTQWL15s9eOPP/bKkTg0YoLq6WiPQA3dKXSus5mbm8uaKwLFxcXssWPHNozlr6ysZAOAnp5ewzX3n3/+OdPExESliO2mYXQdFT4XERFhl5OT03A/Y+XKlblNh5RSFKVAewQU1U60R0D1dHTUEEVRVB9HCwFFUVQfRwsBRVFUH0cLQSc6e/astq+vr7NAIOBHRETYVVZWPvN7Jycn6yclJRm0J6SutRC65rQUTPeidZd2UFRfRwtBJxGLxUx8fLzNDz/88NfFixdFL7/8csWSJUssGq8jk8kQERFR8Tzj4htrKYSuOeoG03Wm7tIOiurraCFQkbqhcydPntQJDAystLKykgLA7NmzSy5duqQLKOYOzJkzZ8CUKVPsExMTTRISEkwBID8/n9s0mK25MLj169ebeXh4uPj6+jpfu3aN1ziErqqqipk4caK9n5+f8/jx4x3q6uqYloLplP9Wtkn577ffftvKzc1t0Keffmo6adIkez6f7/r999/rA8Dy5cvNs7KyuOnp6bxhw4a5+Pr6Oi9dutS8rd+voKCALRAI+L6+vs6xsbHWjfcfHh5uN2PGDGsvLy/+okWLLJTtkEgkT7Vt/vz5Vsogvezs7O418YKierBeM49g16WcDgmdm+Fj0yGhc7m5uVxLS8unZv2yWKyGsboRERGlo0aNqm5cUJoLZmsaBhcVFVV28OBBoytXrtzhcDiQyWRoHEL3wQcf9JswYUL57NmzS9avX2+2e/duIwBoLpiuJVFRUSUbNmx4ZGFh4XHjxo10mUzGxMXF2YaHh1esW7cuHwASEhJMZ8yYURgfH1+sShDd+fPntQMDAysTEhIeNbd+SEhI5a5dux64u7u7AMhrbhsXLlzQvXTp0h02mw1V9klRlGpoj0BF6obOWVlZSR49evTUdXu5XN6QDxQYGPjMNpoGszUXBicSiXju7u41HI6ihisD5ZRu376tuXXr1n4CgYD/9ddfmygzeVoKpnvSrqdee3t7i7W0tIiDg4PY2tpaamdnJykvL39qR7GxsSXXr1/XCgsLs09OTtZv6/cYO3ZslVwux8SJE+3/85//GDd9f9iwYbUAoKmpSYDmc5EWLVqUHx4ebh8XF2ddVVVF/3YpqoP0mh5BS2fyHUXd0LkRI0ZUr1ixYkBubi7HyspK+vnnnxt7eXk1fEaZh99Y02C25sLgXFxc6m7cuKEtk8nAZrMhk8nA5XIbQuj4fL541KhRFbGxsWUAUFdXx3z++efGzQXTEUKY2tpa5vr160+lgSp7Lo0Pxkq5ubkcU1NTmYaGBvnyyy8fisVixtvb2+W1116raO33kMlk2LRp0yMAcHFxcZ03b15J4/eb5v/r6enJ7t+/z9XR0SHKGOoJEyZUREVFlS9dutT8wIEDBjNnzixtbZ8URamm1xSCzqZu6JyWlhbZvHlzziuvvDKQEAJra+u63bt3txpe1jSYbcCAAZKmYXCWlpbSSZMmlXp6erpoamrKv/jii/uNQ+hWrVqVHx0dbff555/3I4Rg3bp1uS3tLzo6usjHx8dl5MiRKt+snjt3rvWmTZsenjhxQnf79u39xGIx67XXXmszufPMmTM6//rXv6wkEgkTHBzcatEAgJkzZxZNmDDBSSAQVJmYmEgAIDQ01LG2tpb1pB33VG0zRVGtoxETauiM0LmeJiYmxmb//v00jbMRGjFB9XS0EKghMDCwIXSuqKiI2zhkrreFzqmracgbwzCEENJwvac3h77RQkD1dPTSkBpSU1Mz2l6rb+qIxFCKoroGHXlBURTVx9FCQFEU1cfRQkBRFNXH0ULQyYRCofPy5ctbjGAoKipi79mzx1D5mgaxPatxFEZKSopefHy8JQAwDOO1a9cuI0AxX8LAwGBoYmKiSXBwsJNAIODr6ekNFQgE/ODgYKem20xLS9NydnZ2tbKyGqJclpycrC8QCPgCgYBvZmbmvm/fPsM7d+5oeHl58b29vfkTJ060l0pbfPQzRfVYtBB0ooKCArahoaE0LS1Nt6V1iouL2T/88IOR8jUNYlOdk5NT7dGjRw0A4OjRo3rW1tZ1APDbb79lXLx4UeTs7Fx78eJF0W+//fbMTX5XV9e6y5cv3zE3N2+IAYmIiKi4ePGi6OLFiyILC4v6sLCwChMTE9mxY8cyL1++LLKzs6v/9ttvDV7cN6SoF4MWAhWpGzoHAElJSYYRERGlFhYWkqysLC4AXL16VVN51rl27dp+n332mVlqaqq+QCDgP3r0iKM8+xUIBPx58+ZZubm5Ddq4caMpABw6dEjPw8PDxcPDw+XQoUN6yvVaCom7f/8+d8KECQ6AYg6En5+fM/BsaN2jR484I0aMcPT19XWOiYmxAYAPP/ywYZ3U1FTtc+fOaSnbMWXKFDsfHx++r6+vs0zW9qOJm4bFLVy40DI8PNxOKBQ6v/rqq3bP8/8DUMw+FovFLLFYzPzwww+G48ePV3mmsZGRkVxfX7/ZwKJbt25pmJqaSg0MDORmZmYyU1NTGQBwOBzSNNKDonqDXjN89HxmWYeEzvk7GnZI6BwAHD9+3ODrr7/ONjAwkCUlJRkuW7ascPHixVaff/75/WHDhollMhkyMzM1Hjx4oPHjjz9mNf381KlTSz755JNHQUFBzgsWLChau3at5enTp+8CwMiRI50nT558B2g5JO7UqVOZNTU1rNLSUtapU6d0hw8fXpGbm8tpGlo3Z86cAcuXL88bNWpU9VtvvWV14sQJnZSUFMO0tDSRrq4ukcvlYLFYEAqFtXV1dUxeXp7GpUuXRMrlbWkuLM7R0VH8/fffZ0dHR9ucPHlSZ+TIkS1OzisuLuYqE0grKirYo0ePLlO+N3z48Ioff/xRv6CggCsQCDpkgl9SUpLRpEmTnioq2dnZ3DNnzuhv2LDhUUfsg6K6E9ojUJG6oXPl5eWsq1ev6owfP37gp59+2v/YsWMGAFBSUsIZNmyYGHg2MK4pZfhb44OtsbGx3NjYWN44ybS1kLhJkyaVfvPNN4ZJSUnG06ZNK2kutO7u3buay5YtGyAQCPhnz57Vf/jwIXf16tWPpk+fbvv666/b5ubmNpww8Hg8Eh0dXTxp0iT7//u//7NSpUfQXFicp6dnw28pEolafZaCiYmJRHnJ5pNPPnnq0llkZGTZ6tWrLX19fTtsstrPP/9sGBkZ2VBsamtrmalTp9rt2LEjm8ul6ddU79NregQtncl3FHVD57777juDf/7zn3kLFiwoAoDJkyfbFxQUsE1MTKTXr1/nubu718lkMmhoaJCWIpUbH+wBRUic8pkEjZNMWwuJe+ONN0ojIiIcJBIJ4+rqWm9oaChvGlrn6OhYFxsbWxwUFFQDKC4j1dXVMaNGjcrevn278fbt202XLVtWUF1dzTIxMZHNnDmzZP78+cWvv/667ZkzZ3RGjBjR6m/RNCwOAP7880+tyMjI8mvXrmnHxsa2mVXUEmtra2lwcHBlTExM6a+//qr3vNtRysnJ4XC5XLm5uXlDhYuOjradM2dOoZeXl7i926eo7oj2CFQUExNTdv78ed3AwECnx48fc9sKnfvxxx8NR48eXal8LRQKq7799lvD9evX57755pt2AoGAv27dun7W1taS0tJSTmhoqENBQUGrXYQVK1Y8CgkJcR4+fLjzv/71L5UuURgbG8t5PJ589OjR5QDQOLTO19fXOT09nbd27dq8NWvWWPj5+TkLhULne/fuaUybNs3W29ubv23btn5TpkwpO3HihG5CQoJZWVkZOyAgwNnT09MlNzdXw8fHp7atNoSGhjp6eXnxjx8/bjBmzJhKALhz546mv7+/c21tLWvUqFHVje9BqGv79u0PXVxc6tte82+ZmZlcoVDofPfuXS2hUOgsEok0AODAgQOGEyZMaOgNnDhxQufYsWNGW7Zs6S8QCPh79+41bHmrFNUz0awhNfTU0LmJEyfab9my5aGtra3kebexevXqfuHh4eVubm7tzk9auHChpfJBOu3dVndAs4aonq7XXBp6EUJCQhpC59asWWO5fPnybh869/rrr9uamZlJ21MEAGDlypWPm1u+Y8cOox07dvRTvi4qKuKYmpo2DLb38fGp2rp1a4tR2E3NmzfPSvlITwCYOXPm4/Y+d+DatWu8mTNn2ilfa2pqypsbUkpRfRXtEVBUO9EeAdXT0XsEFEVRfRwtBBRFUX0cLQQURVF9HC0EnUQkEmlMmjTJ/nk+m5iYaKLKRK2W5OTkcJYsWdJi0F1nUWe/CxcutFTGZCglJyfrJyUlGTT+7WgIH0V1PloIuqH9+/ebNp1kpk5hsLGxka5fvz6/o9v1PPtVp90REREVUVFR5Y2X0RA+iup8tBCo6HlC55patmyZuY+PD9/d3d0lLS1NC3g6YlkgEPBPnz6tffv2bS2hUMjfunWrcXh4uN20adNsgoODnSQSCSZOnGivjESWSCRISUnRGzFihGNwcLCTn5+fc0FBAbvxGfWMGTOsfXx8+F5eXvyMjAwNAIiLi7OWSqU4fvy4jru7u4ufn5/zpk2b2vw+IpFIw8/Pz3ncuHEOgwYNct2zZ49hQECAk5eXF7+iooLVeL8eHh4u0dHRNrNnzx6QnJys7+Pjw3dzcxu0ZcuWhv3s27fPpHHwXGJioklCQsJTk8qUv0/T4Dp1f3uKolrWa+YR5L2m3yGhcxbfVnRY6FxTy5cvf/zhhx/mp6en85YvX255+PDhZ4LmQkJCagYNGlSblpYm4nK5OHXqlH5AQEDV3r17c3bu3Gnk4uIiPnLkSNbixYstdu/ebWRhYSEVi8XMuXPnMnbs2GG0efNms+nTp5cot5eYmJirp6cnP3TokF5iYqLZZ599lrtz584HAHDkyBGDdevW5U6YMKGypZiLpsrLyzlpaWl3d+zYYbx3716TtLS0jKVLl5ofPHhQ38/PryGDqbS0lLNmzZq8gQMHSiorK1kRERGiJwmo/Pnz5xcDzwbPtbbf5oLrKIrqGLRHoCJ1Q+eas337dmNvb29+XFycXUFBwVNnta0d3Pz9/asBIDMzk+fl5VUDAAKBoDozM5MHAO7u7rVPltXcu3fvqQC3lStXmnt5efFXrlxplZeX99Q+Fy5cWJiUlGQ0efJk+zNnzmir8h2cnJxq2Ww2BgwYIHF1dRUDgKWlpaSkpOSpkwoTExPpwIEDJQCQlpamLRQKnYVCIT8jI0NLuY46wXPNBddRFNUxek2PoKUz+Y6ibuhcc7788st+txR4b775ph2gCJKrra1lrl+/rqlcj8PhyKVSKaPMM1Kmjzo4ONRdvnxZOyoqqvzixYs6jo6OYgBIT0/XAoBLly5pOzg4NMxmzs/PZ6empupduXJFdPDgQf39+/cbA0BWVhbX3t5eYmpqKt2/f39OdnY2d9q0aXYhIccgLtgAACAASURBVCFtzrZlGKbxfzfMRmw6MbFxYN6GDRvMd+3alW1nZydxcHBwUy5vGjx369YtTbSgaXBde2cbUxT1t15TCDpbTExM2fjx4wcGBgY66ejoyNsKnQOAc+fO6QmFQmcAiI2NLRw2bFi1QCDg+/v7N0QmR0dHF/n4+LiMHDmy4SZpaGho+ZgxYxzj4uIKG2/vjTfeKJsyZYq9t7c3v3///pK1a9fmHz9+XJfL5ZKgoCCnuro61uHDhzOVEdRmZmYyHR0dmZ+fn7Orq2tDONzUqVPtU1NT7yYkJJgdPnzYqKamhrVgwYJOu7kcFhZWGhYW5jh48OAafX39hrvHyuA5S0vL+lGjRlW3VghCQ0Mda2trWQAwd+7ce53VVorqi2jEhBq6Y+hcSkqK3q+//qqXmJioUhqpTCZDXFycNR2N03FoxATV09EegRp6YuhcU2w2u8UhmcHBwU6Nv0NP+U4URbUP7RFQVDvRHgHV09HRFxRFUX0cLQQURVF9HC0EFEVRfRwtBJ1EJBJpsFgsrxs3bvCA5kPWlNrzvN6mlMFtHbGtjqCMiFDGWlAU1f3QUUOdyN7eXvzxxx/337t3b05r6wmFwlqhUNjmQ+DbIpPJEBERUdHe7TzvvtlsdovvK2MtKIrqfmiPQEXPEzrn5uZWk5WVxSsqKmo4QsrlckRHR9v4+fk5Dx8+3LGwsJCdkpKiFx8fb1lQUMAWCAR8X19f59jYWGuxWMyEhIQ4+vr6Oo8bN84BAE6dOqUjEAj4np6eLps3bzYBFGF1c+bMGTBlyhR7ZXBbXV0d4+/v7+zt7c0fM2bMQOXZuDLW+cMPPzTz8PBw8fX1dU5NTW0zXqK4uJgdEhLi6O3tzY+NjbUGFCFx48ePdxg+fLjjhQsXtMaNG+fg4+PDDwgIcCopKXnqb0sgEPAlEona+6UoqvP1mh6BPCGoQ0LnWAvPdmjoXFxcXNHmzZsbLvskJSUZ2NjY1H/11Vc53377rX5CQoJZQEBANQCcP39eOzAwsDIhIeGRXC7HrVu3eMbGxtLTp09nK7OI3n//fctffvklw9DQUB4YGOg8e/bsEgCIiIgoHTVqVLWyQHG5XHL8+PEMXV1dEh8fb3nkyBH9V155pUI5hyAlJcUwLS1NpKurS1QJcdu0aZNpREREybx580oiIyNtT506pQMABgYGsqNHj/4FAAcOHMjW09OTJyQkmO7atct40aJFzwzrVXe/FEV1vl5TCDrb84bOvfHGG6V+fn58oVBYCQC3bt3S/OGHH4xPnjxpIJPJ4O3t3TA7eezYsVX/+9//9CZOnGgfGhpaPm/evBI3N7fasLAwe09Pz5pVq1YV3LlzRzs0NNQJUCR85uXlcQAgMDDwqTZVVlayYmJibPPz8zWKioo4Tk5OT00EW7169aPp06fbamhokA0bNuRaW1u3egH/3r17vEmTJpUDgJeXV0NInKenZzUASKVSzJ07d8Dt27e1qqqq2OPHj282C0jd/VIU1fl6TSFo6Uy+ozxv6ByHw8G4cePKvvrqK9MRI0ZUuri41EVFRRWvXr26AADq6uqY48eP6wKK6+ybNm16BAAuLi6ucXFxpe+//34Bm81GQECA05tvvlns6upak5KSck9fX19eV1fH8Hg8AgBsNvupmYEHDx7Ud3Jyqjty5EjW22+/baWcOKgMnPPz86sZNWpU9vbt2423b99u+sEHH7SaNeTg4FD3+++/63h7e4uvXLmiPWvWrKKbN29qKgPxzp8/r11TU8O6fPmy6NNPPzXNzc3VaG476u6XoqjOR+8RqCgmJqbs/PnzuoGBgU6PHz/mqhI6pxQfH1+Ul5enAQBTp04ty87O1vDz83P28/NzTk5O1leud+bMGR0vLy++u7u7S3BwcEVGRoaGQCDgDx061MXExERqZWUlXb16de6YMWMcfX19ncPCwhxa2mdQUFD1L7/8YhASEuJ4//59DUCRlRQdHW0HANOmTbP19vbmb9u2rd+UKVPK2voO77zzTtF3331n5OXlxdfQ0CAjR458qhC6u7uLs7OzeUFBQU6XLl1q8dkC6u6XoqjORyMm1NAdQ+fUcfbsWe0rV65ovfPOO8Vd3ZbehEZMUD0dLQRqCAwMbAidKyoq4vaWgLZr167xZs6caad83fS7aGpqyn/77bc2n1XQV9FCQPV0veYewYuQmpraKw+GHh4edRcvXhR1dTsoiuoa9B4BRVFUH0cLAUVRVB9HCwFFUVQfRwtBJxGJRBomJiYevr6+zt7e3vz09HReS+vu27fPsKCggA00H06XmJhocvbs2W4Rx5CTk8NZsmSJeVe3g6KojkMLQScKDAysuHDhwt0PPvggNzEx0ayl9Q4dOmRYWFjY4o37+Pj44qCgIJVnM3cWuVyOAQMGSNevX08ngVFUL0ILgYqeJ3ROqaSkpCF07vz581qenp4uHh4eLtu2bTPOyMjQOHPmjMHUqVPt33vvvf4A8M033xgLhULnyMhIW+DvXkJ2djbX19fX2cvLix8TE2MDPH2GPn/+fCtPT08XX19f5+zsbG5b7UpOTtb38fHhu7m5DdqyZYsJADx69IgzYsQIR19fX2flPhYuXGgZERFhFxQU5HT27FntSZMm2QPAjBkzrH18fPheXl78jIyMZmcSUxTV/fWa4aPyg/M7JnTulS0dFjqXmpqq7+bmNqi4uJjz22+/iQBgxYoVVt98881fdnZ2Eh8fH/6bb75Z+tJLL5WvXr06z83NrW7hwoWWrq6utQcOHLgfEBDg1Di51NzcXJqamnqXy+Vi0qRJ9jdu3OANGTKkTnmGfuHCBd1Lly7dYbPZUCXQbcyYMVUREREiiUQCPz8//vz584tXrlxpvnz58rxRo0ZVv/XWW1YnTpzQAQAnJydxcnJytkgkajjgJyYm5urp6ckPHTqkl5iYaPbZZ5/lqvK7UBTVvdAegYqeJ3QuMDCwIj09/XZERETJH3/8oQUA5eXlbD6fX8/j8YiNjU39o0ePninG7u7uYgAwNzeXNO5NFBQUcMaOHTtQIBDwL1++rJuTk/PUWf+iRYvyw8PD7ePi4qyrqqra/H+blpamLRQKnYVCIT8jI0MLAO7evau5bNmyAQKBgH/27Fn9hw8fcgHAx8fnme+8cuVKcy8vL/7KlSut8vLy2uyBUBTVPfWaHkFLZ/Id5XlD5wBg5cqV+ePHjx8YGRlZbmBgIBOJRBp2dnaS+/fv8ywtLaVcLpfIZLKG9RmGaZju3Xjm986dO43DwsLK4uPji8PCwuwJIYxEIkF+fj7H2tpaOmHChIqoqKjypUuXmh84cMBg5syZzSaAKm3YsMF8165d2XZ2dhIHBwc3AHB0dKyLjY1tuCchkUhw/fp1bRaL9dQU9Pz8fHZqaqrelStXRAcPHtTfv3+/saq/B0VR3UuvKQSdLSYmpmz8+PEDAwMDnXR0dOTqhM6ZmprKLC0t61NTU7XXrl2bGxUV5SCTyTBr1qzHPB6PjBkzpmLOnDm2r7zySqsH7tGjR1fMmDHD/vDhw4bKZSKRiPfvf//bPCkp6X5oaKhjbW0tCwDmzp17r612hYWFlYaFhTkOHjy4Rl9fXwYAa9euzZs+fbptZWUlm8ViYefOndnNfdbMzEymo6Mj8/Pzc3Z1dW3309Uoiuo6NGtIDd0xdG737t2GxsbGsrCwsMqubktfRbOGqJ6OFgI19KTQuaNHj+quXr3aSvmaYRhCCGGUr21sbOqSk5Ozu6RxvQwtBFRPRwsBRbUTLQRUT0dHDVEURfVxtBBQFEX1cbQQUBRF9XG0EHQSkUikoYxiUFdiYqJJ43kF6uqqYDh19ttcuF5ycrJ+UlKSQePfbvr06dad0VaKov5GC0E3tH//ftOmERHqFAYbG5suCYZrbr/qtDsiIqIiKiqqvPGyPXv2POig5lEU1QJaCFTUntA5pWXLlpn7+Pjw3d3dXdLS0rQAwMvLi698XyAQ8E+fPq19+/ZtLaFQyN+6datxeHi43bRp02yCg4OdJBIJJk6caO/t7c2fOHGivUQiQUpKit6IESMcg4ODnfz8/JwLCgrYjc+omwuGi4uLs5ZKpTh+/LiOu7u7i5+fn/OmTZva/D4ikUjDz8/Pedy4cQ6DBg1y3bNnj2FAQICTl5cXv6KigtV4vx4eHi7R0dE2s2fPHtBcuB0A7Nu3z0QoFDq/+uqrdoCiJ5SQkGDaeJ/K30fdQD2KolTXa2YWy6/t65jQOY83Oix0rqnly5c//vDDD/PT09N5y5cvtzx8+HBW03VCQkJqBg0aVJuWlibicrk4deqUfkBAQNXevXtzdu7caeTi4iI+cuRI1uLFiy12795tZGFhIRWLxcy5c+cyduzYYbR582az6dOnlyi311ww3M6dOx8AwJEjRwzWrVuXO2HChEpVQuoAoLy8nJOWlnZ3x44dxnv37jVJS0vLWLp0qfnBgwf1/fz8GvKISktLOWvWrMkbOHCgpLKyktU03A4AHB0dxd9//312dHS0zcmTJ3Va26+6gXoURamO9ghU9Dyhc01t377d2Nvbmx8XF2dXUFDw1Fltawc3f3//agDIzMzkeXl51QCAQCCozszM5AGAu7t77ZNlNffu3XvqATitBcMtXLiwMCkpyWjy5Mn2Z86cUenBN05OTrVsNhsDBgyQuLq6igHA0tJSUlJS8tRJhYmJiXTgwIESoPlwOwDw9PRs+D1FIlGLD+4B1A/UoyhKdb2mR9DSmXxHaU/onNKXX37Z75YC780337QDAEIIU1tby1y/fl1TuR6Hw5FLpVJGmWfEYimOew4ODnWXL1/WjoqKKr948aKOo6OjGADS09O1AODSpUvaDg4ODbOZWwqGy8rK4trb20tMTU2l+/fvz8nOzuZOmzbNLiQkJKOt78AwTOP/bjYc70mbGxY0F24HAH/++adWZGRk+bVr17RjY2OLb926pYkWqBuoR1GU6npNIehszxM6d+7cOT2hUOgMALGxsYXDhg2rFggEfH9//yrlOtHR0UU+Pj4uI0eObLhJGhoaWj5mzBjHuLi4wsbbe+ONN8qmTJli7+3tze/fv79k7dq1+cePH9flcrkkKCjIqa6ujnX48OHM8vJyNtByMNzUqVPtU1NT7yYkJJgdPnzYqKamhrVgwYJOu7ncXLgdANy5c0fT39/f2dLSsn7UqFHVrRUCdQP1KIpSHY2YUEN3DJ1LSUnR+/XXX/USExMfqbK+TCZDXFycNR2N03FoxATV09EegRpCQkIaQufWrFljuXz58m4bOtcSNpvd4pDM4OBgp+4cpEdRVOegPQKKaifaI6B6ut4w+kIul8uZtlejqI735ESqZ59NUX1ebygE6YWFhQa0GFBdoaamRpNhmIKubgdFtUePv0cglUr/kZ+f/2V+fr4bekdho3oWmUwm+1dXN4Ki2qPH3yOgKIqi2oeeQVMURfVxtBBQFEX1cT3uHoGpqSmxs7Pr6mZQFEX1KFeuXCkihJg1916PKwR2dna4fPlyVzeDoiiqR2EY5n5L79FLQxRFUX0cLQQURVF9nMqFgGGYVh8cQlEURfVMbRYChmGEDMPcAnD7yWsPhmG2dXrLKIqiqBdClR7BRgBjABQDACHkGoDgzmwURVEdj1Q+hvzntZAfXgFS1+UJ6lQ3otKlIUJI09hiWbMrUhTV7RCZFORyEsiRpYCJGWDlAHJyA4hc2tVNo7oJVQrBA4ZhhAAIwzAaDMP8E08uE1EU1b2R3OsgSXNAim8Dg0MAng7A0QD624KkfdHVzaO6CVUKwRwA8wBYAXgIYOiT1xRFdVOktgzyYx+BnP4UcBYA5o4AGDCW3gDPCOBwAR0tyK9919VNpbqBNieUEUKKAES/gLZQFNVOhMiB9J9Aft8F2HkAQ0Yp3tDuB8ZGCIanB5jyQW58D7DFgKwK8sz/geU4vEvbTXUtVUYN7WEYxrDRayOGYXZ2brMoilIXKcwEOTAP5MYPwLBQoJ89wLDBDPAF4zhaUQQAMAwLzJAIoF4KsDlAZTbkuX90ceuprqRKxIQ7IaRM+YIQUsowzLBObBNFUWog9TUg53cCN38C+AGAqY3iDT1LMAP8wGg8OwWIYRgQr2nAua2AnhHw+DrkGjpgmfFfcOup7kCVewQshmGMlC8YhjFGD8wooqjehhACcvc0yO4YIO8a4BuuKAIsLhhrIRj7Ec0WASUWiw0IZgJFDwEWG8j9HfKSrBf4DajuQpUD+qcAzjEMk/zk9asAPui8JlEU1RZSlgtyaiOQfwsYFAQYD1C8oT9AcSmIq63SdlgaWpAL/gGkbQEsHIGcs5CzWGAZ2nZi66nups0eASFkL4AIAAUAHgOYQgjZ19kNoyjqWURaD3J+F8ieaYC0CvB7VVEE2DwwNkFg7Ia3WgQIIWj6VEKWnhkwbCrwIB1gGCD7N5CyFoMqqV5I1Us8dwCUKtdnGMaGEJLTaa2iKOoZ5P5lkJOfAnWVwNAxgKGF4g1DWzBWAjAczVY/L68oRsnHr0NWkAWdcXOhPeYfYGkpbiCz+vMh548Dbv4I2A8Duf8bQALBGNl39teiuoE2n1nMMMzbAFZC0SOQAWAAEEKIe+c371ne3t6EPo+A6ktIVRHImS2A6BRgMwQY6KO4ps/RVFwGMrBpcxuykkco+fdkSB/eaVjG6BpBZ9xc6IydDZaOYmCg/NJXwMNLgIOXYh3rADDGDp3zxagXimGYK4QQ72bfU6EQZALwJYQUd0bj1EULAdVXELkU+PMQyLkdAJcHuA4H9J88YMrIAYylNxgOr83tSPP/QsnaSZAV3gfHehB0I5ai+uftkNw5DwBgtPShM3Y2dMbPBaNrDPLreqAmD3BQHDMYayEY44Gd9TWpF6S9heA0gJcJId0imIQWAqovIHk3FZeBCjMB26GKM3SGBXC1FUNC9a1U2o4k5yZK/j0Z8rICcB29YLwsGSw9ExBCUH8rFVXfb0B9+hkAAKOpC+3R/4DOuDlgTnwIsAE4ChTvDfADY+LUWV+XegHaWwj+C4AP4CiAOuVyQkhCRzZSVbQQ9D019VKIJXIY62h0dVM6HRFXgqR+Dlw/DOgaA24jgSeXbWDsBMbSEwxbtd+h/u5FlHwYAVJdBg23l2D07tcN9wSeWk90AVXfr0fdnycUCzS0oDMiBjqyHDBG/QEnXwAAY+ULxtS5Q74n9eK1VghUuVmc8+QfjSf/UNQLc7egCvHfXENxVT1GD+6H6UJbOPfX7epmdThCCHDrF5DftgHiCmCgt6InwDCAhg6YAf5g9CxU3l7d9dMo/XgqSF01eN7jYfTOLjAazd9M1uD7wnj5D6jPvIKqHz5G3eWfUP3LDoh19WDi5gSGyAFnf5DcCwDkYExdOuhbU91Fmz2C594ww1gD2AvAHIAcwBeEkM1N1mEAbAYwDkANgFhCSKtz3WmPoO+4kFWCd7+7gaq6p1PPAx1NMF1oA08bQyj+hHo2UpSluAyUe01xD8B9DMB7MgTU1AWM+VAwbK7K2xNfPILSTTMAaT20gl+HwVtbwbBVnwMqyb6Bqh8+hvjCj9DQ04GhiwMYGzfFrGVAcW/CbJBa35Hqeu3qETAMYwZgMYDBABpOKQghI9r4qBTAIkLIHwzD6AG4wjDMcULIrUbrjAXg9OQfXwD/efJvqo/76UY+Vh2+Damc4GXXfpg73AHfXn6Ig388QmpmMVIzi+E+QB+xQlsEO5uC1QMLApFJQH7fA1zar1jgEgxYPTnb1tBT3KTV7afWNmv+9xXK/zMPIHJoh86Gfux6MCz1Hk3OtRsCo4V7IXl4B1U/fIzKu6egz2JA5HIwg4JAHl0GCAHTz1Wt7VLdlyr3CH4FcADAP6GIpJ4OoJAQskStHTHMjwC2EEKON1r2OYD/EUK+efJaBGA4ISSvpe3QHkHvRgjB7nP38dmpvwAAMX7WeGeUY8OBvrSmHkkXH+Lbyw9RXqsYv2Bvqo3pQluMdesPLlu9g15XIYWZIL98oLgZbNAfGDoe4HAAMIDZIDDmHmBY6iW5VB/dioo9ywAAuhFLoPvq8g7pMUkfZUDy3bvgyYoh7+8MlttwAIBcywoc57bOB6nuor03i68QQrwYhrmunDvAMMwZQshLajTADsBvANwIIRWNlqcA+IgQkvrk9UkASwghLR7paSHovWRygg2/3MV3V3LBAFg42gnRvtbNrltTL8Whq3nY/3sO8isUYxj66/MQ7WuNKZ6W0NbonnFYRCYFLn0F8vtuQC4FHASA/VDFmzwDRVS0tql62yQEVd+tQ1XyegCA/vQPoTO+Yx8ZQuQyyL97B0zun5CZOoDlPgIMw4I4MwNcz3Bw7T06dH9Ux2vvzWLJk3/nMQwzHsAjAAPU2LkugO8BvNO4CCjfbuYjz1QmhmFmAZgFADY2bU+eoXqeWokMK364if/dLYIGm4W1k13xsmvLl0W0NTiY6muNV72tcCy9ALvP5eCvomokHM/El2ezEekzAFE+A2DUjUYakaIskGMfAAUixYKhUwAzC0AuAcwGP+kFsNXbplyOit1LUPPL5wDDgsGcLdAOienwtjMsNlivrAdJmgt20T3IRLpg8QXQdHRCVcqHqKyWQzf8XWg4+XT4vqnOp0qPYAKAswCsAXwGQB/AakLI4TY3zjBcACkAjjU33JReGqIAxeWed5Ku40ZuBfQ1OUiIdIenjWHbH2xETgjOZhRjd9p9XHtYDgDQ5LAwaZgl3vCzhqWhVmc0XSVELgOuJIGc+y8gkwB6/cGMXgrC1AOl9xRx0fYj1L6MQ2RSlP9nHmp/+wbgaMDw/3ZCyzesk77Fk31WFIB8PQuoKQEZOhkwNgPDMKhK/RnVZ1Kg4R4CvfAl0Bgk7NR2UOpr16WhduyUAbAHQAkh5J0W1hkPYD4Uo4Z8ASQSQgStbZcWgt7lQUkN3v7mGnJKamFhoInPXveAg1nL0cmquJpThl3n7iM1QzEZns0wXTb0lJTkgBxbB+TdVCxwmwDmpfmAtBok8xjAsMDwJ4Lh6au33XoxSjfHoe5SChieNoze/Ro89xdzvZ7k3QT5Nh6Q1QPD3wJhy8GAoPriaVQdV4QUawwOgm74YmgMDu4VI7t6g+cqBAzDfIZmLtMoEULi29hpIBQ9iRtQDB8FgOUAbJ58fvuTYrEFQCgUw0dntHZ/AKCFoDe5+agC8d9cQ2mNBC7mutgc5QEzvbYjE1SVUVCFPefv41j6Y8ie/J0HOpogNsAWw6wNOvUARYgcuJoMcvZzxQFTxxTM6CVg7P1AiBzk7k+AuBTo5waWhXrPeZLXVqL046moTz8DRscQxsu+g4bzix1sR0SnQI6uBMAA41cA4kIABJKSCpTu+QikRtEr4/J9oT9tHb1k1A08byGY/uQ/AwC4QjFyCFA8j+AKIWRBRzdUFbQQ9A6/3S3C0h/SIZbI4edgjI8j3KDD65wbvI/KavHVhQc4+McjiKWKc5LOHHpKynJBjn2omBcAAK6hYIbHg9FUzOolRSKQ3IsAVweMS5hao4PkVSUoWRcOSeYVsAz7w/hfh8C1Gdyh7VcV+X234nIXVwuYtBYoywSIHDB0QPWVc6g+uhWksgQMTxvGKw5Bw8WvS9pJKXRE1tBoQojkyWsugF8JISEd3lIV0ELQ833/Ry4+/EkEOQEmepjjX+NdXsiwz84eekqIHLj2o2J2sFQMaBuDefldMAMD/15HKga58yMgqwdj95JKyaFKspI8lHwwGdIHt8E2s4Xxe4fAMe+6MDhCCMjPa4E7xwG9fopikH9VUQxM+IDJIFT8dxFqf0sCo6UH4/cOQ8PRq8va29e1txCIAPgTQkqevDYC8DshpEsebkoLQc9FCMF//peFL1OzAQAzg+ww5yX7F34NuTOGnpKKfJBfPwJyrigW8EeCGbEAjJbBU+vJH5wDStS/QSwtyELJvydBVpANjhUfxu/9CLaxpdrt7GhEWgfy3TtAXjpgPggIXQ48THtSDJwAC2+Uf/YPiM/9AEbHECYrj4JrN6Srm90ntbcQzACwCsDpJ4teArCKELKnIxupKloIeiaJTI5/p9zBkev5YDMMlo1zxhRP1RI0O7NNjYeeAoC+JgeRPgMQ7WsNfa22Yx0IIUB6iuJ5AfU1gJYBmJGLwDg/22Em1YUgmb+ofYNY8uA2Sv49CfLSfHAHDoPxsh/A0jdR78t2IlJTCvL1bKAiD3AeAQTPBLLPKIqBsSNg7oWyjdNRd/koWHomMF79M7gDaF7Ri9buUUMMw5jj7+iHC4SQ/A5sn1poIeh5quukeDf5Bn7/qxSaXBbWh7shyEm1SVNyQnA1txz5lXXQ5rKho6H8hwPtJ/+txWW36zp/c0NPXcz1sCfOq9XLRaTyMcjxDUD2BcUCx5fAjFoERtvo2XWJHCTjZ6C2RK0bxPWZl1GyLhykqhQag4Ng9O43YGmrN8LoRSBFf4EkvaUohn6xYIaMA8k6DRCZohj0H4ayT6ai7s8TYBmZw2T1z116Wasv6ohCYAXAFo0moBFCfuuwFqqBFoKepbCyDvHfXIOooApG2lwkvu6BwZaqHcikcoIz94qQU1bb6noMAG0N9pNCwXmqWOhosKGtRrG4mlOG93+8hdwyMWKFNogf6fjMOoQQ4PYxkNObgboqgKcHZuRCxeWgFvbx9w1ibTD8MJVC5OrSz6B0w+sg4irwvMbCaMFuMBpdNx+iLSTrd5BDSwAiBzP2PcDaHeSvU4pi0N8djLEzSj58FfU3fwPLZICiGPSz7epm9xntvTS0HkAkgJv4exgoIYR07syVFtBC0HP8VViN+d/8ifzyOtgYa+Gz1z1gbdzyg9Ubk8jkOJlZiLyKOmiwGfjbGoMAqKmXobpeiup6GarrZaipl6JWKm9ze8pioaPBhjb32WKho8GG5pNice1hOd7cfQWEAF9MGwYv27/P8El1Mcjxj4G/0hQLHIRgRr0LRrflHs5ToprxbAAAIABJREFUN4htg8EYtn3wE186itJNsYCkDpqBr8Fw7n/AcFRPIO0q5GqyokCyuWBe3QzoGIJknwFAwAzwA9GxQMkHr0AiugB2fzuYrP6lW9zr6As64maxOyGkrtUVXxBaCHqGP3LKsODAdVSKpRhipY9Nke4qxz3USWX49W4hiqrrocFmQSwhKK2VYYiFLjyt9GGm+/R2ZHKCGsmzBUL539X10oZho61hGDT0Kh4W1+Dz01kwN+DhwCxf6PLYgOgkyKmNiucF8HTBDI9XDA1to6chf3AeKMkEdC3AOLTca1Cq+e0blG+bC8hl0B4zE/ozPlY7QbSrEEIUv9G1g4CWIZipXwDSSpCHFwAwYOyHg3D0ULI2DJJ7V8G2dILJqp/BNlQvZZVSX3sLwc8AXiWEVHVG49RFC0H3d/zWY7x36BbqZXIMdzbFB1MGQ4urWoZOTb0Ux+4WoqxWAhYDpOfVoFb69N+otaEmPK30McxKD6YqFheZnCh6E5KnC0TjHkbTYpFVUIWUq/mIcNXGUo0UIEPxSEfYChSTw/TaPng9dYPYeQIYTYNW16/+eTsqdi0GAOhO+Sd0I9/rcTNziVwKcnAxcP8SYGIPJuo/IKUZQMENgGGDGfgyCGGjePUESO+ng2MzGCYrU8DS6z43wHuj9haC7wF4ADiJpx9V2erM4s5CC0H39tWFB0j4NQMEwKveVlg8xhlslmoHskqxFL+IClBVL0OdVA7RYzGkcgI/WwMM7q+La3mVuJFX9dQB29ZIWRT0YazdvksnUjlBbb0MD8prcTGnFAQAI/ofJhR8DWOmGuBqKeIhhkxU6eD81A1is8FgWXq2si5B1fcbUPXtBwAAvTc+gO7Et9v1fboSEVcqbh6X3FfcOxm3EuTBeUW2EpsHxikUcnEdSlaNgzRXBI79UJi8fxgsHfUypijVtbcQTG9uOR0+SjUmJwQbj2fiqwsPAADxIwZiutBG5bPZ0pp6/HTnMeplclTXy5BRJIa1gRZe9egPW6O/b5BKZHLcKqjGH7kVuJFXiXrZ33+/9sZaDT0FQxWGfrbmQV4+6k8kwL7wPADgKnGAdeRq9Btgp/I2SPFdxSWRNm4QE0JQuXc5qo9uVSSIzk6E9ohp7Wp/d0DKckH2TFPcG3ljF2DqoBhJVPkI0NAF4xgKeWUZileGQlaQBa6zAMYrDjb7XGWq/bokdK6z0ELQ/dRJZXj/0G0cv/0YHBaDVWGDMG6Iucqfzy6pwel7RQCACrEMj6skmDS4P7yt9Vsd6VMvleNmQRX+yK1Aen4VJI2KwkATRVEYaqUHA031igK5lwZyYgNQXQIpi4fL9lE4zhUir7QO26YOVWn0kao3iIlMivLP41H7v/0AmwvD//svtPwmq9Xe7kz+v0Tgj+8Ax2Cwwj5QPJXt3q+KXpKWMZiBoxUzpleOhazoATRcA2G8LBkMT7VBBZTq2tsjyEIz4XOEEIeOaZ56aCHoXsprJVj47XVczSmHLo+Nj18dAl97Y5U+K5UTHLtTiLzKWrBZDMprpbAz0kGoiyk0VbynoFQnlSM9X1EUbuVXQSJX/MkyAAaaasPrSVHQayPPiFw7CHLySWK6lQfKg/6Jn/M5qJPJUVAmhqWuJqb5tx0L8fcNYnMwDqNa7BlVfL0a1Yc+VSSILtoP3tBRan3v7o5UFYH8N1JREGN2gunnBCKpVdw3qa96MsM6BLKCbBSvHAt5aR40PEbA+N0kMBqabe+AUll7C0HjOziaUITOGRNC3u+4JqqOFoLu41FZLd7+5hqyimrQT4+HxNc9VI55vlVQhaO3H8NYWzFkUyYHJrj2h3kHpI+KJTLcyK/C1dwK3CqohrRRUXA208YwK30MtdSDbpOiQG7+pAiLA8AEzgJ8osEwLFSIJfjxZj6kcoLS6nq87GiGIQNavulLaooU9wbAKGYQt3CDWFqYg8J3/p+9846Tqjr///vc6Ts72/sCCwssHUQQUFHBbqzYEzXR9Pw0MT3RaIzRVKMm3xhj1ERNTKwxdrE3pEgTWKRu7313er3P74+7LLtsG5ZdWJDP6zWvmblz7znnzu6czzlP+TzzIBIi7ecvY5t58gHf+2iE/t6fYcPTXbsCAAm5kV3LIRaC1ImosccTq91Fy21no7ubjbyJHzx+WITMHi4YdtOQUmqFiCwe/Mzhx1EiGB3YUe/h209sotkbZmKmkz9/fg45yYOv4Jq8YZ7b0kCtO0hBqhWlFNmJNs6ZmjUi0TGBSIwtdV7WV7vZ3uhlj/VIU1CU6WRefhKz81wklL5rCKiJjjr5etT8K3u04w/HeHx9FRazRiAc49LZuWS5et+viHQ6iFsGdRC3/d9XCa54GvuJl5J64z+G9b5HE8TXgjx8eY9dgXG8CSl5szPhbBZazjFEKoppuf1cxNuGfdEyUm78O8o0OsuOHm440B1B9/9kDZgPfEtEDkmR0qNEcOixobKd7zyxCX84xvyCFO6+fBauQezwoajO6zuaeWd3K2kJJsamGCv/2blJHJs/srUB9sAfjrGpzsPGGjfbG310bhSY49/IdQ1/Q0MnsvDL2E68rs/rOwIRHlpVTrrLhq4L58/IISux5w5GWnYh1asHdRBHSjfS/NNTwGwl84/rj/gM2752BQDirkbK3sNIOFuISi8iUrqRltvPRwJuHCdfSfL/e+CwyaMYzTjQmsV3d3sdBcqBy4dhXEdxGKK+I8gPn9mCPxzjrBlZ3H7BdKzmAfR4RFhX7eb54kY6glFykyzkJRmx/wvGpjAj5+Dp5iRYTRxfkMLxBSl4Q1E213lo2rqCc8oeQkPnzeSzea1pEVNXVXHWlAwmpPWUc0h2WDh7ShbPfFLL+Cwnr25r4PSiTMYkG+dJNITUbQBA5c4bMErI/a9bAXCe880jngQA1HFfQDY9D7s/QBp3de0KVNIYGLMQqV6NVH8MZgeWwrmk3fxfWu+8iMAHT4LFTvLX/3TY5VMcThiUZkVkabfHGSLyNRHZcTAGdxSjC+Gozo//W0y7P8KiwlTuvGjGgCRQ1R7k3g8qeGxdLR3BKNOyHOQlWVHA4glpB5UE9kWizczxpnLO3/1HzESpnXg+OyZfRUyguN7LfR9VUt0e7HXdzPxkxiU7+LTajQBv7WyipMVQLpX6jUY1ssQcGEBGIrTxdcJbP0A5U0hc9oORusVRBeVMhzlGNJSseqTnZ+mTIXs2IEjFh4ivCeuUhaT+9Gmw2Am8/SjuR3/C4RbheDhhUCJQSqUrpf5PKbVBKbVeKfWnfRzIR/EZwd1v7KK4xk1Oso1fLZvRb6KYJxTliY11/P7dMkpbA7hsJpZOSiHBqqEpWDopg8kZB7d28L6Q2mJDIC0Whlnnk3/Bj/jOyeP51TmTOTbfRSiq89dVVbT5I72uvfaEAlo6QmwoM5LOPihtobS6Elp2AQqVf1z/4nOxKJ7HjTiLxEt+jJbYW6n0SIU67gtgtkHJh0jjzp6fZc+GtIkgMaTsXSTkxjbjJNJ+9B8wW/G/9gCe//ziM0MGIoKEfUhHDVK/FSlfib7tFaR0ZLQ+4zENPQl8AFzS+f4qjLKVR1ac21EMiJc31fHM+hosJsVdl87qM2ErpgsflrXxyrYmAhEdTcGSwlRcNo0adxCzpjhtciZ5SYc2LFAadiDP/RAiAZh2plE/oHPiTrKbuWZeHu5gFbtb/Px1VRXfO7mgh0SGSVPcceF0rnzwY/yhGIunpuNs2mgsqzKmouz9Z8cG3n2caPV2TFnjcZ71tZG+1VEF5UxH5lwE659CVj2KuvDXez9TCsYsQiJB8NQgpW/DpLOxHXM6qd//J213X43vhXtRNgeuS396CO/iwCGiG3pVgXbjEWxHAu0Q7IBAGwQ6INgOsd6LEEktQBUOf3RZPM7i9SIyb59j6/pzOow0jjqLDz521nv40iPrCUV1fn7+VKZmJ9Hui+K0mUh1Wkh1WqhqD/DfLY3UeQwVkmlZTi6ckcmWejf1nhBWk8aZRZlkJg5fcfqhQJpLkae/bfwQJ5+COvcXfdYM9oVj3PN+OQ3eMFOznHzr+LG9dkAvbarjthe3ccNcnWvH1uEXK8XJS1kwPrPPpDM96KXpO3PR2xtI+e4jOE64pNc5RzrE12LkFURDqKv/jsoq6vl5LGJEEgVauhLOlMlCYOVztP/pyyA6rqvvIPGCGw/RHQwMiYX3TuSdE710n+AD7RByG0V7BoPZDo4UsCcbz44UVGI2atyCIY3tQJ3F7yqlrgSe7nx/KfDKkEZyFIcdPMEIP3y2GLOmuPHMInISE2jzGfV+vaEY3lCMqtYgvnCUMS4byTYTJ09MZXJGAm/taqbZH8ZhMXFWUSapCfEJxI0UpK0SefZ7BglMOB71udv6LRzvtJr41glj+cN75Wxv9PHExjquOja3h8nnvNk5rCup56KsTwHYEJtISXOQQKyZkwszehGH76U/o7c3YJk0D/vxF4/cjQ4ACbRDRzXouhG2KdL5rHce2/PofUz2HOtxvFsbPa7v9pnFAfZko2ynPQVmnAXFryIr/4G66Lc9xqdMFpiwFNn9OgRakYoPYMJSHCdcjESCdNz/LTyP34qy2HCe882D973pMQh5jFV70A3BDmOCD3ZO7oHOyT7ij69Bm6vbJJ+697vpnPCxp6AsB2/n3O+OQCnlwcgoVoCTvbUINMArIofE03d0R3DwoIvww6e3YNI0zp2Zi63TPJKVZCXLZWVVRTtNnjC5ibYephNNgT8apS0YQZTOmUXZuOyHNhZcOmqRp74N3kYYNx910W9R5sF3J+WtAf60ooJITDh3WgbnTM3s8XmwfBXWjt2sb7ay2b4Au8NMOCbkumycNjmzq8JZrL2Bpm8fg4R8pN++HOu0E0bkPvuDhP3IjuVQ+gHo0YPad78QAWsiONO7JsU9E6KYbdC6A/QIpBaixp6AUgr/W4/Q8aCxG0j+xp9JOK1PKbT4h6DHuib2PZO8dHvddTzkpQ+Bhd5QJthnUledk/3e48n9LkBGEke1ho5ivyEiPLOuFrvZTGZntm+yw0xBhoPqjiCPb6ilLWBMKPPHJHFmUTrRqNDijRCM9Nz2JtpNpCUYJiSHVTvoYYDiaUSevgE66iB/NuriP6As8Vf62lzr4aE11QjwxXl5LBhnZAqLvwXZ9SqC4op3M6n0WfnL1cdQ0uYjENFJT7ByRlEmDouJjgdvxP/WI9jmn0vaj58YoTvtDYlFoOxDZPvyvavVjMmG2UFpoGnGszIZBRmUacBjSnX7rOu8fc7XTD2Ph33GCrrLDt5u1DeW2MCDtzkhd5rRnrfN2HHYk4nUlRNcu5xYIIDjnBuwL7oE7Ek9wnUlFjFW8IF2Y0IPdSCBzsk95O4033QYMhdxQYEt0VjB25M6n5M7J/mUvat7W6LxHY1CDEepylRgMobEBHC0VOWRDE8wSnG1B2MzaJDCtLxEUhLMfFTeztOb6tEFxiTbuGx2DhMzDIGwVn+YN3Y0EokJOYl28lx2PIFYj3WUzayR6jST6rSQ5DAfUK3heCC+VoME2qogexrq0ntRNud+t/NeSSvPbm7ApOD6E8cxOSPB0MvxN0PmNO4pTuY/a6ooSEvgwS/O5f2yFjyhKEl2M6c7W/HftBiUIvOejzHnTR6BO+0JEYGaDcjWF8HfYhzMLELNXIZKGTvi/Q8G8bUif7/cMCGd+3NjJ9CDLDrt7OiQWWgQTFMZeBr7b9SaCFanQQDxmmhQYHd1Tex7Jnm1z4SPzYXS9k//arThgHwESqmvAjcCY4BPgEXAKuDUQa77B3Ae0CgiM/v4fAnwAlDWeeg5EfnlYOM5ipFDKKJT0RKgxRsBFJ5ghPZAmMvm5yHAc1saebekFYAzitI5f/pep2iDJ8RbuxoNs0iSjSWT0rGYNGK60O6P0OaL0uaPEIrq1HeEqe8IY1KQnGAhzWkhxWkesFD8UCCBDsMn0FYFmZOMncAQSABgycQ0WnwR3i1p5aE11dw8VyfZ3wxmByp7Nt/OMrGmtJWSJh8Pf1jOjWdM4s0djbQGIlQ9cTPpopNwxlcODgk070aK/wdtFcYBVw5q5jLInj5qkrKUMw2ZswzWPwlbX0dd2Knx1Me5euOnULceMifA+BNRsSgS7CC6ew3SWo2WkICWkIgKe/eu8JUGtiTDHGNLMnYMjpSu18YKPsmY4EfpCv5gIh5D1Y3AccBqEVmqlJoK3B7HdY8C9wH/HOCcD0XkvDjaOooRRFQXalqD1HWEEIGorvP2tkbaAyHuvnw2oajOI2tr2drgxaTg83NzWVSwN0SypiPAO7ubierCuBQHSybudZSaNEV6opX0RKtRfCUYo80fodUXIRDWafUZrwFc9r1RSA7LgZmQJOhB/vt9aCmFtALUJfegHAfm1lo2K4tWf4SdDW2YG4tBA5V3LMpkxQb86qLpXPOPdTyzvobFk9M5Z2o2q998mfTd7xG1JhA4+/sMXJ/swCCeBmTrC1C32ThgS0JNOxcKFo3K1ayRbfw/KFmBNOxAZU/p8zwtazqiR5CGzeBvgIlnoDkzscy+DM+/f47vhT+B2UbaDx/DOmmusSsYxSaa0Yh4vqmgiAQBlFI2EdkO9P0X64ZO01HrAY7vKEYQIkJ9R4iNFW5q2w0SqHcHuP2lT1lb0covLphORzDKvR9UsLXBS4LFxA2Lx/UggfJWP2/taiKqC5MynCyd1DtaZg+UUrgcZsalOzhmXBJzC1yMz3CQ7DCjAE8wRmVLkE2VHjZWeihvCtDhj6Dvpx9Lwn6jVGLjTkjOR136R1TCgSduaUrxpfl5fD6nAacWpSqSRDhxbwZxUY6LG5ZOBOD2l7bhCYSZ9v4fAChd8BXeaICq9sABj2NfSMiD/slTyNu/MkjAZEVN/RzqzNtQE04clSQAGH+TOcsAkNWPDnxy9mxIm9SZcPYOEuww/p+u+iUJZ38doiFa77mOcHUJyp50lAT2E/F8W9VKqRTgeeBNpdQLQO0w9X+8UmqTUuo1pdSMYWrzKAaBiNDmi7CpykNZU4BoTHDZTbQFgtzxyja8oSh3XTaT9mCUu94ro8YdIivRyg+XjGdyxl7Tys4mL++VNKMLTM92sXh82n7Z/O0WE7kpNqbnJzK/MJminAQyXRbMmiIU0anrCPFprY91ZR2UNPqJxAaPvZZICHn+p1BXDK5s1GV/RCVmDOl76guWSDtz7IaP5PH6sTyyrrYHUV21aCzzx6fQ6ovw37/eR7T0E7TUXMxnfYuYLry9q4mSZt+wjEWiYWTHcuSN26HsQyMKZ/wJBgFM+1xcUVGHGnuzjY1dQb/nKYUasxBc+RALI2VvIxE/SimSrv09jlO/BOEAbb+9HN8bf0cioX7bOore2K+oIaXUKUAysFxEwnGcPx54uR8fQRKgi4hXKfU54E8i0qcBVSn1deDrAOPGjZtXUVER95iPoid8oRgVzQE6OiN+bGaNggw7je4gX35sg5E0dt5UxmUm8q/1tUR0oSgzga8uGEOCde/KsrjezdqqdgDm5iczJzdp2OzPIoInGKPNF6HNb5iQACwmxcSsBFKd/Yi5RcPIizdD+RpwpqMuvw+VOmZYxrRnXLL7dfA34U+azC+3puGPxDi5MJXLZmd33X99R5CrH1jB3Zu/QVakieRv3Y9jyVWsr+5gS70bODDBPREdKj9GPn2506EKZM9AzbwIlZQ7LPd6MKG//xfDV1B4Ito+eQX7QmIRpPRNwwHeLeFM9Bgdf/kmgQ+fAkBLzyfxwu+RcOoXjxa46cQhCx8diAj6OLccmC8izQOddzRqaGgIR3WqWoM0ug3+NmkwJtVOTooNbzDKVQ+vpaY9yIVzcllYlMnL25oAOGF8ClfMyelh7tlc52Z9tTEBLRyXyvTska0x6w/HKG304wka4YbZyVYK0h09xiSxKPLyz6HkQyN2+/I/o9LHD+s4pLUEqVoJZjtq6oWUtEW576NKorqwbGYWp03eK8G14cE7yX3r91TYC8i9eyUFmcZ31J1AhyLBLQ3bkOLnwV1jHEgeg5q1DJU5qLV21EL8bUa9gmgQddXD/foKus6PBAxCDnsgMRc1YSlKMyG6TnDNC3if/R3RKiPJT0vNJfGi75Fw2pdQ1vhDho9EDEQEh8yQppTKUZ2/AKXUgs6xtByq8RypiOlCdWuQjRVuGt1hFJCTbGVuQRJ5qcZK6ZbnP6WmPci0XBfjcly8vK0JBSybmcXnj+lJAiUtvi4SOGlC2oiTABjy0TPyExmXbkcBDR1hNld58AaNXY3oMWT5nQYJ2BKNENHhJoFYGKlbD+yRmLYyKSOBq+cZK/D/FTeyscZY7eueFvJW/g2Ah3Ou5ZYXt3eZtWbmJHHShDQUBqGuLG+NywciHTXoH92HrPyLQQKOVNS8L6KW/viwJgHo9BUc0+kr2EeZtM/zLQ5U4WmGSclbZ0hYi6A0Dcfxy8i4ayUp3/8X5oKZ6G11uB/5MY3fnoPvlb8goXjDSj9bGLH0NqXUE8ASIEMpVQ3cBlgAROQBDKmKbymlokAAuFIOt+y2UQwRodkTobI1QDhqfK2pTjMF6Q4c3Uw8D39YzordLaQ5LcwqTGd9jQerSXHdcfnMyu05yde5g6woM7h6wdgUJh1EBVGlFPmpdlISzOxq8BMI6xRXexmTaiVv7R9hx9tgTUBdfDcqc9Kw9y/1n0A0BM4sSJ3QdXz+mGRafRFe/LSJx9bVkmw3k/HSXYi/A9PMJdSmLKK+1sPDH5bzrSVGme9JGYnYzCbe3d3MzmYfwajOKRPTMfdRfEUC7YYJqHINIMZuZMpZMPEUlOnQSnYMJ9T8zyOf/A9KPxowgqjrfJsLJpyGlLwBbaWIJQGVO9f4TNNwLLoQ+4LzCa17Fc9/f0e0bBPux27C+/y9OC+4kYQzvoxmH1oo8ZGIeBPKCoDJIvKWUsoBmEXEM+Kj6wNHTUODwx2IUt4cwBcyTClOm4mCdDvJCT1t6x/tbuE7T2zCYTNz4vRsfBGdFIeZbx4/ljH7lJ1sC4R5dVsD4ZgwPdvFwnGHTj5Z14XKliB17UHGb3mA7PJXEbMd7ZK7Ufmzh70/CbQiO18FQBWdi3L0vHcR4clP6vmovJ28YB3XvXgZ6FEyfreCzWoMX//nRpSCv187jzndah03eIK8tauJcExIdVhYMjGjS9VVIkFk15uw+x1DhVJpUHgyasrZKNuhlfAeKegf3A/rnoDCE9Au+l1c14i7Bil7FxBU/gJURm8CERFC65fj/e9viZRsBEBLyjAI4cyvoNmPzO9zXxyQaUgp9TXgWeBvnYfGYEQQHcUoQyAcY0edj601XnyhWJdzddaYxF4kUNse4Jbnt5KcaOXYokx8EZ2CVDs/WjK+Fwn4wzHe3GlMWAWpDo4b27/M8sGApikKMuzMrfwX2eWvomsWdh73MxoTpwy7Xr2IGJWzEMiY0osEwNitXD4nh+nZThauvw9iEcwnXoll/CzmFaTypRPGoQvc+vxWfKG9Oj/ZLjufm5pNks1MWyDCS5/Ws6uxAyn9EHnzdtjxukECeXNRp9+KNvvSI5YEwNgVYLZD6Uqkfnt81yTlo8YsAkBqPkY6KnufoxT2+eeQ/uv3SP3pM1gmzUN3N+N5/Faarp+F9/l70AOHZF07ahCPj+B64ETADSAiu4CskRzUUewfdBHKmwNsqvTQ6ougKRiTZmduQRJZSdZezshQNMaPni3G4bAyvSANXWBunosbFxeQvE/t4UhM581djfjCMTKdVk4uTB9xWYh4IKv+gXXTM4hmpv6kW2jPPIbSpgA76n1EonFI/MaLtlLwNxkmmZz+y3SbNMUXU+uYXvkWEZONZydeR7jTL/CtJYVMyU6kui3IH97Y1eO61AQrF8zIoTDVQa5vJxlr7kY2PWXIJKRNQJ38fbSFX0ElZvbV7REFw1dgqLLK6sF9BV3XpU/q+ttIxQqkvTcZQCchHHsW6b96h7Sbn8NStADd04LnP7+g8fqZeJ67C93vPvAbOQwRDxGEuoeKKqXMxCXDdxQHAyJCSYOfuvYQAmS6DEfw2DR7v4ldv3ttJwFRFOYlg1KcNSWd6xbk9yo7qYvwbkkzrf4ILpuZ0ydn9mnHPtiQjx+H1Y+C0tA+dxtj5p3K5OwETBq0+aJ8UuXpylY+oH5i4W41iI8d0CYvIgSfMCqPbZ7xBbaGk/nXeiPHwGLS+NWyGdjMGi98Usfb23rq5ZjdVZxU+29ObfkfydFW3OYU1uRcTNv8G1DphQd8H4cT1Pwr93tXAEDWLEifbCScVbyPXrXKEJ7rqw+lsB1zOul3vEnaLc9jmbII8bbhffIOgxCe/S26r32Y7ujwQDy/6veVUjcDDqXUGcAzwEsjO6yjiAciQkVzkGavsQuYmZ/IpOyEAesIP7u+hh1tQfIynGgKrpmXy/nTs3qt8kWEVRWt1HQEsZmNojJ2y6HPUJUNzyIr/gYo1Fk3o4qWAJDhsjJnbBJJDjPRmLCjzkdJo5+YPvQ1i9RvgmgQEjIhdeAJObT2FSLbV6G50pn7lZ9hN2tsrPHwQrEx6RdmOrnxdMOJfecr22nqLOAj25cj790FLbvB6iQ4dRnvjf8m282TeGVbA9saPZ+Z8oxwALsCpVD5C1F58w1/SutuQxnW338golIK2+xTSf/l66T9/CWs005EfO14n/41jdfPwvP0r9G9bQd8T4cD4qlQpgFfAc7E0IR6HXj4UEX4HHUW70VNW5DKliAKmJrnJCWh70SrPVhd1sqDq6tJsFuwaIrrTxzHpE7l0H2xubaD9TUdmJTi7KlZZB3iymIAsvlF5K27AFCn/wg1+4Le54hQ1xGisiWICNgsGpOzEnA59i9ATgJtyE6j/pIq+hzKkdb/udEITT9YSKxuN0lfvgvn2d9ge6OP+1dWogtcPiebkwvTEBG+/cQmVpa0sqgwjfuOa4KtLwAKJp+GKjoTZU0S918eAAAgAElEQVQgquusqWxnZ5MhoFaQ6uDE8enYBiD4Iwk98gq+8CAqZ9r+XR9oQyo/NFRMlYbKOQYy4xPcC326Au8zvyW81RBXVo4knOd8A+d516Ml9v8/cDjgQPMILgT+KSKXicilIvLQ0TDPQ49GtzHZAUzKThiUBDbXunlsfR0Jdgsm4ObTCvslgZIWH+trOgA4uTB9dJDAp68jbxm6PWrJjX2SABirvLwUO7PHuEiwaoQiOsU1XipbAnFrFokIUtPpIE4vGpAEAPxvP0qsbjem3IkknP5lAKZmOfnCXCPH4JlNDWyp86CU4rbzp5HisDC2fW0nCYCadzXazItQVuPvYdY0ThyfximF6Vg0RUVbgBe31tHk/WzIJvTYFcSRV9DrekcqavLnIH2KUVmtbgNS+hYShzS1bfpi0m972SgeNGspEnDjfe4uGv/fTNz/uR3dfWSmOsVDBBcAO5VS/1JKndvpIziKQ4hWX4SSRkO8bHyGgwzXwPHkH1d28ODqakwmjUgkyi/PnkRmYt/X7JsrMD6tb7I4mJCd7yGv/xoQ1OJvoo69dNBrEmwmZo11kZdikFhNW4jiai+B8CDFUADay8DXCGYbKveYAU/V/W68zxgSyq4v3I4y7yXkRQUpnDM1AwEeWVtDZVuATJeN+xYH+EneJgCaJy5DjVvYZ9uF6U4umJFDeoIVbzjGK9sb2FLn/kyYitT8K40Sl2WrkPpt+3+9ZkYbswA1YWln4lk9suOlPqOK+oJ12gmk3/oC6Xe8iW3OaUjQi+/5u2m8fibux39OzD2gAMJhh0GJQESuAyZh+Aa+AJQopR4e6YEdRd9wB6LsrDdEy/JTbeSm9L9aFxFe3dbEP9fXglK0dgS49YxJJDv63j20BcK8s7upS0RuqFo4wwnxtRgkIDosuha14Kq4r9WUoiDDwYz8RGxmhS8UY1OVx3Cs9zOZSiyM1O5xEM8bNGnL9+If0d3NWKYswr7g/F6ff25qBgvGJhOOCQ+sqsJTsoapdYaL7e662dy42kV4gCinJLuFc6dlMz3bhQisq27nzV1NBCNxENphjAPdFXS1kzQGVXQ+uPIMsbrygR3J+8I6ZSFpP/sf6Xe+hW3umUjIh+/FP9J0/Uzc//wZ0foSdF87Ehsl5T+HiLi1hpRSFuBs4DrgJBE5JPFsn2UfgT8Uo7jGS0wXspKsFGY6+rV7RmI6j2+oY321sYIsr3Nzy9mTOX5iep/n+8NRXt7WgC8coyDVqCkwGsJE9Td/D1tegsITUBf+dsjCdtGYUN7sp8ljTADJCWYmZfV2rOuVK6GtBBIyUZPOGrC/WGstjd+ZC+EA6Xe+ibWo75V9VBfu/6gSe3MxX2Y5GkKk6FwuezeZ6rYgXzp+XJcjeSBUtvtZUdpKKKaTYDFxSmE6OUlHrqCa+NuQv18BkQDq839D5U4felsi0LzdiAIT3ajVMG4xKqHv30N/CO9ej/fZ3xHasLz3hxYbmj0RZU9E2Z2dD+O1tue1I7HruLbPOcbD1fO4eWCT7/7ggETnlFJnA1cCS4H3gKeAN0TkkFDgZ5UIQhGd4mqPkYXqtDAlJ6HfScodjPLQmmrKWgPEdJ2dle18fn4+Xz1pfJ/nR2I6r25voNUfIdNp5eypWaMjTLS5FPnXdYBCfekxVFrBoNcMhhZvmNLGAFFdMGuKwiwH6Z1mMnFXG1mqymRkENsHLiPTfv//I/De49gXLSP1+48NeG6wegumtQ9hQmetYxHzzryKrbVuvvLoekTgj1fO5qTJg8tle0NR3i9todEbQgFz8pKZk5c0Kkh7JKB/+ACs/TdMWIS27K4Dbu9AHMndESn9BO9zfyC8ay0S9CFBr0Ewww2ztZNIDLIwj59N6neGZpA5UCJ4EngSeE1EDrm36rNIBJGYoasTjOi47Cam5SX2myNQ6w7ywKpqWv0RYjGdLaUtzB+Xwj1XzOpzstBFeGtXEzUdQVw2M+dNyx4VYaIA+n9/ABUfw5yL0U773rC1G47qlDT6afcba5lMl4XxKSa03S9DNIDKm4fKHHj1GakopvnHJ4JmJvPetZhz+g8vlaYdyMoHQI/wkekYnoqdwIKxKVwzL5eHPizngffLsJo07rliFif0s2PrDl2EjTUdbK4zkp9yXDZOKUwnwXrkue8k0G5EEA3DrqCrTT1q7AyaO+sfJGajxp6Isg5de0hEIBJE7yQFCfqQgBcJense6/a663jAi4R83c7tfAQ8oPc0AVqmLCLjjjeGNMYDqlksIlcOqdejGBbEdGF7rY9gRCfBqjE119kvCWyt9/LI2hqCUR0zsHZXE9kuG3dcNK1PEhARVpWPvlwBAClfY5CA1Yk6/tphbdtqNr7Hho4wFS0BmjwRUlvWkBoNGKJyGVMHbcPz+K0gQsJZXx2YBFpKkVV/Az0CExZTMO4CrCsq+biqg3Snha+dNJ5WX5in19Xw/ae2cO8Vs/o13+2BphTzxqSQ67LzQWkz9Z4QL2yt56TCdMYkH1lSy8qRghxzMaz9N7L6EdQw7AqUZkblL0Bc+YasuLcB2fkyjFmEShnarlMpBVYHJqsDkoanEJKIQDRsEETAIA9GqNpcv/t/pdSKzmePUsrd7eFRSn0287APMnQRdtb78IZi2Mwa0/ISMfdT4P39klYeWFVFMKqTlWBmxdY6TMAfLp2Jy963nXFznZudzT5MSnH65EyS+jnvYEP0GPL+/QCohdcMS5nJfaGUIifFxuyxLnJVPamRGmKYqE08FumzhPpehDa/Q2jT26iEZFyX/Ljf86StEll5P8TCMHYBas7ljEtL4Lrj8lHAa9ubWV3ZwU/OLuKyefmEYzrff3oLq0vjq/Cal2zngpm55CXZCUZ13tzZxNqqNvQDSKIbjdgbQbQaqds6fO0m5aOKztvrSK74YL8cySMNpRTKYkNzpWPOKsAybjqWMSMjOd4vEYjI4s5nl4gkdXu4ROTQh5Mc4RCRLvOFWVNMy3P2mzG8rqqDZzY3IMC8fBevrK9GF7j53CkU5fRdL6Ck2ceGUZYr0IXiV4yi80k5MHfwUNEDgV2FGRcwir1XWqdT5bWwpcrTpdy6L0SP4f7XLQAkLvsBmqvv1bt01CAf3WdkJufPRR17VVcd3Vm5Li6bkw3AExvr2NHk5yfnFHHpvHxCUZ3vPbWZNWXxkUGCxcSZRZlGgRuguN7Dq9sb8IQO7yiW7lCOlG4RRI8Ob9sWB2rCqai84/ZmJO98ZcCM5CMR8aiP/iueY0cxfBARKlqCNHsM6Yhpec4eNQS6oy0Q4alN9QB8bmo6z6yuIBITLpufz3mz+y5bWOcOsqJ8dOUK7IGE/cjKvwOgFn9jROvuGsqiayAWgsQcssbPwGbR8Id1tlR5qGkL9gozDXzwJNGKYkwZY3Ge882+2/U0GCQQ8UPOTNT8L/UqIH9yYRqnTTIE/x5eU82nDV5+ek4RFx+bZ5DBk5v5OE4yUEoxJy+Zc6ZmkWAx0eQL8+LWOirajpwiLF27gvLh3RVA58o7c6qRhGZPgbAH2fUa0lBslAX9DCCe0JAeReU7E8rmjcxwjgKgtj1EXbsRFTIl10mivW9Xji7C4+vrCER0pmc7eXlDDfUdIWbmJ/GDM/os/zwqcwW6Q9b+B/ytkDMdppw2sp21l4G7CjQLauwJuBwW5ox1kZVkRYDKliBba/YmoUnIj+fJOwBIvPLWPmvhiq8ZWfFnQz00aypqwVdQWt9/vwtnZnFsvotgVOeBVdU8tLqab5wygWVzcwlGdb775GbWlcevdZPtsnPhzBzGpjgIx4R3djezuqKV6BFgKlKOFJh7CTD8u4K9fXRmJGdMBQSp34iUvIWEfSPS32jCQD6Cm5RSHmB2d/8A0AC8cNBG+BnD/khHfFDaxo4mH06riWggwurSNlISLPz+kpl9mpH84eioqiuwL8TTaBQxB9SSG4acMxBXXxE/UrPW6CtvflfEiEkzajhMzXViMSk8wRibqzzUd4TwvvIX9NZazBOOwbH48t5t+tsMEgi2Q/ok1MKvo0z9//00pfjS/HwunpWF3ayxpd7Lr98u49iJGVwwxyCD7zy5ifUV8ZOB3WzitEkZLBibgqZgW6OXV7bV0xEcHXbvA4Gad8WI7Qq6+tBMaPnHoSacaqig+gxHsrRXjEh/owUD+Qh+IyIu4K59/APpInLTQRzjZwZt+yEdUe8JdSlbzsxM4NGVFWgKfnPxDHKSe69UjboCTaOurkB3yEcPGeUgJy9B5c0auX5EkKpVhhPXlQ9pE3udk+q0MGeci/REC7pARVkVnv/dC0DSNXeg9smzkKAbWfF/4G+B1PGo47+JMg9eStKkKU6dlM6tZxRy3NgkIrrw2vZmxGbhrNm5BCM6335i/8hAKcWMnCTOnZaNy2am1R/hxa31lLQc3ivbnruCoWcbx9VXlyM5f1Q6kocb8UhM3KSUSlVKLVBKnbzncTAG91mCZz+kI2K68Ni6WiK6MCUjgb++sxuAG06dyMIJvQXSRmtdge6Qhh3w6XIwWVAn9W17Hza07gZPLZisqLGL+t15WEwaRTlOinISSHv3T6iQl8CUU3EXHN/DdyAhr7ET8DVB8hjUCf8PZdm/jN9ku4Uvzc/nuycVkJdko8UfwaPDKTNzQCm+88RmNlTun0Z+htPGBdNzmJCWQFQXPihtYUVZC5HY4Wv3VvP2+ArWILUjsyvo6sviQE1Yisrf15F8ZOkMQXzO4q8CH2DIT9/e+fyLkR3WZwv+UIxtdT50gawkK2PTBp5Elu9opqo9SLLdzPKNNQQjOufNzuFLx4/rde5ozhXYAxFB3v+L8eaYS1Ap+SPXV9iL1BoJiSp/AcoyuKM82V1B4sf/RpRG69k3sbvBz856P5GYbji3P7oPPHXgykGdeH2XiuhQMCkjgZ8sncCls7OxmzUiKI4tyiQj1cF3ntjExv0kA6tZ45TCdE4Yn4ZJKXY1+3jp03rqPcEhj/FQQjmSuyLJ9qdewZD7UwqVsa8jefkR50iOZ1l4I3AcUCEiS4G5QNOIjuozhFBEZ1utoR+U6jQPqB8EUN4a4PUdzSigtslDXUeQOWOSueXcqX1eN1pzBXqgZAVUbwR7EmrhF0esG8MktBL0KCSPg5TxcV3n/s8vIBYlYenVjJt1DCZlKMBuKWsmsuIv0FENzkzU4m+jbH2H6+4PTJpiycQ0fn7GRBaOSwYUY7NcTB2fxi0vbeOT/SQDpRRTMhM5b3o2yXYzHcEor21vZEVZC6Ho4Sdet9dXMPK7gq4+j3BHcjxEEBSRIIBSyiYi24GRyWr4jCESM0ggHBNcdhOTs50DkkAoqvPYulp0AQvCpop2cpJt3H35rD6dw6M6V6ATEosiH/4VALXoOpT9wCfSftGyA7wNhrz0mIVxOaPD21cR+vgllC0B1xW3kJVkY864JJKtOoWl/8LcUUHElop+/PWDahPtL5LsZq6Zl8f3Ti4gP9mG3WpmQl4K93xQzvu79z/OPS3BygUzcjkmLwlNwa5mH89tqWN3s++wkrbusStY9Y+D129fjuQdL6LXrDvsCSEeIqhWSqUAzwNvKqVeAGpHdlhHPvZIRwQiOo5BpCP24PniRpp8Yewmxftb63FYTPzxijmkOXs7JUdzrkAPbH4e2qogZQzMuXDEupGQe2/94fxFKPPgNnwR6Uoec57/bUypOQBYtRhTq54gyVdG2Oxia8G1bGqx0O4fGUfixPQEfrxkApfMygIRkpw2ntrUwMOrqwgNIGHdF8yaYm5+ChfOyCXHZSMY1fmwrIXXdzQeVpFFat4VYE2Aio+R2uKD2/ceR3LSWGN32bwN2fY/9MoVSODwLG0Zj7N4mYi0i8gvgFuBvwMXjfTAjmT0lI5QTB9AOmIPPm3w8mFZGwpYs6MRBH61bDpF2Ym9zh3tuQJ7IEFPV0y4OvlbA4ZaHlA/oiOVKw0Br5QJqJTevpS+EFz9PJFda9GSs3BecKPRlh5D1j4CjdvAmggn3IAlKYtwVNhW66O06cDqJPcHk6ZYOimdO86ZhFl0NE3xSZ2X217fzcaa/S9Wk+KwcPaULBZPSMNm1qjzhHihuI5PajpGZPzDjZ67gpH3FfTq3+JAm7AENfncvSbGtjJk58vopW8jnrrDapc1UB5B2r4PYAuwAug9+/S+/h9KqUalVJ90rQz8n1Jqt1Jqs1Lq2CHfxWGE3tIRiQMWmwfwhWM8vqEOgJomL75glBtOLWTJlN4lIUZ7rkB3yJp/GnLA+XNg4kkj11HTNvA3gdlhRIDEM7ZoGM9/fgGA64qfodkTDUJZ9xjUbQZLAmrxDdjT85k5JpGxaXYU0NARZnOVB09gZCQeUh1W/nDRNOyxGN5ABG84xt8/ruG+j6qo9+yfOLBSiskZiVw8M5dJGU5iAhtrO3hhax117tHvTFbHXn7IdgVdY0hIQys4CTXtIsN/oJnAU2uUxtz1KtJWdlg4lQeagdYD6zqf1+/zPh4d6EcxCtn0h3OAyZ2PrwN/jaPNwx7xSkfsgYjw1Cd1uINRgqEI5Q0ezp2Vw7Un9FZJPBxyBfZAOmrhk/8CoE4ZueQxCbYj9Z8Y/Yw9Pm7JCv8bDxNrKMOcPwXH0msMEtjwb6jZAGa7ER2UPMZoVynGpNmZNTYRh1Uj2FknuaI5/jrJ+wOzpvGbZdPJtCpKajqIxnR2NPn49dulPF/csN/Vy+wWEydNSOfsKVldzuTlOxr5sKxlVFdCO9S7gh5jsSYa/oNpl6By5hg+hEArUrkC2f480rR9VOcgDJRQNkFECjufJ+zzvn/d3b3XfwAMJJZyIfBPMbAaSFFK9S2Oc4Sgti0Yl3REd6yrdrOhxoOIsLWijVn5Sdxy3pReE+fhkCvQHfLh3yAWgWlnonIGl30eUh97TEKiQ9okVFJ8Yam6rx3Ps78DwHX1L0EzIZuegco1Ru7B8d9CpfYmYqfNzOxudZJr20MDCtgdCCwmjV8vm8GMbCcbdjbS2hFAF3hrVyt3vFXKuuqO/TZN5CbZuXBGLnPzkjEp2N3s47niOnY1e0etmWM07Ap6jMdsQ2XPRk27GDVmIVhdEPYhtWuRbc+h132CRAKHepi9EE8egVJKXa2UurXz/Til1IJh6DsfqOr2vrrz2BGJJneYik7piImDSEfsQVsgwtOdgnIltR2k2s3cc/lsbOaeuwgRYVXF6M4V6A6pLYad7xiT6uKvj1xHjcUQaAGLE5UXvzyW93/3IN42rDNOwjr3LKT4f1D2IWhm1KJvoDJ6ZyLvQY86yd0E7KpbewvYHSgsJo3fXDyDkyZlsL2qnZKqNrKcFjqCUR5dW8ufV1RS594/c5FJUxyTn8yFM3PJddkIRXVWlLWyfEcj7YHRt6I9VBFEg0FpJlR6EWrqBajxp0BChpHJ3rjFcCxXr0ZCo0fNP54l4/3A8RiF6wE8wF+Goe++bAF9/lKUUl9XSq1TSq1rajr8UhjafBF2NxpKkOMzHGQOIB2xB4agXC2BiE6rO4jbG+LeK2Z3lVXcAxHh46p2djaN8lyBThjJY/cZb+ZdgXJlj0w/gVakYQvQaRIapAj9HkRrd+F7zbBSJl1zJ2x/FXa/Y5SvXPg1VFZ8kdNJDjNzxrrITjYE7KpagxRX7xWwGy5YTBq/vWQGS4oyaOgI8s6mWpZOSMFpNbGz2c9v3inluS0NBPbTxJNst3DWlCxOnpCOzax1Fr+pY2NN+6gTsVPH7okgWovUbjnUw+kBpTRU8jjUpLNRE8+CpDEgMWjZhWx/Ab38fcR36Oe0eIhgoYhcDwQBRKQNiO9XNTCqgbHd3o+hn7BUEXlQROaLyPzMzN4O0tGM7tIReSkDS0d0hyEo5ycSjVFS28GdF81gyj61BUSENZVtfNrgQVOwZNLozBXogZ3vQt1WSEhFLbh6RLoQPbbXJJQ+BeWKz+KoB7203X01REI4TvkC5kgj7Fhu1LY97lpUzozBG+kGk6YozEwwakmYFN6QIWBX1x4a1t2BxaTxu0tncvLkDNzBKPe9tYurZmezeEIKIvDObsNctLZq/8xFSikmZji5ZFYuRRlOdIFPat28UFxH7ShyJitHEsy9DDj0voL+oJRCJWahTViKmnIBpE0yZCs6KpHdy9F3v450VB0yE1w8RBBRSpnoXK0rpTKB4XCDvwh8sdP0tAjoEJG6YWh31MAf3isdkemyMi49Pv2ZeneI/3UKypXUdPD1k8Zz6tSeBCgirK5sY1ujF03BqZMyGZcySnMFOiHRMLLiAQDUCV89ICmGAftp2AzBNrC6ULlz47tGhI4HbiBatQ1z/hSSTr0I+fQlQKHmXYPKj6+dvpCSYAjYZbgMAbvy5gCf1voIRYYvmsRi0vj9pTM5aXI6HYEo331qM8fmuvjRkvGMT7XjDkZ5bF0tf/qwkqr2/ZvEbWYTJ05I55ypWaTYLbhDUV7f0cgHpc2jxpm811ewFildeaiHMyCUPRlt7PGoacsgawZoFvA1IuXvITteQlp2I/rB/V7jKV5/FXAFcCzwGHApcIuIPDPIdU8AS4AMDOnq2wALgIg8oAxv530YkUV+4DoRGTQa6XApXh+K6BTXeAhHhdQEM1NyB84a3oOYLvzunTJqPSEa2/xMTrVz50XTe1xr+ATa2NHkxaRg6aRMxqaM/lq1su4J5IP7IX0C6pp/9KvTf0B9+JuRXcsBQU08C5WYFdd1vlfvx/3oT1H2RDK+dx9a+TsAqLlfQI0/YdjG1+INU9oYIKoLJm2vqXC4oqbCUZ0fPrOFFbtbSEmw8Ler5zIxy8mayg5eKG7E22mayk60MjvXxew8FwWp9rijy2K6UFzvZlOtm5gIVpPGcWNTmJwR3//3SEJWPWL4CTQz6qybUNPOPKTjiRcSi0DrLqRpm1HMCIxQ58ypkF4Ut1lzMAxUvH5QIuhsYCpwGoZd/20R2TYsIxsCDgciaPNFKG3yE44a0hHT8hIHzRregxe3NvLGzhZC4SghX4iHvji3h3NYRFhZ3tqlH3Ta5AzyD4OC5RJoR/7xeQh5Ucv+gJqwcPj70KPIzlcg5IbM6WhxOojD21fRcvu5EIuS9s3fYOnYDghq9qWoiUuGfZzhqE5pk582n5FrkOo0Mz7DMWwO/lA0xg+fKeajTjJ48Jq5TMpKxB+O8dr2ZtZUduDvtpJPspuZnZvI7FwXRZlOzHH8r7qDEVZVtHWZiLITbZwwPo0Ux6HzT4nohnjhhqcBjECE464+5AQVL0R0aC9HGrcaNS3A2C2kTzaE7zprZgwVQyYCZRRZ3SwiMw9oBMOI0UwEkahOWXOAFq8RXeG0mZie5xw0a3gPylr83P1BOSJQ1+jhwauPIaObzV/vJIFd3UTk8vqoPTAaob9zL3zyHBQch3bJPSPTR+16aPoUbEmoonPj2nHE2upp/slJ6O0NpFz8DWw2I7RPzbgIVXT6iIwTDEJv9kQoa/azRxU6I9FCfqqdBNuBE0IoGuMHT29hZUkrqZ1kMDHLyAON6UJJi5/NdR421Xpo65b8ZjdrzMgxSGF6thPHAOQkIpS1+vm4so1AVEdTMDMniTl5SYc0dFk2PI28dx8gMPtC1KnfHZHd50hBRMBThzRtBW9951EFqRNQmdNRjtQhtXtAOwKl1L+Bm0Skcki9DzNGIxGICE2eMOXNQWK6oCkYm2YnN8UW92okFNW5+dWdhGJCQ6uPO8+d0sM5rIuwoqyVkhYfJq2TBJIODxKQ1krkn18EEdTV/0Bl9h9+OeQ+vA1IyRuAQk0+G5WQMfg10Qitd5xPeNtKkpachyPXqOUwUjuBvhCK6FS1Bmn2hLtC5tKcFvJTbXHlmQzYdjTG957awurSVtKcFh685lgKM3uuKkWE6o4Qm+s8bK71UNMt3NSkYEqWk9m5LmblJpLcTzRaKKqzvrqdHU1eAFw2M8cXpB7SnarsfA957Q4jZLPwRNS5t6Eso3/nvC/E34I0fQrtFYAYJqPpF2Os0fcPB0oE72DIUH8MdEnsicgF+z2SYcBoI4JAOEZpUwB356oq2WGmMGv/t/n3vldGSVsQXzDC52dnc+aMvWGVuggflrZQ2urHrCnOmJxJzmFCAgD6CzcZUtMzz0M78yfD3r7EIoZJKOyBrJlocTqI3Y/dhO+Vv5C4YDHOiYWA6vQJHD/sYxwMoYhOTXuQRneYPT/J5AQzY1LtJDmGTgjBSIzvP72Z1aVtpDut/O2aub3IoDuafeFOUvBS0uLvEc89PtXBnDxjt5Dt6h2d1uAJsbKitSvfoDAtgQXjUgfcVYwkpGYz8sJNEHRDzjTURb9DJQxtNd2jXRHCUekKyXXZzXGbfofcZ8iDNG9DWV2ozGlDauNAieCUPgcm8v6QRnOAGC1EoItQ2xaiui2IiKHqOD7TQUaiZb9tkm9ub+KFbc3oujA1zc53lu5N3NbFqCxVtocEijLJcR0+JCBVG5FnvgMWB+rLT6Cc6cPeh179sSExbU9FTT4HpQ0+8QRWPkf7H68lcd5CnEVTjBDReV9Eje3zd3LQEI7q1LaHaOgIsSdc32U3MSbVTnKCeUj27mAkxvee2syaMoMMHvziXCZkDG5v9oSiFNd72VzrYXujj0i3/IHsRCuz81zMyXUxrpuzWdeF4gYPn9Qa4nVWk+KYvGSmZCUeEnORtFYgz/0I3HWQnIe6+A+o1LGDXydCJCYEIzqBSIxgWCcQ0QlGYgQjOvtOm4k2E0kOc9djpIlhKDhgZ/FowmggAk8wSmmjH3/YMO5muiwUZDiwxOkL6I7qtgB3vl2K2aRhQ+cP3SKEdF14v7SF8jY/Fk1xRlFWnyux0QoRHfn316BxJ+qEr6AWXTv8fXjqkNK3jIl88jkoR+9SnfsiUr2dlpuWkjh7FgmTiowok+OuQ+XNGfbxDRWRmE59e4i6jnCXGqjTZmJMqo1U5/4vNgKRGN99cjNry9tIspu54rgxXHHcmODKAm0AACAASURBVD4lzPtCKKqzvdHLplovxfUe/N1CX/tyNntC0a5sdwCHxcTs3CSKMhPjckYPJ8TXgjz/E2jYAfZkY2eQZ+SERGN7JnidYDjW43VsgKnRYlLYLRoi4O1DQmQ0EsNRIhgmxHShsiVAfUcYAJtZozDLEZdcRF8IRWLc+NynmK1mopEY91w0lQSrYQbQdeG90mYq2gJYTIozi7JGf7LYPpBPlyPLfwWJmajr/rPfdXwHbT8WRna8DBEfKmcOKnv2oNfofjfNNy/FWZCNY3whaBbUoq+hsqcP69iGC1FdaOgIUdseIto5MzmsGvmp9v3efQYiMX70jOFABrCaNM6dnc1VC8cNaC7aF/E6m6dlJdDsC7OxtoPWzloNCRYTsw4yIcR0IejzYl7+CyxVH6ObrFQv+glNmQsHzJI2aQqHRcNu1XBYTNgtGnaLhsNq6jGxx3TBE4jSEYjiDkRHLTEcJYJhQKsvQllnSCgYWcJj0uxD/oOKCDe9sA2vKHRduHHxOKZk743qeK+kmcr2ANZOEsg83EggEkQe+QJ4m1Bn3Yyacc6w96FXrTIK0TvSDQfxIA40EaH97muwWTqwjyvoFJD7JiqzaNjHNtyI6UKjO0xte7Drf9Bm0chPsZGZZI07D0BE2FDZzuOrq/hgZ3OXD+DESelcvWgsC8an7he5dHc2b6r1UNuHs3lalpNUh5nyNj9tgZEjhHBUxxeKEQjHOk06hilnz/eFHmPC5vvJqnwDQVEx6xs0FZ7bNbkbE72pa/I3a2pIprj9JQaXw3xQSPEoERwAwlGd8n1CQidmOXDaDiyi46EV5ayv92E2aSwtTOWSOUb1q5huqIhWdZLAWVOyyHAeXiQAIKsfQ1Y+DFmTUVc9PKQohwHbd1cjZe8aJqGic1H2wesueF+4F1PV+9jyxyAmG9riG1BpE4Z1XCMNXYQmd5ia9lBXZrLVrMhLsZOVZN2vhUlFi59/r6nipU11XZXOpmQnctWisZw1I3tIps4mb6ezuc5DaUugh7PZZTNRlJmAzQShzpjZoRKCiOAP63iCUeMRiPVbrU1hkKajc0WfvvUJnBseNT6c/3nUSd8c9v/P7hgtxDAkIlBKbaEfETgAERl8Hz4COFhEIGKswCpaAsR00BSMS7eTkxx/SGh/eH9HEw99XENKoo2cRAs/O30iSimiuvDu7iaqO4JYTVonCQxPVuHBhPhajOSxSAB16Z9Q44a35pBEQ8iOlyAaQOUei8oaXAMotPltZOVfsWbnIJoV7ZTvoVIGdxqOVogIzd4INW1BAp2+KotJkZtiIzvZtl8TSZs/zLPranhqXTWtPmPBk+my8vnjxnLxsXkkDTFJbI+zeWeTjx2NPtzdJsBku4lxKdauokx2s8acvOR+CSGqC95gFE8giicYwxuM9rLhawoS7WYc3Uw5DquGzaz1+s3+//bOO76u4szfz3t7Uy+2LLnJFQhgTC9h6TGEQEjILgQ2Idl0skv6pu2mbDY/SDaNJQnJJiSBsLApQAiBEDAEAgGbFsAUN8lFtiVfdd3e3t8fcyxLtmRdNduS5vl8js85c86ZOyPfO98zM++8r75yP/rQN0zUumXnml6r5+D81ooRhrDfTdkEC8NYhWCPw/Vrnf1tzv4qIKGqXx13ycbAwRCCZCbP5t0J+lLmP6g85KGxJoTfO/63hk27Y3zq7ldpqC3BLfAfq5ZQGvCQKyiPbIr2u5J+07JaqkJTTwQACg99A17+PSw6A9el/2/i89/6BHQ3Q6gGWXzBiG9z+bbN5O77N7yVlRTUhfu8zyGl0yP0harSGc+yoyvdH/fA7RLqynzMLveP6q0+ncvzwLo2fvnUdprajaV40Ovm0hV1vPPkuTRUjN0OX1Vp68uwPhpnQzTBhnYTr7s84Kau1EvICdAkQH1ZkJVzyigoTsOf6zfMGIjf46Ik4KYk4KEk6Cbkc49uWGvLWvS+f4NMAhpWIJd8HQmUjPzgBDNaYagIj02Yx2s++qSqnj5S2sFiMoWg3yS0M4Vi3rAWVAepGoNJ6FB0xTO899bnmV1bitslvO+kelbUl5IrFFi9sZ2dvUYEVi2rpXKKioC2N6G3vQdEkHfdilQWFx+46Py7t6FbHwOXG1l6MeI/cDzmQryT3N2fxhMOUsjkcF34JVzTRAQGoqr0JHK0dKX6X2BcArPK/Mwp948YDnXfvP66uZNfPr2NNc1d/XmdvayGq0+dx7ENZeMub0GVlu6UEYVonHQ2T0OZnzK/l6DHvZ+pqWAawz2NfknAM6o6DYfu3oje/WmIdxgfWJd9EymdHNfoxZIvKH0pIwo9yRzxVL5/aCboc7Fi3thikB9ICIoZ6A6LyBmq+oST2WnA+JxeHIb0JXNsjib6u9m1pcZb6FjGSYcikyvwyV+/TFlZELdLOKGh1IhAvsDDm6Ls6k0TcESgYoqKAGB8vWgBjn3bxItANonueBrADAmNIAKa6iV//5fwhIPkEwnc531hWooAGDfH5WEv5WEvvUkjCD2JHLu607T2pKkt9VFfHiiqVysinL64itMXV7GhtY/b12zngXVtrH49yurXoxzTUMo/njKPs5bVjMlYIpszY/uah9khH5HZ7v3s8nOFAslcnlgmx/beNNu6kjRUBFlWE2ZpTZiyMVrq7VfX2iVw5Y/Quz8FHc3oHR+Cy75h0g8RbpdQHvL2WyMOFAaPe3ImlYvpERwP3AKUYeYMeoD3qurzk1KiEZjoHkEub0xC23qNSWjA66KxJjhhXzQwb1hf+f3rvLCrj7m1JZQFPHzxvEY8LuHhjVFa+9IEPS5WLZ91SJ12jRdtXmN+UP6IWTwWHHkCt+i8VdGtj0PPNojMQhrPP2AvTRNd5B/6T1yFFLneHjjpA/iOeOOElWcqEEvlaOlK0+WM+wtQXeKjttTHnvbb/AkFZ9cfLWrPn9a5Qmc8wz1/28l9L5r42QrUlfm5fGU9Fx0zm5DPM+CZPXlI/6RubM+kbspY9OxL0Ocyb/sBNxG/mx29KZ7fsdc5XiZXoLUvS3vcfLbPLSyqCrG0JsyymhAN5cV7UB0KTfWh934eWv4GvhDylq8h808cc36HIxNiNSQipc79PRNZuNEykULQEcvQHE2SzSsCzKnw01ARwDXBM/a3PbWN/3lyK0c3ViEiXPfGeSyoCPLQxihtfWmCXjerltVObREo5NDb3gsdzcgbP4yc+M6RHxpN/l3N6LYnjN3/sosRX2T4e2NRCo99C8nEyHZ1km84i9AFkxgS8zAnkc6zoytFe+zwCDW5Z1K3f3w/4B7SMaOqsrUryd929vSbnQrQlczT1JEaZMkS8rpYUh1mSU2IueUB6kv9o3bzorkM+uDXYf1qM/R4/r9OitnzoWJcQ0MiMgv4OjBHVS8UkSOBU1X1pxNczoNGOldgSzTZbyER8btprA0RngCvj/vyl43t3Lh6M8csrkZEOHdxpRGBDVHaYmlCXjerltcO69BryrDufuhohtI6OO7tE5q1ZhPojrUAyJzjDywCfa3oX76HZGJk2qOkclWUnf/+CS3PVCPkd7NkdpiGTJ5d3elBk5F7XgTV+Uf70529k7L3fO99uXyBbF5RFEFwCbhcst+buc8jToNvGv2wv7hJXRFhQWWI+RXBQYJQHnRz+oISykM+2mNZNrYn6UhkeXFXHy/u6ut/virkpb7MT32ZEYb6sgBVYe+wPQfx+OCif0dLauHZO4wo9O2Gk981ZVxZj5Vi5gh+DvwM+IJzvgH4P2DKCYGq0tabYZtjEuoWmFcVZFbZxAUGGcjm3TE+f9crzJ9dQtDvoa7UzwXLqnhw/W6i8Qwhr5sLl9ce1jGGi0EzCbNmAJA3fhDxTNy6BxN28knjRbJkjgnxN9y9PS3oEzdBJkamrZW+rVGqv/Lzaf8jLpagz7zwTDQv7+jhl09tZ/Xru/v9I520sIKrT57LqYsqcY/Tx9BwgpDsTRHyurniuFlUBn1s7kiyuSPBjp4UrX0ZOhJZOhJZXtoV68/L5xbmlAYcgTDiMKfU3+8YT8SFnPkRtGQW+uj3zPe6rw3O/cSUcmU9WoqZI3hGVU8UkRdU9Tgn7W+quuKglHAfxjo0lEjnaYruNQmtCHtYWBPCPwGWB0PRlcjwrp8+SzynHLmgErfAx86czyutvUTjGcI+Nxcum0XJOF0NHw4UnvwfWHMr1B2FXPHDCWt4tZBDtzwGfTvB7TdDQt6hGzLt3IL+9QeQTZDetZOe51+g+uuP4qldMCFlsYzMjq4kdzzTwj0v7CThREJbWB3iLcfWcdLCCpbNKpkQ1wpDDRmFvG6OnVPKkmoTBCpfUNpiaXb0mG1nb4qWnjS9qdyQeQ7sPcwp9dNQFqBq11Nw/1fNS8jCU5A3f2XSwqseDMZrPvpn4O3AQ6q60okvfIOqDumVdLIZqxC8ujNGTyKH1y0srAlSOQbHXcXy/LZubnhgA80dcU5YVovL5eItR1aTyedpj2eI+NysWj6LknGuTj4c0L429JZ3Qj6DXPEDZM7RE5NvPotu+bMJzOH2I4vOHzYgh7ZvQp+6GXIpUi3b6fnr41R85v8IHDc1QhVON/pSWe56fid3rG1hd99elxNlQQ8nzK/gpIUVnLywkrmVwXH9BocThGW1EZZUhwn79v999aVz7OxJs6M35YiE6T0M5XPI5xaOd23nbc3fwZ/tI121BN56A8Gymv3unQqMVwiOB24E3gCsA2qAd6jqixNd0GIYqxAkM3l29aSZVxkoOmLYaNnakeDG1Zt4dH07AMc2VhEO+VhcFaSxyk9HIkvEb3oCkWkgAgCFB74Grz0IS8/GdfHErDHUfBZtfgTiu8ETMCIwjAsJ3f0a+vSPIZ8ltWMHPX95hMjln6XkHZ+bkLJYxk42X+DxDe38dXMHa5q72NmdGnR9dqmfkxZWctJCIw7VY/SnNdykckN5kKU1YRrKgge0KMoXlN2xDDt6HHFwRKLH6T3UZNv4cOv3qM610+6p5s4Fn8Bfs4D6sgCzSnxE/G5KfB4ifjcR/8HxGzQWJiJmsQdYhvn7rlfVQ2Z+cDi4od6XrkSGHz++hd8+t4NcQQl4XVyysoGWeJaQ18XJ80voTeUo8XtYtax22oiAtq1Hb38fuL3Iu3+JlM8Zf575DNr0CCSiJhrTovORwNALmHTXy+jan0IhRzraSffqP+BfcT4V//or5BCGSrQMTUtXkrXNnaxp7uKZLV10JwY3I4tqwo4oVHL8/PJR/05UlZ29KTZEY2ztTvZPcIe8bhZXh1laExlVLzyWzrGjN83OnjTt0TZOe/F66hJNxF1hfjzrWpoDQ89XBb0uIj4jCiV+N5F+kXBT4vf0X4v43UR87glbqzQS4+0RbAa+qao3D0i7T1UvnthiFsfhJATpXJ471rZwyxNbiKXzuAQuObaOI+dV8HhzFyJw8twS0vkCpX4Pq5bXDtldnYqoqgk40/I3OP4KXH937cgPjZRnLo02rYZkB3jDRgT8Qy/515bn0Gd/AVogk4Cu392Ku3YB1df/GVdk5JgElkNLQZWNbTHWNnexprmT57d1D1pf4BbhqPoSTlpQycmNFRxdXzaqlcTJbJ5N7XE2tMcGzQvMKQ2wtCbMvPLQqOcrNJukcN+Xkea/UnB5ef7o61hXcgKxjHEL0ZfOEc/kOYBn6yEJeFyOKDjC4fQsIj5HOPzuQYIyVuEYrxC8DrwIJIAPqmpm4MTxweZwEIKCKg+ua+O/H91Ma48ZAz1tUSXXnL6Ax7d009SZxOOClfUR8qqUBjxcuKy2P9bAdEA3/cUswAmUmcVj4/TRormUIwKd4IsYERjGTFS3Po0+fzug5Px1dPz8BvAGqP7aQ3gXHj7BZSzFk80XeHlHL2uaOnlmSxcvt/SSH9A2BTwuVswr52Snx7BsdqSoBWSqZtJ4QzTGls5kf55+j4vFVaaXMJr1O1rIoY98F176HSBw4juRqgXgC4MvjPpCpFxBYgToVT+xnItYJt8vFOY4N+h8NMIxp9TP589tHPnGIRivEDzvTBJ/BjNp/PfA3ao6sS4li+RQC8FzW7v4zkObeNWxV15SG+G6cxtxeT38+qU2vG6YHfFRFTaO5MoCHlYtm9XvVGs6oPmcCUbftR05+zrkuMvHl182aaKMpbrBV+KIwNBeTLTpcfTFXwFQmH080Zs+A9k0ZR/5IaGzrhpXOSyHD/F0jue3dbOmuYu1zZ1s2h0fdL086OWEBeX9cwxzK0aeeE7nCmzuiLMhGuufSwCYFfGztCbMgopQUfOHqgrP/BJ94scjV8TtA18I/BGz94XBHx60z3pCpCRAUgIkXAH6NEAfAXoLPrryfrryPvqySiydp77Mz0dOG5vrlvEKwUCz0XOB7wOVqlo7ptKMk0MlBFva49y4ejN/3mAmgqsjPq49u5GzltXw65da6YhnqAgN7rbVhH2cu6TmkAXvniz6zUUr5hrHcu6x93Q0m0A3PwzpHvCXGhEYzkR048PounvMyZILaP/hF8hHtxI6772UfeC7Yy6D5fCnI5bhmS2dzlBSF7t69pl4LvNz8sJKTlxQwdJZEeZXhYYdQlFV2uMZNkRjNHUm+i2GfG6h0eklFOP5V5ueQpufgnTMeDDds8/EIO3sC/t7Eh0T3qARktlHjNmj73iF4C2q+vsB5/OBdxfjhlpEVgHfA9zAT1T1+n2unwX8Dmh2ku4aKd+DLQRd8Qw/eryZ3z63k7wqQa+ba06bx0XHzObltj62OaEk91Di97CwMkRjZWhKO48bClVFn/gRPHM7IMhl30AWnjL2/DJxdPNDkOmDQDnSeB7i3d/VsarC+j+ir/3BJBzzDrrv/C7pF1fjXbSSqq8+iHinXvAey9hQVVq6kmbSec/Ec3LwxLPHJcytDLKoJkxjTZjG6jCLasLM20cgsvkCzZ0J1kdjtMcz/enVIR9LasI0VoXxjXFMXlUhl4FMfLA47LPXdNy5Jw7pofYJ+td81x+D6x++P6byjDUewXJVfV1EhhwCGsnpnIi4MauQzwdagGeAK1X11QH3nAV8ajQTzwdLCFLZPHeubeGWJ/dOBF9+Qj1nLK2mNZamL713AqpQUJbUhDlyVglVoclZpXyoUS2YsdEX7zZ+WFZ9EVl+3tjzy8QcEYhBoAJZdB7i2T+msaZ60Nfuhy1PAoKsvIrYmtXEfnMDrpIqqm94HHf11A0wYxk/BVU2tMZYu6WLF7Z1szkaZ0dXcsioWh6XMK8yxKKaMAtrQiyqibCoJszcyiB96RwbojE2d8TJOFFvPC5hYWWIpTURasKH5retWoBM0giD6pjdZI/V19AngfcD3xqqbMA5I3zuScAmVW1yCnEncCnw6gGfOsQUVPnjujZuciaCy0Je3n7ibBbPjpDI5tnYYcYqs3mlO5ljeW2Ei5bXTNrahMMBLeTQB6836wXcPuTiryKLxh6OQtN9RgSycRNvuPHc/dxSaLwD3fgwbH0KCjkTkvKEa0jv3k3sNzeAuCi/7hYrAhZcIiyvK2F5XQnvOtWMnyezeba0J9gcjdEUjdMUTdDUHmNHV4qm9rgJvPPa3jw8LmF+VYjGatODmF0RII/Sk86xsT3OxvY45UEvS6vDLK4O4/ccvOFeEZeZT/BPnvf/YYVAVd/v7M8eY971wPYB5y3AyUPcd6qIvAjsxPQOXhnj542bZ7d08Z2HN7G9K8mS2REuOHZ2vzvqRNb0CqKxHB2JLEGvh3efMIf544jaNBXQXAa9/6uw6THwBpFLv47MG/Klorj80r2OCCRMhLHGcxD33iE07WtDN/wJtj9j4hoAzDkWWXYh+XSW7v9+KwAlV/wb/mPG+tW0THeCXjdH1JVwRN1ga7ZkNk9ze5ym3UYMNkfjNEXj7OhOsTlqzgcKRFXEx8mLK5lfE6I7mWXt9m6eaelmfnmQ5bUlzC4Zf+jaw4FhhUBE3nagB1X1rhHyHuqvs29v7XlgvqrGROQi4B5gv4gQIvIB4AMA8+ZNbLATgOb2OD/4cxPReIYj55Vy3rF7u15elzC7NMD6aILX2xIocMbCci57w6xJ81N0uKDZFHrvF2DrWhNj4LJvInPeMPb8Ut1mYjiXhHAtsvAcxG2EVnta0PV/gh0vAArigrknIksvQErr0HSCrv88H4134z/hzYQv/fgE1dIykwh63RxZV8qRdYODGiUzRiD2CMPmdrPf2Z3i/r+14hJYWBvmqIZS5lWH2NKVZEtXkmyugCiUBTzMrTDiUDEF44wfaI7gZwd4TlX1vQfMWORU4Muq+ibn/HPOg8NOeYvIFuAEVW0f7p6JnCNo7Ulx53Mt9GVyzKncuwzdLcLc8iCNVSGisSz/92IryWyBiM/NVSvrOLru4Mc1PdhoOobe81nY8SIEy5G3f3tcUZs02YU2PQS5NERmIwvOQtxe4yxu/R+hdZ25Udww/xRk6flIuNo8q0rPDz5M8rH/xT27kerrH8MVGn+4RItlJBKZHM3tCSMOjkjs6k1RXerniPoSSvZZg5AvKJ2xDOlMHp/LRXXYx7yKIPOrzTzEobQgHNMcgaq+Z5yf+wywREQWAjuAK4BB0UpEZDbQpqoqIicBLqBjnJ97QDKOLfFTzZ2oQFmJjzJ8qCqzIn6W10aYWx4kV1B+/WIba7ebODxHzYpw1co6SqeBt9CR0GQPetcnoW09RGqRy7+NVM4fe36JDrNOYI8r6flnQmczhfUPQnS9ucnthQWnI0vOHeRcTnNZ4vfdRPKx/wVfkIpP3W5FwHLQCPk8HDWnlKPmDO5BJDI5Nu82PYe2vjSJXB6P20Uk6KGmdPB81+udcf7c1EFrd4pMNk/E72FuRYh5lUHmV4WYXxmirnziwuKOhaJaNRF5M3AU0G/WMZKZp6rmROSjwIMY89FbVPUVEfmQc/1m4HLgwyKSA5LAFVpsyLRRsq07wcZonG3dSVMnt6AFpTeR5cS55Rw/r6J/qGdTe4Jbn9tJZyKL1y287ehZnLGgfFqMBY6ExtrR337CBJkpm4Nc/l2kbOxxfjXRblYM5zNQUg/BanjiRrTTsRj2BKDxTGTx2YPcSRTi3SRW30r8gZspdLQAUP7BG/HOO2pc9bNYJoKQz8PRDWUc3TD4pSSVzdPUEae5I0FHIkO2oEQCHhbPjrB4tlkpny8o7X1pnm3p5r51rbR2p0ik89SXB5lXZcRhXuVeoag9CPMQxawjuBkIAWcDP8E03mtV9Z8mtWTDMNahoXte3kVXKouqsqMzSW88y1tXzOH0RVX99+QKyv2vRXloQwcKzC0PcM0Jc5hVMjNs1LVnF/qbj0PPDqhagLz9O0ikeuz5xXcbB3KFLPhKYddr0GMadbwhZPHZ0Ph3g3y853ZvIX7/D0k+chuaMgFF3PVLiVz2KUJnXjGu+lksBxtVpSeVIxpP09aXZldvilhm/0VmiXSO1u4UrT1pWrtT7O5JkXVMWANeF/MqTc/hiLoSrjl9bL3zcYWqBE5T1WNE5CVV/YqIfAsYaaL4sOOF5i62dSfpimV43+kLuOiY2YN8lbT2pfnFszvZ3p1CgAuWVnHRETWHrUvZiUY7txkRiO2GWcuQt/3XuILPa6zNuJIu5CAVh6a1gJrVw0vOhYVnDDIZzWxYQ/y+75Nac2+/tZDv6LMIv/la/CvOt95ELVMSEaE86KU86GVJtekRZPMF2uMZdsfSRJ09eGicFaFxlrlHVYmn8uzsSrK1PUFrT4oNbTHa+lJjFoIDUYwQJJ19QkTmYMbwF054SSaZ952xgAdf2c1Vp8wdNGGjqjzR3M1d69rI5pXKkJd3HT+HxdVTNxLRaNHoJiMCyW6oPwZ56w2If/i4wCNR6NkOWx4DFPqiEG2CUCWy5HwzEbzHUiifI7X298Tvu4nsxmfMw24vwdP/gfDFH8W7YGKC3FgshxNet4u60gB1pWakXVXpS+f6RSEay9CZzBAJelgaLGHpHDNkKkBwkiwVixGC+0SkHPgmxtxTMUNEU4q5lSHe98YFg9J6Uzluf34Xr7SZIYiT5pbxjmNnTTvfQAdCd65D7/608ZMy/0Tkkq8j3v1X+BaVVy6NbloNqagx/+zdDak4svJqYwrqMn/XQqKXxCO3knjgZvLRbQBIuJzQBe8j/Kb3464c+5yExTLVEBFKA15KA14WVZlFY7l8gfZEhmgsQzSeZncsQzKbp3yS3NaMKASq+h/O4W9F5D4goKo9k1Kag8jLu/q4/fldxDJ5Ql4XV6yoY2VD6cgPTiN023Po7z4H2SQsPhO56EuIZ/RfNM0moekxtOU5qGowIpCMwYKzkIbjzMpIIBfdRuKBH5FY/Qs02QuAu24R4TdfS/DMK3EFJm/lpMUylfC4XcwuCTC7ZG+vIZ7Jk80XRnhyjJ830g2Oz6A3Awv23C8iqOq3J6VEk0w6V+Cul9t4cks3AEtrQvzj8XOoGIVP8umAbn4Sve/fjTXPEW9C3vRZxDXKiFDpGLr5UWh6HHwBqF1sRCBYDUdfhcsZ189sepb4fTeRevp3/d4YfUeeQfjij+JfucqO/1ssIyAikxrZsJicfw+kgJeByZGjg4C+/kcyW9fSnoITCy6Ow0NVJEiVNwyveCm4PMaW3eU1oRfdXtiTNiCdYtLFfVibmurrD6N//JpplI99K3LOx/vf2oe8v5BzPCHGHK+JcWP+ueVJIyThKqhdBCJQcyRStxK0QHLNvWb8f/3TJiO3h8AZf0/k4mvxNh6SuEYWi2UIihGCBlU9ZtJLMslsb91FQ2I3g0afY842BONbzCCo22uCUuzxI+4NDthCxt3ygPPB14PGudskiEnhpXvhke+Ayw3HXALLz4WtTxuX0P0NfWyA69wY5FLDZ1h/LPidOYXao6F8EYk//oj4/T8g37bF/DVCZYTOew/hCz+Iu6p+wutksVjGRzFC8ICIXKCqf5r00kwir1b8Hbd2LeHUhghnzg/jKoFvdwAADgdJREFUIQ/5rDFvzGed42x/mu6bNvDe/rTcENczxvwxnzFbJgbx/cszotCIC91PIPaei2+fNE/A+Tzz5q4DG/Q9DXy8w9wz31mU1dMEa5qKKosJxRdxoipFIFAGlfOgc4OpT+lC4qvvIvHwz9CEmUJyz1pA+KKPEDz7alyBsVshWSyWyaUYIXgauFvM2EEWY8WkqjqlZlbPO3ohi+tnF20WOp53cS3kjSjk0mYiNps03jb7j5NmgnWI9P7zfNY03pmhuyxj7rGIGH8+oQrToPvCThi9iDEZHdjY+8LgLwFvYL+hI21fj+5YC0CyaTO9v7rOCCPgXX4qkYs/iv+Ei/othSwWy+FLMULwLeBU4OXJcv9wMPC45KCtDRCX2wy9ePwQGFovRxIazWfNkExmaKEYLCQpc+z2Oo162Lhr8IVRXxheXw2vPWR6Kmd/DNfRbxlTvTSfgd4WtGMzxFsB6PvTb0g88yi43AROezvhi6/Ft3jsbqotFsvBpxgh2Aism8oiMBWRPRPQ/qE9nRbTY9FCHlZ/G17+vYkqdtGXkKXF+fDXTIpcWzP5XRvQvp243Fk8ZaX9b/hayNP3p1+TevVFwhf/M+GLPmSDxFgsU5RihGAX8GcReQBI70mcquajMwXN59AHvw6vP2Qmnt/yNaTx1P3uyUe3kWvdTH7nJnK7zFbo3IWnupLA8hX4FixHnKXxWiiQ2bKB1MaXyfXGCZ78Vkr/+Q5coSk1SmixWPahGCFodjafs1kOczSXRv/wZdj8BHiD6DmfIZvMkX/4Z+R2biLXupncrk3kW5vNXAQgoQiBZSsIHbUC34LLBrz5F8j1xSgU/EjZPDynnU3ZWxqs7b/FMo04oBA4i8kiqvrpg1QeyxjQfI5CrIt8WzO5ltfwvnoX7kwnhYLSvW4D2b8MHWzOFS4lcPwqAsuPw1Nd22+uqghEZiHlC5CyubiHCCpvsVimDwcUAlXNi8jKg1UYC2gmSaGv02yxzr3HfR3owLQB1zRuVkmL20X58kbcJWHymSzdrzWRS6ZwlVThrluEp24x3oaleObU4wl6kfwAu1ZxQaQOKZ+PlDbsF0zeYrFMX4oZGvqbiNwL/JoBFvFFxCye0WihgKb6BjfqvR3owMZ9UEPficY60XRi9B8mgqu0korFDXi8UHAHyK18J2VvPwFP3SLE54eerWj3NkhEAWeNg7igZA5SNh/KGgYFkbdYLDOHYoSgEuN6+pwBacoUjEmwB83n0HQCzSTNPp1EMwk0lTD7dHLA9eTI9+57LZ2ETHLkggyF24urtApXpBJXSSWuSCVSUomrpMqcO2mukkpcwTCuxG7o2AyvPww9O6G8Affl38XtD0HPNnTXGjQxIAS0uKG0HimbZ/a28bdYZjzFeB8db+ziw4LeWz9P4pHbzBu3M0E62UgggqukEolU4iqp2NuY72nIBzTw4qRJIDKsawlN9cHOl9CWF2HDPbB7fb8TNwDqjoI3vhdtfRaSnXvTXW4oaUDK50FJfX88AIvFYoHivI82AP8NnI7pCTwBXKeqLZNctglFc9l+1weIC/GHEH8Q8YfN3hfcm+bbZ+8P7X/NHxr0DL4QrkCo/zrewLgtazTeCTteNA1/22vQ1wa+4N5t0UlQWgvhSnOuOehwgsG7PFDaYIZ9SueM2rOoxWKZORQTs/gh4H+B25ykq4GrVPX8SS7bkIw1ZnEh0QuFPOIPgWdyHLqNBtWCcUGRS5oVxNkkGouamL6JDrOiWAT8jk+hYsrr8pqx/rL5UFJnG3+LxdLPeGMW16jqzwac/1xEPjYxRTt4TNaiJ9XCAGd0OWdzjnMZyCXRXMo0+NnU3oY/l2ZYj0GBkNn2fogZ2/dF9jqY8wYRj3PsMT0Q/GXWt4/FYhk1xQhBu4hcDdzhnF+JmTyeUqjq4Ea6kNvrPbS/ATfnOky6uX+fxl7HEaIh40wqZxJmn88br55lc6CqEaoXG0dwnv2dvlksFstEUYwQvBe4CfgO5hX2r07alEJbnobOTZOQs5jx+P7ANs4xLmOiGe+Gnl3Q3gzxzr2NfzZpGv2GY5H6FdBwLFQ32jd6i8Vy0CnGamgbcMlBKMvksqeB7t+8gxtwcZtxeFUoFJzNGfLJZcyWTZktk4B0AtJO8JZMYsAWd0xHhxj2idRCw0qkwWn4K+Yd8rkKi8ViGVYIROTfD/CcDghqPywisgr4HuAGfqKq1+9zXZzrFwEJ4BpVfb6Ygo+a7g7Y2Wwa8KzTiGcTe88zifEN8+yLN2hCONYfi8xdYSJ5lc62Db/FYjnsOFCPYIi4WoSBfwKqgAMKgeOn6PvA+UAL8IyI3Kuqrw647UJgibOdDPzQ2U88ra/BKw8c+B6P34nEFXK2gcfmXPrTgwOuO3u/c+wN2iEei8UyZRhWCFT1W3uORaQEuA54D3AnJljNSJwEbFLVJiePO4FLgYFCcClwqxPr4GkRKReROlXdNeqajIAsPxdmLR2ikd+zD1pzS4vFMiMZyftoJfAJ4CrgF8BKVe0qMu96YPuA8xb2f9sf6p56TAyECUXqj4H6YyY6W4vFYpnyDGuTKCLfBJ4B+oCjVfXLoxABGDqI1r4zqMXcg4h8QESeFZFno9HoKIpgsVgslpE4kHH6J4E5wBeBnSLS62x9ItJbRN4twMDYhQ3AzjHcg6r+WFVPUNUTampqivhoi8VisRTLsEKgqi5VDapqiaqWDthKVLWYZbrPAEtEZKGI+IArgHv3uede4F1iOAXomYz5AYvFYrEMz6TNjqpqTkQ+CjyIMR+9RVVfEZEPOddvBu7HmI5uwpiPTgtPpxaLxTKVmFQzGVW9H9PYD0y7ecCxAtdOZhksFovFcmBG9D56uCEiUWDrGB+vBtpHvGvqMp3rZ+s2dZnO9ZtKdZuvqkNOsk45IRgPIvLscG5YpwPTuX62blOX6Vy/6VI369LSYrFYZjhWCCwWi2WGM9OE4MeHugCTzHSun63b1GU6129a1G1GzRFYLBaLZX9mWo/AYrFYLPswY4RARFaJyHoR2SQinz3U5RkPIjJXRB4VkddE5BURuc5JrxSRh0Rko7OvONRlHSsi4haRF0TkPud8OtWtXER+IyKvO/+Hp06X+onIx53v5DoRuUNEAlO5biJyi4jsFpF1A9KGrY+IfM5pY9aLyJsOTalHz4wQggGxES4EjgSuFJEjD22pxkUO+KSqHgGcAlzr1OezwGpVXQKsds6nKtcBrw04n051+x7wR1VdDhyLqeeUr5+I1AP/Apygqm/AeBS4gqldt58Dq/ZJG7I+zm/wCuAo55kfOG3PYc+MEAIGxEZQ1QwmpsKlh7hMY0ZVd+2J5KaqfZiGpB5Tp184t/0CeOuhKeH4EJEG4M3ATwYkT5e6lQJnAj8FUNWMqnYzTeqH8VYQFBEPEMI4kZyydVPVx4HOfZKHq8+lwJ2qmlbVZozrnJMOSkHHyUwRguHiHkx5RGQBcBywBpi1x2mfs689dCUbF98FPgMMjB06XerWCESBnzlDXz8RkTDToH6qugP4L2AbJqZIj6r+iWlQt30Yrj5Ttp2ZKUJQVNyDqYaIRIDfAh9T1WJcgx/2iMjFwG5Vfe5Ql2WS8AArgR+q6nGYkLBTaahkWJyx8kuBhRgX9mERufrQluqgMmXbmZkiBEXFPZhKiIgXIwK3q+pdTnKbiNQ51+uA3YeqfOPgdOASEdmCGcI7R0R+yfSoG5jvYouqrnHOf4MRhulQv/OAZlWNqmoWuAs4jelRt4EMV58p287MFCEoJjbClEFEBDPG/JqqfnvApXuBdzvH7wZ+d7DLNl5U9XOq2qCqCzD/T4+o6tVMg7oBqGorsF1EljlJ52LieE+H+m0DThGRkPMdPRczfzUd6jaQ4epzL3CFiPhFZCGwBFh7CMo3elR1RmyYuAcbgM3AFw51ecZZlzMwXc6XgL8520VAFcaKYaOzrzzUZR1nPc8C7nOOp03dgBXAs87/3z1AxXSpH/AV4HVgHXAb4J/KdQPuwMx3ZDFv/P90oPoAX3DamPXAhYe6/MVudmWxxWKxzHBmytCQxWKxWIbBCoHFYrHMcKwQWCwWywzHCoHFYrHMcKwQWCwWywzHCoFlxiMiKiK3DTj3iEh0gOfTS8bjsVZEPiYioYkoq8UyGVghsFiMm4c3iEjQOT8f2LHnoqreq6rXjyP/j2EcsFkshyVWCCwWwwMYj6cAV2IWEgEgIteIyE3O8c9F5EYR+auINInI5U76WXt6EM75Tc5z/4Lxu/OoiDzqXLtARJ4SkedF5NeOzyhE5HoReVVEXhKR/zootbZYsEJgsezhTox7gABwDMab63DUYVZ3XwwcsKegqjdi/M2crapni0g18EXgPFVdiVlh/AkRqQQuA45S1WOAr423QhZLsXgOdQEslsMBVX3Jcel9JXD/CLffo6oF4FURmTXKjzoFExzpSeOOBx/wFNALpICfiMgfgPuGzcFimWCsEFgse7kX40//LIw/meFIDzje43o4x+AedmCYZwV4SFWv3O+CyEkYR21XAB8Fzimq1BbLOLFDQxbLXm4BvqqqL4/h2a3AkY7nyTJMg76HPqDEOX4aOF1EFgM4njqXOvMEZap6P2ZyecWYa2GxjBLbI7BYHFS1BRNPeCzPbheRX2E8im4EXhhw+cfAAyKyy5knuAa4Q0T8zvUvYsTid84chQAfH2M1LJZRY72PWiwWywzHDg1ZLBbLDMcKgcViscxwrBBYLBbLDMcKgcViscxwrBBYLBbLDMcKgcViscxwrBBYLBbLDMcKgcViscxw/j9GinJAZch+SAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 3))\n",
    "colors = plt.get_cmap('tab20c').colors\n",
    "colors = colors[:len(asvs_tongue)] + colors[4:4+len(asvs_plaque)]\n",
    "\n",
    "for o, color in zip(asvs_tongue + asvs_plaque,\n",
    "                    colors):\n",
    "    series = df_asv_median[o]  \n",
    "    series_m = moving_average(series)\n",
    "    series_m_scale = series_m / series_m.mean()\n",
    "    taxa = ';'.join(saliva_ma.feature_metadata.loc[o, 'taxonomy'].split(';')[-2:])\n",
    "    \n",
    "    ax.plot(series_m_scale, label=taxa, linewidth=2, color=color)\n",
    "    ax.legend(bbox_to_anchor=(0, 1.7), loc='upper left', fontsize=8)\n",
    "\n",
    "ax.set_xlabel('Minutes')\n",
    "ax.set_ylabel('Normalized relative abundance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:qiime2-2019.10]",
   "language": "python",
   "name": "conda-env-qiime2-2019.10-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "462.222px",
    "left": "326px",
    "top": "110.926px",
    "width": "285.74px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}