example_usage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pylazybam example useage\n",
    "\n",
    "Matthew Wakefield\n",
    "\n",
    "<small>The Walter and Eliza Hall Institute and The University of Melbourne</small>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example file is intended to give a flavour of the potential usages of pylazybam. As a library intened for interactive use in notebooks and bespoke scripts its intentionally very open ended and flexable.\n",
    "\n",
    "First you will need to install pylazybam following the instructions in the readme. There are no dependencies, so this should be easy even on Windows and ARM based linux, although you will have issues with anything other than trivial BAM files on a Raspberry Pi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylazybam import bam\n",
    "import gzip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import as `bam` as it gives an easy to type, short, safe and descriptive namespace.\n",
    "You can get help on pylazybam functions by using `help()` or in a notebook prefacing with a `?`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mInit signature:\u001b[0m \u001b[0mbam\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFileReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mubam\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m<\u001b[0m\u001b[0;32mclass\u001b[0m \u001b[0;34m'BinaryIO'\u001b[0m\u001b[0;34m>\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m     \n",
       "A Pure Python Lazy Bam Parser Class\n",
       "\n",
       "Parameters\n",
       "----------\n",
       "    ubam : BinaryIO\n",
       "        An binary (bytes) file or stream containing a valid uncompressed\n",
       "        bam file conforming to the specification.\n",
       "\n",
       "Yields\n",
       "------\n",
       "    align : bytes\n",
       "        A byte string of a bam alignment entry in raw binary format\n",
       "\n",
       "Attributes\n",
       "----------\n",
       "    header : str\n",
       "        The ASCII representation of the header\n",
       "\n",
       "    index_to_ref : Dict[int:str]\n",
       "        A dictionary mapping bam reference numeric identifiers to names\n",
       "\n",
       "    raw_header : bytes\n",
       "        The raw bytestring representing the bam header\n",
       "\n",
       "    raw_refs : bytes\n",
       "        The raw bytestring representing the reference sequences\n",
       "\n",
       "    refs : Dict[str:int]\n",
       "        A dictionary of reference_name keys with reference_length\n",
       "\n",
       "    ref_to_index : Dict[str:int]\n",
       "        A dictionary mapping reference names to the bam numeric identifier\n",
       "\n",
       "    sort_order : str\n",
       "        The value of the SO field indicating sort type. Value is as given in\n",
       "        the BAM file.\n",
       "        Should be one of 'unknown', 'unsorted', 'queryname' or 'coordinate'\n",
       "\n",
       "Notes\n",
       "-----\n",
       "    It is advisable not to call the private functions or operate directly\n",
       "    on the underlying file object.\n",
       "\n",
       "    The detailed specification for the BAM format can be found at\n",
       "    https://samtools.github.io/hts-specs/SAMv1.pdf\n",
       "\n",
       "Example\n",
       "-------\n",
       "    This class requires an uncompressed bam file as input.\n",
       "    For this example we will use gzip to decompress the test file included\n",
       "    as a resource in the package\n",
       "\n",
       "    >>> import gzip\n",
       "    >>> from pkg_resources import resource_stream\n",
       "    >>> ubam = gzip.open('/tests/data/paired_end_testdata_human.bam'))\n",
       "\n",
       "    A bam filereader object can then be created and headers inspected\n",
       "    >>> mybam = bam.FileReader(ubam)\n",
       "    >>> print(mybam.header)\n",
       "\n",
       "    The filereader object is an iterator and yields alignments in raw format\n",
       "    >>> align = next(mybam)\n",
       "\n",
       "    Alignments can be processed using functions from pylazybam.bam\n",
       "    >>> print(mybam.index_to_ref[get_ref_index(align)],\n",
       "    >>>       get_pos(align),\n",
       "    >>>       get_AS(align))\n",
       "\u001b[0;31mInit docstring:\u001b[0m         \n",
       "\u001b[0;31mFile:\u001b[0m           ~/repos/pylazybam/pylazybam/bam.py\n",
       "\u001b[0;31mType:\u001b[0m           type\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "?bam.FileReader"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using data generated by the NYU genome centre from one of the well know genome in a bottle references. This case has been sequenced to high depth, and mapped against a GRCh38 reference with decoys and ALTs\n",
    "It can be downloaded from `lftp ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR323/ERR3239484/NA12778.final.cram`\n",
    "The file we are using has had some manipulations done on it with `samtools view` but the exact manipulations have been lost to the mists of time and advancing years of the operator. Time for some bioinfosleuthing<sup>(tm)</sup> Lets open the file up and have a look at what we can see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@PG\tID:MarkDuplicates\tVN:2.4.1(7c4d36e011df1aec4689b51efcada44e92d1817f_1464389670)\tCL:picard.sam.markduplicates.MarkDuplicates INPUT=[/data/analysis/DarnellR/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/analysis/NA12778.merged.bam] OUTPUT=/scratch/3ad0b0f54a6c428eb90eea99b52683dd/NA12778.bam METRICS_FILE=/data/analysis/DarnellR/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/qc/NA12778.dedup.metrics VALIDATION_STRINGENCY=SILENT COMPRESSION_LEVEL=0 MAX_RECORDS_IN_RAM=2000000    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag REMOVE_DUPLICATES=false ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json\tPN:MarkDuplicates\n",
      "@PG\tID:bwa\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L003\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKV5DSXX.3.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L003_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L003_001.R2.fastq.gz\n",
      "@PG\tID:bwa.1\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L003\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HCKWTDSXX.3.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L003_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L003_001.R2.fastq.gz\n",
      "@PG\tID:bwa.2\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L002\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKV5DSXX.2.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L002_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L002_001.R2.fastq.gz\n",
      "@PG\tID:bwa.3\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HCKWTDSXX.1.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001_001.R2.fastq.gz\n",
      "@PG\tID:bwa.4\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L002\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HCKWTDSXX.2.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L002_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L002_001.R2.fastq.gz\n",
      "@PG\tID:bwa.5\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L003\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKGLDSXX.3.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L003_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L003/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L003_001.R2.fastq.gz\n",
      "@PG\tID:bwa.6\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L004\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKV5DSXX.4.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L004_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L004_001.R2.fastq.gz\n",
      "@PG\tID:bwa.7\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L002\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKGLDSXX.2.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L002_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L002/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L002_001.R2.fastq.gz\n",
      "@PG\tID:bwa.8\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L001\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKGLDSXX.1.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L001_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L001_001.R2.fastq.gz\n",
      "@PG\tID:bwa.9\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L001\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKV5DSXX.1.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L001_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0047_BHFKV5DSXX/01.2018-11-20_11.07.57/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L001/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKV5DSXX_L001_001.R2.fastq.gz\n",
      "@PG\tID:bwa.A\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L004\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HFKGLDSXX.4.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L004_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00404/181118_SNA00404_0046_AHFKGLDSXX/01.2018-11-20_11.05.20/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HFKGLDSXX_L004_001.R2.fastq.gz\n",
      "@PG\tID:bwa.B\tPN:bwa\tVN:0.7.15-r1140\tCL:/data/NYGC/Software/bwa/bwa-0.7.15/bwa mem -Y -K 100000000 -t 16 -R @RG\\tID:NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L004\\tPL:illumina\\tPM:Unknown\\tLB:NA12778\\tDS:GRCh38\\tSM:NA12778\\tCN:NYGenome\\tPU:HCKWTDSXX.4.CTCACCAA /gpfs/internal/sweng/production/Resources/GRCh38_1000genomes/GRCh38_full_analysis_set_plus_decoy_hla.fa /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L004_001.R1.fastq.gz /gpfs/internal/sequence/demux/ILLUMINA_NOVASEQ/SNA00132/181109_SNA00132_0055_BHCKWTDSXX/01.2018-11-11_11.00.38/CCDG_13607_B01_GRM_WGS/NYGC/Project_CCDG_13607_B01_GRM_WGS/Sample_NA12778/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L004/fastq/NA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L004_001.R2.fastq.gz\n",
      "@PG\tID:GATK PrintReads\tVN:3.5-0-g36282e4\tCL:readGroup=null platform=null number=-1 sample_file=[] sample_name=[] simplify=false no_pg_tag=false\n"
     ]
    }
   ],
   "source": [
    "NA12778 = bam.FileReader(gzip.open('Genomes/NA12778_GRCh38_HLA_namesorted.bam'))\n",
    "for headerline in NA12778.header.split('\\n'):\n",
    "    if headerline[:3] == '@PG':\n",
    "        print(headerline)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok so we can see a lot of programs have operated on this file, but there are no entries for samtools so we dont know the basis for the subsetting. Lets see if we can work it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counter({5: 1206714, 3343: 1855, 3326: 391, 3330: 327, 3309: 289, 3307: 284, 3328: 219, 3308: 204, 3329: 179, 3327: 179, 3170: 41, 3164: 38, 3165: 38, 3171: 37, 3167: 36, 3035: 35, 3034: 32, 3169: 31, 3038: 29, 3172: 26, 2877: 24, 3177: 24, 3190: 24, 2861: 24, 3041: 23, 2855: 23, 3166: 23, 2867: 22, 3033: 22, 3188: 21, 3189: 21, 3192: 21, 2863: 21, 2873: 20, 2853: 20, 2875: 20, 3184: 20, 3347: 20, 3168: 20, 2870: 20, 3173: 19, 2866: 19, 2852: 19, 2857: 19, 2871: 19, 3037: 19, 3174: 19, 3183: 18, 2884: 18, 3181: 18, 3070: 18, 3071: 18, 2879: 18, 2856: 18, 3306: 18, 2862: 18, 2885: 18, 3179: 17, 3074: 17, 3191: 17, 2859: 17, 3077: 17, 2878: 17, 2864: 17, 3178: 17, 2869: 16, 2860: 16, 3075: 16, 3079: 16, 3180: 15, 2970: 15, 3080: 15, 2872: 15, 2882: 15, 2956: 14, 3186: 14, 2881: 14, 2876: 14, 3039: 13, 2854: 13, 2865: 13, 2868: 13, 3176: 13, 2883: 13, 2874: 13, 3076: 12, 2967: 12, 3185: 12, 3072: 12, 3107: 12, 3353: 11, 2858: 11, 2962: 11, 3187: 11, 3102: 10, 3082: 10, 2968: 10, 3106: 10, 3040: 10, 3280: 9, 2959: 9, 2961: 9, 3284: 9, 3114: 9, 2960: 9, 2965: 8, 3036: 8, 2958: 8, 3175: 8, 3182: 8, 3355: 8, 3283: 8, 3344: 7, 3088: 7, 3081: 7, 3103: 7, 3354: 7, 3151: 7, 2957: 7, 3279: 7, 2969: 6, 3073: 6, 3113: 6, 3274: 6, 3030: 6, 3078: 5, 2966: 5, 3345: 5, 3155: 5, 3200: 5, 3112: 5, 2955: 5, 3094: 5, 3093: 5, 3286: 5, 3278: 5, 3154: 5, 3351: 5, 3156: 5, 3266: 5, 3357: 5, 2963: 5, 3362: 5, 2880: 5, 3100: 5, 3282: 4, 3207: 4, 2949: 4, 3089: 4, 3216: 4, 3097: 4, 3153: 4, 3157: 4, 3109: 4, 3281: 4, 3031: 4, 3086: 4, 3099: 4, 3268: 4, 2964: 4, 2945: 4, 3285: 4, 3087: 3, 3158: 3, 3161: 3, 3163: 3, 3211: 3, 3123: 3, 3096: 3, 3053: 3, 3110: 3, 3162: 3, 2954: 3, 2992: 3, 3217: 3, 2920: 3, 3108: 3, 3253: 3, 3356: 3, 2948: 3, 3159: 3, 3101: 3, 3026: 3, 3160: 2, 2994: 2, 3130: 2, 3203: 2, 3027: 2, 3202: 2, 3049: 2, 3032: 2, 2930: 2, 3150: 2, 2908: 2, 3196: 2, 3359: 2, 3133: 2, 3029: 2, 2944: 2, 3011: 2, 3270: 2, 3263: 2, 2971: 2, 3132: 2, 2950: 2, 3275: 2, 2926: 2, 2985: 2, 3199: 2, 3142: 2, 3262: 2, 3276: 2, 3264: 2, 3259: 2, 2998: 2, 3052: 2, 3252: 2, 2993: 2, 2918: 2, 3210: 2, 3195: 2, 3204: 2, 3267: 2, 3317: 1, 2995: 1, 2910: 1, 3205: 1, 3140: 1, 2927: 1, 2913: 1, 3117: 1, 3043: 1, 3212: 1, 3213: 1, 3066: 1, 2929: 1, 2986: 1, 2897: 1, 2911: 1, 3004: 1, 3119: 1, 3024: 1, 2915: 1, 2923: 1, 2990: 1, 2947: 1, 3084: 1, 3023: 1, 3289: 1, 2952: 1, 3008: 1, 3193: 1, 2973: 1, 3272: 1, 2849: 1, 3214: 1, 2940: 1, 3256: 1, 3271: 1, 3250: 1, 2953: 1, 3273: 1, 3098: 1, 3005: 1, 3311: 1, 3126: 1, 3006: 1, 3143: 1, 3277: 1, 3002: 1, 3209: 1, 3295: 1, 2933: 1, 3365: 1, 3236: 1, 2894: 1, 3091: 1, 2844: 1, 3141: 1, 3020: 1, 3255: 1, 3292: 1, 3206: 1, 3310: 1, 3363: 1, 3201: 1, 3152: 1, 2914: 1, 3128: 1, 2939: 1, 2899: 1, 3021: 1, 3339: 1, 3358: 1, 3352: 1, 2997: 1, 2912: 1, 2925: 1, 3251: 1, 2928: 1, 3092: 1, 2989: 1, 2987: 1, 2909: 1, 3095: 1, 3194: 1, 3349: 1, 3047: 1, 2943: 1, 3090: 1})\n"
     ]
    }
   ],
   "source": [
    "from collections import Counter\n",
    "id_counts = Counter()\n",
    "for align in NA12778:\n",
    "    id_counts[bam.get_ref_index(align)] += 1\n",
    "print(id_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have over 1.2M reads from index 5, so that must be the chromosome that was specified in the samtools view region. We can look up its name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'chr6'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NA12778.index_to_ref[5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chromosome 6 makes sense at it is MHC related. Lets see if we can work out what the start and end of the selected regions were."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reads from chr6 start at 28509971 and end at 33480576\n"
     ]
    }
   ],
   "source": [
    "with bam.FileReader(gzip.open('Genomes/NA12778_GRCh38_HLA_namesorted.bam')) as NA12778:\n",
    "    pos = []\n",
    "    for align in NA12778:\n",
    "        if bam.get_ref_index(align) == 5:\n",
    "            pos.append(bam.get_pos(align))\n",
    "    print(f'Reads from {NA12778.index_to_ref[5]} start at {min(pos)} and end at {max(pos)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that reads overlapping the region and their pairs will be included in the output, it looks like samtools view - chr:28510000-33480000 is the most likely command used to generate the file, likely converting from CRAM to BAM in the process. Now lets check if it has been sorted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "queryname\n"
     ]
    }
   ],
   "source": [
    "NA12778 = bam.FileReader(gzip.open('Genomes/NA12778_GRCh38_HLA_namesorted.bam'))\n",
    "print(NA12778.sort_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So just like it suggest in the filename this BAM has reads that map to the HLA region on chromosome 6 in a namesorted order. This is just what we need for the next step which is selecting those reads that have mappings on the HLA alt chromosomes, and writing them to a new file.\n",
    "\n",
    "To do this we will set up a bam.FileReader and a bam.FileWriter. We will copy across the header information for the new file from the original, and we will add an @PG line to indicate what we have done. It's a bit hard to make it completely reproducable with an ad hoc script, but we can show we did something and indicate our intent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "NA12778 = bam.FileReader(gzip.open('Genomes/NA12778_GRCh38_HLA_namesorted.bam'))\n",
    "NA12778_any_HLA_hits = bam.FileWriter('Genomes/NA12778_any_HLA_hits.bam')\n",
    "NA12778_any_HLA_hits.raw_header = NA12778.raw_header\n",
    "NA12778_any_HLA_hits.raw_refs = NA12778.raw_refs\n",
    "NA12778_any_HLA_hits.update_header(id = 'pylazybam_any_HLA',\n",
    "                                   program = 'pylazybam',\n",
    "                                   version = '0.1.0',\n",
    "                                   description = \"all alignments for reads where any alignment maps to a HLA alt contig\")\n",
    "NA12778_any_HLA_hits.write_header()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a file to read from, a file to write to, and all we need now are some functions to manipulate the reads along the way. The first function is to get all the alignments that belong to a given read pair at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_alignment_batches(in_bam):\n",
    "    previous_name = None\n",
    "    alignbatch = []\n",
    "    for align in in_bam:\n",
    "        name = bam.get_raw_read_name(align, bam.get_len_read_name(align))\n",
    "        if previous_name and name != previous_name:\n",
    "            yield alignbatch\n",
    "            alignbatch = []\n",
    "        alignbatch.append(align)\n",
    "        previous_name = name\n",
    "    if alignbatch:\n",
    "        yield alignbatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need a function that iterates through the batches of alignments and make a decision on each of the alignments. In our case it is going to find alignment batches where any alignment has a reference assignment that we are interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def any_align_maps_to_ref(in_bam, refs = [0,]):\n",
    "    # Iterate through groups of alignments from the same read\n",
    "    for alignbatch in get_alignment_batches(in_bam):\n",
    "        # Iterate through alignments from the same read\n",
    "        for align in alignbatch:\n",
    "            if bam.get_ref_index(align) in refs:\n",
    "                yield alignbatch\n",
    "                break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok so what alignment ids are we interested in? They are not human friendly so we will need to look them up in the `ref_to_index` dictionary on the `FileReader`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "HLA_refnames = [x for x in NA12778.ref_to_index.keys() if 'HLA' in x]            \n",
    "HLA_refids = [NA12778.ref_to_index[x] for x in HLA_refnames]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now that we have all the pieces we need lets collect out the reads we are interested in looking at. We are going to use timeit so we have an idea of what the perfomance is like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "for alignbatch in any_align_maps_to_ref(NA12778, HLA_refids):\n",
    "    for align in alignbatch:\n",
    "        NA12778_any_HLA_hits.write(align)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we will close the file to ensure the end of file marker is correctly written\n",
    "NA12778_any_HLA_hits.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That took 20s on a seven year old MacBook Pro for 1.2M reads. Reasonable for interactive use. More typical files of 100M reads you will want to go to lunch. Something like the original 700M read file will take 3-4h."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting AS vs XS scores\n",
    "\n",
    "Now that we have the reads we are interested in one file we can ask a slightly more intersting question - what does the distribution of AS (reported primary alignment) and XS (next best alignment) scores look like for these reads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "NA12778_any_HLA_hits = bam.FileReader(gzip.open('Genomes/NA12778_any_HLA_hits.bam'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "primary_AS_XS = []\n",
    "for align in NA12778_any_HLA_hits:\n",
    "    if not bam.is_flag(align, bam.FLAGS['secondary']):\n",
    "        primary_AS_XS.append((bam.get_AS(align,no_tag=0),bam.get_XS(align,no_tag=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 0), (0, 0), (0, 0), (150, 150), (150, 150)] 10079\n"
     ]
    }
   ],
   "source": [
    "print(primary_AS_XS[:5], len(primary_AS_XS))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import hexbin, figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x11f9776d0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAJcCAYAAADKNbH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XmcHFd57//PU73PaCTNaJvRNiMJW9iObRkLgY0Bh4DjgMG5LAGSgAncQG7yyyWEBBIuISQ3yU0cwhIIYYsxm224hOQmAYxjg23AGEc2trEtyTKyJGu3NPtM73V+f1T1qDWapWemq3tm/H2/XvOSuqq66umqU1VPn3PqtDnnEBEREZH68podgIiIiMhipCRLREREJAJKskREREQioCRLREREJAJKskREREQioCRLREREJAJKsmRes8DnzazPzO5rdjzzjZldaWaHIlz/t83suinm32hmfxHV9hvJzO40s//e7Dgaycx+zcxua3YcC8UzsYzA/D3Pzex9Zva5ZscxFSVZUwhPqD4zS42bvt7M/tnMTprZgJk9YmZviTCOJWa238x+rWpam5kdNLPXNiOmBroCeBmw3jm3Y64rM7MeM3NmFh83fewiMtky45Z/S7jM6+ca03zmnPsl59wXYOwz/6DZMUn9OOe+4py7Kop1h9esl0ax7km2t93M/iO8Zveb2WNm9pdm1h7Of4uZlc1s2MwGzewhM7tm3DqSZvZBM9trZiPhZ7jBzHpq2P5KM/uhmZ0Kt/8jM3tB1Xwzs78ws8PhNfpOM7tgivWdtf/Gn4PT7WMz22Rmvpn943Tx1yrqL3Yz2a5z7q+cc/M66VWSNYnwpHoh4IBXjZv9JeApoBtYAbwJOB5VLM65YeAdwEfNbFU4+Xpgp3Pu682IaSpTJSez0A3sd86NNDmO8a4DeoE3z+RNEcckTaBj2nxmdjlwJ/BD4NnOueXA1UAJuLhq0R8555YAy4FPAreY2fKq+V8nuN7/KrAsfO/9wC/UEMYw8FZgFdAO/A3w71Xl43Xh/BcCHcCPCK7bUXoz0Ae8fnxlgTSIc05/E/wBHyA4YT8M/Me4ecPAthrXswu4pup1HHgaeA6QBr4MnAL6gf8C1kyxrhuBm4Erw/d0zjKmlcB/hNvsBb4PeOG8DcA3whhPAZ8Ip3vA+4EDwAngi8CycF4PQTL6NuAgcHc4/fnAPeF2HgKurIrhLcA+YAh4Evi1CeJ8G5ADyuHn+7Nw+m8CT4Sx/xuwtuo9DvgdYC/w5ATrrMQan2Df/sVUy1Qt2w34wGsILuKdEy1X9Tl/CHwk3J+Vbbw1LBt9wHeA7qr3fIwgYR4kuMC/sGpeJoy1D3gM+EPgUNX89wKHw/26B/iFCWLaFB6TyjH/LHCiav6XgN8L/38n8N+B88Ydi/6q/fYPwDfDbf4Y2DLJvpi0vBPc0P4JOBrG/xdArOq9vxnur6Hwcz8nnH5eGGM/8CjwqnHHdNLYCGpIdwMDwCeAu4D/Hs7bAnw3jPUk8BVgedV794f7+mEgHx6Hfx73ef8e+Ngk+2J/+J6HgZHws68Bvh3GejvQXrX8/wWOhbHeDVww7nN+CvjP8L13cWZ5csD/JDjfTgJ/W3Xs3wL8YNyyv0Vw/vSH+8/CeTHg78J1PAn8f0xynhCUIR/IhuXlPeH0V4XHqT88budNce5Meh5MsOwPgI9Pc90b/1lbwvifG75+aRjvhinWcSfwvwnO6SHgNmDlBMt5wCvD9a+uOje/VrXMBUBuim3tB146zWc4a5mqeQb8DPgfBF+4XzvN/pmujP0F0BruIz88rsNUXX/HXaf+juB+MRAen8x0ZSD8PH9AcF4MAF8luG5MuF3gg8CXq97/5nCbp4A/mWr/NOqvaRue738EN/HfBi4FilQlPwQXwB8CbwA2TrOeDwBfqXr9CmBX+P93AP8enuyxcFtLp1hXO8FN6CTwG+PmzSSm/0NwUU6Efy8MT8gYQTL0kbBQp4Erwve8Ndwnm4ElBInYl8J5PQQXky+G78sA68KC/nKCC87LwterwmUGga3h+7uqT+hxsb6FMy8qLwk//3OAFPBxwqQunO8IbjYdlZN63Poqsc4lyfoT4L7w/z8F3j3Fvn4LQSL2uwQJdga4NtyX54XT3g/cU/WeXyeojYwD7ya48KXDeX9NkBR3ECTEjxAmWcBWgpvS2qrPMVnCcxC4NPz/HoIb8HlV8y4J/38npxOPM45F1X47BewI4/0KcMsk25y0vAP/Anw6LBurgfuAd4TzXkeQeD2XoJw+iyDRTYT78X1AMiwbQ1XlatLYCL5oDAGvDdfzrvA4VT7rswjKbIqgzN4NfLTqs+wHHgyPQYagDI8QJmLh9k5U9vEE+2I/cC9BYrUuXPYB4BKC8+67wJ9WLf9WoC2M56PAg+OOwRDwonD+xzg7cfoeQZnZCDw+2TENl/0PgpqejQRftq4O5/0WQYK7nuBadDtTnyf7qbrBAeeG++hl4T5/T3j8kpO8f9LzYNxyrQTJ/5XTXPfGPitB+fsdoMDpJOivgbumWcedBInLueFxvxP463HLPByu1wGfrZreTZAsnht+/uuBf51iW2fsv0mO11nLVM17IcEXgHaC6+S/T/PZpitjlevjlVR9sZtkXf8Q7pt14b6+PFzvlGUg/Dz3ESRQHQRfrH5rsu1SlWQB5xMkX1cQXA8+RHDvVpI13/7Cg1Qk/IZC8G33XVXz28MT8tHw5H6Q8NvQBOt6FsEFsCV8/RXgA1WF+h7gohnEdjswSliLNMuY/hz4f8Czxk2/jOCiOtE30zuA3656vTXcR3FOJyWbq+a/lzAJq5r2HYJmtlaCbzGvYYJEaNx73sKZF5V/Aq6ver0kjKMnfO2Al0yxvkqs/eP+CtSeZO3ldE3PHwMPTRP/wXHTvg28req1Fx7T7knW0QdcHP5/H+FNL3z9dk4nWc8iuFm/FEhMs1+/BPw+0EmQZF1PcBMdX8t1J9MnWZ+rev1yYPck25ywvBMkGvnqsgC8EfheVbl55wTreyHBjdermnYz8MHpYiP4xntv1TwDDlU+6wTb+mXgJ1Wv9wNvneC4/mb4/2uAx6bY//upqr0F/hn4x6rXv8skN2CCBMhxuib5RqoSW4JzokxYIxMuW11mfhu4Y5LzyxF+sQpffw34o/D/3yVMfMPXL2VmSdafcGZNjkeQPF85VVmd6DwYN319GMezq6ZdT1COR4D3V33WUji9SFAz8itV7/ksk3xBqFrmzsr6qvblrRMslw7L8HVV05IECbAL43gS2DRNGRnmzOvUKLUnWZ+rlCGCa3uRMKGsYV9PVMZqSrLC45qd5FhNWQbCz/Pr447jpybbLmcmWR8Abq6a10JwXW9qkqU+WRO7DrjNOXcyfH1TOA0A51yfc+6PnHMXENwgHgT+1cxs/Iqcc08QZOOvNLMWgqrSm8LZXyK4gdxiZkfM7HozS0wWlJn9OkECcDtBe3/1dmqOiaC54AngNjPbZ2Z/FE7fABxwzpUmeM9agmrYigMECdaaqmlPVf2/G3hd2AG038z6CZLXLhf0r3o9wU39qJl908yePdnnnioOF/RXO0XwjWmiOCaz0jm3vPLH6WMypbAj6ybglnDSTcCFZrZtireNj6cb+FjVfukluMmvC7fxB2a2K+wc20/QlLYyfO/aceur3hdPAL9HcOE5YWa3mNnaSWK6i+Ci9SKCWpo7gReHf993zvlTfJ7xjlX9f5TgJj+Rycp7N8G32qNV++TTBDVaEJTLn02wvrXAU+NiPcCZZWGy2M7Yjy64Ko+9NrM14f47bGaDBM2cKznT+OP6BYLaF8J/p+tvU91nMjvB6yVhLDEz+2sz+1kYy/5wmep4qj/LMEGZWjvRfIJ9NFm5gBr3GbWdZ9XGn7t+uI51Ey08zXlQrY+gGamrat3vCc/rfyG4TlXcG05vJ+hq8MKqeaeq1zGFacu7cy7nnLsZ+CMzq/QJ+wBBbewGgiTsz4DvhveFyfzyuOvUb9cQH2aWIagB/koYz48Iaqh/dZLlayljtVpJ8PkmO2enKwO1Xk8mWnf1eTBKcEybSknWOGHh/BXgxWZ2zMyOETQlXFx1sowJE7EPcbp6cyI3E3yruZbg2+0T4XuLzrk/c86dT1Cdeg2TdKQ2s9UEzXi/SdDs8itm9sKJlp0uJufckHPu3c65zQRJ3++b2S8QFNCNk3TkPUJwM6zYSPBtrPrG4Kr+/xRBTdbyqr9W59xfhzF8xzn3MoKL2m6Cb5G1OCMOM2slaFI4PEkc9XYdQUL0YFg2flw1fTLj43mKoEaget9knHP3hMf0PQRlsD28sA6E24SguXhD1bo2nrEh525yzl1BsI8c45LxKncR3GCuDP//A+AFBEnWXTV+jhmZorw/RVCTVZ34Lg2/MBDO3zLBKo8AG8ys+jq2kTPLwmTO2I/hl5Hq/fpXBJ/3QufcUoKkafwXlvH741+Bi8zs58LP9pUa4qjFrxJcO15KkGj0VMKuWqb6sywhOO+PTDSfYB9Vz6vVUYJao4nWOZHx+2f8uVvZ52cdrxrOg9MbCb60/Rh49fQfYew9wwR9ld5kZpeEk28HdpjZ+snfOWMJgi4WANuArzrnDjnnSs65GwmSvfPruL2K/wYsBT5ZdR9bx+TXqVrKWMV014GTBP03JztnayoDs9juGeUzvJevqGG9kVKSdbZfJqhqP5/gpNhG0Hfm+4QJkJn9jZn9nJnFzayN4GR9wjk3WdZ8C3BVuNxYjYmZ/byZXWhmMYI+SkWCb2QT+QRB1e/3nHNHCS5An608MTKTmMzsGjN7VljAB8LP6xO0hR8F/trMWs0sbacfQb4ZeFf4SPASgpvQVyep9YLgm/8rzewXw29J6fAR3PVhLcG1YYKUJ6gSr7Xm5GbgN8xsW/jZ/wr4sXNuf43vn4lUGHflr4Xgov92TpeNbQRNO786g6fMPgX8sYWPb5vZMjN7XTivjSB5fRqIm9kHCC6WFV8L39se3gx+tzLDzLaa2UvC/ZLjdEfRszjn9obzf52gH8ogQcL8GiZPso4D680sWePnPMNk5T0sz7cBf2dmS83MM7MtZvbi8K2fA/7AzC61wLPMrJvgxjoKvMfMEmZ2JUFn41vO2vjZvglcYGavDo/b/yRoOq1oIyiXA2a2jqCT+pScczmCp9NuIuizd7CGOGrRRnCenCJoAvmrCZZ5uZldER6b/01QY1Nd0/SHYZnZALyToEPxTH0NeKeZrbPgibz3TrP8cU4nGJX3v8LMfiGswXx3+LnumeC9050H470HeKuZ/VH4hZTw/Ng02Rucc70EZesD4evbCfpz/ktY1uIWDJXzW2b21mk+K2b2/MoxMLOMmb2XoKa/8kXsvwhq99eEZfxNnO5XOBeJcdepOEEydQNwIaevUy8gqCy4cIJ11FLGKo4DK8xs2UQzw9qpG4APm9na8Pp/WXhdmkkZmNF2Cc69V5rZ5eF58EEmThIbSknW2a4DPu+cO+icO1b5I0hyfi0swC0E1dD9BH1kujl7mIcx4U3kRwTf3qsvbp0EBWOQoEnxLiZoYjCzXyZoahu70DvnPkfwreAD4aSZxHQOwbe24TCuT4bJW5ngJvUsgqrlQwTNehCcNF8iaFp6kuAm/rtMIrzAX0vQKflpgtqIPyQocx5Bf6AjBM0aLyZICqcVXgj/hKAPy1GCb0tvqOW9szBMkIhU/l4d/vvFcWXjBoImiatrWalz7l8IaphusaBq/hHgl8LZ3wFuJeicfIBgP1ffLP8snP4kQWJSXV5SBP3yThJUua8m6DM2mbuAU1U347sILkoPTLL8dwn6/B0zs5OTLDOVqcr7mwn6rDxG0PzzdcKmG+fc/wX+kiB5GSKoMepwzhUIyusvEXzmTwJvds7tni6QsLb3dQT76xTBOfHDqkX+jODhigGChOwbNX7GLxDc2Or5aP4XCY75YYL9c+8Ey9wE/CnB+XQpp5stK/4fQafrBwk+zz/NIo7PEpS5h4GfAN8iSITKkyz/f4D3W9AE/AfOuT1hXB8nOF6vBF4ZHsfxpjsPzuCc+wHBgw8vAh63oHnxVoJm8I9P8Zk+SpCgXhS+fm34ub5KcOwfAbYTXC+nkyLo8H2K4Fi9HHiFc65Sa/g3BA8WPUhwnX4X8BrnXH8N657KtzjzOvVZgiEnPlp9nXLO3U+wTyaqzaqljAEQnl83A/vCYztR0/MfEDwU9F8EZfJvCPpOzqQMzGi7zrlHCe5JtxDcG4YJ+qjmp1t3lCqP5oqIyByZ2UaC5u/OsHawEdu8kaBD8Psnme+AcyrdFOq43V8i6JTcPe3CIg0Wtrj0E5T9J5sVh2qyRETqwIK+Yb9P8IRaQxKsRgqbwF4eNqOtI6g5+5dmxyVSYWavNLOWsCvKhwhq0/Y3MyYlWSIicxRe1AcJxv/50yaHExUjaEbtI2gu3MXp7goi88G1BN1QjhB0AXiDa3JznZoLRURERCKgmiwRERGRCCzoHzZduXKl6+npaXYYIiIi8gxy//33n3TOrZpuuQWdZPX09LBz585mhyEiIiLPIGZ2YPql1FwoIiIiEgklWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIjLvOedwzjU7jBlRkiUiIiLzlnOOb3zjG3R2drJ582buuOOOZodUMyVZIiIiMi898sgjXHbZZbz5zW/mxIkT7N+/n1e96lVcffXVPPnkk80Ob1pKskRERGRe6e3t5e1vfzs7duzgvvvuY2RkZGze6Ogod9xxBxdccAHvec97zpg33yjJEhERkXmhVCrxD//wD/T09PDFL36RbDY7YT+sUqlENpvlE5/4BBs3buSmm26al/21lGSJiIhI0915551s3bqV9773vQwNDZHP56d9TzabHav1uuSSS3jggQcaEGntlGSJiIhI0xw4cIBXvOIVvOIVr2Dfvn2zav4bGRnh4Ycf5oorruBNb3oTJ06ciCDSmVOSJSIiIg03OjrK+973Ps477zxuu+02RkdH57Q+5xzZbJavfe1rbN68mQ996EMUCoU6RTs7SrJERESkYZxzfPWrX6W7u5uPfvSjZLNZSqVS3dZfKBQYGRnhgx/8IFu2bOHWW2+t27pnSkmWiIiINMRPfvITLr30Ut72trdx8uRJstlsZNsaGRnh0KFDvOY1r+ElL3kJe/fujWxbk1GSJSIiIpHq7e3luuuu4wUveAEPPvhgQ4ddGB0d5e677+biiy/mXe96V0O3rSRLREREInX99ddz0003TTokQ9TK5TLZbJZPfepTfOELX2jYdpVkiYiISKTy+Xxd+13Nlu/7De0MryRLREREJAJKskREREQioCRLREREJAJKskREREQioCRLREREJAKRJVlmdoOZnTCzR6qmfdDMDpvZg+Hfy6vm/bGZPWFme8zsF6OKS0RERKQRoqzJuhG4eoLpH3HObQv/vgVgZucDbwAuCN/zSTOLRRibiIiISKQiS7Kcc3cDvTUufi1wi3Mu75x7EngC2BFVbCIiIiJRa0afrP/PzB4OmxPbw2nrgKeqljkUTjuLmb3dzHaa2c6nn3466lhFREREZqXRSdY/AluAbcBR4O9mugLn3Gecc9udc9tXrVpV7/hERERE6qKhSZZz7rhzruyc84HPcrpJ8DCwoWrR9eE0ERERkQWpoUmWmXVVvfxvQOXJw38D3mBmKTPbBJwD3NfI2ERERETqKR7Vis3sZuBKYKWZHQL+FLjSzLYBDtgPvAPAOfeomX0NeAwoAb/jnCtHFZuIiIhI1CJLspxzb5xg8j9NsfxfAn8ZVTwiIiIijaQR30VEREQioCRLREREJAJKskREREQioCRLRERkkRoeHmZ4eLjZYdDR0cHWrVubHQbOuYZuL7KO7yIiItIc+XyegwcPMjQ0BMDSpUvZuHEjyWSyoXFks1n27NnDFVdcwXOf+1zuu+8+/v7v/55Tp041NA6AlpYWli5dyste9rKGbVNJloiIyCJRLpc5evQox48fB07X3PT39zM4OMiaNWvo6urC86JtyCqVSuzfv5/Dhw+PxZBKpXjBC17Ac5/7XG655RZuvvlmisVipHEAJJNJ4vE4/+t//S/e/e53k0qlIt9mhTW66qyetm/f7nbu3NnsMERERJrKOUdfXx8HDx7E9318359wOTMjFovR3d3N8uXLMbO6x3H8+HH27t07ZRyFQoHh4WE+/OEPc88999Q1hmqZTIZrrrmGj33sY3R1dU3/hhqZ2f3Oue3TLqckS0REZOEaHR1l//795HK5SZOa8TzPI51O09PTQ0tLS13iGBoaYteuXWSz2ZrjyOVyPPnkk1x//fUcOHCgLnEAtLa20t3dzQ033MDznve8uq23QkmWiIjIIlYsFnnqqafo6+ubdYduM2PFihWsX7+eeHx2PYgKhQJ79+7l5MmTNSdX1ZxzFAoFbr/9dj796U/PqaN+Op0mnU7zkY98hDe/+c2RNYsqyRIREVmEfN/nxIkTHDlyBOfcnJ+YMzPMjHXr1rF69eqamxB93+epp55i//79dYmjVCqRz+f5zGc+wze/+c0ZJWyxWIxkMsk73vEO/vzP/5y2trY5xTIdJVkiIiKLzMDAAAcOHKBUKs2q1mgqnucRj8fp6elh6dKlUy576tQp9uzZQ7FYrHsc+XyekydP8rd/+7c8/PDD0y7f2trKjh07+PSnP80555xT11gmoyRLRERkkcjlchw8eJDh4eG6JzXjeZ7HkiVL6O7uPutJvNHRUfbs2cPg4GDkceTzeR588EE+8pGPcOLEibPmt7a20tHRwWc+8xmuvvrqSGMZT0mWiIjIAlculzly5AgnTpxo+ECaZsbq1atZu3Ytzjn27dvH0aNHI0+uqvm+T7FY5Otf/zpf/vKXyefzJJNJEokEH/zgB3nnO99JIpFoWDwVSrJEREQWsOHh4bGhEJp1rzYzfN8nm83inGtoglWtUCiQzWZ5//vfz7Zt2/jwhz/M6tWrmxIL1J5kaTBSERGReWhgYIByudzUGJxz5PP5pseRTCZJJpN87nOf4/nPf35TY5kJ/XahiIiIzHtmxtq1a5sdxowoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERERkQgoyRIRERGJgJIsERGRecjMmh0CMH/iAHDO4Zxrdhg1U5IlIiIyjzjnKJfLFItFCoVCU5MKM2Pt2rV0dXXhec1LGTzPo7OzkxUrVpDL5SgWiwsi2Yo3OwAREREJlMtlBgcH2bVrF6Ojo0CQYLS0tBCLxRpWq+R5Hq2trXR3d5NOpwHYsGEDu3fvZnh4GN/3GxZHW1sbmzZtGosDoFQqUS6XSSQSxGKxhsQyG0qyREREmqxcLlMqlXjsscc4derUGfN832d4eJh4PE5LS0ukNUqe5xGLxejp6WHZsmVnzGttbeU5z3kOJ0+eZM+ePZTL5ciSLc/ziMfjbNmy5aw4KpxzFAoFPM8jkUg0taZtMkqyREREmsT3fZxz7Nu3j4MHD07ZBFYqlRgcHCSVSpFOp+taq2VmmBldXV2sWbNm0oTFzFi1ahUdHR089dRTHDhwoK79pCpxbNiwgc7Ozpo+o+/75PN5YrEYiURiXvUhU5IlIiLSYM45fN/n6aef5vHHH6dQKNT83nw+T6FQoKWlhXg8Puekwsxob29nw4YNJBKJmt5Tqe3q6upi7969nDp1as61WmbGihUr6O7urjmOauVy+YwmxPmQbCnJEhERaaByuUw2m+Wxxx5jcHBwVutwzjEyMkIsFhtrQpxpUuF5HqlUip6eHlpbW2cVRyqV4ud+7ucYGBhg9+7d5HK5GSdbnueRyWTYvHnzrOOoViwWKZVKJJPJpjchKskSERFpAN/3KZfL7Nmzh2PHjtVlneVymaGhIRKJBC0tLTUlWmaG53ls3LiRjo6OutT4LFu2jB07dnD06FF+9rOf4fv+tMmW53l4nsemTZvqFkeFc25eNCEqyRIREYlQpd/VoUOH2LdvH+Vyue7bKBaLDAwMkE6nSaVSkyYVZsaaNWvo6uqq+1N5leEeVq9ezb59+zh69OikiZbneXR1dbF27dpInw6sNCHG4/G6NK3OlJIsERGRiJTLZQYGBti1axfZbDby7eVyOQqFAplM5oykwvM8lixZQnd3N6lUKtIY4vE45557LuvXr2fPnj0MDg6OJVue57F06VI2bdoUeRzVSqXSWBNiI4d8UJIlIiJSZ6VSiWKxyK5du+jt7W3otn3fZ2RkhHg8TltbG4lEgp6eHpYuXdrQOFpaWrjkkks4deoUu3fvxszYsmVLw+OoVigUSCQSxOONSX+UZImIiNTZwYMH2bdvX1NjKJVKtLe3s3nz5qbGsWLFCp773OdG0kw6382/kbtEREREFgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiIiIRUJIlIiIiEgElWSIiC5hzDudcs8OQcQqFAsVisdlhAMyL8lEqlSiXy80Oo+HizQ5ARERmzjlHsVgkl8vhnCOVSpFKpTCzZof2jOb7Pt///vf56le/SrlcZtu2bWzZsgXPa3ydxqpVq2hvb2dwcJB0Ok0ymWx4+fB9n/3797Nnzx7MjHPOOYeurq6mltNGbtvmQ4Y7W9u3b3c7d+5sdhgiIg1VLpcZHR3F9/0zppsZmUyGRCLRpMie2fbu3cvnP/95+vv7yefzAMTjcdLpNM973vNYvXp1Q+Joa2ujp6eHZDJ5RnLneR6ZTIZ4vDH1KydPnuShhx6iUCiM1WLFYjFSqRTnn38+y5Yta0gcFWZ21j6Zw7rud85tn3Y5JVkiIguD7/vkcrlpm6FHPjjMAAAgAElEQVRisRiZTIZYLNagyJ7Zent7+cpXvsJjjz1GoVCYcJlYLEZnZyeXXnopra2tkcSRTCbp7u6mra1tymMfj8fJZDKR1a6Njo7y05/+lN7e3kmbCD3PY8WKFWzdupVUKhVJHNUqyVW9arGUZImILBLOOfL5/FjtSK2SySTpdFpNiBEpFAp8+9vf5tZbb6VUKp1Vszhe5Sa/detWLrjggrrVKHmeR1dXF52dnZhZzce73k3MpVKJxx9/nP3799fUV7ASa09PD93d3ZEkffF4nHg8XvdzQEmWiMgiUCwWyWazc+q83Kz+OIuVc47777+fL3/5y+Tz+UlrryYTj8eJxWJs376dDRs2zOm4tLe3093dTSwWm1WSUmlinksi4pzjyJEjPPLII5TL5WmTzfE8zyMej/PsZz+blStX1qWcep4XaZlXkiUisoCVy2Wy2WzdnshqdH+cxerQoUPceOONHDlyZMY1i+PF43GWLl3Kjh07aG9vn9F7M5kMmzZtIp1O16VZeLZNzP39/Tz88MOMjIzMuax6nkdbWxvPfvazWbJkyazWUc9+V9NsR0mWiMhC45wjl8vNuHakVlH3x1mshoeH+frXv869995LqVSq67AIsViM7u5utm3bNm3/pHg8zoYNG+jo6JhR02CtEokE6XR62vKRz+d59NFHOXbs2IxrrqZTaf7csmXLjB7iSCQSxGKxhtTY1ppk6SuNiMg84JyjUCiQy+Ui3U6pVGJoaEhDPtSoXC5z55138s///M+Uy2VKpVIk2zhw4AAHDx7kwgsv5Nxzzz0ryTEzVq1axfr16+vagXu8YrFIsVictInZ93327dvH3r178X0/kjG4fN/n6NGjHDt2jGc961msW7duys8bi8VIJBLzsiyrJktEpMlKpRLZbLbuNQLTqUd/nMVs9+7dfP7zn2doaGjOTYO1isfjpFIpduzYQWdnJwBLly6lp6eHRCLR0BpIM6OlpWWsifnEiRM8/PDDFIvFhg0sGovFSCaTnHfeeWc1qXqe1/B9UqHmQhGRec45x+joaCS1IzMRi8VobW1VohUaGBjghhtu4PHHH4+s2XY6sViMTZs28cY3vpGlS5c2dTiOUqnEnj17GBgYaNqo7Z7n0dHRwXnnnUcymSSZTDZ1n9SaZKlRXkSkSYrFYtMTLEDjaY3zwx/+cMoxrxqhXC6zefNmli1b1vTjc+jQoSnHvGoE3/c5deoUJ06cqFtn/0ZQkiUiIlJlvvwepJnNiwcU5sO+qIiyP1oUmn/0RERERBahyJIsM7vBzE6Y2SMTzHu3mTkzWxm+vtLMBszswfDvA1HFJSIiItIIUQ7hcCPwCeCL1RPNbANwFXBw3PLfd85dE2E8IiIiIg0TWU2Wc+5uoHeCWR8B3gPMn0ZeERERkTpraJ8sM7sWOOyce2iC2ZeZ2UNm9m0zu2CKdbzdzHaa2c6nn346umBFRERE5qBhSZaZtQDvAybqb/UA0O2cuxj4OPCvk63HOfcZ59x259z2VatWRROsiIiIyBw1siZrC7AJeMjM9gPrgQfMrNM5N+icGwZwzn0LSFQ6xYuIiIgsRA377ULn3E+B1ZXXYaK13Tl30sw6gePOOWdmOwiSv1ONik1ERESk3qIcwuFm4EfAVjM7ZGZvm2Lx1wKPmNlDwN8Db3DzafQzERERkRmKrCbLOffGaeb3VP3/EwTDPYiIiIgsChrxXURERCQCSrJEREREIqAkS0RERCQCSrJEREREIqAkS0RERCQCSrJEREREIqAkS0RERCQCSrJERGTe0DjUspgoyRIRaZJ4PI7nNf8yXCwW8X2/qQlOuVwmn8+Tz+cpl8tNi8P3fS666CI6OjpIJBJNiyOVSnHw4EF8329aDBWrV68mk8k0tazGYjGWLFlCV1dX02KYDVvI3xq2b9/udu7c2ewwRERmzTlHqVQim802JcnxPI9YLIaZEYvFxhILM2vI9ivJ3eHDh3n66acB6OjoYMOGDXie17AbeyWOY8eO0dfXh3OORx99lLvvvptyuUypVGpIHIlEgng8zute9zpe+MIXYmYUCgVyuVxDtj9eLBYbOwaHDx/miSeewDnXsOSvUj7PP/981q9f37ByOR0zu985t33a5ZRkiYg0n3NurCanEaqTq/Hi8TjxePCra1Hd1JxzOOfo7e3l0KFDZ9VeeZ5HV1cXq1evxswij6O/v59jx46dlTzk83nuuecefvrTn0Za22dmxONxrrjiCl796lfT0tJyVpy5XI5CoRDJ9sebrHwUi0WeeOKJCfdVPVWOeU9PD+eee+5YeZwvlGSJiCxAvu+TzWYjrTmppZnSzEgkEnieV/cEp1wuk8vl2L9//7Q1NMlkku7ublpbW4nFYnWNw/d9crkchw8fnja57evr4/bbb+f48eMUi8W6xlH5jNddd920zWHlcplsNhtZk2ol2ZvumA8PD7Nr1y5GRkbqHkssFqO9vZ0LL7yQ1tbWuq67XpRkiYgsYJUmxHrWFlSafmaSNHmeRyKRqEttku/7lMtlDh48SH9//4ze29bWRk9PT136sVXiOHLkCENDQzN675NPPsntt99OPp+fc7KVTCbJZDJcd911XHTRRTXv36iamGe6b51znDx5kt27d1MqleZcVivN1RdffDGrVq2a07qipiRLRGSBc87VpT/OVE2DtZpLf61KH57jx49z7NixOSUGq1evZu3atZjZjJOtStPg008/zcmTJ2cdR7lc5oEHHuDHP/7xWMI2E7FYjFgsxjXXXMNVV1016w729Wpinmv5KJfLHDhwgAMHDozt45lu38zYunUrPT098+JhkOkoyRIRWSRm2x+n1qafmUgkEmPNdrWs1/d9BgcHOXjwYN2a2eLxOOvWraOjo6OmGrbKjX9oaIijR4/WrSl2ZGSEu+++myeeeKLmdSYSCbZt28Yb3/hGli1bVpc4ZtvEXO/ykcvlePzxxzl16lTNtVqe57F27VrOO+88UqlUXeJoBCVZIiKLzEz640Q5PEQt/bXK5TLFYpH9+/czMjISSRyZTIaenh5SqdSk/bV836dQKHD48GGy2WwkcRw/fpzbbruNgYGBSRPJVCrFihUr+I3f+A02b94cSRwzaWKOsnz09/eza9euKYfjiMVitLa2cvHFF9ct2WwkJVkiIovQdP1x6tE0WCvP80gmk8DpWq3KE3iHDh3i5MmTkccA0N7ezsaNG88Y8qESx9GjR2fc/2s2nHPs3r2b733ve2cM+ZBIJEgkErzhDW/gsssui7wpbLom5tn0y5ttHEeOHGHv3r1nDPlQ2f6FF15IV1fXvBmSYaaUZImILGLj++M0MrkarzLkg3OOU6dOcfjw4YYPKGpmdHV1sWbNGiB4GvD48eMNH8yzUChw77338tBDD2FmXHnllVx77bVkMpmGxjG+iblZ5aNUKvGzn/2MI0eOALB582bOOeecuj8p2mhKskREngEKhQKlUqnpNQJ9fX2cOnWqYeM4TSadThOPx+s+zMJMJRIJVqxYwcqVK5saRy6Xw/f9ppePQqEw9jTlYlBrkjW/RvcSEZEZaUTTTy0qfZ+arVgszovfP1y+fDkrVqxodhh4njcv9kdLS0tTf6aoWeb/c5IiIiIiC5CSLBEREZEIKMkSERERiYCSLBEREZEIKMkSERERiYCSLBEREZEIKMkSERERiYCSLBEREZEITJtkmdkLapkmIiIiIqfVUpP18RqniYiIiEho0p/VMbPLgMuBVWb2+1WzlgIL+5cdRURERCI21W8XJoEl4TJtVdMHgddGGZSIiIjIQjdpkuWcuwu4y8xudM4daGBMIiIiIgveVDVZFSkz+wzQU728c+4lUQUlIiIistDVkmT9X+BTwOeAcrThiIiIiCwOtSRZJefcP0YeiYiIiMgiUssQDv9uZr9tZl1m1lH5izwyERERkQWslpqs68J//7BqmgM21z8cERGZCc/zSCaTFItFnHNNicH3feLxOKtXr6avr49isdiUOFKpFBs3biQWi3HkyBHy+XxT4ojH4yQSCUZGRshkMsRizRv1KB6PE4vFmlo+zIxYLIZzDjNrSgzNMm2S5Zzb1IhARERk5swMMyOVSlEulxua4DjnyOfzFItFPM8jlUqxZs0aRkdH6e/vx/f9hsTheR5r165l1apVYzfxLVu2MDAwwLFjxyiXG9Od2MxYtmwZra2tmBnlcpnh4WGSySTpdLopCUZ1+SiVSpRKpYZuP5FINDXJbLZpkywzawF+H9jonHu7mZ0DbHXO/Ufk0YmIyJSqb9yxWIxYLBb5zdQ5R7FYnLCmyMxoaWkhk8kwMDDA8PBwZHEArFixgvXr1+N5Hp53ugdMJeFZunQpx48fp7e3N9I4WltbWb58+VhSU61QKFAoFEin0ySTyYYmW9XbisfjxONxisVi5IlnLBYjkUicFcMzTS3NhZ8H7icY/R3gMMETh0qyRETmkcrNrLqJqN61SeVymVwuN+V6K4nGsmXLaGtro7e3t+5Nd62trXR3d5NMJietKakkXZ2dnaxcuZLDhw8zMjJS1ziSySQdHR3EYrEzkryJ5HI5CoUCmUyGeLyW2299VcpHIpEYS7bqXT48zyORSEyYbD4T1XKUtzjnXm9mbwRwzo2a9pyIyLxVucElk0l8369Lfxzf98nn8zOqIavULq1cuZJCoUBvb++ca1ASiQQbNmxg2bJlNd/IK/3Wuru7GRkZ4ciRI3NuVo3FYixfvpx0Oj1tclXN931GRkaIx+NkMpkZvbdeqstHvZqYzYxEIoHneUquqtSSZBXMLEPQ2R0z2wI0pzehiIjUzMzG+krNtgnROTfW3DVblRg6OzsZHh5mcHBwxkmfmbFmzRo6OztnfSP3PI8lS5ZwzjnncOrUKU6cODGrONra2mhra5tTbU2pVGJoaIhUKkUqlWpaf616NDFXmiEr65TTakmy/hS4FdhgZl8BXgC8JcqgRESkPqqbEGfaH6dUKpHL5eryVFolIWlra2PJkiX09fUxOjpa03uXLVtGd3d3TU1ytcaxYsUK2tvbOXr0KAMDAzW9N5PJ0N7eXtfamnw+f0YTYqOTlPHlo1Ao1NyEqH5X07NaTh4zWwE8HzDgXufcyagDq8X27dvdzp07mx2GiMiC4Zwbq52a7Prv+z65XC7SztG+71Mul+nt7Z20liydTtPd3R3pMAiVZtDDhw+Ty+UmXCaRSNDe3j7WHBYVz/NoaWlp6tN4zrlpm5grTY3P5H5XZna/c277dMvV2vNuHRALl3+RmeGc+8ZcAhQRkcar3BQnGvKhekiGqFX6a61evZpsNktfX99YDUosFmPdunWsWLEi8hu553mk02k2b97M4OAgR48eHUsuPc9j2bJltLS0NCSh8H2f4eFhEonEjPt61Ut1E/NE/bWqh2R4piZYM1HLEA43ABcBjwKVOkQHKMkSEVmAKjfHSn+cYrFINpttyuCdZkYmkyGdTjM4OEgqlWL9+vVjN/tGxVD9NOSJEyfI5/Mz6lxfT8VikWKx2JQhH2Di8mFm6nc1C7XUZD3fOXd+5JGIiEhDVW6W5XK5aaOjV+IwMzo7O1myZElTanAqccRiMdrb28nlck1PJho1mOtkqod8qH4ttaulJP/IzJRkiYgsUs36uZXx5lMfn/kQx3zZH/MljoWolpqsLxIkWscIhm4wwDnnLoo0MhEREZEFrJYk65+ANwE/5XSfLBERERGZQi1J1tPOuX+LPBIRERGRRaSWJOsnZnYT8O9UjfSuIRxEREREJldLkpUhSK6uqpqmIRxEREREpjBtkuWc+41GBCIiIiKymNQyGOkq4DeBnurlnXNvjS4sERERkYWtlubC/wd8H7gdiO6HrEREREQWkVqSrBbn3Hsjj0RERERkEallxPf/MLOXRx6JiIiIyCJSS5L1ToJEK2tmg2Y2ZGaDUQcmIiIispDV8nRhWyMCEREREVlManm68DkTTB4ADjjnSvUPSURERGThq6Xj+yeB5xD8diHAhcAjwDIz+x/OuduiCk5ERERkoaqlT9YR4BLn3KXOuUuBbcA+4GXA9VEGJyIiIrJQ1ZJkneuce7Tywjn3GPBs59y+6MISERERWdhqaS581Mz+EbglfP164DEzSwHFyCITERERWcBqqcl6C/AE8Hvh375wWhH4+agCE5HFyTnX7BAAxSEi0atlCIcs8Hfh33jDdY9IRBal6mTCOYeZKY55EIfv+5gZ5XIZM8PzavnuHY1SqfkPrDvn8Dxv7Pg067gAlMvlppYNmbtJkywz+5pz7lfM7KfAWV+1nHMXRRqZiCwazjl836dYLOKcI5FIEIvFgMbexJxzlMtlstkszjnS6TSJRKIpcRSLRYaGhvB9nyVLlpBOpxseR7FY5OjRo/znf/4nuVyO7du3s2HDBmKxWMNv7Mlkkra2tqYmFL7vUygU6O/vx/d9kskk8Xi8KTHF43EymYwSrAVuqpqsd4b/XtOIQERk8anUBhQKBXzfH5teLBYplUokEgk8z4v8RuKcwzlHNps9o7Ykm82Sz+fJZDINSSwqyebQ0BCFQmFs+tDQENlslra2tobc1IvFItlslttuu40DBw6MTf/BD35AR0cHl19++VgsUfM8j6VLl5JMJiPf1mR838f3ffr6+s44Lvl8nmKxSDqdxswakvB4nkcmk2nIvpfo2ULuD7B9+3a3c+fOZochIuNUriulUmnaJiDP88ZusPW+iVXiyOVyZ9w8J1JdcxBFHM45RkZGyGazUy6bSqXGanTqHUe5XKZcLnPPPffwk5/85IzEd7yenh527NhBPB4fq3WsJzOjtbW1qbU1leMyMDDA6OjolMvG43FSqRQQXW1jOp0mmUyq9moBMLP7nXPbp1tuqubCISZoJgQMcM65pXOIT0QWqUqTXLFY28PHvu+Ty+WIx+Nj397nepOpJFfFYpFcLldT5/JSqcTQ0BCpVKpuN9PqJG94eLimOPL5PPl8npaWFlpbW+sWR6lUYu/evdx5553TJnoA+/fv59ChQ1x00UVs3boVz/Pq1l8rnU6zZMmSpvX/qiRXo6OjDA4O1lw+SqUSyWSSRCJR10QomUyO1ZbJ4jJpkqXfLBSRmajcuIrF4pQ1JJOp3MTm2l+r0iSXzWYpl8szfn8+n6dQKIw12cz2xuf7PuVymaGhoVl16B4dHSWXy9HW1jan2o1isUhfXx/f+c53OHHixIzeWyqVeOCBB3j88cd53vOex+rVq+fUjBWPx1m6dGlTm8J836dUKtHX1zer41IoFMaaEOfa1B2LxcaaqmVxqrmkm9lqIF157Zw7GElEIrKgVNcazSapGa/SX2umTYiVJC+Xy9VcizbVukZHR8dugjO5mVbiGBoaIp/PzykO3/cZGBggkUjQ1tY2o35jpVKJYrHIHXfcwZ49e+YUx/DwMHfccQednZ1cdtllpNPpGSVKZkZbWxupVKqpTYO+79Pf308ul5vzurLZLLFYbOwzzeRzmdmck3hZGKbtk2VmryIYvmEtcALoBnY55y6IPrypqU+WSPPMpN/VbMVisWmf/qvuXD/Xm+dkEokEmUympjhGR0cZHR2NZPyrSjPbVDf1Sg3aAw88wL333lv3Y2NmnHvuuVxyySXEYrFpm/wqzZ7N7nc1NDTE8HA0ow4lEomavxRUmqOVXC1sc+6TVeV/A88HbnfOXWJmPw/8+lwDFJGFa/yQDFGpdNSerL9W9ZAMs2mirFWxWBxrIproZuqco1AojA3JEJVcLkc+nx/rMD4+jmKxyKFDh7j99tsZHByMJAbnHHv27OHJJ5/k0ksvpaenZ8IatsqQDM1qCqskV/l8fmxIhqhU175OVjtVebCimeOQSePVkmQVnXOnzMwzM8859z0z+2gtKzezGwiGgDjhnPu5cfPeDXwIWOWcO2lBqfwY8HJgFHiLc+6BGX0aEYlMJZmaS7+r2SqVSpTL5bEhHypxjB+SIWqVpxQr/WgqN/LBwcE5N1HWyjnH8PDwGUM+lMtlRkdH+c53vsNTTz3VkDgKhQI/+tGP2LVrF5dffjnLli0bexKx0o+sGSrHpFwu09/fP+1TpfXcbmXIh1QqNdbErCEZntlqOer9ZrYEuBv4ipmdAEZqXP+NwCeAL1ZPNLMNwFVAdb+uXwLOCf+eB/xj+K+IzAOVZqh69LuajUptUUWjbp7j+b7PyMjI2KjgUTVRTqeSRDz99NMcPnyYRx55pCk/0dPf38+3vvUtXvziF3PeeefR0tLS1KawXC5HNput6QnKKFQeuqgMMKshGZ7Zaqm3vJagZuldwK3Az4BX1rJy59zdQO8Esz4CvIczh4i4FviiC9wLLDezrlq2IyLRq9QONFu5XG5aglWtMjxEs/X29vLYY481/TcQh4aG5sUwBKVSqWkJVjUzU4IlNf12YaXWyge+MNcNmtm1wGHn3EPjCt86oLqe+1A47ei4978deDvAxo0b5xqOiIiISCQa2gPPzFqA9wEfmO06nHOfcc5td85tX7VqVf2CExEREamjRvfE2wJsAiq1WOuBB8xsB3AY2FC17PpwmoiIiMiCM21Nlpm9s5ZptXDO/dQ5t9o51+Oc6yFoEnyOc+4Y8G/Amy3wfGDAOXd0qvWJiIiIzFe1NBdeN8G0t9SycjO7GfgRsNXMDpnZ26ZY/FvAPuAJ4LPAb9eyDREREZH5aKofiH4j8KvAJjP7t6pZS5n4icGzOOfeOM38nqr/O+B3almviIiIyHw3VZ+sewie7FtJ8LM6FUPAw1EGJSIiIrLQTZpkOecOAAfM7KVA1jnnm9m5wLOBnzYqQBEREZGFqJY+WXcDaTNbB9wGvIlgJHcRERERmUQtSZY550aBVwOfdM69Drgg2rBEREREFraakiwzuwz4NeCb4bTm/Ky6iIiIyAJRS5L1e8AfA//inHvUzDYD34s2LBEREZGFrZbfLrwLuCv8SRycc/uA/xl1YCIiIiILWS0jvl9mZo8Bu8PXF5vZJyOPTERERGQBq6W58KPALwKnAJxzDwEvijIoERERkYWuliQL59xT4yaVI4hFREREZNGYtk8W8JSZXQ44M0sA7wR2RRuWiIiIyMJWS03WbxH8puA64AiwDf3GoMgzjud5xOO1fC+LlpkRizV3FBnnHNlsllwuR/Czq83T1dXFjh07mnpsPM9j3bp1mFnTYqhIJBIsWbKk2WGQSCSaHYLMA9bsC8RcbN++3e3cubPZYYg8I1RfK4rFIuVy83oNVGIpl8v4vt/QbRcKBYaGhs74/Ol0mmQy2dAkI5FI0NbWhud5+L5PsVjku9/9Lrt3725YDAAbN27kF3/xF2lpaSEWizU90aqUDd/36e/vJ5fLNXT7qVSK9vZ2PC+ow2j2/pBomNn9zrnt0y43XZIVjov1MeD5gAN+BLwrHMqhqZRkiTSHcw7nHMViseFJzkRxlMvlyGuUyuUyQ0NDFAqFCed7nkcmk4m8RsnzPNra2iZM6orFIv39/dx6662cOHEi0jiWLVvGy172MtauXTtva21836dUKtHX10epVIp0W/F4nOXLl5NIJMYSLFm86plk3Qv8A3BzOOkNwO8655435yjnSEmWSHNVEpxSqdTUZjPnHL7vR1K75pxjZGSE0dHRmpaPx+NkMplIbrStra20tLQAk9eQOOcolUo88cQTfO973yObzdY1hkQiweWXX862bdvwPG/eJxSVcjk6OsrAwEDdy6mZsXTpUlpbW8dey+JXzyTrYefcReOmPeScu3iOMc6ZkiyR5qtcQ0qlUuS1BbXEUq8mROcc+XyeoaGhWd2YU6kUqVSqLjfdVCpFW1sbZlbz+ir74Z577uGBBx6oyz4577zzeMlLXkI8Hp8X/fNmolLrOTg4yMjISF3W2dLSwrJly2Z0XGRxmHOSZWYd4X/fC/QBtxA0F74eaHfO/XGdYp01JVki80flWlIoFJrehAjMqXatVCoxODg456TRzMhkMrNuTovH47S1tRGPx2d9Ey8Wi+RyOW677Tb2798/q3WsWbOGq6++mmXLls3bpsFa+b6P7/v09fVN2vQ7nWQyyfLly4nFYvO+Jk+iUY8k60mCpGqiM9s55zbPLcS5U5IlMv9Umu6KxWLTmxArTWe18n2f4eHhuneWjsViZDKZmp+KNDOWLFlCOp0eez1XxWKRY8eOcdttt9Hf31/Te1paWvj5n/95tmzZMqdEbz7yfZ9CoUB/f3/Nzcye57F8+XJSqZSSq2e4ujUXzmdKskTmp+qn/4rFYtNjma6/VmVIhpGRkUgTw2QySTqdnjJZyWQytLa2RtIEVdkPDz/8MPfcc8+UnfgvvfRSLrvsskVdW1NJxEdGRqZtFl6yZMmMm2xl8ao1yVpYjeoisiBUbkKxWIxYLNbUIR/MbKyD9kT9tQqFAoODg2M33CgVCgUKhcKEQz4kk8mxIRmiuolX9sO2bdu44IILuPPOO3n00UfPWGbTpk1cddVVpFKpBd80OJ1KwlR5oGBgYOCsBwXS6TTLly+P9LjI4qWaLBGJ3Hwa8gFOd9KfakiGqFWGfKh0ak8kEg2/iReLRQYHB7n11lvJ5/NcddVVrFmzZtEnV5OpDPnQ39+Pc4729nbi8fiircmT2VNNlojMG5Uag3g83rSkphIHBM2Yp06daloccLoDdnt7e9OaoBKJBB0dHbz+9a8fi+GZnFB4nkcikWDlypVjx0O1VzIX055NZnZHLdNERKYzX2rOy+XyvLh5VmqMmhlLJfldCGNeNUJ1ojkfyogsbJPWZJlZGmgBVppZO6efMlxK8DuGIiKySCihOJP2h9TDVM2F7wB+D1gL3M/pJGsQ+ETEcYmIiIgsaJMmWc65jwEfM7Pfdc59vIExiYiIiCx4tTTArzazsRH0zGypmX0+wphEREREFrxakqwYcJ+ZXWRmLwP+i6D5UEREREQmMe0QDs6594VPE/6Y4DcMXzxJvG4AACAASURBVOSceyLyyEREREQWsFqGcHgR8PfAnwN3Ah83s7URxyUiIiKyoNUyGOmHgNc55x4DMLNXA98Fnh1lYCIiIiILWS1J1mXOubEfHXPOfcPM7oowJhEREZEFr5aO7yvN7J/M7FYAMzsf+OVowxIRERFZ2GpJsm4EvgN0ha8fJxikVEREREQmUVNNlnPua4AP4JwrAeWp3yIiIiLyzFZLkjViZisAB2BmzwcGIo1KREREZIGrpeP77wP/Bmwxsx8Cq4DXRhqViIiIyAJXy2CkD5jZi4GtBD8Svcc5V4w8MhEREZEFrJaaLIAdQE+4/HPMDOfcFyOLSkRERGSBmzbJMrMvAVuABznd4d0BSrJEREREJlFLTdZ24HznnIs6GBEREZHFopanCx8BOqMORESkUTyvlktf9Mrl+TEaju/76Hu0SP3VUpO1EnjMzO4D8pWJzrlXRRaViCwqlRu47/tNjiSQyWRYsWIFfX19TUswzAwzY3R0lJaWlrHXjVT57MeOHaNcLtPV1UUsFps3SajIQldLkvXBqIMQkcWpkryUy2WKxfn1UHI6naazs5Ph4WEGBwcbmmhlMhlaW1vxPI9yuczQ0BDJZJJ0Og0QebLlnMM5R29vLydPnhxLfvfu3UtHRwerV69uStInstjUMoSDfgxaRGbMOYfv+xSLxXnbFGVmtLW10dLSwsDAAKOjo5FuL5FI0NbWRjx+9qW3UChQLBZJp9MkEomx+OrN932y2SxHjx6lUCicNb+3t5eBgQE6Oztpa2tTsiUyB7UO4SAiUpNKQlUoFOZN8+B0YrEYHR0dtLW10dvbS6lUqmti6HkebW1tpFKpKZdzzpHNZsnn82QyGWKxWN0SHN/3KZfLHDlyhJGRkSmXLZfLHD58mHQ6zdq1a0kmk2pCFJkFJVkiUheVpKRUKlEqlZoczewkEglWr15NNpulv7+/Lklia2vrWJ+rWvm+z8jICIlEgnQ6PafapErT4IkTJ+jt7Z3Re3O5HPv27WPp0qV0dnbieZ6SLZEZqCnJMrMMsNE5tyfieERkAXLO/f/t3XmUXGd95//Pt25VV1V3a2tZm5ElWf4J2QYb2ygOeAjHwByGeBw7M5gtkMQwAZwQEpJfAj/CDGRgmJkMk5NfmAx7OAaGsAUIjgeIzWIgCR6QjbwlBsu2JMtarKX3pdbv/FH3tkrtXqq766ml9X6d00dVt7ZvPVXq++nnee5zO3Le1VKYmXp7e5XP5zUyMqKxsbEl9Wpls1n19/criqIl11IqlVQqlZTNZqd7wRoNW0m4GhkZ0fHjx5d1JOPIyIhGR0e1YcMGDQwMMIQINGjBP0nM7JdUW4j0m/H1K8zsttCFAeh8ybyrZD7RSmJmWrNmjTZt2jTdm9TIY6Io0rp167RmzZplBax6hUJBo6OjDQ9jVqtVTU1N6cCBAzpy5EhTlopIesMeffRRjY+Pd81QMNBOjR5deLWkuyTJ3feZ2YUBawLQ4ZIdfalU6pi1nkJJp9M677zzVCgUNDg4qEqlMmvQMTP19/c3HMgWy901MTGhKIqUz+eVSqWe9jrValXValXHjh3TyMhI02uQap/5oUOH1NfXpy1btiidTjOECMyhkZBVcvfhGf+ZO/NQIQDB1AeLbp53tVTZbFabNm3S+Pi4hoeHz2qP+iUZQqtUKhobGztryQep9vmcOnVKJ0+ebMnRnOPj49q/f7/WrVunjRs3Tr93hhGBMxoJWQ+Z2a9Iisxsl6TfkfSPYcsC0GmSOT6dvCRDaElvVW9vr4aGhlQul9Xf3z/rkgyhJUO0URSpUqno+PHjbRmyHRwc1MjIiLZu3are3t6Wvz7QyRr5s+utkp6l2mrvn5M0IultIYsC0HmSuVfnasCql0qltGbNGq1du7YtASvh7hoaGtLhw4fbOicuWVAVwNkaWYx0QtK74h8AAAA0YMGQZWbPlPQHknbU39/dXxyuLAAAgO7WSD/3lyR9RNInJK3sw4gAAACapJGQVXb3DwevBAAAYAWZM2SZ2UB88W/N7LckfVW1ye+SJHdf3PkZAAAAziHz9WTdo9p6WMmiJ39Yd5tL2hmqKAAAgG43Z8hy9wslycxy7j5Vf5uZ5WZ/FAAAAKTG1smabeFRFiMFAACYx3xzsjZLeoakvJldqTPDhqslsawvAADAPOabk/WvJN0saaukP9WZkDUi6Y/ClgUAANDd5puT9SlJnzKzl7v7l1tYEwAAQNdbcE4WAQsAAGDxGpn4DgAAgEWaM2SZ2Svify9sXTkAAAArw3w9We+M/2W4EAAAYJHmO7rwlJndIelCM7tt5o3ufkO4sgAAALrbfCHrX0u6StJnVFvCAQAAAA2abwmHoqS7zewadz9hZv3x9rGWVQegY7h7u0uQJE1MTKhYLGrt2rVtrcPMFEWRKpVKW+tIp9PK5/OanJxsex0AztbI/4pN8bDhgCQzsxOSft3dHwxbGoBO4O4qlUptDxPlclmPPfaYjh49KklavXq1du/erd7e1p+AIpPJKIqi6bqmpqZUrVZbWoOZKZvNavXq1XJ3TUxM6OjRoyqVSi2tI4oibdq0SatXr27p6wLdoJGQ9TFJv+/u35UkM7s23nZNwLoAtJm7q1KptHynPVsdR48e1f79++Xu02FmaGhIP/7xj7Vlyxbt3LmzJT0pURQpk8lIqoUcqdaD09/fr0KhoEKhELwGSerp6VEul5uuw8zU19eniy66SKdPn9aJEyda0vO4bt06bdy4UalUaro9AJzRyG+lviRgSZK732VmfQFrAtBm1WpVxWKx7UOEw8PDevjhh+fsKapWqzp69KiOHz+uiy66SFu2bAmys0+lUspkMtOBpl5yPZvNqqenR1NTU8GCaRRFyufzs4aapLaBgQGtXbtWx44d08jISJA6+vr6tGXLFqXTaaVSLLcIzKWRkPWYmf0H1SbAS9LrJD0WriQA7eLuKhaLLR/6mqlQKOiRRx7RqVOnFqylWq2qWq3qkUce0RNPPKGLL75Ya9asaUodZqZ0Oq0oihYMb0nIyefzymazmpycbNoQa/K86XR6wTpSqZRSqZTOP/98nXfeeTpy5IimpqaaUkcmk9GWLVvU29tLuAIa0EjIeoOk/yjpK5Jc0g/ibQBWCHdXuVxWuVxuax2VSkWHDh3SoUOH5O6L6kmrVquamJjQvn37tH79eu3atUvZbHbJtaTT6ekhyMX0jpmZUqmU+vr6VCqVNDU1tawewWw2O/0+FlNHKpVSNpvVjh07NDo6qmPHji059JmZNmzYoIGBgVl78wDMbsGQ5e6Dkn5nsU9sZp+UdL2kp9z92fG290m6UVJV0lOSbnb3I/E8r69Jejx++Ffc/b2LfU0Ai5PMcSoWi22v4+TJk/rpT3+qSqWyrJ60arWqkydP6tSpU9q+fbu2bdu2qF6XVCqlnp4eSYsLNfWSx2UyGWUyGU1NTS26jTOZjHK53LJCTfLY1atXa9WqVXrqqad0+vTpRT3H6tWrp4dh6b0CFsdCzbkwsxdKGpP06bqQtdrdR+LLvyPpUne/JQ5Zf+Du1y/mNfbs2eN79+5tcuXAuaFT5l2Nj4/r4Ycf1tjYWNOHKVOplNLptHbv3q3169fPG1bMTJlMJsgk7qRXbnJycsHewlQqpXw+39AQ5WJVq1VVKhUdOXJE4+Pj8943l8vp/PPPV09PD+EKmMHM7nH3PQvdL9jhOO7+fTPbMWNb/SzMPtWGHwG0UKcsyVAqlfToo4/q+PHjweaAJUHyoYceUn9/vy6++GL19T39uJ36JRlCDIUlPUq9vb2qVCqanJx82ns2M+VyuacdvdhMyXytCy64QJOTkzpy5MjTJulHUaTNmzdr1apVDA0Cy9Ty1ePM7P2Sfk3SsKQX1d30fDO7T9IR1Xq1Hprj8W+S9CZJ2rZtW+BqgZWjU+ZdubuefPJJPfbYY2ctyRBStVrVyMiI9u7dq82bN2vnzp3TwSpkqJkpWcC0v79fxWJxekL6zCUZQkulUurt7dVFF12kwcFBnThxQtVqVQMDA9q4cSPhCmiSBYcLzeyZkj4saZO7P9vMLpd0g7v/pwWfvNaTdXsyXDjjtndKyrn7e8xstaSqu4+Z2XWS/tzddy30/AwXAo1xdxUKhbYPDU5NTWnfvn0qFAptO4IxCRjXXHNNkCG5RiWfhbu3NdRUq9Xp4cykpwvA/BodLmzkf9PHJb1TUkmS3P1+Sa9eXnmSpM9Kenn8nCPJ6Xrc/euSMmZ2XhNeA4DO7EjbbXh4uK0BS6q1xZo1a9reW5O8frvrSKVSiqJIURQRsIAma+R/VK+7/2jGtiWNN5hZfe/UjZIejrdvtvi3jJldHdd1aimvAaCzMQx1tk5pj06pA1hJGpmTddLMLlI8Sd3MbpJ0dKEHmdnnJF0r6TwzOyzpPZKuM7Pdqi3hcFDSLfHdb5L0m2ZWljQp6dXeCX92AwAALFEjIestqp2r8GIze1K1taxeu9CD3P01s2z+yznu+xeS/qKBWgAAALrCvCHLzFKS9rj7v4zPV5hy99HWlAYAANC95p2T5e5VSW+PL48TsAAAABrTyMT3b5nZH5jZBWY2kPwErwwAAKCLNTIn61Xxv2+p2+aSdja/HAAAgJWhkRNEX9iKQgAAAFaSBUOWmWUk/aakF8ab7pL0UXcvzfkgAACAc1wjw4UflpSR9KH4+q/G234jVFEAAADdrpGQ9XPu/py669+JT+QMAACAOTRydGElXvFdkmRmOyVVwpUEAADQ/RrpyfpDSd81s8ckmaTtkl4ftCoAAIAu18jRhd+OT+y8O970U3cvhC0LAACguy04XGhmb5GUd/f73f1+Sb1m9lvhSwMAAOhejczJeqO7DyVX3H1Q0hvDlQQAAND9GglZkZlZcsXMIkk94UoCAADofo1MfP+mpC+Y2Ufj62+OtwEAAGAOjYSsd0h6k2qrvkvSnZI+EawiAE3n7u0uQZJkZh1VS7slbdEJtQBovgWHC9296u4fcfebVAtbP3R31skCukS1WlWpVNLU1JTcvW0hx901MDCg888/X6lUqm3BIooilctlVavVtrx+wt1VKpVULBbb+rkACKeRcxfeJemG+L73SHrKzP7R3X8vcG0AlqFaraparWpoaEiFQm3VlcnJSfX19SmTybQs5CThoVwuy921fft2bd68WY8++qiGh4dbFnaiKFI2m9Xll1+ugYGBlrzmbNxd1WpVk5OTqlRqf68Wi0Xlcjml02l6tYAVpJHhwjXuPmJmvyHp0+7+HjO7P3RhAJYm6RUZGRnR+Pj4024bGxtTOp1WX19f0B6lJFxVKpWnBalsNqtLL71Uw8PD2r9/v4rFYrCwlUqllEqltHv3bm3btq1tISZpj8nJSZVKpbNuq1armpiYUDqdVj6fl5kRtoAVoJGQlTazLZJeKeldgesBsERJuJqcnNTIyMi8oaVcLmt4eFjZbDbITj3prUl6auayZs0aXXXVVTp27JgOHjyoarXa1GGzVCqlrVu3avfu3cpkMk173sVI3k+xWNTU1NS89y2XyxodHVVPT49yuZwk5msB3ayRkPVeSX8n6e/d/cfxuQsfCVsWgMWoVqsql8saGhp6Wi/JfAqFgorFovL5vLLZrKTl7dSToFepVBoOS2amLVu2aMOGDTpw4IBOnDix7F6tKIq0atUqXXbZZVq1atWynms5kraYnJxc1HsqFosqlUrK5XLT4ZCwBXQf6+bJlnv27PG9e/e2uwygbZKen+HhYU1OTi7ruVKplPr6+pY0L2jmvKvlmJiY0P79+zU+Pr7osBVFkaIo0mWXXaaNGze2dWgw6VUsl8vLeq5UKqV8Pq8oighaQIcws3vcfc9C92ukJwtAh0l24mNjYxobG2vKEFu1WtXo6KgymYz6+voaHkJMemuaNaeqt7dXl112mU6fPq1HH320oedOat25c6d27typKIqaUstiJZ/D1NSUisViU56zWq1qfHxcmUxGuVyO+VpAFyFkAV2mWq2qWCxqaGhowTlPS1EqlTQ0NKRcLqd8Pi9p9qGqRuddLYWZaf369Vq3bp0OHz6sJ598cs5lDlKplDZu3KhLL710eh5TqyV11S+V0WylUkmlUknZbLYpQ7sAwmtkCYcL3f3xhbYBCCtZkmFwcLBpvSTzmZqaUqFQUG9vr3p6amfSql9MtBlDgwtJpVLatm2bNm3apMcee0xDQ0PTvVpRFCmXy+nyyy/XunXrgtYxn9mWZAipfh4dSz4Ana2RnqwvS7pqxra/lvTc5pcDYKYkVAwPD2tiYqKlr+3uGh8f19TUlPr7+5VKpZo6NNiobDarSy65RCMjI9q/f79KpZIuueQSbd26te3zrqamphZ1sEGzXntiYkJRFKm3t5chRKBDzRmyzOxiSc+StMbM/m3dTasltadPHjgHFQoFDQ4OtnVF8EqlovHxcWWz2bbuzFevXq2rr75a+XxeqVQj57cPp1QqLftgg+WqVCoqFAptGyYFML/5erJ2S7pe0lpJv1S3fVTSG0MWBQBzodcGQLeYM2S5+9ckfc3Mnu/uP2xhTQAAAF2vkTlZ+83sjyTtqL+/u78hVFEAAADdrpGQ9TVJP5D0LUnhD50BAABYARoJWb3u/o7glQAAAKwgjRyec7uZXRe8EgAAgBWkkZD1u6oFrSkzGzGzUTMbCV0YAABAN1twuNDd23cKewAAgC61YE+W1bzOzP5DfP0CM7s6fGkAAADdq5Hhwg9Jer6kX4mvj0n6n8EqAgAAWAEaObrw5939KjP7iSS5+6CZ9QSuCwAAoKs10pNVMrNIkkuSmW2Q1NqzwwIAAHSZRkLWByV9VdJGM3u/pL+X9J+DVgUAANDlGjm68LNmdo+kl0gySb/s7v8cvDIAAIAu1sicLEk6rtqpddKS8mZ2lbvfG64sAACA7rZgyDKz90m6WdKjiudlxf++OFxZAAAA3a2RnqxXSrrI3YuhiwEAAFgpGpn4/qCktaELATA7d1/4Tmg5PhcAC2mkJ+u/SPqJmT0oqZBsdPcbglUFQJJULBY1Pj6uKIpUrVZVrbZn9ZR0Oq3+/n6ZmSqVSlvryOVybXnthLtPfxaZTKbt7ZHNZtvy2gAW1kjI+pSkP5H0gFgfC2iJSqWi0dFRFYu1UXozUxRFSqVSqlQqLetFSaVS6uvrUzqdlplJkqIoUhRFKpfLLa0jn88riqLpOtqhWq2qWCxOv+/kczlX2wPA/BoJWRPu/sHglQCQu2t8fFwTExOz3p7s1N1dlUolaC35fH6616h+R55cTqfTcneVy+WgdeRyOfX09DytjlZyd5VKpVnb/FxsDwCNaSRk/cDM/ouk23T2cCFLOABN4u4qFAoaHR1dsDfEzKZ/Qgwh9vT0qLe3d/o1Fqojk8moWq02PfRlMhnlcrkF6wgpCU2NBKfQ7dHT0zNr6AXQuRoJWVfG/z6vbhtLOABNUiqVNDIysuidcv0QYjPCVhRF6uvrW/QQlJkplUpND2U2o458Pq9UKtXWMFGpVFQqlRY9BLhS2wPA4jWy4vuLWlEIcK6pVqsaHR1VoVBY+M7zWO58LTNTb2/vsoag6udrLaeOXC6nTCaz5Dqawd1VLBaXFY6a1R75fP6s+XAAuksji5FmJb1c0o76+7v7e8OVBaxc7q6JiQmNj4839XmXMl8rl8spn89PP74ZNdQPZTZaRzabnT5KrhuGBhvVze0BYPkaGS78mqRhSfeobk4WgMVL5l2FOuS/0flamUxGfX19weY7zexdm6uOdDqtfD7f9nlXyVGDoTTaHp0wDw1A8zQSsra6+8uCVwKsYOVyWaOjoyqVSi15vbmGEGdbkqEVdcxc4qBTliCYuSRDaJ3eHgCaq5GQ9Y9mdpm7PxC8GmCFcXeNjY1pcnKyLa9fP4TYrqPTZi5xkE6n274EwXxLMoQ2sz0ymUzb56EBCKORkPUCSTeb2eOqDReaJHf3y4NWBqwAxWKxbQErkSwrkAxDtbOOdDqtTCbT9jBRqVTaErDqdVJ7AAijkZD1i8GrAFYwM+M8dwBwDmokZI02uA0AAACxVAP3uVfSCUk/k/RIfPmAmd1rZs8NWRwAAEC3aiRk3SnpOnc/z93XqzZ8eLuk35L0oZDFAQAAdKtGQtbz3P3vkivufoek57v73ZKywSoDAADoYo3MyTpqZu+Q9Pn4+qskHTezSFKYFRUBAAC6XCM9Wb8iaaukv4l/tsXbIkmvDFcaAABA92rkBNEnJb11jpv3N7ccAACAlWHOkGVm/7+7v83M/lbS0xb5cfcbglYGAADQxebryfpM/O9/b0UhAAAAK8mcIcvd74n//V7rygEAAFgZ5hsufECzDBOKcxcCAAAsaL7hwutbVgUAAMAKM99w4cHZtpvZCyS9RtJbQhUFAADQ7RpZjFRmdqVqa2O9QtLjkr4SsigAAIBuN9+crGeq1mP1GkknJX1Bkrn7i1pUGwAAQNearyfrYUk/kHS9u++XJDP7vZZUBQAA0OXmO63Ov5V0VNJ3zezjZvYS1Y4sBAAAwALmDFnu/jfu/mpJF0v6rqS3SdpoZh82s5e2qkBgqdxnW4Gk9TqlDmA+fE+B5lvwBNHuPu7uf+Xuv6TaiaJ/IukdwSsDlsjdVa1WVSqVVK1W21ZHtVqVmalcLre1DkmqVCqqVqtt35FWKpW2vn4iiiKZtb9jPvlM2vm5JK/f7jqAlaihowsT7j4o6WPxz4LM7JOqrbf1lLs/O972Pkk3SqpKekrSze5+xGq/8f5c0nWSJuLt9y6mPpzbkp3E+Pi4hoaG5O7K5XIaGBhQKpVSKrXg3xRNrePEiRN68sknVa1Wlc/np+to9c49lUopnU5PB610uvbfvtV1RFGkTCbTEeHGzJTNZlWpVFQqldpWh7urUChMt01SW6teW5LK5bLK5bIktaUOYCULvde5VdLLZmz7gLtf7u5XSLpd0rvj7b8oaVf88yZJHw5cG1aQarWqQqGgY8eOaXBwcHoHMjU1pSNHjmh4eLglPTmVSkVjY2N66KGH9MQTT0z3YE1OTurJJ59sWR1SbSfZ09NzVrBxd5VKJVUqlZb1WqRSKWWzWfX09HTUjtvMlE6nlcvlFEVRW2upVCqamppSuVwO3qOUPH/9a7ajDuBcsKierMVy9++b2Y4Z20bqrvbpzKl7bpT0aa/9r77bzNaa2RZ3PxqyRnS3arWqarWq06dPa2pqas77jY6Oanx8XOvWrVM+n5eZNXWHX6lUVKlUdODAAY2MjMx5v5GRkek6crlcsN61TCYz73Mn7RZFUbDeNTNTJpNpe4BZSBJGO2GIuVwuq1KpTH9+zf5ckuBULBbnDVCh6wDOFYsKWWZ2vbvfvtwXNbP3S/o1ScOSknW3niHpibq7HY63HZ3x2Dep1tOlbdu2LbcUdKlkZzE8PKzR0dGGHlOtVnXq1Cn19PRoYGBA6XR62SEnmf915MgRPfXUUw395V+pVHTy5En19PRo/fr100GnGaIoWtR8o1BDiOl0Wul0uqt2zkmPWzKE2K5enCQEpVKp6V7I5bZj8l6SXsx21QGcaxb7m/29zXhRd3+Xu18g6bOSfnuRj/2Yu+9x9z0bNmxoRjnoIkmomZiY0JEjRxoOWPWKxaKOHTum06dPT/foLLWO06dP64EHHtDx48cXvVMuFos6evSoBgcHlz2EmEql1NPTs6RgM3MIcTl1RFGkXC7XMXOvliKKImWz2eng2S7JEHgS+JbyuSSPK5fLmpqaWtKBB0kdDCECi7fY3yLN/q35WUlfl/QeSU9KuqDutq3xNkBS7Zd9uVzW6dOnVSwWl/18ExMTmpyc1Jo1a9Tf39/wX+qVSkWFQkEHDhzQxMTEsusYHx/XxMSE1q5dq76+vkX1GCTziprRE7acIcRkyK1VBxeElgx1ptPpRfX+hJAMRdcPvTby2dQfZduMYJRMkF9sHcC5bLEh683LfUEz2+Xuj8RXb1RtZXlJuk3Sb5vZ5yX9vKRh5mNB0nRP0+DgoMbHx5v63O6uoaEhjY2Nad26dcpms3MGhSSEHDp0SKdPn256HYODgxodHdXAwEBDgSUJV83e0SU79aRXbKHnT2pdiTvc+vlaC81jCq1UKqlcLqunp2e6ttkkvU2h5pc1WgeAxS/h8KPF3N/MPifpWknnmdlh1XqsrjOz3aot4XBQ0i3x3b+u2vIN+1VbwuH1i3ktrDzJzmJsbEzDw8NBd3DlclknTpyYdcmHpI7jx4/r6NGjQSdGl8tlPfXUU9N1zDa/KlmSIfTOrVwuT/eUSU/fmXbjvKulSuZrJWGrXeZb8mEp865C1AHgDOvm8fU9e/b43r17210Gmqz+CKjTp0+fdYh5q6xatUpr1qyRVDsy8eDBg23ZuSZ1JKGvXaEmlUpNB75kIvRKGRpcrGSOUzu+lzMlQVdSW9f8qq+DsIVzgZnd4+57Frpfe2d2ArMoFAoaHh5WoVBoWw2jo6MaHBxUoVBo+hDlYutwd23cuLGtSyEkc7V6eno6fkmG0JL5WpLaHrSSpRak9p4WJ5kUn7QLgJo5/xQ1s+1mtqbu+ovM7M/N7PfNrKc15eFc1c4hmUSlUmnKxPblcveO6R04V3uvZtMpnwlH/AGda77fmF9UbbFQmdkVkr4k6ZCk50j6UPjSAAAAutd8w4V5dz8SX36dpE+6+5+aWUrSvvClAQAAdK/5erLq+8JfLOnbkuTu7TvnBAAAQJeYryfrO2b2RdVOa7NO0nckycy2SGr/hBkAAIAONl/IepukV0naIukF7p4cG7xZ0rtCFwYAANDN5gtZF7j752dudPefmNkvBKwJAACg6803J+suM3u7mU0vimNmm8zsf0n6s/ClAQAAdK/5QtZzJV0kaZ+ZvdjMflfSjyT9UNLVrSgOAACgW805XOjug5LeHIerb0k6Iul57n64VcUBAAB0q/lWfF9rZh9V7UTNngDFDgAAFaRJREFUL5P015K+YWYvblVxAAAA3Wq+ie/3qray+1vcvSzpjnjl9w+Z2UF3f01LKgQAAOhC84WsF84cGnT3fZKuMbM3hi0LAACgu805XDjf3Ct3/3iYcgAAAFaG+Y4uBAAAwBIRsgAAAAIgZAEAAARAyAIAAAiAkIWOk06ntWnTJmWz2bbVYGbauHGjLr74YvX397etjiiKtHnzZq1atUpRFC38gEDMTJlMpm2v34lSqZR6enpkZm2tI5PJtL0OM2vr9xPoVPMt4QC0RRRFSqfT2rBhgwqFgk6fPq1KpdKy1+/v79fatWtlZjIz7dq1SyMjIzp06JBKpVLL6tiwYYO2bt06XUdfX59KpZKmpqbk7i2ro6enR7lcrmWv1y2SzyWbzapSqbT0uyHV/hhJp8/8Cu+UOgCcwf8MdJzkL/JUKqVcLqctW7ZodHRUIyMjQcNFNpvVwMCAoihSKnWmkzeKIq1du1Zr1qzR0aNHdezYsaB1rFq1Sjt27FA6nX5a70Amk1Emk1GhUFChUAhWg1TbeeZyOaVSqbb31nSi+jaJokhRFKlUKgX/gyDpQZtZQyfVAaCGkIWOlvQWrFq1Sv39/RocHNTExERTXyOKIg0MDCibzZ4VrmarY8uWLdq4caMOHjyooaGhptbR09Oj7du3q7+/f86hl2Rnls1m1dPTo8nJSZXL5abWkYTbdDrNzrNBSTtlMhml02mVSiVVq9Wmv0Ymk5k39HZKHQBqCFnoCkn4GRgY0OrVq3Xq1KllD4uYmVavXq1Vq1ZNh6hG6kilUrrwwgs1NTWlAwcOaHJycll1pFIpbdmyRZs2bWq4juR+vb29qlQqmpycbMrONJvNTs+FYwe6eMnn0tPTo2q1qlKp1JRez0wmMx28F/P9aHcdwLmOkIWukkqllMlktGnTJk1MTGhoaGhJ4SKfz2tgYEBmNmfv1XyiKFJvb68uueQSnTp1SocPH17S8MzAwIC2bds2Hd4WK5lw3N/fr2KxqKmpqUU/h1TbeeZyuYZDHuaXfK+y2azK5fKSexujKJo+4GApn0un1AGcqwhZ6Dr1E8F7e3s1NDSksbGxhh6byWS0fv16pdPpJYWa2epYv369BgYGdPjwYZ04caKhx/b29mrHjh3KZrPLPior2en19PRMDyE22ssXRZHy+TxDPwEk7ZlMDF/MPKnkj4lmhN5OqQM4FxGy0LWSX/xr166dHkKcazJ4KpXS2rVr1dvb2/QdRhLWtm7dqs2bN+vAgQMaHR2d9b7pdFoXXHCB1q1b1/Q6kufK5/PKZrOanJycc2dqZsrlcvROtMDMeVLFYnHOobuQ8506pQ7gXELIQtdLhtrmWvJh5pIMoSRHdu3atUujo6M6ePCgisWipDPrbp1//vlLHqJsVDKE2NfXp3K5rMnJybN2pvVLMrADbZ36gxZmW2qhfimEkJ9Lp9QBnAsIWVgx6pd8GBsb09TUlAYGBpY832k5daxatUrPfvazdfz4cY2NjWn79u1NGaJcDDNTOp3WqlWrVCgUVKlUlM/nGfppo6Tdk0BeLpfl7i3vUZxZRxK06NkEmouQhRWlfsmH5KjBdkjC1KZNm7R58+a21VHfa1F/He1VP0+q/nq76iBcAWEQsrAidcrOopU9V/PplPbA2Trlc+mUOoCVpjP2AAAAACsMIQsAACAAQhYAAEAAhCwAAIAACFkAAAABELIAAAACIGQBAAAEQMgCAAAIgJAFAAAQACELAAAgAEIWAABAAIQsAACAAAhZAAAAARCyAAAAAiBkAQAABEDIAgAACICQBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAEQsgAAAAIgZAEAAARAyAIAAAiAkAUAABAAIQsAACAAQhYAAEAAhCwAAIAACFkAAAABELIAAAACIGQBAAAEQMgCAAAIgJAFAAAQACELAAAgAEIWAABAAIQsAACAAAhZAAAAARCyAAAAAiBkAQAABEDIAgAACICQBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAEEC1lm9kkze8rMHqzb9gEze9jM7jezr5rZ2nj7DjObNLN98c9HQtUFAADQCiF7sm6V9LIZ2+6U9Gx3v1zSzyS9s+62R939ivjnloB1AQAABBcsZLn79yWdnrHtDncvx1fvlrQ11OsDAAC0UzvnZL1B0jfqrl9oZj8xs++Z2S/M9SAze5OZ7TWzvSdOnAhfJQAAwBK0JWSZ2bsklSV9Nt50VNI2d79S0u9L+iszWz3bY939Y+6+x933bNiwoTUFAwAALFLLQ5aZ3SzpekmvdXeXJHcvuPup+PI9kh6V9MxW1wYAANAsLQ1ZZvYySW+XdIO7T9Rt32BmUXx5p6Rdkh5rZW0AAADNlA71xGb2OUnXSjrPzA5Leo9qRxNmJd1pZpJ0d3wk4QslvdfMSpKqkm5x99OzPjEAAEAXCBay3P01s2z+yznu+2VJXw5VCwAAQKux4jsAAEAAhCwAAIAACFkAAAABELIAAAACIGQBAAAEQMgCAAAIgJAFAAAQACELAAAgAEIWAABAAIQsAACAAAhZAAAAARCyAAAAAiBkAQAABEDIAgAACICQBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAEQsgAAAAIgZAEAAARAyAIAAAiAkAUAABAAIQsAACAAQhYAAEAAhCwAAIAACFkAAAABELIAAAACIGQBAAAEQMgCAAAIgJAFAAAQACELAAAgAEIWAABAAIQsAACAAAhZAAAAARCyAAAAAiBkAQAABEDIAgAACICQBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAEQsgAAAAIgZAEAAARAyAIAAAiAkAUAABAAIQsAACAAQhYAAEAAhCwAAIAACFkAAAABELIAAAACIGQBAAAEQMgCAAAIgJAFAAAQACELAAAgAEIWAABAAIQsAACAAAhZAAAAARCyAAAAAiBkAQAABEDIAgAACICQBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAEQsgAAAAIgZAEAAARAyAIAAAiAkAUAABAAIQsAACAAQhYAAEAAhCwAAIAACFkAAAABELIAAAACCBayzOyTZvaUmT1Yt+0DZvawmd1vZl81s7V1t73TzPab2U/N7F+FqgsAAKAVQvZk3SrpZTO23Snp2e5+uaSfSXqnJJnZpZJeLelZ8WM+ZGZRwNoAAACCChay3P37kk7P2HaHu5fjq3dL2hpfvlHS59294O6PS9ov6epQtQEAAISWbuNrv0HSF+LLz1AtdCUOx9uexszeJOlN8dVC/XAkmuI8SSfbXcQKQ5s2H23afLRp89Gmzdcpbbq9kTu1JWSZ2bsklSV9drGPdfePSfpY/Dx73X1Pk8s7p9GmzUebNh9t2ny0afPRps3XbW3a8pBlZjdLul7SS9zd481PSrqg7m5b420AAABdqaVLOJjZyyS9XdIN7j5Rd9Ntkl5tZlkzu1DSLkk/amVtAAAAzRSsJ8vMPifpWknnmdlhSe9R7WjCrKQ7zUyS7nb3W9z9ITP7oqR/Um0Y8S3uXmngZT4WpPhzG23afLRp89GmzUebNh9t2nxd1aZ2ZsQOAAAAzcKK7wAAAAEQsgAAAALo6JDFqXmabzFtamY7zGzSzPbFPx9pX+Wda442fV/cnvvM7A4zOz/ebmb2wfh7er+ZXdW+yjvXItv0WjMbrvuevrt9lXeu2dq07rb/18zczM6Lr/M9bcAi25TvaQPm+L//x2b2ZF3bXVd3W0fv9zs6ZIlT84Rwqxps09ij7n5F/HNLi2rsNrfq6W36AXe/3N2vkHS7pOQX6i+qdvTsLtUW1f1wq4rsMreq8TaVpB/UfU/f26oiu8ytenqbyswukPRSSYfqNvM9bcytarxNJb6njbhVs7SppD+ra7uvS92x3+/okMWpeZpvkW2KBszRpiN1V/skJUeY3Cjp015zt6S1ZralNZV2j0W2KRowW5vG/ky1pXXq25PvaQMW2aZowDxtOpuO3+93dMhqwBskfSO+/AxJT9TdNuepeTCv+jaVpAvN7Cdm9j0z+4V2FdWNzOz9ZvaEpNfqTK8L39NlmKNNJen5ZnafmX3DzJ7VpvK6jpndKOlJd79vxk18T5donjaV+J4ux2/HQ9efNLN18baO/552bciyZZyaB7ObpU2PStrm7ldK+n1Jf2Vmq9tVX7dx93e5+wWqtedvt7uelWCONr1X0nZ3f46k/yHpb9pVXzcxs15Jf6SzwyqWYYE25Xu6dB+WdJGkK1TbL/1pe8tpXFeGLDtzap7Xcmqe5pitTeMu2FPx5XskPSrpmW0rsnt9VtLL48t8T5tjuk3dfcTdx+LLX5eUSSYbY14XSbpQ0n1mdkC17+K9ZrZZfE+Xas425Xu6dO5+3N0r7l6V9HGdGRLs+O9p14Us49Q8TTdXm5rZhmQSoZntVK1NH2tPld3FzHbVXb1R0sPx5dsk/Vp89NbzJA27+9GWF9iF5mpTM9tsVjuFhJldrdrvtVOtr7C7uPsD7r7R3Xe4+w7Vhlqucvdj4nu6JPO1Kd/TpZsxH/DfSEqOPOz4/X7LTxC9GNaaU/OcUxbTppJeKOm9ZlaSVJV0i7s3OiHxnDFHm15nZrtVa7eDkpIjM78u6TrVJmhOSHp9ywvuAots05sk/aaZlSVNSnp1XQ83YrO1qbv/5Rx353vagEW2Kd/TBszxf/9aM7tCtQMJDkh6syR1w36f0+oAAAAE0HXDhQAAAN2AkAUAABAAIQsAACAAQhYAAEAAhCwAAIAACFlAi5jZL5uZm9nFddtSZvZBM3vQzB4wsx/H670047l3JGeyN7M9ZvbB5ryT5onrvrTddTQiXovnW2a2z8xeNeO2m83s/LrrB1qx0KSZfaJb2m8h9d9XYKUgZAGt8xpJfx//m3iVpPMlXe7ul6m20N5Qk557mrvvdfffWcLzhvbLkloSEpKFdZfhSkly9yvc/QszbrtZtc+xZcwscvffcPd/avLzLnr9xHjRUvYnwAz8pwBawMz6Jb1A0r+T9Oq6m7ZIOhqfLkLuftjdB5v03PX3udbMbo8vbzCzO83sobgn5KCZnRf3JPyzmX08vu0OM8vHj7nLzP7MzPbG9/k5M/uKmT1iZv+p7nVeZ2Y/int7Plp3xoAxq53c+T4zu9vMNpnZNZJukPSB+P4Xzaj5FXEP331m9v14W2Rm/z3efr+ZvTXe/hKrncj8AaudQDYbbz9gZn9iZvdKeoWZXWRm3zSze8zsB/U9f3WvO2BmfxM//91mdrmZbZT0vyT93MxazewmSXskfTa+LR/f9FYzuzeu6eL4vn1xfT+K671xjs/q+2b2v83sp2b2kSTAxO34p2Z2n2onG77LzPbU3faB+LP7lpldHd/+mJndEN9nR/y+741/rql7zR+Y2W2S/snM3mtmb6ur6f1m9rsz6twR1/dp1VbgvsDMXmpmP4yf+0vxd1Nm9m6r9dI+aGYfM5te+fy58ed7n6S3zPbdBbqau/PDDz+BfyS9VtJfxpf/UdJz48tbVVvBeJ9qJz29sonPvUPSg/HlayXdHl/+C0nvjC+/TLVVlM+L71+WdEV82xclvS6+fJekP4kv/66kI6oFxKxqpw5ZL+kSSX8rKRPf70OSfi2+7JJ+Kb783yT9+/jyrZJumuN9PSDpGfHltfG/vynpryWl4+sDknKSnpD0zHjbpyW9Lb58QNLb657z25J2xZd/XtJ3Znnd/6Hayt2S9GJJ+2a24SyPuUvSnrrrByS9Nb78W5I+EV/+z3VtulbSzyT1zXiuayVNSdopKZJ0Z9JGcTu+crbXjW/7xfjyVyXdISkj6Tl176FXUi6+vEvS3rrXHJd0Yd135974ckq185aun1HnDtVW339efP08Sd9P3o+kd0h6d/I51T3uM3XfhfslvTC+/AHF31d++FkpP/RkAa3xGkmfjy9/Pr4udz8sabdqpzaqSvq2mb2kGc89jxck93f3b0qq7zl73N33xZfvUW1Hmrgt/vcBSQ+5+1F3L6h2PssLJL1E0nMl/djM9sXXd8aPKUq6fY7nncs/SLrVzN6oWtiQpH8p6aPuXo7rP61a+z3u7j+L7/Mp1U4JlfiCNN3jd42kL8X1fVS1oDjTC1QLAnL370hab2arG6h3pq/E/9a/35dK+v/i179LtYC4bZbH/sjdH/PaKUI+F9ckSRVJX57j9YqSvhlffkDS99y9FF9OXj8j6eNm9oCkL+nsodofufvjkuTuBySdMrMr45p/4vHJ4mc46O53x5efFz/fP8Tv79clbY9ve5GZ/Z/4dV8s6Vlmtla18Pz9+D6fmeN9AV2ro89dCKwEZjag2o7lMjNz1QKDm9kfek1B0jckfcPMjqs2T+nbdY//edUCgVTrGbitkedeYrmFussVSflZbqvOuF9Vtd8lJulT7v7OWZ635O7JObwqauB3j7vfEr/3fy3pHjN7bmNv4WnG439Tkobc/YolPs9iJW1U/35N0svd/acLPHbm+c6S61M+97nZ6tt4+jNy96qdmWf1e5KOq9a7lVKtxywxrrN9QrW5ZpslfXKO16x/jEm6093PCvlmllOtV3OPuz9hZn+sWrgEVjx6soDwbpL0GXff7u473P0CSY9L+gUzu8rio9LieTeXq3by42nu/n+8Ntn6ivqAtdBzz1PPP0h6ZfyaL5W0rhlvUrVgeFM8fymZ27R9gceMSlo12w1mdlH83t8t6YRqvWV3SnpzEhrikPlTSTvM7P+JH/qrkr438/ncfUTS42b2ivixZmbPmeWlf6DaEKzM7FpJJ+PHLul9zPB3qs3VSuYkXTnH/a42swvj78SrVDuooRnW6MwcwF/VmR7C2XxVteHkn4vrXsjdkv5F8jnE88+eqTOB6mTcm3iTJLn7kKQhM0t66V672DcDdDpCFhDea1TbYdX7crx9o6S/tdqh6/erNifqL5r03HP5j5JeGr/mKyQdUy0kLIvXjnL795LuMLP7VQtEsw3H1fu8pD+MJ4FfNOO2D8STxh9Uba7Zfar1rhySdH88WfpX3H1K0utVGwZ8QLVenI/M8XqvlfTv4sc+JOlpE88l/bGk58bv4b+qNuy1kFslfWTGxPfZvE+1Ibv7zeyh+Ppsfqza9+CfVQvNMz/jpfqQpF+P3//Fenrv1TR3L0r6rqQvztN7Vn//E6r1fH0ubrsfSro4DlMfV21y/N+p9t4Sr5f0P+PhRVvSOwI6mJ3pXQZwLrDakXcVdy+b2fMlfbiFQ2hYQNx79gfufn2b60hJulfSK9z9kXbWAnQr5mQB555tkr4Y70SLkt7Y5nrQYay2wOntkr5KwAKWjp4sAACAAJiTBQAAEAAhCwAAIABCFgAAQACELAAAgAAIWQAAAAH8XwOXskwM8xxgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(1, 1, 1, aspect=1)\n",
    "ax.set_xlim(120,152)\n",
    "ax.set_ylim(120,152)\n",
    "ax.set_ylabel(\"XS - Alignment score of the next best alignment\")\n",
    "ax.set_xlabel(\"AS - Alignment score of the primary read\")\n",
    "ax.set_title(\"AS vs XS scores for HLA reads with secondary mapping to a GRCh38 HLA alt contig\")\n",
    "hexbin(x=[x[0] for x in primary_AS_XS],y=[y[1] for y in primary_AS_XS], gridsize=75, bins='log', cmap='Greys')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a simple example, but by combining the alignbatch approach above more complex examples that look at relationships between different alignments of a read pair and/or factors like the presence and type of variants can easily be constructed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outputing reads as fastq\n",
    "Sometimes it is useful to output a subset of reads for further analysis, and occasionally you will want to do this in fastq format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@A00132:55:HCKWTDSXX:1:1102:27263:27759\n",
      "AGAAAGAAAGAAAAAGGAAGGAAGGAAGCACCGAAGGAAGGCACCCACGAAGGAAGGCAGGAAAAGCAGACATACAGGCTGAAAGTGAAGGAATGAAAACAGCTATTCCATGCAAATGGTAACCAAAATACAGCAAGAGTGATTACACTT\n",
      "+\n",
      "?????+?+?++?++?+?????+?+++55+++%+??++?5+5++++++%+??++?5+?+?++???5?+?++????????????????????+??????????????+?????????+????????????+???????????????????5?\n",
      "\n"
     ]
    }
   ],
   "source": [
    "test = b'\\xcd\\x01\\x00\\x00\\xf5\\x0c\\x00\\x00^\\x08\\x00\\x00\\'\\tI\\x12\\x02\\x00s\\x08\\x96\\x00\\x00\\x00\\x0b\\x00\\x00\\x00h\\x88\\xeb\\x04\\xd0\\x01\\x00\\x00A00132:55:HCKWTDSXX:1:1102:27263:27759\\x00\\xf0\\x01\\x00\\x00t\\x07\\x00\\x00\\x14\\x11\\x14\\x11\\x14\\x11\\x11\\x14A\\x14A\\x14A\\x14!\"A\\x14A\\x14B\\x12\"\\x12A\\x14A\\x14B\\x14A\\x11\\x14!A!\\x81!D(A\\x11HA\\x14A\\x18A\\x11\\x12\\x14(\\x18\\x82!\\x84!\\x11\\x84H\\x11\"\\x11\\x11\\x81!B\\x11AHA\\x88\\x12\\x12\\x88\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\n\\x1e\\n\\n\\x1e\\n\\n\\x1e\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\n\\n\\n\\x14\\x14\\n\\n\\n\\x04\\n\\x1e\\x1e\\n\\n\\x1e\\x14\\n\\x14\\n\\n\\n\\n\\n\\n\\x04\\n\\x1e\\x1e\\n\\n\\x1e\\x14\\n\\x1e\\n\\x1e\\n\\n\\x1e\\x1e\\x1e\\x14\\x1e\\n\\x1e\\n\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\n\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x1e\\x14\\x1eSAZchr12,82544431,+,102M48S,60,3;\\x00XAZchr12,+82544561,122S28M,0;\\x00MCZ150M\\x00PGZMarkDuplicates\\x00ASc\\x1eXSc\\x1cMDZ0G30\\x00NMI\\x01\\x00\\x00\\x00RGZNA12778_CTCACCAA-CTAGGCAA_HCKWTDSXX_L001\\x00'\n",
    "\n",
    "def align_to_fastq(align):\n",
    "    readname_length = bam.get_len_read_name(align)\n",
    "    number_cigar_operations = bam.get_number_cigar_operations(align)\n",
    "    sequence_length = bam.get_len_sequence(align)\n",
    "    readname = bam.get_read_name(align, readname_length)\n",
    "    sequence = bam.decode_sequence(bam.get_raw_sequence(align, readname_length, number_cigar_operations, sequence_length))\n",
    "    qual = bam.decode_base_qual(bam.get_raw_base_qual(align, readname_length, number_cigar_operations, sequence_length))\n",
    "    return f\"@{readname}\\n{sequence}\\n{'+'}\\n{qual}\\n\"\n",
    "\n",
    "\n",
    "print(align_to_fastq(test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastqf = open('NA12778_any_HLA_hits.R1.fastq.gz','w')\n",
    "fastqr = open('NA12778_any_HLA_hits.R2.fastq.gz','w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "NA12778_any_HLA_hits = bam.FileReader(gzip.open('Genomes/NA12778_any_HLA_hits.bam'))\n",
    "for align in NA12778_any_HLA_hits:\n",
    "    if not bam.is_flag(align, bam.FLAGS['secondary']):\n",
    "        if bam.is_flag(align, bam.FLAGS['forward']):\n",
    "            fastqf.write(align_to_fastq(align))\n",
    "        else:\n",
    "            fastqr.write(align_to_fastq(align))\n",
    "fastqf.close()\n",
    "fastqr.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### your uses\n",
    "If you have an interesting use of pylazybam feel free to send a pull request against the readme and I will add it to a list of contributed example scripts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}