000_wasabi_plots_initial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_02_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam'"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import glob\n",
    "bam_filenames = glob.glob('/home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/*.bam')\n",
    "\n",
    "sample_id = 'M2_02'\n",
    "\n",
    "bam_filename = [x for x in bam_filenames if sample_id in x][0]\n",
    "bam_filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# figwidth = 4\n",
    "# figheight = len(bam_filenames) * 1\n",
    "# figsize = figwidth, figheight\n",
    "# print figsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pysam\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(context='talk')\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "HTSeq._HTSeq.SAM_Alignment"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(read)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.8.3'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pysam.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from outrigger.region import Region\n",
    "\n",
    "snap25_mxe_event = 'exon:chr20:10265372-10265420:+@exon:chr20:10273530-10273647:+@exon:chr20:10273809-10273926:+@exon:chr20:10277573-10277698:+'\n",
    "\n",
    "snap25_exons = snap25_mxe_event.split('@')\n",
    "snap25_regions = map(Region, snap25_exons)\n",
    "snap25_regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import HTSeq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bam_htseq = HTSeq.BAM_Reader(bam_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bam_htseq?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!samtools index $bam_filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# exon1 = snap25_regions[0]\n",
    "\n",
    "# list(bam_htseq[HTSeq.GenomicInterval(snap25_regions[0].chrom, snap25_regions[0].start, \n",
    "#                                      snap25_regions[-1].stop, snap25_regions[0].strand)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'D': 'deleted',\n",
       " 'H': 'hard-clipped',\n",
       " 'I': 'inserted',\n",
       " 'M': 'matched',\n",
       " 'N': 'skipped',\n",
       " 'P': 'padded',\n",
       " 'S': 'soft-clipped'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTSeq.cigar_operation_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "INSERTION_DELETIONS = ('I', 'D')\n",
    "\n",
    "\n",
    "def skip_bad_cigar(read, bam_filename, bad_cigar=INSERTION_DELETIONS):\n",
    "    # Skip reads with no CIGAR string\n",
    "    if read.cigar is None:\n",
    "        print(\"Skipping read with no CIGAR string: {read_name} (from {bam})\".format(\n",
    "                read_name=read.read.name, bam=bam_filename))\n",
    "        return\n",
    "\n",
    "    # Check if the read contains an insertion (I)\n",
    "    # or deletion (D) -- if so, skip it\n",
    "    for cigar_operation in read.cigar:\n",
    "        cigar = cigar_operation.type\n",
    "        if cigar in bad_cigar:\n",
    "            print(\"Skipping read with CIGAR string {abbrev} (a base in the read was {full}): {read_name} (from {bam})\".format(\n",
    "                read_name=read.read.name, bam=bam_filename, abbrev=cigar, \n",
    "                    full=HTSeq.cigar_operation_names[cigar]))\n",
    "            return\n",
    "    return read\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start: 10265372\n",
      "Skipping read with CIGAR string D (a base in the read was deleted): HWI-ST1001:211:C1PUCACXX:8:2204:7788:53829 (from /home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_02_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([80, 80, 80, ..., 89, 87,  5])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "def count_region_reads(bam_htseq, bam_filename, chrom, start, stop, strand,\n",
    "                       allowed_cigar=('M'), bad_cigar=('I', 'D'), offset=True):\n",
    "    \"\"\"Get the number of reads that matched to the reference sequence\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    bam_htseq : HTSeq.BAM_Reader\n",
    "        Bam file object to get reads from\n",
    "    bam_filename : str\n",
    "        Name of the bam filename for logging purposes\n",
    "    chrom : str\n",
    "        Name of the reference chromosome\n",
    "    start, stop : int\n",
    "        Genome-based locations of the start and stop regions\n",
    "    pad : int, optional\n",
    "        Add a few nucleotides to the left and right of the array for\n",
    "        visually pleasing padding. (default=10)\n",
    "    allowed_cigar : tuple of str, optional\n",
    "        Which CIGAR string flags are allowed. (default=('M') aka match)\n",
    "    bad_cigar : tuple of str, optional\n",
    "        Which CIGAR string flags are not allowed. (default=('I', 'D') aka \n",
    "        insertion and deletion)\n",
    "    offset : bool, optional\n",
    "        If True, offset the region counts so that the array starts at zero \n",
    "        and not the whole chromosome. Useful for plotting just the one region.\n",
    "        (default=True)\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    counts : numpy.array\n",
    "        Number of reads that matched to the genome at every location\n",
    "    \"\"\"\n",
    "\n",
    "    # Add a small amount of nucleotides to left and right for making the plot look nice\n",
    "#     start = start - pad\n",
    "#     stop = stop + pad\n",
    "    \n",
    "    length = stop - start + 1 if offset else stop\n",
    "\n",
    "    counts = np.zeros(shape=(length), dtype=int)\n",
    "    \n",
    "    region_reads = bam_htseq[HTSeq.GenomicInterval(chrom, start, stop, strand)]\n",
    "\n",
    "    print('start:', start)\n",
    "    \n",
    "    for read in region_reads:\n",
    "        read = skip_bad_cigar(read, bam_filename, bad_cigar)\n",
    "        if read is None:\n",
    "            continue\n",
    "        for cigar_operation in read.cigar:\n",
    "            # Only count where the read matched to the genome\n",
    "            if cigar_operation.type not in allowed_cigar:\n",
    "                continue\n",
    "            match_start = cigar_operation.ref_iv.start\n",
    "            match_stop = cigar_operation.ref_iv.end\n",
    "            \n",
    "            if offset:\n",
    "                match_start = match_start - start\n",
    "                match_stop = match_stop - start\n",
    "                \n",
    "            if match_stop < 0:\n",
    "                # If the match_stop is negative, that means we have the other read of the paired end read\n",
    "                # that mapped to somewhere else in the genome\n",
    "                continue\n",
    "            match_start = max(match_start, 0)\n",
    "            match_stop = min(match_stop, length)\n",
    "\n",
    "            counts[match_start:match_stop] += 1\n",
    "    return counts\n",
    "\n",
    "\n",
    "chrom = snap25_regions[0].chrom\n",
    "start = snap25_regions[0].start\n",
    "stop = snap25_regions[-1].stop\n",
    "strand = snap25_regions[0].strand\n",
    "\n",
    "snap25_counts = count_region_reads(bam_htseq, bam_filename, chrom, start, stop, strand)\n",
    "snap25_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 12327)]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAABhCAYAAAAgEOx6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE3xJREFUeJzt3XuQXGWZx/Fvz+QCmWgmI6BUuNeSh1stSJXcClCwKMWt\nGBUVsSxYcMFVoLRAdtEqgiKbWFwKAVcQSSgQBGstI4QCkVqB2gtKBARBeEBRTJaQDDOZyfR0Zrqn\nu/ePc4Y0nemenjndPeed/D5VqZo+l/SZ/HLmPHP6Pe+TKZfLiIiIiIjIzjpm+gBERERERNJKxbKI\niIiISA0qlkVEREREalCxLCIiIiJSg4plEREREZEaVCyLiIiIiNQwZyobm9kxwFp3XxK/XgJ8HzgJ\nyAM/Ay5190K8fhXwRaATuAu4xN01V52IiIiIBKHhO8tmdh7wCDC3YvHdwAZgb+Ao4APAFfH2FwGn\nA0cAhwInApc25ahFRERERNqgoWLZzL4JXAxcXbFsLpAFrnb3grtvAe4BTog3+QLwPXffEq9bBZzb\nzIMXEREREWmlRodhrHb3lWb2wfEF8VCLZVXbLQN+H399CPDHinUOLJ3ugYqIiIiItFtDd5bdffNk\n25jZTYAR3UEG6AJyFZvkgA4zmzfVgxQRERERmQlTesBvIma2G9HY5cOBk929L16VA3av2HQBMObu\n+aTvKSIiIiLSDomKZTNbDPwS2AYc5+6DFatfIrrTvD5+fUi8rCHlcrmcyWSSHJ6IiIiISCNqFp1J\n7yyvBTYBZ7h7sWrd3cBlZvYYMAZcTjR9XEMymQwDA8OUSpppLgQdHRm6u7uUWSCUV3iUWXiUWViU\nV3iamVlPz8Ka66ZdLJvZ8UTzK48AA2Y2fpTPuPuHgB8AewFPAfOAHwM3TOU9SqUyxaL+w4ZEmYVF\neYVHmYVHmYVFeYWn1ZlNtYPfdqAA4O5PAu8hmnu5CPQDX4sLZdy9RDQn89z4fTSmQkRERESCkrQp\nye3AELAn8BngmrjLn5qSiIiIiEjqjBaKjOarRw/XlqQpSRewHFgRNyVZD/wEODveRE1JRERERCRV\nhobzDA6PNjzOudE7y6vd/f3A7yqWLQXy7v56xTInmvUC1JRERERERFKmd3CEbcMF+odGGto+SVOS\nBURjmCvl4uWgpiQiIiIikjLD2ws8+6deerduZ6xYmnT7JFPH5YDdqpYtALIV6xM1Jeno0DOBoRjP\nSpmFQXmFR5mFR5mFRXmFZ7qZvdE3zH8/v4nBbJ5F75rPvnvVnjYOkhXLrwLzzGwfd98YLzN2DL1I\n1JQEoLu7K8HhyUxQZmFRXuFRZuFRZmFRXuGZamajY2WGcgVe/Gs/Hz3hwLpzLEOCYtnds2Z2P7DK\nzC4gmvXiLKIZMCBhUxJAE4MHRJO5h0V5hUeZhUeZhUV5hWe6mW0fiQY5bMvm+eWTf2GPhXPZd8ni\nmtsn7eB3AXArsJFoCrmvu/v4Q4BqSrILUmZhUV7hUWbhaSSzsWKJOZ1TbX0graBzLDxTzWy8sC4D\nT3svpx+zH/suqb39lIpld3+CqAAef70VOLPGtiVgRfxHREREanj59a0ccdB7ZvowRHYJlTehR/NF\nnnjuDY46bO+a2ye9s4yZnQDcSDQt3BvAVe5+r5l1A2uAU4GBePmapO8nIiIym5TLZd7sz6lYFmmT\ncvmdd6H/9w9v8tWzam+f6DMfM+sA1gIr3X0RcD5wp5ntR53ufiIiIhIZyObZ3B/NxDpWLDGSH6NU\n1jAAkWYbL5Krh2yMFup380s6QKob2IMdLbDLwChQon53PxEREQGGcnk2b80xWihSLJXJjYzxwmt9\nM31YIrPOa5u2AfDKxoF3LC+Wyiy79P75tfZLVCy7ez9wC3CfmRWAJ4CLiAroet39GqbfrkVEZLba\nvDXH628OsaE3S/+2EUqlMqOFIg8++frkO4vIlGzYnGWsWKJ/2+hEq2sOTU46DCND1HzkDKIGJB8n\nGr/8bup395vU+JOKm/tzDA433MdEREQkGIPZPL986m8MZvMMj4yRLxQZyhXYsCU7+c4iMiWDw3n+\n+ubQlPdL+oDfp4Bj3P1f4tcPmdmDwLeo391vUttHx+joyPDW4Aj77z6Xzk511EkzdT4Ki/IKjzIL\nTyOZlcplNvXlAOgd2E5hrMgbfTny+aKue22mcyw8U83slQ0DPPtqb63VNf+SpMXyfkD1GI8x4Gng\nxDrd/Sb1p40DHHnwnnTO3Up39wJ6Fu0++U4y49T5KCzKKzzKLDz1Mtt98/DbXw/nSyykg18//X+Q\nYdKuYtIaOsfC02hmHZ0dbOwdnnzDKkmL5UeBlWZ2jrvfaWYfBD5BNF3cAezc3e9jjf7F2VyBgYFh\n/vjnXpb07EamWP9JRZlZ6nwUFuUVHmUWnkYy2/LWjo+EB7Zt57WNW9myNUe5DP39GorRTjrHwjPV\nzAqF4rSyTVQsu/sLZvZp4GozuxH4G3C2uz8TF8nV3f3WN/p3l0plSqUyW4dGKZV2nuZD0kmdj8Ki\nvMKjzMJTL7N8ofT210PDebK5AsXSxNNbSXvoHAtPo5lVz69cpWXDMAB+T1QQH0Q0ldz+48dENKVc\nZ/y6c+ddaxufBWMgm5/smxMREQneM6/0MmeOWl6LpE0zzspfAC8Ci4GPAFea2XEkbEpSKpf5w5/7\n9ESwiIjsEoZHxhjMavYnkbRJdGfZzI4F9ga+4e5l4CUzO56oMcly4O/cvQCsN7PxpiRPNfJ3l0pl\nbl/3ItnthSSHKCIiIiLCJOMUag7DSHpn+WiiGS6uNbNNZvYycDzQQ8KmJMPbC2wbVqEsIiKzV71R\nhhqCKJIOSYvlHuAUYAuwL3AucBOwkIRNSYp6ElVEREREmmWapWXSB/xGgT53vyZ+/aSZ/Rz4Ngmb\nklTq6MxocvaU02TuYVFe4VFm4Wkks446t6w6OzNkMsq7XXSOhWfKmdXfrGWzYTgwx8wy8ZhliGa9\neAY4KUlTkkqLFi2gZ3HDN6VlBmky97Aor/Aos/DUy6yrq7qv1w6LFy9U4TYDdI6Fp9HM5s6Z0sRs\nb2tGU5Ic0QwY3wGOJWpKchoJm5JUGhzIMadcmnxDmTGazD0syis8yiw8jWSWHR6tuX//1iwdurPc\nNjrHwjOdpiTTkbQpyYiZfQj4d6Jxy4PAxe7+VNKmJJWKRU0QHgpN5h4W5RUeZRaeepmV6mRZLJYp\nq1ZuO51j4Wm4KUn9Qcstmw0Dd3/N3U8HDgO6gLfePqYETUkqTfLNiYiIzD669ImkQjNbBa0mmh1j\nXKKmJCIiIiIiTTNDs2EAYGZfIiqMN8Svu0jYlOQd9Nu1iIjMQvUub9GnqhqHIdImrRuGYWZLgUuA\nL1e80cEkbEoiIiIiIjLTEhXLZtYJ3EX0UN9AxaouEjYlERERERFplukOVEg6DGMF8Ky7/6pqeQ41\nJdmlaDL3sCiv8Ciz8DSSWb2Z4To6dO1rJ51j4ZlqZpNs1bKmJJ8F3mdmZ8avFwH3AdcA85ralKRH\nk4SHQJO5h0V5hUeZhad+U5Lq+0o79PQsZE5nM5/Dl0boHAtPo5nNmYmmJO5+aOVrM/sL8BV3f9jM\njqRZTUkGc8zL6Cm/NNNk7mFRXuFRZuFpJLPh4ZGa+/f3Z1Ust5HOsfBMuSnJ2Aw0JZlA5aO7TWtK\nMqamJMHQZO5hUV7hUWbhqduUpM4Fvlgsk9F0UG2ncyw8DTclqb9Jy4ZhvIO7H1Tx9VbgzDqbi4iI\n7NJUkomkX+Ji2cxOBK4jmhauF7jW3W8zs25gDXAqMABc5e5rkr6fiIiIiMjUTe/X06RTx3UD9wM3\nuHs30QN/K83sw8CPaFYHv0num4uIiMw2uvSJtFXLmpLsDzzo7j8FcPdngceBE4g6+K1w90I8Vnm8\ng5+IiIiAxmGIBCDpbBjPAeeMvzazxcBJwHNAYYIOfp9M8n4iIiK7DlXSIs003U9rmvaAn5ktAh4A\n1hPdXf5q1SbT7uCnidnTT5O5h0V5hUeZhSdxUxI15GornWPhmXJTkvqbtXY2DDM7EFgHvAp8DjiM\nJnbwW9S9gJ6ehYmOUdpDk7mHRXmFR5mFp35Tkvk11/UsXsi8udNroiDTp3MsPI1m1jkTTUkAzOxo\n4GHgLne/LF72Kk3s4DcwkGM3zcueaprMPSzKKzzKLDyNNSUZrbl/f39WxXIb6RwLz1QzGyvMQFMS\nM3svUaF8nbtfO77c3bNmdj9N6uBXUlOSYGgy97Aor/Aos/AkaUpS7FDW7aZzLDwNNyWpv7plwzDO\nA/YArjCzFRXHciNwPvBDmtDBT0RERERkJiSdDWMVsKrWejP7LnAAcDjwVTN72d1/O9X30e93IiKy\nq9G1T6TJptnuumUjgc1sPtHsGKuBRcDNwANmNq0ZMURERGYbFcQi6dfKx+ZOAYrufpu7F939DmAz\n0xy3LCIiIiLSbq0slg9h59kvPF4+JWX1/BQRkVmo7uVNlz6RpirXP6naPwwD6CJqRFJp2o1JRERE\nRETarWkd/CaQA3avWjatxiSd6mKUeup8FBblFR5lFp5GMuuoc8sq00Hda1/vwHb27K6+zMKbfTk6\nOzMTrpPadI6FZ8od/GrfPK6rlcXyS8CFVcsMuKeRndese/F9wD4Al9/6mxfXXb98pLmHJ62gzkdh\nUV7hUWbhqZfZjx95ZS9g34nW/fN1Tzy/7vrlhVr7nvNv//kuILvu+uXlquULgfK665cPT++Id206\nx8LTaGavbx46jJ1v5I7rq7VfplXjgc1sHvBn4LtE8y2fDawEDnT37S15UxERERGRJmrZmGV3zwOn\nA58nqtYvBJapUBYRERGRULTszrKIiIiISOhaORuGiIiIiEjQVCyLiIiIiNSgYllEREREpAYVyyIi\nIiIiNahYFhERERGpoZVNSabFzN4P3AocDrwCfNndfzuzR7XrMbMTgeuAQ4Be4Fp3v83MuoE1wKnA\nAHCVu6+p2G8V8EWgE7gLuMTdy/G6s4Crgb2Ax4B/cvct7fuuZj8zey/wPHCuuz+kvNLLzJYQ/aw7\nGRgkOsduVmbpZWYnADcCS4E3iLK5V5mlj5kdA6x19yXx65ZkpJqlOSbIawnwfeAkIA/8DLjU3Qvx\n+rbmlao7y2Y2H3gAWA0sAm4GHjCzBTN6YLuY+IfK/cAN7t4NfBZYaWYfBn4EDAF7Ap8Bron/k2Nm\nFxHNrX0EcChwInBpvO7vgVuAM4E9gM3AHW38tnYVq4Geite3o7zS6hfAi8Bi4CPAlWZ2HMoslcys\nA1gLrHT3RcD5wJ1mth/KLFXM7DzgEWBuxeKmZ6SapTlq5HU3sAHYGzgK+ABwRbx92/NKVbEMnAIU\n3f02dy+6+x1E3+jHZvi4djX7Aw+6+08B3P1Z4HHgBGA5sMLdC+6+HvgJUXdGgC8A33P3LfFvcauA\nf4zXfR74hbv/zt1HgX8FPmpme7bpe5r1zOxLRBeDDfHrLpRXKpnZsUQXgW+4e8ndXwKOJ7pbqczS\nqZvo4jt+QS8Do0AJZZYaZvZN4GKiO4vjy1r1s/BUVLMkUiOvuUAWuDrOawtwD1ENAjOQV9qK5UOA\nP1Yt83i5tIm7P+fu54y/NrPFRB+FABTc/fXKzdmRT3V+DthE69y9H+ivWC8JmNlS4BLgy0AmXnww\nkFdeqXQ00b/vtWa2ycxeJiqWe1BmqRT/e94C3GdmBeAJ4CKiAlqZpcdqd38/8LuKZUtpbkZ98XpD\nNUtSO+UVF8jLqoYjLQN+H3/d9rzSVix3AbmqZTlAH2nMEDNbRPSxxXqiu8vV7cor86nOLwd0mNm8\nCdZV7yvTZGbjY7YudveBilVdKK+06iH6JG0LsC9wLnATsBBllkpmliH69zwD2B34ONH45XejzFLD\n3TdPsHgBzc1oe7yv8kuoRl7vYGY3ERW6q+JFbc8rbcVyjuiHUKUFRLfjpc3M7EDgf4C3iC4QWWC3\nqs0q86nObwEw5u75CdZV7yvTtwJ41t1/VbU8h/JKq1Ggz92vcfcxd38S+DnwbZRZWn0KOMbd18aZ\nPQQ8CHwLZZZ2rfpZqPxayMx2M7OfAacBJ7t7X7yq7XmlrVh+iZ0/fprotrm0mJkdDfwGeNjdPxmP\n/XkVmGdm+1Ruyo58qvM7JF620zoz24PowaaXkKQ+C3zOzPrNrB/YD7gP+AeUV1o5MCe+WzmuE3gG\nZZZW+wHzq5aNAU+jzNKuVdcu1SwtEg//fILoQbzj3P1vFavbnlfapo77NTDfzC4Efkg0+H4voqck\npU3i6cceBq5z92vHl7t71szuB1aZ2QVET6KeRfRUKkRPr15mZo8RXUQuJxoeAHAv8LiZrSEqCFYB\nD7n71nZ8T7OZux9a+drM/gJ8xd0fNrMjUV5p9CjRXY4rzew7wLHAJ4juoByAMkujR4lmBTrH3e80\nsw8SZXYqyizVWnXtMjPVLK2zFtgEnOHuxap1bc8rVXeW41vopxM9zdgHXAgsc/fqsUbSWucRPbRy\nhZkNxX+2xRf184F5wEbgP4Cvu/v4wPwfEE059xTwAvBfwA0QPTQY73sH8Cbwvvh9pPnK7HjI7wKU\nV+q4+wjwIaIieQvRD/+L3f0plFkqufsLwKeBr5nZANG0U2e7+zMosxA0PSPVLK1hZscTTSpwGjAQ\n1x/bzOzxeJO255Upl8vN+v5ERERERGaVVN1ZFhERERFJExXLIiIiIiI1qFgWEREREalBxbKIiIiI\nSA0qlkVEREREalCxLCIiIiJSg4plEREREZEaVCyLiIiIiNSgYllEREREpIb/B2OPfmek6eJKAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b8e24b134d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 1))\n",
    "\n",
    "xmax = snap25_counts.shape[0]\n",
    "xvalues = np.arange(0, xmax)\n",
    "\n",
    "ax.fill_between(xvalues, snap25_counts, y2=0, linewidth=0);\n",
    "ax.set(xlim=(0, xmax))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict, Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pad_location(start, stop, pad):\n",
    "    start = start - pad\n",
    "    stop = stop + pad\n",
    "    return start, stop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({(48, 8157): 78, (48, 8436): 2, (8275, 12200): 82, (8554, 12200): 2})"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "region_reads = bam_htseq[HTSeq.GenomicInterval(chrom, start, stop, strand)]\n",
    "bad_cigar = INSERTION_DELETIONS\n",
    "\n",
    "junctions = Counter()\n",
    "\n",
    "# pad = 10\n",
    "# start, stop = pad_location(start, stop, pad)\n",
    "\n",
    "length = stop - start + 1\n",
    "\n",
    "for read in region_reads:\n",
    "#     read = skip_bad_cigar(read, bam_filename, bad_cigar)\n",
    "    if read is None:\n",
    "        continue\n",
    "    for cigar_operation in read.cigar:\n",
    "        # N = did not match to genome and is an insertion\n",
    "        if cigar_operation.type == 'N':\n",
    "            junction_start = cigar_operation.ref_iv.start - start\n",
    "            junction_stop = junction_start + cigar_operation.ref_iv.length #- start\n",
    "            if (junction_stop < 0) or (junction_stop > length):\n",
    "                continue\n",
    "            \n",
    "#             print(junction_start, junction_stop, cigar_operation.ref_iv.length, \n",
    "#                   cigar_operation.ref_iv.start,\n",
    "#                   cigar_operation.ref_iv.end)\n",
    "            junctions[(junction_start, junction_stop)] += 1\n",
    "junctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def count_junctions(bam, chrom, start, stop, strand, bad_cigar=INSERTION_DELETIONS, offset=True):\n",
    "    region_reads = bam[HTSeq.GenomicInterval(chrom, start, stop, strand)]\n",
    "    junctions = Counter()\n",
    "\n",
    "    length = stop - start + 1\n",
    "\n",
    "    for read in region_reads:\n",
    "        if read is None:\n",
    "            continue\n",
    "        for cigar_operation in read.cigar:\n",
    "            # N = did not match to genome and is an insertion\n",
    "            if cigar_operation.type == 'N':\n",
    "                junction_start = cigar_operation.ref_iv.start\n",
    "                \n",
    "                if offset:\n",
    "                    junction_start -= start\n",
    "                    \n",
    "                junction_stop = junction_start + cigar_operation.ref_iv.length\n",
    "                \n",
    "                if (junction_stop < 0) or (junction_stop > length):\n",
    "                    # If any of the junctions start or end outside of the region, skip it\n",
    "                    continue\n",
    "                junctions[(junction_start, junction_stop)] += 1\n",
    "    return junctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def cubic_bezier(points, t):\n",
    "    \"\"\"\n",
    "    Get points in a cubic bezier.\n",
    "    \"\"\"\n",
    "    p0, p1, p2, p3 = points\n",
    "    p0 = np.array(p0)\n",
    "    p1 = np.array(p1)\n",
    "    p2 = np.array(p2)\n",
    "    p3 = np.array(p3)\n",
    "    return p0 * (1 - t)**3 + 3 * t * p1 * (1 - t) ** 2 + \\\n",
    "        3 * t**2 * (1 - t) * p2 + t**3 * p3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matplotlib.path.Path"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpl.artist.Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ymin = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bezier curve is specified by four points:\n",
    "\n",
    "![](http://matplotlib.org/_images/path_tutorial-2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8275 12200\n",
      "\t 83 88\n",
      "\t\t [(8275, 83), (8275, 93), (12200, 98), (12200, 88)]\n",
      "8554 12200\n",
      "\t 2 88\n",
      "\t\t [(8554, 2), (8554, 12), (12200, 98), (12200, 88)]\n",
      "48 8157\n",
      "\t 84 77\n",
      "\t\t [(48, 84), (48, 94), (8157, 87), (8157, 77)]\n",
      "48 8436\n",
      "\t 84 2\n",
      "\t\t [(48, 84), (48, 94), (8436, 12), (8436, 2)]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FNX6+PHPbkJIJSF06UE59KZUwcKVr4IFEAVB6Yjd\na+NaEAhFQcELcrGgFAUFxEJTuV69gvq7ClKVegQVBKkhvWd35/fHbkISigk7yexunvfrlddmZ2Zn\nnyeb3WfPzJlzbIZhIIQQQgDYrQ5ACCGE75CiIIQQooAUBSGEEAWkKAghhCggRUEIIUQBKQpCCCEK\nBJdmY6VUJ2CV1rqu535dYB7QA8gFPgKe1FrnedZPB0YDQcAS4AmttfSBFUIIH1XiloJSahTwBVCp\n0OL3gCNAHaAd0BGY4Nn+YaA30ApoDnQHnjQlaiGEEGWiREVBKfUc8AgwrdCySkA6ME1rnae1PgW8\nD3TzbHIPMEdrfcqzbjow0szghRBCmKukh48Waq1fVEpdm7/Ac4jo1mLb3Qrs9PzeDNhbaJ0Gml5q\noEIIIcpeiVoKWuuTf7WNUmouoHC3CAAigMxCm2QCdqVUSGmDFEIIUT5KdaL5fJRSobjPLbQErtFa\nn/GsygTCCm0aDji01rnePqcQQoiy4VVRUEpVBf4NpAJdtNYphVbvw91y2OK538yzrEQMwzBsNps3\n4QkhREXk1Qenty2FVcBxYIDW2lls3XvAOKXUBsABPIO7W2qJ2Gw2kpMzcLkCpwer3W4jJiYi4PKC\nwM1N8vI/gZpbSfOKjY306nkuuSgopbrivj4hG0hWSuVHuV1rfR3wOlAT+BEIAZYCs0vzHC6XgdMZ\nOC9qvkDNCwI3N8nL/wRqbmWdV6mKgtb6G9wf9Gitf8B9UdqFtnUBEz0/Qggh/IAMcyGEEKKAFAUh\nhBAFpCgIIYQoIEVBCCFEASkKQgghCkhREEIIUUCKghBCiAJSFIQQQhSQoiCEEKKAFAUhhBAFpCgI\nIUQAMAxzxkOSoiCEEAHApJogRUEIIQKBw+kyZT9SFIQQIgDk5BWf0ubSSFEQQogA4DBpjgUpCkII\nEQDMmmVOioIQQgQAp0vOKQghhPBwSktBCCFEvtw8aSkIIYTwyHNIURBCCOGRnp1nyrUKUhSEECIA\nOBwuU84rSFEQQogA4HQZpox/JEVBCCECQG6e05Txj6QoCCFEADiZlCXnFIQQQriHzc5zuEwZ6kKK\nghBC+DmXYZCamWvKoHjBpdlYKdUJWKW1ruu5HwMsAnoCycAUrfWiQttPB0YDQcAS4AmttUmjfgsh\nhMh3JjWbjOw8r/dT4paCUmoU8AVQqdDiBUAaUAO4E3jZUzhQSj0M9AZaAc2B7sCTXkcshBDiHEmp\n2eTket9SKFFRUEo9BzwCTCu0LALoC0zUWudprbcAy4Bhnk3uAeZorU9prU8B04GRXkcshBCiCMOA\nzBwHuSZc1VzSlsJCrXV7YGuhZU2BXK314ULLNNDM83szYG+xdU0vNVAhhBAX5nAaOMqrKGitT55n\ncTiQVWxZpmc5QITnfuF1dqVUSGmDFEIIcXGGYeAwYfjsUp1oLiYTCC22LBxIL7Q+rNg6h9Y6t6RP\nYLfbvAjP9+TnE2h5QeDmJnn5n0DN7WJ5GdhMGzrbm6JwAAhRStXTWh/1LFOcPWS0z3N/i+d+M8+y\nEouJifAiPN8VqHlB4OYmefmfQM3tfHnlOVwYBtiCgrze/yUXBa11ulJqDTBdKTUWdy+jwbh7HAG8\nB4xTSm0AHMAzuLulllhycoZpU8z5ArvdRkxMRMDlBYGbm+TlfwI1t4vllX8lc2pattfP401LAWAs\n8CZwFHfX1Ke01vkno18HagI/AiHAUmB2aXbuchk4TZqM2pcEal4QuLlJXv4nUHM7X175948lZHi9\n/1IVBa31N7g/6PPvJwGDLrCtC5jo+RFCCFHGssvrOgUhhBC+K3901NzyHuZCCCF8XVCQveDWZiv5\n4SMz+vhbzYyL16QoCCECRnCwneho96VS+bcllZSU4feFISvH4fU+5PCREKIIM2bvEtbIzJaiIIQw\nmRlj8ovy5n7NzBg6W4qCEKIIMz5YhDXMuKpZioIQoog8Pz+ubjaXy0V6ejopKcmkpaWSnp5GRkYG\nmZmZZGdnk5OTQ15eHg6Hw/JDbzm5Tm59co1Xn+tyolkIUUR2rgOobHUYpvn4449JTk4mKyuLhIQE\nQkND6d+/P0qpIts9+eTf2blzB5mZmZ4f9wd/VlbxcT8vzGazUaVKNNHR0URHxxTcxsTEUKVKNDEx\nMdSpcxmNGzehSZPLqVatGjabeWM0eVp51YFTl7oPKQpCiCL8raWQkZHBH38c5vDhQ+TmZjF69Igi\n6w3DIDU1ldzcXMaNG8emTZs4ffr0OUVh587t7Nixo8iy0NBQYmNjCQ+PICgoyHM9gIHL5cIwzt7m\n/56Xl0tqaiopKcnAYf5KlSrRxMXFERd3OXFxTbjiiqZ069adWrVql+pvUKyBEoMUBSFEIMvLy+Pg\nwQPs3bubAwc0hw4d4vBh909CwumC7dq3b39OUTh69CiPP/44s2fPZvPmzZw4cYK77777nOd4++13\ncLkgLCyciIhwwsLCCbqEAeYMwyA9PY3k5GRSUlJISUkmOTmZ1NQUkpKS+PPPI/z660F+++1Xjhz5\ng507d7BzZ9Fi1Lx5C669tifXXXc9XbpcTXh4qbrXVil10IVIURBCFGF1j9STJ0+yd+9u9u7dU3D7\nyy/7ycs7//zDtWvXoWHDRjRo0JBu3bqdsz42NpZ33nmHkJAQFi1aRPfu3dmzZw8tW7Yssl1cXBNT\nrlOw2WxERVUhKqoK9etffNvc3FyOHDlcUCR2797FN99sYN++vezbt5c335xHSEgInTt3Y9SoMQwb\nNrgkIXh17M9m9YmRizASE9MDakCroCAbsbGRBEJewcFFz2XZ7Taio8NJScn8y5Ep/ekCoUB6zQq7\nUF5Ol4vDJ9KJu8yrL5sllp2dzY4d29i06Xs2bfqeXbt+LvLNP194eATNm7egRYtWNGvWjMaN42jY\nsDH16tUnLOzstC3BwXaqVr20IbN95eI1wzDQej/ffPM133yzge+//39kZrrnK2vevDmPPvoEffsO\nIDj47Hf63Dwn97/yTf7d69e90nfjpT6/tBREqV3sjVeSq0h95c0nzmUYYFB2xS8tLZUtWzbzww/u\nIrBjxzZyc4vOu9WoUWNatmxNixYtadGiFS1atKRhw0bY7RWjs6TNZqNZs+Y0a9ac++57iJycHP7z\nn/XMmTOLXbt+5oEH7mXevLksWPAOTZpcAVD8FfPqzLUUBSFEAbNbQ4Zh8PPPO/nss3X8979fsmfP\nLlyFpowMDw+nc+dudOnSla5dr6Zdu/ZERkaZGoO/q1y5Mrfe2o++ffuxefN3PPvsc/z880/ccMO1\nLFq0lOuv/1vxh0hREEKYI8fhBANchoH9ErtKOp1OtmzZzGefreWzz9Zx9OiRgnVVq1alc+eudO7c\nja5du9G6dVsqVapkVvgBzWaz0adPHzp06EJ8/ETmz3+NMWOGs379f2nY+PIim3rzPFIUSsnpdJKR\nkU5aWhppaWmkp7tv8/JycTpdOBwOnE4HTqfT8/vZW5fLSeXKQaSlZZKX58DhcP+4XE4cDidOp5Og\noCCCg4Ow24MIDg4mKCiIoKDgc5afXRfkWRdc6DF2goKCCQ0NJTQ0lLCwcMLCwgr9hMsbUZxX/nj8\neQ4XlSuVvOeNYRh8++1G1q5dzfr1nxY5L9CiRStuvvlWbrrpZlq2bFWmh4EcDhcpKZklPr9V/LH+\nICQkhKlTp+NyOXn77Td59tmnWLZideFNpChciszMTI4f/5OjR49y7NifnD59yvNBn0p6enqhD/yz\n99PS0sjM9H5mI18QFBR0nmLhLhgXKiTh4eFUqVKFuLg4Bg0aUGR/n3/+OUePHmXlypVce+21hISE\nMGrUKGrUqGFRhuJSOBwuDDyHkUrwvSErK4uPPvqAt956Ha33Fyy/6qpO3HzzbfTpcwuNG8eVXcDn\n4fRMTel0ugKqc0BhhmHQqlUbwsLC+O67bzh06FDh1VIUisvLy+P48WMcO/Ynf/55tODHff9Pjh07\nSmJiYqn3GxERSe3adYiKiiIqKoqICPdtZGQklStXLvjGnv/tvvC39uDgYCpVCiYqKpzcXCd2e+EW\ngPvWbrfjcrmKtC6cTofnd1exFoijUMvk3O0cDge5uTlkZWWRmZlFdnYWWVlZZGW5L83Pysr0XOF5\nutSX5rdv3/6cotCnTx/+/e9/88QTT7BhwwYqVapUpFdIYatWfURmZrbnSs+qVK9enWrVqhEdHVNh\nTib6qvy5fvMcTi728XDq1CkWL36bd95ZwJkzZwB3i2Do0BH06XMLdepcVh7hVkipqSnMnPky/frd\nzq239mPlyuXs2LkDiM3fpGIXBcMw+P33X9m2bSvbtm1h+/at7N69C4fjwkPIBgcHU79+Ay67rC51\n69blssvqUatWLaKjY4iIiCz40I+KqkJkZKSnAER6/YHli90bDcMgJyenSLE4XxHJzMwkIyODtLRU\nqlQ5f3fFrVu30q1bN4YPH05wcDDffvstffr0OWe7efNePefKUXC3XqpWjaV69erExlajWjV3sYiN\nrVZsmXt5tWrVCQkJMf1vUpE5XQYYkJnjIDry3O7up06dYsaMqaxcubyg19ANN/wf99//MD16XGvq\nkA3iXCkpKcyY8SLjxj1LlSrRzJw5hxenzyYt28V/fv0xf7OKN/aR1vv57LO1bNmyme3bt5KUlFRk\nfXR0DHFxcVx2WT3q1avHZZfV83z416Vu3XrUrFnrkq5UDEQ2m63g3ENJFb9GAdyH42JiYmjQoAEf\nffQRoaGh3Hnnned9/KhR97J/vyY1NYXExEQSE89w5kxCwc/5+qlfSFRUFWJjYz2tjepUr16D2rXr\nULt2HerUuYw6depQq1YdatSoIa2QEnA6DZwuFxnFxuU3DINPPvmQ8eP/QWJiIqGhoQwdOpL77nuQ\npk3VBfYmzJSamsKcOTN55pnniIhwfzELCwsjNz2H3KJfgitGUTh+/BgrVy5n1aqP2bt3d8HyoKAg\n2rRpR4cOV3LllR258sqOxMU1kQ+AchYeHs7DDz8MwDPPPHPRbe++e9gFT+o5nU6Sk5M5cyaBxMQz\nJCQkFCsaZwpu85fnD3dwMcHBwdSqVZvatWtTu7a7WBQuHu7bOhW+O6TD6cJpGLgMF1k5DsIqB3Py\n5AmeeurvfPHFegDuvPMu4uNfkPNF5ch9yGg6L730Ik5nUJEjDS4DsnMqUFHIyMhg3rw5vP763ILR\nCps0uZx+/QZw7bU9adOmbWnHBRE+LCgoyHNoqFqJtjcMg4yMjIICcerUKU6cOF7wc/z4MU6cOMGJ\nE8cKzi1dTGRkFLVr1y4oFI0aNaJly2ZUr16b+vUbUatW7YD+wpHrcGGzgWEzyM51cuzo7wwc2J8j\nR/6gTp3LmDVrDr163WR1mBWK0+lk1qyXePrpZ4mOjiYxMf2cbXLyzOs55dNF4bfffuW22/pw/Pgx\nKlWqxL333s9dd91Dq1at5dilANyHvyIjI4mMjKRBg4YX3TY7O5uTJ09w/PhxTpw45ika+b+f8BSQ\n4xw8eICDBw+cdx+VK1emQYOGNGzYqNBP44KxdyIjI8sizXJzKimTGjFhOF0GBw8eYPjdt5GQkECv\nXjfyxhsLqFIl2uoQK5zXXnuV4cNHEh0dc8FtzJhcJ5/PFoW0tDSGDh3M8ePHaNu2HQ888AgdO3am\nfv0GVocm/FRoaGjBB/mFGIZBSkoyx48f5/jxPzly5A9OnvyT/ft/4dChQxw69DsHDvzCgQO/nPfx\n1avXKHiORo3cBaNx4ziUakbVqrHnfYwvScvMo1p0KDk5OTzzxP0kJCQwYMBA5s59Q65tscCaNZ/Q\nqlXrguEsLsTpqgAthddee439+/fRtGkzZs58lczMDP73v+84cuQPbDYbbdu245prrqdy5cCZDERY\nz2azERNTlZiYqjRv3uKcHmOGYZCUlFhwHuPwYXehyP/9zz+PkpBwmm3btpyz71q1anvGtGlB8+Yt\naNasOU2bNvOp1kVWjgOH02DNR++yb8/PtGnTjldffV0KggV++mkHCQmn6dv39r/c1szejD5bFK6/\n/nqGDhvJ8GEjadOmHQDdunUH3MfYfvppB/PmzcFut3Pbbf3+spIK8zgcLpKSil7EF6ijpBZns9mI\njXV3k23f/spz1ufm5nL06JEiRePXXw+wf/8+Dh8+xMmTJ/jmmw1FHtOgQSOaN3cXi/yicfnlV1jy\nhScrx0FOroMv138MwIwZs6TbrwVOnTrFunVrGD9+Uom2d1SElkLnzp2pXvtyqoSf+w8ZFBREhw5X\n0aHDVWRmZrJ27SpWrlxOixat6N37FvknLgfFP9iDgtzneAL5KtKSCAkJIS6uCXFxTc5Zl5GRwYED\nmv3797Fv31603sf+/fv4449D/PHHoYLePeD+H2/S5HKaNWtBu3Yd6Ny5C23atCvzQpGcnkPC6ZMc\nOXSQ2nXqceWVHcv0+cS5cnNzmTv3FcaPjy/xudPcinCi2TCMEk0LGB4ezl13uWdR2r17F3Pn/pNK\nlSpx++13yvkH4VMiIiJo164D7dp1KLI8JSWZ/fv3s3//Xs/PPvbt28Mvv2h++UWzdu0qwH2Su337\nK+nUqQudO3ehY8fOxMRUNTXG7Dwn4RExPDP1dWpUjZIOHRaYPXsmDzzwyAVHBDifXPfczKbwuigo\npboBrwJNgWPAFK31cqVUDLAI6Akke5YvKul+DaNgEuoSa9WqNa1atSYjI4OVK5dz+vQpBg++R4qD\n8GnR0TF07uz+oM9nGAanT59m797dbN36I5s3/8DWrVsKJqPJ16xZczp2dD+2U6cuNGzYyKsP8uwc\nJ6Fh4bRsfw31avjOuY6KYvXqj7n66h7UrVuvVI9LSssxLQavioJSyg6sAu7XWq9SSnUHvlZK/Q/4\nJ5AG1ADaAeuVUru11j9eeI9nuSfCvrTDEBEREYwcOYasrCyWL3+PhITTDBkylHr1/mJuPCF8hM1m\no2bNmtSs2ZPrrusJgMPhYN++Pfz44yY2b/6BzZs3sX+/+xDU0qWLAffJ7K5du3Hzzbdxww03EhFR\nulnIHE4XTqeLXD8+7+OvDh36nSNHjtCv34C/3riYhJRs0+LwtqUQA1Tn7HiKBpADuIC+wOVa6zxg\ni1JqGTAMKFFRcHnGdPdGWFgYo0bdS2ZmJitWvM+ZMwkMHnyPFAfhl4KDg2ndui2tW7dl9Oj7MAyD\no0ePsHnzD55CsYn9+/eyevUnrF79CWFhYdxww4307dufv/3t/0pUIAwgz2kUjDQqykdeXh4LF85n\n4sSpl/T4wyfTTIvFq6KgtU5USr0BrFBKvY97dL7RuAtFrtb6cOHNgf4l3bfhnhfQFOHh4QXFYfny\n90hMPMOQIUNL3UQTwpfYbDbq129A/foNuOOOQYD7/MTGjV+zdu1qvvrqC9atW826dasJDw/3FIjb\n6d37ZoKCLtDF1HC3FtKy8soxE/HGG/MYM+b+S+76W5LzryXl7eEjG5AJDADWAf8HLANuA7KKbZ4J\nlHg8CpdhYLOf7dVihqioCMaOvY+MjAyWLXuP5OQkT3Goa9pzXIzdbityG0gCNTd/yys2tiq33z6A\n228fQEZGBl9++QWrV3/CV1/9h7VrV7F27SoaNGjIQw89zEMPPXBOXgbgchmkZ+VhN/n9V1787TX7\n9ttvaNy4EXFxjS+63YXyCgqyFQx5bgZvDx/dDnTSWv/Dc/9zpdSnQDxQfNjNcODcQTsuICfXSURE\nKLGx5p/sio2N5OmnnyQjI4NFixaRlJTE6NGjy604xMSU7jivPwnU3Pwxr9jYSEaNGsaoUcNIT0/n\ns88+Y/78+WzYsIGnnx7Hyy/PYPz48Tz66KMFowbb7TZCKlciK9dJVFRYmbz/yos/vGYJCQls3foD\n06ZNK/FjiudlBAVd8vnX8/G2KDQAinecdgDbgO5KqXpa6/wRyBSwt6Q7djhdpKVln3fwJzPdffdI\nMjIyWLJkCWlpaYwePabMhiOw223ExESQnJxh6ovoCwI1t0DKq1evm+nV62a2b9/GvHlzWLt2DU88\n8QQfffQxr702nwYNGuJ0usjMzCE1LYe0tKwyf/+VBX95zQzDYPLkqYwb92yJ/s4Xyis5LQeHD13R\n/CXwolJquNb6XaXUtUA/3N1QGwHTlVJjgVbAYODcGVcuIM/hwukqnwuhQkPDGTPmfpKTk1i40D3o\n19ChI8rsIjiXywjYC7wCNbdAyqtt2w4sWrSUw4cPMHjwEL7//n/06NGVf//7a8A9uFqe04XLZe7w\nCeXN11+z9957l759BxAWFlGqOIvn5XQaXnfKKcyrMYC11ruBO4DHlFLJwL+AYVrr7cBYIAQ4CnwI\nPKW1PndAmAtwOF3lXuVjYqry+OPj6NnzBmbOnM5nn60r9VSVQviL9u3b8/XX3zF06AjS09MYO3Yk\nebnu/u5mjropzrV37x5yc3Np27a91aGcw+uL17TWnwGfnWd5EjDoUvfrcl36dQreatw4jvHjJ7Fp\n0w9MnjyBfv1uP+cqVCECQWhoKC+99E/27dvL1q0/EtX4P7RvNoLM7AtPZyu8k99FPj6+5OcRypPP\nDnMB1n9b6dKlK506dWbNmk9Yv/5Thg0bJd1YRcAJDg5m0qRpfPfdRnYn1cTlKtkQM+LSzJs3h4ce\n+rvPTtbk00UhJSPX6hCw2+30738HWVlZLF26mJycXEaOHF3hp20UgSV/mI3H5n6HYRimdnEUZ61f\n/xlXXdWRWrVqWR3KBflmqcK069ZMExYWxtixDzJo0BDeeGMey5e/h9Np3iBUQvgKA+tb6YHozz+P\nsn//Xnr27GV1KBfls0UB3IPi+ZqaNWsybtyztG9/JS+8MJkNG/5rdUhCmMqQgmA6p9PJm2++xoMP\nPmp1KH/Jp4uCL2vWrDkTJ07BMAymTJnI/v37rA5JCFNITTDfwoXzGT58VHlNnOTVpdw+XRT8oTto\nz5438NxzE9m+fSszZ04nKSnR6pCE8Io/vO/8yZ49uwkLC+fyy/1jdkjfPdFs3nh4ZS44OJghQ4aS\nkpLMggXzqVu3HgMHDvbZ3gVCXEx2rpwrM0teXh4ffLDMZ7ufno9vf2r5S1XwiI6O4cknn6Zly9bE\nxz/Pzp3brQ5JiFLLzJFrFMyyePHbjBgx2q++IPp0pH5WEwq0bt2GyZNf4ODBA7z88oskJp6xOiQh\nSiwx1bwJWyqy/fv3UalSyHnn6/Zlvnv4CHC5/LevtM1m4447BpGamsKCBfOpU+cyhgy52+qwhPhL\nGXI1s9ccDgfLli1l0qRLmzTHSj7bUjCA9Cz//+esUiWaJ574B23atGPy5Als2VLi4Z+EsESOnFPw\n2rvvLmTYsJEFQ5L7E58tCmDubEJWa9myFfHx0zh48CAzZrzAmTNySEn4pvwRNw2/PYBrrYMHDwD4\nTW+j4ny6KOQ6Ausbi81mY/DgwTz44CMsWbKI999fIldFCxFAnE4nS5YsZsSIMVaHcsl8uigEUkuh\nsCpVqvD44+Po0OEqpkyZyLZtckhJiECwZMli7r57mF8eNsrns0UhITmLI6f8b9an0mjevAXx8dP4\n44/DzJgxjYSEBKtDEkJcot9++5W8vFyUamZ1KF7x2d5HTqerQpzwstls9O9/B2lpqSxc+BbVqlVn\nyJChfv1NQ4iKxuVy8c47C5kwYbLVoXjNZ1sKOXlOcgP08NH5REVV4bHHnqJTpy5MmTKRHTu2WR2S\nEKKE3n9/CYMGDaFSpUpWh+I1ny0KgXo+4a8o1Yz4+GkcOPALs2bNIC0t1eqQhBAXcfjwIdLS0mjZ\nspXVoZjCZw8f5eYF/qGjC7HZbAwcOJgzZ87w2muv0rp1O/r0uQWbzavBD4UQJjMMg0WL3mb8+ElW\nh2Ian20pWDU/sy+pVq0azzwzgejoaKZMmcjRo0esDkkIUciKFe8zYMCdhISEWB2KaXy2pSAzP53V\nvfs1XHVVJxYufIuwsDCGDRtJcLDPvnRCVAhHjx7h9OnTtGnTzupQTOWzLQUpCkWFhoby0EOP0qPH\ntUyZMpGfftphdUgiAMm7rmQMw+Dtt9/kvvsetDoU0/lsURDnd8UVTZk8+QX279/HK6+8RHp6mtUh\nCVHhfPjhCvr27V9eM6mVKykKfshmszFo0BCGDx/Nv/41m/XrP7M6JCEqjBMnjnP8+DE6dLjK6lDK\nhBQFP1a9enWefXYikZGRTJkykT//PGp1SEIENMMwmD//dcaODbzDRvnkbGUA6NHjWjp27MyCBfOJ\niIjw2yF7hfB1q1Z9RO/etxAWFmZ1KBfjVd91aSkEiNDQUB5++O9cfXUPJk+ewM8/77Q6JCECysmT\nJzl06Hc6depsdShlyuuWglKqLvAmcA2QAszUWv9LKRUDLAJ6AsnAFK31Im+fT1xc06aK+PhprFy5\nnK+++g9jxz5IZGSk1WEJ4dfch41e46mnnrE6lHMYhrl9xsxoKawG9gBVgRuBSUqpLsACIA2oAdwJ\nvKyU6mTC84m/YLfbueuuuxk6dCRz5/6TL75Yb3VIQvi1Tz9dww03/B/h4eFWh1LmvCoKSqnOQB3g\nWa21S2u9D+gKHAP6AhO11nla6y3AMmCYtwGLkqtRowbPPTeRypUrM21aPKdPn7Y6JCH8TkJCAlrv\np1u37laHUi68bSl0APYCM5VSx5VS+3EXhVggV2t9uNC2GvDvgcb91HXX9eSxx55i6dLFfPzxStOb\nm0IEsjffnMf99z9sdRjlxttzCrHA9cBXQH2gI7AeuAXIKrZtJlCqtpcNCAoKnEHg7HZbkdvyFB0d\nxbhxT7Nly2amTp3AvffeT7169U3bv5W5laWKnJfdbvPL95+Zr9kXX6znmmuuJTo6yut9eetCeZn9\nGnlbFHKAM1rrlz33f1BKfQJMBkKLbRsOlGoqNXuQjdjYwDtJGhMTYdlz33jj37juuu7MnTuXmJgY\nRo8ejd2nVQOHAAAZiUlEQVRuXic0K3MrSxUlr8IfOFFRYX79/vP2NcvIyGDXrh28+OKLJkVkjuJ5\nuezmdj/3tihoIFgpZdNa5x+TCAK2Az2UUvW01vlXVCnch5pKzOk0SEwMnCk57XYbMTERJCdnWD4K\n7OjRD7Bnz24eeeTvDBs2kiuuaOrV/nwpNzNVtLwK/56WmuWX7z+zXrPZs2cxdOhon/kbXCiv5NRs\nU5/H26LwJe7DQpOUUlOBzkA/oBfQCJiulBoLtAIGA31Kt3sDpzNw3oj5XC7fyKtZs5Y8//wU3n13\nIV999SUjR97r9cxRvpKb2SpiXv6eszfx7927h2rValC9ek2f+xsUz8vs+Lw6bqC1zgauw10MTgHv\nAY9orX8ExgIhwFHgQ+ApTy8k4UOCg4MZPfo+brjhRqZOncSuXT9ZHZIQlnK5XKxY8T6DB99jdSiW\n8PriNa31b0Dv8yxPAgZ5u39RPuLimhAfP40VK95n48YN3Hvv/YSGFj8tJETg++CDZQwcONjUc23+\npGJmLc7LbrczZMhQBgy4k5deeoFNm36wOiQhytXp06c5ceI4rVq1tjoUy0hREOe47LK6TJw4hRMn\njjFr1gyZs0FUGAsWvMG99z5gdRiWklFSxXnZbDb69RvAmTNnmDt3Nl26dKVnz15WhyXEJQkO/uvv\nvzt2bKdnzxuIiamCw+Eqh6h8kxQFcVHVqlXjuecm8uWX/2b69Cncd99DxMZWszosIUosONhO1ap/\nfc1Cz549Cn5PSsqosIVBioIokV69bqJLl2688cY8mjdvwS239MVm87+rXcXFyQgoQs4piBKLiqrC\nP/7xHNWqVWfq1EmcOnXqnG2CguwEB5fsRwjhe6SlIEqtW7futGvXgfnzX6NJk8u57bb+Beuio0s+\nvFVFbqIL4avk65q4JOHh4Tz++Dhq1KjJtGnxJCTIsNzCf6xevZo33niDqVOnWh2Kz5GiILzStevV\nPPbYUyxevJCvv/7a6nCEKJFu3brxwAMPyHmx85CiILwWGRnJuHHPyLSfwm/ExsYyd+5c7rjjDqtD\n8TlSFIRpOnWS2VaFf5g6dSqpqals3LgRp9NpdThm86r5IyeaRZlavXo1x48fJyEhgQkTJlgdjhAA\nTJ482eoQfJa0FESZkmO3QvgXKQqiTMmxWyH8ixQFUaYC/NitEAFHzimIMiXHboXwL1IUhGV27NhO\n69btrA5DBDiHw8WxYwksX/4eo0bdW+LHVFQ+ffhIBucKbOnpacyaNYOMjAyrQxEB7vXX53HTTbfi\ncLhK9FOR+XRREIGtR49rGT58NK+88hKbNn1vdTiimED5TrZ37x6qV69BzZo1rQ7FL0hREKZKSckk\nKSmjRD8Oh4saNWowYcJkTpw4zj//+TKZmZlWpyACiMvl4oMPljF48D1Wh+I35JyCMJXT6cLpLN13\nzPxZ3k6ePMHMmdPp0+cWOnbsXEYRiorkgw+Wceedd2G3y/ffkpK/lPAZtWrVZuLEKfzxx2HmzJlF\ndna21SEJP3b69GlOnjxBq1atrQ7Fr0hRED7FZrMxYMBABg0awvTpU9mxY5vVIQk/9dZbbzBmzP1W\nh+F3pCgIn1SnzmXEx0/jl180c+fOJicnx+qQhB/ZsGED7dt3kJF7L4EUBeGzbDYbgwYN4fbb7+DF\nF6fw8887rQ5J+IG8vDzWr19Pnz63WB2KX5KiIHxevXr1iY+fxu7du3jttbnk5uZaHZLwYUuXvsu9\n95bsIjVxLikKwi/YbDaGDBnKLbfcxrRp8ezevcvqkIQPOn36NCkpyVxxxRVWh+K3TCsKSqlaSqmT\nSqk+nvsxSqlPlFLJSqlDSqlRZj2XqLgaNmxEfPw0tm3bwoIFb+JyVeyrT0VRixa9xZgx91kdhl8z\ns6WwEIgtdH8BkAbUAO4EXlZKydRcwmt2u53hw0fRvfu1TJo0nmPH/rQ6JOEDfvppB3FxTYiKirI6\nFL9mSlFQSt2HuwAc8dyPAPoCE7XWeVrrLcAyYJgZzycEQLNmzXn++XhWrlzOunVrrA5HWMgwDFat\n+pgBAwZaHYrf87ooKKWaAk8AD3B2btArgFyt9eFCm2qgmbfPJ0RhlStX5rHHnqJKlSpMnz6F9PQ0\nq0MSFli16iNuu62fXLlsAq+GuVBKBQFLgEe01slKqfxVEUBWsc0zgfDS7N8GBAUFzjSOdrutyG0g\nsTq3nj170q5dW2bPnsnNN99Kx47mHKm0Oq+yUpK87HabX7z/0tPTOXBAc+ed7lZCRXvN7Ca/Rt6O\nfTQR2KG1/k+x5ZlAaLFl4UB6aXZuD7ITGxt4F5/ExERYHUKZsTK32NhIXn31nyxdupR33tnJo48+\nSnCwOcN7BeprVjyvwnNpR0WF+sX7b8GC13nqqcfPibWivGZOW1DxTbyqEt6+YwYCtZVSgzz3o4EV\nwMtAiFKqntb6qGedAvaWZucup4vExFLVEZ9mt9uIiYkgOTkDlytQBiZ286Xcbrnldn799SCPPvoY\nY8bcR6NGjS95X76Ul5kulJdRaBKTtLRsn3//HTr0O4YRREhIZEGsFe01S041d4wwr4qC1rp54ftK\nqd+BB7XW65VSbYHpSqmxQCtgMNCnNPs3oNQjbvoDl8sIyLzAd3Jr1KgJzz8/hfnzX6d27doMGDCw\nyLfg0vKVvMx2sbz8Ied33lnM00+PP2+c/hD/pSiel8vkHM0+K2NwtukyFggBjgIfAk95eiEJUS4q\nVarEww//nXr1GvDCC5NJTk6yOiRhoo0bv6Zr126EhIRYHUpAMXU+Ba11XKHfk4BBF9lciHLRpUtX\nWrRowb/+NYfrruvJ1Vf3sDok4aW8vDy++WYDkyZNtTqUgCP9t0SFUKVKNOPHT+L06VPMnftPGT/J\nz73//hLuuUcuewLz57KXoiAqlH79BtC//x1MnTqJgwcPWB2OuAQJCQmkpCTTpImMb1QWpCiICqd+\n/QbEx0/j66+/ZNmypUV63IizfPXPsmjRW4waJaOglhUpCqJCCgoKYuzYB2nevAWTJ08gISHB6pBE\nCeza9RONGjUmKqqK1aEELCkKokJr3/5KnnrqGRYunM/XX39ldTjiIgzD4MMPP+COO6T/SlmSoiAq\nvMjISJ5+ejw5OTm88spLZGebezGQMMeaNZ/Qt29/Gd+ojMlfVwiP3r1v5u67hzF9+lT27t1jdTii\nkIyMDPbv38uVV3a0OpSAJ0VBiEJq167DpElT+fHHTSxevEAm8fERixcvYORIOblcHqQoCFGM3W5n\nxIjRdOrUxTOJzzGrQyo3vtgT6/DhQ4SGVqZWrdpWh1IhSFEQ4gJatmzF+PGTWLRoEZ9//qnV4VRY\nS5YsZujQkVaHUWFIURDiIkJDQ3n++ecJDQ1lxoxppKf79qihgebbbzfSpUtXKleubHUoFYYUBSFK\noGfPGxgz5n5mzZrBzp3brQ6nQnA4HGzY8F969brJ6lAqFCkKQpRQ9erVmTRpKrt2/cw77yz0yePv\ngeS9997l7rtlfKPyJkVBiFKw2WwMHTqCK6+8ismTJ5CUlGh1SAEpISGBpKRELr9cxjcqb6YOnS1E\nRdG6dVsaN45jzpxX6NXrJjp37mJ1SAFl0aK3uP/+h6wOo0KSloIQlygyMorx4ydx6NBvvPXW63JN\ng0l27fqZBg0aUqVKtNWhVEhSFITwgs1mY9CgIVxzzfXExz/PqVOnrA7JrxmGwccfr2TgwMFWh1Jh\nSVEQwgTNmjXnmWeeZ9Git/j2241Wh+O3Pv/8U3r3vkXGN/LOpU9GjhQFIUwTHh7OM888T2LiGebN\nexWHw2F1SH7F4XCwfftWOT9jMSkKQpisX78B9OlzM5MnT+D48YozRIa3PvhgGYMGDbE6jApPioIQ\nZSAu7nKefz6eZcuW8tVXX1gdjs/LyMjg6NEjNG2qrA6lwpOiIEQZqVy5Mk8++TR5eQ7mzJlFbm6u\n1SH5rKVLFzNsmIxv5AukKAhRxnr3vpkBAwYydepEDh8+ZHU4PichIYHc3Dzq1LnM6lAEUhSEKBf1\n6zdg4sSprFnzCZ9+utbqcHzKu+8uZMSIUVaH4bcMzB1uRYqCEOWkUqVKPProE4SHhzFz5nQ/mPaz\n7Md2Onz4EFWrxsqFaj7Ep4uCjDcmAlHPnr0YOnQEL7wQz6FDv1sdjqWWLVsig975GJ8uCkIEKve0\nn9NYs+aTCts7adeun2jatJnMleBjvC4KSqnuSqlNSqlkpdQBpdRYz/IYpdQnnuWHlFJy0FCIQoKD\ng/n7358kOzuH11//F06n0+qQytWqVR/Tv/8dVochivGqKCilYoA1wGytdQwwEHhRKfU34G0gDagB\n3Am8rJTq5GW8QgScW265jV69bmTy5AkkJp6xOpxysXHj13Tvfo0MZ+GDvH1FGgKfaq0/ANBa7wA2\nAt2AvsBErXWe1noLsAyQg4dCnMcVVzTlH/94jtdemxvwM7u5XC6+/XYj11//N6tDEefhVVHQWv+k\ntR6ef18pVRXo4bmbp7U+XHhzoJk3zydEIIuMjOT55+PZuXMHK1a8b3U4ZWbNmk/o27c/NptX47aJ\nMmJa200pFQ2sBbbgbi1kFdskEwg36/mECEQ2m40RI0bTqFEcM2ZMJSur+NvIv+Xk5LBv317atm1v\ndSjiAkyZeU0p1RhYBxwA7gJaAKHFNgsH0ku776CgwPk2YbfbitwGkkDNzaq8rr66G02axDF9+hTG\njr2fhg0bmbr/kuRlt9tMf/+tWPEe99wztEzf1xXtf9Hsv6XXRUEp1QFYDyzRWo/zLDsAhCil6mmt\nj+ZvCuwtzb6DgmzExkZ6G6LPiYmJsDqEMhOouVmRV2zs5cydO5uXX36ZDh060Lt3b9Ofo3hehQ/p\nREWFmfr+S01NJTMzlSuvbGPaPi+movwvOky+ssCroqCUqoW7IMzSWs/MX661TldKrQGme7qotgIG\nA31Ks3+n0yAxsdSNC59lt9uIiYkgOTkDlyuwrswL1Nx8Ia8HHvg7a9euYerUF3nggYdN6bFzobyM\nQleMpqVlmfr+mzv3VQYPHlrm72lfeM3KwoXySk4x9xCjty2FUUB1YIJSaqJnmQG8CtwLzAeO4u6a\n+pSnF1KpOJ2B86Lmc7mMgMwLAjc3q/O6+ebb+OUXzYQJ43n88aeoWjXWlP1eLC8zcz558gQ2WxCx\nsdXL7e9o9WtWVornZXaOXhUFrfV0YPpFNhnkzf6FEGc1bar4xz+eZfbsWfTt2582bdpZHVKJvfvu\nIh588FGrwxAlIFeOCOFHIiOjeP75eLZt28rKlcutDqdEDh8+RI0aNYmMDLzzg4FIioIQfsZmszFy\n5Bjq1LmM2bNn+vzwGCtWvM9dd91tdRiihKQoCOGnevS4lr59+zN58gTS0lKtDue8Dh8+RM2atQgL\nC7M6lIrEqz6qUhSE8GNxcZfzxBPjePnl6fzxx+G/fkA5k1aC/5GiIISfi4mpysSJU/jgg2Vs3rzJ\n6nAKSCvBP0lRECIAVKpUiXHjnuWXX/bz8ccrrQ4HkFaCv5KiIEQAGTp0BLGx1Zg7dzYul8uyOKSV\n4L+kKAgRYK6//m/cdFMfpkyZSHq6NSMCSCvBf0lRECIANW2qePTRx3nppRf488+jf/0AE0krwb9J\nURAiQMXGVmPChMksXfoO27aVeoQZjEscPUFaCeXM5JE8fLooGGZnK0QFExISwtNPj2fXrp9Zs+aT\nMn8+aSX4P58uCkII7+VP3BMREcGbb84rMhKq2T78cIW0EvycFAUhKogbbriRzp27Mnv2zAsWBm/q\nRVJSIqGhYdJK8HNSFISoQNq3v5LevW9h+vRppo+ZtHLlcgYOHGzqPkX5k6IgRAXTvHkLBg0azPjx\n48nLyzNln3l5eSQlJVGzZk1T9iesI0VBiAqoSZPLefDBB5k6dRLZ2dle72/dutXcdlt/EyITVpOi\nIEQF1aBBAx588FFeeCHeq4vcDMNg3769tGjR0sTohFWkKAhRgdWuXZvHHx/HjBlTSUlJvqR9bN78\nA507dzE5MmEVKQpCVHCxsdV4+unxzJw5nezMlFI//uuvv6Jnz15lEJmwghQFIQRRUVV47rlJ7P9+\nJdnpiSV+3O+//0aDBg2x2+WjJFDIKymEACA8PJzW143gt22rycko2aGkNWs+4fbb7yzjyER5kqIg\nhCgQFFwJ1e1uDv74EXl5uRfd1jAMcnJyCA8PL6foRHmQoiCEKCKoUmXirurL+++8dtEhMXbu3E67\ndh3KMTJRHqQoCCHOERZVg6s692DZsqUX3Gbjxq+57rqe5RiVKCGbNw+WoiCEOK9Wba4iKCiI//3v\nu3PWGYZBbm4ulStXtiAyUZakKAghLuiuu+7mu+82kph4psjyn3/eSdu27S2KSpQlKQpCiIu6776H\neO+9d4ssk0NHgSu4LHeulGoPvAm0BH4BHtBaby7L5xRCmKtq1ViCgytx8uRJatWqBYDT6SQ0NNTi\nyERZKLOWglKqMrAWWAhEA/8C1iqlpP+aEH6mT59b+P77s+cWgoKCLIxGlKWybClcDzi11m957i9W\nSj0O9AE+KsPnFUKYrGHDRtSt16Dg/ojRD1sYjShLZXlOoRmwt9gy7VkuhPAjNpuNnDxXwf20LHPm\nYRC+pyyLQgSQWWxZJiCHj4TwA8WvW8tznC0KLlfZzfMsrFWWh48ygeKTtYYDJR643YaNoCCvrsPw\nKXa7rchtIAnU3CpcXraivxd5/9mMIvd99b1Z0V4zu8mvQ1kWhX3AQ8WWKeD9kjx40bo9tYEasbGR\nu80OzGoxMRFWh1BmAjW3ipJXZrajHRAEMPejn39b90rfpPx1w1/4b+V1r/TN8fwese6VvhnlGWtp\nVZTXbPgL/60MtCq06Ddv9m+72Ngm3lBKhQC/AjOA+cAw4EWgsdY6q0yeVAghhFfK7JyC1joX6A0M\nAc7gbjXcKgVBCCF8V5m1FIQQQvgfGeZCCCFEASkKQgghCkhREEIIUUCKghBCiAJSFIQQQhSQoiCE\nEKJAmc6ncCn8aQ4GpVR3YBbuQf5OAzO11m8ppWKARUBPIBmYorVeVOhx04HRuK8cXQI8obU2POsG\nA9OAmsAGYIzW+lT5ZXWWUqoW8DMwUmv9eSDkpZSqi/v/6xogBfdr9i9/z00p1Q14FWgKHPPEv9yf\n81JKdQJWaa3reu6XSS7l/ZlznrzqAvOAHkAu7lGkn9Ra51mRl0+1FPxpDgbPP+gaYLbWOgYYCLyo\nlPob8DaQBtQA7gRe9vwjoJR6GPdFfa2A5kB34EnPujbAG8AgoDpwElhcjmkVtxCILXR/Af6f12pg\nD1AVuBGYpJTqgh/nppSyA6uAF7XW0cC9wLtKqQb4aV5KqVHAF0ClQotNz6W8P3MukNd7wBGgDtAO\n6AhMsCovnyoKFJqDQWvt1Fovxp1oH4vjOp+GwKda6w8AtNY7gI1AN6AvMFFrnae13gIswz3MB8A9\nwByt9SlPRZ8OjPCsGwKs1lpv1VrnAE8DNymlapRTTgWUUvfhfgMe8dyPwM/zUkp1xv3Ge1Zr7dJa\n7wO64v5m7c+5xeD+UMj/oDGAHMCFH+allHoOeAT3N+D8ZWX1/9eTcvrMuUBelXAPEjrNk9cp3OPD\ndbMqL18rCn4zB4PW+iet9fD8+0qpqribfwB5WuvDhTfnbA7Fc9S4Bwo8Z53WOhFILLS+XCilmgJP\nAA9wdtzMK4Bcf84L6OCJY6ZS6rhSaj/uohCLH+fmec43gBVKqTzgG+Bh3IXCH/NaqLVuD2wttKwp\n5uZyxrNeUX6fOefk5SkEtxY7LHcrsPN8sVMOeflaUfDLORiUUtG4m2pbcLcWio/vVDiH4jlmAnbP\nAIKW56+Uyj9u+YjWOrnQqgj8OC+PWNyt0VNAfWAkMBeIxI9zU0rZPM85APdw9bfhPr9QBT/MS2t9\n8jyLwzE3lyzPY8stzwvkVYRSai7uD/TpnkXlnpevFQWv52Aob0qpxsD/gATcb8p0oPiM5oVzKJ5j\nOODwDCDoC/lPBHZorf9TbHkm/p0XuA+pnNFav6y1dmitfwA+ASbj37ndDnTSWq/y5PU58CkQj3/n\nVVhZ/f/5RJ5KqVCl1EdAL+AarfUZz6pyz8vXisI+zm2enq8Z5BOUUh2ATcB6rXV/z3G9A0CIUqpe\n4U05m0PxHJt5lp2zTilVHfcJ0X2Un4HAXUqpRKVUItAAWAHcjH/nBe7mc7Dnm3W+IGA7/p1bA6By\nsWUOYBv+nVdhZfW+svwzx3Po+RvcJ4S7aK3/KLS63PPytS6pXwOVlVIPcXYOhpq4z9b7FE93zfXA\nLK31zPzlWut0pdQaYLpSaizuXgODcfcgAHdPg3FKqQ2437jP4D5cA7Ac2KiUWoT7g2o68LnWumCi\nk7KmtW5e+L5S6nfgQa31eqVUW/w0L48vcX+DmqSUmgp0Bvrh/nbWCP/N7UvcPd+Ga63fVUpdizuv\nnvh3XgXK6n2llPKFz5xVwHFggNbaWWxdueflUy0F7V9zMIzCfSJvglIqzfOT6vmwuRcIAY4CHwJP\naa3zTy69jrsr64/AbuA7YDa4T157HrsYOAHU9jyPlQzOnmweix/npbXOBq7DXQxO4X7DPaK1/hE/\nzk1rvRu4A3hMKZWMu/vhMK31dvw4r/MwPRerP3OUUl1xd1DpBSR7PkNSlVIbrcpL5lMQQghRwKda\nCkIIIawlRUEIIUQBKQpCCCEKSFEQQghRQIqCEEKIAlIUhBBCFJCiIIQQooAUBSGEEAX+PyRIhgpf\nK5axAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b8e26f4fb10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.artist import Path\n",
    "from matplotlib.patches import PathPatch\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "xmax = snap25_counts.shape[0]\n",
    "xvalues = np.arange(0, xmax)\n",
    "\n",
    "ax.fill_between(xvalues, snap25_counts, y2=0, linewidth=0);\n",
    "ax.set(xlim=(0, xmax))\n",
    "\n",
    "\n",
    "for (left, right), n_junction_reads in junctions.items():\n",
    "    print(left, right)\n",
    "#     midpoint = (right - left)/2\n",
    "    curve_height = 10 # 3 * float(ymin) / 4\n",
    "    \n",
    "    \n",
    "    left_height = snap25_counts[left-1]\n",
    "    right_height = snap25_counts[right+1]\n",
    "    print('\\t', left_height, right_height)\n",
    "    vertices = [(left, left_height), \n",
    "                (left, left_height + curve_height), \n",
    "                (right, right_height + curve_height), \n",
    "                (right, right_height)]\n",
    "    print('\\t\\t', vertices)\n",
    "    codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]\n",
    "    midpoint = cubic_bezier(vertices, 0.5)\n",
    "    \n",
    "    if n_junction_reads:\n",
    "        plt.text(midpoint[0], midpoint[1], '{}'.format(n_junction_reads),\n",
    "             fontsize=6, ha='center', va='center', backgroundcolor='w')\n",
    "\n",
    "    path = Path(vertices, codes)\n",
    "    patch = PathPatch(path, ec='black', lw=np.log10(n_junction_reads + 1), fc='none')\n",
    "    ax.add_patch(patch)\n",
    "    \n",
    "#     xs, ys = zip(*vertices)\n",
    "#     ax.plot(xs, ys, 'x--', lw=2, color='black', ms=10)\n",
    "\n",
    "#     for i, (x, y) in enumerate(zip(xs, ys)):\n",
    "#         ax.text(x, y, 'P{}'.format(i))\n",
    "#     ax.text(0.15, 1.05, 'P1')\n",
    "#     ax.text(1.05, 0.85, 'P2')\n",
    "#     ax.text(0.85, -0.05, 'P3')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def junctionplot(left, right, n_junction_reads, read_counts, curve_height_multiplier=0.1, #curve_height=10, \n",
    "                 color=None,\n",
    "                 text_kws=dict(fontsize=6, horizontalalignment='center', verticalalignment='center', backgroundcolor='w'),\n",
    "                 linewidth_function=lambda x: np.log10(x + 1),\n",
    "                 patch_kws=dict(facecolor='none',)):\n",
    "    \"\"\"\n",
    "    \n",
    "    Uses the y-axis limits to determine the curve height so make sure to use this function AFTER \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    patch_kws : dict\n",
    "        To change the color of the line, use the \"edgecolor\" keyword, e.g. dict(edgecolor=\"black\")\n",
    "    \"\"\"\n",
    "    ymin, ymax = ax.get_ylim()\n",
    "    yrange = ymax - ymin\n",
    "    curve_height = yrange * curve_height_multiplier\n",
    "    \n",
    "    left_height = read_counts[left-1]\n",
    "    right_height = read_counts[right+1]\n",
    "#     print('\\t', left_height, right_height)\n",
    "    vertices = [(left, left_height), \n",
    "                (left, left_height + curve_height), \n",
    "                (right, right_height + curve_height), \n",
    "                (right, right_height)]\n",
    "#     print('\\t\\t', vertices)\n",
    "    codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]\n",
    "    midpoint = cubic_bezier(vertices, 0.5)\n",
    "    \n",
    "    if n_junction_reads:\n",
    "        plt.text(midpoint[0], midpoint[1], '{}'.format(n_junction_reads), **text_kws)\n",
    "\n",
    "    path = Path(vertices, codes)\n",
    "    \n",
    "    patch_kws['linewidth'] = linewidth_function(n_junction_reads)\n",
    "    if color is not None:\n",
    "        patch_kws['edgecolor'] = color\n",
    "    patch = PathPatch(path, **patch_kws)\n",
    "    return ax.add_patch(patch)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "color = 'steelblue'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.90848502,  1.90848502,  1.90848502, ...,  1.95424251,\n",
       "        1.94448267,  0.77815125])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snap25_counts_log = np.log10(snap25_counts+1)\n",
    "snap25_counts_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def coverageplot(counts, color, ax=None, **kwargs):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    xmax = counts.shape[0]\n",
    "    xvalues = np.arange(0, xmax)\n",
    "#     ax.set(xlim=(0, xmax))\n",
    "    \n",
    "    return ax.fill_between(xvalues, counts, y2=0, linewidth=0, color=color, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[<matplotlib.axis.YTick at 0x2b8e27320990>,\n",
       "  <matplotlib.axis.YTick at 0x2b8e26f66310>,\n",
       "  <matplotlib.axis.YTick at 0x2b8e273b1390>,\n",
       "  <matplotlib.axis.YTick at 0x2b8e273b1a10>,\n",
       "  <matplotlib.axis.YTick at 0x2b8e273bb1d0>,\n",
       "  <matplotlib.axis.YTick at 0x2b8e273bb950>],\n",
       " [<matplotlib.text.Text at 0x2b8e26fc8b10>,\n",
       "  <matplotlib.text.Text at 0x2b8e2734b310>,\n",
       "  <matplotlib.text.Text at 0x2b8e273b1890>,\n",
       "  <matplotlib.text.Text at 0x2b8e273bb050>,\n",
       "  <matplotlib.text.Text at 0x2b8e273bb7d0>,\n",
       "  <matplotlib.text.Text at 0x2b8e273bbf50>],\n",
       " <matplotlib.text.Text at 0x2b8e272c5410>]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAECCAYAAADdD/HDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WeYXFed5/Fvxc5VnaM6Sa0+rSwrSw5ywNkYjG1sk4cF\nBtae3R0YdpeZMR4GFliYgR0YYJZkYGewDTYOOGOwbGNbVrCSFY5S55xDVXflfVElubvVkjpUrv/n\nefpR162q2+eo7q3fveeee44hEAgghBBChIsx1gUQQgiRXCRYhBBChJUEixBCiLCSYBFCCBFWEixC\nCCHCSoJFCCFEWJljXYBICgQCgf5+B8nUpdpgMFBQkEWy1QuSt25Sr8STrHWbbb2KinIMC/k7SX3G\nYjAYMCZZDY3G5KwXJG/dpF6JJ1nrFq16Jdl/mxBCiFiTYBFCCBFWEixCCCHCSoJFCCFEWEmwCCGE\nCCsJFiGEEGElwSKEECKskvoGSSGEmA+TyXj2X4Nh9jdIer3+SBUpoUiwCCHEJGazEbs9E+Dsv7M1\nOOiQcEGCRcSQ2Ty1JdZoDI4iMZujRNl5hYhfEiwiJsxmI3l5WTM+N5ujRDkyFCJ+ycV7IYQQYSVn\nLEIIcRGPPfYYQ0NDjI+P09fXR3p6OrfddhtKqVgXLS5JsAghRIjH52do3HNOM20gEGBkZAS3280X\nv/hFdu7cSW9v7znBsuNoFx6PD1uGBVuGlZwMC7YMCxkWEwbDgkaiTygSLEKIlOJweWjuc9A7OhH8\nGZk4+/uQw82SUhs/+PTlU97T1tbGX//1X/Pd736Xt956i66uLj784Q+fs+6XDndwsmvknOUWkwFb\nhpViWzq1RTnUFmVTW5RDiS09KQNHgkUsmD8QwO314/H6cXl9uH2Tfvf68fj8uDx+3D4ffn8Ao9FA\nfnY6V007Knz22Wdpa2vjN7/5Ddu3b8dqtfLJT36SoqKic/7mkMONx+vDaDSQbjZhNRuTcgcVCzc8\n7kZ3DnOsY5hjncO09js4X5/DnHQLFTN0KsnPz+cXv/gFVquVn//851x22WUcPnyYFStWTHndDasX\n0VQ6ysi4m5FxT+jHzciEh/4xF/1jLo52DJ99fabVfDZkFof+LcxJS/ht2ZBMs6PNIDAwMIbPlzx1\nNJkM5OdnE656BQIBHC4vw+Mehp1uhpxuhp3us4/HPV7cXv+7Pz7flMcenw/PPMpRN8NRIcDzzz+P\n3+/n5ZdfxmKx8Ld/+7dkZ2ef87p7f/LalCNDs9FAVrqFrDQzWVZz8N9pP5lpZrLSLO8us5rJTjdj\nNZvmXP65CPdnFi/itV4DYy6OdQ6f/ekYdE55PifdQl1JDiX2DIpy0oM/tnQKs9PIsJov2GPxYi7U\nW3HC46N9wMHp3jEae0dp7B2jfdCBf9p/nS3DwsbaQrbWFVFfZscYxpCZ7We20Bkk5YwlSU14fAw5\n3Yw43QyNuxl2eiaFRvBx8F833ulb9ixZTAYsJhPZaUYsZiNWs5E0swmLadLvZiPW0GOr2YTJCP4A\n5Oekz7jOPXv2sG3bNj7+8Y9jNpt59dVXuemmm8553aL8LLxePz5/gAmPNxiOofrNpx6Zae+Gki3D\nQmFO8IumMPSFU5STTlaaOeGPJJPVsNPNa8e7eU130z4tSHIzrTSU22kos7OszE55XmZMPsd0i4kl\nJTaWlNjOLnN5fLT0O2jsHeX02bBx8scjnfzxSCd5WVY2Lyli65IiFhfnJMz2J8GSoAKBAENONx2D\nTjqHnXQOjtMx5KR7eJwhpxvXLO/xyEozU5xpxZ5hwZ5pJTfDin3S40yrORQKoR+T6WyILORIavrN\nkQBOp5Pc3Fyqqqp49NFHSU9P584775zx/Z+9Wp1zZOj3B3C6gyEz5cftweHy4py+fNLPyCxCKd1i\nojAnGDKF2ekU5qQFAygn+Lst3ZIwO34y8PkDHGwdYMfRLva3DOALHSAVZKexLBQkDeW5cX0dI81i\nYmmpjaWl74ZN/5iLXad6efNkL6d7R3n+YDvPH2ynKCedLXVFbK0rojI/K27rBNIUFve8Pj/dIxN0\nDDnpHHTSNTJOz6iLlt5RnG7fjO9JMxuDIZE5NSSmh4Y904rFFJtbmRbS3ADhv0HSHwgwHgqlIaeb\nvlEXfaMT9I256B2doH90gt5RFx7f+f+m1WwMnuVMCpuinHSqC7NZuaSI4SFnQm+L08WqKaxraJxX\ndBev6W6GQgcD6RYTW+uK2N5QypIFHtlHqilsPrqHx9l5qpedJ3tpHXCcXb64KIc7NlWzalHenOoa\nraYwCZY4MjbhoXFS+2vrgIOekfFz2mAheE2h1J5BWV4m5bkZlOVmBn/sGWSmxf+JaLwFy2wEAgFG\nJjz0jUwOnOC/Z0JowjNz2FtMRhblZ1JVkEV1QTZVBVlUFWQnxGd1PtEMlkAgwP6WAZ7Z38axzncv\nfteX2riyoZRNS4pIt4TnWlk8BctkbQMOdp7q5Y0TPfSMTACgymzcuamWhjL7rNYhwRIecRssTreX\npkkhcrp39OzGMpkt3UJpbgbluZmU52VSkZ/J8ppCrH4/ifzRJWKwXEwgEGDM5Q2GTOiMp3tknNYB\nB639Tpxu7znvKcpJoyoUNNWFwX+LcuK36WayaAVLc98Yv37zNIfbhwCwZ1i4XJVwRUMp5blzGyRy\ntqxWE3Z7JsPDTvxzuAYZjW3S5w/w+oluHt/TTO+oC4BVi/K4Y1MNS4pzLvheCZbwiJtgmfD40J3D\nHGkf4kjHEE19Y+cEQ1FOWqiPe7Cfe3VhNjnplimvideeOHOVjMFyPiaTgdy8LHRTH409Y7T0j9Hc\n76Clb4y+Mdc5r8+0mqgpykGV2qgvtVNXkkOGNf7ObCK9LQ46XPx2VxOv6W4CQGF2Gh/cXMumxYWY\nI9yEmwj7mdfn55VjXTyxt4XBUJPg5aqEj2xbQtZ5zoQlWMIjZsHi9vo43jXC0Y4hjrQPc7p39OzF\nRQieidSV2s72Xa8tysaWYb3oehNhg5+tmUY3nu1RYqKEClz4M3O4PLT0O2juc9DcP0ZLv4P2AceU\nnnpGA1QXZlNfakeFLvTmZaVFuxrniNS26Pb6eXp/K8/sb8Xl9ZNuMfG+dVVcv6oC6wydPiIhkfYz\nt9fHHw938sTbLThcXvKyrHxqez1rqvLPea0ES3hENVgcLi/7m/vZ09jPwdaBKT2zstLMLCu3s7w8\nl+UVuVTMs8tjIm3wc5WsdZtrvbw+P019YxzvGkF3DnO8a4TRCc+U15TY0lldmc/qqjyWl+eSFqbr\nC3MRic+rdcDBD186RuuAA6MBrlpWxgc2VmOfxUFXOCXitjjkdPPgqyfY29QPwFXLSvnYZXVTOuhI\nsIRHxINlyOlmb1Mfexr7OdI+dPasxGo2ng2R5RW5VBVkheVGp0Tc4GcrWeu20HoFAgE6h8c5HgqZ\n410jdA2Pn33eYjLQUJ7Lmsp81lTlUWrPiMo1mnB+XoFAgD+808FDO0/j8QVYUpzDZ66spyJ//s2l\nC5Go22IgEOCNEz386vVTOFxe6ktt/Nfrl58NZgmW8IhIsLi9fnad7mXH0S505/DZ4SEyrSYuqS5g\nw+JCVi/Ki8hRZKJu8LORrHWLRL36Ric42DrIgZYBDrcPTemNVmxL55LqArbWFS246+2FhKteLo+P\nH/3pGHsa+zEY4H3rqnj/uqqIX0e5kETfFntHJvjn59+hbcBJUU4an79xJZX5WRIsYRLWYGkfdPLy\nkU7+fLybMVewh489w8L62kI21hawrDxXLiouQLLWLdL18vr8HO8a4UDLAAdaB2gbePfO82JbOluW\nFLF1aTGVYT76D0e9Rsbd/PNzhznVM0p+lpV737MMNcuus5GUDNui0+3lhy8dY3/LAJlWE3936xoW\nl+RIsITBgoPF6/Oz63QffzrSebb/vMEA66oLuGpZGasr885OqRsNybDBn0+y1i3a9eobnWDX6T7e\nPNlDY+/Y2eWL8jPZVlfMtqXFFJ5nSJ25WGi9uobH+dYzh+gZmaC6IIu/uWllXHRKgOTZFv3+AL/8\n80n+eKQTW7qFf7h9LSsWF0uwLNC8g8XvD/DGyR5+t6f57P0lBdlpXNlQyvaGUvKzY7MDJMsGP5Nk\nrVss69U55GTnyV7eONlD51DwusyZA6PrV1WwrNw+76ayhdSrb3SCrzy+n0Gnm1WL8vgv1y2Lqy7V\nybQt+gMBfvyy5s/HeyjITuP/fPJSzD6fBMsCzDlY/IEAu0718bu9zWdHRW0os3Pz2kWsqcyP6tnJ\nTJJpg58uWesWD/UKBAI09zt440QPr+ouxiaCTbmV+Vlct6qcbXXFc74mON96jU14+McnD9Ax6GTj\n4kLuvaYhptdTZhIPn1k4+fwBvvfiEfY29VNfbuf+W1dj4PzfZRIsFzanYDnUOshDO0/T0h8ck2dp\niY3bN1azoiI3bu6ETrYNfrJkrVu81cvt9fHmyV5eONR+dlvPTjNz1fIyblq9iJwMy0XWEDSferm9\nPr7x+0Oc6B6hoczOf795VdTuTZmLePvMwsHt9fOVx/fR3O/gtg1V3L6h5ryvlWHzw2Dc7eXXb57m\n5aNdANQUZnPHxhrWVM1tgDchEoHVbGJ7QylXqBKOdQ7z4qEO9jT18ft9rbx0uIP3rq3k+lUVEenV\n+O9vnOZE9wiV+Vn89Q0r4jJUkpXVbOTea5fx94++zRN7W1hZkRexjhIpHyxHO4b48cua3lEX2Wlm\nPnZZHVvriiRQRNIzGAwsK89lWXkufaMTPPl2CzuOdfGbXU384XAHt2+o5gpVGrbm3wMtA/zpSCeZ\nVhNfuHHFeYcdEZFTWZDFp65p4IcvHOEnO47zv+/agCkCzfspe7jg9fn599dP8b+eOkjvqItLqvP5\n5l0b2La0WEJFpJzCnHT+0/Z6vnnnetbVFDDocPPTV07wpd/u5Z22wQWvf2zCw092HAfgY5fVhaVX\nmpifWzcGB6vsGh7njRM9EfkbcX/IoJTaAtwNFAF+rfVHF7pOt9fHv7xwhAOtg2RYTXx02xIuVyUS\nKCLlVeRn8fkbVqA7h3lo52lOdo/yzacPcePqCj64uXbe8/c8vLORIaebjbWFXLq0OMylFnNhMBi4\nY1MN//vpQzyxt5ltS4vDftaSCGcs9wEOrfWHgduUUgv6Hxh3e/n2s+9woHWQUnsGX79jPVc0lEqo\nCDGJKrPzwPvX8tmrFRlWE88dbOeB3+2jbdJkU7PVMzLOa8e7STMb+cTldbKvxYE1VXnUleTQPTIR\nkbOWmJyxKKWWAy9qrRdNWnYfcA+ggXXAVq31uNb6I6Hn0wCn1nreXTQcLi/ffvYQJ7tHWZSfyZdu\nWY09M7qD2wmRKAwGA5fVl1BfauPf/qQ53jXC/Y/t456ttVy7ohwu0F11st/va8XnD3D9qgrZ32LA\n6/NzsHWQtgEHAQOkp1sx+nxc2VDGye5RXj/ezeWqBAh2Sw9H8Ec9WJRSvwDWACWTltmB7wB2rfW4\nUuoZ4BHg1klvfRD42Hz/rt8f4J+fe4eT3aPUFmXz329edc5cJwAenx+z0SBHVUKEFNsy+Ltb1/DU\nvhYe39PMr/58it6RCT562ZKLvrdvdIJXdTdWs5Gb1yy66OtF+PSPuXj9RDcer5+11QWsqcrHajGS\nn59NS8cgz+5vw2o2cqRjiLEJD9npFkYmPGEZSTrqwaK1/oRSahvwyqTF9wEjWuszQ7Y+DXzzzJNK\nqS8CPwLmPdjR84faOd4V7Ob4pVtW4/H5+cM7HThdU2f1M5kMU+ZNOTPCZIAARoOBYlsGi/IzKbVn\nxN1NXUJEislo4Lb1wXu6/unZwzx3sJ0A8N9uXXPB971+ogefP8A1y8vkbCVK3F4fLx3uxAC8Z0U5\nmTOMaJCdbuH2jTXsbeyjZcDJ2039XNFQSufQODlps7uP6ULi5eJ9AzA+6XEnkAOglPpb4O8BJ+AH\n5nzlr2PIyW93NWE0wKX1xTx7sA1bhoXNi4vmtLH7/AG6R8ZpH3Syr3kAn99PAMi0mqkvtVFdkB3z\nO/OFiKT6Ujtfeu8qvvn0IZ4/2I7FaubuTdXnff2exj4AttbJBfto2Nfcz4muEa5ZUU7BLIadunlt\nJT/6k2Zf8wCXqxKOtA9RlpvxbnPSPMVLsBQB7kmPnYQacLXWXwe+Pt8VG4wGfrrjOB6fn/pSG6sq\n86grtc1rXRaguiib6qLsKcsdLg/HOoZ5+mArgQBkWEyocju1RdmYjOE9qzkTXGazEaMxOe4IPiNZ\n65Zs9aovt/Pl29bwtScO8vs9zeD385FLz20W6xudoLF3jLwsKw2L7GGZjyhaEu0z8/j8PLm3hSXF\nOXz4Ak2U0+u1ujo4y2TrwBhGk4E/HulgaamNuqqCBZUnXoJliKlXAu2A7zyvnZPjvcGZ+GwZFh64\nawOFtoxwrHaKXKCiJPfsY4fLw5HWQV480o3X56c0N5PN9cWzmnp4tmy2zLCtK94ka92SqV5rc7P4\n9scy+cIv3+T3b7eyfWUFq6qnfhm9eqIXgG0NpeTnZc+0mriXCJ9Zx4CDp/a2cvuWxRTN8vvtTL3s\n9kzSzEZ6RibIyslg2Omhc9S14DLFS7DsBd4z6XE1wbOWBXv5YDsAt22owuz3MzQ09+6S87G0MIul\nhcFLQl3D4zz9ViMOlwd7ppUNtYXzHh3ZaDRgs2UyMnLxeeETTbLWLVnrVZBh5i+vW853fn+Qf3ry\nAN+6Z8PZYWACgQBHWoJT5C4tyo7afhcuifKZvXWyl67hce7cUIlpFt9vM9Wr1J5Bc7+D46HPa6fu\n4sNXLF1QuWIVLNPbh34IfF0pVaS17gXex9SL+/Pi8/t561QPBmDL4mI8Hv9F3xMJBZlpoe6ZwamM\nd5/qY9DhIsNqZm11PuW5sz8qMpmCJ3Zerz9pBsc7I1nrlsz1um7NIl7a38rB1kEefrORD21dDASb\nwdpDE44V56THbN+br3j/zLw+P7/f18qSkhxuWr0Ivw/8vov/H89Ur/K8TJr7HbT3Bz+vpr4xrv/q\nM/YX7r95eL7li3q3JqXUD4BfAWal1C6l1Ge11g6C11GeUUo9SLB16QML/VuHWgZwunwsLbXNesTW\nSMvNtHL18jJu31jD9oZSGntGeWx3E0+93UJj7yhJPtq0SDIGg4FPXVlPusXEcwfb6Boax+8PMDzu\noXNoHAPBI2IRPm6vn0feauTS+mJWV+YveH1neo25vb6z6wdWL2SdsehufO95lj8APBDOv3W4JTjG\n0YbawnCuNmyy0sxcWh/sf+H2+jjUNsTepmay0yxsXVoUlv7kQkRakS2d96wo5+n9rexq7GXLkiIc\nEx5GJzwUZqdhNYd/lORUNe728ujuZm5Zuyhss20uLs6htd+BberB94JWHi/XWCLC4fIAwXnp453V\nbGJ9TQHrawoYGXez82Qvw+MeGsrsrFiUm1A9akTq2VhbwNP7W9nb2E99iY2x0P1h2TPchCzmZ3TC\nwxN7W/jAhiqywnCvyRnbG0q5pDqfruHJd3xIsJzXVlXK6JiL2qKcWBdlTmwZVq5bVUEgEEB3DvO7\nPc1kWc1sXVpMQU58zAkuxGS1xTnkZ1k51TPKGyd7qczLJCvNLPOthMnAmItnD7Zx56Ya0iMwT47P\nH6Bvam8wmejrfFZU5oHHS2kEuhhHg8FgoKE8l4byXEYnPOw82cvIhJsN9aXU5Mqw4yJ+GA0G1tcW\n8od3OtjfMkBtaLK8stzE3PfiSffwOH880sldCxhdejZO9YyGbV1JHSwuj48JT1huh4m5nHQL164s\nx2iEToeHx3Y1k24xsbWuKGxtrUIsxJkQcXl8jIx76Bx2UmqXA6CFaB1wsPNkD3dtro3IhFyTvaa7\nw7aupA6WgTEX3ll0wUskBoOBlVX5lGdbGXa4efNkL0NOF0tL7KyuzJMhZUTMFGSns6Q4h+7hcRxu\nLye7R2XulQU42T3CO22D3L6xJirXWJ1u78VfNEtJHSwer59Bh/viL0xQ2aGzmEAgwInuEX6398xZ\nTPGsxgkSIpzOdD6591c7cbg89I5OxLpICet0zyjHOoZ537qqhBxpPamDBcDlTY6msAsxGAzUl9qp\nL7XjcHnZebKH/jEXq6vyUaW2hNwwReIKBAKMu31TRwkXs9Yx6ORAywDvX5+YoQJJHiyBQODMzT4p\nIyvNzDUrgmcxB1oG+c2uJuqKc7ikpkC6LIuoOdE1EusiJKT+seD8NR/cXJOwoQKJMTXxvPkDAbqn\n9s1OGQaDgbXV+dy1uRZ7ppVHdzXx6rGupLvmJOLToDN5m6AjZXTcw/MH27ljY3XCHwQm9RlL56CT\nI+1DsS5GzNWV2KgrsdEx6OTJt1vISbdwmSqZcQIgIUT0TXh8PPl2C3dsqkmKCQST+ptl0OHCI0fo\nZ5XnZXL7xhoGHS5eeqcDDLBdlcrMfkLEkMfn59HdTbx/fVVEbn6MhaQOFleCjagaLXlZady6rgqH\ny8truguHy8ul9SUyWKAQUeb3B3hsdxM3rV6UVMPfJHWwiAvLSjNzw+pFuL1+Xj/RzavHuti4uDDh\nhsARIhEFAgEe39vMlQ1l856fKV5JsAisZiNXLSvD7w+wu7GPt071saLCzspFeQndM0WIePbsgTbW\n1RRQmoTD3iT+VSIRNkajgc1Lirhrcw1Go4Hf7Gpi58keuR9BiDB7+UgntUU5Sds6IGcs4hwGg4EV\nFXmsqMijqW+M3+1pojA7nctUSUQHwRMiFRxsHSDdamJ5RW6si3IhMrqxiJyawmxqCrPpGh7nib3N\nlNgyuLS+JOID4gmRjDqGnDT1jXHrJVWxLkpEyeGnmJVSewZ3bqqltjiHR3cHm8j8Mo2yELPmcHnY\ncbSLW9ZUxrooESfBIuakMj+LuzbXUmzL4DdvNbK3qZ+ABIwQF+TzB3h8bwvvX1+VEiOQS7CIeVlc\nnMPdWxaTlWbm4Z2NvNM2KAEjxHn8fl8L160sT5nRLiRYxII0lNm5e0st/kCAh3c2crxrONZFEiKu\nvKa7qS+1U5ygM9nOhwSLWDCDwcDqynzu3lLLyLiHh3c20tQ3FutiCRFzxzqH8QcC8d4DLOwkWETY\nGAwGNtQW8sFNNXQMOnnkrUbaB52xLpYQMdE7MsGR9iG2N5TGuihRlxoNfiKqjEYD25YW4/X5ef1E\nD2+c6GF7Q0lKNQWIoFS96jbh8fHCO+3cvbk21kWJCQkWETFmk5HtDaW4vX5e010MObu5alnyjYsk\nUoPZPLsGHn8gwGtHenjfuqqkGAJ/PiRYRMRZzUauWVHOhMfHjqNdTHi8XL28DFvGzMP1m0xGDIbZ\nHet6U2yGUBEbZrORvLysWb/+Q9vrGRx0JMz2Ge4OnUkfLNIDNn6kW0zcsLoCh8vLy0c7IQDXrCgj\nY1oXTLs9c9brTKSdV4hUkfTBIuJPVpqZW9ZWMjzu5tkDbRTlpHNZfQkmU/LfOCZEKkjNBkARF+wZ\nVm7fWENNUTYPv9Uo00gLkSQkWETMVRVkc8+WWpxuL0/saox1ccRCSfNzypOmMBEXztwDkzOH6ytC\nxNITTzxBZ2cnfX193H///bEuTrjJsPkiecw030uS78AiQW3bto3i4mK+9rWvxboocUeawkTc27Zt\nG5/73OdkmmQRV/Lz8/ne977HHXfcEeuixB0JFhH3ZAcW8eirX/0qIyMj7NixA5/PF+vixBVpChNx\n76tf/SoWi4UdO3awdOlSTCZTrIskBF/5yldiXYS4JcEi4t6FdmCfX7ogCRFvpClMJLSXDndwsHUg\n1sUQQkwiwSIS2vWrKgB46M3TMkS/EHFCmsJEwltdmc+KijxePdbFrlO9XLuynOx0S6yLJZKI1+vn\nzSMdQHDW1Nm+J1VJsIi4MzzsxD/Laydndl6T0cBVy8twuDw8f7CdUnsG25YWSxdlERYj4272N/Vz\nx8aalA6M2ZKmMBF3fD4/Xu/sfqbLSrNw+8YaSuwZ/PrN03QMSfOYWJhAIMAz+9u4ZW1lrIuSMCRY\nRFKqK7Fx95bFHG4b4un9rbjlKFPM0+snelhfW0C6Rbq5z5YEi0haJqOBa1eWs7WuiN/taeZQ62Cs\niyQSTO/IBIMON/Wls7uuIoIkWETSK8hO5+4ttXj9fh55q5EhpzvWRUoZiTzRnj8Q4PmDbdy4uiLW\nRUk4cvFepIxLqgtYVp7Li4faycmwsF2VYjTKxX0xs6f3NHPlsrKUnbd+IeR/TKSUdIuJW9dVUVds\n46Gdp2nqG4t1kUQcaukfw2AwUFkw+3nuk8yCjrgkWERKqizI4p6ti2nuG+PJvS2Mu72xLpKIEx6f\nn1eOdnPz+qpYFyVhSVOYSFlGg4HtDaWMjLt5en8bNUXZbKgpkHtfUtxzB9u4cU0FRtkO5m3WwaKU\nsgEPAD8FNPDvwJ3A28BdWuumSBRQiEizZVi5c1MNR9qHeHhnI9euLKcwJz3WxRIxoDuHKcxOl89/\ngebSFPavwPWAH/gQ8D7gI0BL6DkhEtryilzu2FTD7tN9vHioHa9P7n1JJeNuL/uaB9haVxTrokRd\nuDvvzSVYbgY+orXWwAeAZ7XWjwB/B2wPc7mEiAmLyciNaxaxuiqf3+xq4njncKyLJKLkmQNt3Lx2\nkTSFhsFcgsUMjCmlrMB7gOdCy7OAiXAXTIhYKrVncM+WWgYcbh7d3YTD5Yl1kUQEHWgZoK7YRo4M\nXhoWc7l4/2fgu8AwYAGeVEqtJ9gM9ocIlE2ImDIYDGypK8LpyuPZg23UFdtYW50f62KJMJvw+NCd\nw3xwc22si5I05nLG8mnAB6wE/kJr3Q+8H2gD7otA2YSIC5lpZu7YWAPAb3c14XRJ1+QLSbSb7Z8/\n2MYNcnd9WM36jEVr3UEwSCYvuz/sJRIiTq2tzmdpaQ5P729lWbmdVZVy9pLoTnWPUJSTji3DGuui\nJJULBotS6jOzXZHW+scLL44Q8S0rzcIHN9eyt7GPx3Y3cdOaRWRY5XawROT1+XnrdB/3bJEmsHC7\n2B7xpWmPqwAXcBrwAHVAOrAPkGARKWN9bSH1pXae2tfKqkV5LK/IjXWRxBy9dLiT96wok15gEXDB\nayxa69pdwNCPAAAVO0lEQVQzPwSD42lgkdZ6pdb6EqACeBx4NfJFFSK+5GRYuGtzLWMTHh7f08yE\nxxfrIolZ6hxyYjIaKLZlxLooSWkuF++/CPxPrfXAmQVa6xHgywQv7AuRkjYtKeLq5WU8sbcZLfe9\nxL1AIMDLR7u4enlZrIuStOYSLONA/QzL1xPsgixEyrJnWrlrcy0DDhdP7m3B7ZWzl3j12vFuttYV\nYZIpEyJmLlcd/wX4hVLq+8ABgsMqbwY+B3whAmUTIqEYDAa21hUz6HDx2O5mNi4upK7EFutiiUmG\nnG6GHG6uUKWxLkq8i86w+VrrbxFsDrsB+CXwIHAp8Emt9f9dSCGESCZ5WWncvaWWruFxfr+vFY+M\nORY3XjjUzvWr5J6VSJtTP0mt9U8Jjm4shLgAg8HAZfUl9I+5eGRnIzduqKEgzRTrYqW0t5v6WVGR\nS5pFPodIm8uw+SaCw+SvAM58MgYgDVivtZaBKIWYpiA7jQ9vW8zetiHe6B3jupXlMtVtDDjdXk71\njHLnpppYFyUlzOWM5V+BjwN7ga3AG8ASoBT4QfiLJkRyMBgMXLemEt3UxyNvNXJZfQnVhdmxLlZK\neeFQuwzbEkVzOXS6A/io1vpy4BTwWaAa+C0g4yEIcRFFtnTu2bqY0z2jvHioHX8g0UbVSkxNfWMU\n5aTLyMVRNJdgsQG7Qr8fAjZqrb3A1wnO1SKEuAijwcBVy8tYVpHLQ2+eZnjcHesiJbVAIMAbJ3rY\ntrQ41kVJKXMJlhagIfT7UWBd6HcXkBfOQgmR7Crzs7h9Yw0vHurgSPtQrIuTtN461cemxYUyf32U\nzSVYfgw8pJS6BXgC+JRS6u+AHxKc914IMQfpFhN3bqphZNzNswfa8PulaSycJjw+WgfG5F6iGJjL\nfSzfBv4rMKa13ktwDpa7ADfwqcgUT4jkt6WumHXVBfx652kGxlyxLk7SeOlwB+9ZUR7rYqSkud7H\n8h8ASikz8Aut9YMRKZUQKaY0N4MPbqrlmf2tLC7OYU2VzPWyEH2jE5iNRvKy0mJdlJQ0pw71SqnP\nKqVOAE6gVin1Y6XU15RS0oApxAJZzUZu21CNx+fnybdb8Mod+/P2xyOdXLNCBpmMlVkHi1LqvwD3\nA98CzszN+ieC3Y6/Ev6iCZGaNtQWsq2umId3NtIzMh7r4iScI+1D1JXYsMiNqDEzl//5zwF/qbX+\nCeAD0Fo/DHwM+ET4iyZE6iqypXP3llp2nupld2NfrIuTMHz+APtbBlhXLU2JsTSXYKkGDs+w/CRQ\nGJ7iCCHOMJuM3HpJFWlmI7/b04zbK01jF/Oq7mK7KpFZIefsnB6J0RndGNgP3DZDST5LcGriuCQd\nOEWiW12Zz1XLSvnNrkY6h5yxLk7cGpvwMOx0U5GfFeuipLy59Ar7PPC8UuoqggNPfkUp1QAsIziU\nvhAiQvKy0vjQlsU8f6idxt4xuZN8Bn94p4PrVkn34ngwlzOWzwO3ErwZ8kkgB3gJaNBavx6Bsgkh\nJjEaDdy0ZhF5WVZ+u6uRCU98zlIZiEE7QeuAg7wsK1lpMh5YPJjLGcs1BOe8fyBShZmJUqqCYOeA\n+7TW0n9QpLxl5blU5GXy2O4mrlClVBZI08+fdTd3ba6NdTFEyFyC5TvAg0qp7wKNwJR+kFrr4+Es\n2KT1tgP/Sykl0x8LEWLLsHLP1sW8eKiDrpFxNtambv+Zw+2DLK/IxShz2MeNuQTLV0P/Xj5pWYBg\n74EA707+dVFKqeXAi1rrRZOW3QfcA2iCA1xu1VpLJ34hzsNoMHDD6gr2NvXz7IE2blhdkXKDLQYC\nAQ62DHL3FjlbiSdzCZawfHJKqV8Aa4CSScvsBM+I7FrrcaXUM8AjBK/pCCEuYH1NAa39Dh7Z2cht\nG6pJT6Gpd/e3DLC2Ol+6F8eZWQeL1ro5HH9Qa/0JpdQ24JVJi+8DRiadoTwNfHPaW2XLEeI8Kguy\nuGVtJY/uauKG1RUU5qTHukgR5w8E0J0jcrYSh+Y0CGUENTD1mk0nwV5nKKUKCE4mlqWU+rLW+h/n\nsmIDYLEkz9AOZ9qRzWYjRmNy3aWTrHWLVr3yLWl89PIlPL6nmVWVeTSU50bsb8H56zX5CNBkMkRs\n/9t5soetS4sisv5U2xbN5vD+H8ZLsBQRHH7/DCeh7VNr3Q/8ZehnzoxGA7m5yddrxmbLjHURIiZZ\n6xaten36+pW8sL+V3c2DXLtm0cXfsEDT6zW5WSo7Oz0i+5/PH6Bj2MUNGyJ7tpIq26LbkJzBMsTU\nAx07ofHIFsrvDzA05AjHquKC0WjAZstkZMSZdBNDJWvdYlGvzTX5HG4b5GcvHub9G6owGSNzVD9T\nvQKBd38fG5uIyP736rEuNlTnRWzfTrVtcWR0Iqx/J16CZS/wnkmPqwmetSxYAPB4kmeMJZMpmL9e\nrx+fL3k2eEjeusWqXvUldnIzrPzqtVO8f10V2enhvXnwfPWaXEOvNxD2/c/r89Pa52DrkuKI7dup\nti16wzwOXawuPkz/uz8E7EqpotDj9zH14r4QYh6KbRncvqGap/a10j6YHOOMvXa8m8tVycVfKGIm\n6sGilPoB8CvArJTapZT6rNbaQfAC/TNKqQeBXOAD0S6bEMkow2rm7i217GvqZ3/zQKyLsyBur4+B\nMRdlucl57SNZRL0pTGt973mWPwBEdbgYIVKF0WDglksqefNkD68c62J7Q2msizQvrxzrTtiyJ5io\nDZsvhEhwW+uKKchO4+n9rVMusieCcbeXsQlPStyjk+gkWIRIMSsX5bFqUR6P7m7C60ucji07jnZx\n5TI5W0kEEixCpKDqwmyubCjjkbfid/j9ySY8Pia8PvKy0mJdFDELEixCpKgiWzq3rqvit7uaGB53\nX/wNMfTGiR4ulcnNEoYEixApLCfdwl2ba3hmfxtdw/E5mLg/EKB/bIJiW0asiyJmSYJFiBRnNZu4\ne3Mtrx/vpqlvLNbFOceBlgHWVOXHuhhiDiRYhBAYjQY+sKGagy0DnOoZjXVxpjjeNcLSElusiyHm\nQIJFCAEEB4987yWVHO0Y4njXcKyLA0Bz3xjVBVky30qCkWARQpxlMBi4ZW0lp3tGOdoxFOvisLux\nj42Liy7+QhFXJFiEEOe4YfUi2gacHGodjFkZBh0uctItmGQu+4gL972yyR8sCXZ3sRDx4tqV5fSO\nTnCgJTbji/35eA+X1ctgk4ko+YNFCDFvVy8vY8jp5u2m/tm/KQzHci6PD6/fT1ZavMzsIeZCgkUI\ncUHbG0oZGXfzTlv0msXkhsjEJsEihLioK5eV0drv4FT3SMT/lj8QoHdUbohMZBIsQohZuWF1BQfb\nBmkfiOxU3wdbBuWGyNiTYfOFEJFnMBh437oqXj/RQ98s50gPzOOCy8meEepL5YbIRCbBIoSYNaMh\neIf+i+90MOQM/8CVTreXdItJbohMcBIsQog5MZuM3LGxmqfebsEx4Qnrut9u6mdDTWFY1ymiT4JF\nCDFnVrOJ2zfW8B+vncAfxnvFuofHKc2Vi/aJToJFCDEvWWlmblhbyTP72sKyvpFxN9nplrCsS8SW\nBIsQYt6qinKoyM9kd2Pfgte1u7GPDbXSDJYMJFiEEAuyrqaAvpEJWvsX1g150OGmIFumHk4GEixC\niAW7fnUFrx3vZmyeF/P7xybIl/nsk4YEixBiwYwGA7etr+LJt1sIzONi/u7T/WyUZrCkIcEihAiL\nDKuZy1UJXv/cg8Xh8pCTIRfuk4UEixAibKoKsuc8U0XnkJMSu3QxTiYSLEKIsLKY5nbX/L7mAdbX\nFESoNCIWJFiEEGE11+FYXF4fGVaZdyWZSLAIISJm+CLjifkDgYUNoyvikgSLECJiDrZeeHKw9kEn\nFflZUSqNiBYJFiFExORnpdF6gflbTnSNsLREhshPNhIsQoiIaSi38+aJnvM+PzzuJjfTGsUSiWiQ\nYBFCRIzJaKDEnkHHoDPWRRFRJMEihIioS5cW8/qJ7nOWj4y7yZHRjJNS0gdL+GaKEELMh9lkJCfd\nyuj41HHEmvrGqC3MjlGpRCQlfbAIIaJrpoO51VV5HGqb2kOsd3SCIlt6dAolokqCRQgRcWX2DLqG\nx6csGx33SFNYkpJgEUJE3Jm78aePfDzXu/RFYpBgEUJERWV+1gXvaRHJQ4JFCBEVFfmZdA292xwm\nZyvJS4JFCBEVGRYT4x5frIshokCCRQgRFfYMKxsmDY9/WX1xDEsjIkmCRQgRFUajYUrz17hbzl6S\nlQSLECJq5Ibl1CDBIoQQIqwkWIQQQoSVBIsQIoqkMSwVSLAIIaJm8o33EjHJS4JFCBE1gfM+EMlE\ngkUIIURYSbAIIaJHzlJSggSLEEKIsJJgEUJETUBOWVKCBIsQImoCU36XkElWEixCiJgISK4kLQkW\nIUT0SJikBAkWIYQQYSXBIoSImOnNXXLCkhokWIQQQoRV0geLXCAUIn4EZIeMS+H+WJI+WIQQ0TXb\nLymJmOQlwSKEiBgJj9QkwSKEECKsJFiEEEKElQSLEEKIsJJgEUJEzrQr+VNmkJQeYklLgkUIETUS\nJalBgkUIIURYSbAIISLm3DMUOWdJBRIsQgghwkqCRQghRFhJsAghhAgrCRYhRMRMv6IytbtxVIsi\nokiCRQgRE5IryUuCRQghRFhJsAghhAgrCRYhRORMn5pY2r9SgjnWBbgYpVQZ8P+APwPlWuvPxLhI\nQgghLiARzli+B/xJa/0PwK1KKUuMyyOEmKcA0i0sFcTkjEUptRx4UWu9aNKy+4B7AA2sA7ZqrceB\npcBDoZc5gY3AG9EtsRAi3CRWklfUz1iUUr8A/gMombTMDnwHeI/W+pNAO/DwpDJ6Qr+bAG/UCiuE\nWJCAxEdKinqwaK0/Adw7bfF9wEjoDAXgaeDK0O9HABX6PQvYF+EiCiGEWIB4uXjfAIxPetwJ5IR+\n/yLwM6VUNvCC1toz/c1CCCHiR7wESxHgnvTYCRgAtNatwHXzXbHBABZLIvRRmB2j0QCA2WzEaEyu\nZoZkrVuq1ctgePc1ZpNxyv5nNr/7ePLv8SbVPjOzxXC+t8xLvATLEKEgCbEDvoWu9CcvHS0FinJz\ns95Z6Lrijc2WGesiREyy1i1V6uVwedcSvB7Kd547fPqF+28ePPPcXd/fkfbC/Te7AP7++zuyX7j/\n5rFolnWuUuUzu+v7O9KAlZMWnV7I+g2xmHdaKbUNeEVrbQk9/iLwP7TWhaHHfwN8WWtti3rhhBBC\nLEiszkOn/90fAnalVFHo8fuAV6JbJCGEEOEQ9TMWpdQPgBuBWmA38HOt9b8ppb4SWn4Y2ACskwv1\nQgiReGLSFCaEECJ5xWeXDCGEEAlLgkUIIURYSbAIIYQIKwkWIYQQYSXBIoQQIqzi5c77sFJKXQL8\nG7ACOA58Tmv9VmxLNTOl1GXAPxEcL60X+LbW+sdKqVzg58DVBEcm+Eet9c8nve8bwH8ieIfzr4DP\na60DoefuAb4GFAMvA5/SWvdEr1bvUkqVAAeBv9BaP5sM9VJKVRDcvq4Ahgl+Zt9P9LqFblz+F6Ae\n6AiV/6FErpdSahPwuNa6IvQ4InWJ9nfODPWqAP4VuJzg8FiPAl84c8tGtOuVdGcsSqk04CngZwSH\nhvk+8JRSKu7GZght5E8C39Va5wIfBL6ulLoG+AkwSnActTuBb4U2pjNz19xIcAiGZcBlwBdCz60G\nfgTcBRQC3cCDUazWdD8D8ic9/imJX68nCN5vlQdcDzyglNpCAtdNKWUEHge+rrW2A58GfqmUqiJB\n66WU+iTwAjB5csCw1yXa3znnqde/A61AGbCW4LxV98eqXkkXLMBVgE9r/WOttU9r/SDB/6ybYlyu\nmVQDT2utHwHQWu8DdgDbCI4+8GWttUdrvRv4NfCx0Ps+AvwfrXVP6MjiG8AnQs99CHhCa71Ha+0C\n/gdww6RRDaJGKfWXBHfi1tDjLBK8XkqpzQR33i9prf1a66PAVoJH+Ilct1yCXyxnvqwCgAvwk4D1\nUkr9LfBXBI/EzyyL1PZ3NVH6zjlPvSzAGPC1UL16CM55tS1W9UrGYGkgOIfLZDq0PK5orQ9orT9+\n5rFSKo/gqSyAR2vdPPnlvFuH6XXUvDtnzZTntNYDwMCk56NCKVUPfB74HO8OMLoUcCdyvQjObnoE\n+LZSqlMpdYxgsOSTwHUL/c0fAQ8rpTwEh1S6j2DYJGK9fqa1vgTYM2lZPeGtS3/oeUX0vnPOqVco\nTN47rYnxvcD+mcpOFOqVjMGSRXDY/cmcQNw1hU0WmkXzKYLD3Oxg6vw0MLUO0+voBIxKKesMz01/\nb8Qppc604/6V1npo0lNZJHC9QvIJnhX3AJXAXwDfA7JJ4LoppQyhv3k7kAHcSvB6i40ErJfWunuG\nxZmEty7jofdGrZ7nqdcUSqnvEQyFb4QWRb1eyRgsToI7xmSZBE8V45JSqhZ4HegjuGOPAenTXja5\nDtPrmAl4tdbuGZ6b/t5o+DKwT2v94rTlThK7XhBsHurXWn9La+3VWr8J/A74Coldtw8Am7TWj4fq\n9SzBmVz/gcSu12SR2v7iop5KqXSl1KPAtcAVWuv+0FNRr1cyBstRzj3VnumULi4opdYBO4HntNa3\nhdo5TwBWpdSiyS/l3TpMr2NDaNk5zymlCgleZD5K9HwQuFspNaCUGgCqgIeBm0nsekGwKcAcOsI/\nwwS8TWLXrQpIm7bMC+wlses1WaT2q5h/54Sa0V8heJF9i9a6ZdLTUa9XMnY3/hOQppS6F/i/BC/M\nFRPsRRFXQl1xnwP+SWv97TPLtdZjSqkngW8opT5DsDfHPQR7dkCwB8gXlVIvE9z5/yfBpieAh4Ad\nSqmfE/yy+wbwrNb67GRLkaa1Xjb5sVKqEfjPWuvnlFJrSNB6hfyB4JHcA0qprwKbgfcTPEqsIXHr\n9geCPRI/rrX+pVJqO8F6XU1i1+usSO1XSql4+M55nOCU7rdrradPkhj1eiXdGUvo9O5Ggr0d+oF7\ngfdqrae3rcaDTxK8OHq/Umo09DMS+sL6NGAF2oDfAn+jtT5zwe6HBLsp7wLeAV4DvgvBDgGh9z4I\ndAGlob8TSwHevYD/GRK4XlrrCeBKgoHSQ3Cn/Sut9S4SuG5a63eAO4D/ppQaIti19GNa67dJ4HrN\nIOx1ifV3jlJqK8FOP9cCQ6HvkBGl1I5Y1UuGzRdCCBFWSXfGIoQQIrYkWIQQQoSVBIsQQoiwkmAR\nQggRVhIsQgghwkqCRQghRFhJsAghhAgrCRYhhBBh9f8BrjkAfKp1/LIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b8e272591d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "xmax = snap25_counts.shape[0]\n",
    "xvalues = np.arange(0, xmax)\n",
    "\n",
    "ax.fill_between(xvalues, snap25_counts_log, y2=0, linewidth=0, color=color);\n",
    "ax.set(xlim=(0, xmax))\n",
    "\n",
    "ymin, ymax = ax.get_ylim()\n",
    "yrange = ymax - ymin\n",
    "\n",
    "for (left, right), n_junction_reads in junctions.items():\n",
    "    junctionplot(left, right, n_junction_reads, snap25_counts_log, color=color, curve_height_multiplier=0.2)\n",
    "yticks = map(int, ax.get_yticks())\n",
    "yticklabels = ('$10^{}$'.format(y) for y in yticks)\n",
    "ax.set(yticks=yticks, yticklabels=yticklabels, ylabel='reads')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# class WasabiPlotter(object):\n",
    "    \n",
    "#     def __init__(bam_filename, chrom, start, stop, log_base, ax):\n",
    "#         self.\n",
    "\n",
    "def wasabiplot(bam_filename, chrom, start, stop, strand, log_base=10, ax=None):\n",
    "    bam = HTSeq.BAM_Reader(bam_filename)\n",
    "    counts = count_region_reads(bam, bam_filename, chrom, start, stop, strand)\n",
    "    counts = np.log(counts)/np.log(log_base)\n",
    "    \n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "        \n",
    "    coverageplot(counts, color, ax)\n",
    "\n",
    "    ymin, ymax = ax.get_ylim()\n",
    "    yrange = ymax - ymin\n",
    "\n",
    "    junctions = count_junctions(bam, chrom, start, stop, strand)\n",
    "    \n",
    "    for (left, right), n_junction_reads in junctions.items():\n",
    "        junctionplot(left, right, n_junction_reads, counts, color=color, curve_height_multiplier=0.2)\n",
    "    yticks = map(int, ax.get_yticks())\n",
    "    yticklabels = ('$10^{}$'.format(y) for y in yticks)\n",
    "    ax.set(yticks=yticks, yticklabels=yticklabels, ylabel='reads')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bam_filename2 = [x for x in bam_filenames if 'M2_03' in x][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_03_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam\r\n",
      "/home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_03_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam.bai\r\n"
     ]
    }
   ],
   "source": [
    "! ls $bam_filename2*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start: 10265372\n",
      "Skipping read with CIGAR string I (a base in the read was inserted): HWI-ST1001:211:C1PUCACXX:8:2302:5344:69707 (from /home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_03_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam)\n",
      "Skipping read with CIGAR string D (a base in the read was deleted): HWI-ST1001:211:C1PUCACXX:8:2312:7987:8949 (from /home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_03_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam)\n",
      "Skipping read with CIGAR string D (a base in the read was deleted): HWI-ST1001:211:C1PUCACXX:8:2205:8970:2276 (from /home/obotvinnik/projects/singlecell_pnms/analysis/bams_from_aws/M2_03_R1.fastq.gz.polyATrim.adapterTrim.rmRep.sorted.rg.bam)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/obotvinnik/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:9: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEKCAYAAADNSVhkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUJGWd5vFvZlU1faG7i9YGGVQG2eUFZ8Z1cGyBoyLI\n7IgsAg6Oo6IiCyJrO7tHh13BVUQdOEfHGyo6qKDuILozR9GDOt4QdHAWFBgEbH4it6bv3XXJqsrM\nyoyMiP0jMqsys7O6K6syIyIzn8850JWRl/plVEY8+Ua88b6ZMAwRERGJUzbpAkREZPAofEREJHYK\nHxERiZ3CR0REYqfwERGR2Cl8REQkdsNJF3AwzrlnAK+q3jzGzD6SZD0iIrJ8vdDyeRFwKPAT4KUJ\n1yIiIh2QSeIiU+fc4cD1ZnZ+3bKXE7Vwfg+8EHiPmXnV+54DXAp8w8weir1gERHpqNhbPs65i4Er\ngRPrlq0EvghcZWY3AtuAv6ndb2ZPAx8Cbom3WhER6YbYz/mY2Zecc0cDr6lb/ArgqVpLB/glcC3w\ncefc+4BfmdmPnHPPcs6tM7OpeKuWBYwAf9Kh13oQ8A76KBHpC2npcPBcYLLu9iRwfPXn7wNHO+fe\nBrxPwZMqfwLc26HXehFwX4deS0RSLi3h8wxgtu52GVgLYGb3A/cnUZSIiHRHWnq75YBM3e1VwL6E\nahERkS5LS8vnd8C5dbefCWxd7ouGYRhmMpmDP1DSoFOH70Rk+bq+40wqfDI0vrmfA190zq0xszxw\nGvCtZf+STIbJyTxB0HtzFmWzGUZH16S6/qGhLOvXr25Y9r3vfY9cLsf4+DhbtmzhhS98Ieeccw5f\n//rXGR0dZXR0lHPPPXe/18rlCvh+EFfpB9UL638hvVw7qP6k1ervttjDxzn3JuBM4Ajn3MeAW83s\nLufcu4APOuceAQ4DrurE7wuCEN/vvQ9ATZrrz2T2r+uss87i9ttv52c/+xkbN27E8zxWr17NzMwM\nAEcccUTL10rr+0xrXYvRa7UPD0dnAWpHKzKZDEs5cFGppONLTK+t/7gl0dX6ZuDmFstvA26r3vxy\nrEVJxzz22GOcfvrp/PSnP+WKK65g27Zt/OQnP2F0dJTNmzfz8Y9/nDPPPDPpMiVlhoezHHZY47ft\n5lb1Yk1M5FMTQLKwtHQ4kD7x4IMPcsMNN3DUUUfx6U9/mh/+8IecdNJJTE1N8bWvfY1jjz026RJF\nJAXS0uFA+kSr8zkAV155ZcyViEiaqeUjIiKxU8tHRFLp1ltvZefOnWzfvp1Vq1axceNGzj//fL71\nrW+Rz+fZtWsX1157bdJlyhKp5SMiqXTKKadw2WWX8cADD/DqV7+ak08+mTvuuIOLL76YMAx1KLfH\nqeUjIqm0YcMGrrvuOj70oQ9x1113sWrVKg4//HByuRzDw8OsXbs26RJlGRQ+smSVSsDERL7lfdls\nhvXrV5PLFRZ1oZ26xkqzD3/4w4yMjHD33XfjeR4jIyOcdtpp3HvvvWzatCnp8mSZFD6yLAuFxtBQ\ndHWg7we60E6W5Oqrr265/NRTT425EukGnfMREZHYKXxERCR2Ch8REYmdwkdERGKnDgcikrj6npPt\n9pRs9VqSfgofEUmFWmiop+Rg0GE3ERGJncJHRERip/AREZHYKXxERCR2Ch8REYmdwkdERGKn8BER\nkdgpfEREJHYKHxERiZ3CR2QAhKFGCpB06evwGZueTboEkVTYMVFIugSRBn0dPvlZL+kSRFLhyX0z\nlDXgpqRIX4fPk3tnGM+Xki5DJHFbx/LM6MtYomY9P+kSUqWvw+c3T41pgxMBdkwWmNK2kKiSwqdB\nX4fPL7bsTLoEkVTYOpbn8T3TOvSWIM8PqPha/zV9HT76piEyb8/ULE/tm0m6jIFV9Hz2Tes0QE1f\nh4+IzCuWK9z16O6kyxhYfhAyU9KhzxrNZCoyIH63a4qhbCbpMgbSdNEjDEM8HXabo5aPyIB4ejxP\nrlhOuoyBND3r4QchutZ3nlo+IgPED7T3S0Kp4jNZKOs8dB21fEQGiL55J+e32ycZz6vlWaOWj4hI\nDHZMFFgxPJR0Gamhlo+ISAy2juWZ0jm3OQofEZEYTM16zJQqSZeRGgofEZGYBOrwMUfhIyISE82r\nNE/hIyISE3V1n6fwERGJSUkDu85R+IiIxERz+sxT+IiIxCTQOZ85Ch8RkZgUy+pqXaPwERGJic75\nzFP4iAwAdfFNh3ypomt9qjS2m3Tc8HCWbHXemKGhLJnM4je2ir4ZSp+b9XxWH6Jdr9aAdNTwcJbD\nDlszd3v9+tVtPX9iIq8Akr6mTgcRhY+IpNbw8NLPDKT1S0yhXOHQlSNJl5E4hY+IpNLQULbtlnO9\ntLaiyymsKQnqcCAiEiPPV/iAWj4SE8/z+NSnPsVxxx3H/fffz8qVKznvvPNwziVd2kAJ0fmGpO2d\nnuWYjWuTLiNxavlILPbt28cJJ5zA+vXrufzyy3nxi1/M3r17ky5LJHaTmkobUMtHYnLkkUdy+OGH\nUygUuPvuu9m1axdvetObki5LelCtFe2cY+vWreTzeS655BI2bNiQdGmLMj3rJV1CKqjlI7F66KGH\nuPHGG5menubhhx9OuhzpQfWt6PPOO4/x8XFGRnqn95jO+UTU8pHYbNq0iU2bNrF58+akS5EeVmtF\n79q1i7Vr1/LmN7+Zn//855x11llJl7Yo6u0WUfiISE/KZDJcd911rFq1igsuuCDpchatqGkVAIWP\niPSgTZs2AXD22WcnXEn7CiWNbA065yMiEquSWj6AwkdEJFZeoHM+oMNu0mGVSsDERJ5sNsP69avJ\n5QptDSGfxuFQRDrJq+hCX1D4SBdUKgFDQ9GUCr4f4Pva2ERqNMpEJPXh45w7CfhrYCMQmNmbEy5J\nRGLg+1EreqnUik63XjjnsxnIm9mbgPOcc5mkCxKReFQqwZL/k3RLpOXjnHs+8CMze3bdss3AGwAD\nTgRONrOimV1Qvf8QoGBmarOKtEkbjaRN7C0f59xXgJuBI+qWrQc+AZxhZhcB24FvNj31JuAtMZUp\nIiJdFHv4mNmFwDubFm8GpsysWL19G3Bq7U7n3OXA52MpUEREui4tHQ6OB4p1t3cCawGcc1cC/xso\nAAFweOzViYhIR6UlfDYC9ZNcFIAMgJldA1yz1BceGsoyMtIL/SoaZbNRv4rh4SzZbO8dsVf9yWlV\ne62XTgZSvz308rqH1vUPD8+v80wm3X+DWv3dlpbwmWR++wBYD3RkDIpD16xkdHRNJ14qEevWLX0O\n+zRQ/cmprz2bre7sMpme2R56ed1DY/1ri/PjuQ0NDfXM36Cb0hI+9wJn1N0+mqj1s2wz+VkmJ4c6\n8VKxymYzrFu3mqmp9kYISAvVn5xWtQfVIV3CMGRycunXzsShl9c9tK5/enr+rILv+6n+G9Tq77ak\nwqe5zXk9cI1zbqOZ7QXOAe7sxC/y/QDP670+/7URAiqV3hwhQPUnp1Xtc+8gJPXbQy+ve2hdf/11\nR2HK/wa1+rst9vBxzn0OOBMYds7dA9xoZl9wzl0DfM859zAwCpwed20iIhKP2MPHzJq7WdeWXwVc\nFXM5IiKSgPR2uRAR6UNh7x1J7AqFj4iIxE7hIyIisVP4iIhI7BQ+IiISO4WPiIjETuEjIiKxU/iI\niEjsFD4iIhI7hY/IINCFjZIyCh8REYmdwkdERGKn8BERkdgpfEREJHYKH5EBon4HkhYKHxGRGOkL\nQEThIyIisVv0TKbOuXVEM41+CTDgH4HXAfcBrzezJ7tRoIiI9J92Wj6fBf4CCIA3AucAFwBbq/eJ\niIgsSjvhcxZwgZkZ8Frg+2b2TeB9wKndKE5ERPpTO+EzDMw451YAZwA/qC5fA8x2ujAREelfiz7n\nA/wr8EkgB4wA33HOvYjokNuPu1CbiIj0qXZaPpcAPvDHwNvMbAw4F9gGbO5CbSIi0qcW3fIxsx1E\nYVO/7P0dr0hERPreAcPHOff2xb6Qmd2w/HJk0FT8gCAMWTE8lHQpIhKjg7V8rmi6/VygBDwOeMB/\nAFYC9wMKH2lbsewDKHxEBswBw8fMjqn97Jy7AjiJ6HzPeHXZOqKLTrd2s0jpX0WvouARGUDtdDi4\nHHhvLXgAzGwK+ABRZwSRtnmVgDDUaFcig6ad8CkCx7VY/iKi7tciixaGIUEY8tjeacIQAgVQV2nt\nStq0c53Pp4GvOOc+AzwAZICXAJcB7+lCbdLH8qUKITCZL+MHIV4l4JARHX6TAaAvWkB7Xa0/6pwb\nB94OvJvoy9SDwEVm9n+7VJ/0IT8I2TZRYCSbIVf0mPV8VgxnFT4iA6Sdlg9m9iWiDgYiS1au+JQ8\nn0IQMjYzy66pIjMlj/WrVyRdmojEpJ0pFYaIplD4I6D2FTUDHAK8yMw0uKgsShhCyfPZOVnE80O2\nbJ/keYevZdbzWanWT1fpiI+kRTstn88CbwXuBU4GfgkcCzwL+FznS5N+5YchhXKFXz+5j7UrR3hk\nZ46jNqxm52SBYzauTbo8EYlBO73dzgfebGYvAx4D3gEcDfwToOMlsmj5UoXH984wPlOiXAkolCps\nHy9EnRCqX839QF/RRfpZO+GzDrin+vODwIvNrAJcQzTXj8ii+EHIE3unyRU9CuUK07MeT+ydYWym\nRKkSALA7V0y4ShHppnbCZytwfPXnLcCJ1Z9LwGGdLEr6W7ni41VDZleuSKkSsGOywFP7ZihVouF2\nJgtlto7NJFmmiHRRO+FzA3CLc+6/ALcCFzvn3gdcD9zXjeKkP816PlNFD4CgengtX6pw16N7yBXK\nAEzkS+zOaY5CkX616PAxs48B/x2YMbN7iebweT1QBi7uTnnSj8ZnSuTLFaDxyvt8qULFj5bsmIy6\nX4tIf2r3Op+bAZxzw8BXzOymrlQlfW2iOqoB7D+sTq3DwWShxIqhdhrm84IwZLqo64ZE0qytrds5\n9w7n3KNAATjGOXeDc+4jzrlMd8qTflT0KnM/N193UrsZhPDkEs/5BEHI73ZPLbE6EYnDosPHOfc3\nwPuBjwK1vcftRF2ur+58adKv9k4tfC6nUD0cZztz7JwsLPl3FEuVgz9IRBLTTsvnMuBSM/si4AOY\n2TeAtwAXdr406TfFcoWxmRJ7p0sLPsarBFT8gJLns286ug6oXUGoUZxF0q6d8DkaeLjF8t8Dz+xM\nOdLPZj2fselZcsXygo/x/IBi2ccPQ2Y9n0rQfvgUy5W5LtsiaaMvRpF2wuffgfPqbtfW4TuIptEW\nOSA/CHlybGaum3Urnh/y+N5pZqvTay9lLLIHt02wfXzph+xEpPva6e32buBfnHOnEQ0merVz7njg\nBOBV3ShO+ssDW8d58OkJSt7CrZIwDNk+UaCyQG+4xSh5PmW//RaTiMSn3fB5DfBKYBZYC/wEOM/M\ntnehNukz2yYKjM2UDnjYIQjDxo4GS2j5lCqBxoYTSbl2wueVwHvN7KpuFSP97bHd0wftQFD2A57c\nN9/FOlhC+kzPeszM6gLVeopiSZt2wucTwE3OuU8CTwANIz+a2e86WZj0n0K5wuwBDrkBTBc9vLpD\nZlMFj/Wr2rtYNAzDuWF7RCSd2gmfD1f/fVndspBoQrmQ+QnmRFoqef5Bw2d3rkixPP+YsZkSz3nG\nmrZ+jx+E+qYvknLthM8xXatCBkKpEhywswFEo1nXBhcFWEqMlCvB3DA90kzrRdJh0eFjZk91sxDp\nf35w8DM4nh/O9XSD6ILRdpUq/pKeJyLxWdrIjSJt2l2dt+dg9kw3Dr3TbscBzw/YOVnksT3TTB/g\neiIRSZbCR2Kx2OtuvKaRCdodqcAPopERPD9Y0ugIIhIPhY/E4pEduUU9rlBuDJt7Htu36N+RK5R5\nYOv4fgOSakZUkfRR+EgsvCWOOFAb5Xqxv2Oq6M2d76n9u5TBSUWkuxQ+kmq7c9HlZIsZZuebdz/B\njon5Vs/TY3mK5Qq3/3Zn1+oTkaVR+Eiqlf2o2/Riek7vyhUbWkqzXjQt93h+4SkcRCQZCh+JxZKv\nu6nOzfPb7RMNi1tdrBqEsG9mvrecH4TYrhz7DjB/kIgkQ+EjqRZW//fo7umG5a0uVg3DsKF7dRDC\ntvE8UweYP0hEkqHwkVgs55rPohfNgFp/3qdVB4YwpOFaonL1YlNdcLo4Gglc4qTwkdR7ePskv3x0\nD15dsNQPwVMThiFjM/OH2KZnvSX3shtE7fQsFFmudsZ2S4Rz7ijgQmCzmR2ZcDmyVMv4Uu1VAvwg\naGj5tHq55mV+EEaDjGqcN5HUSX3Lx8y2m9nfEc2eKj1qOS2QiUKZIIRHd09Ft/MlHto2sd/jmofi\nKZR99k3P4vkB05rfpyVdAyVJSaTl45x7PvAjM3t23bLNwBsAA04ETjaz4gIvIT1mMeO6LaQ2xcK2\n8QIveM4GZmYrPLUvv9/jmkdH8PyAyWpwFcsV1q4cWXIN/arh+ik1ECVGsbd8nHNfAW4Gjqhbtp5o\nsrozzOwiYDvwzbhrk+5ZzEWiC4ueW6m2noIwxG/xes2H18IwXObvFZFuiT18zOxC4J1NizcDU3Ut\nnduAU5sek+lyadJFD7c4TLZYtU5Yv3oiGuctZP8OB//ym20NUzEAFD2fXCE63KYMEkmXtJzzOZ7G\nabl3AmsBnHPPcM79A7DGOfeBJIqT5ZspLb0nVS04pqrX8ARBONcKqtmV2/8IbcWfn7xO2dNa2NCJ\nQ2tJ4pOW3m4bgfqvsgWqLR0zGwMurf7XtqGhLCMjacnYxctmo4be8HCWbLb3dgrN9S+nx1mm7s83\nMpJlaCgDGRr+rq0uUQlC5qZVGB7OtPU56OX136r2+sMG9euhEs5vHyOVdGwrvbzuoXX9w8ON6zUN\n63khtfq7LS3hM0nj9rEeaG8ilwUcumYlo6NrOvFSiVi3bnXSJSxLrf7p2aW3fB7fG02JkMlmGB1d\nw+7H9lEJaPi7TrSYOK5Y8eeG4VmxcsWSPge9vP7ra5/boWQyDethRanC6kOqu4GRYUbXpKdTaS+v\ne2isf21x/vM/NJTt6X1Sp6QlfO4Fzqi7fTRR62fZZvKzTE4OdeKlYpXNZli3bjVTUwWCHrzyvL5+\n3w+W1dV6T3V+nsAPmJzMs2tshnLFZ3JyvsdbrsXgoYXZ+ekVJnIFJlcu/uPey+u/Ve21f8MwbFhv\nhXKF8opoveQKZfCSv9C0l9c9tK5/enr+sLBf/RynVa3+bksqfJrbnNcD1zjnNprZXuAc4M5O/CLf\nD/C83ruWYWgo+qZaqQT4fu9tgPX1e8u8lsSf23GC50XrIwzChr+r3yLcKnXrreK19zno5fXfqvb6\nd1C/HjwvwMtEt71KOraVXl730Lr+StM2kIb1vJBa/d2WRFfrzwFfA4adc/c4595hZnngGuB7zrmb\ngFHgtXHXJt1RXOawLbPVDdWvdjTYPVWkVGm8cLTVdUT5uk4Oy7nOSKST1PMyEnvLx8yau1nXll8F\nXBVzORKD5U5pUBtzrOwHFMoV7ntynEK5wrbxPCf8wWj0mBa96SbrumPnl9HbTkQ6L71dLqRvdPKL\nXlCbY4HGb5CtLjptqEFfN1sLF/hZpMsUPtJ1ndzxB0Hrq1EONh2A9qsi6aLwkZ6yK1ec60H0++pA\nowB+cOBzOs0DaNrOHI/szHW+QBFZFIWP9IwwDNk9VZwbr60+PJqH1mlWaQqnR3bk2LJjsvNF9jC1\nDiVOCh/puo6ebqmbmbR+iJ2D/o6m+wNC9TpCgSPJUfhI13U4e+ZaPn64+GF79ntUqE4IIklS+EjX\neZWOjJQ0p5YZxbJPsezPDR7ajny5MneoruIH5EuabK4XRxOQ3qXwka7bPtmRkZIIaRx5eetYnkd3\nT7FjEa/f3EfuN1vHmagOyTNV9Pj1E2MdqbGXaVZTiZPCR7qv08fd6m8u9rWbn7fwXQNmsN+9JEfh\nIz2jVdCES+w4EIaN0zBoxlOReCl8pOu62PBpOY/PomuYC5xQDQCRmCl8pGcEYbjfGG6ThcWNG9d8\nPiMIw4ZrgyaL5eanDAT1MZCkKHyk6zp1RMsPQrbsaByV4OeP7F5Uo2XPVOM0234QUiz7c/X9wnZ3\npsge03CtlJp/EiOFj/SU5h3kwcZ0m3vefh0VQu1sRRKk8JGe0hwineoooP4GIvFS+EjXdbKF0fxK\nlSBcdnKEddM0iEg8FD7SW5qCZqlX5TfnzaC2fAb1fUvyFD7SdV4Hr5yvTaldUyhX2JUrLvDoxdk2\nkadUCXSFv0iMFD7SdZ2cwjrX1CV6quhxz+P7Dvq8lgOLVn+885FdTM96Azm+m0Z6kKQofKSnaCQC\nkf6g8JGu6+gIBy2H2FnM8xofFdY9b7DzbKDfvCRI4SM9pVUHg07NyzPYISQSL4WP9JRWh92WNLBo\niycOYvYocCUpCh/pKX6r8OnYq2tPLBIXhY/0lFaH3Tz/4F2k61tMs55PGESXvta/XmnQu1oreyVG\nCh/puk4e2mkVEL/dPnnQ500V57tR37FlF1Oz0e1csTy3z/23R/d0pEYROTiFjwygKG7C5Y/M0/MG\n/f1LchQ+MtC07xVJhsJHYpD8Lr7+G37jt/2wY121e5GmlZCkKHxk4MxPnh3qsFMdrQqJk8JHBprC\nRyQZCh/purTu4MNQh53S+reR/qfwkYHQMHpz3R634VxQfOWkxiC+Z0kHhY8MNO18RZKh8JGuK5Q7\nN59PJ9QC56l9M2zZMckjO3OJ1pMopa8kROEjXecvcarrbitVAvKlCsWyn3Qp6aATQBIjhY90XSp2\naQvsWBe+/kdEuknhIwNtkC8wFUmSwke6L8X7d2WPSDIUPjIQFsqYoLETdhyliAgKH4lB6i7k1Hke\nkcQpfGSg6ZzPPK0JiZPCRwZaSnuBi/Q9hY8MhIUyZqGhdkSkuxQ+IiISO4WPdF2aWxQpLk2kryl8\nZOAsNMK1iMRH4SNdl4bdu+8HLZff8/i+uZ+3TRTiKkdk4Cl8ZCAs1Kttz9Ts3M8zs15M1YiIwkdE\nRGKn8BERkdgpfKT7UnZSP3XD/YgMIIWPiIjETuEjIiKxU/hI16XhIJcOtYmki8JHRERip/AREZHY\nKXxk8OgInEjiFD7SdWnY16est7fIwFP4iIhI7BQ+IiISO4WPiIjETuEj3afzLSLSZDjpAg7GOXck\n8H+AfwX+wMzennBJ0uOUhSLJ64WWz3XA7Wb2QeA1zrmRhOsREZFlSqTl45x7PvAjM3t23bLNwBsA\nA04ETjazIvAfgVuqDysALwZ+GW/FshxpaGmoq7VIusTe8nHOfQW4GTiibtl64BPAGWZ2EbAd+EZd\njbUpJoeASmzFiohIV8QePmZ2IfDOpsWbgalqSwfgNuAV1Z9/C7jqz2uA+7tcooiIdFlaOhwcDxTr\nbu8E1lZ/vhz4snPuUOCHZuY1P1lERHpLWsJnI1Cuu10AMgBm9jTwn5fyon927OEcumqEkZFe6FfR\nKJvNADA8nCWb7b0TFvX1H3vEWjw/SLSeow5bPfc5OGrDak78w2fs/5gN84/p5fXfqvYXPOcw8qUK\nhwxnG7aHNauG526nZVvp5XUPretft2bF3Gdu47qVqVjPC6nV322ZMIEzsc65U4A7zWykevsbwCYz\ne1719uuAW8wsLeEoIiIdlJb4vRdYV3f7aKLWj4iI9KGkwqf5914PrHfObazePge4M96SREQkLrEf\ndnPOfQ44EzgG+BVwo5l9wTl3dXX5w8CfASeqc4GISH9K5JyPiIgMtrSc8xERkQGi8BERkdgpfERE\nJHYKHxERiV1fXsTpnPtT4AvAHwG/Ay4zs7sTrumlwN8TDSW0F/iYmd3gnBsFbgROByaBD5nZjXXP\nuxb4r0SDqn4NeLeZhdX73gB8BDgc+BlwsZnt6eJ7OAL4DfA2M/t+j9V+FNFn4uVAjmj9f6YX3kP1\nouxPA8cBO6o13tIjtW8Cvm1mR1Vvd6XmbmzzLWo/Cvgs8DKiEVn+GXhPrVdummpvVX/d8gxwO/Ar\nM/ufdctjrb/vWj7OuUOA7wJfBtYDnwG+65xbnWBNo8B3gE+a2SjwV8A1zrlXAl8EpomGGHod8NHq\nh6Y2zcSZwB8DJwAvBd5Tve8FwOeB1wPPBHYDN3X5rXwZ2FB3+0s9VPutRN34DwP+ArjKOXdS2t+D\ncy4LfBu4xszWA5cAX3XOPbcHar8I+CFQPwdXx2vuxja/QO3/CDwNHAm8kGh6l/enrfYD1F/zt9X6\n6h8fe/19Fz7AaYBvZjeYmW9mNxGtrFcnWNPRwG1m9k0AM7sfuAM4heiC2g+YmWdmvwK+Dryl+rwL\ngE+Z2Z7qt4xrgQur970RuNXMfm1mJeB/Aa+qu1C3o5xzlxLtNJ6u3l7TQ7W/hGiHcYWZBWa2BTiZ\nqBWR9vcwSrTB13YiIVACgjTX7py7EngX0bfl2rJufWZOp4Pb/AK1jwAzwEeqte8hmhrmlDTVvlD9\ndfe9oFrXt5vuir3+fgyf44mmYahn1eWJMLMHzOyttdvOucOImu4Anpk9Vf9w5mttfi/G/PQSDfeZ\n2TgwXnd/xzjnjgPeDVxGdcBXokn+ymmvverE6u/7mHNup3PuEaLw2UDK30P1dT8PfMM55xGN/LGZ\nKJDSXPuXzexPgV/XLTuuwzWPVe93dHab36/2auCc3XRo8mzg31NWe8v6AZxzK4CvAhcD+abnxF5/\nP4bPGvYfF64AJHbYrV514rzvEo3ucAeNU0lAY63N76UAZKsfoljep3Oudvz3XWY2WXfXGlJee50N\nRC3iPcBzgLcRTc9+KCl/D9Xj8wXgL4FVwGuIzv+sS3PtZra7xeLVdLbmYvW5HX0/C9TewDl3HdGO\n99rqolTUDges/1rgB2b2by3ui73+fgyfAtFGWm81UZM5Uc65Y4C7gH1EO5MZYGXTw+prbX4vq4GK\nmZVb3Nf83E75AHC/mf2oaXmB9NdeUwLGzOyjZlapbnzfAq4m/e/htUQjvn+7Wvv3iSZb/GAP1N6s\nW5+Z2N6Pc26lc+6fgT8HXm5mY71Qu3PudKJDZB9Y4CGx19+P4bOF/Q8ftGoaxso5dyLw/4i+eZxX\nPXb6KLAjrujJAAACAklEQVTCOffs+ocyX2vzezm+umy/+5xzzyQ6mb6Fzvor4K+dc+POuXHguURT\nnJ/VA7XXGDBcbUXUDAH3kf738FzgkKZlFaKR4NNee7Nufd5j2earh8vvJDqxfpKZba27O9W1E3UW\neB6wp7odvxHY7Jz7blL192NX69uBQ5xz7wT+gehk5uFEPT8SUe2i/APg783sY7XlZjbjnPsOcK1z\n7u1EPU3eQNTrBKLeNZc7535GtMN5L9EhMIBbgDucczcS7USvBb5vZhOdrN3MTmh6L08A/83MfuCc\n+09prr3Oj4m+oV3lnPsw8BLgXKJvr3+Y8vfwY6KekW81s686506t1n56D9TeoFufd+dcXNv8t4lm\nWf5LM/Ob7kt17WZ2KXBp7bZz7iZgb11X69jr77uWT7WZeCZRso8B7wTONrPmY81xuojoBPH7nXPT\n1f+mqjvCS4AVwDbgn4C/NbPaicLribpo3wM8BPwC+CREnRiqz70J2AU8q/p7ui1kvtPB23uhdjOb\nBV5BFDp7iDa0d5nZPWl/D2b2EHA+8D+cc5NEXVnfYmb3pb32BXS85ji2eefcyUSdhP4cmKxuv1PO\nuTvSXvsixV6/RrUWEZHY9V3LR0RE0k/hIyIisVP4iIhI7BQ+IiISO4WPiIjETuEjIiKxU/iIiEjs\nFD4iIhI7hY+IiMTu/wP9f/VchL4aqAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2b8e275eaa10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wasabiplot(bam_filename2, chrom, start, stop, strand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation = read.cigar[1]\n",
    "cigar_operation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.end_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.ref_iv.start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cigar_operation.check()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "read.iv.start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# subset_reads = bam.fetch(reference=exon1.chrom, start=snap25_regions[0].start, end=snap25_regions[-1].stop)\n",
    "# subset_reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Junction plotting code:\n",
    "\n",
    "```python\n",
    "for jxn in jxns:\n",
    "    leftss, rightss = map(int, jxn.split(\":\"))\n",
    "\n",
    "    ss1, ss2 = [graphcoords[leftss - tx_start - 1],\\\n",
    "        graphcoords[rightss - tx_start]]\n",
    "\n",
    "    mid = (ss1 + ss2) / 2\n",
    "    h = -3 * ymin / 4\n",
    "\n",
    "    numisoforms = 0\n",
    "    for i in range(len(mRNAs)):\n",
    "        if leftss in sslists[i] and \\\n",
    "            rightss in sslists[i]:\n",
    "            numisoforms += 1\n",
    "    if numisoforms > 0:\n",
    "        if numisoforms % 2 == 0: # put on bottom \n",
    "            pts = [(ss1, 0), (ss1, -h), (ss2, -h), (ss2, 0)]\n",
    "            midpt = cubic_bezier(pts, .5)\n",
    "        else:                         # put on top \n",
    "            leftdens = wiggle[leftss - tx_start - 1]\n",
    "            rightdens = wiggle[rightss - tx_start]\n",
    "\n",
    "            pts = [(ss1, leftdens),\n",
    "                   (ss1, leftdens + h),\n",
    "                   (ss2, rightdens + h),\n",
    "                   (ss2, rightdens)]\n",
    "            midpt = cubic_bezier(pts, .5)\n",
    "\n",
    "        if number_junctions:\n",
    "            text(midpt[0], midpt[1], '%s'%(jxns[jxn]),\n",
    "                 fontsize=6, ha='center', va='center', backgroundcolor='w')\n",
    "\n",
    "        a = Path(pts, [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4])\n",
    "        p = PathPatch(a, ec=color, lw=log(jxns[jxn] + 1) /\\\n",
    "            log(junction_log_base), fc='none')\n",
    "        axvar.add_patch(p) \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Code to get read counts and junctions\n",
    "\n",
    "```python\n",
    "def readsToWiggle_pysam(reads, tx_start, tx_end):\n",
    "    \"\"\"\n",
    "    Convert reads to wiggles; uses pysam.\n",
    "    \"\"\"\n",
    "    wiggle = zeros((tx_end - tx_start + 1), dtype='f')\n",
    "    jxns = {}\n",
    "    for read in reads:\n",
    "        # Skip reads with no CIGAR string\n",
    "        if read.cigar is None:\n",
    "            print \"Skipping read with no CIGAR string: %s\" %(read.cigar)\n",
    "            continue\n",
    "        cigar_str = sam_utils.sam_cigar_to_str(read.cigar)\n",
    "\n",
    "        if (\"N\" in cigar_str) and (cigar_str.count(\"N\") > 1):\n",
    "            print \"Skipping read with multiple junctions crossed: %s\" \\\n",
    "                  %(cigar_str)\n",
    "            continue\n",
    "\n",
    "        # Check if the read contains an insertion (I)\n",
    "        # or deletion (D) -- if so, skip it\n",
    "        for cigar_part in read.cigar:\n",
    "            if cigar_part[0] == 1 or \\\n",
    "               cigar_part[1] == 2:\n",
    "                print \"Skipping read with CIGAR %s\" \\\n",
    "                      %(cigar_str)\n",
    "        aligned_positions = read.positions\n",
    "        for i, pos in enumerate(aligned_positions):\n",
    "            if pos < tx_start or pos > tx_end:\n",
    "#                print \"=>\",pos\n",
    "                continue\n",
    "            wig_index = pos-tx_start\n",
    "            wiggle[wig_index] += 1./read.qlen\n",
    "            try:\n",
    "                # if there is a junction coming up                \n",
    "                if aligned_positions[i+1] > pos + 1: \n",
    "                    leftss = pos+1\n",
    "                    rightss= aligned_positions[i+1]+1\n",
    "                    if leftss > tx_start and leftss < tx_end \\\n",
    "                           and rightss > tx_start and rightss < tx_end:                      \n",
    "                        jxn = \":\".join(map(str, [leftss, rightss]))\n",
    "                        try:\n",
    "                            jxns[jxn] += 1 \n",
    "                        except:\n",
    "                            jxns[jxn] = 1\n",
    "            except:\n",
    "                pass\n",
    "    return wiggle, jxns\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}