work/9_coi_its/2_phylo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trees and diversity estimates for molecular markers. Env from 62_phylo_reduced.\n",
    "\n",
    "After first run and evaluation, manually resolve problems with selected sequences (on plate fasta level):\n",
    "- VBS00055 (aconitus) its2 trimmed 3' - weak signal, multipeaks\n",
    "- VBS00021,VBS00022,VBS00023 (barbirostris) its2 re-trimmed to the same length - sequences too long, weak signal after 1000 bases\n",
    "- VBS00024 (barbirostris) its2 removed - weak signal, multipeaks\n",
    "- marshalli mis-identifications - remove all samples from analysis\n",
    "- VBS00059, VBS00061 (minimus) coi trimmed 5' - retained variation is true\n",
    "- VBS00145 (sundaicus) is rampae according to its2 and ampseq - removed from analysis\n",
    "- vin.M0004, vin.B0009, vin.M0009 its2 removed - problem with amplification after ca 60 bases\n",
    "\n",
    "After finalisation of sequences, need to re-run BLAST and 1_blast.ipynb\n",
    "\n",
    "Higher than expexted variation left without changes:\n",
    "- carnevalei its2 - legitimately highly variable\n",
    "- coustani, tenebrosus, ziemannii - highly variable unresolved group with two branches according to all markers\n",
    "- hyrcanus - two distinct lineages - VBS00082, VBS00083 and VBS00085,VBS00086 according to all markers\n",
    "- nili - two distinct lineages - nil.237, nil.239 and nil.233, nil.236, nil.238 according to all markers \n",
    "- paludis: Apal-81 within coustani group, while Apal-257 is an outgroup according to all markers\n",
    "- VBS00113 (rampae) coi groups with maculatus (closest relative), while its2 and ampseq group it with rampae\n",
    "- sundaicus coi - legitimately highly variable\n",
    "- theileri its - multipeaks at the middle of all samples result in several bases mis-called, do not change \n",
    "- vagus coi - legitimately highly variable\n",
    "\n",
    "Other:\n",
    "- brohieri, hancocki, demeilloni - single group of closely related individuals according to all markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from Bio import AlignIO, SeqIO\n",
    "from Bio.Phylo.TreeConstruction import DistanceCalculator\n",
    "from collections import defaultdict\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import ete3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in\n",
    "WD = '../../../data/phylo_ampl_dada2/coi_its2/work/'\n",
    "PRD = '../../../data/phylo_ampl_dada2/phylo_reduced'\n",
    "OD = 'data'\n",
    "META = os.path.join(OD, 'species_predictions.csv')\n",
    "FA = os.path.join(WD, 'seqman_fa/plate{}.fas')\n",
    "AMPL_ALN = os.path.join(PRD,'aln_all/{}.fa')\n",
    "AMPL_HAP = os.path.join(PRD,'0_haplotypes.csv')\n",
    "TAXONOMY = '../7_species_id/data/0_taxonomy.csv'\n",
    "# out\n",
    "ALN = os.path.join(WD, 'phylo/{}.aln.fas')\n",
    "TREE = os.path.join(OD, '{}.nwk')\n",
    "TREE_FIG = os.path.join(OD, '{}.pdf')\n",
    "TREE_FIG_UNFILTERED = os.path.join(WD, 'phylo/{}.all.pdf')\n",
    "DIVERSITY = os.path.join(OD, 'diversity.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p {os.path.join(WD, 'phylo')}\n",
    "! mkdir -p {OD}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping of markers used and plates\n",
    "MARKER_PLATES = {\n",
    "    'coi':[1,2],\n",
    "    'its2':[3,4]\n",
    "}\n",
    "# samples with conflicting phylogenetic positions - to be removed from diversity estimates\n",
    "CONFLICT_SAMPLES = ['Amar-3-1','Amar-42','Amar-5','VBS00145']\n",
    "# minimum number of samples per species to be used in diversity plotting\n",
    "MIN_SAMPLES_DIV = 2\n",
    "# species labels containing multiple - hard filtering\n",
    "MULTISP_LABELS = ['brohieri', 'hancocki', 'demeilloni', \n",
    "                  'coustani','tenebrosus','ziemanni',\n",
    "                  'hyrcanus','nili','paludis']\n",
    "# species that are overly diverged - softer filtering\n",
    "BAD_SPP = ['hyrcanus','nili','paludis']\n",
    "MULTISP_LABELS = ['Anopheles_'+l for l in MULTISP_LABELS]\n",
    "BAD_SPP = ['Anopheles_'+l for l in BAD_SPP]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metadata "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(144, 44)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meta = pd.read_csv(META, index_col=0)\n",
    "meta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(123, 44)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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>partner_species</th>\n",
       "      <th>amplicons_recovered</th>\n",
       "      <th>ampseq_species</th>\n",
       "      <th>bold_best_id</th>\n",
       "      <th>bold_top_%</th>\n",
       "      <th>bold_low_%</th>\n",
       "      <th>COI_length</th>\n",
       "      <th>COI_num_alignments</th>\n",
       "      <th>COI_species_predictions</th>\n",
       "      <th>COI_top10_species</th>\n",
       "      <th>...</th>\n",
       "      <th>ITS2_first_hit_score</th>\n",
       "      <th>ITS2_second_hit</th>\n",
       "      <th>ITS2_second_hit_species</th>\n",
       "      <th>ITS2_second_hit_hsps</th>\n",
       "      <th>ITS2_second_hit_length</th>\n",
       "      <th>ITS2_second_hit_identities</th>\n",
       "      <th>ITS2_second_hit_identity</th>\n",
       "      <th>ITS2_second_hit_score</th>\n",
       "      <th>COI_seqid</th>\n",
       "      <th>ITS2_seqid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Abro-21</th>\n",
       "      <td>Anopheles_brohieri</td>\n",
       "      <td>21.0</td>\n",
       "      <td>Anopheles_brohieri, Anopheles_demeilloni, Anop...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>355.0</td>\n",
       "      <td>gi|374676298|gb|JN994151.1|</td>\n",
       "      <td>Anopheles_theileri</td>\n",
       "      <td>1.0</td>\n",
       "      <td>481.0</td>\n",
       "      <td>438.0</td>\n",
       "      <td>91.060291</td>\n",
       "      <td>344.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>bro.21_A10-ITS2A.ab1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 44 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            partner_species  amplicons_recovered  \\\n",
       "Abro-21  Anopheles_brohieri                 21.0   \n",
       "\n",
       "                                            ampseq_species bold_best_id  \\\n",
       "Abro-21  Anopheles_brohieri, Anopheles_demeilloni, Anop...          NaN   \n",
       "\n",
       "         bold_top_%  bold_low_%  COI_length  COI_num_alignments  \\\n",
       "Abro-21         NaN         NaN         NaN                 NaN   \n",
       "\n",
       "        COI_species_predictions COI_top10_species  ... ITS2_first_hit_score  \\\n",
       "Abro-21                     NaN               NaN  ...                355.0   \n",
       "\n",
       "                     ITS2_second_hit  ITS2_second_hit_species  \\\n",
       "Abro-21  gi|374676298|gb|JN994151.1|       Anopheles_theileri   \n",
       "\n",
       "         ITS2_second_hit_hsps  ITS2_second_hit_length  \\\n",
       "Abro-21                   1.0                   481.0   \n",
       "\n",
       "         ITS2_second_hit_identities  ITS2_second_hit_identity  \\\n",
       "Abro-21                       438.0                 91.060291   \n",
       "\n",
       "        ITS2_second_hit_score COI_seqid            ITS2_seqid  \n",
       "Abro-21                 344.0       NaN  bro.21_A10-ITS2A.ab1  \n",
       "\n",
       "[1 rows x 44 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# metadata table for use in diversity estimates - only include ampseq'd and non-conflicting samples\n",
    "# remove non-ampseq'd\n",
    "ampseq_meta = meta[~meta.ampseq_species.isna()]\n",
    "# remove conflicting samples \n",
    "# ampseq_meta = ampseq_meta[~ampseq_meta.index.isin(CONFLICT_SAMPLES)]\n",
    "display(ampseq_meta.shape)\n",
    "ampseq_meta.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(584.2857142857143, 87.54966966782146)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# basic stats\n",
    "meta[~meta.ampseq_species.isna()].COI_length.mean(), meta[~meta.ampseq_species.isna()].COI_length.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment and phylogeny"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nthread = 0\n",
      "nthreadpair = 0\n",
      "nthreadtb = 0\n",
      "ppenalty_ex = 0\n",
      "stacksize: 8192 kb\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "Gap Penalty = -1.53, +0.00, +0.00\n",
      "\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 3 ambiguous characters.\n",
      "  101 / 127\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree (efffree=0) ... \n",
      "  120 / 127\n",
      "done.\n",
      "\n",
      "Progressive alignment 1/2... \n",
      "STEP   126 / 126 \n",
      "done.\n",
      "\n",
      "Making a distance matrix from msa.. \n",
      "  100 / 127\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree (efffree=1) ... \n",
      "  120 / 127\n",
      "done.\n",
      "\n",
      "Progressive alignment 2/2... \n",
      "STEP   126 / 126  h\n",
      "done.\n",
      "\n",
      "disttbfast (nuc) Version 7.407\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "0 thread(s)\n",
      "\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-2 (Fast but rough)\n",
      " Progressive method (guide trees were built 2 times.)\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "nthread = 0\n",
      "nthreadpair = 0\n",
      "nthreadtb = 0\n",
      "ppenalty_ex = 0\n",
      "stacksize: 8192 kb\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "Gap Penalty = -1.53, +0.00, +0.00\n",
      "\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "  101 / 126\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree (efffree=0) ... \n",
      "  120 / 126\n",
      "done.\n",
      "\n",
      "Progressive alignment 1/2... \n",
      "STEP   123 / 125  h\n",
      "Reallocating..done. *alloclen = 3196\n",
      "STEP   125 / 125  h\n",
      "done.\n",
      "\n",
      "Making a distance matrix from msa.. \n",
      "  100 / 126\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree (efffree=1) ... \n",
      "  120 / 126\n",
      "done.\n",
      "\n",
      "Progressive alignment 2/2... \n",
      "STEP   123 / 125  h\n",
      "Reallocating..done. *alloclen = 3375\n",
      "STEP   125 / 125  h\n",
      "done.\n",
      "\n",
      "disttbfast (nuc) Version 7.407\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "0 thread(s)\n",
      "\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-2 (Fast but rough)\n",
      " Progressive method (guide trees were built 2 times.)\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for marker, plates in MARKER_PLATES.items():\n",
    "    ! cat {FA.format(plates[0])} {FA.format(plates[1])} | mafft - > {ALN.format(marker)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FastTree Version 2.1.10 Double precision (No SSE3)\n",
      "Alignment: ../../../data/phylo_ampl_dada2/coi_its2/work/phylo/coi.aln.fas\n",
      "Nucleotide distances: Jukes-Cantor Joins: balanced Support: SH-like 1000\n",
      "Search: Normal +NNI +SPR (2 rounds range 10) +ML-NNI opt-each=1\n",
      "TopHits: 1.00*sqrtN close=default refresh=0.80\n",
      "ML Model: Jukes-Cantor, CAT approximation with 20 rate categories\n",
      "Ignored unknown character n (seen 2 times)\n",
      "Ignored unknown character w (seen 1 times)\n",
      "Initial topology in 0.03 seconds\n",
      "Refining topology: 27 rounds ME-NNIs, 2 rounds ME-SPRs, 13 rounds ML-NNIs\n",
      "Total branch-length 1.828 after 0.33 sec, 101 of 200 nodes   \n",
      "ML-NNI round 1: LogLk = -7673.123 NNIs 21 max delta 9.36 Time 0.50\n",
      "Switched to using 20 rate categories (CAT approximation)1 of 20   \n",
      "Rate categories were divided by 0.851 so that average rate = 1.0\n",
      "CAT-based log-likelihoods may not be comparable across runs\n",
      "Use -gamma for approximate but comparable Gamma(20) log-likelihoods\n",
      "ML-NNI round 2: LogLk = -6618.283 NNIs 8 max delta 1.73 Time 0.61\n",
      "ML-NNI round 3: LogLk = -6617.109 NNIs 2 max delta 0.49 Time 0.65\n",
      "ML-NNI round 4: LogLk = -6617.108 NNIs 0 max delta 0.00 Time 0.68\n",
      "Turning off heuristics for final round of ML NNIs (converged)\n",
      "ML-NNI round 5: LogLk = -6616.090 NNIs 0 max delta 0.00 Time 0.81 (final)\n",
      "Optimize all lengths: LogLk = -6616.087 Time 0.83\n",
      "Total time: 1.02 seconds Unique: 101/127 Bad splits: 0/98\n",
      "FastTree Version 2.1.10 Double precision (No SSE3)\n",
      "Alignment: ../../../data/phylo_ampl_dada2/coi_its2/work/phylo/its2.aln.fas\n",
      "Nucleotide distances: Jukes-Cantor Joins: balanced Support: SH-like 1000\n",
      "Search: Normal +NNI +SPR (2 rounds range 10) +ML-NNI opt-each=1\n",
      "TopHits: 1.00*sqrtN close=default refresh=0.80\n",
      "ML Model: Jukes-Cantor, CAT approximation with 20 rate categories\n",
      "Initial topology in 0.08 seconds\n",
      "Refining topology: 27 rounds ME-NNIs, 2 rounds ME-SPRs, 14 rounds ML-NNIs\n",
      "Total branch-length 7.889 after 0.81 sec, 201 of 214 nodes   \n",
      "ML-NNI round 1: LogLk = -16941.037 NNIs 33 max delta 3.76 Time 1.31nges (max delta 3.758)   \n",
      "Switched to using 20 rate categories (CAT approximation)\n",
      "Rate categories were divided by 0.770 so that average rate = 1.0\n",
      "CAT-based log-likelihoods may not be comparable across runs\n",
      "Use -gamma for approximate but comparable Gamma(20) log-likelihoods\n",
      "ML-NNI round 2: LogLk = -16096.686 NNIs 18 max delta 6.66 Time 1.68nges (max delta 6.661)   \n",
      "ML-NNI round 3: LogLk = -16096.211 NNIs 4 max delta 0.00 Time 1.85\n",
      "Turning off heuristics for final round of ML NNIs (converged)\n",
      "ML-NNI round 4: LogLk = -16095.651 NNIs 0 max delta 0.00 Time 2.25 (final)\n",
      "Optimize all lengths: LogLk = -16095.616 Time 2.34 \n",
      "Total time: 2.93 seconds Unique: 108/126 Bad splits: 0/105ternal splits   \n"
     ]
    }
   ],
   "source": [
    "for marker in MARKER_PLATES:\n",
    "    ! fasttree -nt {ALN.format(marker)} > {TREE.format(marker)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['zie.1032_A04-HCO2198.ab1', 'zie.70_C04-HCO2198.ab1', 'ten.79_F04-HCO2198.ab1', 'pal.81_C09-HCO2198.ab1', 'ten.333_A05-HCO2198.ab1', 'zie.1055_E04-HCO2198.ab1', 'zie.77_D04-HCO2198.ab1', 'ten.954_B05-HCO2198.ab1', 'cou.80_A03-HCO2198.ab1', 'cou.71_H02-HCO2198.ab1', 'zie.334_B04-HCO2198.ab1', 'cou.962_C03-HCO2198.ab1', 'cou.959_D03-HCO2198.ab1', 'ten.185_G04-HCO2198.ab1', 'cou.956_B03-HCO2198.ab1', 'ten.191_H04-HCO2198.ab1', 'pal.10_B09-HCO2198.ab1', 'pal.4_A09-HCO2198.ab1', 'pal.257_D09-HCO2198.ab1', 'VBS00084_D03-HCO2198.ab1', 'VBS00083_C03-HCO2198.ab1', 'VBS00082_B03-HCO2198.ab1', 'VBS00086_F03-HCO2198.ab1', 'VBS00085_E03-HCO2198.ab1', 'nils.7_H12-HCO2198.ab1', 'nil.239_B02-HCO2198.ab1', 'nil.237_H01-HCO2198.ab1', 'nil.236_G01-HCO2198.ab1', 'nil.233_F01-HCO2198.ab1', 'nil.238_A02-HCO2198.ab1', 'car.191_A01-HCO2198.ab1', 'car.192_D02-HCO2198.ab1', 'car.195_A06-HCO2198.ab1', 'car.193_E02-HCO2198.ab1', 'car.194_B05-HCO2198.ab1', 'VBS00025_A02-HCO2198.ab1', 'VBS00023_G01-HCO2198.ab1', 'VBS00024_H01-HCO2198.ab1', 'VBS00022_F01-HCO2198.ab1', 'VBS00150_G05-HCO2198.ab1', 'VBS00049_C02-HCO2198.ab1', 'VBS00048_B02-HCO2198.ab1', 'VBS00122_E04-HCO2198.ab1', 'vin.M0012_H08-HCO2198.ab1', 'vin.M0004_H06-HCO2198.ab1', 'vin.B0004_G11-HCO2198.ab1', 'vin.B0009_H11-HCO2198.ab1', 'vin.M0009_A07-HCO2198.ab1', 'dur.97.6_F02-HCO2198.ab1', 'dur.97.5_E02-HCO2198.ab1', 'dur.97.4_D02-HCO2198.ab1', 'dur.97.3_C02-HCO2198.ab1', 'gam.38_E12-HCO2198.ab1', 'gam.37_D12-HCO2198.ab1', 'col.558_E08-HCO2198.ab1', 'col.570_F08-HCO2198.ab1', 'mar.5_E03-HCO2198.ab1', 'col.645_G08-HCO2198.ab1', 'VBS00154_H05-HCO2198.ab1', 'VBS00158_B11-HCO2198.ab1', 'VBS00157_H09-HCO2198.ab1', 'VBS00156_E09-HCO2198.ab1', 'han.615_A08-HCO2198.ab1', 'han.6.1_B08-HCO2198.ab1', 'dem.32_G09-HCO2198.ab1', 'bro.30_C10-HCO2198.ab1', 'dem.19_F09-HCO2198.ab1', 'bro.33_D10-HCO2198.ab1', 'bro.22_B10-HCO2198.ab1', 'mar.3.1_H03-HCO2198.ab1', 'the.6.2_E05-HCO2198.ab1', 'the.6.1_D05-HCO2198.ab1', 'the.2_C05-HCO2198.ab1', 'the.1_F05-HCO2198.ab1', 'rho.7.4_G06-HCO2198.ab1', 'rho.7.3_F06-HCO2198.ab1', 'rho.7.2_E06-HCO2198.ab1', 'rho.7.1_D06-HCO2198.ab1', 'gab.M192_E11-HCO2198.ab1', 'mou.2.3_C01-HCO2198.ab1', 'mou.2.2_B01-HCO2198.ab1', 'mou.2.4_D01-HCO2198.ab1', 'mou.2.1_A01-HCO2198.ab1', 'mar.42_G03-HCO2198.ab1', 'jab.B620_G10-HCO2198.ab1', 'fun.13_A06-HCO2198.ab1', 'fun.3_C06-HCO2198.ab1', 'fun.9_G05-HCO2198.ab1', 'fun.1_B06-HCO2198.ab1', 'VBS00004_D01-HCO2198.ab1', 'VBS00003_C01-HCO2198.ab1', 'VBS00002_B01-HCO2198.ab1', 'VBS00061_A03-HCO2198.ab1', 'VBS00059_H02-HCO2198.ab1', 'VBS00055_G02-HCO2198.ab1', 'VBS00053_F02-HCO2198.ab1', 'VBS00115_D04-HCO2198.ab1', 'VBS00114_C04-HCO2198.ab1', 'VBS00112_A04-HCO2198.ab1', 'VBS00108_H03-HCO2198.ab1', 'VBS00107_G03-HCO2198.ab1', 'VBS00113_B04-HCO2198.ab1', 'mac.50_E07-HCO2198.ab1', 'mac.44_B07-HCO2198.ab1', 'VBS00128_H04-HCO2198.ab1', 'VBS00127_G04-HCO2198.ab1', 'VBS00125_F04-HCO2198.ab1', 'VBS00147_D05-HCO2198.ab1', 'VBS00144_A05-HCO2198.ab1', 'VBS00146_C05-HCO2198.ab1', 'VBS00148_E05-HCO2198.ab1', 'imp.M0030_C11-HCO2198.ab1']\n",
      "['han.615_A08-ITS2A.ab1', 'han.6.1_B08-ITS2A.ab1', 'dem.19_F09-ITS2A.ab1', 'bro.30_C10-ITS2A.ab1', 'bro.22_B10-ITS2A.ab1', 'bro.21_A10-ITS2A.ab1', 'dem.32_G09-ITS2A.ab1', 'bro.33_D10-ITS2A.ab1', 'dem.34_H09-ITS2A.ab1', 'mar.3.1_H03-ITS2A.ab1', 'the.6.1_D05-ITS2A.ab1', 'the.1_F05-ITS2A.ab1', 'the.2_C05-ITS2A.ab1', 'the.6.2_E05-ITS2A.ab1', 'rho.7.3_F06-ITS2A.ab1', 'rho.7.1_D06-ITS2A.ab1', 'rho.7.4_G06-ITS2A.ab1', 'gab.M192_E11-ITS2A.ab1', 'VBS00055_G02-ITS2A.ab1', 'VBS00053_F02-ITS2A.ab1', 'VBS00061_A03-ITS2A.ab1', 'VBS00059_H02-ITS2A.ab1', 'mou.2.5_E01-ITS2A.ab1', 'mou.2.4_D01-ITS2A.ab1', 'mou.2.3_C01-ITS2A.ab1', 'mou.2.2_B01-ITS2A.ab1', 'mou.2.1_A01-ITS2A.ab1', 'mar.42_G03-ITS2A.ab1', 'jab.B620_G10-ITS2A.ab1', 'fun.3_C06-ITS2A.ab1', 'fun.1_B06-ITS2A.ab1', 'gam.38_B09-ITS2A.ab1', 'col.570_F08-ITS2A.ab1', 'gam.35_C08-ITS2A.ab1', 'col.558_E08-ITS2A.ab1', 'gam.32_C07-ITS2A.ab1', 'mar.5_E03-ITS2A.ab1', 'col.645_G08-ITS2A.ab1', 'gam.30_A07-ITS2A.ab1', 'gam.37_A09-ITS2A.ab1', 'VBS00157_E06-ITS2A.ab1', 'VBS00156_A06-ITS2A.ab1', 'VBS00158_H06-ITS2A.ab1', 'VBS00154_H12-ITS2A.ab1', 'VBS00148_E12-ITS2A.ab1', 'VBS00147_D12-ITS2A.ab1', 'VBS00146_C12-ITS2A.ab1', 'VBS00144_A12-ITS2A.ab1', 'VBS00003_C01-ITS2A.ab1', 'VBS00001_A01-ITS2A.ab1', 'VBS00002_B01-ITS2A.ab1', 'VBS00004_D01-ITS2A.ab1', 'VBS00128_H04-ITS2A.ab1', 'VBS00125_F04-ITS2A.ab1', 'VBS00127_G04-ITS2A.ab1', 'VBS00112_A04-ITS2A.ab1', 'VBS00113_B04-ITS2A.ab1', 'VBS00115_D04-ITS2A.ab1', 'VBS00114_C04-ITS2A.ab1', 'VBS00145_B12-ITS2A.ab1', 'VBS00108_H03-ITS2A.ab1', 'VBS00107_G03-ITS2A.ab1', 'mac.50_E07-ITS2A.ab1', 'mac.44_B07-ITS2A.ab1', 'car.194_E02-ITS2A.ab1', 'car.193_A02-ITS2A.ab1', 'car.195_A05-ITS2A.ab1', 'car.191_B11-ITS2A.ab1', 'car.192_D11-ITS2A.ab1', 'nil.233_F01-ITS2A.ab1', 'nil.238_A02-ITS2A.ab1', 'nil.236_G01-ITS2A.ab1', 'nils.7_H10-ITS2A.ab1', 'nil.239_B02-ITS2A.ab1', 'nil.237_H01-ITS2A.ab1', 'VBS00049_C02-ITS2A.ab1', 'VBS00050_D02-ITS2A.ab1', 'VBS00048_B02-ITS2A.ab1', 'VBS00122_E04-ITS2A.ab1', 'VBS00150_G12-ITS2A.ab1', 'VBS00149_F12-ITS2A.ab1', 'dur.97.4_D02-ITS2A.ab1', 'dur.97.5_E02-ITS2A.ab1', 'dur.97.3_C02-ITS2A.ab1', 'vin.M0012_H05-ITS2A.ab1', 'vin.B0004_G11-ITS2A.ab1', 'VBS00021_E01-ITS2A.ab1', 'VBS00022_F01-ITS2A.ab1', 'VBS00023_G01-ITS2A.ab1', 'zie.1032_A04-ITS2A.ab1', 'zie.77_D04-ITS2A.ab1', 'zie.1055_E04-ITS2A.ab1', 'ten.79_F04-ITS2A.ab1', 'cou.71_H02-ITS2A.ab1', 'zie.334_B04-ITS2A.ab1', 'ten.185_G04-ITS2A.ab1', 'cou.962_C03-ITS2A.ab1', 'cou.959_D03-ITS2A.ab1', 'cou.956_B03-ITS2A.ab1', 'ten.191_H04-ITS2A.ab1', 'pal.81_C09-ITS2A.ab1', 'ten.333_A05-ITS2A.ab1', 'zie.70_C04-ITS2A.ab1', 'cou.80_A03-ITS2A.ab1', 'ten.954_B05-ITS2A.ab1', 'pal.257_D09-ITS2A.ab1', 'VBS00082_B03-ITS2A.ab1', 'VBS00083_C03-ITS2A.ab1', 'VBS00085_E03-ITS2A.ab1', 'VBS00086_F03-ITS2A.ab1', 'imp.M0030_C11-ITS2A.ab1']\n"
     ]
    }
   ],
   "source": [
    "# draw trees\n",
    "for marker in MARKER_PLATES:\n",
    "    t = ete3.Tree(TREE.format(marker))\n",
    "    # set outgroup to implexus \n",
    "    outgroups = []\n",
    "    for leaf in t:\n",
    "        if leaf.name.startswith('imp'):# or leaf.name.startswith('cou'):\n",
    "            outgroups.append(leaf.name)            \n",
    "#     print(outgroups)\n",
    "    if len(outgroups) > 1:\n",
    "        t.set_outgroup(t.get_common_ancestor(*outgroups))\n",
    "    elif len(outgroups) == 1:\n",
    "        t.set_outgroup(outgroups[0])\n",
    "    \n",
    "    # style\n",
    "    t.ladderize(direction=1)\n",
    "    ns = ete3.NodeStyle(size=0)\n",
    "    for n in t.traverse():\n",
    "        n.set_style(ns)\n",
    "        \n",
    "    t.render(TREE_FIG_UNFILTERED.format(marker)); \n",
    "    \n",
    "    # remove non-ampseq samples from tree\n",
    "    pruned_taxa = [leaf.name for leaf in t if leaf.name in ampseq_meta[marker.upper() + '_seqid'].to_list()]\n",
    "    print(pruned_taxa)\n",
    "    t.prune(pruned_taxa)\n",
    "            \n",
    "    t.render(TREE_FIG.format(marker)); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'col.554_D08-ITS2A.ab1' in ampseq_meta.COI_seqid.to_list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Diversity estimates for Sanger sequencing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "(0, 0)\n"
     ]
    }
   ],
   "source": [
    "# variable sites\n",
    "def distances(aln, seqids):\n",
    "    s_aln = [seq for seq in aln if seq.name in seqids]\n",
    "    s_aln = AlignIO.MultipleSeqAlignment(s_aln)\n",
    "    aln_len = len(s_aln[0].seq)\n",
    "    dist_matrix = DistanceCalculator('identity').get_distance(s_aln)\n",
    "    dists = []\n",
    "    for i, d in enumerate(dist_matrix):\n",
    "        dists.extend(d[:i])\n",
    "    dists = [int(d * aln_len) for d in dists]\n",
    "    return dists\n",
    "\n",
    "def aln_stats(aln, seqids, verbose=False):\n",
    "    \n",
    "    var_sites = 0\n",
    "    aln_len = 0\n",
    "    if len(seqids) == 0:\n",
    "        return var_sites, aln_len\n",
    "    # subset alignment to seqids\n",
    "    s_aln = [seq for seq in aln if seq.name in seqids]\n",
    "    s_aln = AlignIO.MultipleSeqAlignment(s_aln)\n",
    "    # iterate over alignment columns\n",
    "    for i in range(len(s_aln[0].seq)):\n",
    "        chars = s_aln[:,i]\n",
    "        charset = set(chars)\n",
    "        # any aligned bases\n",
    "        if charset != set('-'):\n",
    "            aln_len += 1\n",
    "        # any variable bases\n",
    "        if len(charset - set('-')) > 1:\n",
    "            if verbose:\n",
    "                print(i, chars)\n",
    "            var_sites += 1\n",
    "            \n",
    "    return var_sites, aln_len\n",
    "\n",
    "aln_div = defaultdict(dict)\n",
    "for marker in MARKER_PLATES:\n",
    "    marker_aln = AlignIO.read(ALN.format(marker), format='fasta')\n",
    "    marker_aln = list(marker_aln)\n",
    "\n",
    "    for species, sp_meta in ampseq_meta.groupby('partner_species'):\n",
    "        # remove conflicting samples\n",
    "        sp_meta = sp_meta[~sp_meta.index.isin(CONFLICT_SAMPLES)]\n",
    "        # subset samples\n",
    "        sp_marker_samples = sp_meta[marker.upper()+'_seqid'].dropna().to_list()\n",
    "        # debug\n",
    "        if species == 'Anopheles_marshallii' and marker=='its2':\n",
    "            print(sp_marker_samples)\n",
    "            print(aln_stats(marker_aln, sp_marker_samples, verbose=True))\n",
    "        # skip small datasets\n",
    "        #if len(sp_marker_samples) < MIN_SAMPLES_DIV:\n",
    "        #    continue\n",
    "        v, l = aln_stats(marker_aln, sp_marker_samples)\n",
    "        aln_div[marker + '_len'][species] = l\n",
    "        aln_div[marker + '_var'][species] = v\n",
    "        aln_div[marker + '_nseq'][species] = len(sp_marker_samples)\n",
    "        #aln_div[marker + '_samples'][species] = sp_marker_samples\n",
    "        #aln_div[marker + '_dists'][species] = distances(marker_aln, sp_marker_samples)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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>coi_len</th>\n",
       "      <th>coi_var</th>\n",
       "      <th>coi_nseq</th>\n",
       "      <th>its2_len</th>\n",
       "      <th>its2_var</th>\n",
       "      <th>its2_nseq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Anopheles_aconitus</th>\n",
       "      <td>619</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>475</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_annularis</th>\n",
       "      <td>619</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>430</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_balabacensis</th>\n",
       "      <td>619</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>458</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_barbirostris</th>\n",
       "      <td>619</td>\n",
       "      <td>10</td>\n",
       "      <td>4</td>\n",
       "      <td>1024</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_brohieri</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>470</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        coi_len  coi_var  coi_nseq  its2_len  its2_var  \\\n",
       "Anopheles_aconitus          619        1         2       475         1   \n",
       "Anopheles_annularis         619        6         3       430         2   \n",
       "Anopheles_balabacensis      619        0         1       458         0   \n",
       "Anopheles_barbirostris      619       10         4      1024         3   \n",
       "Anopheles_brohieri          616        1         3       470         0   \n",
       "\n",
       "                        its2_nseq  \n",
       "Anopheles_aconitus              2  \n",
       "Anopheles_annularis             4  \n",
       "Anopheles_balabacensis          1  \n",
       "Anopheles_barbirostris          3  \n",
       "Anopheles_brohieri              4  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aln_div = pd.DataFrame(aln_div)\n",
    "aln_div.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diversity estimates for ampseq\n",
    "\n",
    "Use alignments generated for 62_phylo_reduced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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>c_Consensus</th>\n",
       "      <th>s_Sample</th>\n",
       "      <th>s_Species</th>\n",
       "      <th>source_id</th>\n",
       "      <th>target</th>\n",
       "      <th>combUID</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TCTCCATTTCTATCTGATGAAACTTGTGCCCGTGGTTTAATTAATG...</td>\n",
       "      <td>Abro-21</td>\n",
       "      <td>Anopheles_brohieri</td>\n",
       "      <td>run1</td>\n",
       "      <td>13</td>\n",
       "      <td>13-0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                         c_Consensus s_Sample  \\\n",
       "0  TCTCCATTTCTATCTGATGAAACTTGTGCCCGTGGTTTAATTAATG...  Abro-21   \n",
       "\n",
       "            s_Species source_id target combUID  \n",
       "0  Anopheles_brohieri      run1     13    13-0  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ampl_haps = pd.read_csv(AMPL_HAP)\n",
    "ampl_haps.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ampseq_var = defaultdict(dict)\n",
    "ampseq_len = defaultdict(dict)\n",
    "ampseq_nseq = dict()\n",
    "\n",
    "for species, sp_meta in ampseq_meta.groupby('partner_species'):\n",
    "    # remove conflicting samples\n",
    "    sp_meta = sp_meta[~sp_meta.index.isin(CONFLICT_SAMPLES)]\n",
    "    # subset samples\n",
    "    sp_marker_samples = sp_meta.index.dropna().to_list()\n",
    "    # skip small datasets\n",
    "    #if len(sp_marker_aln) < MIN_SAMPLES_DIV:\n",
    "    #    continue\n",
    "    ampseq_nseq[species] = len(sp_marker_samples)\n",
    "    for target in range(62):\n",
    "        target = str(target)\n",
    "        st_uids = ampl_haps.loc[ampl_haps.s_Sample.isin(sp_marker_samples) \n",
    "                                    & (ampl_haps.target == target),\n",
    "                                'combUID'].to_list()\n",
    "        # no sequences - no divergent sites\n",
    "        if len(st_uids) == 0:\n",
    "            ampseq_var[target][species] = np.nan\n",
    "        # estimate divergent sites\n",
    "        else:\n",
    "            taln = AlignIO.read(AMPL_ALN.format(target), format='fasta')\n",
    "            ampseq_var[target][species], ampseq_len[target][species] = aln_stats(taln, st_uids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Anopheles_aconitus</th>\n",
       "      <td>131.0</td>\n",
       "      <td>153.0</td>\n",
       "      <td>140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_annularis</th>\n",
       "      <td>138.0</td>\n",
       "      <td>153.0</td>\n",
       "      <td>153</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_balabacensis</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>171</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            0      1    3\n",
       "Anopheles_aconitus      131.0  153.0  140\n",
       "Anopheles_annularis     138.0  153.0  153\n",
       "Anopheles_balabacensis    NaN    NaN  171"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ampseq_var = pd.DataFrame(ampseq_var)\n",
    "ampseq_len = pd.DataFrame(ampseq_len)\n",
    "ampseq_len.iloc[:3,:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "comb_div = aln_div.copy()\n",
    "comb_div['total_ampseq_len'] = ampseq_len.sum(axis=1).astype(int)\n",
    "comb_div['total_ampseq_var'] = ampseq_var.sum(axis=1).astype(int)\n",
    "comb_div['total_ampseq_nsamples'] = pd.Series(ampseq_nseq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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>coi_len</th>\n",
       "      <th>coi_var</th>\n",
       "      <th>coi_nseq</th>\n",
       "      <th>its2_len</th>\n",
       "      <th>its2_var</th>\n",
       "      <th>its2_nseq</th>\n",
       "      <th>total_ampseq_len</th>\n",
       "      <th>total_ampseq_var</th>\n",
       "      <th>total_ampseq_nsamples</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Anopheles_aconitus</th>\n",
       "      <td>619</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>475</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>7712.0</td>\n",
       "      <td>45</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_annularis</th>\n",
       "      <td>619</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>430</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>9068.0</td>\n",
       "      <td>245</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_balabacensis</th>\n",
       "      <td>619</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>458</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>8047.0</td>\n",
       "      <td>28</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_barbirostris</th>\n",
       "      <td>619</td>\n",
       "      <td>10</td>\n",
       "      <td>4</td>\n",
       "      <td>1024</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>8054.0</td>\n",
       "      <td>255</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_brohieri</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>470</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9739.0</td>\n",
       "      <td>62</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_carnevalei</th>\n",
       "      <td>620</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>480</td>\n",
       "      <td>18</td>\n",
       "      <td>5</td>\n",
       "      <td>8773.0</td>\n",
       "      <td>203</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_coluzzii</th>\n",
       "      <td>616</td>\n",
       "      <td>8</td>\n",
       "      <td>3</td>\n",
       "      <td>526</td>\n",
       "      <td>8</td>\n",
       "      <td>3</td>\n",
       "      <td>9336.0</td>\n",
       "      <td>278</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_coustani</th>\n",
       "      <td>616</td>\n",
       "      <td>11</td>\n",
       "      <td>5</td>\n",
       "      <td>586</td>\n",
       "      <td>52</td>\n",
       "      <td>5</td>\n",
       "      <td>8716.0</td>\n",
       "      <td>335</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_demeilloni</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>469</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>9286.0</td>\n",
       "      <td>55</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_dirus_A</th>\n",
       "      <td>619</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>813</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>8218.0</td>\n",
       "      <td>100</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_dureni</th>\n",
       "      <td>616</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>385</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>7196.0</td>\n",
       "      <td>26</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_funestus</th>\n",
       "      <td>616</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>803</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>9552.0</td>\n",
       "      <td>142</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_gabonensis</th>\n",
       "      <td>616</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>461</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>8577.0</td>\n",
       "      <td>50</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_gambiae</th>\n",
       "      <td>402</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>527</td>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>9556.0</td>\n",
       "      <td>382</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_hancocki</th>\n",
       "      <td>616</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>469</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>9385.0</td>\n",
       "      <td>35</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_hyrcanus</th>\n",
       "      <td>619</td>\n",
       "      <td>40</td>\n",
       "      <td>5</td>\n",
       "      <td>615</td>\n",
       "      <td>148</td>\n",
       "      <td>4</td>\n",
       "      <td>8389.0</td>\n",
       "      <td>635</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_implexus</th>\n",
       "      <td>616</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>534</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7659.0</td>\n",
       "      <td>54</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_jamesii</th>\n",
       "      <td>619</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>474</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>9578.0</td>\n",
       "      <td>142</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_jebudensis</th>\n",
       "      <td>213</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>502</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>9154.0</td>\n",
       "      <td>25</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_maculatus_A</th>\n",
       "      <td>619</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>426</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>8563.0</td>\n",
       "      <td>142</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_maculipalpis</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>449</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>8998.0</td>\n",
       "      <td>69</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_marshallii</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_minimus_A</th>\n",
       "      <td>576</td>\n",
       "      <td>16</td>\n",
       "      <td>2</td>\n",
       "      <td>467</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>9018.0</td>\n",
       "      <td>104</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_moucheti</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>491</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>9283.0</td>\n",
       "      <td>82</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_nili</th>\n",
       "      <td>616</td>\n",
       "      <td>23</td>\n",
       "      <td>5</td>\n",
       "      <td>472</td>\n",
       "      <td>48</td>\n",
       "      <td>5</td>\n",
       "      <td>8401.0</td>\n",
       "      <td>353</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_nili_ss</th>\n",
       "      <td>343</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>450</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>4087.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_paludis</th>\n",
       "      <td>616</td>\n",
       "      <td>36</td>\n",
       "      <td>4</td>\n",
       "      <td>605</td>\n",
       "      <td>147</td>\n",
       "      <td>2</td>\n",
       "      <td>5490.0</td>\n",
       "      <td>222</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_rampae</th>\n",
       "      <td>619</td>\n",
       "      <td>32</td>\n",
       "      <td>4</td>\n",
       "      <td>421</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9396.0</td>\n",
       "      <td>335</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_rhodesiensis</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>496</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>9658.0</td>\n",
       "      <td>141</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_sundaicus</th>\n",
       "      <td>619</td>\n",
       "      <td>18</td>\n",
       "      <td>4</td>\n",
       "      <td>561</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9303.0</td>\n",
       "      <td>90</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_tenebrosus</th>\n",
       "      <td>616</td>\n",
       "      <td>11</td>\n",
       "      <td>5</td>\n",
       "      <td>585</td>\n",
       "      <td>56</td>\n",
       "      <td>5</td>\n",
       "      <td>8506.0</td>\n",
       "      <td>356</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_tessellatus</th>\n",
       "      <td>452</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>541</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>8776.0</td>\n",
       "      <td>306</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_theileri</th>\n",
       "      <td>616</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>501</td>\n",
       "      <td>35</td>\n",
       "      <td>4</td>\n",
       "      <td>9729.0</td>\n",
       "      <td>197</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_vagus</th>\n",
       "      <td>619</td>\n",
       "      <td>22</td>\n",
       "      <td>4</td>\n",
       "      <td>656</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9104.0</td>\n",
       "      <td>122</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_vinckei</th>\n",
       "      <td>616</td>\n",
       "      <td>11</td>\n",
       "      <td>5</td>\n",
       "      <td>385</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>9137.0</td>\n",
       "      <td>264</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Anopheles_ziemanni</th>\n",
       "      <td>616</td>\n",
       "      <td>9</td>\n",
       "      <td>5</td>\n",
       "      <td>586</td>\n",
       "      <td>55</td>\n",
       "      <td>5</td>\n",
       "      <td>7824.0</td>\n",
       "      <td>305</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        coi_len  coi_var  coi_nseq  its2_len  its2_var  \\\n",
       "Anopheles_aconitus          619        1         2       475         1   \n",
       "Anopheles_annularis         619        6         3       430         2   \n",
       "Anopheles_balabacensis      619        0         1       458         0   \n",
       "Anopheles_barbirostris      619       10         4      1024         3   \n",
       "Anopheles_brohieri          616        1         3       470         0   \n",
       "Anopheles_carnevalei        620        5         5       480        18   \n",
       "Anopheles_coluzzii          616        8         3       526         8   \n",
       "Anopheles_coustani          616       11         5       586        52   \n",
       "Anopheles_demeilloni        616        1         2       469         1   \n",
       "Anopheles_dirus_A           619        1         2       813         3   \n",
       "Anopheles_dureni            616        0         4       385         2   \n",
       "Anopheles_funestus          616        3         4       803         0   \n",
       "Anopheles_gabonensis        616        0         1       461         0   \n",
       "Anopheles_gambiae           402        0         2       527         4   \n",
       "Anopheles_hancocki          616        0         2       469         0   \n",
       "Anopheles_hyrcanus          619       40         5       615       148   \n",
       "Anopheles_implexus          616        0         1       534         0   \n",
       "Anopheles_jamesii           619        2         3       474         3   \n",
       "Anopheles_jebudensis        213        0         1       502         0   \n",
       "Anopheles_maculatus_A       619        1         2       426         1   \n",
       "Anopheles_maculipalpis      616        1         2       449         0   \n",
       "Anopheles_marshallii          0        0         0         0         0   \n",
       "Anopheles_minimus_A         576       16         2       467         1   \n",
       "Anopheles_moucheti          616        1         4       491         1   \n",
       "Anopheles_nili              616       23         5       472        48   \n",
       "Anopheles_nili_ss           343        0         1       450         0   \n",
       "Anopheles_paludis           616       36         4       605       147   \n",
       "Anopheles_rampae            619       32         4       421         0   \n",
       "Anopheles_rhodesiensis      616        1         4       496         0   \n",
       "Anopheles_sundaicus         619       18         4       561         0   \n",
       "Anopheles_tenebrosus        616       11         5       585        56   \n",
       "Anopheles_tessellatus       452        0         1       541         3   \n",
       "Anopheles_theileri          616        0         4       501        35   \n",
       "Anopheles_vagus             619       22         4       656         0   \n",
       "Anopheles_vinckei           616       11         5       385         1   \n",
       "Anopheles_ziemanni          616        9         5       586        55   \n",
       "\n",
       "                        its2_nseq  total_ampseq_len  total_ampseq_var  \\\n",
       "Anopheles_aconitus              2            7712.0                45   \n",
       "Anopheles_annularis             4            9068.0               245   \n",
       "Anopheles_balabacensis          1            8047.0                28   \n",
       "Anopheles_barbirostris          3            8054.0               255   \n",
       "Anopheles_brohieri              4            9739.0                62   \n",
       "Anopheles_carnevalei            5            8773.0               203   \n",
       "Anopheles_coluzzii              3            9336.0               278   \n",
       "Anopheles_coustani              5            8716.0               335   \n",
       "Anopheles_demeilloni            3            9286.0                55   \n",
       "Anopheles_dirus_A               3            8218.0               100   \n",
       "Anopheles_dureni                3            7196.0                26   \n",
       "Anopheles_funestus              2            9552.0               142   \n",
       "Anopheles_gabonensis            1            8577.0                50   \n",
       "Anopheles_gambiae               5            9556.0               382   \n",
       "Anopheles_hancocki              2            9385.0                35   \n",
       "Anopheles_hyrcanus              4            8389.0               635   \n",
       "Anopheles_implexus              1            7659.0                54   \n",
       "Anopheles_jamesii               3            9578.0               142   \n",
       "Anopheles_jebudensis            1            9154.0                25   \n",
       "Anopheles_maculatus_A           2            8563.0               142   \n",
       "Anopheles_maculipalpis          2            8998.0                69   \n",
       "Anopheles_marshallii            0               NaN                 0   \n",
       "Anopheles_minimus_A             2            9018.0               104   \n",
       "Anopheles_moucheti              5            9283.0                82   \n",
       "Anopheles_nili                  5            8401.0               353   \n",
       "Anopheles_nili_ss               1            4087.0                 3   \n",
       "Anopheles_paludis               2            5490.0               222   \n",
       "Anopheles_rampae                4            9396.0               335   \n",
       "Anopheles_rhodesiensis          3            9658.0               141   \n",
       "Anopheles_sundaicus             4            9303.0                90   \n",
       "Anopheles_tenebrosus            5            8506.0               356   \n",
       "Anopheles_tessellatus           2            8776.0               306   \n",
       "Anopheles_theileri              4            9729.0               197   \n",
       "Anopheles_vagus                 4            9104.0               122   \n",
       "Anopheles_vinckei               2            9137.0               264   \n",
       "Anopheles_ziemanni              5            7824.0               305   \n",
       "\n",
       "                        total_ampseq_nsamples  \n",
       "Anopheles_aconitus                          2  \n",
       "Anopheles_annularis                         4  \n",
       "Anopheles_balabacensis                      1  \n",
       "Anopheles_barbirostris                      5  \n",
       "Anopheles_brohieri                          4  \n",
       "Anopheles_carnevalei                        5  \n",
       "Anopheles_coluzzii                          3  \n",
       "Anopheles_coustani                          5  \n",
       "Anopheles_demeilloni                        3  \n",
       "Anopheles_dirus_A                           3  \n",
       "Anopheles_dureni                            4  \n",
       "Anopheles_funestus                          4  \n",
       "Anopheles_gabonensis                        1  \n",
       "Anopheles_gambiae                           5  \n",
       "Anopheles_hancocki                          2  \n",
       "Anopheles_hyrcanus                          5  \n",
       "Anopheles_implexus                          1  \n",
       "Anopheles_jamesii                           3  \n",
       "Anopheles_jebudensis                        1  \n",
       "Anopheles_maculatus_A                       2  \n",
       "Anopheles_maculipalpis                      2  \n",
       "Anopheles_marshallii                        0  \n",
       "Anopheles_minimus_A                         2  \n",
       "Anopheles_moucheti                          5  \n",
       "Anopheles_nili                              5  \n",
       "Anopheles_nili_ss                           1  \n",
       "Anopheles_paludis                           4  \n",
       "Anopheles_rampae                            4  \n",
       "Anopheles_rhodesiensis                      4  \n",
       "Anopheles_sundaicus                         4  \n",
       "Anopheles_tenebrosus                        5  \n",
       "Anopheles_tessellatus                       2  \n",
       "Anopheles_theileri                          4  \n",
       "Anopheles_vagus                             4  \n",
       "Anopheles_vinckei                           5  \n",
       "Anopheles_ziemanni                          5  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "comb_div"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "comb_div.to_csv(DIVERSITY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>subgenus</th>\n",
       "      <th>section</th>\n",
       "      <th>series</th>\n",
       "      <th>group</th>\n",
       "      <th>subgroup</th>\n",
       "      <th>complex</th>\n",
       "      <th>species_ncbi</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>species</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>Anopheles_aconitus</th>\n",
       "      <td>Cellia</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Myzomyia</td>\n",
       "      <td>Minimus</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Anopheles aconitus</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   subgenus section    series    group subgroup complex  \\\n",
       "species                                                                   \n",
       "Anopheles_aconitus   Cellia     NaN  Myzomyia  Minimus      NaN     NaN   \n",
       "\n",
       "                          species_ncbi  \n",
       "species                                 \n",
       "Anopheles_aconitus  Anopheles aconitus  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "taxonomy = pd.read_csv(TAXONOMY, index_col=0)\n",
    "taxonomy.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "comb_div['series'] = taxonomy['series']\n",
    "comb_div['subgenus'] = taxonomy['subgenus']\n",
    "def region(species):\n",
    "    sp_data = ampseq_meta[ampseq_meta.partner_species == species]\n",
    "    if sp_data.index[0].startswith('VBS'):\n",
    "        return 'SE Asia'\n",
    "    elif sp_data.index[0].startswith('A'):\n",
    "        return 'Africa'\n",
    "    return 'Unknown'\n",
    "comb_div['region'] = comb_div.index.map(region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index(['Anopheles_minimus_A', 'Anopheles_rampae', 'Anopheles_sundaicus',\n",
      "       'Anopheles_vagus'],\n",
      "      dtype='object')\n",
      "Index(['Anopheles_carnevalei', 'Anopheles_coustani', 'Anopheles_tenebrosus',\n",
      "       'Anopheles_theileri', 'Anopheles_ziemanni'],\n",
      "      dtype='object')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwYAAAE+CAYAAAA3cAc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3xb5dn/8c9ledtxHGcDCRmEEQIJQexSwiqzrEIh7FWgpfNp+5TSRZ/+6IAyCi1tmYGWvVOgjFLChuCwQsJISAKE7O0Vz+v3xzkmSiLHsi1Zlvx9v156SbrPObcumXBLl+5l7o6IiIiIiPRuOekOQERERERE0k+JgYiIiIiIKDEQERERERElBiIiIiIighIDERERERFBiYGIiIiIiAC56Q6gqwYMGOAjRoxIdxgi0o1mzJixwt0HpjuOZFAbJtK7ZFP7Jdkn4xODESNGUFlZme4wRKQbmdkn6Y4hWdSGifQu2dR+SfbRUCIREREREVFiICIiIiIiSgxERERERAQlBiIiIiIighIDERERERFBiYGIiIiIiJAFy5UmpL4K1q+DZbMhtxAGjIGi8uCxiEi2q10JdWth5RwoHQRlW0HJIDBLd2QiItKDZH9iULsKXvsrvHQ1tDQFZfmlcMKNMGoS5JekMzoRkdSqXgaPfAvmPrOhrHw4nPYg9N8OctRxLCIigez/RJj/ArxwxYakAKChGu49HaoWpy8uEZFUq6+GZ365cVIAsOZTuP0oqFmWnrhERKRHyu7EoHo5PP+H+Me8JehJaKrv3phERLpL/TqYeX/8Y9XLYOms7o1HRER6tOxODLwZVi9o+/iKOdC4vtvCERHpVo11G/eWbmrFR90Xi4iI9HjZnRjk5MLAHdo+PnRXyCvuvnhERLpTXhHkFrR9fMgu3ReLiIj0eNmdGJQMgAN/Hv9YJB/2+Abk5nVvTCIi3aWwHCaeFf9Y+fBghTYREZFQdicGANtE4cirNl59qM8QOOMRKB2cvrhERFItvxgO+AmMPxUsprkfsguc9ZjaQBER2Uj2L1daVA67nQ47HR2sQhQpgOL+UDwAIpF0RyciklolA+CIP8BBP4fqJUEvQmHfoFxERCRG9icGAHmFkDck6CkQEeltCsuCW9+t0x2JiIj0YNk/lEhERERERNqlxEBERERERJQYiIiIiIiIEgMREREREUGJgYiIiIiIoMRARERERERQYiAiIiIiIigxEBERERERlBhIe2pWQO3KDc/r1kDtqvTFIyIikm7NjdBQE9y3qlsN69emLyaRJFBiIG2rWQEPXwTPXBYkA3Vr4M3b4favQtWSdEcnIiLS/ZobYdls+NOusHRW8HzdIrj5YJj1iJIDyWi56Q5AeqjmRljwEsx9JnjuzVA+HKb9Lnj+5h2wz7cgvzR9MYqIiHS3xjp49OLgx7PbjoBj/wLP/T9Y+TE8eQnseFS6IxTpNPUYSHyRPBg1CQ7+VfD87Ts3JAUTToc9zlNSICIivU9hGUy+FypGQWMtPHBOkBTkFcHZT0BBWbojFOk0JQbStqJy2PMbMGSXmLJ+cOQVUNw/fXGJiIikU9+t4aTbNy778o9h8E6Qm5+emESSIKWJgZkVmtl0M3vHzGaZ2a/D8svM7HMzezu8HRlzzU/NbK6ZfWhmh6UyPmlH3Rp44xZYMjOmbDX8+5KNJySLZCG1XyLSpnWL4IGzNy574Y+w7IONJySLZJhU9xjUAwe5+3hgAnC4me0dHrvG3SeEtycAzGwscAqwM3A4cIOZRVIco8TT3AjzpsF/wqFEu06GA34SPH7rjiBhaKhOW3gi3UDtl4hsbv06uOvrG4YPfe2WDcOKbjsC1q9Jd4QinZbSycfu7kDrt8e88OZbuORY4B53rwfmm9lcYE/g1VTGKXFE8mDEfrDdIVA6BL7yf2ARyC2E9x6AiWdqjoFkNbVfIhJXXhEc82f45wlw2oMwZBwM3wduPxr2/Q5ENJRIMlfK5xiYWcTM3gaWAc+4++vhoW+b2btmdquZ9QvLtgY+i7l8YVgm6VAyEI7/e5AUFPcP5hxEz4Wz/gV9hqQ7OpGUU/slIpuJ5MHgneH77wb3kbxgzsE3/gvjvgaFfdMdoUinpTwxcPdmd58AbAPsaWbjgL8Cowm65xcDV4WnW7wqNi0wswvMrNLMKpcvX56iyAWAkgEbTzQuKtfEY+k1UtF+gdowkYwXyQt6zWMnGhf1U1IgGa/bViVy9zXANOBwd18afuC2ADcRdLdD8AvbsJjLtgEWxanrRnePunt04MCBKY5cRHq7ZLZfYX1qw0REpMdJ9apEA82sPHxcBBwCfGBmQ2NOOx54L3w8FTjFzArMbCQwBpieyhhFROJR+yUiIr1Nqnc+HgrcHq7MkQPc5+6Pmdk/zGwCQTf7AuBCAHefZWb3AbOBJuBid29OcYwiIvGo/RIRkV7FgoU3Mlc0GvXKysp0hyEi3cjMZrh7NN1xJIPaMJHeJZvaL8k+2vlYtqy5CTZNHrV5i4iIZDL34PMt1qbPRXqhVA8lkkzW3ATLP4D8Yug3EsyC3R5XzYetdgvKRbpgbW0DtQ3NNDS3UJwfoaKkgEhOvMV9erc14d+psbmF4vxcKkry9XcS6Sx3WDUPGutg4I4QyYXqZbD4HRi2NxT2SXeEImmjxEDia26C5e/DrYcHm7mc9zTkFsEdxwQN6mkPBA2okgPppBVV9Vz68EyeeX8p7jC4rIDrJ+/G+GHlFORqw+BWy6vW878PvMtzHwbLmm7Vt5C/nDaRnbcqI19/J5GOaU0Kbv0KNNXDOU9Cn8Fwz6mw8A044UbY/ggoLEt3pCJpoaFEEp/lBIlAbgHULIebD4ZbDoUVcyAnDwrKIEdfSqRzahuauPLpD3l69tIvRqotXVfPGbdMZ3WNhqq1ql7fyG8em/1FUgCwaO16Tr3pdVbX6u8k0mFmwc7EecVQXwW3HhbcFr4RfKYV9ddnm/RqSgwkvpwcqBgF5z4dbNpSuwrWfhY0qGc/DkN2CZIGkU6oXt/Ew29+vll5fVMLL81dkYaIeqbq+maemLlks/K6xmYqF6xKQ0QiWaB8GJzzb+i7DTRUw8qPg/LJ98G2+0B+SXrjE0kjJQbStpycoIHML91QFsmH4orgXqSTHGhqaYl7rLZBEwA3cJpa4q8ct74x/t9PRBIQyd94l+KcCJQNhYh+8JLeTYmBtG3dIrjjq2FPQV6QJDRUwy1fgdXzN1+tSCRBRXkRDtxh0GblZsQt760K8yLsO7r/ZuU5BnvHKReRBFQvg7snw9JZwbDZgjJoaQ6GFC3/QKsTSa+mxEDia26ElXODSVq5hcEErQtfhOL+wZyDBS9DQ026o5QMVVaUx/87bhyjB27ojcqLGFeeuCv9itUb1aq8OJ8rTtyVEf03TPIvyM3h2lN2o19RXhojE8lQ7lC9FJa+F/QSnHo/XPw6lA8P5hx88FjwA5hIL6UNzqRtDTXwySvBHIMhuwSTjlfNg09ehp2P27gbVqQTVlTVs7y6ntW1DYwaUEpZUS7F+e0vlpZNGwQl0oYtr6pnWdV61tU1MmpgKX0KE/s7iUgcTfWw+F1Yv3bDnII1n8Gsh2DimcFnXgplU/sl2UefLNK2/BLYdr/gV5XWicYVo4Kl3Qq0zrN03YA+BQzoozG97RnYp4CB+juJJEduAQzdNRg+1LrkdvkwiJ6rzzbp9ZQYyJZtuk9BTo4aThERyWzxVtXTZ5uI5hiIiIiIiIgSAxERERERQYmBiIiIiIigxEBERERERFBiICIiIiIiKDEQERERERGUGIiIiIiICEoMREREREQEJQYiIiIiIoISAxERERERQYlB71Fftfnzlpb0xCIiIpJp4n2OimQZJQa9Qc0KeP4PUL00eF63Bt66E1bNU3IgIiLSnqol8OLVwecpQO1KeOkaWLcovXGJJJkSg2xXuxqmfgdeuR7uOjlICmZMgSd/ArceClWL0x2hiCRZY3ML85ZXc/OL8/hwSRX1jc3pDkkkc9WthjuOhZeuhif+Fxpq4MlL4cWr4I5joHZVuiMUSZrcdAcgKZZXAHteAHOehkVvwbXjNnR/bn8k5BWmNz4RSbpV1Q189fqXqGlopiD3Q57/8YEM6RtJd1gimclyYI/z4YkfwawHYe7TGz5Hdz8nOC6SJVL6r9nMCs1supm9Y2azzOzXYXmFmT1jZnPC+34x1/zUzOaa2Ydmdlgq4+sV8oph+F4w+e7geWtjNvY4+MpvoLh/+mIT6cEyuf1qaG6hpiHoJahvaqG2oSldoYhkvsK+sOvX4ZD/C563fo4e+DOYcDoUlacvNpEkS3WaWw8c5O7jgQnA4Wa2N3AJ8Ky7jwGeDZ9jZmOBU4CdgcOBG8xMP3N1VVMDLJ65cdnqBdDSmJZwRDJExrZffQpz+faB2zGgNJ8z996W8uL8dIQhkj1ammDpexuXLZ2tz1HJOilNDDxQHT7NC28OHAvcHpbfDhwXPj4WuMfd6919PjAX2DOVMWa91jkF/w1/6dhqIuTkwuK3gzkHrROSRWQjmdx+lRfnc9EBo3jiu/vzo8N3oKJEiYFIp9WuDOYUzLwveL7NHsH97IeDOQetE5JFskDKB8aZWcTM3gaWAc+4++vAYHdfDBDeDwpP3xr4LObyhWGZdJY7NITfbXY7A854BE67P0gOmhuC4yISVya3X6WFeQwqK6SsMC9dIYhkB3dorA0eH/ZbOOMhOOLK4HlTXfriEkmBlE8+dvdmYIKZlQMPm9m4LZxu8arY7CSzC4ALAIYPH56UOLNWcT/Y52IYOgG23QeK+sLwveGsx6BiJPQZku4IRXqsVLRfoDZMJKOUDICjr4EJpwWfnwVlMP5k6LctbL17cFwkS3TbVHp3XwNMIxh7u9TMhgKE98vC0xYCw2Iu2wbYbJFgd7/R3aPuHh04cGBK484KRf1gh8M3TDTOKw66QpUUiCQkme1XWJ/aMJFMUjIAxhyyYaJxYV/Y7mAlBZJ1Ur0q0cDwlzbMrAg4BPgAmAqcFZ52FvBo+HgqcIqZFZjZSGAMMD2VMfYaOZt0DkW0Uq3Ilqj9EpGNbPo5uulzkSyQ6n/VQ4Hbw5U5coD73P0xM3sVuM/MzgM+BU4CcPdZZnYfMBtoAi4Ou/JFRLqb2i8REelVzDN88mk0GvXKysp0hyEi3cjMZrh7NN1xJIPaMJHeJZvaL8k+2q5PRERERESUGIiIiIiIiBIDERERERFBiYGIiIiIiKDEQEREREREUGIgIiIiIiIoMRAREREREVK/wZl0t9pV0NwAuYUbtm4XEemqlmaoXQHuUFgOeYXpjkhERJJMPQbZZPUCuPc0+NOucP9ZsOazdEckItlg/Tp4fyrcdBD8ZS948SqoWZnuqEREJMmUGGSL6mVw54nwySvQVA/zpsF9Z0LNinRHJiKZrmoJ3H82rF0I69fAC1fAx8+mOyoREUkyJQbZomk9rJizcdmiN4MkQUSkKz7+7+ZlM++D+uruj0VERFJGiUG2iORDUb+Ny/oMgRxNIxGRLhqyy+ZlW00M5jKJiEjWUGKQLYoq4IQbIa8oeJ5fAl+7BYoHpDcuEcl8A3eAcSdueD5kF4ieBxH98CAikk3UqmeL3HwY8WX47ltQXwUFZUEPQiSS7shEJNOVDIAj/wiHXAYtjVDQB0oGpjsqERFJMiUG2SSvEPKGQp+h6Y5ERLJNcb/gJiIiWUtDiURERERERD0G0o76KqhbEyx/2twAow8MhihtOtFZRLpHUz3UrYbPXod1i2H4XtB3Gw3tERGRLlNiIG2rWwNv3gH/+RV4y4by3c+Fg38BxRXpi02kN2qsh0Uzgj1LGmo2lG81EU69B0oHpy82EenVzOwYYKy7/z7dsUjnaSiRtG3VPHjmFxsnBQAzbo2/rrmIpFbdKvjnCRsnBRDsWfLUL7SvgIgkhQU69B3R3acqKch87f5HN7OImV3ZHcFID9JQAy9f2/bxl6+FmuXdF49IJ2VVGzZvGjTWxT82+yGoX9et4YhI9jCzEWb2vpndALwJnGFmr5rZm2Z2v5mVhucdaWYfmNlLZnadmT0Wlp9tZn8OH29rZs+a2bvh/fCwfEp4zStmNs/MTmwrHkmPdhMDd28Gdjcz64Z4pKdoqod1i9o+XrUYWlraPi7SQ2RVG7ZmQdvHmhuDeUAiIp23A3AHcChwHnCIu08EKoH/MbNC4O/AEe7+JaCtyU1/Bu5w912BO4HrYo4NBb4EHA2oh6GHSXSOwVvAo2Z2P/BFH7a7P5SSqCT98ktg2N6w8I34x4fupl1PJZNkRxs2fN+2jxVX6P9JEemqT9z9NTM7GhgLvBz+ppIPvArsCMxz9/nh+XcDF8SpZx/ghPDxP4ArYo494u4twGwz08SoHibRxKACWAkcFFPmQGZ9qEpca+saycsxiguCfw6raxoAo99eF0LlzZsPXbAcOPiXUNS3+4MV6ZzsaMMG7QQVo4L5P5va/0dQ3J/G5hbqG5spLcz74lDV+kbyc3MoyNWGhyKyRa0/nBjwjLtPjj1oZrt1sl6PeVwfW2Un65MUSWhiibufE+d2bqqDk9RbXlXP2bdN5+WPV1Jb38Tq2gau+c9H/OW5udTk94dzngy+jLQqHw6nPQAVo9MXtEgHZU0bVjoIzvoXjDpwQ1lhXzjk1zB+Mo1E+HBJFcf+5RUWrq4FYMnaOk7866vMXLiW+qbmNAUuIhnmNWA/M9sOwMyKzWx74ANglJmNCM87uY3rXwFOCR+fBryUulAlmRLqMQjHlJ0H7Ax80VedkR+s8oXq9Y1c/985vPXpGi78RyV/OmUC0xes5h+vfgLA/tsPZP/txpNz5r+gsQbcgyFGxQMgRwtaSebIqjas7zZw0hRoqA7mAuWXBPuK5BbQUN/Eebe/wdJ19Zz0t1e54bSJ/M997zB/RQ3n3v4Gr/zkYAq0SLWItMPdl5vZ2cDdZlYQFv/c3T8ys28BT5rZCmB6G1V8F7jVzH4MLAfOSXnQkhTm7u2fFIzL/QA4Ffg/guzvfXf/XmrDa180GvXKysp0h5GxVlTX861/zmD6gtUblX/noO0470sjKS/OT1NkIm0zsxnuHu3A+b2iDWtubuHj5TV87a+vUFXf9EV5YV4O916wDzsN7UO+hhOJpFVH26+exsxK3b06XNDhL8Acd78m3XFJciT6s+927v4LoMbdbweOAnZp7yIzG2Zmz4XLX80ys++F5ZeZ2edm9nZ4OzLmmp+a2Vwz+9DMDuvMm5LEDSgt4Naz96Awb8M/hT1HVnDBl0cpKZBs0uE2LBPbr0gkh9GDSvjTKRM2Kv9/x41j7FZlSgpEJBm+YWZvA7OAvgSrFEmWSLRTuTG8X2Nm44AlwIgErmsCfujub5pZH2CGmT0THrvG3f8Ye7KZjSUYk7YzsBXwHzPbPlxuUFJgdW0D1zzzEesbNyw9WrlgFa/PW8W+o/t/MSFZJMN1pg3LyPZreVU9v3n8/Y3Krn76I/Ye1Z9t+hV3ZygikoXC3gH1EGSpRHsMbjSzfsAvgKnAbOAP7V3k7ovd/c3wcRXwPrD1Fi45FrjH3evDpbDmAnsmGKN00Lq6Rm57eT53hHMKzt53W6Lb9qPF4YJ/VDJr8TqaW9ofaiaSATrchmVi+7WurpFTb3qd+StqKMzL4edH7URZYS6L1q7n6397lTW12udARETalmhicJu7r3b35919lLsPcvcOdR2FM9h3A14Pi74d7oh3a/iBDcGH7mcxly1kyx/EWau2oYnq9Y3tn9gFfQpzOXHiNgzsU8B3DtqOHxy6PX87Y3f2HNGPPUdWMHJACZEcrSQmWaFLbVgmtF+NTS00Nbfwg0O3pygvwr0X7MOZ+4zggYv2pawwl29O2o6cLNjjTUREUifRcSLzzexJ4F7gv57IjOUY4TbaDwLfd/d1ZvZX4DcE69r+BrgKOJf469lu9lpmdgHhhhrDhw/vSCg93traBj5aVs3NL86jvrGFU/caTnREBRUlyR/vb2YMqyjm39/bn9wco29R8Bo3nL47EMw/EMkSnW7Dkt1+hXUmrQ1rbG5h2bp6bnlpHh8sqeIbXx7Fy5ccRGlBhPzcHEYPLGHajw8kN8coK8prv0IREem1Ek0MdgC+ClxMsPzUvwi6zNtdl9bM8gg+VO9s3WXU3ZfGHL8JeCx8uhAYFnP5NsCiTet09xuBGyFY0SPB99DjrW9s5qG3PufX/5r9Rdm0j5ZzcnQYlx65E32Lk/+hbmabJQBKCCQLdaoNS0X7FdaRtDbs89V1HHndi9Q2BFMZXvl4JbtsXcYtZ+3BoLIIkUhOSn5YEBGR7JPoBmd17n6fu58ATADKgOfbuy5cyuoWgmUBr44pHxpz2vHAe+HjqcApZlZgZiOBMbS9Rm7WWVvXyJVPfbhZ+b2Vn7E2xcOKRLJZZ9qwTGi/quubuOqZD79IClrN/Hwdj7z9ueYIiUjCzOxn4Qps74Yrru0Vlk8LV1prXYntgS3U8aiZvZrAa0XN7Lpkxi/JkfCSM2Z2AMEOd0cAbwBfT+Cy/YAzgJnh0lYAlwKTzWwCQTf7AuBCAHefZWb3EUwMbAIu7k0rEtU1NG/2Ad/q01U1DK/QiiIindWJNqzHt191Dc28u3Bt3GOvz1/F6Xttq5XFRKRdZrYPcDQw0d3rzWwAENvVeJq7b3HDFTMrByYC1WY2MlyEIa6wLm1C1QMluvPxfOBt4D7gx+5ek8h1YTd9vHG3T2zhmsuByxOpP9sU5kUoyM2hvqlls2NblyspEOmszrRhmdB+FeblMGZQKZ+srN3s2K7blFOQqx3KRbLRiEsePxX4LTAc+BS4dMHvj7qrC1UOBVa4ez2Au6/oRB1fA/4FLCVYuvl3AGZ2EvAroBlY6+5fNrNJwI/c/Wgz2xO4FigC6oBz3H3z4RPSLRL91Bjv7se7+93xPlDN7KdJjqtX6luUxzcnjd6s/LCdh1CegvkFIr1IVrZhfQrz+MnhO5If2bgp71ecx8nRYUQiSgxEsk2YFNwEbEvw48W2wE1heWc9DQwzs4/M7IawhzXWnTFDia5so47JwN3hbXJM+S+Bw9x9PHBMnOs+AL7s7ruF5/62C+9DuiihHgN3X9fOKScRZobSeUX5Ec7edwQ7DSnj5pfmsb6xhVP2GMbh44bQT7sQi3RaNrdhwyqKefL7+3PNMx8xZ1k1e4yo4JuTRjOojxYREMlSvwU2HUZQHJZ3qtfA3avNbHdgf+BA4F4zu8Tdp4SnbHEokZkNBrYDXnJ3N7MmMxvn7u8BLwNTwqGWD8W5vC9wu5mNIRiiqV9C0yhZg0+1OHaSlBfnc9i4Iew1qoIWd8qL8snRXgIiqZax/5MV5kUYNbCU339tV9Y3NlNSkEthXiTdYYlI6rS1xnGX1j4O50RNA6aZ2UzgLGBKgpefDPQjWBoaggUeTgF+7u4XhROZjwLeDudoxfoN8Jy7Hx/uGTOtK+9DuiZZ/cxa+iLJyovzqSgpUFIg0j0yvg0rKcilf2mBkgKR7PdpB8vbZWY7hL/Yt5oAfNKBKiYDh7v7CHcfAexOkBhgZqPd/XV3/yWwgo2XdYagx+Dz8PHZnQhfkihZiYG+vYpIJlMbJiKZ4lJg0xUHasPyziolGM4z28zeBcYCl8Ucj51j8J/YC8Nf+YcDr7WWhSsSrQt7Cq40s5lm9h7wAvDOJq99BfA7M3sZ0C8baZasoUT3J6keEZF0UBsmIhlhwe+PumvEJY9DElclcvcZwL5tHJvUzrULgK3jlE8MH54Q57Jp4Q13fxXYPubYL9oJV1Io0eVKt7gJhbt/NznhiIgkn9owEckmYRLQleVJReJKdChRIcGmFXPC2wSC9WhnhDcRkZ5MbZiIiEg7Eh1KNAY40N0bAczsb8DT7v6DlEUmIpI8asNERETakWiPwVZAn5jnpWGZiEgmUBsmIiLSjkR7DH4PvGVmz4XPD2Dj2eoiIj2Z2jAREZF2JLrz8W1m9m9gr7DoEndfkrqwRESSR22YiIhI+xIaSmTBNnaHAOPd/VEg38z2TGlkIiJJojZMRKR9Zna8mbmZ7RhTdqWZzTKzK+Ocf4yZXdK9UUoqJTrH4AZgH4Kd7QCqgL+kJCIRkeRTGyYi0r7JwEuEuxaHLgQmuvuPY080s1x3n+ruv+/OACW1Ep1jsJe7TzSztwDcfbWZ5acwLhGRZFIbJiLZ47K+p7LJBmdctrZL+xqYWSmwH3AgMBW4zMymAiXA62b2O+AIYBWwG/Cmmc0Eou7+bTMbDPwNGBVW+U13f8XMHgGGESwb/Sd3v7ErcUpqJdpj0GhmEcABzGwg0JKyqEREkkttmIhkhyApuAnYFrDw/qawvCuOA55094+AVWY20d2PAercfYK73xuetz1wiLv/cJPrrwOed/fxBPvGzArLz3X33YEo8F0z69/FOCWFEk0MrgMeBgaZ2eUE3Uy/TVlUIiLJpTZMRLLFb4HiTcqK6XqbNhm4J3x8DxuGXm7qfndvjlN+EPBXAHdvdve1Yfl3zewd4DWCnoMxXYxTUijRVYnuNLMZwMEE2elx7v5+SiMTEUkStWEikkWGd7C8XeGv+AcB48zMgQjgZva/cU6v6UC9kwgWftjH3WvNbBrBkCLpoRJdlWg0MN/d/wK8BxxqZuUpjUxEJEnUholIFvm0g+WJOBG4w923dfcR7j4MmA98qQN1PAt8E8DMImZWBvQFVodJwY7A3l2IUbpBokOJHgSazWw74GZgJNClSS4iIt1IbZiIZItLgdpNymrD8s6aTDDcMtaDQEfmLXwPODCckDwD2Bl4Esg1s3eB3xAMJ5IeLNFViVrcvcnMTiCYUX596+oeIiIZQG2YiGSHy9bexWV9IYmrErn7pDhl14UPvxlTdvYm50wBpoSPlwLHxqn+iM7GJd0v0cSg0cwmA2cCXw3L8lITkohI0qkNE5HsESQB6vWUpEt0KNE5BJsDXe7u881sJPDP1IUlIpJUasNERETakeiqRLPN7EfAjma2C/ChdrrrAneoWQaNdZCTC3nFUFyR7qhEspbasF6ibi00VENLI+QWQvFAiETSHZWISMZIKDEws6MIdrP7mGCpv5FmdqG7/zuVwWWl+mpYOB0e/yGsmp0hD+EAACAASURBVBeUDdsLjv0LVIyCHH2IiSSb2rBeYO1CeOx/YO7TwY8vpYPh4F/BjkdBkRagEhFJRKJDia4CDnT3Se5+AMF22dekLqwstuJD+OcJG5ICgM9eh1sOheql6YtLJLupDctm1ctgylEw56kgKYCgPX30WzDnaWjRJtciIolINDFY5u5zY57PA5a1d5GZDTOz58zsfTObZWbfC8srzOwZM5sT3veLueanZjbXzD40s8M69G56utrV8MwvN3xwxapbDe/eC83xNhMUkS7qcBum9iuDLHoLVi+If+zZ/4Pa5d0ajohIpko0MZhlZk+Y2dlmdhbwL+ANMzshXP6vLU3AD919J4JNLS42s7HAJcCz7j6GYEOMSwDCY6cQrH17OHCDmWXP2JrmelhY2fbx+S9A46ZLE4tIEnSmDVP7lSnmPtv2sbWfQVND98UiIpLBEk0MCoGlwAHAJGA5UEGw7N/RbV3k7ovd/c3wcRXwPrA1wTq3t4en3Q4cFz4+FrjH3evdfT4wF9izA+8nvoYaWLcIZk8NPkCqlqTnl3kzKB3U9vHSIRDRCooiKdDhNqzHtF/SvvLhbR/LLQgWeeit3KFqKSx4CWY9DGs+hfqqdEclPYyZlZvZt9o5Z4SZtbvhWXjee8mLTlqZ2QIzG5DK10h0VaJzuvpCZjYC2A14HRjs7ovDuhebWeu35a3ZeFe8hWFZ5zXUBmNMHzwfWpqCssJyOOdxGDyuS1V3WPEA2Ptb8OQl8Y/vczHkFXZvTCK9QFfbsLS1X5KYnY+D//wSWuL84LPryVDYt/tj6ilWL4DbjoCqxcFzM/jK5bDb6b377yKbKge+BdywhXNGEOyErP0TslhCPQZmNtLMrjazh8xsaust0Rcxs1KCrbW/7+7rtnRqnLLNBuSb2QVmVmlmlcuXtzN2dP0aeOgbG5KC1rL7zgomrHWnnAjsclKwSkYsMzj0N1A+rHvjEeklutKGJbv9CutMvA2T9hVVwNf/sXnPwFa7wUE/h/zi9MSVbrWrYep3NiQFEPQgPHUp1K5KX1zSE/0eGG1mb5vZleHtPTObaWYnx5yzf3jOD8KegRfN7M3wtm8iL9TWdWY2ycyeN7P7zOwjM/u9mZ1mZtPDOEaH500xs7+FdXxkZkeH5TuH575tZu+a2Ziw/PSY8r+3DvE0s3PC6583s5vM7M9biHmKmV1nZq+Y2TwzOzEsH2pmL4R1v2dm+4flfw3b+Flm9uuYehaY2W/N7NXw+EQze8rMPjazi2L+Di+Y2cNmNjt8r5t9X4/3vsLblJj/dj9I5L9JrET7Vx8BbiEYl9uh5R3MLI/gQ/VOd38oLF5qZkPDX9uGsmES4EIg9tvxNsCiTet09xuBGwGi0WjcD94vLJsNzY2bl6+cm57x/CUD4Jjr4cCfwUdPQV4J7HhE0ItRWNb98Yj0Dp1qw1LRfkEH2zBpX34xjD4Qvj8T5v4nGDozahJUjICSgWkOLo2a6mDBi/GPffxfqDive+ORnuwSYJy7TzCzrwEXAeOBAQTzsV4Iz/mRu7d+ES8GDnX39eGX8LuBaAKvtWwL140HdgJWESwScbO772nB4g/fAb4fnjeCYGjoaOA5M9sujPlP7n6nmeUDETPbCTgZ2M/dG83sBuA0M3sG+DWwO7AWeA54q524hwJfAnYEpgIPEPSgPOXul4cJR+uvED9z91Vh2bNmtqu7vxse+8zd9zGza4ApwH4Ew11nESyrDcEw1LHAJ8CTwAnh6wHQ1vsK69ja3ceF53V4reZEE4P17n5dRys3MyP4MH7f3a+OOTQVOIsg+zwLeDSm/C4zuxrYChgDTO/o624kv09bwUFOmsbzF/cPboN3Ts/ri/Q+HW7DekT7JYnLKw5uE89MdyQ9h+UEG701rd/8WKH2dpA2fQm4292bCX4IeR7YA9i0xzQP+LOZTQCage0TrH9L173ROlTTzD4Gng7LZxIsM93qPndvAeaY2TyCL+uvAj8zs22Ah9x9jpkdTPDl/42gSaeIIDHZC5jm7svD17o3gfgfCV9ztpkNbo0XuDX8EekRd387LP+6mV1A8D17KMGX/NbEoLW3eiZQGs5hqzKz9TFf5Ke7+7wwtrsJ/pt8kRgAbb2vfwGjzOx64PGYv1/CEk0M/mRmvwpfoL61sHVi3hbsB5wBzDSz1j/WpQQfqPeZ2XnAp8BJYX2zzOw+YDbBiiAXh/8wO6/fCOgzJJhwHGu7r/Te7mWR3qczbVj62y+Rrigog93OgDdu2rg8vxS2TWjUh/RO8YZFxvMDgkUdxhMMTY+TgXb4uvqYxy0xz1vY+Dvrpj2t7u53mdnrwFHAU2Z2PsF7ud3dfxp7spkdF6eO9sTGZuGLvmBmXw5f8x9mdiXwIvAjYA93X21mUwh6BDatJ/b9bfoeN3t/mzyP+74AzGw8cBhwMfB14NyE3l0o0cRgF4IPyIPY0A3v4fM2uftLtP0P7OA2rrkcuDzBuNpXOgjOfiKYfLzozeAXlB2OgqOvgqJ+7V8vItmgw21Yj2i/RLoivxgm/QSa6uHdu4NhtQN3hBNu6t1DrCSeKqB1iMULwIVmdjvB6m1fBn5MsJhC7DCMvsBCd2+xYBnoRJdn7ux1sU4K4xsJjAI+NLNRwDx3vy58vCvBj0GPmtk17r7MzCrC9/A6wQ9G/Ql6Qk4C3uloEGa2LfC5u99kZiXAxLCeGmBt2LNwBDCtg1XvaWYjCYYSnUw49DTGs228rxqgwd0fDHtcpnT0PSWaGBwPjHL3zFsM2gz6j4bT7ofGumACcH6JVmMQ6V0ytw0T6YqSgXD47+DAS6GlEXKLoFRJgWzM3Vea2csWLDP6b4JhL+8Q/IDyv+6+xMxWAk1m9g7BF84bgAfN7CSCMfo1Cb5cZ6+L9SHwPDAYuCicr3AycLqZNQJLgP8Lx/n/HHg6nMDbSNCb+5qZXUYw/Ggx8CadS1AmAT8OX7MaONPd55vZWwTj/ecBL3ei3lcJeqd3IUjUHo496O6z470voA64LWay8mY9Cu0xj7cL76YnBWOvvuPu3byMT/ui0ahXVm5h0zARyTpmNsPdE5nk1nq+2jAR6RE62n7JxsKhOY+5+wPtnduBOs8Gou7+7WTV2YVYJhEzybu7JdpjMBj4wMzeYOPxucekJCoRkeRSGyYiItKORBODX6U0ChGR1FIbJiLSjczsMOAPmxTPd/fju1Kvu5/dlevbqHMKMMXMfka4oESM+8P5Y93C3afR8TkJSZPozsfPpzoQEZFUURsmItK93P0p4Kl0x9ERWkAi8Z2P9zazN8ys2swazKzZzLa0A6iISI+hNkxERKR9CSUGwJ+BycAcgk0Uzg/LREQygdowERGRdiQ6xwB3n2tmkXDDntvM7JUUxiUiklRqw0RERLYs0cSg1szygbfN7AqCNV9LUheWiEhSqQ0TEdkCMxsG3AEMIdgI8kZ3/1N6o5LuluhQojPCc79NsBHFMOBrqQpKRCTJ1IaJiGxZE/BDd98J2Bu42MzGpjkm6WaJrkr0SfhwPfDrTY+b2YPurg9ZEemR1IaJSLaJRqO5wABgRWVlZVNX63P3xQS9qbh7lZm9D2wNzO5q3ZI5Eu0xaM+oJNUjIpIOasNEJGNEo9F9geXAfGB5+DxpzGwEsBvwejLrlZ4vWYmBJ6keEZF0UBsmIhkh7Cl4HCgHCsP7x6PRaCQZ9ZtZKfAg8H1317LOvUyyEgMRERERSb0BBAlBrEJgYFcrNrM8gqTgTnd/qKv1SeZJVmJgSapHRCQd1IaJSKZYQTBfKtZ6gqFFnWZmBtwCvO/uV3elLslcyUoMfpKkekRE0kFtmIhkhHCi8VHAGoKEYA1wVGVlZXMXq96PYAW3g8zs7fB2ZBfrlAyzxVWJzGwm8cfeGuDuvivBg6dTEJuISJeoDRORbFRZWflKNBodQDB8aHkSkgLc/SXUe9rrtbdc6dHdEoWISGqoDRORrBQmA0vSHYdkly0mBjFrf4uIZBy1YSIiIolLaI6Bme1tZm+YWbWZNZhZs5lpCSsRyQhqw0RERNqX6OTjPwOTgTlAEXA+cH2qghIRSTK1YSIiIu1ob47BF9x9rplF3L0ZuM3MXklhXCIiSaU2TEREZMsSTQxqzSwfeNvMrgAWAyWpC0tEJKnUhomIiLQj0aFEZ4TnfhuoAYYBJ6QqKBGRJFMbJiKyBWZWaGbTzewdM5tlZr9Od0zS/RJNDI5z9/Xuvs7df+3u/4OWARSRzKE2TERky+qBg9x9PDABONzM9k5zTNLNEk0MzopTdnYS4xARSSW1YSKSNaLR6KnRaHRBNBptCe9P7WqdHqgOn+aFt3gbREoW22JiYGaTzexfwEgzmxpzmwasbK9yM7vVzJaZ2XsxZZeZ2efxtts2s5+a2Vwz+9DMDuvC+xIRURsmIlknTAJuArYl2Kl4W+CmZCQHZhYxs7eBZcAz7v56V+uUzNLe5ONXCCbpDQCuiimvAt5NoP4pBMsE3rFJ+TXu/sfYAjMbC5wC7AxsBfzHzLYPVxAREekMtWEikm1+CxRvUlYclt/VlYrD9mqCmZUDD5vZOHd/r73rJHtsscfA3T9x92nuvg/wAdAnvC1096b2Knf3F4BVCcZyLHCPu9e7+3xgLrBngteKiGxGbZiIZKHhHSzvMHdfA0wDDk9WnZIZEt35+CRgOnAS8HXgdTM7sQuv+20zezfspu8Xlm0NfBZzzsKwLF48F5hZpZlVLl++vAthiEhvoDZMRLLIpx0sT4iZDQx7CjCzIuAQgh9UpBdJdPLxz4E93P0sdz+T4FewX3TyNf8KjCaY8b6YDd37FufcuJNe3P1Gd4+6e3TgwIGdDENEehG1YSKSLS4Fajcpqw3Lu2Io8JyZvQu8QTDH4LEu1ikZJtENznLcfVnM85UknlRsxN2Xtj42s5uA1n90CwnWFm+1DbCoM6/Rbeqrgtv6dVBYBgVlUFCa7qhEZHOZ3YbVrABvhuYmiORDUTlE8pJStYhklsrKyrui0SgEcwqGE/QUXFpZWdnV+QXvArt1PULJZIkmBv82s6eAu8PnJwNPdOYFzWyouy8Onx4PtE5qmQrcZWZXE0zcG0PQ9Z9ajbVQswKf/SjU12DjjofSIVDUF4AV1fVM+3AZHy2p5tgJWzGsopiyojyoXQUvXQOv3QAtTcGH9L7fg30uhuKKlIctIh2SuW1Y1VKoXgpznoa61TD6IOi7DS3lI2kgQmFeZLNLVtU0MHvROp77YBn7jenPhGHlVJQUdDkUEekZwiSgS4mASDyJJgYO/B34EkF3+Y1Au5temNndwCRggJktBH4FTDKzCWGdC4ALAdx9lpndB8wGmoCLU76aR3MTfPoa3HkS1hLOQ3z+d3DYb2G3M1nZlM85t73BzM/XAnDji/P42+kTOWznIdiCl+CV62LqaoQX/wjD94Yxh6Y0bBHpsMxsw+qrYMVHcNdJ0FgXlL36Z1r2vJA5Y7/N1S8u53sHj2GHwaVEIkEHSG19Eze/OI8bpn0MwC0vz2fynsP46ZE7UVaoXgYREWmbube/d4WZvenuEzcpe9fdd01ZZAmKRqNeWVnZuYurlsLfvwTVyzYuz4nA92byYV0Zh137wkaHth9cyl3n7s6Ah0+GBS9tXufog+GkKcHQIhFJCTOb4e7RDpyfmW3Y2oVw54mw7P3NDi0762UOnvI5lgPP/s8BDOxTCMDSdev58hXPUd/U8sW5OQYvX3IQQ/sWpeQ9iEjiOtp+iXSn9jY4+6aZzQR2CFfgaL3NJ7E1wHu2xprNkwKAlmZY/iGROH+d5tbPWm/Z/OAX5dooUKQnyPg2rKUpblIAkLfwNbYbVEpLy+YtTssmP/ioRRIRkUS0N/nuLuCrBGNnvxpz293dT09xbKkXKQh6B+IpHUS/4nx2GNxno+LvHrwdFX1KIHp+/Ov2OB8K+yY5UBHppMxuw3LyID/+ggb5A7Zlu0Gl3Hn+XvQrzv+ivLQgl9P32najc48dvxXF+YmOHBURkd5qi58U7r4WWAtM7p5wullBH9jlJHjnno3LB+8MpYPpX1rAP8/fi8feXcT7i9dx8h7D2G5gKTk5BqMmwe7nwpu3gTtYDuxxAWy7TzreiYjEkfFtWHF/2OuiYP5SrLKtKdp6F341sj8lBbmYbVgptaQgl+8cPIa9R/fn6VlLOGD7Qey3XX/6Fml+gYiIbFlCcwx6si7NMQCoWQ6v3gCVt0DTehh7PBxyGZQN3ei0lhYPEoJYdWuhoSoYjlQ6OFiqVL0FIimXTWN0223DalfCm3fAq3+GujWw/eHBAgn9tm37mlDcdktE0qqnt19mFgEqgc/d/eh0xyPdS33LJQNh0iXBr3IQfLnPL9nstLgfrkV9g1vfbVIcpIj0WsX9YZ/vwIRTg8kCeUUJL26gpEAkO0Wj0e2Bi4AdCXYn/ltlZeVHSar+e8D7gFZR6YU6tcFP1sktgD6Dg1ucpEBEJK0iuUGvZJ/BWvFMpJeLRqPnA28D3waOCO/fDsu7xMy2AY4Cbu5qXZKZlBiIiIiIZICwp+A6oAhonTiUFz6/LhqNjuniS1wL/C/QxtKLku2UGIiIiIhkhotoexh4bni8U8zsaGCZu8/obB2S+ZQYiIiIiGSGHdnQU7CpPGCnLtS9H3CMmS0A7gEOMrN/dqE+yUBKDEREREQywwdAYxvHGgkmDXeKu//U3bdx9xHAKcB/M2K/F0kqrUokkkIrqupZsm49K2vq2X5wH/oU5lFaoP/tpGOWV9WzeG0da+saGTOolL5F+RTlt7E5o4hks78RDBeK12vQFB4X6TR9QxFJkcVr6jj9ltf5eHkNAJEc47KvjuW43bamT6E2m5LELFxVyyk3vcbC1XUA5EWMK762K1/ZeQglSjJFepXKysqPotHodwkmIOcSJAiNBEnBdysrK+ck43XcfRowLRl1SWbRUCKRFKha38hl/5r1RVIA0Nzi/OLRWaysaUhjZJJJ1tQ28OMH3v0iKQBobHZ+eP87rKlrazSBiGSzysrKm4HxwPXAv8P78WG5SJfo5yaRFKhtaOaZ2UvjHnt29lLO239UN0ckmWh9YzOvzlu5WXmLw2sfr+Rru2tzRZHeKOwZ+GG645Dsox4DkRQwID83/v9exRr+IQkyjNw2di/WHAMREUk2JQYiKVBWlMfJ0eGbledHcpi0/cA0RCSZqKQwwtG7Dt2svDg/wu7b9ktDRCIiks2UGIikQGFehO8evB3HTtiK1h98h/Yt5J/n70lFaX56g5OMUVqQx8+PGstXxg7Gwn9HwyqKuPsbe1NRognsIiKSXBrTIJIi/UsL+M1x47jkiB1paGqhOD9CRUkBkTaGhojEM6BPAX88aTzV9U00NbdQlJ/LgNJ8zPTvSEREkkuJgUgKlRXmUaalSaWLyoryKCvSvyMR2SAajY4EtgIWVVZWzk9GneGux1VAM9Dk7tFk1CuZQ4mBiIiISIaIRqNR4O/ATkADkB+NRt8HLqysrKxMwksc6O4rklCPZCDNMRARERHJAGFSMA2YCBQBfcP7icC08LhIpykxEBEREckMfwdK2jhWAvyti/U78LSZzTCzC7pYl2QgJQYiIiIiPVw4p2Cndk4bG57XWfu5+0TgCOBiM/tyF+qSDKTEQERERKTn24pgTsGWNITndYq7LwrvlwEPA3t2ti7JTEoMJOOsrWtk3vJqPliyjjW17bWRItmnvrGZxWvqePuzNayork93OCLSPRYB7W2Ekx+e12FmVmJmfVofA18B3utMXZK5UpoYmNmtZrbMzN6LKasws2fMbE543y/m2E/NbK6ZfWhmh6UyNslcL85ZzkFXPc/h177I3dM/pb6pOd0hSZbqqW3YknXrmfTHaRz3l5e5+M43WVWjBFkk24VLkr7fzmmzu7B06WDgJTN7B5gOPO7uT3ayLslQqe4xmAIcvknZJcCz7j4GeDZ8jpmNBU4Bdg6vucHMIimOD4DlVfW8PHcFT89awrJ162lp8e54WemExuYWnp619Ivn//1gGXUNSgwkZabQA9uwD5ZUUd/UAsD0Batoam7Z6Pi6sFft0bc+Z87SKtbWNqYiDBHpfhcCNW0cqwEu6mzF7j7P3ceHt53d/fLO1iWZK6WJgbu/AKzapPhY4Pbw8e3AcTHl97h7vbvPB+bSDWPbllfVM/mm1zjt5te54B8zOOzaF1i8bn2qX1Y6KS+Sw3lfGklhXg6RHOOiA0ZTWqDtOCQ1emobNmFYOVuXFwFwzr4jKMjb0JQ3NLXwxMzFHHTV83zv3rc59JoXuGv6J9TWN6UiFBHpRuE+BZOAGUAdsDa8nwFMStI+BtKLpeMb1WB3Xwzg7ovNbFBYvjXwWsx5C8OylJo+fxVzl1V/8Xx1bSM3vziPS4/cibyIpmD0RDsN7cPzPz4Qd6dvUT65+u8k3SvtbdjgskIeuXhfGpud4vwIfYs2DDteXdvA5Y9vPNrg2v/M4fjdtqFYSbRIxgu//EdTsfOxSE/6lLA4ZXHH9IRr614AMHz48C696LKqzXsHlqxdT3OLk9ctA5mko/JzIwwu038c6XG6tQ0b2Kcwbrk7VDds3DtQ39RCs4ZIimSVMBlQQiBJlY6fWpea2VCA8H5ZWL4QGBZz3ja0MbPe3W9096i7RwcOHNilYA4dO5j8TX5xDoaq6IuniMTVo9qwTRXnRzhs7JCNyvYeVUFRvnrWRERky9LxSTEVOCt8fBbwaEz5KWZWYGYjgTEEs+JTakBpAVO/vR+Hjh3EvqP784/z9mT7IX1S/bIikrl6VBu2qbKiPP7f8eP4/iFjmDi8nG8eMIrrJ+9GRUlBd4ciIiIZJqVDiczsboJJMgPMbCHwK+D3wH1mdh7wKXASgLvPMrP7gNlAE3Cxu6d8uZnCvAg7Di3jmpMn0NwCfYvyUv2SkibrG5tZW9dIfWMzBXkR+hXnk5+rX1GlbT2lDVtX10hNQxONzU5RXoQBpfmYxRu5FBhQWsDFB27HWfuOoDQ/lzz9OxcRkQSYe2aPO41Go65J+NKeNbUN3D39U67/71xqG5opK8rlZ0fuxBHjhlKmZDDjmNkMd4+mO45kaK8NW1FVzy+mvsdT7y2hxWFE/2KuP3UiOw7powUSRDJQNrRfZnYZUO3uf0x3LJJc+lSRXmH6/FX84ckPqQ33PFhX18RPHpzJvBVtLQctkn419U389on3+ffMICkAWLCyllP+/qo2NRMRkaRTYiBZb01tAze+MC/usZte+FgbpEmPVV3fxNR3Np+/XNPQzJufrk5DRCKSrcysxMweN7N3zOw9MzvZzBaY2YDweNTMpsVcMt7M/hvuAv+N8JwcM7vBzGaZ2WNm9oSZnRge293MnjezGWb2VMwiDtPM7A9mNt3MPjKz/cPys83szzHxPWZmk8wsYmZTwhhnmtkPuutv1Bv0pOVKRVKiucWpa4z/5b+mvplmbwG0CpX0PC3uNLWxzGjNem1YJiJJdTiwyN2PAjCzvsAftnD+rsDeQAnwlpk9DuwLjAB2AQYB7wO3mlkecD1wrLsvN7OTgcuBc8O6ct19TzM7kmAu1yFbeN0JwNbuPi6Ms7wzb1biU48BQM0KWPMprJoP1UuhpSXdEUkSlRfnc+Lu28Q9NnmvYZQWaI6B9ExFeRH2GlmxWbkZ7DOqH6yaB1VLoEnDikSky2YCh4S/3u/v7mvbOf9Rd69z9xXAcwQ7vX8JuN/dW9x9SVgOsAMwDnjGzN4Gfk6wpHOrh8L7GQSJxZbMA0aZ2fVmdjiwLsH3Jwno3YlBcyMsnQX/PAGu3QWumwA3HwrznoP66vavl4wQyTGOGb8V+48ZsFH5cRO2Yo8R/dMUlUj7yovzufKk8WzVd8NmZrk5xhXHbEff6VfBdbvBDXvDGzdB7ao0Rioimc7dPwJ2J0gQfmdmvyRYYa31u+Kmuypu2p3pxN/okbB8lrtPCG+7uPtXYo7Xh/fNbBjNEvvaX7y+u68GxgPTgIuBm9t/d5Ko3j2U6P+3d+dRUpTnHse/P4cZZpgBEUREFHHBKG6gE664i2sSEzHuu8YTE5NczU00Md6TuN1ovN7E5ZrEJa5Ro0Zj9KpxQ3GJRBwVAZdEUVQCiojsyDbP/aNqQjv0DAMz09U9/fucM6er3qqufii63+qn36UWfAQ3HABLc5KAOe/B7YfBN8fCRsMyC806Vt+67lx59HA+WbCE92YvYot+dazXo5LePaqyDs2sVYP69OD+7+3OtE8XMWfhErbptYRer1xLj4ZfJzss/hQePRfW6Qb1p0JFeVfrZrZ2JG0EzI6I2yQtAE4GppIkC38BDmv2lEMkXULSlWhv4BygCjhJ0i1Av7T8DuDvQD9JIyNiXNq1aKuIeK2VkKYC35G0DjCQpEWCdMzD0oi4V9IU4Ob2/cstV/leQVYsg5du+nxS0CQCnjgfjrgJatYreGjWOfrUVtGntooh/X0DOyst/Xp2p1/P7vDBeLjhIGjMM2bm6Uth6CHQc8NVt5mZrd72wGWSGoFlwOlADXCDpHOBF5rtPx54CBgEXBQR0yXdC+wLTAb+kT5nbkQsTQchX5WOXegGXAG0lhj8FXiXpAVjMvByWj4QuClNGAB+0o5/szVTvonB0gXw3vMtb5/+Mixf0vJ2M7NCmzImf1IAsOgTWP5ZYeMxsy4jIh4FHs2zaas8+57fwjEaJZ0VEQsk9SVJHial2yYAe+Z5zt45y7NIxxhEcqOt41oId6dW/inWDuWbGHTrDj03anl7XX+QZ6oxsyLSe1DL29bpBhXuGmdmmXswnSmoiqQl4cOsA7K2K9/EoLIH7PrvMPme/Nt3PQNq18+/zcwsC5vvA92q87cMDB0N3d1NzsyyldsCYKWnvGclWm8w7H9RMvdfrh2Ohq2/vGq5mVmWavrAcX+EyprPlw/YEQ78uRMDMzNrl/JtMQCo6Q07nwzbjoa3xyS/wm25L/RYH+F0jQAAEzRJREFUH3qsOne4mVmmKqth43+DM16Bqc/D/BkwaBfovSnU9cs6OjMzK3HlnRgAVPdK/upPyToSM7PVq+wOlQNg++YzB5qZmbVPeXclMjMzMzMzwIlBorExudnZ/A893Z+ZFacFM5M6aunCrCMxsy5K0oaS7pQ0RdLrkh6WtMp0pTn7L0gfB0uanC7XS7qqUDFbx3JiADB7CvxmJFy5I7z/N1ixPOuIzMxWmvMBXD8KLh8Kbz4MSxdnHZGZdTGSBNwHjI2ILSJiKHAu0H9NjhMRDRFxRmfEaJ3PiUEEPP+/K28O9MQF8NncrKMyM1tp8r0w94Pk5mZP/AyWuI4yK2f19fXH1tfXT62vr29MH4/tgMPuAyyLiGuaCiJiQkQ8K+lsSS9KmijpgtYOImlvSQ+myyMkPS/plfTxCx0Qp3UiJwYSDN5j5fomI1adCtDMLEubjFi5vNFOyb0MzKwspUnA9cCmgNLH6zsgOdgOeKl5oaQDgCHACGAYsLOkVe5g3II3gT0jYjjwM+DidsZoncyzEi2eAxsNgxPvT5b7bwvLFkNVj6wjMzODJfOhtj9849Gk1WDA8KTlwMzK1cVA8y8pPdLyOzrh9Q5I/15J1+tIEoVn2vDcdYFbJA0BAqjshPisA5V3YrDss+TOxw/9EHpumPwK9+lUGHY8HPhfULNe1hGaWTlrbIR3n4W7joUefaF7r6SO2mwvOOx3vju7WXkatIblbfUacHiecgGXRMS1a3HMi4CnIuJQSYOBsWsdnRVEeXcl+mwOPHF+sjz/w+SCCzDhtqT1wMwsSws/hofPSsZCLZwFs9+BaIR3noK507KOzsyy8f4alrfVk0B3Sd9sKpD0RWAe8A1JdWnZQEkbtPGY6wL/TJdPbmd8VgDlnRgsW5Q00+cz+53CxmJm1lzjMpj3z/zbZkwobCxmVizOBRY1K1uUlq+1iAjgUGD/dLrS14DzSbon3QGMkzQJuAfo2cbD/jdwiaS/AhXtic8Ko7y7EnWrSQYaL8sz9V/v9rbImZm10zrdoLZf0nLQ3AZDCx+PmWWuoaHhjvr6ekjGFAwiaSk4t6Ghod3jCyJiOnBknk1Xpn/N969LH6eSDF4mIsaSdhmKiHFA7n0QftreGK1zlXdiUL0u7HE2VFTCpiOTJvpPpsA7Tyf9ec3MstRjfdjvArj/O58v3/c86LVR0rJZVQu1GyQzrJlZWUiTgM4YaGxlrrwTg6oeMOxYePJCeDxNYgeNhMNvgh59so3NzKyiG2z9Fej9f/DMZbD40yQpmPMeXDUcVixNWjdPfAD6bJZ1tGZmVuLKe4wBJHc9npCTdL8/Dibe5ekAzaw41PSGzfaEI38PJ/wZNtgGHj47SQoA5rwPD/0gSRrMzMzaIbPEQNJUSZMkTZDUkJb1kfS4pLfSx86fL3TaKvfygA9egOVLOv2lzax0FbwOq+mdTE+6YGbS7THXR5NdZ5mZWbtl3WKwT0QMi4j6dP0cYExEDAHGpOuda8t9Vy0bOto3ODOztih8HdZrwKp3Z99yf+je1klCzMzM8ss6MWjuEOCWdPkWYHSnv+K6G8NXr0rGFFTWwK5nwpD9Ov1lzaxL6vw6rLoPHH8f9NkctE4yBmHf85JByGZmZu2Q5eDjAB6TFMC1EXEd0D8iZgBExIw1uIHG2qvpnQxA/sJBSUTde7q1wMzaIps6rLI7DNoFvvFocuOzyupkhjUzs3aSdCjwJ2CbiHizg499M/BgRNyzBs9Z0DQlqhVGlonBbhExPb1wPi6pzW9ASacBpwEMGtQB9xuoqIS6/u0/jpmVk+zqMAnqOv93EzMrO8cAzwFHk9zczMpMZl2J0ptoEBEzgfuAEcBHkgYApI8zW3judRFRHxH1/fr1a3cscxcv5Y0Z83hx6mxmL1za7uOZWdeXZR22eOly3p+9kOfe+phZ8z3o2Kyc1NfXV9XX1/8l/avLWa5qz3El1QG7AaeSJAZI2lvSWEn3SHpT0u1SctMUSftKeiWdhOFGSd3T8qmSLpU0Pv3bMudl9pT0vKR3JB2e89pnS3pR0kRJF7QQ3yr7SKqV9JCkVyVNlnRUe86BZZQYpP+RPZuWgQOAycADwEnpbicB9xcinufemsWXrnyWI64ZxzVjp7Bo6fJCvKyZlais67AP5y1h1P88zfE3jOfUW1/kkwVODszKyP3AXunftJzl9tY3o4FHIuIfwGxJO6Xlw4HvA0OBzYHdJFUDNwNHRcT2JD1QTs851ryIGAFcDVyRUz4A2B04GPgFgKQDgCEkP64MA3aWtGduYK3scxAwPSJ2jIjtgEfaeQ7KXlYtBv2B5yS9CowHHoqIR0jeJPtLegvYP13vVBHBs2/N+tf63979hM+WNbbyDDOzbOuwtz6az/LGAGDitLmsSJfNrKzUAOumjx3hGODOdPnOdB1gfERMi4hGYAIwGPgC8G6aREAy2ULul/k/5DyOzCn/c0Q0RsTrJPUoJD+sHAC8ArwMbE2SBORqaZ9JwH5pC8UeETF3bf7htlImYwwi4h1gxzzlnwB55g/tPJI4dffN+MvkD1m0dDk/2H8relaX9w2hzax1WddhO27Smy361TLl44Wcue8QqqsqOvslzax4HEHSUpCbECwFDs+/++pJ6guMArZLJ1SoIJlg4WEgt0lyBcl3R63mkNHCcu6xlPN4SURc21qILe0jaWfgy8Alkh6LiAtXE5u1wt+Agc3Xr+XxH+xJBPSq6UZlRbHN4mpmtlL/XtXcedpIGhuDmqoKelVXZh2SmRXOH4Hm4wmqgHuAL63lMQ8Hbo2IbzUVSHqapNtPPm8CgyVtGRFvAycAT+dsP4qkxfQoYNxqXvtR4CJJt0fEAkkDgWXp+K1W9yH5Hjs7Im6TtAA4ua3/YMvPiQFQUbEOG/SszjoMM7M269eze9YhmFm2FpO0FLRr0HHqGFbt+ngvybiBKc13jojPJJ0C/FFSN+BF4JqcXbpLeoGky/oxzZ/f7FiPSdoGGJeOa14AHE/O5A2t7LMlcJmkRpJE4XSsXRRR2n1T6+vro6GhIeswzKyAJL2Uc7fhkuY6zKy8tLf+SmcfahpofARJCwLAIQ0NDZlPrShpKlAfEbNWt68Vn/JpMVg4Cz4YDzXrwQZbJ49mZqXgs7nwyRSYOw0GjYS69k/TbGalKf3yn9tlaG27D5mtojwSgyXzYcwF8PKtyfrXfwc7HJFtTGZmbTXjVbjlq8nykAPh69f6xw0zK0oRMTjrGGztlcco2+VLYPqElesfjM8uFjOzNZVbf300KanTzMzMOlh5JAbVveHAi6GqFnoNhJHfyToiM7O22+4w6LsFdOsOB12a1GlmZmYdrDy6ElV0g01GwPcaAEHPDbOOyMys7dYdCKc8AhFQ3RMq88yitnwJLFsMVXVJnWdmZraGyuPqsXAWvPtMMsagsgfsdgb02xpq/KubmZWIug3yl69YBvOmw7irYebrMHh32Olk6DWgoOGZmVnp6/qJwaLZ8PCP4LV7V5b9/SH4yuUw/Likad7MrFTNeR+u2R2WLUrWpz6X/AjyzafcOmpmZmuk648xWDz780lBkzHnw+JPCx6OmVmHWTIfnrxoZVLQZN50mHg3NK7IJi4zMytJXT8xmPV2/vLP5ib9cc3MStWyRclUpvm891fPXmRmZmuk6ycGfTbPX15Vm38An5lZqehWA/22yb9t4M5Q4a6SZmbWdl0/MejRF4YcsGr5Xj+Gat8gyMxKWHUv2O+8VcdK1a4Pw0+Aiops4jIzs5LU9Qcf1/aF0b+F1++HCbdBZS3sdiYMrHeLgZmVvt6bwunPw9hL4eM3YdNdYdczoK5/1pGZmVmJ6fqJASS/nu18CgwdDeusAzVuKTCzLqKyGvpuCQdfAcvT+xj4Rw8zM1sL5ZEYQJIQ1PbNOgozs87RvTb5MzMzW0tdf4yBmZmZmZmtlhMDMzMzMzNzYmBmZmZmZk4MzMzMzMwMJwZmZmZmZoYTAzMzMzMzw4mBmZmZmZkBioisY2gXSR8D763BU9YHZnVSOGvD8bTO8bSuXOPZNCL6FeB1Ot0a1mHF9v+dj2PsOKUQp2Ncc12m/rKup+QTgzUlqSEi6rOOo4njaZ3jaZ3jKS+lcH4dY8cphTgdo1nX4q5EZmZmZmbmxMDMzMzMzMozMbgu6wCacTytczytczzlpRTOr2PsOKUQp2M060LKboyBmZmZmZmtqhxbDMzMzMzMrJmySQwkHSTp75LelnRO1vEASJoqaZKkCZIaMnj9GyXNlDQ5p6yPpMclvZU+rpdxPOdL+md6jiZI+nKBYtlE0lOS3pD0mqQz0/JMzk8r8WR1fqoljZf0ahrPBWl5Zu+frq5I67CiqkNaiLGoPsstxFgynydJFZJekfRgMcaY77pabDGaFbOySAwkVQC/Br4EDAWOkTQ026j+ZZ+IGJbRVGo3Awc1KzsHGBMRQ4Ax6XqW8QBcnp6jYRHxcIFiWQ78MCK2AXYBvpu+Z7I6Py3FA9mcnyXAqIjYERgGHCRpF7J9/3RZRVyH3Uxx1SH5FNtnOZ9S+jydCbyRs16MMTa/rhZjjGZFqSwSA2AE8HZEvBMRS4E7gUMyjilzEfEMMLtZ8SHALenyLcDojOPJRETMiIiX0+X5JBfCgWR0flqJJxORWJCuVqZ/QYbvny6uKOuwYqtD8im2z3I+pfJ5krQx8BXgdznFRRVjC0ohRrOiUC6JwUDgg5z1aWT4pSpHAI9JeknSaVkHk+ofETMguaACG2QcD8D3JE1Muy0UvAlY0mBgOPACRXB+msUDGZ2ftEvBBGAm8HhEFMX56aKKtQ7Lp2jfA8X2Wc5VIp+nK4AfAY05ZcUWY77rarHFaFa0yiUxUJ6yYpiOabeI2Imke8B3Je2ZdUBF6LfAFiTN6zOAXxbyxSXVAfcC34+IeYV87TbGk9n5iYgVETEM2BgYIWm7Qr12GSrWOqxkFNtnubli/zxJOhiYGREvZR3Lavi6atYO5ZIYTAM2yVnfGJieUSz/EhHT08eZwH0k3QWy9pGkAQDp48wsg4mIj9ILZiNwPQU8R5IqSb5I3B4Rf0qLMzs/+eLJ8vw0iYg5wFiSvuZF9f7pQoqyDmtB0b0Hiu2z3Joi/jztBnxN0lSSrmyjJN1GccXY0nW1qGI0K2blkhi8CAyRtJmkKuBo4IEsA5JUK6ln0zJwADC59WcVxAPASenyScD9GcbSVIk3OZQCnSNJAm4A3oiIX+VsyuT8tBRPhuenn6Te6XINsB/wJkX2/ulCiq4Oa0VRvQeK7bOcTyl8niLiJxGxcUQMJnn/PRkRx1NEMbZyXS2aGM2KXdnc4EzJNI5XABXAjRHx84zj2Zzk1wyAbsAdhY5J0h+AvYH1gY+A84A/A3cDg4D3gSMioiADgluIZ2+SbjIBTAW+1dRXtJNj2R14FpjEyv6055L0TS74+WklnmPI5vzsQDKIr4LkB4a7I+JCSX3J6P3T1RVbHQbFV4e0EGNRfZZbiLGkPk+S9gbOioiDiynGlq6rxRSjWbErm8TAzMzMzMxaVi5diczMzMzMrBVODMzMzMzMzImBmZmZmZk5MTAzMzMzM5wYmJmZmZkZTgzMzMzMzAwnBlYCJH1b0olZx2FmJun59HGwpGNXs+8mkp6S9Iak1ySdWZgozczWju9jYGVNUreIWJ51HGZWWnJv8tXKPgOAARHxcnpH3peA0RHxeifE47rMzNrNLQaWGUknSpoo6VVJv5e0qaQxadkYSYPS/c6XdFYLx9hG0vic9cGSJqbLP5P0oqTJkq6TpLR8rKSLJT0N+Bc8M2szSQvSxV8Ae0iaIOk/JG0raXy6PlHSkIiYEREvA0TEfOANYGALx3VdZmaZc2JgmZC0LfCfwKiI2JHkonY1cGtE7ADcDly1uuNExBtAlaTN06KjgLvT5asj4osRsR1QA+T+stc7IvaKiF92zL/IzMrMOcCzETEsIi4Hvg1cGRHDgHpgWu7OkgYDw4EX8h3MdZmZFQMnBpaVUcA9ETELICJmAyOBO9Ltvwd2b+Ox7gaOTJePAu5Kl/eR9IKkSenrbZvznLswM+s444BzJf0Y2DQiFjdtkFQH3At8PyLmtXIM12VmliknBpYVAasb4NLWATB3AUdK2gqIiHhLUjXwG+DwiNgeuB6oznnOwjUN2MysJRFxB/A1YDHwqKRRAJIqSZKC2yPiT6s5jOsyM8uUEwPLyhiSC2BfAEl9gOeBo9PtxwHPteVAETEFWAH8lJW/njVdOGelv9Yd3kFxm5kBzAd6Nq2kXYDeiYirgAeAHdKxADcAb0TEr1Z3QNdlZpa1blkHYOUpIl6T9HPgaUkrgFeAM4AbJZ0NfAycsgaHvAu4DNgsPf4cSdcDk4CpwIsdGL6Z2URguaRXgZtJvsAfL2kZ8CFwIbAbcAIwSdKE9HnnRsTDrRzXdZmZZcbTlZqZmZmZmbsSmZmZmZmZuxJZCZH0a5Km+VxXRsRNWcRjZrY2XJeZWbFyVyIzMzMzM3NXIjMzMzMzc2JgZmZmZmY4MTAzMzMzM5wYmJmZmZkZTgzMzMzMzAz4f83/MYZE2r1uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1,2,figsize=(10,5))\n",
    "for ax, marker in zip(axs, ('coi','its2')):\n",
    "    # exclude compromised species labels and species with insufficient sample size\n",
    "    d = comb_div[~comb_div.index.isin(BAD_SPP) & # harder filter -  MULTISP_LABELS\n",
    "                 (comb_div.total_ampseq_nsamples >= MIN_SAMPLES_DIV) &\n",
    "                 (comb_div[marker+'_nseq'] >= MIN_SAMPLES_DIV)]\n",
    "    print(d[d[marker+'_var'] > 12].index)\n",
    "    sns.scatterplot(x=marker+'_var',y='total_ampseq_var',\n",
    "                    hue='region',\n",
    "                    size='total_ampseq_nsamples',\n",
    "                    style='subgenus',\n",
    "                    data=d,\n",
    "                    ax=ax)\n",
    "    ax.legend().set_visible(False)\n",
    "ax.legend(bbox_to_anchor=(1, 0.7), frameon=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU1fn48c+THZKwJiwSdgFDQogQ2QUUEGoVQRBFK6TaYlUatd8frbU/rVqt/NRqoWotLtWiIgqoiNtXEYpYK0KACMiiECQBIawhkD3P748ZQmIWMkkmsz3v1yuvmTn33nOfuUmeuXPuueeIqmKMMSawBHk6AGOMMU3Pkr8xxgQgS/7GGBOALPkbY0wAsuRvjDEBKMTTAdRVTEyMduvWzdNhGGOMT9mwYcNhVY39cbnPJP9u3bqxfv16T4dhjDE+RUT2VlduzT7GGBOALPkbY0wAsuRvjDEByGfa/I3xNcXFxWRlZVFQUODpUEwAiIiIIC4ujtDQ0Dqt79/Jv+g0FOaClkFwGETGeDoiE0CysrKIjo6mW7duiIinwzF+TFU5cuQIWVlZdO/evU7buLXZR0QiRGSdiGwWka0i8oCz/H4RyRaRTc6fyxt953k58PF9MD8ZnoiH166BAxlQUtjouzKmOgUFBbRt29YSv3E7EaFt27Yufct095l/IXCpquaJSCiwVkQ+cC57UlUfd8teTx+Dd26HXR+dLctOhxfGweyvoFUXt+zWmB+zxG+aiqt/a24981eHPOfLUOeP+8eQLjpZOfGfUVIA656D0mK3h2CMMd7M7b19RCRYRDYBh4CPVfVL56LZIpIhIi+KSOsatp0lIutFZH1OTk7dd5pXy7oHtzo+BIwx57R8+XLmzp3r6TCMG7g9+atqqaomA3HAIBFJBP4O9ASSgQPAX2rYdoGqpqhqSmxslbuTaxbdHmr6CtR5EIQ2c+k9GOMPVJWysjKXtpk4cSJ33323myIyntRk/fxV9TiwGpigqgedHwplwHPAoEbdWVgU9L+hanlEKxgwA4L8u5OTMWdkZmYSHx/PbbfdxoABA1i4cCFDhw5lwIABXHPNNeTlOVpl33//fS644AJGjBhBWloaV1xxBQAvvfQSs2fPBmDv3r2MGTOGpKQkxowZw/fffw9AamoqaWlpDBs2jB49erBkyRLPvFnjEnf39okVkVbO582AscB2EelYYbXJwJZG3XGzVjDuARj3J2jRCUKbQ/xEmLUKoto36q6M8XY7duxgxowZfPzxx7zwwgt88sknpKenk5KSwhNPPEFBQQG33HILH3zwAWvXrqWmJtbZs2czY8YMMjIyuOGGG0hLSytfduDAAdauXcuKFSvsm4KPcPcpcEfgZREJxvFB84aqrhCRhSKSjOPibyZwS6PvOTIGhtwG/a8FVQiLhPDoRt+NMd6ua9euDBkyhBUrVrBt2zaGDx8OQFFREUOHDmX79u306NGjvH/49OnTWbBgQZV6vvjiC5YtWwbAjTfeyG9/+9vyZZMmTSIoKIi+ffty8ODBJnhXpqHcmvxVNQO4sJryG92533LBIXambwJeZGQk4GjzHzduHIsWLaq0fOPGjfWqt2LXwvDw8PLnqu7v0Gcazsb2MSZADBkyhM8//5xvv/0WgNOnT7Nz504uuOACdu/eTWZmJgCLFy+udvthw4bx+uuvA/Dqq68yYsSIJonbuIdd+TQmQMTGxvLSSy8xffp0Cgsdd7o/9NBD9O7dm2eeeYYJEyYQExPDoEHV97+YP38+N910E4899hixsbH885//bMrwTSMTX/mKlpKSojaZi/El33zzDfHx8Z4Oo07y8vKIiopCVbn99tvp1asXd911l6fDMi6q7m9ORDaoasqP17VmH2MMzz33HMnJySQkJHDixAluuaXx+2AY72LNPsYY7rrrLjvTDzB25m+MMQHIkr8xxgQgS/7GGBOALPkbY0wAsuRvjB97+OGHSUhIICkpieTkZL780jGi+ujRo+nTpw/JyckkJyczderUGuu46qqrGDp06Dn3tX79+krj/RjvZr19jPFTX3zxBStWrCA9PZ3w8HAOHz5MUVFR+fJXX32VlJQq3b8rOX78OOnp6URFRbFnz55a54dNSUk5Z33Ge9iZvzFe4u2N2Qyf+ynd736P4XM/5e2N2Q2q78CBA8TExJSPuxMTE8N5553nUh1Lly7lyiuv5Lrrrisf2gHgzTffJDExkf79+zNy5EgAVq9eXT4U9Lp16xg2bBgXXnghw4YNY8eOHQ16L6bxWfI3xgu8vTGb3y/7muzj+SiQfTyf3y/7ukEfAJdddhn79u2jd+/e3Hbbbfz73/+utPyGG24ob/aZM2dOtXUsWrSI6dOnM3369EoDwj344IN89NFHbN68meXLl1fZ7oILLmDNmjVs3LiRBx98kHvuuafe78O4hzX7GOMFHvtoB/nFpZXK8otLeeyjHUy6sFO96oyKimLDhg189tlnrFq1imuvvZa5c+eSmpoKnLvZ5+DBg3z77beMGDECESEkJIQtW7aQmJjI8OHDSU1NZdq0aVx99dVVtj1x4gQzZ85k165diAjFxTZvtrexM39jvMD+4/kulddVcHAwo0eP5oEHHuCpp55i6dKldd528eLFHDt2jO7du9OtWzcyMzPLm36effZZHnroIfbt20dycjJHjhyptO29997LJZdcwpYtW3j33XcpKLB5s72NJX9jvMB5raqfV7qm8rrYsWMHu3btKn+9adMmunbtWuftFy1axIcffkhmZiaZmZls2LChPPl/9913DB48mAcffJCYmBj27dtXadsTJ07QqZPjG8tLL71U7/dg3MeSvzFeYM74PjQLDa5U1iw0mDnj+9S7zry8PGbOnEnfvn1JSkpi27Zt3H///eXLK7b5jx07ttK2mZmZfP/99wwZMqS8rHv37rRo0YIvv/ySOXPm0K9fPxITExk5ciT9+/evtP1vf/tbfv/73zN8+HBKSys3ZxnvYEM6G+Mmrg7p/PbGbB77aAf7j+dzXqtmzBnfp97t/SYwuTKks13wNcZLTLqwkyV702Tc2uwjIhEisk5ENovIVhF5wFneRkQ+FpFdzsfW7ozDGGNMZe5u8y8ELlXV/kAyMEFEhgB3AytVtRew0vnaGGNME3Fr8leHPOfLUOePAlcBLzvLXwYmuTMOY4wxlbm9t4+IBIvIJuAQ8LGqfgm0V9UDAM7Hdu6OwxhjzFluT/6qWqqqyUAcMEhEEuu6rYjMEpH1IrI+JyfHfUEaY0yAabJ+/qp6HFgNTAAOikhHAOfjoRq2WaCqKaqaEhsb21ShGuNX3nrrLUSE7du3l5fNmTOHhISEasf0Wb58OXPnzm3KEI0HuLu3T6yItHI+bwaMBbYDy4GZztVmAu+4Mw5jAtmiRYsYMWJEpVE5//GPf5Cens5jjz1Wad2SkhImTpzI3XdbHwx/5+4z/47AKhHJAL7C0ea/ApgLjBORXcA452tjAlvGG/BkItzfyvGY8UaDq8zLy+Pzzz/nhRdeKE/+EydO5NSpUwwePJjFixeTmprKb37zGy655BJ+97vf8dJLLzF79mzAMbjb5MmT6d+/P/379+c///kPAJMmTWLgwIEkJCSwYMGCBsdpmp5bb/JS1QzgwmrKjwBj3LlvY3xKxhvwbhoUOwdyO7HP8RogaVq9q3377beZMGECvXv3pk2bNqSnp7N8+XKioqLYtGkTAB988AE7d+7kk08+ITg4uNJYPGlpaYwaNYq33nqL0tJS8vIcnfdefPFF2rRpQ35+PhdddBFTpkyhbdu29Y7TND0b28cYb7DywbOJ/4zifEd5AyxatIjrrrsOgOuuu67SmPwVXXPNNQQHB1cp//TTT7n11lsBxwihLVu2BGD+/Pn079+fIUOGsG/fvkoDyBnfYMM7GOMNTmS5Vl4HR44c4dNPP2XLli2ICKWlpYgIjz76aJV1IyMj61zv6tWr+eSTT/jiiy9o3rw5o0ePtiGbfZCd+RvjDVrGuVZeB0uWLGHGjBns3buXzMxM9u3bR/fu3Vm7dm2d6xgzZgx///vfASgtLSU3N5cTJ07QunVrmjdvzvbt2/nvf/9b7xiN51jyN8YbjLkPQn80dn9oM0d5PS1atIjJkydXKpsyZQqvvfZaneuYN28eq1atol+/fgwcOJCtW7cyYcIESkpKSEpK4t5776007LPxHTakszFu4uqQzmS84WjjP5HlOOMfc1+DLvaawGNDOhvji5KmWbI3TcaafYwxJgBZ8jfGmABkyd8YYwKQJX9jjAlAlvyNMSYAWfI3xpgAZMnfGD92/PhxnnnmmVrXyczMrNONX5mZmSQm1nkuJuOCbt26cfjw4SbdpyV/Y/xYYyZ/418s+Rvjx+6++26+++47kpOTmTNnDnPmzCExMZF+/fqxePHi8nU+++wzkpOTefLJJ8nMzOTiiy9mwIABDBgwoHwM/3OpabvVq1czatQopk2bRu/evbn77rt59dVXGTRoEP369eO7774DIDU1lV/96ldcfPHF9O7dmxUrVgCwdetWBg0aRHJyMklJSeUjiL7yyivl5bfccgulpaUA/POf/6R3796MGjWKX/7yl+VzE1QnNTWVtLQ0hg0bRo8ePViyZAkABw4cYOTIkSQnJ5OYmMhnn30GwK233kpKSgoJCQn88Y9/LK+nW7du3HPPPQwdOpSUlBTS09MZP348PXv25Nlnny0/DiNHjmTy5Mn07duXX/3qV5SVlVWJqbr3VVpaSmpqavnv7sknn6zT76RWquoTPwMHDlRjfMm2bds8HYLu2bNHExISVFV1yZIlOnbsWC0pKdEffvhBO3furPv379dVq1bpT3/60/JtTp06pfn5+aqqunPnTj3zv1exrurUtN2qVau0ZcuWun//fi0oKNDzzjtP77vvPlVV/etf/6p33HGHqqrOnDlTx48fr6Wlpbpz507t1KmT5ufn6+zZs/WVV15RVdXCwkI9ffq0btu2Ta+44gotKipSVdVbb71VX375Zd2/f7927txZDx06pIWFhTps2DC9/fbba4x55syZOnXqVC0tLdWtW7dqz549VVX18ccf14ceekhVVUtKSjQ3N1dVVY8cOVJeNmrUKN28ebOqqnbt2lWfeeYZVVW98847tV+/fpqbm6uHDh3S2NjY8uMQHh6u3333nZaUlOjYsWP1zTffLN8+Jyenxve1fv16HTt2bHncx44dq/b9VPc3B6zXanKqDe9gTIBYu3Yt06dPJzg4mPbt2zNq1Ci++uorWrRoUWm94uJiZs+ezaZNmwgODmbnzp11qr+27S666CI6duwIQM+ePbnssssA6NevH6tWrSpfb9q0aQQFBdGrVy969OjB9u3bGTp0KA8//DBZWVlcffXV9OrVi5UrV7JhwwYuuugiAPLz82nXrh1ffvklo0eP5syc39dee+054580aRJBQUH07duXgwcPlsd70003UVxczKRJk0hOTgbgjTfeYMGCBZSUlHDgwAG2bdtGUlIS4Jgh7cx7ysvLIzo6mujoaCIiIjh+/DgAgwYNokePHgBMnz6dtWvXMnXq1PJYanpfV155Jbt37+bXv/41P/3pT8uPX0NY8jcmQGgdB3F88sknad++PZs3b6asrIyIiIgGbxceHl7+PCgoqPx1UFAQJSUl5ctEpFKdIsL111/P4MGDee+99xg/fjzPP/88qsrMmTN55JFHKq3/9ttvV6njXCrGduYYjRw5kjVr1vDee+9x4403MmfOHC6++GIef/xxvvrqK1q3bk1qamqleQwqvqcfv98z77G691dRTe8LYPPmzXz00Uc8/fTTvPHGG7z44osuvc8fszZ/Y/xYdHQ0J0+eBBwJbfHixZSWlpKTk8OaNWsYNGhQpXUATpw4QceOHQkKCmLhwoXlbennUt/tKnrzzTcpKyvju+++Y/fu3fTp04fdu3fTo0cP0tLSmDhxIhkZGYwZM4YlS5Zw6NAhAI4ePcrevXsZPHgwq1ev5siRIxQXF/Pmm2+6HAPA3r17adeuHb/85S+5+eabSU9PJzc3l8jISFq2bMnBgwf54IMPXK533bp17Nmzh7KyMhYvXsyIESMqLa/pfR0+fJiysjKmTJnCn/70J9LT0+v1viqyM39j/Fjbtm0ZPnw4iYmJ/OQnPyEpKYn+/fuXz+jVoUMH2rZtS0hICP379yc1NZXbbruNKVOm8Oabb3LJJZfUeZav+m5XUZ8+fRg1ahQHDx7k2WefJSIigsWLF/PKK68QGhpKhw4duO+++2jTpg0PPfQQl112GWVlZYSGhvL0008zZMgQ7r//foYOHUrHjh0ZMGBAvT6EVq9ezWOPPUZoaChRUVH861//onv37lx44YUkJCTQo0cPhg8f7nK9Q4cO5e677+brr78uv/hbUd++fat9X82aNePnP/95+QXi6r4ZuMqt4/mLSGfgX0AHoAxYoKrzROR+4JdAjnPVe1T1/drqsvH8ja9xeTz/AJeamsoVV1xRqQ28oV566SXWr1/PU0891Wh11tfq1at5/PHHy3sxuYM3jedfAvyPqqaLSDSwQUQ+di57UlUfd/P+jTHGVKNJZ/ISkXeAp4DhQJ4ryd/O/I2v8dcz/48++ojf/e53lcq6d+/OW2+95aGIzu3hhx+u0v5/zTXX8Ic//MFDEbmHK2f+TZb8RaQbsAZIBH4DpAK5wHoc3w6OVbPNLGAWQJcuXQbu3bu3SWI1pjH4a/I33suV5N8kvX1EJApYCtypqrnA34GeQDJwAPhLddup6gJVTVHVlDP9do0xxjSc25O/iITiSPyvquoyAFU9qKqlqloGPAcMcnccxhhjznJr8hfHHQwvAN+o6hMVyjtWWG0ysMWdcRhjjKnM3Wf+w4EbgUtFZJPz53LgURH5WkQygEuAu9wchzEB68MPP6RPnz6cf/75zJ0719PhGC/h1q6eqroWqO5e61r79BtjGkdpaSm33347H3/8MXFxcVx00UVMnDiRvn37ejo042E2vIMxXqSkpITDhw9XGu+mIdatW8f5559Pjx49CAsL47rrruOdd95plLqNb7Pkb4yX2Lx5M2PHjmXixImMHTuWzZs3N7jO7OxsOnfuXP46Li6O7OzsBtdrfJ8lf2O8QElJCXfccQd5eXkUFRWRl5fHHXfcUa9xaSqq7j4eV0e9NP7Jkr8xXuD48eMUFRVVKisqKuLYsSr3ProkLi6Offv2lb/OysrivPPOa1Cdxj9Y8jfGC7Rq1YqwsLBKZWFhYbRu3bpB9V500UXs2rWLPXv2UFRUxOuvv14+6YgJbJb8jfECISEhzJs3j6ioKMLCwoiKimLevHkEBwc3uN6nnnqK8ePHEx8fz7Rp00hISGikqI0vs/H8jfES/fv3Z+XKlRw7dozWrVs3OPGfcfnll3P55Zc3Sl3Gf1jyN8aLBAcHExMT4+kwTACwZh9jjAlAlvyNMSYAWfI3xpgAZMnfGGMCkCV/Y4wJQJb8jfFj+/bt45JLLiE+Pp6EhATmzZvn6ZCMl7CunsZ4kYKCAg4fPkxMTAwRERENri8kJIS//OUvDBgwgJMnTzJw4EDGjRtnQzobS/7GeIOSkhLmz5/P0qVLCQoKoqysjClTppCWlkZISP3/TTt27EjHjo6J86Kjo4mPjyc7O9uSv7Hkb4w3mD9/PsuWLaOwsLC8bNmyZYgId93VOBPdZWZmsnHjRgYPHtwo9RnfZm3+xnhYQUEBS5cupaCgoEr5kiVLqpTXR15eHlOmTOGvf/0rLVq0aHB9xvdZ8jfGww4fPkxQUPX/ikFBQRw+fLhB9RcXFzNlyhRuuOEGrr766gbVZfxHnZK/iASLiMvfPUWks4isEpFvRGSriNzhLG8jIh+LyC7nY8PGrTXGh8XExFBWVlbtMlVt0Fg/qsrNN99MfHw8v/nNb+pdj/E/dUr+qloKXFWP+kuA/1HVeGAIcLuI9AXuBlaqai9gpfO1MQEpIiKCKVOmVOndU1O5Kz7//HMWLlzIp59+SnJyMsnJybz//vsNDdn4AVcu+H4uIk8Bi4FTZwpVNb2mDVT1AHDA+fykiHwDdMLxQTLaudrLwGrgd64Ebow/SUtLQ0RYsmRJpd4+v/71rxtU74gRI6qdytEYqesfhoisqqZYVfXSOm7fDVgDJALfq2qrCsuOqWqVph8RmQXMAujSpcvAvXv31ilWY7zBN998Q3x8vEvbNHY/fxNYqvubE5ENqpry43XrfOavqpfUNyARiQKWAneqam5dJ5BW1QXAAoCUlBQ7fTF+LyIigri4OE+HYQKAS/38ReSnQAJQfkqiqg+eY5tQHIn/VVVd5iw+KCIdVfWAiHQEDrkWtjHGmIaoc1dPEXkWuBb4NSDANUDXc2wjwAvAN6r6RIVFy4GZzuczgXdciNkYY0wDudLPf5iqzgCOqeoDwFCg8zm2GQ7cCFwqIpucP5cDc4FxIrILGOd8bYwxpom40uyT73w8LSLnAUeA7rVtoKprcXxLqM4YF/ZtjDGmEbly5r9CRFoBjwHpQCawyB1BGWMaV2lpKRdeeCFXXHGFp0MxXsKV3j5/cj5dKiIrgAhVPeGesAynDkNuNuz9D0R3gM5DoFkbCA33dGTGDUpKSli+fDmLFi0iJyeH2NhYpk+fzsSJExs0qucZ8+bNIz4+ntzc3EaI1viDOv9VichmHDd4LVbV74DCc2xi6ivvILz+M8had7YstBlMXwydB0Oo9f/2JyUlJaSlpZGRkVE+iFteXh5PPPEEn3zyCfPnz2/QB0BWVhbvvfcef/jDH3jiiSfOvYEJCK40+0zEMVzDGyLylYj8HxHp4qa4AlfRafj0ocqJH6A4H16dCvlHPROXcZvly5dXSvxnFBQUkJGRwbvvvtug+u+8804effTRGgePM4Gpzn8NqrpXVR9V1YHA9UASsMdtkQWqghOQsbj6ZaVFsPOjpo3HuN2iRYtqHLa5oKCA1157rd51r1ixgnbt2jFw4MB612H8k6s3eXUDpuHo718K/LbxQwpwZcVQUkuL2rHMJgvFNI2cnJwGLa/N559/zvLly3n//fcpKCggNzeXn/3sZ7zyyiv1rtP4B1du8voSWObc5hpVHaSqf3FbZIEqOAyi2tW8vOvwpovFNInY2NgGLa/NI488QlZWFpmZmbz++utceumllvgN4Fqb/0xVHaCqc1V1948XisjM6jYyLmoeA6N+X/2yVl3hvOSmjce43fTp02scxC0iIoLrr7++iSMygcCVNv/t51jljgbGYgCCQyBhEoyfCxGtzpZ3HwWp79X+rcD4pIkTJ5KUlFTteP5JSUlceeWVjbKf0aNHs2LFikapy/i+xpzAvW5DdZpza94GLroZEic7ev+EhEFYJDSzCc/8UUhICPPnz+fdd9/ltddeK+/nf8MNN3DFFVc0Sj9/Y36sMf+qbMjlxhQS5ri5ywSEkJAQJk+ezOTJkz0digkQjdnx1878jTHGRzRm8v+8EesyxhjjRnVK/iJygYiMcc7IVbF8wpnnqjq7sYMzxhjjHudM/iKShmOylV8DW0TkqgqL/+yuwIwxxrhPXc78fwkMVNVJwGjgXhE5063T2vmNaUTZ2dls2rSJ7OzsRqvz+PHjTJ06lQsuuID4+Hi++OKLRqvb+K669PYJVtU8AFXNFJHRwBIR6Yolf2MaxbZt2/jzn//Mnj17CA0Npbi4mO7du3PPPffQt2/fBtV9xx13MGHCBJYsWUJRURGnT59upKiNL6vLmf8PIlJ+W6nzg+AKIAbo567AjAkU27ZtY9asWWzfvp3CwkLy8vIoLCxk+/btzJo1i23bttW77tzcXNasWcPNN98MQFhYGK1atTrHViYQ1CX5zwB+qFigqiXO+XxHuiUqYwLIn//851pH9XzkkUfqXffu3buJjY3l5z//ORdeeCG/+MUvOHXqVL3rM/7jnMlfVbNU9QcR6Ski4QAiMtp5IXhrbduKyIsickhEtlQou19Esn80obsxASk7O5s9e2ofGX337t31vgZQUlJCeno6t956Kxs3biQyMpK5c+fWqy7jX1zp578UKBWR84EXcEzefq6Bxl8CJlRT/qSqJjt/3nchBmP8Sk5ODqGhobWuExoaWu9hnePi4oiLi2Pw4MEATJ06lfT09HrVZfyLK8m/TFVLgMnAX1X1LqBjbRuo6hrApp4ypgaxsbEUFxfXuk5xcXG9h3Xu0KEDnTt3ZseOHQCsXLmywReQjX9wJfkXi8h0YCZwZmjA2k9ZajZbRDKczUI1jlYmIrNEZL2IrG/IhBbGeKtOnTrRvXv3Wtfp0aMHnTp1qvc+/va3v3HDDTeQlJTEpk2buOeee+pdl/EfriT/nwNDgYdVdY+IdAfqMyvE34GeQDJwAKhxQhhVXaCqKaqa0pAJLYzxZvfcc0+t4/n//vc1zO9QR8nJyaxfv56MjAzefvttWre20WGNa8l/nKqmqeoiAFXdA+S7ukNVPaiqpapaBjwHDHK1DmP8Sd++fVmwYAHx8fGEh4cTFRVFeHg48fHxLFiwwJppjFu4MqTzTGDej8pSqymrlYh0VNUDzpeTgS21rW9MIOjbty8LFy4kOzu7fDz/hjT1GHMu50z+znb+64HuIrK8wqJo4Mg5tl2EY0iIGBHJAv4IjHbeNKZAJnBLvSI3xgeoKiJ1vxG+U6dOlvRNvai6NqVKXc78/4OjbT6Gyu3zJ4GMcwQzvZriF+ocnTE+LCIigiNHjtC2bVuXPgCMcZWqcuTIkRqvHVXnnMlfVfcCe3Fc7DXG1FFcXBxZWVn17qNvjCsiIiKIi4ur8/p1afZZq6ojROQkladqFEBVtYXrYRrj/0JDQ8/ZjdMYT6nLmf8I52O0+8MxxhjTFBpzGkdjjDE+wpK/McYEIEv+xhgTgCz5G2NMALLkb4wxAciSvzHGBCBL/sYYE4As+RtjTACy5G+MMQHIkr8xxgQgS/7GGBOALPkbY0wAsuRvjDEByJK/McYEIEv+xhgTgCz5G2NMAHJr8heRF0XkkIhsqVDWRkQ+FpFdzsfW7ozBGGNMVe4+838JmPCjsruBlaraC1jpfG2MMaYJuTX5q+oa4OiPiq8CXnY+fxmY5M4YjDHGVGD9aEkAABK3SURBVOWJNv/2qnoAwPnYrqYVRWSWiKwXkfU5OTlNFqAxxvg7r77gq6oLVDVFVVNiY2M9HY4xxvgNTyT/gyLSEcD5eMgDMRhjTEDzRPJfDsx0Pp8JvOOBGIwxJqC5u6vnIuALoI+IZInIzcBcYJyI7ALGOV8bY4xpQiHurFxVp9ewaIw792uMMaZ2Xn3B1xhjjHtY8jfGmADk1mYf42PyDsGJLDh1GNonQERLCI/ydFTGGDew5G8ccg/AvybC4Z2O10HBMPFpiL/SPgCM8UPW7GOg6BSsfuRs4gcoK4Xlt0NhrufiMsa4jSV/A4UnYcd7VcvLSuGHLVXLjTE+z5K/AQmGyBqGWIq0YTWM8UeW/A1ExcIl91Qtj+kNLTs1fTzGGLez5G8cul0M177i6OUTGQMDUmHGcoiqcdBVY4wPs94+xqFZK0fPni5DHG394S0grLmnozLGuIklf1OZtfEbExCs2ccYYwKQJX9jjAlAlvyNMSYAWfI3xpgAZMnfGGMCkCV/Y4wJQJb8jTEmAFnyN8aYAOSxm7xEJBM4CZQCJaqa4qlYjDEm0Hj6Dt9LVPWwh2MwxpiAY80+xhgTgDyZ/BX4XxHZICKzPBiHMcYEHE82+wxX1f0i0g74WES2q+qaiis4PxRmAXTp0sUTMRpjjF/y2Jm/qu53Ph4C3gIGVbPOAlVNUdWU2FgXR5ssOgW52bDxVfh8PhzIgFNHGiN0Y4zxeR458xeRSCBIVU86n18GPNhoOyg8BTs/hLdmQVnJ2fJuI2Hqi46Zq4wxJoB56sy/PbBWRDYD64D3VPXDRqv9dA4s+0XlxA+QuQa+eBpKChttV8YY44s8cuavqruB/m7bwcZXQcuqX7b+BRhyK0S3d9vujTHG2/lfV8+yUjieWfPywlzQ0iYLxxhjvJH/Jf+gYOhxac3LY3pBkKfvbQsMR04VknOygPyiknOvbIxpUv6ZBXteCtEdAOF0vxspioih5e73kMx/w9gHIapd+aqlZcqRvELyi0sJDwkiKjyEqIjQKlXm5hdzqqiE4pIymoWF0DYyjKAgacI35TuOnipi9Y5DvPj5Ho6fLmZYz7bMvrQXHVqEExYS7OnwjDGAqKqnY6iTlJQUXb9+fd1WVqUo9yBZuSX85d8HOJRXzHVJLbm0Zwtat4iGZq0AKCwuJSPrBLMXpXMwt5CQIGH6oC7cNa43bSLDyqs7klfIn1ZsY/nm/ZQpxLVuxj9uHEif9tGEBPvfl6eGOHqqiHuWfc2HW3+oVB4RGsQ7t4+gT4doD0VmTGASkQ3VjZ3mn5lLhMPakssXfM17Ww7xVeYx/md5Jou/KaIotEX5asdOF/OzF77kYK6j909JmbLwv3t5L2M/ZWWOD8XikjL++Xkmb29yJH6ArGP5TH/uvxw9VdTkb83bZR87XSXxAxQUl3HvO1s4ftqOmTHewD+TP/DBlgMUFFfu8fP82t2cyD+bfDZ+f4zCkqq9ghb+d295Yj+WX8Tir/ZVWSc3v4TMI6cbOWrft2xjdo3L1u05Sn6xXWw3xhv4bfIPlqpvLUgqt9HX1GYfHCQgFderYR9+e/TqL1hqvg4ilQ+rMcaD/DZ9jU9sT1R45evZt43uSevmZ9vykzu3onlY1QuQqcO60dbZ5t+6eRg/G9K1yjptI8Po3Lp5I0ftOlXlcF4hh3K9o1fNlIFxNS4b3rMtzUL9s4+BMb7Gb5N/bGgRH9x+EdendGR0n1iev+4CrkpsQ0iFs/02kaG8PmsI3WMiAcdFydtG9+Syvh0Q5xlsaHAQNwzuwowhXQlznur3aR/N67OG0DYqvOnfWAXHThXx9qZsfvb8l1zxt7X8cflW9h/PL79e4QkdWkQwZUDVD4Do8BAeuCqRls2r9qQyxjQ9/+ztU5ALq/4MGa9T1Pcaipq1JWr3R5CbBbesgeiOlVY/fLKQgpJSQoODaBERSrNqvg2cKizhZEEJJWVlNAsN9njizyss5m8rv+Ufa3ZXKm8bGcb7aRfTvmWEhyKDo6cK2fj9cZ7/bA/H84u4uFcsqcO60T46nGBrKzOmSdXU28c/v4MXnYSvFkBZKWEbFhBWcdnm12FoGgSfTfAx0edO5JHhIUSGe8/hyiso5YW1e6qUHzlVxOtffc/tl55PSE0XK9ysTWQ4Y+Lbk9KtNSWlSlRECOHWv98Yr+Kfp2G5PziGeajOns+gxPd76RzJK6Skhuad/+45yulCz/eqadksjLZR4Zb4jfFC/pn8I9vWvCymFwR7tsmmMdTWdt69bSThIf75qzXGNA7vacdoTOHR0O1iyPyscnlwKAy+BUIqNATlHYIDm2HXx9AyDhImQ/M2EBbZtDG7KCo8hJG9Yliz63Cl8pAg4ZcjexAe6gNn23mH4OAW2PEhtDgPEq+GZm0gPMrTkRnj9/zzgi9A3kH44HfwzXJHE1Db8+Gqp6BjMoQ2c6yTux9e+ikcrXDRVIJg6j/h/HEQ7t0fADknC3novW28l3GAkjKlW9vm/L8pSSR1bun9XSpP/gD/mgg5O86WicBVf4f4K33iA+BwXiFbsk+w98hpRvWJJSYqvEr3YmM8raYLvv6b/AEKTjou/paVQEhEpQHdKMhF374N2f5u1e2CguGODMc3AS+XV1BMXmEJJWVKREgwbaPCyrupeq3CPHjvN5CxuOoyEUjbDK2r3lvhTQ7nFTLzxXVs3Z8LOMJecGMKl17QznGToDFeIrDG9jkjItrRnNCqS+XED1CUh+x8v/rtykphRw3LvExURCgdWjYjrnVzYqLDvT/xAxSehK3Lql+mWvMyL7I1+0R54gdH2Pe9s4Wjp2yWOOMb/Dv516a0uOYeQQAnDzZdLIGmrMRx/GtysurAcN5m79GqPcYO5hbgI1+kjQng5B8SAS06OZ63T4T+06HPTyDE2RPo/DGei83fhYRDmx41L+81vuliqaeRvWP58ZessfHtiajmBkFjvFHgXp2KjIGfPOpIRKcOwd4voF0CXPw/sO1daNvL9TrzcmD/Rtj+LrTqBknTIDIWQj13t61XioyFyx6C16+vuiymN3RIbPqYXBQbFc7zM1K49+0t/JBbwLi+7fnTVYm0qGYioHMpKC7lSF4hb23MJutYPhMSO5DYqSUxHr6L3BFXEW9tzCLrWD7jEzrQr1PLOt0UaRqmsLiUI6eKeGfTfjIP5zGubweSu7Rq1L8Jj13wFZEJwDwgGHheVefWtn59Lvjqif2OxJ53CNonQPN2SGiFf86TP8Ci6xwJ+4zgMLjxLegy1HHht67ycuC1abA/veF1BYKCE5C5Fv73/zp6WwWHQsLVMO6BKsNveKuyMuXIqSIUJSI0uF6Jv6S0jK8yjzLjxXUUl579X7ywSyuem5HisQ+A0jJlfeZRbnxhHUWlZ4c979+5JS/MvMjjH0z+rKxM2bTvONOf+2+lIefjO0bzr5sGE+vih69XXfAVkWDgaeAnQF9guoj0bcx9lJ38Af77FLJgFPLaNcjzl8LpQxVWKIOtb1VO/AClRbBsFpyu3H++9p2dqSu9cnlpESz9hWt1BYqIltDncrjpQ7hzC9zxNfz0CZ9J/OAYEjw2Opx20RH1SvwAR08XkbZoU6XED7Dx++N8vO0gnjo5O3qqkLTXN1ZK/ACb953g/a8PeHTwQH935FQRaa9vrDLXyDcHTvLWxixKS6vOQVIfnmrzHwR8q6q7VbUIeB24qjF3IKXFyJf/OFtw8gf49pOzrwuOwddvVL9xbjYUnqr7zvKPwdfVdFsEOHnA0bXRVCUCUe2hVWdo0dEn+vY3tryCEnLyqu8htHRDFifya7kw7kZ5haXlM9z92Fvp2R6LKxCcLioh61h+tcve3rif44107D2V/DsBFafHynKWVSIis0RkvYisz8nJcW0PIo6LuhVFtKywPAhCmtW8vSvNNEFBVfdVaXngXloxtavtnoCI0OAaJxxyt5BzxeULXYp9VO1/E0GNduw9lfyri77K90hVXaCqKaqaEhsb69IONKQ5evlfHO3ugHYZBp0HnV2hWWsY/KvqN26f6NrwDs1aw5DbaqgrAcI8P+mL8U6RYSGc3676bzw3De9W7+akhmoeFkyf9tHVLrtpRDebl8GNmoeF0D+uZbXLUod1o3VkWLXLXOWp5J8FdK7wOg7Y35g7CIpsg/a8FE3biKZtginPIy3Oq7xSt+EwMLVyWcs4uHahozeQK7oMgZSbKpe16ATTFjp6txhTjZjocJ6bkUL7FpUv4s0c2pULu7T2UFTQNiqcZ28cSMcfzQtxw+AuDOzaxkNRBYY2kWH87foBxLWu3DIxdWAcI3o1Xi7xSG8fEQkBdgJjgGzgK+B6Vd1a0zb1Gt6hLvKPO35+2Oy42NiqK0S3r2ddxyD/BPyQ4WjLbt3NcWexfUU2tVBVcvIK+f7oaXJOFpJ4XktaNAuhZbPGOcNrrLgSzmtJy4gQWjb3bFyBIudkAVnH8vkht4D4Di1o1TyUVvU49l41mYuqlojIbOAjHF09X6wt8btVs1aOnzbdGqGu1o6fxqjLBAwRoV10BO2ivet+EG+NK1DERkcQ68Zj77Erkar6PuAbA+gYY4yfCdzhHYwxJoBZ8jfGmABkyd8YYwKQJX9jjAlAPjOTl4jkAHvruXkM4MsD7Phy/L4cO/h2/L4cO/h2/N4Ue1dVrXKDgM8k/4YQkfXV9XP1Fb4cvy/HDr4dvy/HDr4dvy/Ebs0+xhgTgCz5G2NMAAqU5L/A0wE0kC/H78uxg2/H78uxg2/H7/WxB0SbvzHGmMoC5czfGGNMBZb8jTEmAPl98heRCSKyQ0S+FZG7PR2PK0QkU0S+FpFNIuKG8awbl4i8KCKHRGRLhbI2IvKxiOxyPnpukPpzqCH++0Uk2/k72CQil3syxpqISGcRWSUi34jIVhG5w1nu9ce/lth95dhHiMg6EdnsjP8BZ7lXH3u/bvN3ThS/ExiHYwKZr4DpqrrNo4HVkYhkAimq6i03i9RKREYCecC/VDXRWfYocFRV5zo/fFur6u88GWdNaoj/fiBPVR/3ZGznIiIdgY6qmi4i0cAGYBKQipcf/1pin4ZvHHsBIlU1T0RCgbXAHcDVePGx9/czf7dPFG/OUtU1wNEfFV8FvOx8/jKOf2qvVEP8PkFVD6hquvP5SeAbHPNie/3xryV2n6AOec6Xoc4fxcuPvb8n/zpNFO/FFPhfEdkgIrM8HUw9tVfVA+D4JwfaeTie+pgtIhnOZiGv+upeHRHpBlwIfImPHf8fxQ4+cuxFJFhENgGHgI9V1euPvb8n/zpNFO/FhqvqAOAnwO3OZgnTtP4O9ASSgQPAXzwbTu1EJApYCtypqrmejscV1cTuM8deVUtVNRnHfOSDRCTR0zGdi78nf7dPFO9Oqrrf+XgIeAtHM5avOehs0z3TtnvIw/G4RFUPOv+xy4Dn8OLfgbO9eSnwqqoucxb7xPGvLnZfOvZnqOpxYDUwAS8/9v6e/L8CeolIdxEJA64Dlns4pjoRkUjnxS9EJBK4DNhS+1ZeaTkw0/l8JvCOB2Nx2Zl/XqfJeOnvwHnR8QXgG1V9osIirz/+NcXuQ8c+VkRaOZ83A8YC2/HyY+/XvX0AnN3D/srZieIf9nBIdSIiPXCc7YNjruXXvD12EVkEjMYxnO1B4I/A28AbQBfge+AaVfXKi6o1xD8aR7ODApnALWfacb2JiIwAPgO+BsqcxffgaDv36uNfS+zT8Y1jn4Tjgm4wjhPqN1T1QRFpixcfe79P/sYYY6ry92YfY4wx1bDkb4wxAciSvzHGBCBL/sYYE4As+RtjTACy5G+MMQHIkr8xdSQivxKRGZ6Ow5jGYP38jfECIhKiqiWejsMEDjvzNwFPRGY4R47cLCILRaSriKx0lq0UkS7O9e4Xkf9TQx3xIrKuwutuIpLhfH6fiHwlIltEZIFzOANEZLWI/FlE/o1j/HdjmowlfxPQRCQB+ANwqar2x5GEn8IxoUsS8Cow/1z1qOo3QJhzWA6Aa3Hc2g/wlKpe5JwgphlwRYVNW6nqKFX12hErjX+y5G8C3aXAkjOzpTnHXhkKvOZcvhAYUce63sAx+xQ4kv9i5/NLRORLEfnaub+ECtssxhgPsORvAp1w7jke6nphbDEwTUR645jgaZeIRADPAFNVtR+OoYkjKmxzytWAjWkMlvxNoFuJI2G3Bcek28B/cAz/DXADjjlZz0lVvwNKgXs5e0Z/JtEfdk5WMrWR4jamQUI8HYAxnqSqW0XkYeDfIlIKbATSgBdFZA6QA/zchSoXA48B3Z31HxeR53AMV5yJY44JYzzOunoaY0wAsmYfY4wJQNbsY4yLRORpYPiPiuep6j89EY8x9WHNPsYYE4Cs2ccYYwKQJX9jjAlAlvyNMSYAWfI3xpgA9P8B771ngFquHP4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(data=comb_div[~comb_div.index.isin(MULTISP_LABELS)],\n",
    "                x='coi_var',y='its2_var',hue='region',size='total_ampseq_nsamples');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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>coi_len</th>\n",
       "      <th>coi_var</th>\n",
       "      <th>coi_nseq</th>\n",
       "      <th>its2_len</th>\n",
       "      <th>its2_var</th>\n",
       "      <th>its2_nseq</th>\n",
       "      <th>total_ampseq_len</th>\n",
       "      <th>total_ampseq_var</th>\n",
       "      <th>total_ampseq_nsamples</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "      <td>27.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>601.666667</td>\n",
       "      <td>6.370370</td>\n",
       "      <td>3.333333</td>\n",
       "      <td>536.888889</td>\n",
       "      <td>9.259259</td>\n",
       "      <td>3.370370</td>\n",
       "      <td>8943.111111</td>\n",
       "      <td>178.444444</td>\n",
       "      <td>3.703704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>51.456479</td>\n",
       "      <td>8.043861</td>\n",
       "      <td>1.208941</td>\n",
       "      <td>143.064088</td>\n",
       "      <td>17.790647</td>\n",
       "      <td>1.148516</td>\n",
       "      <td>664.986370</td>\n",
       "      <td>111.377160</td>\n",
       "      <td>1.137298</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>402.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>385.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>7196.000000</td>\n",
       "      <td>26.000000</td>\n",
       "      <td>2.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>616.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>468.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>8639.500000</td>\n",
       "      <td>86.000000</td>\n",
       "      <td>3.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>616.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>491.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>9104.000000</td>\n",
       "      <td>142.000000</td>\n",
       "      <td>4.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>619.000000</td>\n",
       "      <td>10.500000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>573.000000</td>\n",
       "      <td>3.500000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>9390.500000</td>\n",
       "      <td>271.000000</td>\n",
       "      <td>5.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>620.000000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>1024.000000</td>\n",
       "      <td>56.000000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>9739.000000</td>\n",
       "      <td>382.000000</td>\n",
       "      <td>5.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          coi_len    coi_var   coi_nseq     its2_len   its2_var  its2_nseq  \\\n",
       "count   27.000000  27.000000  27.000000    27.000000  27.000000  27.000000   \n",
       "mean   601.666667   6.370370   3.333333   536.888889   9.259259   3.370370   \n",
       "std     51.456479   8.043861   1.208941   143.064088  17.790647   1.148516   \n",
       "min    402.000000   0.000000   1.000000   385.000000   0.000000   2.000000   \n",
       "25%    616.000000   1.000000   2.000000   468.000000   0.000000   2.000000   \n",
       "50%    616.000000   2.000000   4.000000   491.000000   1.000000   3.000000   \n",
       "75%    619.000000  10.500000   4.000000   573.000000   3.500000   4.000000   \n",
       "max    620.000000  32.000000   5.000000  1024.000000  56.000000   5.000000   \n",
       "\n",
       "       total_ampseq_len  total_ampseq_var  total_ampseq_nsamples  \n",
       "count         27.000000         27.000000              27.000000  \n",
       "mean        8943.111111        178.444444               3.703704  \n",
       "std          664.986370        111.377160               1.137298  \n",
       "min         7196.000000         26.000000               2.000000  \n",
       "25%         8639.500000         86.000000               3.000000  \n",
       "50%         9104.000000        142.000000               4.000000  \n",
       "75%         9390.500000        271.000000               5.000000  \n",
       "max         9739.000000        382.000000               5.000000  "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# stats for good species (not multiple species with same label)\n",
    "comb_div[~comb_div.index.isin(BAD_SPP) & \n",
    "         (comb_div.total_ampseq_nsamples >= MIN_SAMPLES_DIV)] \\\n",
    "    .describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}