{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lassa virus glycoprotein PacBio sequencing \n",
    "This example shows how to use [alignparse](https://jbloomlab.github.io/alignparse/) to process circular consensus sequences that map to multiple different targets. Here, our two targets are a wildtype and a codon optimized sequence for the Lassa virus (LASV) glycoprotein from the Josiah strain. For such targets we do not expect large internal deletions, so we use alignment settings optimized for codon-level mutations as these will be the settings used when analyzing mutant LASV GP sequences in later experiments.\n",
    "\n",
    "Here we analyze a snippet of the full data set of circular consensus sequences so that the example is small and fast.\n",
    "\n",
    "In the other included example notebooks ([RecA deep mutational scanning libraries](https://jbloomlab.github.io/alignparse/recA_DMS.html) and [Single-cell virus sequencing](https://jbloomlab.github.io/alignparse/flu_virus_seq_example.html)), we align and parse the PacBio reads using the single [align_and_parse](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.align_and_parse) function. Here we use separate functions to do the aligning and parsing steps. This illustrates an additional use case and shows how this package could be used to parse alignments generated elsewhere as long as they have a [cs tag](https://lh3.github.io/minimap2/minimap2.html#10) and are in the SAM file format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up for analysis\n",
    "Import necessary Python modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:27.680590Z",
     "iopub.status.busy": "2024-05-23T18:06:27.680149Z",
     "iopub.status.idle": "2024-05-23T18:06:31.717771Z",
     "shell.execute_reply": "2024-05-23T18:06:31.716937Z",
     "shell.execute_reply.started": "2024-05-23T18:06:27.680557Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import tempfile\n",
    "import warnings\n",
    "\n",
    "import Bio.SeqIO\n",
    "\n",
    "import pandas as pd\n",
    "import numpy\n",
    "\n",
    "from plotnine import *\n",
    "\n",
    "import pysam\n",
    "\n",
    "import alignparse.ccs\n",
    "import alignparse.minimap2\n",
    "import alignparse.targets\n",
    "import alignparse.cs_tag\n",
    "import alignparse.consensus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppress warnings that clutter output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.720334Z",
     "iopub.status.busy": "2024-05-23T18:06:31.719872Z",
     "iopub.status.idle": "2024-05-23T18:06:31.723583Z",
     "shell.execute_reply": "2024-05-23T18:06:31.722917Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.720306Z"
    }
   },
   "outputs": [],
   "source": [
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directory for output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.724731Z",
     "iopub.status.busy": "2024-05-23T18:06:31.724436Z",
     "iopub.status.idle": "2024-05-23T18:06:31.728224Z",
     "shell.execute_reply": "2024-05-23T18:06:31.727532Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.724705Z"
    }
   },
   "outputs": [],
   "source": [
    "outdir = \"./output_files/\"\n",
    "os.makedirs(outdir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Color palette for plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.729604Z",
     "iopub.status.busy": "2024-05-23T18:06:31.729285Z",
     "iopub.status.idle": "2024-05-23T18:06:31.733335Z",
     "shell.execute_reply": "2024-05-23T18:06:31.732450Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.729573Z"
    }
   },
   "outputs": [],
   "source": [
    "CBPALETTE = (\"#999999\", \"#E69F00\", \"#56B4E9\", \"#009E73\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Target amplicons\n",
    "We have performed sequencing of several LASV GP amplicons that include the glycoprotein sequence along with a PacBio index and several other features. Here we analyze reads mapping to two of these amplicons.\n",
    "The amplicons are defined in [Genbank Flat File format](https://www.ncbi.nlm.nih.gov/genbank/samplerecord/).\n",
    "If there are multiple targets, they can all be defined in a single Genbank file (as we did in the [Single-cell virus sequencing](https://jbloomlab.github.io/alignparse/flu_virus_seq_example.html) example) or each target can be defined in its own Genbank file (as we've done here).\n",
    "\n",
    "First, let's just look at the files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.735845Z",
     "iopub.status.busy": "2024-05-23T18:06:31.735032Z",
     "iopub.status.idle": "2024-05-23T18:06:31.742332Z",
     "shell.execute_reply": "2024-05-23T18:06:31.741750Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.735790Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LOCUS       LASV_Josiah_WT          1730 bp ds-DNA     linear       14-JUN-2019\n",
      "DEFINITION  .\n",
      "ACCESSION   \n",
      "VERSION     \n",
      "SOURCE      Kate Crawford\n",
      "  ORGANISM  .\n",
      "COMMENT     PacBio amplicon for LASV Josiah WT sequence\n",
      "FEATURES             Location/Qualifiers\n",
      "     T2A             85..147\n",
      "                     /label=\"T2A\"\n",
      "     WPRE            1639..1730\n",
      "                     /label=\"WPRE\"\n",
      "     ZsGreen         15..84\n",
      "                     /label=\"ZsGreen\"\n",
      "     termini3        1639..1730\n",
      "                     /label=\"3'Termini\"\n",
      "     index           9..14\n",
      "                     /label=\"index\"\n",
      "     leader5         1..8\n",
      "                     /label=\"5' leader\"\n",
      "     termini5        1..147\n",
      "                     /label=\"5'Termini\"\n",
      "     variant_tag5    34..34\n",
      "                     /variant_1=T\n",
      "                     /variant_2=C\n",
      "                     /label=\"5'VariantTag\"\n",
      "     variant_tag3    1702..1702\n",
      "                     /variant_1=G\n",
      "                     /variant_2=A\n",
      "                     /label=\"3'VariantTag\"\n",
      "     spacer          1624..1638\n",
      "                     /label=\"3'Spacer\"\n",
      "     gene            148..1623\n",
      "                     /label=\"LASV_Josiah_WT\"\n",
      "\n",
      "ORIGIN\n",
      "        1 GACTGATANN NNNNcagcga cgccaagaac cagYagtggc acctgaccga gcacgccatc\n",
      "       61 gcctccggcT CCGCCTTGCC CGCTGGATCC GGCGAGGGCA GAGGAAGTCT GCTAACATGC\n",
      "      121 GGTGACGTCG AGGAGAATCC TGGCCCAATG GGACAAATAG TGACATTCTT CCAGGAAGTG\n",
      "      181 CCTCATGTAA TAGAAGAGGT GATGAACATT GTTCTCATTG CACTGTCTGT ACTAGCAGTG\n",
      "      241 CTGAAAGGTC TGTACAATTT TGCAACGTGT GGCCTTGTTG GTTTGGTCAC TTTCCTCCTG\n",
      "      301 TTGTGTGGTA GGTCTTGCAC AACCAGTCTT TATAAAGGGG TTTATGAGCT TCAGACTCTG\n",
      "      361 GAACTAAACA TGGAGACACT CAATATGACC ATGCCTCTCT CCTGCACAAA GAACAACAGT\n",
      "      421 CATCATTATA TAATGGTGGG CAATGAGACA GGACTAGAAC TGACCTTGAC CAACACGAGC\n",
      "      481 ATTATTAATC ACAAATTTTG CAATCTGTCT GATGCCCACA AAAAGAACCT CTATGACCAC\n",
      "      541 GCTCTTATGA GCATAATCTC AACTTtccac ttgtccatcc ccaacTTCAA TCAGTATGAG\n",
      "      601 GCAATGAGCT GCGATTTTAA TGGGGGAAAG ATTAGTGTGC AGTACAACCT GAGTCACAGC\n",
      "      661 TATGCTGGGG ATGCAGCCAA CCATTGTGGT ACTGTTGCAA ATGGTGTGTT ACAGACTTTT\n",
      "      721 ATGAGGATGG CTTGGGGTGG GAGCTACATT GCTCTTGACT CAGGCCGTGG CAACTGGGAC\n",
      "      781 TGTATTATGA CTAGTTATCA ATATCTGATA ATCCAAAATA CAACCTGGGA AGATCACTGC\n",
      "      841 CAATTCTCGA GACCATCTCC CATCGGTTAT CTCGGGCTCC TCTCACAAAG GACTAGAGAT\n",
      "      901 ATTTATATTA GTAGAAGATT GCTAGGCACA TTCACATGGA CACTGTCAGA TTCTGAAGGT\n",
      "      961 AAAGACACAC CAGGGGGATA TTGTCTGACC AGGTGGATGC TAATTGAGGC TGAACTAAAA\n",
      "     1021 TGCTTCGGGA ACACAGCTGT GGCAAAATGT AATGAGAAGC ATGATGAgga attttgtgac\n",
      "     1081 atgctgaggc TGTTTGACTT CAACAAACAA GCCATTCAAA GGTTGAAAGC TGAAGCACAA\n",
      "     1141 ATGAGCATTC AGTTGATCAA CAAAGCAGTA AATGCTTTGA TAAATGACCA ACTTATAATG\n",
      "     1201 AAGAACCATC TACGGGACAT CATGGGAATT CCATACTGTA ATTACAGCAA GTATTGGTAC\n",
      "     1261 CTCAACCACA CAACTACTGG GAGAACATCA CTGCCCAAAT GTTGGCTTGT ATCAAATGGT\n",
      "     1321 TCATACTTGA ACGAGACCCA CTTTTCTGAT GATATTGAAC AACAAGCTGA CAATATGATC\n",
      "     1381 ACTGAGATGT TACAGAAGGA GTATATGGAG AGGCAGGGGA AGACACCATT GGGTCTAGTT\n",
      "     1441 GACCTCTTTG TGTTCAGCAC AAGTTTCTAT CTTATTAGCA TCTTCCTTCA CCTAGTCAAA\n",
      "     1501 ATACCAACTC ATAGGCATAT TGTAGGCAAG TCGTGTCCCA AACCTCACAG ATTGAATCAT\n",
      "     1561 ATGGGCATTT GTTCCTGTGG ACTCTACAAA CAGCCTGGTG TGCCTGTGAA ATGGAAGAGA\n",
      "     1621 TGAGCTAGCT AAACGCGTTG ATCCtaatca acctctggat tacaaaattt gtgaaagatt\n",
      "     1681 gactggtatt cttaactatg tRgctccttt tacgctatgt ggatacgctg \n",
      "//\n",
      "\n",
      "LOCUS       LASV_Josiah_OPT          1730 bp ds-DNA     linear       14-JUN-2019\n",
      "DEFINITION  .\n",
      "ACCESSION   \n",
      "VERSION     \n",
      "SOURCE      Kate Crawford\n",
      "  ORGANISM  .\n",
      "COMMENT     PacBio amplicon for LASV Josiah OPT sequence\n",
      "FEATURES             Location/Qualifiers\n",
      "     T2A             85..147\n",
      "                     /label=\"T2A\"\n",
      "     WPRE            1639..1730\n",
      "                     /label=\"WPRE\"\n",
      "     ZsGreen         15..84\n",
      "                     /label=\"ZsGreen\"\n",
      "     termini3        1639..1730\n",
      "                     /label=\"3'Termini\"\n",
      "     index           9..14\n",
      "                     /label=\"index\"\n",
      "     leader5         1..8\n",
      "                     /label=\"5' leader\"\n",
      "     termini5        1..147\n",
      "                     /label=\"5'Termini\"\n",
      "     variant_tag5     34..34\n",
      "                     /variant_1=T\n",
      "                     /variant_2=C\n",
      "                     /label=\"5'VariantTag\"\n",
      "     variant_tag3     1702..1702\n",
      "                     /variant_1=G\n",
      "                     /variant_2=A\n",
      "                     /label=\"3'VariantTag\"\n",
      "     spacer          1624..1638\n",
      "                     /label=\"3'Spacer\"\n",
      "     gene            148..1623\n",
      "                     /label=\"LASV_Josiah_OPT\"\n",
      "\n",
      "ORIGIN\n",
      "        1 GACTGATANN NNNNcagcga cgccaagaac cagYagtggc acctgaccga gcacgccatc\n",
      "       61 gcctccggcT CCGCCTTGCC CGCTGGATCC GGCGAGGGCA GAGGAAGTCT GCTAACATGC\n",
      "      121 GGTGACGTCG AGGAGAATCC TGGCCCAATG GGCCAGATCG TGACCTTCTT CCAAGAAGTG\n",
      "      181 CCTCATGTGA TTGAGGAGGT GATGAATATC GTGCTGATCG CTTTAAGCGT GCTGGCCGTT\n",
      "      241 CTTAAGGGCC TCTATAACTT CGCCACTTGT GGTTTAGTCG GACTGGTGAC ATTTCTGCTG\n",
      "      301 CTGTGTGGCA GATCTTGTAC CACATCTTTA TACAAGGGCG TGTACGAGCT GCAGACTTTA\n",
      "      361 GAACTGAACA TGGAGACTTT AAACATGACC ATGCCTTTAA GCTGTACCAA GAACAATAGC\n",
      "      421 CACCACTACA TCATGGTGGG CAACGAGACC GGTTTAGAAC TGACACTCAC CAACACCAGC\n",
      "      481 ATTATCAACC ATAAGTTCTG CAACCTCTCC GACGCTCACA AGAAGAATTT ATACGACCAC\n",
      "      541 GCTTTAATGA GCATCATCTC CACCTTCCAT CTCTCCATTC CTAATttcaa ccagtacgag\n",
      "      601 gccatgAGCT GCGACTTTAA CGGCGGCAAG ATCTCCGTGC AGTACAATTT ATCCCATAGC\n",
      "      661 TACGCCGGCG ATGCCGCCAA TCACTGCGGA ACCGTGGCCA ACGGCGTGCT GCAGACATTC\n",
      "      721 ATGAGGATGG CTTGGGGCGG CTCCTATATC GCTTTAGACT CCGGCAGAGG AAACTGGGAC\n",
      "      781 TGTATCATGA CCAGCTACCA ATATTTAATC ATTCAGAACA CCACATGGGA GGACCACTGC\n",
      "      841 CAATTCTCTC GTCCCTCTCC TATCGGCTAT CTGGGACTGC TGTCCCAGAG GACCAGAGAC\n",
      "      901 ATCTACATCT CTCGTAGGCT GCTGGGCACA TTCACTTGGA CTTTAAGCGA CAGCGAAGGC\n",
      "      961 AAAGATACTC CCGGTGGCTA CTGTTTAACA AGATGGATGC TGATCGAGGC CGAGCTCAAG\n",
      "     1021 TGCTTCGGAA ATACCGCCGT GGCCAAATGC AACGAGAAAC ACGACGAGGA GTTCTGCGAC\n",
      "     1081 ATGCTGAGGC TCTTCGACTT CAacaagcaa gccattcaga ggcTGAAGGC CGAAGCCCAG\n",
      "     1141 ATGTCCATCC AGCTGATTAA TAAGGCCGTG AATGCCCTCA TTAACGACCA GCTGATCATG\n",
      "     1201 AAGAACCATT TAAGGGACAT CATGGGCATC CCTTATTGCA ACTACAGCAA ATACTGGTAT\n",
      "     1261 TTAAATCATA CCACCACCGG TCGTACATCC TTACCTAAGT GCTGGCTGGT CAGCAATGGC\n",
      "     1321 TCCTATTTAA ACGAGACACA CTTCTCCGAC GACATCGAGC AGCAAGCCGA CAACATGATC\n",
      "     1381 ACCGAAATGC TCCAGAAGGA GTACATGGAG AGGCAAGGTA AGACTCCTCT GGGTTTAGTG\n",
      "     1441 GATTTATTCG TCTTCAGCAC CTCCTTCTAT TTAATCTCCA TCTTTCTTCA TCTGGTGAAG\n",
      "     1501 ATTCCTACCC ACAGACACAT TGTGGGCAAG AGCTGTCCTA AGCCTCATAG ACTGAACCAC\n",
      "     1561 ATGGGCATCT GTAGCTGCGG TTTATATAAA CAGCCCGGTG TTCCCGTTAA GTGGAAGAGG\n",
      "     1621 TGAGCTAGCT AAACGCGTTG ATCCtaatca acctctggat tacaaaattt gtgaaagatt\n",
      "     1681 gactggtatt cttaactatg tRgctccttt tacgctatgt ggatacgctg \n",
      "//\n",
      "\n"
     ]
    }
   ],
   "source": [
    "target_file_names = [\"LASV_Josiah_WT\", \"LASV_Josiah_OPT\"]\n",
    "\n",
    "targetfiles = [\n",
    "    f\"input_files/{target_file_name}.gb\" for target_file_name in target_file_names\n",
    "]\n",
    "\n",
    "for targetfile in targetfiles:\n",
    "    with open(targetfile) as f:\n",
    "        print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Along with the Genbank files giving the sequences of the amplicons, we have a YAML file specifying how to filter and parse alignments to these amplicons.\n",
    "\n",
    "Below is the text of the YAML file.\n",
    "\n",
    "Additional information about these filters can be found in the [RecA deep mutational scanning libraries](https://jbloomlab.github.io/alignparse/recA_DMS.html) example notebook or the [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets) documentation. \n",
    "\n",
    "A filter setting of `null` indicates this filter is not applied. When filters are missing for a feature, they are automatically set to zero.\n",
    "\n",
    "Here we filter the `gene` based on `mutation_op_counts` by setting the `mutation_nt_counts` filter for the `gene` to `null`. Although we do not expect these sequences to have large indels, we want to confirm this. Filtering on mutation \"operations\" allows us to retain sequences with large indels by only filtering on the number of indels, not the number of nucleotides inserted or deleted. Overall, this example uses very loose filters to allow us to do further analyses of the types of mutations that are arising in these samples.\n",
    "\n",
    "The YAML file also specifies what information is parsed from alignments that are not filtered out.\n",
    "As you can see, for some features we parse the mutations or the full sequence of the feature, along with the accuracy of that feature in the sequencing query (computed from the Q-values). \n",
    "\n",
    "As seen below, we can use YAML syntax to apply one set of filters to multiple targets. Here, we apply the same filters to both targets, but this is not necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.743703Z",
     "iopub.status.busy": "2024-05-23T18:06:31.743415Z",
     "iopub.status.idle": "2024-05-23T18:06:31.749052Z",
     "shell.execute_reply": "2024-05-23T18:06:31.748095Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.743677Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_WT: &LASV_target_parse_specs\n",
      "  query_clip5: 10\n",
      "  query_clip3: 10\n",
      "  termini5:\n",
      "    filter:\n",
      "      clip5: 10\n",
      "      mutation_nt_count: 5\n",
      "      mutation_op_count: null\n",
      "  termini3:\n",
      "    filter:\n",
      "      clip3: 10\n",
      "      mutation_nt_count: 5\n",
      "      mutation_op_count: null\n",
      "  gene:\n",
      "    filter:\n",
      "      mutation_nt_count: null\n",
      "      mutation_op_count: 30\n",
      "    return: [mutations, accuracy]\n",
      "  spacer:\n",
      "    filter:\n",
      "      mutation_nt_count: 1\n",
      "      mutation_op_count: null\n",
      "  index:\n",
      "    return: sequence\n",
      "  variant_tag5:\n",
      "    return: sequence\n",
      "  variant_tag3:\n",
      "    return: sequence\n",
      "\n",
      "LASV_Josiah_OPT: *LASV_target_parse_specs\n",
      "\n"
     ]
    }
   ],
   "source": [
    "lasv_parse_specs_file = \"input_files/lasv_feature_parse_specs.yaml\"\n",
    "with open(lasv_parse_specs_file) as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the amplicons in `feature_parse_specs` into a [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets) object, specifying the features that we require the target to contain. The [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets) in this example have more features specified in their Genbank files than we want to parse, so we set `allow_extra_features` to `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.755591Z",
     "iopub.status.busy": "2024-05-23T18:06:31.755090Z",
     "iopub.status.idle": "2024-05-23T18:06:31.893571Z",
     "shell.execute_reply": "2024-05-23T18:06:31.892700Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.755546Z"
    }
   },
   "outputs": [],
   "source": [
    "targets = alignparse.targets.Targets(\n",
    "    seqsfile=targetfiles,\n",
    "    feature_parse_specs=lasv_parse_specs_file,\n",
    "    allow_extra_features=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we look at the [targets.feature_parse_specs](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.feature_parse_specs), we now see that the previously unspecified specs have been filled in with the defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.895056Z",
     "iopub.status.busy": "2024-05-23T18:06:31.894755Z",
     "iopub.status.idle": "2024-05-23T18:06:31.906140Z",
     "shell.execute_reply": "2024-05-23T18:06:31.905413Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.895023Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_WT: &id001\n",
      "  query_clip5: 10\n",
      "  query_clip3: 10\n",
      "  termini5:\n",
      "    filter:\n",
      "      clip5: 10\n",
      "      mutation_nt_count: 5\n",
      "      mutation_op_count: null\n",
      "      clip3: 0\n",
      "    return: []\n",
      "  termini3:\n",
      "    filter:\n",
      "      clip3: 10\n",
      "      mutation_nt_count: 5\n",
      "      mutation_op_count: null\n",
      "      clip5: 0\n",
      "    return: []\n",
      "  gene:\n",
      "    filter:\n",
      "      mutation_nt_count: null\n",
      "      mutation_op_count: 30\n",
      "      clip5: 0\n",
      "      clip3: 0\n",
      "    return:\n",
      "    - mutations\n",
      "    - accuracy\n",
      "  spacer:\n",
      "    filter:\n",
      "      mutation_nt_count: 1\n",
      "      mutation_op_count: null\n",
      "      clip5: 0\n",
      "      clip3: 0\n",
      "    return: []\n",
      "  index:\n",
      "    return:\n",
      "    - sequence\n",
      "    filter:\n",
      "      clip5: 0\n",
      "      clip3: 0\n",
      "      mutation_nt_count: 0\n",
      "      mutation_op_count: 0\n",
      "  variant_tag5:\n",
      "    return:\n",
      "    - sequence\n",
      "    filter:\n",
      "      clip5: 0\n",
      "      clip3: 0\n",
      "      mutation_nt_count: 0\n",
      "      mutation_op_count: 0\n",
      "  variant_tag3:\n",
      "    return:\n",
      "    - sequence\n",
      "    filter:\n",
      "      clip5: 0\n",
      "      clip3: 0\n",
      "      mutation_nt_count: 0\n",
      "      mutation_op_count: 0\n",
      "LASV_Josiah_OPT: *id001\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(targets.feature_parse_specs(\"yaml\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets). All features specified in the targets' Genbank files will be annotated, even if they are not in `feature_parse_specs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:31.907942Z",
     "iopub.status.busy": "2024-05-23T18:06:31.907578Z",
     "iopub.status.idle": "2024-05-23T18:06:32.445595Z",
     "shell.execute_reply": "2024-05-23T18:06:32.444747Z",
     "shell.execute_reply.started": "2024-05-23T18:06:31.907912Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA18AAAIQCAYAAAB+P9RkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAADXGklEQVR4nOzdd3RURRvH8e+m9wZJCAkkEHrvJfQuvYs0AUEsCKKiomIXlSYqKIq8SC8C0quoKL2JqID0Tmihpie7+/6BrKwJkECym8Dvc86ek7137sxzBw7Mk5k712A2m82IiIiIiIhItnKwdwAiIiIiIiIPAyVfIiIiIiIiNqDkS0RERERExAaUfImIiIiIiNiAki8REREREREbUPIlIiIiIiJiA0q+REREREREbEDJl4iIiIiIiA0o+RIREREREbEBJV8iIvJAmTJlCgaDgWPHjmX62vr161OmTJmsD0pERAQlXyIiucbNpGLHjh0ZKl+tWjUMBgMTJky4bZk///yTTp06ER4ejpubG6GhoTRp0oRx48YB8Ntvv2EwGBg2bNht6zh48CAGg4EXX3wxQ3EdO3YMg8HA6NGjM1Q+N2rRogX+/v6YzWar47t27cJgMBAeHp7mmp9++gmDwcDEiROJiIjAYDDc9TNlyhQb3ZGIiGQFJ3sHICIiWe/gwYNs376diIgIZs6cyTPPPJOmzKZNm2jQoAEFCxbkySefJF++fJw8eZItW7bw2WefMXDgQCpVqkSJEiWYPXs2H3zwQbptzZo1C4AePXpk6z1lVM+ePXnsscdwdXW1Wwy1a9dm5cqV/PXXX5QtW9ZyfOPGjTg5OXHixAlOnTpFWFiY1bmb13766afExsZazq1YsYLZs2czduxY8ubNazkeFRVlg7sREZGsouRLROQBNGPGDIKCghgzZgydOnXi2LFjREREWJUZPnw4vr6+bN++HT8/P6tz58+ft/zcvXt33nzzTbZs2UKNGjXStDV79mxKlChBpUqVsuNWMs3R0RFHR0e7xlC7dm0ANmzYkCb5atGiBT/99BMbNmzgscces5zbsGEDefLkoWTJkpQqVcqqvrNnzzJ79mzatWuX5s9RRERyDy07FBF5AM2aNYtOnTrRqlUrfH19LbNTtzp8+DClS5dOk3gBBAUFWX7u3r27pc7/2rlzJ/v377eUuR9Hjhyhc+fOBAQE4OHhQY0aNVi+fHmacuPGjaN06dJ4eHjg7+9PlSpVrGJL75mvxYsX07JlS/Lnz4+rqyuRkZG8//77GI3GdGPZu3cvDRo0wMPDg9DQUEaOHJmpe6lWrRouLi6W2aybNm7cSN26dalWrZrVOZPJxJYtW4iKisJgMGSqLRERyT2UfImIPGC2bt3KoUOH6Nq1Ky4uLnTo0IGZM2emKRceHs7OnTv566+/7lhfoUKFiIqK4rvvvkuTrNxMerp163ZfMZ87d46oqChWr17Ns88+y/Dhw0lMTKRNmzYsXLjQUu6bb75h0KBBlCpVik8//ZR3332XChUqsHXr1jvWP2XKFLy8vHjxxRf57LPPqFy5Mm+99RZDhw5NU/by5cs88sgjlC9fnjFjxlCiRAleffVVVq5cmeH7cXNzo3LlymzYsMFy7OTJk5w8eZKoqCiioqKskq8///yTa9euWWbMRETkAWUWEZFc4dtvvzUD5u3bt9+x3HPPPWcuUKCA2WQymc1ms3nNmjVmwLxr1y6rcmvWrDE7OjqaHR0dzTVr1jS/8sor5tWrV5uTk5PT1PnFF1+YAfPq1astx4xGozk0NNRcs2bNTN3H0aNHzYB51KhRlmODBw82A+b169dbjl2/ft1cqFAhc0REhNloNJrNZrO5bdu25tKlS9+x/pv9dPToUcux+Pj4NOWeeuops4eHhzkxMdFyrF69embAPG3aNMuxpKQkc758+cwdO3bM1H2+/PLLZsB86tQps9lsNs+ePdvs5uZmTkpKMq9YscLs6OhovnbtmtlsNpvHjx9vBswbN25Mt65Ro0aluScREcl9NPMlIvIASU1NZe7cuXTp0sWyfK1hw4YEBQWlmf1q0qQJmzdvpk2bNuzevZuRI0fSrFkzQkNDWbJkiVXZLl264OzsbLW875dffuH06dNZsuRwxYoVVKtWzWrmx8vLi/79+3Ps2DH27t0LgJ+fH6dOnWL79u2Zqt/d3d3y8/Xr17l48SJ16tQhPj6ev//+26qsl5eX1eYhLi4uVKtWjSNHjmSqzZv3sn79euDGksPKlSvj4uJCzZo1LUsNb55zc3OjSpUqmWpDRERyFyVfIiIPkDVr1nDhwgWqVavGoUOHOHToEEePHqVBgwbMnj0bk8lkVb5q1ap8//33XL58mW3btvHaa69x/fp1OnXqZEl4APLkyUOzZs1YuHAhiYmJwI0lh05OTjz66KP3Hffx48cpXrx4muMlS5a0nAd49dVX8fLyolq1ahQtWpQBAwakea4qPXv27KF9+/b4+vri4+NDYGCgJcG6evWqVdmwsLA0z135+/tz+fLlTN1TrVq1MBgMlvg2btxIrVq1gBtJZKlSpazOVa1aFRcXl0y1ISIiuYuSLxGRB8jN2a1HH32UokWLWj5z587l9OnT/PLLL+le5+LiQtWqVfnwww+ZMGECKSkpzJs3z6pMjx49uHbtGsuWLSM5OZkFCxbQtGlTAgMDs/2+bipZsiT79+9nzpw51K5dmwULFlC7dm3efvvt215z5coV6tWrx+7du3nvvfdYunQpP/zwAyNGjABIk5DebqdE83/e2XU3efLkoUSJEmzYsIHY2Fj++OMPq63ho6Ki2LBhA6dOneLEiRN63ktE5CGgreZFRB4QcXFxLF68mC5dutCpU6c05wcNGsTMmTNp0KDBHeu5ufQtOjra6nibNm3w9vZm1qxZODs7c/ny5SxZcgg3Nv/Yv39/muM3lwTe+lJiT09PunTpQpcuXUhOTqZDhw4MHz6c1157DTc3tzR1rFu3jpiYGL7//nvq1q1rOX706NEsif1OateuzeTJk1mzZg1GozFN8jV79mzWrVtnKSsiIg82zXyJiDwgFi5cSFxcHAMGDKBTp05pPq1atWLBggUkJSUB8PPPP6c7m7NixQqANMsA3d3dad++PStWrGDChAl4enrStm3bLIm9RYsWbNu2jc2bN1uOxcXFMXHiRCIiIizvvYqJibG6zsXFhVKlSmE2m0lJSUm37pszWbfea3JyMl9++WWWxH4ntWvXxmg0Mnr0aIoWLWo1SxgVFUVsbCxffvklDg4OemGyiMhDQDNfIiK5zOTJk1m1alWa4z/++CN58uS57SC+TZs2fPPNNyxfvpwOHTowcOBA4uPjad++PSVKlCA5OZlNmzYxd+5cIiIi6NOnT5o6evTowbRp01i9ejXdu3fH09MzS+5p6NChzJ49m+bNmzNo0CACAgKYOnUqR48eZcGCBTg43PhdYdOmTcmXLx+1atUiODiYffv2MX78eFq2bIm3t3e6dUdFReHv70+vXr0YNGgQBoOB6dOnZ3oZ4b24OZu1efNmevfubXWuWLFi5M2bl82bN1O2bNl037cmIiIPFiVfIiK5zIQJE257rmfPnrd9ZqlRo0Z4eHgwY8YMOnTowOjRo5k3bx4rVqxg4sSJJCcnU7BgQZ599lmGDRuWbjLQsGFDQkJCiI6OvuclhzeTnlvjDA4OZtOmTbz66quMGzeOxMREypUrx9KlS2nZsqWl3FNPPcXMmTP55JNPiI2NJSwsjEGDBjFs2LDbtpcnTx6WLVvGSy+9xLBhw/D396dHjx40atSIZs2a3dM9ZFThwoXJnz8/Z86cSTcpjoqKYsmSJVpyKCLykDCYbfGrPxERkX/88ccflC9fnkmTJtG3b197hyMiImIzeuZLRERs6uY7um4+xyUiIvKw0MyXiIhkieTkZC5dunTb80ePHuWnn35i1KhR5MuXj71791qe5cpNLly4gNFovO15FxcXAgICbBiRiIjkFnrmS0REssSmTZvuuo29h4cHderUYdy4cbky8YIbL6a++dLn9NSrV8+yfbyIiMitNPMlIiJZ4vLly+zcufOOZUqXLk1ISIiNIsoeGzduJCEh4bbn/f39qVy5sg0jEhGR3ELJl4iIiIiIiA3kzjUfIiIiIiIiuYySLxERERERERtQ8iUiIiIiImIDSr5ERERERERsQMmXiIiIiIiIDSj5EhERERERsQElXyIiIiIiIjag5EtERERERMQGlHyJiIiIiIjYgJIvERERERERG1DyJSIiIiIiYgNKvkRERERERGzAKTsrN5vNHDp0iD179nD27FlSUlKys7ks4+joSEBAACVLlqRkyZK4uLjYOyQREREREcnlDGaz2ZwdFZvNZpYuXcquXbvwzRNAaEQhnF1cwJAdrWUtY2oqF85Ec+7UKUJDQ+nRowdubm72DktERERERHKxbEu+tmzZwurVq6nftg3FypXDYMgFWdd/nD99hhUzZ1IkMpJOnTrZOxwREREREcnFsu2Zrz///JOIEsUpXr58rky8AIJC81Ohdi32799PcnKyvcMREREREZFcLFuSr+TkZM6cOUN4sWLZUb1NhRctSmpqKqdOnbJ3KCIiIiIikotlS/KVlJQEgLuHZ7rnn+/Th/ply913O3/9/jshDo5sWrfuvuu6HXcvLwASExOzrQ0REREREXnwZctuhzcfI7vdcsMXhg0jPi4uO5rOcrl1yaSIiIiIiOQs2brV/O1EREbao1kRERERERG7sctLlm9ddjh3yhRCHBz5c9cuurVoQWEvb6KKFee7adPSXDf2g+GUC8lPpLcPT3TsyMXz59OUMZvNTBg9hlrFSxDu5k71yCJ8PfZTy/mTx49TzM+fd4e8bHVdtxYtqFm0WK6ZkRMRERERkdzFLslXegb06Em9Jk35duH3lKlYgcF9nuDAvn2W85PHf8HIt96iU48eTJo/j/BChXmp35Np6nnz+cGMfPttOj/+ONOXLaVLr14MHzqUqV99BUCB8HDeG/sJEz/9lE2//ALA1AkT+PWHtXw+dQoenuk/pyYiIiIiInI/7LLsMD1PDBhA72efAaBqVBRrl69g+YLvKTbsDYxGI+M+/phOPXvw1qiRADRo1oyLF84zf/oMSx3HDh9m8hdfMGLCl/Ts3x+Auo0bkxAfzyfvvU/P/v1xcHDgsT59WLV4Mc/37sP/FsznvVde5dmXX6ZqVJTtb1xERERERB4KOWbmq17TJpafPTw9CQsPJ/qf7d3PnDrF2TNnaN6undU1rTp2tPr+69q1ALTs2JHU1FTLp07jRpw/e5bTJ09ayo6eOJGE+Hja1KpNoSJFGPLO29l0ZyIiIiIiIjlo5svHz8/qu4uLC0lJN7Z3Px8dDUDeoCCrMoHBwVbfL12MwWw2UzrQutxNZ06epEB4uKWuOo0asmjOXHo8+SQuLi5ZcRsiIiIiIiLpyjHJ150EhYQApNlg48K5c1bf/QP8MRgMLF7/K87pJFNFihe3/PzTqlUsmjOXMhUrMvrdd2nVqWOa5E5ERERERCSr5Jhlh3eSPyyM4JAQVi5aZHV82YIFVt9rN2oEwOWYGCpUqZLm4+XtfeP8pUu81O9J2nd9jO9//gk3d3eGPPWUTe5FREREREQeTrli5svR0ZHnXn2VNwcPJjAomLpNGvPLmh/Y9PM6q3KRxYrR59lnGfh4L54ZMoRK1auRkpLCkQMH2bjuZ6YsXAjAawMGAPDh+PF4+/jw2beT6dy4CXOnTKFL7942vjsREREREXkY5IrkC6DvwOe4duUK3375JVMmTKBu40aM/mYi3Zq3sCr3weefEVm8ONMnTmTs++/j6eVFZPHitOrUCYBFc+aweO53zFyxHD9/fwBqNWhA34EDeXPwC9Rq2JCwggVtfn8iIiIiIvJgM5jNZnNWV3rt2jXGjh1L865dKVi0SFZXb1NJiYlMGTmKzp07U6pUKXuHIyIiIiIiuVSueOZLREREREQkt1PyJSIiIiIiYgPZknwZDAYAzGZTdlRvU2bTjXu4eU8iIiIiIiL3IluSL3d3dwwGA7HXrmVH9TYVe/UqAB4eHnaOREREREREcrN73u0wLi6Oa9euYTKlP7sVEhLCkb37KF2lyj0HlxMc2fc3Li4uODk5ER8fryRMRERERB44SUlJXLlyhdTUVHuHkikGgwFXV1f8/f1xcMj5T1RlerfDffv2sWnzJk6dPJWh8hVq1aJSndo4u7jcU4D2YjKZOLxnDz8vXozZ9G8XFQwvSO1atSlatKgdoxMRERERuX/R0dGsW7eOw4cPYzQa7R3OPfPy8qJ06dI0atQIZ2dne4dzW5lKvnbv3s2iRYvIX6gARcuXwj9vAA6Ojrctv//3PezZugsnZyeCQkNxdnElNzw6lZqaysXosyTGx5O/UEGqNaqNGTOXz1/kwO97OXfyDJ06ddLW8yIiIiKSa505c4Zp06bh4eFDiRKVCAwMw8kp17wGGACz2UxSUgInThxk//6dhIaG0a1b1xybgGU4+UpNTWX06NGEFStEg/aPZHgDimuXrnBk7wFizl4gJTnlvoK1FUcnR3z8/Shcqih58wdb3avJZGLtd8u4ePocL77wQq6Y3hQRERER+a8pU6YSG5tAq1Z9cHFxtXc49y06+jjLlk2hbdu2VKhQwd7hpCvDqe3Ro0dJSkqiQq2qmdr5zyfAjwq1q91TcDmRg4MD5WtVYcnkuZw8eZLw8HB7hyQiIiIikilxcXGcOHGc2rVbPRCJF0BISDj58hVkz569OTb5yvC0zfnz53FxdcE/KE92xpMrBIWFYDAYuHDhgr1DERERERHJtIsXL2I2mwkOLmDvULJUcHCBHD1Gz3DylZqaipOzc7qzXoN7P0PDMjWyNLA7KelXkDHvfJRl9Q3u/QyhBt80n59XrU23vMFgwMnZiZSU3LGMUkRERETkVjd3NXRySvts1Jw5kwkNNXDkyEGr45MnjyM01MDo0W9bHb98+RJhYQ58+eVI5s6dQmiowfIpWdKP1q1rsnr1Ykv5kyePWZUpXNiNunVLMHr02yQkJFjVPWbMO1Zlb/2MH/9xmtidnFxITc25Y/Tc9URdNgovHMG4mZOsjhUtWewOV+SCnUNERERERO4gvYmVKlVqAbBjxyYKF/53h+/t2zfi7u7Bjh2brMrv2LEJs9lM1aq1OXLkAAAzZ67C29uXmJgLTJz4CU880Y6ZM1dRv34zy3VDh35IVFQD4uPj+OGHJYwd+x4XLpxjxIivrOp3c3Pnu+9+ShNnaGjBdO4nEzdvBw918pWQkIC7uzsAbu7uVK5R1c4RiYiIiIjYV5EixcmTJ5AdOzby6KO9LMe3b9/Io4/2Zv78aRiNRhz/2fV8x46NuLm5Ub58FUvyVa5cZQIC8gIQFVWfqlULMHnyOKvkq1CholSufGP1XJ06jTh4cB/z50/jo4++tNrUzsHBwVIut8uWrfrOnDrNwB5PUiZvISLdg+lQtzl/7NxlVWbetNm0q92M0gHhlPIvSKf6Ldm1bWeaulYvXk7dElUo7BZEy2oN+H172jIAa5evplX1hkS6B1M2sDBDn3mB+Lg4y/lN69YTavBl7fLVPNmpJ8V9wniqc6906xIREREReZhVrVqL7ds3Wr6fPn2C6OhT9O37PElJiezb94fl3PbtGylXrgout3mvr5eXN5GRxTl58ugd2yxTpiKJiQnExOTcZ7buV5YnX1cuX6Z97UfY8/uffDBuFBMXTMfD04NHG7bh4vl/O/LUsRN0evwxvpo3lfGz/kf+gmF0rNucwwcOWcr89fsfPNmxJ4WKRvLN9zPo3KsbTz/am+SkJKs2l81fRJ82j1GibCkmLZzJsJHvsfL7pbzU97k08b3a/3kiIgsxaeEMnh4y0HL82KEjlPAtQIRLXh6pXJdVi5ZlddeIiIiIiOQKVavW4uDBfVy5chm4kWDlz1+AyMhilCpV3pKYpaSksHv3dqpVq33buoxGI2fOnCQ4OP8d2zx16jheXt6WGbNbpaampvnkRlm+7HDSpxO4duUqy7f9RN6gQABqN6pHnWKV+Wr05wwb+T4AL7z1quUak8lE3SYN+H3bTr6bMpPXPrzxEN8XH48ltGAYkxfNskxrurm7WSVVZrOZ94e8SZsuHRg9abzleFBIPnq26MTgN1+heOmSluNN2jTnjRHvWcVcpmI5KlStRLHSJbh25SrTJvyPvu278/W8qbTq1C5rO0hEREREJIerWrU2ZrOZnTs306hRC3bs2ESVKlEAVKkSxY4dm+jT5zn+/PM3EhMTqVrVOvkyGo2kpqYSE3OBzz77gHPnonnxxXesyphMJlJTU4mPj2PNmiWsWLGAV18dbhn33xQfH0d4eNqNQRYuXH/HpC8nyvLk65c1PxHVoA5+Af6WjNTR0ZEa9Wrx+/Z/lx4e3Lefj19/jx2btlrNiB25ZeZr19YdNG3TwuoPoGWntlbJ1+EDhzh1/ATvfvqRVQZco14tHBwc2L1jl1Xy1ajlv+tMb+r3/DNW35u2aUGbqCaMfutDJV8iIiIi8tApV64ybm7ubN++kUaNWlie9wKoXLkmH344FLgxI2YwGCyJ2U0VKuSz/Ozm5s7zzw+je/cnrco880wXq+9t2z7Gs8++kiYWNzd3vv/+1zTHixQpcU/3Zk9ZnnxduhjDb1u2E+6c9n1gEZGFAIi9fp2uTduTJzAPb3/yIWHhBXB1c2VIv0EkJf67pPB89DnyBFlPO3r7+ODm5mb5fvliDAB923dPN54zJ09bfQ8MDrzrPTg4ONCyYxs+eOUtq005REREREQeBs7OzlSoUJUdOzYSFxfLvn1/WM18nT59gjNnTrFjx0aKFSuFn5+/1fVz5qzFx8cXX19/wsLCcXJKm3a88cYIatVqyNWrV5gyZTyLF8+hZs369Oz5lFU5BwcHypevkn03a0NZnnz5B/hT6JHGvPz+G2nOubjeeHv2zs3biD51mqnL5lK6fFnL+etXrxIS9u9a0KCQYGLOX7Sq4/q1ayQmJlq++wXc+IMePn40FatXTtNmcP4Qq+/pbacpIiIiIiLWqlatzaRJn7Jt2wZcXFwpXboCAGFh4QQHh7B9+0Z27NhE06Zt01xbunT5dJ/dulXBgoUtSVWtWg1o0aIqI0cOo2PHHnh4eGb5/eQEWb7hRu3G9Tmwdz9FSxanfJVKVp+SZUsDkJhwI3m6dUeU7Zu2cvLYCau6KlSrzA9LV2I0Gi3Hls9fbFWmSIlihISFcvzIsTTtla9SiXz/Sb4ywmQysWzeIoqXLqlZLxERERF5KFWrVpuEhHi+/XY8FSpUtZq9qlw5igULpnP+/FmqVq113205OjryxhsjuXTpIjNmTLzv+nKqLJ/56v/icyycOY+O9VrQ9/lnCC0YRsyFi+zauoPg/CH0f2EAlWpUxdPLi9cHvMRzQ1/g7OloRr/9EflCrXdAeW7oC7So2oAn2nWj17P9OHHkGF+N/txq2aHBYODtT4bzXLd+xMfF0ahlMzw8PTh9/CRrl69m6IdvE1msyG3jPXX8BIN7PUPbrh2JKFKYq5evMG3C/9i9YxffLJie1d0jIiIiIpIrVK5cEwcHB376aQUDBgxNc+6DD14GyLJNL+rWbUy1arX55pux9OnzHM7ONzbZMJlM7Ny5JU35vHmDCA8vnCVt20qWJ18BeQJYumUtI4e9z4evvs3lmEvkCQqkUo0qPNK+NQCBwUF8PW8q7w8ZxhNtu1GoWBFGfP0pX4741KquMhXL8/W8qXw49B36te9O8TIl+XLOt3Rv1sGqXOvO7fH18+Xz4WP4fsZ3ABSIKEj9Rxrd9RkvT28vvH19+OyD0cScv4CziwvlqlRkxsr51G/WOOs6RkREREQkF/H19aN48dLs2/dnmg01qlSJwmw2ky9ffgoWLJRlbb7wwtt07dqE77+fSZcuvQFITEygTZuaacp27dqX0aMnZVnbtmAwm83mjBT85Zdf2LZjO91f6p/dMeUKUz76ggb161OzZtq/CCIiIiIiOdnhw4eZMWMGXbsOxsvL197hZJldu35l375tDBkyxN6hpCvLn/kSERERERGRtJR8iYiIiIiI2ECGky9HR0eMqca7F3wImM1mjKmp6b6vQEREREQkp3N0dAQgNTXVzpFkrdTUVBwdc+4YPcPJV0BAAEmJiVy7fCUbw8kdLp2/iMlkwt/f/+6FRURERERymJvj2JiYaDtHkrUuXozO0WP0DKeFRYoUwcnJib3bdlOjWb3sjClHM5vN7N22GycnJw4cOMChQ4ds2r6DgwM+Pj4UL148R//FEhEREZF7YzabOXfuHIcOHSIuLo4M7o+XaZ6enmzd+gPnzp3EYDBkSxv/5eDgiLe3HxERJfDw8M7Suq9cuciZM0do1qxZltablTK82yHAr7/+ys8//0yJSmUoVqE0foEBlinLB50x1cil8xf5+7e/OPTHPjw83XH3cLv7hVnMZDRx7WosRqORkiVL0qFDBy1/FBEREXlAxMfHM2vmTE6fOYOLizPenp44OGRPYmQ2mzGZTBgw3Ei+bJCAGVONXI2NxWQyUaJEJWrUaHZfiZ/ZbCYxMZ6TJw+ya9eveHi488QTfXB3d8/CqLNOppIvgE2bNrFp82biYmOzK6Yczcvbk2p1KlC+aim7xZCcnMKBPUf4aflGSpQsQaeOnewWi4iIiIhkDZPJxDcTJ3Lt2lXaNq5LZMGwB3KiIyExkV17D/DDhq1ZVqeDgwORkZG0bdsWT0/PLKs3q2U6+YIbfzFOnz7NtWvXMJlM2REXZrOZlatWkjfIn7KVS2RLG5nh4OiAj68XwfkDbTYteze7d+zlp+UbefHFF/H2ztppWxERERGxraNHjzJt2jR6d2xFeGiIvcPJdj9u2s7W3/+iRcuW95Vkurm5ERYWlmNnu251T+vVHBwcKFCgQFbHYiU6OprEhERq1KtEgUL5s7Wt3Kp46Uh+XrGJgwcPUqlSJXuHIyIiIiL3Yf/+/fh4e1Ewfz57h2ITZYtHsmHH73h7exMZGWnvcGwix77nK/afZY3+eR6cN25nNTd3V9w93Lh+/bq9QxERERGR+xQbG0seX58cs8oqu+XxuzHOf5jGsjk2+bq5nNHBMW2Iyxat4psvp9g4ovQdP3YSX0Moi+Yvy9R169dtwtcQym87dluOtazfCV9DaJrPgb9vv6Oig6NDti39FBERERHbMZvNODjk2OF5lrt5rw/TWDZXbpO3fNEqdu34gyef7W3vUMgXEsTazUuILFY4U9eVr1SWtZuXULxkUavjNWpV5YPRb1odKxgRdt9xioiIiIiIfeXK5CsrJSQk3NfDea6urlStUTnT1/n4eKd7na+fzz3VJyIiIiKSk93vuPtBkOvmNZ/pPZhZU+exb89+y7K8Z3oPBmDb5h20atiZEM8iFPAtQd9uA7hw/qLl2ptLBGdOmcvAJ18mIk9pGlZrBYCvIZSxI77gvTc+JjKoHAX9SvLmKx9gNptZ9+N6aldoQn6vorRu9CinTp5OU+etyw7LRlRnyHNv8M0XUygTXo0CviXo1u4JLl6IsZRJb9mhiIiIiEh69uzbT4tHe5InsjQeoZEUr1aHkZ9/CUDvAYMpE9WQlT/8RJmohriFFKZyg0fYsn2nVR3T5syjdvN2BBQujX+hUtRv3YltO3elaWvf/oN0eLwfAYVvtFW+TmNmL1hkOW82mxk97iuKVa2Na75CFK5Yk7FfTrSq452Px+BVoCjbdu6iZtPWuIUU5otJU7O+Y3KZXDfz9fKbg7l4IYYDfx9m0sxxAOQJzMO2zTtoWb8zTVs05Nu5E4iLi+eDYSPp2rYPazcvtarj3dc+pmnLRvxv9peYb1lj+s34b6ldvyZfT/+cnVt38eHbozEZjfz8w3peemMgLi4uvDroTZ7rO4RFa2bfMc6VS9Zw+OBRRn8xnJiLl3jthXd5eeAwvp0z4Y7XbfxlCyGeRTAaTVSpXpE33n+ZWnVr3GNviYiIiMiDoHW33gQH5uV/n4/B18ebQ0eOcepMtOV89LlzPPvy67zz6ov4+/nx8afjadapOwd3bCAoMC8Ax06c4vEunYgsFE5ycgqzv19E3VYd+WP9DxQrcmO3wYOHj1CzWRsKhObn84/fI19QIH/t28+JU/9OPjz/2ltMmj6LN14cRPXKFdm0bSevvvsh7u5uPN3ncUu55OQUuvV/jheeeZIP3xxKngB/G/VWzpXrkq/CkRHkDczDyeOnrZbnPffES1SsUo4Z30+y7BBTumxJapRpyJoVP9K0RSNL2bIVSjN+0ug0defLn4+J028kdI2b1WfFkjV8MfYbtu752fJsVvTps7w8cBhXrlzFz+/2OzGazWbmLPkWV1dXAE4cO8WYD8dhMplu+yBlrXo1eOzxTkQWLUT0mbOMG/01bRs/xopf5lOtZpVM9pSIiIiIPAguxlzi6PETfPbRu7R+pCkADerUsipz6fIV5n37NQ3r1gagXq0aFChTlbETvuGjt14D4K1XXrCUN5lMNGlQl22//c6U2d/x4Zs3yrwz4hNcXJzZuHIRPj433iPbuH5dy3WHjx5j/Dff8tWYj+nfu4flfHxCAu+OHEv/Xj0sY92UlBSGv/EKXTq0zY5uyZVy3bLD9MTHJ7Bl43badW6F0WgkNTWV1NRUihQrTFiB/Py23XppX7OWjdKtp0GTOlbfixQrTEj+YKtNMW5urHHmVDR3UqteTUviBVC8VFFSUlKslkH+1+vvDqHnE48RVac6Hbu0Zfm6+YTkD2bk+5/dsS0REREReXDlCfAnvEAYr733MVNnf8ep02fSlPH18bEkXje/N65fh607f7Mc27f/IO179iW4eHkc8xbAOSic/QcPc+DQEUuZH3/dQKc2LS2J13+t/WU9AB3btLCMuVNTU2lcrzZnz53n5H9ia9m08X3d+4PmgUi+rly+gtFo5LUX3iGPc7jV5+SJ05w6af2XIDA4b7r1+P5nJsvFxTndYwCJiUl3jMnXz+c/17lk6LpbeXp60LRlI37f+UeGrxERERGRB4vBYGDNglmULFaEAa+8QYGyVanSsDm/btpiKROYNyDNdcGBeYk+ex6A69djadqxK8dPnuKTD95m/YqFbP9xBeXLlCIx6d/xacyly+TPF3zbWC7GXMJsNpO3SFmcg8ItnyYdugJYJV8eHu54eXne9/0/SHLdssP0+Pr5YjAYeOn1gbRs90ia83n+85fxYXlxnYiIiIg8GIoViWTelImkpKSwadsOXn//Y1p37c3pPTc21bhw8VKaa85duEhIviAANm/fyakz0SybM5XyZUpbyly9dp2w/CGW73kC/Dlz9txt4wjw98NgMLBhxULL5MKtiv/z7BhozJ2eXJl8Obu4WM0geXp6UK1mZfbvO8SbH5S3Y2RZKy4untXL1lKp6oNzTyIiIiJy75ydnalXqyZDBw+gTbc+nDl7FoCr167x068bLEsPr167xtp16xnQrzcACYmJALg4/5swbdq6nWMnTlK6RDHLscb16jB/yXJGvP0G3t5eadpv9E/9MZcvW54/k4zLlclX8ZJFmDF5DvNnL6Jw0ULkyRvA+6OG0bphF3p3eZqOj7XFz9+XM6ei+fmHX+nepwt16kfZO+w72rR+K5+PmkCr9s0pGBHG2TPnGDfma86dvcDUeV/bOzwRERERsZM/9uzlpWHv0aV9GyILhXP12nU+GjueiIIFiCwUAdyYkeo7aAjvDn0JP19fPv50PGbMDH66HwA1qlTCy8uTAa+8ztDnn+N09Fne/ng0oSH5rNp6+5UXWLZ6LbVbtOOVQc8SEhzE3v0HiU9I4JVBz1KsSCQD+vWm59PP8/LAp6leuSIpKakcOHyEnzdsYtGMybbunlwlVyZfPft2Zee233l54DAuxVymW6/OTJjyKas2LOSjt0czoM+LJCcnkz8shHqNalO4SIS9Q76r4JAgkpNTeO/1j7kUcxkPTw+qR1Xm068+pnK1ivYOT0RERETsJF9QEPmCA/no0/Gcjj6Lr483dWpUY8bXn+Po6AhASHAwI955nZff+oDDx45TukQxVs+bSXBQIADBQYHMm/w1Q956n7Y9nqBYZCG+/mQEI/55V9hNRSMLs2nVYl57/yOeHfI6qcZUikUWZujzAyxlPv/4fYoXieTrKTN4b9SneHl6ULxIJJ3btrJdp+RSBrPZbLZ3EOnZv38/c+bM4akhPfDwfLjfhH0n34ydRaUKlWnQoIG9QxERERGR+zBv3jySYq/Ro13zTF3Xe8Bgduz6g782/ZRNkWUPs9nMe+Mm0bp1aypVqmTvcGzigdjtUEREREREJKfL8clXDp2YyznUPSIiIiIPEA3uHmQ59pkvy3uxEpLw9PKwczQ5k9lsJikxyeplziIiIiKSO7m4uHAlMTnT10354tOsD8YGEv7ZvfxhGsvm2JmvkJAQHBwcOH74tL1DybHOnDxHSkoqYWFh9g5FRERERO5TaGgo0RcuEhefYO9QbOLwiVPAjft+WNhlw41Lly6xe/duDhw4QGxsLCZT+iEkJydjMptwcXF+oF/SZjCAi6sLBQqFUKJsEcLCQ+56TXJyCgtnrOT6lXheeOGFB7p/RERERB4GcXFxjBkzhjLFImnXpB4ODjl2nuS+xcbH8+28pVyNjcPVzS3b2zMAbm7uFC5ciHLlytlt8sLmydepU6eYMWMGZjMUjCiOj68/Dg6OtgwhxzFjJiE+luNH9xMXe43iZSIpUCh/umVNRhOXLlzm4L5jpKYY6dGjh2a+RERERB4Qf/31F99//z3+vj6UjIzAz8cLB8ODk4SlGo1EX7jIvkPHMGOgWNmKuLln/yNGZrOJ+NjrnDh8gLjY67Rr147y5ctne7v/ZdPky2g0MmbMGHz88vJIq+64uDw86zszwmw2s3XjGv7cveW2ZRwcHPD28aZE8RJUqlSJoKAgG0YoIiIiItnt+PHj/P777xw8eJD4+PgHagM6B0dHfHz8iChWkhLlKuPt62fT9s1mM+vXLGX/n7t47rnnyJMnj03bt+mGG4cPHyYhIYEWbVsq8UqHwWCgWlQTDh34kwoVytG0aVN7hyQiIiIiNhYeHk54ePg9XWs2mxk3fjyB+QtS75G2WRxZ7mcwGIhq2Jwjf+9hz5491K1b16bt23QO88SJE3h6+eAfoNma23FwcCC0YCTHj5+wdygiIiIiksvExcVx+dIlChYuau9QciwnZ2dCCkbYZbxt0+QrKSkJNzePNJtDbNq4jtBAwx0/Y0a+k+F2Tp86weuvDCCqaiSFw9woFuHNI40qM3rE28RcvJDFd5X13N08SEpKsncYIiIiIpLL3BxDurql/xzVqqWLmDLxS1uGdFsnjx8j1MvAsoXzM3Xdpl/XEeplYPdvOyzH3nt9CA2qlKZYPm+Kh/jQom5VFs+bc9s63NztM962+Xu+0tuVr2y5SixZuTnd8qM+fpPNG9dRr0GzDNX/286t9HysOX5+AfR98nlKlCpLakoKO7ZvYvqUrzhy+ABfTpx9X/eQ3W700YOztldEREREbOt2O2GvWraIP37bQe/+z9o4orSC8oWw5KfNFC5SLFPXla1QiSU/baZo8ZKWY3GxsXTr8yRFipXAYDCwfNF8nu3TFZPZRPtHu6Wpw17j7RzxkmVvbx8qV6mR5viaVUtY/8taXh76PlWq1rxrPYmJiTzVtzMh+cNYuGwD3t4+lnP1GjTlqWdf4ofVS+9YR0JCAu7u7pm/CRERERGRh8j9jptdXV2pXC1tDnA33j4+aa4b8flXVt/rN27Ggb/38t2MKekmX/aSY/etPHv2DC8+/wQ1a9Vn0AuvA5CSksL777xM1QoFKRTqSsXSIfTq3ppr164CsGzJPM6cPslrwz62Srxu8vLypn3Hfzt/7uwphAYa2LF9M491akKRcE8+eOdlAM6cOcXAZ3pQpnheIgu406F1Xf7YvTNNnXNnT6FxvXIUDnOjctlQPh7+BkajMU0bf/2xix5dmlMk3JNa1Yoyb+60LO0vEREREZE7GfxUb+bNnMr+fXsI9TIQ6mVg8FO9AdixdTOdWzSkSJAnJfL7MqBPNy6eP2+59uYSwbkzpvDyc09SumAeWtWrBkCol4EvPhnBx++8QbmIIEqG+vHBsFdu7Cz48480qVmBosFePNqyEadPnUxT563LDquXiuCNF59jytdfUK1kOCXy+/LEY+2IufDvo0PpLTtMj39AHlKSk7Oi67JMjpj5+i+TycTAZ3oAMP6rmZYXzI3/7COmT/mK198aQfHipbl06SK/rFtD8j/rNTdvXIeTkxO16jTMVHvPPd2N7j37M3Dw67i7e3DlymXat6qNp6cXH3w0Dm9vX76dNI5H2zdkw9aD5A28sWHI1xM+Yfi7r/Dk0y/w1rtjOHhgHyM+fAOT0cjrb31s3cYz3enW80n6P/MiM6d/wwsDe1OhYlWKFiuZXkgiIiIiIllq8KtvEnPxAocP/M24/80EIE/ewBuJV/P6NGzagglT5xIfF8fI94fR57G2LP3J+tGgj99+jUaPtOTLb2djMpksx7/9ejw169Tn80nT2bV9K6OHv43RaGT9Tz8w8OU3cHFx4c2XBzHk2b7MXrLmjnGuWbGEo4cPMvyTL7gUc5F3h77AsCEDmTD19s9wwY2dHo1GI3Gxsfywcim//riGzyfNuMfeyh45Mvn64vMRbNrwM9/OWEK+fP++bHjXb9uoV78pvZ/4d41qy9YdLT+fO3uGgIC8uP3nLdlGo9HyfgSDwYCjo/VLnXv2epoBg161fB894m2uXb3C8tXbLIlW7bqNqFOjGF99OZphb48kNvY6Y0a8zTPPvcJrwz4EoG79Jji7uPDuWy/y9HMvExDw73sDevd9zhJ3lapR/Lh2OcuXLmDwS8Puq69ERERERDIionAkefIGcvrEcatley898wTlKlVh0uzvLc+KlSxdlobVyvDj6hU0atbCUrZ0uQqM/mJSmrrzheRn3KTpwI0lf2tWLOGb8WP5efseipa4Mdlw9sxphg0ZyNUrV/D187ttnGazmW+/W4Kr641XU506foxxoz/EZDJZJmXSs37dj3Rt3QQAJycnPhgznlbtO2Wwd2wjxy07/G3nVkaPeIsn+g2kabPWVufKlqvETz+uYMzId/h913arbPum9B4uLFHYl/AQZ8JDnClR2DfN+UZNWlp9/2XdGqJqN8DPP4DU1FRSU1NxdHSkRlQ9ft+1HYAd2zYRFxdL6zadLWVSU1OpU7cxiQkJ7N/3l1Wd9er/+84uD09PwsLCiY4+lfGOERERERHJYgnx8WzfspFW7TtjNBotY9rCRYuRP6wAu3dutyrfqFnLdOup06CJ1ffCRYoRHJLfkngBFC56Y2ON6NN3HgPXrF3PkngBFC1RipSUFC5eOH+Hq6BSleqs+HU7c5aupd+Awbw5ZCCzp/7vjtfYWo6a+bp+/RoDnupKseKlGfbOqDTnB73wBg4ODsybO5VPRr1LnryB9H5iAC8MeQuDwUBwvvys/3UtSUlJVn9gC5eux2gyMnPaRBYumJWm3sDAYKvvl2Iu8tuOLYSHOKcpGxEReaPMpYsANGtUKd17OXPmpNV33/+8vdvZ2YWkxMR0rxURERERsYUrVy5jNBp559UXeOfVF9KcP3PaekybNyg4TRkgzUyWs4tLuuNfgKSkO4+Bff5znYvLP9fdZezs5e1N+UpVAKjToBGpqam8+9qLPNqjd5qVb/aSo5Kv115+hgsXzjF99gqr5OkmV1dXXnrlHV565R2OHjnEnFmTGTPyHQqGF6bToz2pWas+c2ZNZtOGn2nQ6BHLdWXKVQRg7Zpl6bb739kyf/8AChV6hJdfez9NWReXG3H5+QUAMGnK9+QPLZCmXIGChTJ41yIiIiIi9uHr64fBYGDgkNd5pHW7NOcD8uS1+n67LexzonIVKzPpi0+JuXiBoOB89g4HyEHJ17y501i4YBZjPv0fRYqWuGv5QoWL8NqwD5kx7WsOHtgHQKs2nRnx4Rt89MFrVK1eCy8v73uKpXbdxnw/fwZFi5bEw9Mz3TKVq9bE3cOD6DOnaN6y/T21IyIiIiJiSy7OLlYzTx6enlSuXpND+/dR/u0P7BhZ1tu2aQPePj5pEkh7yhHJ17Gjh3nj1QGUr1CFosVLsXPHljRlvL19+Hj465QtX5kyZSvi4eHJD6uXcvXKZcvuhm5ubnz9v3n06PIIzRpW4ol+AylRqixGo5GjRw6yZNHcDCVk/Z95kYULZtKxbT369n+e0NCCxMRcYNdvWwnOl5/+T7+Ar68fQ159j+HvvUJ09ClqRtXH0dGR48ePsGblYr75dgHuHum/WVxERERExB6KlCjJnOmTWfTdbAoVKUpAnrwM+2AUXVo25OnHu9C202P4+vkTfeYUv/70A1169CGqbn17h31He//6gw/ffJVW7TsTVjCC+LhY1q5cxqypk3jt3Y9wcsoRKQ+QQ5KvrVvWExcXy+7fd9CmefovU64ZVY9GTVqydPF3TPxyDKmpqUQWKc74r2ZSt15jS7lKlavzw7rdjP/8YyZN/JSz0adxcnKmcGQxWrXpTJ++z901noCAPCxduYWRHw3jw/de5fLlGPLkDaJS5Ro80uLfWa6nn32JkJBQJk74hMmTxuHs5Ex4RCSNm7bC+Z+1qSIiIiIiOUXXx/vy+45tDBsykMuXYujcvReffj2FhT9sYPTwt3nxmT4kJycTEhpG7XqNiIgsYu+Q7yowMBgfXz/GfvweF86dxdvHlyLFSvC/2Qtp1qqtvcOzYjDf3IPdBpYvX86x46dp/2h/WzWZK23duIbTJw/w3HN3TxRFRERERG6KiYlh/PjxtOrSm5AC4fYOJ8f6dfUSYi+dp1+/fjZtN8dtNS8iIiIiIvIgUvIlIiIiIiJiAzZNvgwGA+Z0Xows1kxmEwaD8mIRERERyZybW8GbzBpz34nZZLbLeNumLXp6ehIbexWTErA7un71Cp632eJeREREROR2PP7ZbTv26hX7BpLDXb962S7j7dtuuJGSksK+ffvYu3cvMTEXSU013ndjRqORuLg43Nw9cHTMERst5jxmM/Hx13F1dbW80BnAwcGAh4cnRYsWpWzZsvj7+9sxSBERERGxp8TERPbu3cu+vXu5fOUKRuO/Y/W4uDgwGHBz12uP0mM2m4iPjcXdzS3Ldig3GAx4enhQ5J+xekBAQPrl0ku+kpOTmTFzBidPnCQsLJjQ0ECcnJUs2ZPJZOLq1VgOHTwJGOjevTsFCxa0d1giIiIiYmNxcXFMmzaNCxcuEF4ghHzBeXBydLR3WA81k8nE1WuxHDxyErMZunXrRkRERJpy6SZfixcvZu/ePXTr0ZwCBYJtEa9kUFJSMnPnrCH6TAwvvfQSLnqfmIiIiMhDZebMmURHn+bxLi0ICkx/hkXsIzk5hbkLf+Dk6fO8+OKLuLm5WZ1P88xXamoq+/bto0bNskq8ciBXVxdat6lLcnIyBw4csHc4IiIiImJD8fHxHD58mPq1KyvxyoFcXJxp26IeKSkp7N+/P835NMnX2bNnSUpKomgxLWnLqfz9fcgb6M+xY8fsHYqIiIiI2NCJEycwm80UK1zA3qHIbfh4exISnDfdsXqa5CspKQkAT0/3bA9M7p2Xp7vlz0pEREREHg6JiYmAxuo5naeHm+XP6lZpkq+bj4D984oAmyocUZ2Bz72RLXUvWrSKCV9OyfR1774zhk2btmd9QLfo03swjobQNJ9Vq36+7TUGe/wBiYiIiEiOYI+xYETx6jw3OJvG6ktW8eXXUzJ93TsfjGHT5uwdq4/6ZAIVazTFL19JPPMUoWyVRoyf8C232TQeuP2fT47awnDBwv/h7++bLXUvXrSKnTv+4Jlne2fquvfe/QQvL0+ioqpmS1w3FS4czvSZ46yOlSxZNFvbFBERERHJqIVzs2+svmjpKnb89gfPPtU7U9e9O/wTvDw9iaqZfWP1K1ev0qVjG8qULo6bmys//ryRQS+9ybXr13n9lUGZqitHJF8JCQm4u7tTsWIZe4diN+7ubtSoUdneYYiIiIiIWLGM1Ss8nGP14e8OtfreuGFdTpw8zZTp32U6+Uqz7PB2pkyZi4tTQc6du2B1/NKly7i5RPD119PZvHkHbdv0Jix/Jbw9i1CpQhOmT59vVX7duk04GkJZvnwtnTs9iZ9PcR7t/BSQdtlhZur74Ydf6d5tAL7exSgUXo1RI7+0lOnTezDTps5jz579liV9fXoPvus9OxpCAXjl5fct161btwmAT8Z8RfWqLfD3LUG+oHK0bvU4Bw4cTlPH119Pp1B4Nbw8Imna5DF27foLR0MoU6bMvWv7IiIiIiIZMWX6XJy80h+ru/hE8PWk6WzesoM2nXqTv1AlPPMUoUL1Jkyf9Z+x9a+bMLiHsnzlWjp1fRKfoOJ07n5jrP7fZYeZqe+HH3+lW68BeAcWI7xYNUaO+Xes3vvJwUydMY89e/djcA/F4B5K7ycH3/WeDe43xuovv/6+5bp1v94Yq4/59Cuq1mqBb3AJggqWo1WHxzlwMJ2x+qTphBerhkdAJE1aPsau3//C4B7KlOl3HqvnCfAnOTnlrjH+V4Znvtq3b86zT7/G/HnLGPBcH8vxBQtWANC5cyvWrPmVWrWq8tTTPXFzc2XTxu082XcIJpOJXr0etarv6f6v0r1HB55e+DiOt3kp3PHjpzNc37NPD6VHz44sWDiJxYtWM/TV4ZQtV5JHHmnAsDcHc+FCDPv/PmxZ2hcYmOeu97xx8xJq1WzDcwOfoGu3dgCUKlUMgFOnonn2ud6Eh4dx7VosX381ndpRbfn7wHoCAvwBWLJkDc8+PZS+/brRsVNLdv++h8cefSrdtg4dOoa/bwkSEhIpW7YEb7w5mHbtHrlrjCIiIiIi7ds05+mBrzHv+2U898wtY/VF/4zVO7RizdpfqVWjKk/3uzG23rh5O32f/mds3cN6bN1/wKv06NqBhf3vMFY/cTrD9T09cCg9u3Vk4dxJLFqymleHDadc2ZI80rQBb742mAsXY/h7/2FmfpvxsfrmdUuoWb8NA595gm5d2gFQquQ/Y/XT0Tz3TG/CC94Yq381aTpRDdpy4I9bxurL1vD0wKH069ONTu1b8vvuPTzaI/2xOtx4JVdCQiK/btjCtFnzefv1F+4a439lOPny9fWheYuGzJm9yCr5mjt7EU2a1iUgwJ/HHmtrOW42m6lbtwanTkXzzdcz0iRLrds04eMRd35gLzP1dejYgrffeQmARo3qsGL5jyyYv5xHHmlAZGQEgYF5OHH8dKaW9t0sW7BgaJrrPhn7ruVno9FIkyZ1yBdUnvnzl9O/fw8APvzgMxo2rMXEb0YB0KxZfVJSUnjrzVFWdVWsWIaqVStQqnQxrly5xlcTptGxfV/mzvuaTp1aZTheEREREXk4+fr60KJZQ2Z/t8gq+Zr93SKaNv5nrP7of8bWtWtw6nQ0X0+akSZZatOqCSOG32Wsnon6OrZrwTvD/hmrN6jD8lU/Mn/hch5p2oDIwhEE5s3D8ROnqVE9E2P1f8oWLBCa5rqxo/4zVm9Uh6CC5Zm/cDn9+94Yq3/w8Wc0rF+Lb778Z6zepD4pqSm8+a71WB3g0OGjFC1T2/J92NDneWFQ/wzHelOGlx0CPNa1LZs37+TEidMAREef45dftvBY13YAXL58hecHvUmh8Gq4Oofj6hzONxNncuDAkTR1tWjZ6K7tZaa+Jk3rWn42GAyUKFmUU6eiM3N7mbJly06aNnmMwDylcXEqiJdHEWJj4zj4T2xGo5Fdu/6idZumVte1adssTV2Dnu/HswN6U79+FO3aPcKKlTOoXr0i77w1OtviFxEREZEHS9dH27J563/G6uu30LVzO+DG2HrQi28SXqwazt7hOHuHM/F/MzlwKO3YuuUjGRurZ7S+po2tx+olSxTl1OlsHKtv3UmTlo+RJ7Q0Tl4F8Qi4MVY/cPCWsfruv2jT0nqs3rZV2rE6QIGw/GzfsIKfV8/jvbeG8MnnE3n7/cyP1TOVfLVq1RhPTw/mzlkMwLzvluLm5mpZHten9wvMmb2Il4Y8zao1s9i6fQV9nniMxMS076MKDg68a3uZqc/Pz3rnFRcXZ5LSKZcVTpw4zSNNu2EyGpnw9QjWb1zE1u0rCArKa4ntwoUYUlNTyfufKdOgoLx3rd/BwYEOHVuyb99BEhISsuUeREREROTB0qrFjbH6nHk3xurfLfhnrN7mxli9d/8XmD1vEUMGP82aZbPYvmEFT/S6zVg96O5j9czU5+f7n7G6s3O65bLCiROnadq6G0ajka/HjWDjT4vYviH9sXpg3v+M1QPTH6u7urpSpXJ56teN4s3XXuDDd4cyfMTnnD17PlOxZWq3Q3d3d9q2a8bcOYt5+ZVnmTtnMa1aN8HT04PExESWL1vLmE/e5rmBT1iuMZlM6dZ1t3cTZLY+W1q16mdiY+OY//0kS9KXmprKpUtXLGUCA/Pg5OTExQsxVteeP3/RlqGKiIiIyEPC3d2ddq2bMWf+Yl556VnmzFtM6xb/jtWXrVjLJyPeZuCzt4ytv773sXpm6rOlVT/cGKt/P+fuY/ULF/8zVr+QsbF65YrlMBqNHDt+knz5gjIcW6ZmvgAe69qOXbv+YvXqdWzZ8ptlyWFSUjImkwlnF2dL2evXY1m6ZE1mm8iW+lxcXO4pu3Z2dk7zdurEhEQMBgPOzv/G9t13S0lNTbV8d3R0pGLFMixZvNrq2sWLVt21TZPJxPx5yyhdujju7np7uYiIiIhkTNdH27Hr979Y/cM6tmz7ja6PtgP+HVu7/GdsvWT5/Y3Vs6q++xqrJ1mP1RPSG6vPT2esXr4Mi5dZj9UXLb37WB1gw6ZtGAwGCkUUzFS8mX7PV5MmdcmTx59+T7yEn58vzZs3AG485Fe1agVGfvyFJZMc+fF4fH197mm2J6vrK1myCN9OnsPs2YsoWrQQefMGEBFRIEPXLVm8htp1quPp6UHx4pE0aFgLgCf6vED/p3qwd88BPhnzdZqlj68Pe572bfvQ/8mX6dS5Fb/v+otpU29sv+ngcCPvPX78FH16DaZL17YUKRLB5ctX+XrCNHbs2M28Bd9k+j5FRERE5OHVpNGNsfoTT/0zVm92y1i9cgU+Hv0FgXlvjK0/Hj0eXx+fDM/23Cqr6ytZvAiTp85h9txFFC3yz1g9PANj9RJFWLx0DXVqVcfTw4PixSJpWP/GWL1P/xd4ql8P9uw9wJjP0o7Vhw19nrad+/Dksy/TuUMrdv3+F1NnWI/Vr169Rot2PenRtQNFIguRkpLCul8389kX/+Opfj0y9CjVrTI98+Xs7EzHTi05c+YsHTq2wMXFxXJuxqzxFCkSQZ9egxk86E06dmpFz8c7ZbaJbKnvib5d6dS5Fc8PHEb1qi14950xGbpu3BcfYjKZaNm8B9WrtmDnzj8oW7Ykk6eM5bedf9KmVW/mzF7Ed/Mn4uvrbXVtmzZN+WLCR6xZvY72bZ9g1cqf+WLCR8CNv7AA3t6e+Ph68+EHn9GqxeP07fMiJpOJ5Stn0L5983u6VxERERF5ODk7O9OpfUvORJ+lYzvrsfqsqeMpEhlBrycHM+ilN+nUvhWPd7/3sXpW1te3d1c6d2jFwJeGUbV2C975IGNj9S/GfojJbKJ52x5Urd2Cnbv+oGyZkkz5Ziw7d/1Jqw69mf3dIubPmoivz3/G6q2aMuHzj1j9wzradn6ClWt+ZsLn/4zVfW6M1d3cXClWtDCffD6Rtp370LPvIH5Zv4Wvxn3M+LHDM32fBrPZbL71wKFDh5g5cybPv9AVHx+vTFcod/a//82mf78hHD66JUMzb7czY9oKPD396dixYxZGJyIiIiI52e+//87ixYt58+W+ltkZyTr/mzKbfs8M4ejfWzI083Y7s+atwtHFky5dulgdz/SyQ8m4S5cu8967Y2nQsBbe3p7s2L6bD4d/Tpu2ze4r8RIRERERkftz6dJl3h0+lob1b4zVt+/YzfCRn9O2VbP7Srzu5KFOvm596O6/DAbDbd/mnVHOzs4cPnyM2bMWcuXKNQID89CjZ8e7vlxaRERERORhZ5Ox+tFjzPrun7F63jz07NaRER9k31g9TfJ1c1tJ68WID55jx04SWajGbc/Xq1eTn9bNv682vL29WLps2n3VcTvmB/0PSERERERu60EfCx47fpJCJe4wVq9Tk3Vr7n+svux7247V0yRfrq6uAMTHJeDr++A+85U/fzBbt6+47Xlvb08bRpN5cfGJ5Mnjau8wRERERMSG3NzcAIiPT8zx49X7kT8kmO0bcvFYPSGRYN+0L2xOk3zly5cPFxcXDh48SUj+zG2dmJu4uLhQpUp5e4dxT65cuc6F85eILFyM1NRUnJwe6tWjIiIiIg+NAgUKYDAYOHD4JJUrlLB3ONnGxcWFKpVz51j9+vU4os9epFr1WmnOpdkixcnJiRIlS7Bl85+cOX3BJgFKxqWkpLJ08a8YDAa2bNnCmDFjOHXqlL3DEhEREREb8PT0pHChQvyy8TdiLl2xdzjyH6mpqSxZtR4nJyeKFy+e5nyareYBkpKSmD5jOqdPnaZgeAihoYE4O2t2xZ5MJjNXrlzn4IETGI1mmjXrhoeHF7/+uoQrV87z9NNP4+fnZ+8wRURERCSbxcbGMnXqVC5dukShgiHky5cXp/vcfELuj8lk5uq16xw4fJLUVCOPPfYYkZGRacqlm3wBJCcns2fPHvbu3UtMTMwddxuxtatXr+Dr65fm59uJiYnBZDISGBh01/ruJCkpEaPJhJurRyYjvn8GgwE3N08KFChKkSJl8fO7sYY0OTmJmTPHULduHerUqWPzuERERETE9hISEvjrr7/Yt28fly9fxmg03ld9GR0P36uLFy+SN2/aZ6AAYmIuYjKZbjtWzw0MBgOenh4UKVKUsmXLEhiY/uNbt02+crI2bdqwZMmSND/fzr87OKZ/qxmpw2w289FHH1GuXG0qVap7D1Fnn7Vr55GSco2+ffvaOxQRERERyYUyMh6+H/ny5ePs2bPpnrvbWP1BotdiZ1BKSgopKSn4+PjbO5Q0vL39iI2Ns3cYIiIiIiJyB0q+MuhmJu7gYN1loaGGu37mzp3CoUP7eeON56hfvxSRkR7UqFGIoUOf4dKli7dts3Hj8oSGGti6df0dY3NwcMRsNt3/TYqIiIiISLbRLhr3acmSzVbf27SpyRNPDKRdu26WYxERkSxZMpetW9fTo8dTlCpVnlOnjjN69Fts3ryONWt+t7xf7ab9+/ewb98fACxcOIvq1fU8l4iIiIhIbqbk6z5Vrpz2zduhoQXTHG/btiu9ew+wrGkFKFSoKO3a1WLt2mW0bNnRqvz338/EwcGBGjXqsWzZPN5//3OcnZ2z5yZERERERCTbPVTLDu35EF9AQB6rxAugTJmKAJw7d8bquNlsZvHi2dSq1ZD+/V/k8uUYfv55lc1iFRERERGRrPdQJV85abt8gO3bNwBQpEhJq+M7dmzi5MljtG/fjfr1m+Hvn4dFi2bZI0QREREREckiD1XylZiYaO8QLBITE3nvvSGUKVOROnUaWZ1buHAWbm5uNG/eAWdnZ1q27MSaNUuIi4u1U7QiIiIiItkrJSXF3iFkOyVfdjJ06NOcPHmUzz6bZrUcMTU1lWXL5tGwYQt8fHwBaN++GwkJ8axcudBe4YqIiIiIZKucNFbPLg9V8nW/b/7OKiNGDGPhwpl8/fU8SpQoY3Xul1/WEBNzgcaNW3P16hWuXr1CiRJlCQ4O0dJDEREREXlg5ZSxenZ6qHY7dHNzs3cITJ48jnHjPmTs2CnUr98szfmbCdaLL/bhxRf7WJ2LibnAxYvnyZs3yCaxioiIiIjYiru7u71DyHYPVfJl7z/QRYtm89Zbz/Paax/RufPjac4nJMSzevViHnmkHX37Pm917sKFszz7bFeWLJnLE08MtFXIIiIiIiI24eLiYu8Qst1DlXzZ8w908+ZfGDy4F7VqNaRGjXrs3LnFci4kJIz8+cNYvXoxcXGxPPHEIKKi6qep48svR7Jw4SwlXyIiIiLywPnva5keRA/VM1/2/APdtOlnUlJS2LDhR9q0qWn1mT17EnBjl8PQ0ILpJl4AnTv34rfftnDs2GEbRi4iIiIiIlnhoZr5soXTp9N/kfNLL73DSy+9c8drp05desfz/fo9T79+z9+xjIiIiIiI5EwP1cyXiIiIiIiIvSj5yiSzOf2ZLfvKiTGJiIiIiMitlHxlkLOzMwBJSQl2jiStxMR4XF1d7R2GiIiIiIjcgZKvDHJwcCAsrADHjx+wdyhWzGYzp04dpmDBgvYORURERERE7uCh2HBjxYoVxMfH33c9ZcqUZtWqVRw5sofChUtnQWT3x2w289tvvxAbe5UyZcrYOxwRERERkUxbvHixvUOwmYci+WrevHmW1FO1alVmzZrFjz/O588/t5A/fyFcXd2ypO7MMJvNJCbGceLEQa5cuUjDhg0JDw+3eRwiIiIiIverTZs29g7BZnJV8vXkk08SGRkJwLhx4/j777+Jj49n0KBB9zWztXv3bvr162f5XrZsWQYOHIiDg/WqTAcHB44ePUrfvn3Zu3cvhw//TlxcHCaT6Z7aTU014uTkiMFgwNnZhYy/hsyAm5sr4eERtGzZjCJFitxT+yIiIiIi/zVq1Cj279+fpXVevXrVarydnvDwcF599VVcXFyytO2cxGDOmdv3pevmS5IbN27ML7/+SkpyMkWKFOHQoUM4epfD4OB+T/UajUYcHR3/+WYi9doO2rZpx8yZ0/H09LQq26ZNG5YsWWL5Hhycn4sX43Fw8L2HdlNxdHTCbD5H8eJFWbFimWawRERERMTmbh3jOru4YPINxcEnOMvqT01Jxcn5zvM+phO7qFmzBou+X0DevHmzrO2cJFfNfDm5u5CakExMTAzuXl6kXLpEXFwcAMbC88C92D3XnXrrl0tLWLq8GzVq1mbF8iUUKFDgDleaMZkqYjLVv7d2UwHOsX//XCpVqsKyZUuoWbPmPdUlIiIiIpIVTI1ewVTnmSytM/VuBY5sYsukdlSqUo2Vy5dSurT991jIarlqt0P3QF8Cyhbk4sWLlmOxCdmw9XtAG0wlN7Hv0CUqVqrK1q1bs74NK8EYjX25csWTunXrMX369GxuT0REREQkhykchXHIds6kelGtRk1WrFhh74iyXK5KvgDyNChOzJVLN74E+FtmvrKcZzmMJbdxOSWSOnXrMWvWrOxp598GMZl6kppamscff5yhQ4fe87NkIiIiIiK5UkA4xhc2klC4Aa1at2bs2LHkoqek7irXJV9BDUsSf/3GJheO5cpiSknJvsZcgjEV/4kU3y50796dN94Yls1/+E5AW6ApI0aMpF279sTGxmZjeyIiIiIiOYybN+Z+CzE3epkXX3yRfv2eJDk52d5RZYlcl3wF1i0OBgNGYyoO5ctmf4MOrlB4ChQcwYcffcj2HTuzb7YNAAMQBTzG8uWrqV69JsePH8/G9kREREREchgHB2j7MfScypRp02nYuInVo0e5Va5Lvlz8PfEvX4DUVCOGvHlwCArK/kYNBsj/ChRdxLnzV6hRszYnT57M5kaLYzI9wf79Z6hUqQqbN2/O5vZERERERHKY6o9jGvQzW3bvo1KVauzZs8feEd2XXJd8AQQ1LkXqjW0CMVQqTyZekHV/AtpA6S3aiENERERExFb+sxHH8uXL7R3RPcudyVfDUpjNNzajcChfDhwdwVYP4nmWw1hqu2UjjsTEpOxuUBtxiIiIiMjD7ZaNOFq3bs0nn3xi74juSa5MvvLWLsqNZ6PAsWzpGy/LSrXhGlDnoBsbcXi14tq1qzZo8OZGHA0YMWIEb7zxhg3aFBERERHJQdy8Mff7HnOdAbz00kv2juae5KqXLN+UdDEWuDHTZb5y5cZBRy/bBhH3O47xGzE7OGKbiajrODoewNnZnQYNGtiiQRERERGRnOXCIRz3r8bLP8DekdyTXDnzdf6nvdyc+TLt/hOcncDB3XYBXJyFw9/1qFK+MAEBtviDP42j4ySCgmDz5k00bdrUBm2KiIiIiOQg+37A8ZMaRPo789uO7faO5p7kzuTr579xcPwn9N1/QkqqbRo2m+DkMDjUne7duvDLLz/h4JDdm338hYPDFCpWLM5vv+2gQoUK2dyeiIiIiEgO88sXGCY0p1GdmmzbsonChQvbO6J7kuuSL7PZzKWf9+Pk6ARGI8Y/bbTdpDEOw6HOcPpDRowYwdSpU3B1dc3GBk3Az8B8unTpzPr1v5IvX75sbE9EREREJIcxpsDcZ2Heczw/aBArli3F19fX3lHds1z3zNf1A2eJPXMJDx8fkg8fwZSQkP2NJp3E8VBbnFMPULZqFV555ZVsbjAZg2ExZvMehg8fzmuvvYbBVtvpi4iIiIjkBPGXcZjcGcOhX/jqm2/o16+fvSO6b7lu5uv8z/twdHLC0dHxxvNe2e36Vhz3VSXYJ4YtmzfaYPbpGo6OU3F1PcKCBQt4/fXXlXiJiIiIyMPl3H4cx1TH+9wu1v7wwwOReEEuTL4u/LSPylWrYDAYMP3xF3lDQrKvMcvGGpH8tnMb5cuXz762gP9urNGhQ4dsbk9EREREJIf5Z2ONwn5O7Ny+jfr169s7oiyTu5Ivk4mYdQdo0rCR5VCJIkWyvp1bNtbo1vXGxhrBwcFZ344VbawhIiIiIg+5WzbW2L51M5GRkfaOKEvlqme+rp+KAaBhw4Z88umnABQrWpQN69dD4iEwZ8WuhyYMp96Cy4v4eMQIXn755Qws+4sDzt9Hm38Bv9KlSzcmT/4fbm5u91GXiIiIiMh9uhoN0Xtt2+avX8D6L3l+8AuMHj0KR0dH27ZvA7kq+QLw8PSkZs2aeHl7kxAXR8WKFW+c2N8yy9pw8/BizqJFtGnT5q5lvby8OH9+O3B/7xrQxhoiIiIikhN4eHhybdX7sOp9m7br6OT0wGyscTu5KvnavHkzwcHBuLu7s3PbNi5cuECFChWoW7cucXFxWdZOREQEIRl8luzXX9dx4sSJ+2ovb968FC1a9L7qEBERERHJCrt/30V0dHSW1vnyyy8zatSoO5YJCQkhIiIiS9vNaXJV8lWjRg3LzwUKFKBAgQIAlCtXzl4hERoaSmhoqN3aFxERERHJShEREVmeBAUEBFCzZs0srTM3ylXJl4enB2az2a4xpKSk4O7hnqV1BgbmYcb0WdStW/euZYcNG8Yno0dnafsiIiIikvsFBwYyY84catWqddeyr732Gp99OhaA9u3aWZ0rX74y+/fvy9LYUlJScHPzuGOZAgXCWbBgrl0nVrJbrkq+HFxMtHuvlL3DyHK/fX+Gxo0bMWHCV/Tt2/eOZffu3UtYUhK9bBSbiIiIiOQOy86coWH9+kycNIleve48WtyzZw+FXZJoGATj5sy1mgTYu/dPUlMfAapnaXxJSXc+f/ToXKpXj2LOnJm0bds2S9vOKXJV8uXh40KTgQ/es1ENno5kxsDf6NevH3v27mHUyDvv7lIAeHAfQxQRERGRe/G4ycRQk4nevXuzZ88ePvroozuOKQt7wiflYfl5Jw4dPPCfs7WA3tkZbhpGY1dMpkG0b9+ejz76iFdeeeWB24wud73n6wHl5OxArwmV6f55RT779FNat2nFtWvX7B2WiIiIiOQiLsAY4F1gzOjRtG/bluvXr9/xGicHGFo0lejos+zbl7VLDTPPE7P5G8zmwQwdOpTHH+9N0t2my3IZJV85hMFgoMnAorywog6/bvyZmlHVOXLkiL3DEhEREZFcxAA8BUwzm/l51SpqVa/OsWPH7njN4+EQ7G7g448+skWId+EAvApMYNasudSr15Dz5+/nfbo5i5KvHKZss3y8saUBMYlnqFqtCr/++qu9QxIRERGRXKYxsMxo5OqBA1SrVImNGzfetqyrIwwtbmbmrJkcPXrUdkHeUXtMpoXs2HGIihWrsnv3bnsHlCWUfOVA+Uv48ObWBgSXdaFx40b873//s3dIIiI25d6/P3mPHSMoLg6/NWtwqlCBYLMZt1seIHfr1YuA3bsJSkgg76lTeH7wATg4WJ0PNptxqlABvxUrCIqNJc+BA7j17JmmPZcWLQjYsoWg+HgCz5/H+8svwePOu3KJiOR0JYCVRiNFrl6lYf36TJ069bZl+xUCfxcDI0eMsF2Ad1UJo3EV5875UqNGLRYvXmzvgO6bkq8cyiuPKy+tqUOtPgXp168fLw15CaPRaO+wRESynWvr1vh8/TVJa9ZwpX17kteuxfe776zKeLzwAj6TJpG8ejVXWrcmbsQIPAYNwmv48DT1+c6cSfKaNVxp147UXbvwmTIFxxIl/m2vY0f8liwh9c8/udK+PddfeQXXDh3w1S++ROQBkAeYazLRKTWV3r1788orr6Q7pvR0ghcijUye/D+7v9rJWn6MxkUkJdWnffv2jBgxIofFlzm5arfDh42TswO9vqpMaGkfPn1hLHv2/IVB+bKIPOA8hw0j6ccfud6/PwDJa9ZgcHbG64MPADB4eeH57rvEjxxJ7Btv3Lho7VpITsb7k0+IGzUK86VLlvrix48nYcKEG3Vt2kRQy5a4dexI3D+Jmvfo0STOncu1J5+0XGOKjsZvxQoc338f4969trhtEZFsc3MjjuLAu6NHs+fPP0k1mXD9T7kBRWDkQTPXUk22D/KObmzEAaMYOnQov/++my++GE9AQIC9A8u0XDWST05JsXcINmcwGGgyqBiDl9dm3fofWfPDauLtHZSISHZxcMCpYkWSliyxOpx0y1IT56goHLy9SZw3DxwdLZ/ktWsxeHjgVKaM1bXJa9b8+yU+HuPx4ziEhQHgWKwYjhERJH33nXVdv/wCJhPOVapk372KiNjQfzfiWHPrv43/8HWG5wobMZjNQJytQ7yLmxtxfMmcOfPJkyfQ3gHdk1yVfD3Myj0SQo/xlTAZzZywdzAiItnEITAQg7Mz5gsXrI6bbtnpyiFvXgDy7NpFcGqq5ZP30CEAHAsUsL72yhWr7+bkZAxublZ1+S1aZFVXcEICBienNHWJiOR2jYEP7nD+0TC4Me+VU3fd7gC8zM0oc5tctezQxdnZ3iHYzaYZx5n21G8E5PWn8MXL9g5HRCRbmC5cwJySgiHQ+jeaDkFB/5b5Z0nhlfbtMZ48maYOYyZ26rpZ17UBA0jZujXt+TNnMlyXiEhuMBd4zcEBf19fIO2Y8ttj4GgAo7morUPLADMwEYPhY5o0ecTewdyTXJV8PYxMJjML3viT5R//Ta/evbh86TKx/1mOIyLywDCZSN21C7e2bUn4/HPLYdd27Sw/p2zejDkuDoewMJIWLbqv5ox//43x5EkcCxcm4csv76suEZGczAgMB74Enujdm3PnzsFvy63KXEiCr485/DOn5GbzGO8sGXgNmMmQIS/zUY54J1nmKfnKwRJjU5jYYxu7lpxh1KhRvPTSS3Ts2NHeYYmIZKu4Dz7Ab8kSvCdOJGnePJwqVvx3i3mTCfPVq8S+9RbeI0fiGBZG8rp1YDTiWLgwrm3bcqVjR0hIyHB71198Ed9ZszB4epK0fDnmuDgcw8NxbdmS2Ndfx3jwYPbcqIiIjcQCzxoMrAXGfvIJzz//PG3btk1T7tMD4ODkgoMplZy1yXYMDg79cHDYyTfffEvv3r3tHdA9U/KVQ8WciOPzNpu5eDiRJUuW0KpVK3uHJCJiE0lLl3Lt6afxfP113Hv0IGXrVq4/8wz+P/yA+epVAOI/+QTj6dN4vvgiHgMHYk5JwXj4MEnLlkFycubamz+fK1eu4PnGG7j16AGA6dgxklatwnTuXJbfn4iILZ0AHnd05LSrK8vmz6d58+bplruSDOOPOvLMcwP49LPP0y1jH/txdHwcH594li79iVq1atk7oPtiMOeijfKDwn0YdTz9vzAPkkObLzKu3RZ83QNYtnQFZcuWtZzr0KEDVxcuZJYd4xMRsTW3J57A93//40JEBKbjx+0djohIrrAV6OvoiHf+/CxduZLSpUtbzrVp0wZ2LGVJ7Rvfh++D9w84c/TYcQoWDCc19X2gtz3CvsWPODo+TdGi4axcuZSIiAg7x3P/NPOVw2yacZxv++6gatVqLPx+EUG3PGQuIvIwMPj74/X22yT/9BOm69dxrloVzzfeIHHRIiVeIiIZNBd42cGBGjVqsGDRIvL+s7treuJS4dPDjvTt24+QkBDbBXlbNzfWeJdHHmnJ7Nkz8fb2tndQWULJVw5hMpn5ftifLPvoxsYaX3/1Na6u/331nYjIQyAlBcfISHy6dcPg54fpwgUSp08n9tVX7R2ZiEiO99+NNSZMmICLi8sdr/nmCFxONvPyK6/YIsS7SLuxhqOjo72DyjJKvnKAxNgUvum5nd8Wn2bkyJEMGTIEg8Fg77BEROzCHBvLldat7R2GiEiuEws8YzDwI/9urHG3MWWSEUYdcqJH9245YFnfzY01duT6jTVuJ1clX7GXk5j/xp/2DiPL/bHsHDFHElm8eDGtMzDgOAzkzs01RURERCS7rM7Axhq32ncdem+H6Hgjr73++n9rA6KzJc7bcXJahLd3HEuW/ETt2rVt2rat5KrkKyhvPv6Ydd3eYWS5wMAIlm2aYrWxxu1ERUWxa+tWFtsgLhERERHJPYJDQpg/darVxhq3U6tWLf7ctZMtZhj0XEeKFy9uOdegQSP27PkbyPhL67NCoUIRzJgxJQfMwGWfXLXboYiIiIiISG7lYO8AMio1NZV169aRmppq71AeaOpn21A/Zz/1sW2on7Of+tg21M/ZT31sG+rnnC3XJF9Go5FffvkFY8563fYDR/1sG+rn7Kc+tg31c/ZTH9uG+jn7qY9tQ/2cs+Wa5EtERERERCQ3U/IlIiIiIiJiA0q+REREREREbCDXJF+Ojo7Uq1fvgXrDdU6kfrYN9XP2Ux/bhvo5+6mPbUP9nP3Ux7ahfs7ZtNW8iIiIiIiIDeSamS8REREREZHcTMmXiIiIiIiIDSj5EhERERERsQElXyIiIiIiIjbgZO8AMuL48eNs2rSJM2fOEBsbS5cuXShRooS9w8o11q9fz99//83FixdxcnKiQIECNG7cmLx581rKpKamsnr1avbs2UNqaipFihShRYsWeHl5WcqsXLmSkydPcv78efLmzcvTTz9tj9vJFTZs2MCPP/5I9erVeeSRRwD1cVa5du0aa9eu5dChQ6SkpBAQEEDbtm3Jnz8/AGazmXXr1vHbb7+RmJhIgQIFaNmyJXny5LHU8euvv3Lw4EHOnj2Lo6MjQ4cOtdft5Dgmk4l169bx559/Ehsbi7e3N+XLl6du3boYDAZAfXwv7vb/WEb6NCEhgZUrV7J//34MBgMlS5akefPmuLi4ADf+jVm2bBnR0dFcuHCBYsWK8dhjj9n8Xu3pTv1sNBr56aefOHToEJcvX8bV1ZXChQvTuHFjvL29LXWon+8sM2OyZcuWsXPnTpo1a0aNGjUsx9XHd5eRfr5w4QJr167l+PHjmEwmAgMDefTRR/H19QU07sipcsXMV3JyMsHBwbRo0cLeoeRKx48fp2rVqvTt25eePXtiMpmYMWMGycnJljKrVq3iwIEDdO7cmd69e3P9+nW+++67NHVVqFCB0qVL2zL8XOf06dPs3LmT4OBgq+Pq4/uXkJDA5MmTcXR0pHv37jz77LM0bdoUNzc3S5mNGzeydetWWrZsSb9+/XBxcWHGjBmkpqZayhiNRkqVKkWVKlXscRs52saNG9mxYwfNmzdnwIABNG7cmE2bNrFt2zarMurjzLnb/2MZ6dPvv/+e8+fP07NnT7p168aJEydYunSp5bzJZMLJyYlq1apRuHDhbL+nnOhO/ZySksLZs2epW7cu/fv3p0uXLsTExDB79myrcurnO8vomGzfvn2cOnXKKrG9SX18d3fr50uXLvHtt9+SN29eevXqxdNPP03dunVxcvp3XkXjjpwpVyRfRYsWpWHDhpQsWdLeoeRKPXr0oEKFCgQFBZEvXz7atm3L1atXiY6OBiAxMZFdu3bRrFkzChUqRP78+Wnbti0nT57k1KlTlnqaN29OtWrV8Pf3t9et5HjJycl8//33tG7d2iohUB9njY0bN+Lr60vbtm0JDQ3F39+fyMhIAgICgBuzB1u3bqVu3bqUKFGC4OBg2rVrx/Xr1/n7778t9TRo0ICaNWumSZAFTp48SfHixSlWrBh+fn6UKlWKyMhITp8+DaiP79Wd/h/LSJ9euHCBQ4cO0aZNG8LCwihYsCDNmzfnr7/+4vr16wC4uLjQqlUrKleubPWb7YfJnfrZzc2Nnj17Urp0afLmzUtYWBjNmzcnOjqaq1evAurnjMjImOzatWusXLmSDh064OBgPdRUH2fM3fr5p59+omjRojRp0oSQkBACAgIoXrw4np6egMYdOVmuSL4kayUlJQHg7u4OQHR0NCaTyeq3S3nz5sXX15eTJ0/aJcbcasWKFRQtWjTNb+rUx1lj//79hISEMG/ePEaNGsXXX3/Nzp07LeevXLlCbGysVT+7ubkRFhamfs6gAgUKcPToUWJiYgA4e/YsJ06coEiRIoD6ODtkpE9PnTqFm5ubZXktQOHChTEYDFYDKcmcm/8f3vxlmfr5/pnNZhYuXEhUVBRBQUFpzquP75/ZbObgwYMEBAQwY8YMRo0axaRJk6x+AaZxR86VK575kqxjNptZtWoVBQoUsPyjGBsbi6Ojo9VMDYCnpyexsbH2CDNX+uuvv4iOjubJJ59Mc059nDUuX77Mjh07qFmzJrVr1+bMmTOsWrUKR0dHKlSoYOnLm7/5u8nT05O4uDh7hJzr1K5dm6SkJMaPH4+DgwMmk4mGDRtSrlw5APVxNshIn8bGxqY57+DggLu7u/4NuUepqamsXbuWsmXL4urqCqifs8KGDRtwcHCgevXq6Z5XH9+/uLg4kpOT2bhxIw0aNKBx48YcOnSIuXPn0qtXLyIiIjTuyMGUfD1kli9fzvnz53niiSfsHcoD5erVq6xatYqePXtarbeWrGU2m8mfPz+NGjUCICQkhPPnz7Nz504qVKhg3+AeEHv27OHPP/+kY8eOBAYGcvbsWVavXo23t7f6WB4YRqORefPmYTabadmypb3DeWCcOXOGrVu38tRTT1k26JGsZzabAShevDg1a9YEIF++fJw8eZKdO3cSERFhx+jkbjRKfIisWLGCgwcP0rt3b3x8fCzHvby8MBqNJCYmWv2GJC4u7qFda51Z0dHRxMXF8fXXX1uOmc1mjh8/zrZt2+jRo4f6OAt4e3sTGBhodSxv3rzs27cPwNKXcXFxVg95x8XF6dmjDPrhhx+oVasWZcqUASA4OJirV6+yYcMGKlSooD7OBhnpUy8vrzQziyaTiYSEBP0bkklGo5H58+dz9epVHn/8ccusF6if79eJEyeIi4tj7NixlmNms5k1a9awZcsWBg8erD7OAh4eHjg4OKT7/+HNJYUa2+VcSr4eAmazmZUrV/L333/Tq1evNA9VhoSE4ODgwJEjRyhVqhQAFy9e5OrVqxQoUMAeIec6hQoV4plnnrE6tnjxYvLmzUutWrXw8fFRH2eBAgUKWJ5FuikmJsayra6fnx9eXl4cOXKEfPnyATee6Th16pR23cuglJSUNL+xNhgMlt+0qo+zXkb6NCwsjMTERM6cOWN5Vubo0aOYzWbCwsLsFntuczPxiomJoVevXnh4eFidVz/fn3LlyqV55nnGjBmUK1fOMnOuPr5/jo6O5M+fP83/h5cuXbL8f6ixXc6VK5Kv5ORkLl26ZPl++fJlzp49i7u7u+UvmdzeihUr+PPPP3nsscdwdXW1rPV1dXXF2dkZNzc3KlasyJo1a3B3d8fV1ZWVK1cSFhZm9Q/hpUuXSE5OJjY2ltTUVM6ePQtAYGAgjo6Odrm3nMLV1TXNg8XOzs64u7tbjquP71+NGjWYPHky69evp3Tp0pw+fZrffvuNVq1aATeShOrVq7N+/Xry5MmDn58fP//8M97e3lbvR7l69SoJCQlcvXoVs9ls6eeAgADLe2YeVsWKFWP9+vX4+voSFBREdHQ0W7ZssQyc1Mf35m7/j92tTwMDAylSpAhLly6lVatWGI1GVqxYQZkyZaxmyy5cuIDRaCQhIYHk5GRLv99M6h50d+pnLy8v5s2bR3R0NF27dsVsNlv+P3R3d8fR0VH9nAF3+7v834TWwcEBLy8vy7tF1ccZc7d+joqKYv78+RQsWJBChQpx6NAh9u/fT+/evQE0tsvBDOabv87MwY4dO8bUqVPTHC9fvjzt2rWzfUC5zLvvvpvu8bZt21oGVDdfxPfXX39hNBqJjIykZcuWVlPTU6ZM4fjx42nqef755/Hz88uO0HO1KVOmkC9fvjQvWVYf358DBw7w448/EhMTg7+/PzVq1KBy5cqW8zdfVrtz504SExMpWLBgmpfVLlq0iN27d6ep++aDyg+zpKQkfv75Z/7++2/LMrgyZcpQr149y3/E6uPMu9v/Yxnp04SEBFasWMGBAwfSfTEtwKeffmrZNv1Wb7/9dvbcWA5zp36uX78+n332WbrX3fr3Uv18Z5kdk3366afUqFEjzUuW1cd3lpF+3rVrFxs2bODatWvkyZOH+vXrW/0STOOOnClXJF8iIiIiIiK5nd7zJSIiIiIiYgNKvkRERERERGxAyZeIiIiIiIgNKPkSERERERGxASVfIiIiIiIiNqDkS0RERERExAaUfImIiIiIiNiAki8REREREREbUPIlIiIiIiJiA0q+REREREREbEDJl4iIiIiIiA0o+RIREREREbEBJV8iIiIiIiI2oORLRERERETEBpR8iYiIiIiI2ICSLxERERERERtQ8iUiIiIiImIDSr5ERERERERsQMmXiIiIiIiIDSj5EhERERERsQElXyIiIiIiIjag5EtERERERMQGlHyJiIiIiIjYgJIvERERERERG1DyJSIiIiIiYgNKvkRERERERGxAyZeIiIiIiIgNKPkSERERERGxASVfIiIiIiIiNqDkS0RERERExAaUfImIiIiIiNiAki8REREREREbUPIlIiIiIiJiA0q+REQk15syZQoGg4Fjx45l+tr69etTpkyZrA9KRETkP5R8iYjkYDeTih07dmSofLVq1TAYDEyYMOG2Zf788086depEeHg4bm5uhIaG0qRJE8aNGwfAb7/9hsFgYNiwYbet4+DBgxgMBl588cUMxXXs2DEMBgOjR4/OUPnczGw2M336dOrWrYufnx8eHh6ULVuW9957j7i4uDTl69evj8FgsHwCAgKoWrUqkydPxmQysW7dOqvzd/qIiEjO5mTvAEREJGscPHiQ7du3ExERwcyZM3nmmWfSlNm0aRMNGjSgYMGCPPnkk+TLl4+TJ0+yZcsWPvvsMwYOHEilSpUoUaIEs2fP5oMPPki3rVmzZgHQo0ePbL2njOrZsyePPfYYrq6udo3DaDTSrVs3vvvuO+rUqcM777yDh4cH69ev591332XevHmsXbuW4OBgq+vCwsL46KOPALhw4QLTpk2jb9++HDhwgBdeeIHp06dblX/ttdfw8vLijTfesNm9iYjI/VPyJSLygJgxYwZBQUGMGTOGTp06cezYMSIiIqzKDB8+HF9fX7Zv346fn5/VufPnz1t+7t69O2+++SZbtmyhRo0aadqaPXs2JUqUoFKlStlxK5nm6OiIo6OjvcNg5MiRfPfddwwZMoRRo0ZZjvfv359HH32Udu3a0bt3b1auXGl1na+vr1Ui+9RTT1G8eHHGjx/P+++/nybJ/fjjj8mbN2+OSX5FRCRjtOxQROQBMWvWLDp16kSrVq3w9fW1zE7d6vDhw5QuXTpN4gUQFBRk+bl79+6WOv9r586d7N+/31Lmfhw5coTOnTsTEBCAh4cHNWrUYPny5WnKjRs3jtKlS+Ph4YG/vz9VqlSxii29Z74WL15My5YtyZ8/P66urkRGRvL+++9jNBrTjWXv3r00aNAADw8PQkNDGTlyZKbuJSEhgVGjRlGsWDHLLNatWrduTa9evVi1ahVbtmy5Y103+yIuLo4LFy5kKg4REcm5lHyJiDwAtm7dyqFDh+jatSsuLi506NCBmTNnpikXHh7Ozp07+euvv+5YX6FChYiKiuK7775Lk6zcTHq6det2XzGfO3eOqKgoVq9ezbPPPsvw4cNJTEykTZs2LFy40FLum2++YdCgQZQqVYpPP/2Ud999lwoVKrB169Y71j9lyhS8vLx48cUX+eyzz6hcuTJvvfUWQ4cOTVP28uXLPPLII5QvX54xY8ZQokQJXn311TQzVHeyYcMGLl++TLdu3XBySn9hyeOPPw7AsmXL7lrfkSNHcHR0TDdRFhGR3EnLDkVEHgAzZsygQIEC1KpVC4DHHnuMyZMn8/vvv1OhQgVLuSFDhtC8eXMqVKhAtWrVqFOnDo0aNaJBgwY4Oztb1dm9e3cGDBjAjz/+SNOmTQEwmUzMnTuXmjVrUrhw4fuK+eOPP+bcuXOsX7+e2rVrA/Dkk09Srlw5XnzxRdq2bYuDgwPLly+ndOnSzJs3L1P1z5o1C3d3d8v3p59+mqeffpovv/ySDz74wOr5sDNnzjBt2jR69uwJQN++fQkPD+d///sfzZs3z1B7e/fuBaB8+fK3LXPz3L59+6yOG41GLl68CMDFixeZMGECv/32G61bt8bDwyND7YuISM6nmS8RkVwuNTWVuXPn0qVLF8uOdw0bNiQoKCjN7FeTJk3YvHkzbdq0Yffu3YwcOZJmzZoRGhrKkiVLrMp26dIFZ2dnq+V9v/zyC6dPn86SJYcrVqygWrVqlsQLwMvLi/79+3Ps2DFLMuPn58epU6fYvn17puq/NfG6fv06Fy9epE6dOsTHx/P3339blfXy8rJ6fsrFxYVq1apx5MiRDLd3/fp1ALy9vW9b5ua5a9euWR3/+++/CQwMJDAwkJIlSzJu3DhatmzJ5MmTM9y+iIjkfEq+RERyuTVr1nDhwgWqVavGoUOHOHToEEePHqVBgwbMnj0bk8lkVb5q1ap8//33XL58mW3btvHaa69x/fp1OnXqZEl4APLkyUOzZs1YuHAhiYmJwI3ZJCcnJx599NH7jvv48eMUL148zfGSJUtazgO8+uqreHl5Ua1aNYoWLcqAAQPYuHHjXevfs2cP7du3x9fXFx8fHwIDAy0J1tWrV63KhoWFpdmq3d/fn8uXL2f4fm4mVjeTsPTcLkGLiIjghx9+YO3atWzYsIGzZ8+ybNky8ubNm+H2RUQk51PyJSKSy92c3Xr00UcpWrSo5TN37lxOnz7NL7/8ku51Li4uVK1alQ8//JAJEyaQkpKSZmlfjx49uHbtGsuWLSM5OZkFCxbQtGlTAgMDs/2+bipZsiT79+9nzpw51K5dmwULFlC7dm3efvvt215z5coV6tWrx+7du3nvvfdYunQpP/zwAyNGjABIk5DebqdEs9mcqTgB/vjjj9uWuXmuVKlSVsc9PT1p3LgxjRo1olatWlabn4iIyINDz3yJiORicXFxLF68mC5dutCpU6c05wcNGsTMmTNp0KDBHeupUqUKANHR0VbH27Rpg7e3N7NmzcLZ2ZnLly9nyZJDuLH5x/79+9Mcv7kkMDw83HLM09OTLl260KVLF5KTk+nQoQPDhw/ntddew83NLU0d69atIyYmhu+//566detajh89ejRLYk9P7dq18fPzY9asWbzxxhvpJnTTpk0DoFWrVtkWh4iI5Fya+RIRycUWLlxIXFwcAwYMoFOnTmk+rVq1YsGCBSQlJQHw888/pzubs2LFCoA0ywDd3d1p3749K1asYMKECXh6etK2bdssib1FixZs27aNzZs3W47FxcUxceJEIiIiLLNDMTExVte5uLhQqlQpzGYzKSkp6dZ9M/G59V6Tk5P58ssvsyT29Hh4eDBkyBD279+f7suPly9fzpQpU2jWrFm6704TEZEHn2a+RERygcmTJ7Nq1ao0x3/88Ufy5MlDVFRUute1adOGb775huXLl9OhQwcGDhxIfHw87du3p0SJEiQnJ7Np0ybmzp1LREQEffr0SVNHjx49mDZtGqtXr6Z79+54enpmyT0NHTqU2bNn07x5cwYNGkRAQABTp07l6NGjLFiwAAeHG78fbNq0Kfny5aNWrVoEBwezb98+xo8fT8uWLW+7uUVUVBT+/v706tWLQYMGYTAYmD59eqaWEd7rPe3atYsRI0awefNmOnbsiLu7Oxs2bGDGjBmULFmSqVOnZmsMIiKScyn5EhHJBSZMmHDbcz179rztM0uNGjXCw8ODGTNm0KFDB0aPHs28efNYsWIFEydOJDk5mYIFC/Lss88ybNiwdN8p1bBhQ0JCQoiOjr7nJYc3k55b4wwODmbTpk28+uqrjBs3jsTERMqVK8fSpUtp2bKlpdxTTz3FzJkz+eSTT4iNjSUsLIxBgwYxbNiw27aXJ08eli1bxksvvcSwYcPw9/enR48eNGrUiGbNmt3TPWSEo6Mj3333HdOmTWPSpEm8+eabJCcnExkZydtvv81LL72UZcmriIjkPgZzdv8aUEREHnp//PEH5cuXZ9KkSfTt29fe4YiIiNiFnvkSEZFsd/MdXf/d5U9ERORhopkvERG5Z8nJyVy6dOm2548ePcpPP/3EqFGjyJcvH3v37rU8y5WbXLhwAaPReNvzLi4uBAQE2DAiERHJjfTMl4iI3LNNmzbddRt7Dw8P6tSpw7hx43Jl4gU3Xkx986XP6alXrx7r1q2zXUAiIpIraeZLRETu2eXLl9m5c+cdy5QuXZqQkBAbRZQ9Nm7cSEJCwm3P+/v7U7lyZRtGJCIiuZGSLxERERERERvInes/REREREREchklXyIiIiIiIjag5EtERERERMQGlHyJiIiIiIjYgJIvERERERERG1DyJSIiIiIiYgNKvkRERERERGxAyZeIiIiIiIgNKPkSERERERGxASVfIiIiIiIiNqDkS0RERERExAacbNHIlStXOHv2LCkpKbZo7r45OjoSEBBAcHAwBoPB3uGIiIiIiMgDIFuTr+joaJYtW8aZM2eys5ls4+fvT+NGjShdurS9QxERERERkVwu25Kv8+fPM336dDx8fGjUoT35IyJwdnUlN8wjpaamcvFMNH9t386CBQsAlICJiIiIiMh9ybbka9OmTTi5utK61+O4urllVzPZwsnZmbDIwuQvFMHqud/x87p1lCpVSksQRURERETknmXLhhsmk4n9+/dTtGyZXJd43crBwYHSVSoTc/EiFy9etHc4IiIiIiKSi2VL8pWQkEBiYiJ5Q0Kyo3qbunkPly5dsnMkIiIiIiKSm2VL8mU0GgFwdEx/VePzffpQv2y5+27nr99/J8TBkU3r1t13Xbfj6HTjHm7ek4iIiIiIyL2wyVbz//XCsGHEx8XZo2kRERERERG7sEvyFREZaY9mRURERERE7CZblh3eza3LDudOmUKIgyN/7tpFtxYtKOzlTVSx4nw3bVqa68Z+MJxyIfmJ9PbhiY4duXj+fJoyZrOZCaPHUKt4CcLd3KkeWYSvx35qOX/y+HGK+fnz7pCXra7r1qIFNYsW04yciIiIiIhkC7skX+kZ0KMn9Zo05duF31OmYgUG93mCA/v2Wc5PHv8FI996i049ejBp/jzCCxXmpX5PpqnnzecHM/Ltt+n8+ONMX7aULr16MXzoUKZ+9RUABcLDeW/sJ0z89FM2/fILAFMnTODXH9by+dQpeHh62uaGRURERETkoWKXZYfpeWLAAHo/+wwAVaOiWLt8BcsXfE+xYW9gNBoZ9/HHdOrZg7dGjQSgQbNmXLxwnvnTZ1jqOHb4MJO/+IIRE76kZ//+ANRt3JiE+Hg+ee99evbvj4ODA4/16cOqxYt5vncf/rdgPu+98irPvvwyVaOibH/jIiIiIiLyUMgxM1/1mjax/Ozh6UlYeDjRp04BcObUKc6eOUPzdu2srmnVsaPV91/XrgWgZceOpKamWj51Gjfi/NmznD550lJ29MSJJMTH06ZWbQoVKcKQd97OpjsTERERERHJQTNfPn5+Vt9dXFxISkoE4Hx0NAB5g4KsygQGB1t9v3QxBrPZTOlA63I3nTl5kgLh4Za66jRqyKI5c+nx5JO4uLhkxW2IiIiIiIikK8ckX3cS9M+Ljv+7wcaFc+esvvsH+GMwGFi8/lec00mmihQvbvn5p1WrWDRnLmUqVmT0u+/SqlPHNMmdiIiIiIhIVskxyw7vJH9YGMEhIaxctMjq+LIFC6y+127UCIDLMTFUqFIlzcfL2/vG+UuXeKnfk7Tv+hjf//wTbu7uDHnqKZvci4iIiIiIPJxyxcyXo6Mjz736Km8OHkxgUDB1mzTmlzU/sOnndVblIosVo8+zzzLw8V48M2QIlapXIyUlhSMHDrJx3c9MWbgQgNcGDADgw/Hj8fbx4bNvJ9O5cRPmTplCl969bXx3IiIiIiLyMMgVyRdA34HPce3KFb798kumTJhA3caNGP3NRLo1b2FV7oPPPyOyeHGmT5zI2Pffx9PLi8jixWnVqRMAi+bMYfHc75i5Yjl+/v4A1GrQgL4DB/Lm4Beo1bAhYQUL2vz+RERERETkwWYwm83mrK702rVrjB07luZdu1KwaJGsrt6mkhITmTJyFJ07d6ZUqVL2DkdERERERHKpXPHMl4iIiIiISG6XzclXlk+q2V7WTwyKiIiIiMhDKFuSL2dnZwBSU1Kyo3qbunkPN+9JRERERETkXmRL8uXm5oaXlxdnT57KjuptKvrkSQACAwPtHImIiIiIiORmmd7t0GQycfToUfbv38+1a9cwmozplnN0dOTAH39Qtno1vP387jdOu0hJSeGPzVtwc3dj1apV+Pr6UqJECcLDw3Fw0ONyIiIiIpK7mc1mzp49y969e4mJiSElJdXeIWWKwQCurq6Eh4dTsmRJPD097R3SHWVqt8OEhASmz5hB9JkzePv54hcYgINj+klIakoq509G4+TsTKkqlckfEYGLi8uNHsrhjKmpXDgTzb7ffuPKxYsEhuXDxcWFyxdiiL16nYIFC9K9e/cb9yMiIiIikguZTCYWL17CH3/sxs3Nnbx58+PomNsetTGTmBjP+fOnMBgMdOnShWLFitk7qNvKVPI1ddpUos+epUmX1uQrGIrhLolU7LXrbF+7geMHjpCSlHzfwdqSwcGBsMLhVKxbjeAC+YEbvxk4feQEP85bRkREBF0f62rnKEVERERE7s2PP/7Ixo0bqVOnNUWLls/VK7sSE+NZv34ZJ07sp3///gQHB9s7pHRleNnh5cuXOXb0GA06NCckPCxD13j5eNOgQ3NSU1O5fvlqrtmAw9HRCS9fb1zcXK2OGwwGwiLDqdakDhuX/0RsbCxeXl52ilJERERE5N6YTCZ27fqdUqWqUrx4RXuHc9/c3Dxo2LAjs2d/wh9//EGTJk3sHVK6Mpx8HT58GAcHB8KLF858I05O+AfmyfR1OVWhkkXZsOxHjhw5Qrly5ewdjoiIiIhIppw7d464uFgKFSpl71CyjKOjIwULFufQoUM5NvnK8NxifHw8ru5uOOs5J9w83HFyciI+Pt7eoYiIiIiIZFpCQgIAnp4+do4ka3l6+uToMXqGky+z2XzbZ7wG936GhmVqZFlQd1PSryBj3vkoy+ob3PsZQg2+aT4/r1p722sMDg5k4nE5EREREZEc4+Y4Nr3nvObMmUxoqIEjRw5aHZ88eRyhoQZGj37b6vjly5cIC3Pgyy9HMnfuFEJDDZZPyZJ+tG5dk9WrF1vKnzx5zKpM4cJu1K1bgtGj37YkhTeNGfOOVdlbP+PHf5wmdoccPkbP9FbzD6rwwhGMmznJ6ljRkjl3pxQRERERkexQpUotAHbs2EThwkUtx7dv34i7uwc7dmyyKr9jxybMZjNVq9bmyJEDAMycuQpvb19iYi4wceInPPFEO2bOXEX9+s0s1w0d+iFRUQ2Ij4/jhx+WMHbse1y4cI4RI76yqt/NzZ3vvvspTZyhoQWz7J5t5aFOvhISEnB3dwfAzd2dyjWq2jkiERERERH7KlKkOHnyBLJjx0YefbSX5fj27Rt59NHezJ8/DaPRiKOjIwA7dmzEzc2N8uWrWJKvcuUqExCQF4CoqPpUrVqAyZPHWSVfhQoVpXLlG6vn6tRpxMGD+5g/fxofffSl1Yycg4ODpVxuly37SZ45dZqBPZ6kTN5CRLoH06Fuc/7YucuqzLxps2lXuxmlA8Ip5V+QTvVbsmvbzjR1rV68nLolqlDYLYiW1Rrw+/a0ZQDWLl9Nq+oNiXQPpmxgYYY+8wLxcXGW85vWrSfU4Mva5at5slNPivuE8VTnXunWJSIiIiLyMKtatRbbt2+0fD99+gTR0afo2/d5kpIS2bfvD8u57ds3Uq5cldu+A9fLy5vIyOKcPHn0jm2WKVORxMQEYmL+3959x8d8/wEcf91d9pa9EEnE3mLvPWLvVbOqRlGjWlrVoqjq0qHVllJKlNjErr1H+dk7EkQWicy7+/2RujoJEpK7C+/n43GPR+47Pp/3faTN+32f7/fzjc6bD2GC8rz4io+Lo0OdFpw58Q9Tv/2cn/5ahI2tDV0bteXe3f8GMuLaDTq/0Z0fQxcyd8kveBfxpVO9lly+cEl3zOkTp3izUx+KFQ/g55WL6dK3J0O69iMtNVWvz3Urwujftjsly5Vm/qo/mDTrEzauXMuYgcOzxPfe4JH4BRRj/qrFDBk7Qrf92qUrlHQsjJ+FKy2q1GNT2Lq8HhohhBBCCCEKhODg2ly8eJb4+Dggs8Dy9i5MQEAQpUtX0BVm6enpnDx5mGrV6jy1LbVaTWTkTTw8vJ/ZZ0TEdezs7HUzZo/LyMjI8iqI8vyyw/lf/cD9+ATWH9qOq7sbAHUa16duUBV+nP0Nk2Z9CsDoj97TnaPRaKjXtCEnDh1l+YI/eH965k183834Ep8ivvwatkQ3rWllbaVXVGm1Wj4d+yFtu3Vk9vy5uu3uXp70adWZUR+Op0SZUrrtTdu2ZOLMT/RiLlupPBWDKxNUpiT34xP4/YdfGNihF/NCFxLSuX3eDpAQQgghhBAmLji4DlqtlqNH99O4cSuOHNlH1aq1AKhatRZHjuyjf//h/PPPMVJSUggO1i++1Go1GRkZxMRE8/XXU7lzJ4p33/1Y7xiNRkNGRgYPHyYRHr6GDRv+4r33puny/kcePkyiaFHzLDGuWrX7mUWfKcrz4mtX+HZqNayLk3MhXUWqUqmoUb82Jw7/d+nhxbPnmfHBJxzZd1BvRuzKYzNfxw8eoVnbVnr/AK07t9Mrvi5fuETE9RtM+eozvQq4Rv3aKJVKTh45rld8NW7933Wmjwwa+bbe+2ZtW9G2VlNmfzRdii8hhBBCCPHaKV++ClZW1hw+vJfGjVvp7vcCqFKlJtOnTwAyZ8QUCoWuMHukYkVP3c9WVtaMHDmJXr3e1Dvm7be76b1v1647Q4eOzxKLlZU1K1f+nWV7YGDJF/psxpTnxVfsvRiOHThMUfOsD1X2CygGQOKDB/Ro1gEXNxcmz5mOb9HCWFpZMnbQO6Sm/HdJ4d2oO7i460872js4YGVlpXsfdy8GgIEdemUbT+TNW3rv3TzcnvsZlEolrTu1Zer4j/QW5RBCCCGEEOJ1YG5uTsWKwRw5spekpETOnj2lN/N169YNIiMjOHJkL0FBpXFyKqR3/p9/bsXBwRFHx0L4+hbFzCxr2TFx4kxq125EQkI8CxbMZfXqP6lZswF9+ryld5xSqaRChar592ENKM+Lr0LOhSjWognjPp2YZZ+FpSUAR/cfIiriFgvXLaNMhXK6/Q8SEvDy/e9aUHcvD2Lu3tNr48H9+6SkpOjeOzln/kNPmzubStWrZOnTw9tL7/3TnlUmhBBCCCGE+E9wcB3mz/+KQ4f2YGFhSZkyFQHw9S2Kh4cXhw/v5ciRfTRr1i7LuWXKVMj23q3HFSniryuqatduSKtWwcyaNYlOnXpjY2Ob55/HFOT5ght1mjTgwv/OU7xUCSpUraz3KlWuDAApyZnF0+Mrohzed5Cb127otVWxWhW2rN2IWq3WbVu/YrXeMYElg/Dy9eH6lWtZ+qtQtTKeTxRfOaHRaFgXGkaJMqVk1ksIIYQQQryWqlWrQ3LyQ377bS4VKwbrzV5VqVKLv/5axN27twkOrv3SfalUKiZOnEVs7D0WL/7ppdszVXk+8zX43eGs+iOUTvVbMXDk2/gU8SUm+h7HDx7Bw9uLwaOHUblGMLZ2dnwwbAzDJ4zm9q0oZk/+DE8f/RVQhk8YTavghgxo35O+Qwdx48o1fpz9jd5lhwqFgslzpjG85yAeJiXRuHVzbGxtuHX9JlvXb2bC9MkEBAU+Nd6I6zcY1fdt2vXohF+gPwlx8fz+wy+cPHKcn/9alNfDI4QQQgghRIFQpUpNlEol27dvYNiwCVn2TZ06DiDPFr2oV68J1arV4eefv6R//+GYm2cusqHRaDh69ECW411d3Sla1D9P+jaUPC++nF2cWXtgK7Mmfcr09yYTFxOLi7sblWtUpUWHNgC4ebgzL3Qhn46dxIB2PSkWFMjMeV/x/cyv9NoqW6kC80IXMn3Cxwzq0IsSZUvx/Z+/0at5R73j2nTpgKOTI99M+4KVi5cDUNivCA1aNH7uPV629nbYOzrw9dTZxNyNxtzCgvJVK7F44woaNG+SdwMjhBBCCCFEAeLo6ESJEmU4e/afLAtqVK1aC61Wi6enN0WKFMuzPkePnkyPHk1ZufIPunXrB0BKSjJt29bMcmyPHgOZPXt+nvVtCAqtVqvNyYG7du3i0JHD9BozOL9jKhAWfPYdDRs0oGbNrL8IQgghhBBCmLLLly+zePFievQYhZ2do7HDyTPHj//N2bOHGDt2rLFDyVae3/MlhBBCCCGEECKrHBdfCoWCHE6SvRa0Wq2snCiEEEIIIQq0Vy2/z/w8ppuj57j4srS0JDUlFfVjDzJ+XaWnpZGRnq638IcQQgghhBAFxaM8Njk5yciR5K3k5CSTztFzvOBGsWLF0KjVRFy5QdGggrWqSF67ceEqAJ6ennrPHDMEpVKJubm5zLoJIYQQQrzitFotaWlp+TI75eTkhIWFBVevnsHJySXP238apVKFSmWWL7msVqvlxo0LlCjx9JXOjS3HxZebmxvuHu4c2LQTZ3dX7J0c8jMuk5UQE8eBzbtQKBTMmzfPKDFYWVkRVCKIasHV8PHxMUoMQgghhBAif1y4cIHjx45x6fJlMvL5qrNTp/Zz6tT+fO3jSdbWNvj5laJMmeoUKvTslclzSqPRcOjQFhITEyhXrlyetJkfcrzaIUB8fDwLFizg/v37+PgXwcnVGaXq9VizQ52hJvbuPW5fv4WNnTU16lfGytrS4HFo1Bpio+M5f+YKifeT6NGjB/7+r/dMpBBCCCHEq+Lo0aOsW7cOL3dXSgcWw8nBHuUrdMVThlrN7egY/jl/mZS0NIoVK4O1te0Lt6fVaklJeUhExCWSk5No2bIl1apVy8OI81auii+Ahw8fcubMGc6dP8f9+/fRqDX5FRuJSYloNBps7KxRKoxb5KlUSuwd7Qgs6Ufx0sWMUng9LiMjg9VLw4m5E8e7745BqXw9imAhhBBCiFfVgwcPmDNnDlXKlaJ1g9qv9G0mySmpLPhrLbEJD3BweJkr6hRYWVlSpEgRypYta/JXheX6Ics2NjYEBwcTHBycH/HoPHz4kNmzZ9OwZS0qBJfO174KIjMzM2o3qsrS+au5fv06xYrl3cPthBBCCCGE4Z09exalUknjmsGvdOEFYG1lSf1qlQnduI1evXrh7Oxs7JAMwmSnS+7cuYNWq6WIv7exQzFZHt5uWFiYExkZaexQhBBCCCHES4qMjMTL3RVrK+NeYWUoxQpnzlK9TrmsyRZfaWlpAFhYWmTZty5sEz9/v8DAEWXv+rWbOCp8CFuxLlfn7d65D0eFD8eOnNRta92gM44KnyyvC+cuZduGQqHAwtJCN1ZCCCGEEKLgSk9Px8rC3NhhGIzVv3n+65TL5vqyQ0PLbsp1fdgmjh85xZtD+xk+oCd4ermzdf8aAnK5/H6FyuXYun8NJUoV19teo3YwU2d/qLetiJ/v0xt6tWekhRBCCCFeM5LcvcpMvvjKb8nJyVhbW7/w+ZaWlgTXqJLr8xwc7LM9z9HJ4YXaE0IIIYQQwpS9bN79KjDZyw6f5u1+o1iyMJSzZ87rLst7u98oAA7tP0JIoy542QZS2LEkA3sOI/ruPd25jy4R/GPBMka8OQ4/lzI0qhYCgKPChy9nfscnE2cQ4F6eIk6l+HD8VLRaLTu37aZOxaZ42xWnTeOuRNy8laXNxy87LOdXnbHDJ/LzdwsoW7QahR1L0rP9AO5Fx+iOye6yQyGEEEIIIbJz5ux5WnXtg0tAGWx8AihRrS6zvvkegH7DRlG2ViM2btlO2VqNsPLyp0rDFhw4fFSvjd//DKVOy/Y4+5ehULHSNGjTmUNHj2fp6+z5i3R8YxDO/pl9VajbhKV/hen2a7VaZn/7I0HBdbD0LIZ/pZp8+f1Pem18POML7AoX59DR49Rs1gYrL3++m78w7wemgClwM1/jPhzFvegYLpy7zPw/vgXAxc2FQ/uP0LpBF5q1asRvy34gKekhUyfNoke7/mzdv1avjSnvz6BZ68b8svR7tJr/lsr/ee5v1GlQk3mLvuHoweNMnzwbjVrNji27GTNxBBYWFrz3zocMHziWsPClz4xz45pwLl+8yuzvphFzL5b3R09h3IhJ/PbnD888b++uA3jZBqJWa6havRITPx1H7Xo1XnC0hBBCCCHEq6BNz354uLnyyzdf4Ohgz6Ur14iIjNLtj7pzh6HjPuDj996lkJMTM76aS/POvbh4ZA/ubq4AXLsRwRvdOhNQrChpaeksXRlGvZBOnNq9haDAAAAuXr5CzeZtKezjzTczPsHT3Y3TZ89zI+K/yYeR73/E/EVLmPjuO1SvUol9h47y3pTpWFtbMaT/G7rj0tLS6Tl4OKPffpPpH07AxbmQgUbLdBW44ss/wA9XNxduXr+ld3ne8AFjqFS1PItXztfdJ1amXClqlG1E+IZtNGvVWHdsuYplmDt/dpa2Pb09+WlRZkHXpHkDNqwJ57svf+bgmR26e7Oibt1m3IhJxMcn4OTk+NQ4tVotf675DUvLzNVqblyL4Ivp36LRaJ76TK7a9WvQ/Y3OBBQvRlTkbb6dPY92TbqzYdcKqtWsmsuREkIIIYQQr4J7MbFcvX6Drz+bQpsWzQBoWLe23jGxcfGE/jaPRvXqAFC/dg0Klw3myx9+5rOP3gfgo/GjdcdrNBqaNqzHoWMnWLB0OdM/zDzm45lzsLAwZ+/GMBwc7AFo0qCe7rzLV68x9+ff+PGLGQzu11u3/2FyMlNmfcngvr11uW56ejrTJo6nW8d2+TEsBVKBu+wwOw8fJnNg72HadwlBrVaTkZFBRkYGgUH++Bb25thh/Uv7mrdunG07DZvW1XsfGOSPl7eH3qIYjxbWiIyI4llq16+pK7wASpQuTnp6ut5lkE/6YMpY+gzoTq261enUrR3rd67Ay9uDWZ9+/cy+hBBCCCHEq8vFuRBFC/vy/iczWLh0ORG3si7N7ujgoCu8Hr1v0qAuB48e0207e/4iHfoMxKNEBVSuhTF3L8r5i5e5cOmK7phtf++hc9vWusLrSVt37QagU9tWupw7IyODJvXrcPvOXW4+EVvrZk1e6rO/al6J4is+Lh61Ws37oz/Gxbyo3uvmjVtE3NT/JXDzcM22HccnZrIsLMyz3QaQkpL6zJgcnfSf1G1hYZGj8x5na2tDs9aNOXH0VI7PEUIIIYQQrxaFQkH4X0soFRTIsPETKVwumKqNWvL3vgO6Y9xcsz6k2MPNlajbdwF48CCRZp16cP1mBHOmTmb3hlUc3raBCmVLk5L6X34aExuHt6fHU2O5FxOLVqvFNbAc5u5Fda+mHXsA6BVfNjbW2NnZvvTnf5UUuMsOs+Po5IhCoWDMByNo3b5Flv0uT/wyvupPDBdCCCGEEK+WoMAAQhf8RHp6OvsOHeGDT2fQpkc/bp3JXFQj+l5slnPuRN/Dy9MdgP2HjxIRGcW6PxdSoWwZ3TEJ9x/g6+2le+/iXIjI23eeGodzIScUCgV7NqzSTS48rsS/946B5NzZKZDFl7mFhd4Mkq2tDdVqVuH82Ut8OLWCESPLW0lJD9m8biuVg1+dzySEEEIIIV6cubk59WvXZMKoYbTt2Z/I27cBSLh/n+1/79Fdephw/z5bd+5m2KB+ACSnpABgYf5fwbTv4GGu3bhJmZJBum1N6tdlxZr1zJw8EXt7uyz9N/63/Zi4ON39ZyLnCmTxVaJUIIt//ZMVS8PwL14MF1dnPv18Em0adaNftyF06t4Op0KOREZEsWPL3/Tq3426DWoZO+xn2rf7IN98/gMhHVpSxM+X25F3+PaLedy5Hc3C0HnGDk8IIYQQQhjJqTP/Y8ykT+jWoS0BxYqScP8Bn305F78ihQko5gdkzkgNfGcsUyaMwcnRkRlfzUWLllFDBgFQo2pl7OxsGTb+AyaMHM6tqNtMnjEbHy9Pvb4mjx/Nus1bqdOqPePfGYqXhzv/O3+Rh8nJjH9nKEGBAQwb1I8+Q0YybsQQqlepRHp6BhcuX2HHnn2ELf7V0MNToBTI4qvPwB4cPXSCcSMmERsTR8++XfhhwVds2rOKzybPZlj/d0lLS8Pb14v6jevgH+hn7JCfy8PLnbS0dD75YAaxMXHY2NpQvVYVvvpxBlWqVTJ2eEIIIYQQwkg83d3x9HDjs6/mcivqNo4O9tStUY3F875BpVIB4OXhwcyPP2DcR1O5fO06ZUoGsTn0Dzzc3QDwcHcj9Nd5jP3oU9r1HkBQQDHmzZnJzH+fFfZI8QB/9m1azfuffsbQsR+Qoc4gKMCfCSOH6Y75ZsanlAgMYN6CxXzy+VfY2dpQIjCALu1CDDcoBZRCq9VqjR1Eds6fP8+ff/7JW2N7Y2P7ej8J+1l+/nIJlStWoWHDhsYORQghhBBCvITQ0FBSE+/Tu33LXJ3Xb9gojhw/xel92/Mpsvyh1Wr55Nv5tGnThsqVKxs7HIN4JVY7FEIIIYQQQghTJ8WXEEIIIYQQQhiAyd7z9ejJ2Gq12siRmDaNWqMbKyGEEEIIUXAplUrUGk2uz1vw3Vd5H4wBqNWZn/V1ymVN9pM6OGQ+pDg2Ot64gZiwh0nJJD9M0Y2VEEIIIYQouOzt7bkXF4+JLsmQ5+7FxQO8Vrms0Wa+0tLSuHr1KomJiWieUuFbW1uzd/th4mITUPAKP6RNAZaWFvj6eWFnn/OngJ87dQmAoKCg5xwphBBCCCFMXcmSJdm/fz+Xb0QQWLSwscPJdyfOXsDc3Jzo6GhiYmLytS+FQoGVlRXFihXD1jbn+Xaex2Ho1Q41Gg0bN27kxIkTZGRkAM+eatRq4VV/OLZWq9V9w1HYz5tWnRs9c4XHh0nJnD11kd1bDlGxYkXatm1rqFCFEEIIIUQ+0Wq1/Pbrr9yNvkur+rUo4V8USwuL559YwCQ8SOTo6bPsPnwChUKBwgDJ/qN8W6FQEBgYSMeOHbGyssr3fp9k8OJr1aow/vnnFJWD6xMYVA57h0IGGXBTptVqSU1J5sa1Cxzct4X09FSsrC2zHRe1WkPyw2QUCgUVK1YkJCTktbpOVgghhBDiVZaamkro8uVcvnIFpVKJtZVl7nM9rTbzyrJ/CxsFCkzlIrKMDDXJKSkoVSrKV61J1TqNDFZ8JT9M4vql8xzevQ1n50IM6N8fMzPDXgho0OIrLi6Ob775htr1W1G6bLChui1Q4mKjWbH0e0qUKIGnp2eW/UqlEgcHB4KCgrCxsTFChEIIIYQQIr/Fx8dz6dIlkpKScn0P2OnTp3mYnEKJ8pVN7kt6lUqFvWMhChcrjoWlpVFiuBMZwZolv9C1a1dKlSpl0L4NWuqdO3cOlcqM4kEVDNltgVLI2Q0v76JoNFoaNGhg7HCEEEIIIYQRODk5UbVq1Vyfl56ezp49e6hapxHlg2vlQ2QFn4e3L85uHpw9e87gxZdBS+H4+HgcHAth/gpeu5qXCrm4k5AQb+wwhBBCCCFEAZOYmIharcbZzcPYoZg0ZzcPo+TbBi2+NBoNKlXWybZ9e3fi46Z45uuLWR/nuJ9bETf4YPwwagUH4O9rRZCfPS0aV2H2zMnE3IvOw0+UP8xUZvJ8MyGEEEIIkWuPVhHPLucG2LQ2jAU/fW/IkJ7q5vVr+NgpWLdqRa7O2/f3TnzsFJw8dkS37ZMPxtKwahmCPO0p4eVAq3rBrA7986ltqFQqo+TbJvGQ5XLlK7Nm4/5s930+40P2791J/YbNc9TWsaMH6dO9JU5Ozgx8cyQlS5cjIz2dI4f3sWjBj1y5fIHvf1qal+ELIYQQQghRIGxaF8apY0foN3iosUPB3dOLNdv34x+Yu8cmlatYmTXb91O8xH+XDCYlJtKz/5sEBpVEoVCwPmwFQ/v3QKPV0KFrz7wO/YWZRPFlb+9Alao1smwP37SG3bu2Mm7Cp1QNrvncdlJSUnhrYBe8vH1ZtW4P9vb/PbCtfsNmvDV0DFs2r31mG8nJyVhbP32ZdyGEEEIIIcTL582WlpZUqZa1BngeeweHLOfN/OZHvfcNmjTnwrn/sXzxApMqvkxr+ZPH3L4dybsjB1CzdgPeGf0BkHkD4acfjyO4YhGK+VhSqYwXfXu14f79BADWrQkl8tZN3p80Q6/wesTOzp4Onf4b/GVLF+DjpuDI4f1079yUwKK2TP14HACRkRGMeLs3ZUu4ElDYmo5t6nHq5NEsbS5buoAm9cvj72tFlXI+zJg2UW8K81Efp08dp3e3lgQWtaV2teKELvs9T8dLCCGEEEKIZxn1Vj9C/1jI+bNn8LFT4GOnYNRb/QA4cnA/XVo1ItDdlpLejgzr35N7d+/qzn10ieCyxQsYN/xNyhRxIaR+NQB87BR8N2cmMz6eSHk/d0r5ODF10ni0Wi27d2yjac2KFPewo2vrxtyKuJmlzccvO6xe2o+J7w5nwbzvqFaqKCW9HRnQvT0x0f/dOpTdZYfZKeTsQnpaWl4MXZ4xiZmvJ2k0Gka83RuAuT/+oVsic+7Xn7FowY988NFMSpQoQ2zsPXbtDCctNRWA/Xt3YmZmRu26jXLV3/AhPenVZzAjRn2AtbUN8fFxdAipg62tHVM/+xZ7e0d+m/8tXTs0Ys/Bi7i6uQMw74c5TJsynjeHjOajKV9w8cJZZk6fiEat5oOPZuj38XYvevZ5k8Fvv8sfi35m9Ih+VKwUTPEgw66wIoQQQgghXk+j3vuQmHvRXL5wjm9/+QMAF1e3zMKrZQMaNWvFDwuX8TApiVmfTqJ/93as3a5/a9CMye/TuEVrvv9tqe7+MoDf5s2lZt0GfDN/EccPH2T2tMmo1Wp2b9/CiHETsbCw4MNx7zB26ECWrgl/ZpzhG9Zw9fJFps35jtiYe0yZMJpJY0fww8Kn38MFmc/yUqvVJCUmsmXjWv7eFs438xe/4GjlD5Msvr77Zib79uzgt8Vr8PT01m0/fuwQ9Rs0o9+A/65Rbd2mk+7nO7cjcXZ2zfK0arVarXs+gkKhQKVS6e3v03cIw955T/d+9szJ3E+IZ/3mQ7pCq069xtStEcSP389m0uRZJCY+4IuZk3l7+HjenzQdgHoNmmJuYcGUj95lyPBxODu76NrsN3C4Lu6qwbXYtnU969f+xagxk15qrIQQQgghhMgJP/8AXFzduHXjut5le2PeHkD5ylWZv3Sl7oHHpcqUo1G1smzbvIHGzVvpji1TviKzv5ufpW1PL2++nb8IyLzkL3zDGn6e+yU7Dp+heMnMyYbbkbeYNHYECfHxODo5PTVOrVbLb8vXYPnvc8Airl/j29nT0Wg0z3xu2e6d2+jRpikAZmZmTP1iLiEdOudwdAzD5C47PHb0ILNnfsSAQSNo1ryN3r5y5SuzfdsGvpj1MSeOH9arth/J7gnZJf0dKeplTlEvc0r6O2bZ37hpa733u3aGU6tOQ5wKOZORkUFGRgYqlYoatepz4vhhAI4c2kdSUiJt2nbRHZORkUHdek1ISU7m/NnTem3Wb9BM97ONrS2+vkWJiorI+cAIIYQQQgiRx5IfPuTwgb2EdOiCWq3W5bT+xYPw9i3MyaOH9Y5v3Lx1tu3UbdhU771/YBAeXt66wgvAv3jmwhpRt56dA9esU19XeAEUL1ma9PR07kXffcZZULlqdTb8fZg/125l0LBRfDh2BEsX/vLMcwzNpGa+Hjy4z7C3ehBUogyTPv48y/53Rk9EqVQSumwhcz6fgourG/0GDGP02I9QKBR4eHqz+++tpKam6v2DrVq7G7VGzR+//8Sqv5ZkadftiecgxMbc49iRAxT1Ms9yrJ9fQOYxsfcAaN64crafJTLypt57R0cnvffm5hakpqRke64QQgghhBCGEB8fh1qt5uP3RvPxe6Oz7I+8pZ/Turpn//ywJ2eyzC0sss1/AVJTn50DOzxxnsW/zwh+Xu5sZ29PhcqZD6au27AxGRkZTHn/Xbr27pflyjdjMani6/1xbxMdfYdFSzfoFU+PWFpaMmb8x4wZ/zFXr1zizyW/8sWsjylS1J/OXftQs3YD/lzyK/v27KBh4xa688qWrwTA1vB12fb75GxZoULOFCvWgnHvf5rlWAuLzLicnJwBmL9gJd4+hbMcV7hIsRx+aiGEEEIIIYzD0dEJhULBiLEf0KJN+yz7nV1c9d5nd5WZqSpfqQrzv/uKmHvRuHt4GjscwISKr9Blv7PqryV88dUvBBYv+dzji/kH8v6k6Sz+fR4XL5wFIKRtF2ZOn8hnU98nuHpt7OzsXyiWOvWasHLFYooXL4WNrW22x1QJrom1jQ1RkRG0bN3hhfoRQgghhBDCkCzMLfRmnmxsbalSvSaXzp+lwuSpRows7x3atwd7B4csBaQxmUTxde3qZSa+N4wKFatSvERpjh45kOUYe3sHZkz7gHIVqlC2XCVsbGzZsnktCfFxutUNraysmPdLKL27taB5o8oMGDSCkqXLoVaruXrlImvCluWoIBv89rus+usPOrWrz8DBI/HxKUJMTDTHjx3Ew9ObwUNG4+joxNj3PmHaJ+OJioqgZq0GqFQqrl+/QvjG1fz8219Y29jk+VgJIYQQQgjxogJLluLPRb8StnwpxQKL4+ziyqSpn9OtdSOGvNGNdp274+hUiKjICP7evoVuvftTq14DY4f9TP87fYrpH75HSIcu+Bbx42FSIls3rmPJwvm8P+UzzMxMouQBTKT4OnhgN0lJiZw8cYS2LbN/mHLNWvVp3LQ1a1cv56fvvyAjI4OAwBLM/fEP6tVvojuucpXqbNl5krnfzGD+T19xO+oWZmbm+AcEEdK2C/0HDn9uPM7OLqzdeIBZn01i+ifvERcXg4urO5Wr1KBFq/9muYYMHYOXlw8//TCHX+d/i7mZOUX9AmjSLATzf69NFUIIIYQQwlT0eGMgJ44cYtLYEcTFxtClV1++mreAVVv2MHvaZN59uz9paWl4+fhSp35j/AICjR3yc7m5eeDg6MSXMz4h+s5t7B0cCQwqyS9LV9E8pJ2xw9Oj0D5ag90A1q9fz7Xrt+jQdbChuiyQDu4N59bNCwwf/vxCUQghhBBCiEdiYmKYO3cuId364VW4qLHDMVl/b15DYuxdBg0aZNB+TW6peSGEEEIIIYR4FRmh+DLYRFuBJSMkhBBCCCFejmSUpsigxZeZmRnp6WmG7LJASk9Lxdw86zPGhBBCCCGEeJZHi0ukp6cbORLTlp6WZpR8+5kLbmi1Wu7cucO9e/dQq9Uv3VlaWhoJ8bH8c/IAlpbWL93eK0mr5ca1C7i4FOLkyZO6zQqFAltbW4oUKSKFmRBCCCGEQKvVcuvWLeLj43W5ukajwdzCgv8dP0RK8kMjR2iitFoirl7C19dHL99+GQqFAhsbG4oWLfrMXP2pC26cPXuWLVu2EBcXlycBibxhYWFB+fLladmyJUql3LInhBBCCPE6On78ODt37uD+/QfGDkU8xtzcnHLlytGqVStUKlWW/dnOfJ0/f57Q0FACAn1p2boGPj7umJtnPVkYjkajJT7+Af87c4U9u4+RnJxMp06dCtRTxoUQQgghxMs7fvw4a9asoWypAKqG1MfT3QUzM8nVjUmj0ZJwP5H/nb/C3/uO8zApia7dumXJ1bMtvnbt2omfnxfdujeT2RUToVKBm1sh6jeogpOTPWtW76JevXq4u7sbOzQhhBBCCGEgWq2WHTt2ULZUAB3bNJQv4k2ESgWuLk7Uq1WZQk4OrFy7g6ioKLy9vfWOy1JZxcfHExV1m8pVSknhZaLKlgvA0tKC//3vf8YORQghhBBCGNCtW7d48OABwZVLS+FlosqU9MfG2oqzZ89m2Zelunp0j5enl0v+RyZeiEqlwt29EPHx8cYORQghhBBCGJAuV3eXXN1UKZVKPNyds107I0vx9WilFGNcN+rvV50RwyfmS9thYZv44fsFuT5vysdfsG/f4bwP6DH9+41CpfDJ8tq0acdTzzEzM8uTFSiFEEIIIUTBYcxc3a9EdYaPyqdcfc0mvp+3INfnfTz1C/btz99c/fM5P1CpRjOcPEth6xJIuaqNmfvDbzxl3UIAzFSqbHP1Zy41b2h/rfqFQoUc86Xt1WGbOHrkFG8P7Zer8z6ZMgc7O1tq1QrOl7ge8fcvyqI/vtXbVqpU8XztUwghhBBCiJxatSz/cvWwtZs4cuwUQ9/ql6vzpkybg52tLbVq5l+uHp+QQLdObSlbpgRWVpZs27GXd8Z8yP0HD/hg/Du5asskiq/k5GSsra2pVKmssUMxGmtrK2rUqGLsMIQQQgghhNCjy9Urvp65+rQpE/TeN2lUjxs3b7Fg0fJcF185XlFjwYJlWJgV4c6daL3tsbFxWFn4MW/eIvbvP0K7tv3w9a6MvW0glSs2ZdGiFXrH79y5D5XCh/Xrt9Kl85s4OZSga5e3gKyXHeamvS1b/qZXz2E42gdRrGg1Pp/1ve6Y/v1G8fvCUM6cOa+7pK9/v1HP/cwqhQ8A48d9qjtv5859AMz54keqB7eikGNJPN3L0ybkDS5cuJyljXnzFlGsaDXsbAJo1rQ7x4+fRqXwYcGCZc/tXwghhBBCiJxYsGgZZnbZ5+oWDn7Mm7+I/QeO0LZzP7yLVcbWJZCK1ZuyaMkTufXf+1BY+7B+41Y693gTB/cSdOmVmas/edlhbtrbsu1vevYdhr1bEEWDqjHri/9y9X5vjmLh4lDO/O88CmsfFNY+9Htz1HM/s8I6M1cf98GnuvN2/p2Zq3/x1Y8E126Fo0dJ3IuUJ6TjG1y4mE2uPn8RRYOqYeMcQNPW3Tl+4jQKax8WLHp2ru7iXIi0tPTnxvikHM98dejQkqFD3mdF6DqGDe+v2/7XXxsA6NIlhPDwv6ldO5i3hvTBysqSfXsP8+bAsWg0Gvr27arX3pDB79Grd0eGrHoj2weQAVy/fivH7Q0dMoHefTrx16r5rA7bzIT3plGufClatGjIpA9HER0dw/lzl3WX9rm5Pf8mxb3711C7ZluGjxhAj57tAShdOgiAiIgohg7vR9Givty/n8i8HxdRp1Y7zl3YjbNzIQDWrAln6JAJDBzUk06dW3PyxBm6d30r274uXbpGIceSJCenUK5cSSZ+OIr27Vs8N0YhhBBCCCE6tG3JkBHvE7pyHcPffixXD/s3V+8YQvjWv6ldI5ghgzJz6737DzNwyL+5dW/93HrwsPfo3aMjqwY/I1e/cSvH7Q0ZMYE+PTuxatl8wtZs5r1J0yhfrhQtmjXkw/dHEX0vhnPnL/PHbznP1ffvXEPNBm0Z8fYAenZrD0DpUv/m6reiGP52P4oWyczVf5y/iFoN23Hh1GO5+rpwhoyYwKD+PencoTUnTp6ha+/sc3WAjIwMkpNT+HvPAX5fsoLJH4x+boxPynHx5ejoQMtWjfhzaZhe8bVsaRhNm9XD2bkQ3bu3023XarXUq1eDiIgofp63OEux1KZtU2bMfPYNe7lpr2OnVkz+eAwAjRvXZcP6bfy1Yj0tWjQkIMAPNzcXbly/latL+x4dW6SIT5bz5nw5RfezWq2madO6eLpXYMWK9Qwe3BuA6VO/plGj2vz08+cANG/egPT0dD768HO9tipVKktwcEVKlwkiPv4+P/7wO506DGRZ6Dw6dw7JcbxCCCGEEOL15OjoQKvmjVi6PEyv+Fq6PIxmTf7N1bs+kVvXqUHErSjmzV+cpVhqG9KUmdOek6vnor1O7Vvx8aR/c/WGdVm/aRsrVq2nRbOGBPj74ebqwvUbt6hRPRe5+r/HFinsk+W8Lz9/IldvXBf3IhVYsWo9gwdm5upTZ3xNowa1+fn7f3P1pg1Iz0jnwyn6uTrApctXKV62ju79pAkjGf3O4BzH+kiuHuTVvUc79u8/yo0btwCIirrDrl0H6N6jPQBxcfGMfOdDihWthqV5USzNi/LzT39w4cKVLG21at34uf3lpr2mzerpflYoFJQsVZyIiKjcfLxcOXDgKM2adsfNpQwWZkWwswkkMTGJi//GplarOX78NG3aNtM7r2275lnaemfkIIYO60eDBrVo374FGzYupnr1Snz80ex8i18IIYQQQrxaenRtx/6DT+Tquw/Qo0t7IDO3fufdDykaVA1z+6KY2xflp1/+4MKlrLl16xY5y9Vz2l6zJvq5eqmSxYm4lY+5+sGjNG3dHRefMpjZFcHGOTNXv3DxsVz95GnattbP1duFZM3VAQr7enN4zwZ2bA7lk4/GMuebn5j8ae5z9VwVXyEhTbC1tWHZn6sBCF2+FisrS93lcf37jebPpWGMGTuETeFLOHh4A/0HdCclJTVLWx4ebs/tLzftOTnpr7xiYWFOajbH5YUbN27RollPNGo1P8ybye69YRw8vAF3d1ddbNHRMWRkZOD6xJSpu7vrc9tXKpV07NSas2cvkpycnC+fQQghhBBCvFpCWmXm6n+GZubqy//6N1dvm5mr9xs8mqWhYYwdNYTwdUs4vGcDA/o+JVd3f36unpv2nByfyNXNzbM9Li/cuHGLZm16olarmfftTPZuD+PwnuxzdTfXJ3J1t+xzdUtLS6pWqUCDerX48P3RTJ8ygWkzv+H27bu5ii1Xqx1aW1vTrn1zlv25mnHjh7Lsz9WEtGmKra0NKSkprF+3lS/mTGb4iAG6czQaTbZtPe+J3Lltz5A2bdpBYmISK1bO1xV9GRkZxMbG645xc3PBzMyMe9ExeufevXvPkKEKIYQQQojXhLW1Ne3bNOfPFasZP2Yof4aupk2r/3L1dRu2MmfmZEYMfSy3nvfiuXpu2jOkTVsyc/WVfz4/V4++90SuHp2zXL1KpfKo1WquXb+Jp6d7jmPL1cwXQPce7Tl+/DSbN+/kwIFjuksOU1PT0Gg0mFuY64598CCRtWvCc9tFvrRnYWHxQtW1ubk5KSkpettSklNQKBSYm/8X2/Lla8nIyNC9V6lUVKpUljWrN+uduzps03P71Gg0rAhdR5kyJbC2ts51zEIIIYQQ4vXUo2t7jp84zeYtOzlw6Bg9urYH/sutLZ7Irdesf7lcPa/ae6lcPVU/V0/OLldfkU2uXqEsq9fp5+pha5+fqwPs2XcIhUJBMb8iuYo318/5atq0Hi4uhRg0YAxOTo60bNkQyLzJLzi4IrNmfKerJGfNmIujo8MLzfbkdXulSgXy269/snRpGMWLF8PV1Rk/v8I5Om/N6nDq1K2Ora0NJUoE0LBRbQAG9B/N4Ld6878zF5jzxbwslz5+MGkkHdr1Z/Cb4+jcJYQTx0/z+8LM5TeVysy69/r1CPr3HUW3Hu0IDPQjLi6BeT/8zpEjJwn96+dcf04hhBBCCPH6ato4M1cf8Na/uXrzx3L1KhWZMfs73Fwzc+sZs+fi6OCQ49mex+V1e6VKBPLrwj9ZuiyM4oH/5upFc5Crlwxk9dpw6taujq2NDSWCAmjUIDNX7z94NG8N6s2Z/13gi6+z5uqTJoykXZf+vDl0HF06hnD8xGkWLtbP1RMS7tOqfR969+hIYEAx0tPT2fn3fr7+7hfeGtQ7R7dSPS7XM1/m5uZ06tyayMjbdOzUCgsLC92+xUvmEhjoR/++oxj1zod06hxCnzc657aLfGlvwMAedO4SwsgRk6ge3IopH3+Ro/O+/W46Go2G1i17Uz24FUePnqJcuVL8uuBLjh39h7Yh/fhzaRjLV/yEo6O93rlt2zbjux8+I3zzTjq0G8CmjTv47ofPgMxfWAB7e1scHO2ZPvVrQlq9wcD+76LRaFi/cTEdOrR8oc8qhBBCCCFeT+bm5nTu0JrIqNt0aq+fqy9ZOJfAAD/6vjmKd8Z8SOcOIbzR68Vz9bxsb2C/HnTpGMKIMZMIrtOKj6fmLFf/7svpaLQaWrbrTXCdVhw9fopyZUux4OcvOXr8H0I69mPp8jBWLPkJR4cncvWQZvzwzWds3rKTdl0GsDF8Bz9882+u7pCZq1tZWRJU3J853/xEuy796TPwHXbtPsCP385g7pfTcv05FVqtVvv4hkuXLvHHH38wcnQPHBzsct2geLZfflnK4EFjuXz1QI5m3p5m8e8bsLUtRKdOnfIwOiGEEEIIYcpOnDjB6tWr+XDcQN3sjMg7vyxYyqC3x3L13IEczbw9zZLQTagsbOnWrZve9lxfdihyLjY2jk+mfEnDRrWxt7flyOGTTJ/2DW3bNX+pwksIIYQQQgjxcmJj45gy7UsaNcjM1Q8fOcm0Wd/QLqT5SxVez/JaF1+P33T3JIVC8dSneeeUubk5ly9fY+mSVcTH38fNzYXefTo99+HSQgghhBBCvO4MkqtfvcaS5f/m6q4u9OnZiZlT8y9Xf22Lr2vXbhJQrMZT99evX5PtO1e8VB/29nasXff7S7UhhBBCCCHE6+ba9ZsUK/mMXL1uTXaGv3yuvm6lYXP1LMXXowoyI0Nt0EAMzdvbg4OHNzx1v729rQGjyb2MjIyXrvaFEEIIIUTB8niubmHx6t7z5e3lweE9BThXV6uxyCZXz1J8OTk5AXDndizOzo5P7n5lWFhYULVqBWOH8ULUajXR0XEUK1bC2KEIIYQQQggD0uXqd2Mp7Oth3GDykYWFBVWrFMxcXaPRcDc6joq+xbLsy1IuFypUCE9PD44dO4tGY/wnVIuszpy5QkpKGqdPn2bv3r3y7ySEEEII8Zrw9fXFzs6OIyf+xxOLlgsTcfbCNZIeJlOqVKks+7Kdq6xXrz5Xr0QSunwr165Gkp7+9JvdhGFoNBpiYxPYs/sEa1f/jaurF46Onmzbto2wsDApwIQQQgghXgMKhYL69etz6swlVm/Yxc1bd17524UKAo1GQ1z8ffYePEnY+p0UL14cb2/vLMdlec7XI2fOnGHLlnASEu7ne7Ai51QqMwIDy1G7dmtUKhWXL59h+/YVdOvWjZIlSxo7PCGEEEIIYQBHjhxh586dJCUlGTsU8RiVSkXZsmUJCQnBzCzr2oZPLb4AtFotkZGRxMTEPHOpR0P79ttvGTFiRJafn+bNN98E4Oeff35ue8+yY8cOrK0dKVmyci4jfnkKhQIrK1u8vf0wN7fQ2/fXXz/g7e0uD1wWQgghhHiNaDQaIiIiiIuLQ61+udmvnObDL2rMmDF88cUX2e57Xq5eECgUCmxsbPDz88PS0vKpxz1zqXmFQoGPjw8+Pj55HuDLiImJoXLlyll+fppjx44BPPW4nLShVqtZu3YtlSo1NErx9SxeXn7cvn3V2GEIIYQQQggDUiqVFClShCJFirx0WznJh1/G2bNnn9r+83L1V8mruz5lHns082durl/J+vgonvtatmwBly6dZ+LE4TRoUJqAABtq1CjGhAlvExt776l9NmlSAR8fBQcP7n5mbObmFqSnp738hxRCCCGEEELkm9f2IcsvSqHQf79mzX69923b1mTAgBG0b99Tt83PL4A1a5Zx8OBuevd+i9KlKxARcZ3Zsz9i//6dhIefyDI9ef78Gc6ePQXAqlVLqF697rOieqnPJIQQQgghhMh/Uny9pCpVsj5528enSJbt7dr1oF+/YSgeq96KFStO+/a12bp1Ha1b69+vtXLlHyiVSmrUqM+6daF8+uk3mJub58+HEEIIIYQQQuQ7uezQQJydXfQKL4CyZSsBcOdOpN52rVbL6tVLqV27EYMHv0tcXAw7dmwyWKxCCCGEEEKIvPdaFV+m9iysw4f3ABAYqP8AtiNH9nHz5jU6dOhJgwbNKVTIhbCwJcYIUQghhBBCCJFHXqviKzU11dgh6KSkpPDJJ2MpW7YSdes21tu3atUSrKysaNmyI+bm5rRu3Znw8DUkJSUaKVohhBBCCCHyl6lNlOSH16r4SklJMXYIOhMmDOHmzat8/fXvepcjZmRksG5dKI0atcLBwRGADh16kpz8kI0bVxkrXCGEEEIIIfKVKU2U5Bcpvoxg5sxJrFr1B/PmhVKyZFm9fbt2hRMTE02TJm1ISIgnISGekiXL4eHhJZceCiGEEEKIV5ap5Or56bVa7dAUVgv89ddv+fbb6Xz55QIaNGieZf+jAuvdd/vz7rv99fbFxERz795dXF3dDRKrEEIIIYQQhmIKuXp+e62KLysrK6P2Hxa2lI8+Gsn7739Gly5vZNmfnPyQzZtX06JFewYOHKm3Lzr6NkOH9mDNmmUMGDDCUCELIYQQQghhEMbO1Q1Bii8D2b9/F6NG9aV27UbUqFGfo0cP6PZ5efni7e3L5s2rSUpKZMCAd6hVq0GWNr7/fharVi2R4ksIIYQQQrxyzMxe/dLktbrny5j/oPv27SA9PZ09e7bRtm1NvdfSpfOBzFUOfXyKZFt4AXTp0pdjxw5w7dplA0YuhBBCCCGEyAuvfnlpYLduabPdPmbMx4wZ8/Ezz124cO0z9w8aNJJBg0Y+8xghhBBCCCGEaXqtZr6EEEIIIYQQwlik+BJCCCGEEEIIA5DiK4eUysyhUqvVRo4kK40mA6VSZewwhBBCCCGEEM8gxVcOmZmZYW1tQ2zsHWOHkkVsbDSOjo7GDkMIIYQQQgjxDK/FghutW7ciJubeS7WhUCgoVaokZ8+eolKlelhYWOZRdC/n/v04IiOv0Lx51gc2CyGEEEIIYeoaNWyARqMxdhgG8VoUX6tXr0GrzX4VwtyoVq0aR44cZePGxVSuXB8fn2JGu9wvPT2NGzcucOjQVhwdHSlbtqxR4hBCCCGEEOJlbA7fYuwQDKZAFV+3bt3CwcEBgKSkJB4+fAjAvXv3SElJeeF2k5OTiYiI0L13cXHB2to6y3EeHh6cP38OFxcXNm36AwAzM3MUCsUL9atWq1Gpcl+8abVaMjLSAfDx8aVLl87Y2Ni8UAxCCCGEEEI8LjExkfj4+DxtU61W6+Xb2XFycsLOzi5P+zU1Bar48vX1pZCLM7Vq1KRM+fJcv3KFihUr4ubm9tJtFy5cWPezq5sn69aGUb169SzHPXz4kBEjhnPnzh0iIyP59NOp3L4d9VJ9K5UqunfvRnBwcI7PsbKyomjRohQqVOil+hZCCCGEEOJxAcWDuPuS+W12Hs+3s2NrZ8+yP5fSunXrPO/bVBSo4gsgLiaW+Ph4bkdl/kLExsZm7vD/GSx8X74DrYbY29OoW68+C377lZ49e2Y5RKFQ4OnpiaenJ/v27efuXTeg/Et0eoZ9+95h/PjxfPbZZ7qVFYUQQgghhDC02Jh7UG8YlA0xaL8P9/xAmzZtmD17NqNHj37hq8tMWYEqvux8XEhPTOHevXtYWluTmpxMUlJS5k77emAdlCf9aBwbo7k6mF69enHmzP/49NNPnlMQuQDFX6LHQMCdWbM+5+zZcyxZ8scrP+UqhBBCCCFMmFcZKN3CoF1qSzaFNR8wZswYTp8+w48//oCFhYVBY8hvBWqKRaFS4lo/iHuxMbptiamped+R0hL8F0CRmUz/bDodO3b+r8jLFwqgFtCd9es3U716Ta5fv56P/QkhhBBCCGFilCpoPxN6L2DhosU0bNyE6OhoY0eVpwpU8QXg1rAkcbFxAChLliA1MTF/OlIowHs8FA9j7fpwatSsw82bN/OnL50SaDQDOH8+ksqVq7Jv37587k8IIYQQQggTU6Mvmnd2cPDUOSpXrcbp06eNHVGeKXDFl3uj0mjUajIyMlBWKJf/HTq3RVNqH2cvxVKpcjBxcXH53KEHavVA4uPtqF+/AYsWLcrn/oQQQgghhDAx/rVQjz1MlNqe6jVrsX79emNHlCcKXPHlWNYHCydb1Go1yhLFwdw8/zu1LY+69GHi0gPYu+8AS5Ysye8O0Wh6k5FRhjfeeIMJEya8Ng+eE0IIIYQQAgDnoqhH7yXZvyFt2rRhzpw5efLsXmMqcMWXQqnErVEpMjIywNwcVakShunY3B1Nie1onXvQq1cvJk6clM8FkRnQDmjGzJmzaN++A4n5dYmlEEIIIYQQpsjKHu2glWgbj2PMmDEMHDiItLQ0Y0f1wgpc8QXg0bgUarUaAEXlimBmoEUblZYQsFBvIY78rb5lIQ4hhBBCCPGae2IhjkZNmnLv3j1jR/VCCmTx5d6wFJBZ9KjKlwO1GjQGqoCfWIgjNja/7wGDJxfiOHz4sAH6FEIIIYQQwoT8uxDHgZNnqVy1mrGjeSEFsviy8XMlc1YIFF6eoNWC2hBF0GMKhaBx6oxanWGgDj1Qq2sQG3uP33//3UB9CiGEEEIIYUL8a6GuOZib168aO5IXUqAesvxIzP5LPJr50vzvLKhUYOZuuADUSSiuvAGxq7CxtSdfHwEGgAbYBeyiZ89efP755/ndoRBCCCGEEKZFnQ4rRsLuHxg5cpSxo3khBXLm6+6OsygUmaFrTpzKnPlSKAzTeepNVOfqYpUcTlhYGLa2NvncYRqwAtjF9OnTWbx4EVZWVvncpxBCCCGEECYkKRblDy1R7f+Zn3/+ma+++tLYEb2Qgjnztf0cKpUKAMXJf8BQy7A/OIjiQhs83K3ZuGEv5cuX580338rHDu+jUi3D3DyWJUtW0qFDh3zsSwghhBBCCBN05zyqn9pglx7L6q1bqV+/vrEjemEFrvhKT0zh3qErWFhZo46PJ/3GTcN0fG8pyqv9cXCw4/ixw7i75/dljrdQqZbh7m7Phg37qFixYj73J4QQQgghhIk5uwXVgq4EFPZm4/pD+Pv7Gzuil1LgLju8t+cCmgw1ZmZmaE6dzv8OtRq4OQku9aRnj27UqlndAIXXaZTKBVSqVIJjx45I4SWEEEIIIV4/u75D8UNLGtetyaED+wp84QUFsPi6u+Mc7l4eKJVKNCf/wdzKMv86UyehuNQFReR0Zs6cye+/L9Bd7pg/NMAOYAXdunVh9+6/8fT0zMf+hBBCCCGEMDHqdFg2FEKHM/Kdd9iwbi2Ojo7GjipPFLjiK2b7ORo3bAyA9m40/v4B+dNR6k1U5+pg+XAzYWFhjB8/HkW+Lurx38Ia06ZN448/FsvCGkIIIYQQ4vXyxMIaX345J58nPwyrQBVfqfFJxBy7RuNGjXTbShYvnvcdPTiI6mwwHg6xHDywj7Zt2+Z9H3ruo1ItxMrqKitXruSDDz7I50JPCCGEEEIIE3PnPKo5NbC/e4JtW7cyaNAgY0eU5wrUghtp9x8C0LBhQ5ITEwEICgrK3HltJKjs86AXDcqEdQQHV2X16pU5vL/rDHD3hXtUqSJkYQ0hhBBCCGE69v4MF3YYtEvVhS2vzMIaT1Ogiq9mzZsREBBAsWLFGD1qFGfOnqV3795cuHiJpMQkICFP+ilffjjTp0/D0vL595O99dab7N9/4KX6c3evyBdffCH3dwkhhBBCCKN7e8jbnD13jrzKrQFOnDjx3EkGvypdmD3781fm/q7sKLRardbYQRQkbdu2Zc2aNcYOQwghhBBCiDyT3zmu5NCZCtTM1+eff46xa8VLly4xa9asPG3Tzc2NN954I0c3E544cYLw8PA87V8IIYQQQhR8Hh4e9O7dO0c55fHjx9myZQsAtWvX1tv3119/cfny5TyNLSc5dNGiRenatesrvfZBgZr5UiqV2BUy7gqAarU6z1dcSYxLoUXL5vy5dBkODg7PPLZz586sXrkS+1do1RchhBBCCPHy4jMyaBMSwuIlS7C3f/ZaCO3bt2f92jVYqBTYOjhStVoNNmzYAIClpRUZGSqUyrx7pFNOcuiMjDh6936D+fN/ytHtPwVRgZr5ci1sx+fXWxo7jDz3z+bb/NhtBzVrVWftmvXPvMFQo9FQT6tlSUaGASMUQgghhBCmbiswZONGalevztqNGylatOhTj9VoNLT00PJVRS1Bm+O5ceOG3j6NZgoaTb88je/56esqliwZxcWLl1izZlUOF74rWArUUvOvqnLNPZl4oCExKZEEV6vK33//beyQhBBCCCFEAdMEWKdWk3DhAsGVKrF3797nnuNvBz0Ka7l2+RJpaWn5H+QzdUCjWcWRI5eoVCmYU6dOGTmevCfFl4nwLunAhwcb4lnegiZNGvPLL78YOyQhhBBCCFHAlAQ2qtUEJiTQqEEDFi5c+Nxz3i8JSSmpLFq0KP8DfK7KqNWbuHPHkerVa7F69WpjB5SnpPgyIXYulry7uS61BxRh0KBBvDvmXdRqtbHDEkIIIYQQBYgLsEyjoXNGBv369WP8+PHPzClLO0AHH5gxbSoZJnFrizdqdRipqQ3o0KEDM2fONPqie3lFii8TY2aupO8PVej1TSW+/uor2rQN4f79+8YOSwghhBBCFCAWwBfAFOCL2bPp2L49Dx48eOrxk0rBpavXCA0NNVSIz2GLVvszWu0oJkyYwBtv9CM1NdXYQb00Kb5MkEKhoOmI4ozeUJe/92YuxHHlyhVjhyWEEAZjPXgwrteu4Z6UhFN4OGYVK+Kh1WLVt6/uGKu+fXE+eRL35GRcIyKwnToVlEq9/R5aLWYVK+K0YQPuiYm4XLiAVZ8+WfqzaNUK5wMHcH/4ELe7d7H//nuwsTHIZxVCiPyiAN4CFmm1bP93IY7r169ne2zlQtDCS8H0Tz8xoVkmJfAe8ANLliyjfv1G3L1719hBvRQpvkzY4wtxVKsezO7du40dkhBC5DvLNm1wmDeP1PBw4jt0IG3rVhyXL9c7xmb0aBzmzydt82bi27QhaeZMbN55B7tp07K05/jHH6SFhxPfvj0Zx4/jsGABqpIl/+uvUyec1qwh459/iO/QgQfjx2PZsSOOcu+tEOIV0ZicLcQxsaSW02fPmVDx9cirsxCHFF8m7tFCHO5lzWncuJHMgAkhXnm2kyaRum0bDwYPJi08nIezZpHy2A3jCjs7bKdM4eGsWSSOH0/a1q0kf/stiePGYfPOOyicnfXaezh3Lg+/+oq0rVtJ6N8fkpOx6tRJt99+9mxSli3j/ptvkrZ5MykLFnC/b18su3ZFVbq0wT63EELkp0cLcQT8uxBHdjNgdVyhnrsSpVYNmFoBpr8Qx/InvpQrKApU8ZWWnm7sEIzCzsWSMeF1qd2/CCdPnuTG808RQoiCSanErFIlUtes0duc+thqV+a1aqG0tyclNBRUKt0rbetWFDY2mJUtq3duWnj4f28ePkR9/TpKX18AVEFBqPz8SF2+XL+tXbtAo8G8atX8+6xCCGFgLsByjYaOGRlPnT36oISGDC3AVUOGlkOZC3GkpNSjW7duxg7mhRSohyy/zjJS1cRHJQPyjyaEeHUp3dxQmJujjY7W26557Bp/pasrAC7Hj2fbhqpwYR7/qk4TH6+3X5uWhsLKSq8tp7Cwp7YlhBCvkjQg+hn7k3SLHZrnfzAvJAWIIfOOtoKnQOXxFuam+kuQv2JuJPFN2/3cu5xCtWrVsDt0yNghCSFEvtBER6NNT0fh5qa3Xenu/t8xsbEAxHfogPrmzSxtqK/m/NvaR23dHzaM9IMHs+6PjMxxW0IIYepuAG+oVNyytKRq6dJw64jefq0Wpl1QYabUkKExxS+fzqNSvYGDQxLLl4c//3ATVKAuO3wdXdp/j0+Cd0CCLfv3HcDHx8fYIQkhRP7RaMg4fhyrdu30Nlu2b6/7OX3/frRJSSh9fck4ejTLS/tvQZUT6nPnUN+8icrfP9u2NFFRefXJhBDCqA4CrVQq0nx8OHD4MF5eXlmOCb8Dx2LUaEyyRNiGShVCUJADx44dpkmTJsYO6IUUqJmv182+xdf5beARgoOrsWplGO6PffMrhBCvqqSpU3Faswb7n34iNTQUs0qV/ltiXqNBm5BA4kcfYT9rFipfX9J27gS1GpW/P5bt2hHfqRMkJ+e4vwfvvovjkiUobG1JXb8ebVISqqJFsWzdmsQPPkB98WL+fFAhhDCQZcA4pZKaNWuyYtUqXP+95PpJ084rqVa1IkePm9JqglrgJxSKKbRo0ZqlS//A3t7e2EG9MFMsa197Go2WFR+c4qc+B+nZszfbt+2QwksI8dpIXbuW+0OGYNm8OU6rV2PZsiUP3n4bAG1CAgAP58whoX9/LBo2xOmvv3AMDcV68GDSDx+GtLTc9bdiBfGtWmFWsiSOS5fitGYNtmPGoL52Dc2dO3n++YQQwlDUwKfASKBP//6Eb9v21MJrdzTsvqth4oeTUShM5X6qNGAsMJmxY8ewevWqAl14gcx8mZyUxHR+7nOYY6tv8fnnnzNmzBgT+g9ACCEMI3nePJLnzdO9txowAID0kyd121KXLSN12bKntpGycKHeEvWPxFaqlGVb2tatpG3d+jIhCyGESUkEhioUbAW+nDOHkSNHPjOnnHZeSbnSJQgJCTFYjM8Wg1L5JkrlEebPX0DfR1dAFHBSfJmQxxfWWLNmjQn98gshhOEoChXCbvJk0rZvR/PgAebBwdhOnEhKWBiabJ5LI4QQQt+jhTUiraxYFxpKy5Ytn3n8kVjYHKVh6ZyPUCpN4cK4RwtrPGTt2u3Url3b2AHlGSm+TMSl/feY2+EAjtYu7N+3k3Llyhk7JCGEMI70dFQBATj07InCyQlNdDQpixaR+N57xo5MCCFM3kFgoEqFg48P+zdupHQOHhY//byC4v5+dOnSJf8DfK5tqFRDCAryY8OGnfj5+Rk7oDxVoIovrRYy0jXGDiPPHVp2k98GHSW4anCOFtbQAq/n46aFEK+FxESi27QxdhRCCFHgrCRnC2s8ogFOxcOqCC2//DIJlUr12F41hs84f0Oh+PiVWFjjqbQFCJl1xyv56tuvrzYlJeW5Y9C1a1ejxyoveclLXvKSl7zkJS/TfA0cOFCbmpr63JyyY8eOunMKe3vpnWNlZWO0+MeNG6fNyMh4qZrBlCm0Wq2WAuL333+nAIWbY25ubrRs2TJHC2tcuXKF3bt3GyAqIYQQQghRkHh4eNC8efMc5ZSXLl1i7969AFStWpUyZcro9u3atYtr167lV5hPVbRoURo0aGDwfg2pwBRfGRkZ7Nmzhzp16mBmVqCulixQZJwNQ8Y5/8kYG4aMc/6TMTYMGef8J2NsGDLOps0UljPJEbVaza5du1Cr1cYO5ZUm42wYMs75T8bYMGSc85+MsWHIOOc/GWPDkHE2bQWm+BJCCCGEEEKIgkyKLyGEEEIIIYQwACm+hBBCCCGEEMIACkzxpVKpqF+//hPPHxB5TcbZMGSc85+MsWHIOOc/GWPDkHHOfzLGhiHjbNoKzGqHQgghhBBCCFGQFZiZLyGEEEIIIYQoyKT4EkIIIYQQQggDkOJLCCGEEEIIIQxAii8hhBBCCCGEMAAzYweQE9evX2ffvn1ERkaSmJhIt27dKFmypLHDKjB2797NuXPnuHfvHmZmZhQuXJgmTZrg6uqqOyYjI4PNmzdz5swZMjIyCAwMpFWrVtjZ2emO2bhxIzdv3uTu3bu4uroyZMgQY3ycAmHPnj1s27aN6tWr06JFC0DGOK/cv3+frVu3cunSJdLT03F2dqZdu3Z4e3sDoNVq2blzJ8eOHSMlJYXChQvTunVrXFxcdG38/fffXLx4kdu3b6NSqZgwYYKxPo7J0Wg07Ny5k3/++YfExETs7e2pUKEC9erVQ6FQADLGL+J5f8dyMqbJycls3LiR8+fPo1AoKFWqFC1btsTCwgLI/H/MunXriIqKIjo6mqCgILp3727wz2pMzxpntVrN9u3buXTpEnFxcVhaWuLv70+TJk2wt7fXtSHj/Gy5ycnWrVvH0aNHad68OTVq1NBtlzF+vpyMc3R0NFu3buX69etoNBrc3Nzo2rUrjo6OgOQdpqpAzHylpaXh4eFBq1atjB1KgXT9+nWCg4MZOHAgffr0QaPRsHjxYtLS0nTHbNq0iQsXLtClSxf69evHgwcPWL58eZa2KlasSJkyZQwZfoFz69Ytjh49ioeHh952GeOXl5yczK+//opKpaJXr14MHTqUZs2aYWVlpTtm7969HDx4kNatWzNo0CAsLCxYvHgxGRkZumPUajWlS5ematWqxvgYJm3v3r0cOXKEli1bMmzYMJo0acK+ffs4dOiQ3jEyxrnzvL9jORnTlStXcvfuXfr06UPPnj25ceMGa9eu1e3XaDSYmZlRrVo1/P398/0zmaJnjXN6ejq3b9+mXr16DB48mG7duhETE8PSpUv1jpNxfrac5mRnz54lIiJCr7B9RMb4+Z43zrGxsfz222+4urrSt29fhgwZQr169TAz+29eRfIO01Qgiq/ixYvTqFEjSpUqZexQCqTevXtTsWJF3N3d8fT0pF27diQkJBAVFQVASkoKx48fp3nz5hQrVgxvb2/atWvHzZs3iYiI0LXTsmVLqlWrRqFChYz1UUxeWloaK1eupE2bNnoFgYxx3ti7dy+Ojo60a9cOHx8fChUqREBAAM7OzkDm7MHBgwepV68eJUuWxMPDg/bt2/PgwQPOnTuna6dhw4bUrFkzS4Es4ObNm5QoUYKgoCCcnJwoXbo0AQEB3Lp1C5AxflHP+juWkzGNjo7m0qVLtG3bFl9fX4oUKULLli05ffo0Dx48AMDCwoKQkBCqVKmi98326+RZ42xlZUWfPn0oU6YMrq6u+Pr60rJlS6KiokhISABknHMiJznZ/fv32bhxIx07dkSp1E81ZYxz5nnjvH37dooXL07Tpk3x8vLC2dmZEiVKYGtrC0jeYcoKRPEl8lZqaioA1tbWAERFRaHRaPS+XXJ1dcXR0ZGbN28aJcaCasOGDRQvXjzLN3Uyxnnj/PnzeHl5ERoayueff868efM4evSobn98fDyJiYl642xlZYWvr6+Mcw4VLlyYq1evEhMTA8Dt27e5ceMGgYGBgIxxfsjJmEZERGBlZaW7vBbA398fhUKhl0iJ3Hn09/DRl2Uyzi9Pq9WyatUqatWqhbu7e5b9MsYvT6vVcvHiRZydnVm8eDGff/458+fP1/sCTPIO01Ug7vkSeUer1bJp0yYKFy6s+59iYmIiKpVKb6YGwNbWlsTERGOEWSCdPn2aqKgo3nzzzSz7ZIzzRlxcHEeOHKFmzZrUqVOHyMhINm3ahEqlomLFirqxfPTN3yO2trYkJSUZI+QCp06dOqSmpjJ37lyUSiUajYZGjRpRvnx5ABnjfJCTMU1MTMyyX6lUYm1tLf8PeUEZGRls3bqVcuXKYWlpCcg454U9e/agVCqpXr16tvtljF9eUlISaWlp7N27l4YNG9KkSRMuXbrEsmXL6Nu3L35+fpJ3mDApvl4z69ev5+7duwwYMMDYobxSEhIS2LRpE3369NG73lrkLa1Wi7e3N40bNwbAy8uLu3fvcvToUSpWrGjc4F4RZ86c4Z9//qFTp064ublx+/ZtNm/ejL29vYyxeGWo1WpCQ0PRarW0bt3a2OG8MiIjIzl48CBvvfWWboEekfe0Wi0AJUqUoGbNmgB4enpy8+ZNjh49ip+fnxGjE88jWeJrZMOGDVy8eJF+/frh4OCg225nZ4darSYlJUXvG5KkpKTX9lrr3IqKiiIpKYl58+bptmm1Wq5fv86hQ4fo3bu3jHEesLe3x83NTW+bq6srZ8+eBdCNZVJSkt5N3klJSXLvUQ5t2bKF2rVrU7ZsWQA8PDxISEhgz549VKxYUcY4H+RkTO3s7LLMLGo0GpKTk+X/IbmkVqtZsWIFCQkJvPHGG7pZL5Bxflk3btwgKSmJL7/8UrdNq9USHh7OgQMHGDVqlIxxHrCxsUGpVGb79/DRJYWS25kuKb5eA1qtlo0bN3Lu3Dn69u2b5aZKLy8vlEolV65coXTp0gDcu3ePhIQEChcubIyQC5xixYrx9ttv621bvXo1rq6u1K5dGwcHBxnjPFC4cGHdvUiPxMTE6JbVdXJyws7OjitXruDp6Qlk3tMREREhq+7lUHp6epZvrBUKhe6bVhnjvJeTMfX19SUlJYXIyEjdvTJXr15Fq9Xi6+trtNgLmkeFV0xMDH379sXGxkZvv4zzyylfvnyWe54XL15M+fLldTPnMsYvT6VS4e3tneXvYWxsrO7voeR2pqtAFF9paWnExsbq3sfFxXH79m2sra11v2Ti6TZs2MA///xD9+7dsbS01F3ra2lpibm5OVZWVlSqVInw8HCsra2xtLRk48aN+Pr66v2PMDY2lrS0NBITE8nIyOD27dsAuLm5oVKpjPLZTIWlpWWWG4vNzc2xtrbWbZcxfnk1atTg119/Zffu3ZQpU4Zbt25x7NgxQkJCgMwioXr16uzevRsXFxecnJzYsWMH9vb2es9HSUhIIDk5mYSEBLRarW6cnZ2ddc+ZeV0FBQWxe/duHB0dcXd3JyoqigMHDugSJxnjF/O8v2PPG1M3NzcCAwNZu3YtISEhqNVqNmzYQNmyZfVmy6Kjo1Gr1SQnJ5OWlqYb90dF3avuWeNsZ2dHaGgoUVFR9OjRA61Wq/t7aG1tjUqlknHOgef9Lj9Z0CqVSuzs7HTPFpUxzpnnjXOtWrVYsWIFRYoUoVixYly6dInz58/Tr18/AMntTJhC++jrTBN27do1Fi5cmGV7hQoVaN++veEDKmCmTJmS7fZ27drpEqpHD+I7ffo0arWagIAAWrdurTc1vWDBAq5fv56lnZEjR+Lk5JQfoRdoCxYswNPTM8tDlmWMX86FCxfYtm0bMTExFCpUiBo1alClShXd/kcPqz169CgpKSkUKVIky8Nqw8LCOHnyZJa2H92o/DpLTU1lx44dnDt3TncZXNmyZalfv77uD7GMce497+9YTsY0OTmZDRs2cOHChWwfTAvw1Vdf6ZZNf9zkyZPz54OZmGeNc4MGDfj666+zPe/x30sZ52fLbU721VdfUaNGjSwPWZYxfracjPPx48fZs2cP9+/fx8XFhQYNGuh9CSZ5h2kqEMWXEEIIIYQQQhR08pwvIYQQQgghhDAAKb6EEEIIIYQQwgCk+BJCCCGEEEIIA5DiSwghhBBCCCEMQIovIYQQQgghhDAAKb6EEEIIIYQQwgCk+BJCCCGEEEIIA5DiSwghhBBCCCEMQIovIYQQQgghhDAAKb6EEEIIIYQQwgCk+BJCCCGEEEIIA5DiSwghhBBCCCEM4P/bnyKeL2zm4wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = targets.plot(ax_width=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if needed, it is also possible to get a [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets) object for just some of the sequences specified in `seqsfile` or `feature_parse_specs`.\n",
    "This is most commonly useful when `seqsfile` contains additional sequences that are not of interest.\n",
    "Below we illustrate how to do this by:\n",
    " - Setting `select_target_names` to only keep the *LASV_Josiah_WT* sequence in `seqsfile`\n",
    " - Setting `ingore_feature_parse_specs` to ignore the other targets (in this case, *LASV_Josiah_OPT*) in `feature_parse_specs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:32.446884Z",
     "iopub.status.busy": "2024-05-23T18:06:32.446568Z",
     "iopub.status.idle": "2024-05-23T18:06:32.459985Z",
     "shell.execute_reply": "2024-05-23T18:06:32.459342Z",
     "shell.execute_reply.started": "2024-05-23T18:06:32.446847Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are the names of the retained targets: ['LASV_Josiah_WT']\n"
     ]
    }
   ],
   "source": [
    "targets_subset = alignparse.targets.Targets(\n",
    "    seqsfile=targetfiles,\n",
    "    feature_parse_specs=lasv_parse_specs_file,\n",
    "    allow_extra_features=True,\n",
    "    select_target_names=[\"LASV_Josiah_WT\"],\n",
    "    ignore_feature_parse_specs_keys=[\"LASV_Josiah_OPT\"],\n",
    ")\n",
    "\n",
    "print(f\"Here are the names of the retained targets: {targets_subset.target_names}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PacBio CCSs\n",
    "We will align PacBio circular consensus sequences (CCSs) to the target.\n",
    "First, we want to look at the CCSs.\n",
    "A FASTQ file with these CCSs along with an associated report file were generated using the PacBio `ccs` program (see [here](https://github.com/PacificBiosciences/ccs) for details on `ccs`) using commands like the following (generates report file and BAM of CCSs):\n",
    "\n",
    "    ccs --minLength 50 --maxLength 5000 \\\n",
    "        --minPasses 3  --minPredictedAccuracy 0.999 \\\n",
    "        --reportFile lasv_pilot_report.txt \\\n",
    "        --polish --numThreads 16 \\\n",
    "        lasv_pilot_subreads.bam lasv_pilot_ccs.bam\n",
    "        \n",
    "The BAM file was then converted to a FASTQ file using [samtools](http://www.htslib.org/) with flags to retain the number of passes (`np`) and read quality (`rq`). Here we convert the BAM file to a gzipped FASTQ file to demonstrate the use of compressed files:\n",
    "\n",
    "    samtools bam2fq -T np,rq lasv_pilot_ccs.bam | gzip > lasv_pilot_ccs.fastq.gz\n",
    "    \n",
    "Here is a data frame with the resulting FASTQ and BAM files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:32.461260Z",
     "iopub.status.busy": "2024-05-23T18:06:32.460833Z",
     "iopub.status.idle": "2024-05-23T18:06:32.474256Z",
     "shell.execute_reply": "2024-05-23T18:06:32.473547Z",
     "shell.execute_reply.started": "2024-05-23T18:06:32.461235Z"
    }
   },
   "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>name</th>\n",
       "      <th>report</th>\n",
       "      <th>fastq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>input_files/lasv_pilot_report.txt</td>\n",
       "      <td>input_files/lasv_example_ccs.fastq.gz</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         name                             report  \\\n",
       "0  lasv_pilot  input_files/lasv_pilot_report.txt   \n",
       "\n",
       "                                   fastq  \n",
       "0  input_files/lasv_example_ccs.fastq.gz  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_names = [\"lasv_pilot\"]\n",
    "ccs_dir = \"input_files\"\n",
    "file_name = \"lasv_example\"\n",
    "\n",
    "pacbio_runs = pd.DataFrame(\n",
    "    {\n",
    "        \"name\": run_names,\n",
    "        \"report\": [f\"{ccs_dir}/{name}_report.txt\" for name in run_names],\n",
    "        \"fastq\": [f\"{ccs_dir}/{file_name}_ccs.fastq.gz\"],\n",
    "    }\n",
    ")\n",
    "\n",
    "pacbio_runs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a [Summaries](https://jbloomlab.github.io/alignparse/alignparse.ccs.html#alignparse.ccs.Summaries) object for these CCSs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:32.475522Z",
     "iopub.status.busy": "2024-05-23T18:06:32.475184Z",
     "iopub.status.idle": "2024-05-23T18:06:33.310276Z",
     "shell.execute_reply": "2024-05-23T18:06:33.309033Z",
     "shell.execute_reply.started": "2024-05-23T18:06:32.475494Z"
    }
   },
   "outputs": [],
   "source": [
    "ccs_summaries = alignparse.ccs.Summaries(pacbio_runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot how many ZMWs yielded CCSs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:33.312544Z",
     "iopub.status.busy": "2024-05-23T18:06:33.312072Z",
     "iopub.status.idle": "2024-05-23T18:06:34.138029Z",
     "shell.execute_reply": "2024-05-23T18:06:34.137019Z",
     "shell.execute_reply.started": "2024-05-23T18:06:33.312500Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "p = ccs_summaries.plot_zmw_stats()\n",
    "_ = p.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.139792Z",
     "iopub.status.busy": "2024-05-23T18:06:34.139349Z",
     "iopub.status.idle": "2024-05-23T18:06:34.194531Z",
     "shell.execute_reply": "2024-05-23T18:06:34.193912Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.139764Z"
    }
   },
   "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>name</th>\n",
       "      <th>status</th>\n",
       "      <th>number</th>\n",
       "      <th>fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>Success -- CCS generated</td>\n",
       "      <td>250</td>\n",
       "      <td>0.4996</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>Failed -- Not enough full passes</td>\n",
       "      <td>157</td>\n",
       "      <td>0.3144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>Failed -- CCS below minimum predicted accuracy</td>\n",
       "      <td>86</td>\n",
       "      <td>0.1720</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>Failed -- No usable subreads</td>\n",
       "      <td>6</td>\n",
       "      <td>0.0137</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>Failed -- Other reason</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         name                                          status  number  \\\n",
       "0  lasv_pilot                        Success -- CCS generated     250   \n",
       "1  lasv_pilot                Failed -- Not enough full passes     157   \n",
       "2  lasv_pilot  Failed -- CCS below minimum predicted accuracy      86   \n",
       "3  lasv_pilot                    Failed -- No usable subreads       6   \n",
       "4  lasv_pilot                          Failed -- Other reason       0   \n",
       "\n",
       "   fraction  \n",
       "0    0.4996  \n",
       "1    0.3144  \n",
       "2    0.1720  \n",
       "3    0.0137  \n",
       "4    0.0000  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ccs_summaries.zmw_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistics on the CCSs (length, number of subread passes, accuracy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.195919Z",
     "iopub.status.busy": "2024-05-23T18:06:34.195449Z",
     "iopub.status.idle": "2024-05-23T18:06:34.797352Z",
     "shell.execute_reply": "2024-05-23T18:06:34.796682Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.195892Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "for stat in [\"length\", \"passes\", \"accuracy\"]:\n",
    "    if ccs_summaries.has_stat(stat):\n",
    "        p = ccs_summaries.plot_ccs_stats(stat)\n",
    "        _ = p.draw()\n",
    "    else:\n",
    "        print(f\"No information available on CCS {stat}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Align CCSs to target\n",
    "Now we use [minimap2](https://github.com/lh3/minimap2) to align the CCSs to the target.\n",
    "\n",
    "First, we create a [Mapper](https://jbloomlab.github.io/alignparse/alignparse.minimap2.html#alignparse.minimap2.Mapper) object to run [minimap2](https://github.com/lh3/minimap2), using the options for codon-level deep mutational scanning (specified by [alignparse.minimap2.OPTIONS_CODON_DMS](https://jbloomlab.github.io/alignparse/alignparse.minimap2.html#alignparse.minimap2.OPTIONS_CODON_DMS)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.798626Z",
     "iopub.status.busy": "2024-05-23T18:06:34.798266Z",
     "iopub.status.idle": "2024-05-23T18:06:34.809961Z",
     "shell.execute_reply": "2024-05-23T18:06:34.809078Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.798601Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using `minimap2` 2.22-r1101 with these options:\n",
      "-A2 -B4 -O12 -E2 --end-bonus=13 --secondary=no --cs\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "mapper = alignparse.minimap2.Mapper(alignparse.minimap2.OPTIONS_CODON_DMS)\n",
    "\n",
    "print(\n",
    "    f\"Using `minimap2` {mapper.version} with these options:\\n\"\n",
    "    + \" \".join(mapper.options)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use this mapper to do the alignments to a SAM file.\n",
    "First, add the names of the desired alignment files to our data frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.811699Z",
     "iopub.status.busy": "2024-05-23T18:06:34.811270Z",
     "iopub.status.idle": "2024-05-23T18:06:34.823979Z",
     "shell.execute_reply": "2024-05-23T18:06:34.823306Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.811661Z"
    },
    "scrolled": true
   },
   "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>name</th>\n",
       "      <th>report</th>\n",
       "      <th>fastq</th>\n",
       "      <th>alignments</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>lasv_pilot</td>\n",
       "      <td>input_files/lasv_pilot_report.txt</td>\n",
       "      <td>input_files/lasv_example_ccs.fastq.gz</td>\n",
       "      <td>./output_files/lasv_pilot_alignments.sam</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         name                             report  \\\n",
       "0  lasv_pilot  input_files/lasv_pilot_report.txt   \n",
       "\n",
       "                                   fastq  \\\n",
       "0  input_files/lasv_example_ccs.fastq.gz   \n",
       "\n",
       "                                 alignments  \n",
       "0  ./output_files/lasv_pilot_alignments.sam  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pacbio_runs = pacbio_runs.assign(\n",
    "    alignments=lambda x: outdir + x[\"name\"] + \"_alignments.sam\"\n",
    ")\n",
    "\n",
    "pacbio_runs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run [targets.align](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.align) using the mapper to actually align the FASTQ queries to the target:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.825383Z",
     "iopub.status.busy": "2024-05-23T18:06:34.825027Z",
     "iopub.status.idle": "2024-05-23T18:06:34.930651Z",
     "shell.execute_reply": "2024-05-23T18:06:34.929339Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.825354Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aligning input_files/lasv_example_ccs.fastq.gz to create ./output_files/lasv_pilot_alignments.sam...\n"
     ]
    }
   ],
   "source": [
    "for tup in pacbio_runs.itertuples(index=False):\n",
    "    print(f\"Aligning {tup.fastq} to create {tup.alignments}...\")\n",
    "    targets.align(queryfile=tup.fastq, alignmentfile=tup.alignments, mapper=mapper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These SAM files now contain the alignments along with the [cs tag](https://github.com/lh3/minimap2#cs). \n",
    "\n",
    "An example [cs tag](https://github.com/lh3/minimap2#cs) is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.932964Z",
     "iopub.status.busy": "2024-05-23T18:06:34.932389Z",
     "iopub.status.idle": "2024-05-23T18:06:34.940875Z",
     "shell.execute_reply": "2024-05-23T18:06:34.939933Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.932907Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First alignment in ./output_files/lasv_pilot_alignments.sam has `cs` tag:\n",
      ":8*ng*na*ng*na*nc*ng:19*nc:1667*na:28\n"
     ]
    }
   ],
   "source": [
    "for fname in pacbio_runs[\"alignments\"][:1]:\n",
    "    with pysam.AlignmentFile(fname) as f:\n",
    "        a = next(f)\n",
    "        print(f\"First alignment in {fname} has `cs` tag:\\n\" + a.get_tag(\"cs\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parse the alignments\n",
    "Now we use [Targets.parse_alignments](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.parse_alignment) to parse the SAM files to get the information we specified for return.\n",
    "This function returns a data frame (`readstats`) on the overall parsing stats, plus dicts keyed by the names of each target in [Targets](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets) giving data frames of the aligned and filtered reads. Here we return the `aligned` and `filtered` reads as dictionaries of data frames. \n",
    "\n",
    "The `aligned` and `filtered` read information can instead be returned as CSV files containing these data frames by setting the `to_csv` argument to `True`. Note: When using the combined [Targets.align_and_parse](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.align_and_parse) function as in the other example notebooks ([RecA DMS libraries](https://jbloomlab.github.io/alignparse/recA_DMS.html) and [Single-cell virus sequencing](https://jbloomlab.github.io/alignparse/flu_virus_seq_example.html)), these CSV files will always be created, even if `to_csv` is `False` and the [align_and_parse](https://jbloomlab.github.io/alignparse/alignparse.targets.html#alignparse.targets.Targets.align_and_parse) function is returning `aligned` and `filtered` as dictionaries of data frames in its final output.\n",
    "\n",
    "Here we only have one PacBio run, but in practice we will often have multiple. As such, our example shows how to concatenate the `readstats`, `aligned`, and `filtered` data frames for each PacBio run and then look at the data frames together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:34.942818Z",
     "iopub.status.busy": "2024-05-23T18:06:34.942189Z",
     "iopub.status.idle": "2024-05-23T18:06:35.502285Z",
     "shell.execute_reply": "2024-05-23T18:06:35.501601Z",
     "shell.execute_reply.started": "2024-05-23T18:06:34.942783Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parsing PacBio run lasv_pilot\n"
     ]
    }
   ],
   "source": [
    "readstats = []\n",
    "aligned = {targetname: [] for targetname in targets.target_names}\n",
    "filtered = {targetname: [] for targetname in targets.target_names}\n",
    "\n",
    "for run in pacbio_runs.itertuples():\n",
    "    print(f\"Parsing PacBio run {run.name}\")\n",
    "    run_readstats, run_aligned, run_filtered = targets.parse_alignment(\n",
    "        run.alignments, filtered_cs=True\n",
    "    )\n",
    "\n",
    "    # when concatenating add the run name to keep track of runs for results\n",
    "    readstats.append(run_readstats.assign(run_name=run.name))\n",
    "    for targetname in targets.target_names:\n",
    "        aligned[targetname].append(run_aligned[targetname].assign(run_name=run.name))\n",
    "        filtered[targetname].append(run_filtered[targetname].assign(run_name=run.name))\n",
    "\n",
    "# now concatenate the data frames for each run\n",
    "readstats = pd.concat(readstats, ignore_index=True, sort=False)\n",
    "for targetname in targets.target_names:\n",
    "    aligned[targetname] = pd.concat(aligned[targetname], ignore_index=True, sort=False)\n",
    "    filtered[targetname] = pd.concat(\n",
    "        filtered[targetname], ignore_index=True, sort=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First lets look at the read stats:\n",
    "\n",
    "From the known composition of the library, there should be more `LASV_Josiah_OPT` reads than `LASV_Josiah_WT`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:35.503468Z",
     "iopub.status.busy": "2024-05-23T18:06:35.503245Z",
     "iopub.status.idle": "2024-05-23T18:06:35.512662Z",
     "shell.execute_reply": "2024-05-23T18:06:35.511886Z",
     "shell.execute_reply.started": "2024-05-23T18:06:35.503446Z"
    },
    "scrolled": true
   },
   "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>category</th>\n",
       "      <th>count</th>\n",
       "      <th>run_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>filtered LASV_Josiah_WT</td>\n",
       "      <td>23</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>aligned LASV_Josiah_WT</td>\n",
       "      <td>84</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>filtered LASV_Josiah_OPT</td>\n",
       "      <td>32</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>aligned LASV_Josiah_OPT</td>\n",
       "      <td>111</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>unmapped</td>\n",
       "      <td>0</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   category  count    run_name\n",
       "0   filtered LASV_Josiah_WT     23  lasv_pilot\n",
       "1    aligned LASV_Josiah_WT     84  lasv_pilot\n",
       "2  filtered LASV_Josiah_OPT     32  lasv_pilot\n",
       "3   aligned LASV_Josiah_OPT    111  lasv_pilot\n",
       "4                  unmapped      0  lasv_pilot"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "readstats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:35.518503Z",
     "iopub.status.busy": "2024-05-23T18:06:35.518113Z",
     "iopub.status.idle": "2024-05-23T18:06:35.649601Z",
     "shell.execute_reply": "2024-05-23T18:06:35.648283Z",
     "shell.execute_reply.started": "2024-05-23T18:06:35.518472Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "p = (\n",
    "    ggplot(readstats, aes(\"category\", \"count\"))\n",
    "    + geom_bar(stat=\"identity\")\n",
    "    + facet_wrap(\"~ run_name\", nrow=1)\n",
    "    + theme(\n",
    "        axis_text_x=element_text(angle=90), figure_size=(1.5 * len(pacbio_runs), 2.5)\n",
    "    )\n",
    ")\n",
    "_ = p.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look at the information on the filtered reads.\n",
    "This is a bigger data frame, so we just look at the first few lines for the first target (of which there is only one anyway):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:35.651362Z",
     "iopub.status.busy": "2024-05-23T18:06:35.651000Z",
     "iopub.status.idle": "2024-05-23T18:06:35.663816Z",
     "shell.execute_reply": "2024-05-23T18:06:35.662995Z",
     "shell.execute_reply.started": "2024-05-23T18:06:35.651338Z"
    }
   },
   "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>query_name</th>\n",
       "      <th>filter_reason</th>\n",
       "      <th>filter_cs</th>\n",
       "      <th>run_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194989/ccs</td>\n",
       "      <td>termini5 clip5</td>\n",
       "      <td>*nc*na*nc:19*nc:113</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4260241/ccs</td>\n",
       "      <td>index clip5</td>\n",
       "      <td>*ng*na*nc*na*nc</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4260334/ccs</td>\n",
       "      <td>query_clip5</td>\n",
       "      <td>None</td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>m54228_190605_190010/4391712/ccs</td>\n",
       "      <td>termini5 clip5</td>\n",
       "      <td></td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>m54228_190605_190010/4456843/ccs</td>\n",
       "      <td>termini3 clip3</td>\n",
       "      <td></td>\n",
       "      <td>lasv_pilot</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name   filter_reason            filter_cs  \\\n",
       "0  m54228_190605_190010/4194989/ccs  termini5 clip5  *nc*na*nc:19*nc:113   \n",
       "1  m54228_190605_190010/4260241/ccs     index clip5      *ng*na*nc*na*nc   \n",
       "2  m54228_190605_190010/4260334/ccs     query_clip5                 None   \n",
       "3  m54228_190605_190010/4391712/ccs  termini5 clip5                        \n",
       "4  m54228_190605_190010/4456843/ccs  termini3 clip3                        \n",
       "\n",
       "     run_name  \n",
       "0  lasv_pilot  \n",
       "1  lasv_pilot  \n",
       "2  lasv_pilot  \n",
       "3  lasv_pilot  \n",
       "4  lasv_pilot  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filtered[targets.target_names[0]].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:35.665706Z",
     "iopub.status.busy": "2024-05-23T18:06:35.665307Z",
     "iopub.status.idle": "2024-05-23T18:06:36.006445Z",
     "shell.execute_reply": "2024-05-23T18:06:36.005674Z",
     "shell.execute_reply.started": "2024-05-23T18:06:35.665669Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "for targetname in targets.target_names:\n",
    "    target_filtered = filtered[targetname]\n",
    "    nreasons = target_filtered[\"filter_reason\"].nunique()\n",
    "    p = (\n",
    "        ggplot(target_filtered, aes(\"filter_reason\"))\n",
    "        + geom_bar()\n",
    "        + facet_wrap(\"~ run_name\", nrow=1)\n",
    "        + labs(title=targetname)\n",
    "        + theme(\n",
    "            axis_text_x=element_text(angle=90),\n",
    "            figure_size=(0.3 * nreasons * len(pacbio_runs), 2.5),\n",
    "        )\n",
    "    )\n",
    "    _ = p.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error filtering\n",
    "\n",
    "Before looking at the information for the validly aligned (not filtered) reads, it is important to get a sense of the error rate for these sequencing reads. \n",
    "\n",
    "These reads do not have random barcodes on the initial viral entry protein plasmids, but we can use the `gene_accuracy` information output from constructing the `ccs`s to examine accuracy. \n",
    "\n",
    "We will do this using a similar method to that implemented in the [RecA DMS libraries](https://jbloomlab.github.io/alignparse/recA_DMS.html) example notebook. However, here we will plot the graphs for each target. \n",
    "\n",
    "We anticipate excluding all CCSs for which the error rate for either the gene or barcode is $>10^{-4}$.\n",
    "We specify this cutoff below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.007729Z",
     "iopub.status.busy": "2024-05-23T18:06:36.007417Z",
     "iopub.status.idle": "2024-05-23T18:06:36.675127Z",
     "shell.execute_reply": "2024-05-23T18:06:36.674361Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.007703Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "error_rate_floor = 1e-7  # error rates < this set to this\n",
    "error_cutoff = 1e-4\n",
    "\n",
    "for targetname in targets.target_names:\n",
    "    aligned[targetname] = aligned[targetname].assign(\n",
    "        gene_error=lambda x: numpy.clip(1 - x[\"gene_accuracy\"], error_rate_floor, None)\n",
    "    )\n",
    "    p = (\n",
    "        ggplot(\n",
    "            aligned[targetname].melt(\n",
    "                id_vars=[\"run_name\"],\n",
    "                value_vars=[\"gene_error\"],\n",
    "                var_name=\"feature_type\",\n",
    "                value_name=\"error rate\",\n",
    "            ),\n",
    "            aes(\"error rate\"),\n",
    "        )\n",
    "        + geom_histogram(bins=25)\n",
    "        + geom_vline(xintercept=error_cutoff, linetype=\"dashed\", color=CBPALETTE[1])\n",
    "        + facet_grid(\"~ feature_type\")\n",
    "        + theme(figure_size=(3, 3))\n",
    "        + labs(y=(\"number of CCSs\"), title=(targetname))\n",
    "        + scale_x_log10()\n",
    "    )\n",
    "\n",
    "    _ = p.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we store all reads with an error rate $<10^{-4}$ in new data frames for retained sequences. \n",
    "\n",
    "We will use these retained sequences for further analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.676680Z",
     "iopub.status.busy": "2024-05-23T18:06:36.676344Z",
     "iopub.status.idle": "2024-05-23T18:06:36.683034Z",
     "shell.execute_reply": "2024-05-23T18:06:36.682323Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.676656Z"
    }
   },
   "outputs": [],
   "source": [
    "retained = {targetname: [] for targetname in targets.target_names}\n",
    "for targetname in targets.target_names:\n",
    "    target_retained = aligned[targetname][\n",
    "        aligned[targetname][\"gene_error\"] <= error_cutoff\n",
    "    ].reset_index(drop=True)\n",
    "    retained[targetname] = target_retained"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data analysis\n",
    "\n",
    "Now we can examine our data in more detail. \n",
    "\n",
    "We will only display a few columns so the example renders well in the documentation. First, let's look at the `query_name`, `gene_mutations`, and `index_sequence` columns for the first few entries in the `retained` data frame for each target:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.684798Z",
     "iopub.status.busy": "2024-05-23T18:06:36.684431Z",
     "iopub.status.idle": "2024-05-23T18:06:36.702566Z",
     "shell.execute_reply": "2024-05-23T18:06:36.701654Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.684771Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_WT\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194436/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194472/ccs</td>\n",
       "      <td>T572C</td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194509/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>m54228_190605_190010/4194553/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>m54228_190605_190010/4194576/ccs</td>\n",
       "      <td></td>\n",
       "      <td>AGACAC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194436/ccs                        GGTATG\n",
       "1  m54228_190605_190010/4194472/ccs          T572C         GGTATG\n",
       "2  m54228_190605_190010/4194509/ccs                        GGTATG\n",
       "3  m54228_190605_190010/4194553/ccs                        GGTATG\n",
       "4  m54228_190605_190010/4194576/ccs                        AGACAC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_OPT\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194382/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GAGACG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194390/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194399/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>m54228_190605_190010/4194439/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>m54228_190605_190010/4194445/ccs</td>\n",
       "      <td></td>\n",
       "      <td>CTTCAC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194382/ccs                        GAGACG\n",
       "1  m54228_190605_190010/4194390/ccs                        ACGACC\n",
       "2  m54228_190605_190010/4194399/ccs                        ACGACC\n",
       "3  m54228_190605_190010/4194439/ccs                        ACGACC\n",
       "4  m54228_190605_190010/4194445/ccs                        CTTCAC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_columns = [\"query_name\", \"gene_mutations\", \"index_sequence\"]\n",
    "for target_name in targets.target_names:\n",
    "    print(target_name)\n",
    "    display(retained[target_name][display_columns].head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen by the `index_sequence` columns, each target here has 2 or 3 different indices that map to it. These indicies indicate different samples. \n",
    "\n",
    "We can then split these data frames into sample-specific data frames based on the known starting index sequences for each target. We will store these data frames in a new dictionary, keyed by target and index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.704097Z",
     "iopub.status.busy": "2024-05-23T18:06:36.703709Z",
     "iopub.status.idle": "2024-05-23T18:06:36.717796Z",
     "shell.execute_reply": "2024-05-23T18:06:36.716977Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.704065Z"
    }
   },
   "outputs": [],
   "source": [
    "indices = {\n",
    "    \"LASV_Josiah_WT\": [\"AGACAC\", \"GGTATG\"],\n",
    "    \"LASV_Josiah_OPT\": [\"ACGACC\", \"CTTCAC\", \"GAGACG\"],\n",
    "}\n",
    "target_idx_retained = {}\n",
    "index_counts = {target: [] for target in indices}\n",
    "index_counts_dfs = {}\n",
    "for target_name in targets.target_names:\n",
    "    for index in indices[target_name]:\n",
    "        target_idx_retained[f\"{target_name}_{index}\"] = retained[target_name][\n",
    "            retained[target_name][\"index_sequence\"] == index\n",
    "        ].reset_index(drop=True)\n",
    "        index_counts[f\"{target_name}\"].append(\n",
    "            (index, len(target_idx_retained[f\"{target_name}_{index}\"]))\n",
    "        )\n",
    "    index_counts[f\"{target_name}\"].append(\n",
    "        (\n",
    "            \"invalid\",\n",
    "            (\n",
    "                len(retained[target_name])\n",
    "                - sum(idx_tup[1] for idx_tup in index_counts[target_name])\n",
    "            ),\n",
    "        )\n",
    "    )\n",
    "    index_counts_dfs[target_name] = pd.DataFrame(\n",
    "        index_counts[target_name], columns=[\"index\", \"count\"]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the process of making these separate data frames, we also kept track of how many reads for each target map to each index or don't map to an index (so have an 'invalid' index) and can now plot these counts. As seen below, only one sequence in this example has an invalid index and it maps to the `LASV_Josiah_OPT` target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.719460Z",
     "iopub.status.busy": "2024-05-23T18:06:36.719131Z",
     "iopub.status.idle": "2024-05-23T18:06:36.988751Z",
     "shell.execute_reply": "2024-05-23T18:06:36.987969Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.719436Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "for target_name in targets.target_names:\n",
    "    df = index_counts_dfs[target_name]\n",
    "    df[\"target\"] = [target_name] * len(df)\n",
    "    id_order = [\"invalid\"] + df[\"index\"][:-1].to_list()\n",
    "    df[\"index\"] = pd.Categorical(df[\"index\"], categories=id_order, ordered=True)\n",
    "    index_count_plot = (\n",
    "        ggplot(df, aes(x=\"target\", y=\"count\", fill=\"index\"))\n",
    "        + geom_bar(stat=\"identity\", position=\"stack\")\n",
    "        + scale_fill_manual(values=CBPALETTE)\n",
    "        + theme(\n",
    "            axis_text_x=element_text(angle=90, vjust=1, hjust=0.5), figure_size=(1, 2)\n",
    "        )\n",
    "        + ylab(\"Reads\")\n",
    "        + xlab(\"Target\")\n",
    "        + ggtitle(\"Reads per Sample Index\")\n",
    "    )\n",
    "\n",
    "    _ = index_count_plot.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two target-specific data frames are now five index-specific data frames. Again we will display the `query_name`, `gene_mutations`, and `index_sequence` columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:36.990634Z",
     "iopub.status.busy": "2024-05-23T18:06:36.990275Z",
     "iopub.status.idle": "2024-05-23T18:06:37.029580Z",
     "shell.execute_reply": "2024-05-23T18:06:37.028748Z",
     "shell.execute_reply.started": "2024-05-23T18:06:36.990612Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_WT_AGACAC\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194576/ccs</td>\n",
       "      <td></td>\n",
       "      <td>AGACAC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194612/ccs</td>\n",
       "      <td></td>\n",
       "      <td>AGACAC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194613/ccs</td>\n",
       "      <td></td>\n",
       "      <td>AGACAC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194576/ccs                        AGACAC\n",
       "1  m54228_190605_190010/4194612/ccs                        AGACAC\n",
       "2  m54228_190605_190010/4194613/ccs                        AGACAC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_WT_GGTATG\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194436/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194472/ccs</td>\n",
       "      <td>T572C</td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194509/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GGTATG</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194436/ccs                        GGTATG\n",
       "1  m54228_190605_190010/4194472/ccs          T572C         GGTATG\n",
       "2  m54228_190605_190010/4194509/ccs                        GGTATG"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_OPT_ACGACC\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194390/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194399/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194439/ccs</td>\n",
       "      <td></td>\n",
       "      <td>ACGACC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194390/ccs                        ACGACC\n",
       "1  m54228_190605_190010/4194399/ccs                        ACGACC\n",
       "2  m54228_190605_190010/4194439/ccs                        ACGACC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_OPT_CTTCAC\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194445/ccs</td>\n",
       "      <td></td>\n",
       "      <td>CTTCAC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194501/ccs</td>\n",
       "      <td></td>\n",
       "      <td>CTTCAC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194549/ccs</td>\n",
       "      <td></td>\n",
       "      <td>CTTCAC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194445/ccs                        CTTCAC\n",
       "1  m54228_190605_190010/4194501/ccs                        CTTCAC\n",
       "2  m54228_190605_190010/4194549/ccs                        CTTCAC"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LASV_Josiah_OPT_GAGACG\n"
     ]
    },
    {
     "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>query_name</th>\n",
       "      <th>gene_mutations</th>\n",
       "      <th>index_sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>m54228_190605_190010/4194382/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GAGACG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>m54228_190605_190010/4194467/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GAGACG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>m54228_190605_190010/4194491/ccs</td>\n",
       "      <td></td>\n",
       "      <td>GAGACG</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         query_name gene_mutations index_sequence\n",
       "0  m54228_190605_190010/4194382/ccs                        GAGACG\n",
       "1  m54228_190605_190010/4194467/ccs                        GAGACG\n",
       "2  m54228_190605_190010/4194491/ccs                        GAGACG"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for target_idx in target_idx_retained:\n",
    "    print(target_idx)\n",
    "    display(target_idx_retained[target_idx][display_columns].head(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add mutation info for the genes in these targets using [alignparse.consensus.add_mut_info_cols](https://jbloomlab.github.io/alignparse/alignparse.consensus.html#alignparse.consensus.add_mut_info_cols)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:37.030955Z",
     "iopub.status.busy": "2024-05-23T18:06:37.030612Z",
     "iopub.status.idle": "2024-05-23T18:06:37.046260Z",
     "shell.execute_reply": "2024-05-23T18:06:37.045472Z",
     "shell.execute_reply.started": "2024-05-23T18:06:37.030927Z"
    }
   },
   "outputs": [],
   "source": [
    "for target_idx in target_idx_retained:\n",
    "    target_idx_retained[target_idx] = alignparse.consensus.add_mut_info_cols(\n",
    "        target_idx_retained[target_idx],\n",
    "        mutation_col=\"gene_mutations\",\n",
    "        n_sub_col=\"n_gene_subs\",\n",
    "        n_indel_col=\"n_gene_indels\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then plot the number of reads for each sample that has each number of substitutions or indels to get a sense of how many mutations are in these reads, if some smaples have more mutations than others, and if substitutions or indels are more prevalent in these sequences. \n",
    "\n",
    "With such a small sample snippet of the data, it is difficult to make any conclusions, but with a full dataset, this can provide important information about the mutational processes at work for each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-23T18:06:37.047681Z",
     "iopub.status.busy": "2024-05-23T18:06:37.047263Z",
     "iopub.status.idle": "2024-05-23T18:06:38.053149Z",
     "shell.execute_reply": "2024-05-23T18:06:38.051878Z",
     "shell.execute_reply.started": "2024-05-23T18:06:37.047655Z"
    }
   },
   "outputs": [],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "mut_cols = [\"n_gene_subs\", \"n_gene_indels\"]\n",
    "for target_idx in target_idx_retained:\n",
    "    target_idx_plot_muts = target_idx_retained[target_idx][mut_cols].melt()\n",
    "    mut_counts_plot = (\n",
    "        ggplot(target_idx_plot_muts, aes(\"value\"))\n",
    "        + geom_bar()\n",
    "        + facet_wrap(\"~ variable\", ncol=2)\n",
    "        + xlim(-0.5, 2.5)\n",
    "        + labs(title=f\"{target_idx}\")\n",
    "        + theme(\n",
    "            axis_text_x=element_text(angle=90),\n",
    "            figure_size=(3, 2.5),\n",
    "        )\n",
    "    )\n",
    "    _ = mut_counts_plot.draw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}