20 Biopython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biopython\n",
    "\n",
    "Večina primerov je povzetih iz [tutoriala o Biopythonu](http://biopython.org/DIST/docs/tutorial/Tutorial.html).\n",
    "\n",
    "Biopython nudi vse glavne funkcije, ki jih potrebuje bioinformatik:\n",
    "\n",
    "* [Branje in pisanje datotek (FASTA, FASTQ, GENBANK)](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc47)\n",
    "* [Poravnava zaporedij](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc69)\n",
    "* [Uporaba programa BLAST](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc87)\n",
    "* [Dostop do zbirke NCBI](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc111)\n",
    "* [Dostop do zbirke Swiss-Prot in ExPASy](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc137)\n",
    "* [Delo s 3D strukturami PDB](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc149)\n",
    "* [Filogenetika](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc198)\n",
    "* [Analiza kratkih zaporedij - motivov](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc211)\n",
    "* [Gručenje](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc231)\n",
    "* [Nekaj nadzorovanega strojnega učenja](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc241)\n",
    "* [Izris genomskih diagramov](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc254)\n",
    "* [Dostop do zbirke KEGG](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc272)\n",
    "\n",
    "V vajah se bomo osredotočili na:\n",
    "\n",
    "* [Delo z datotekami](#Delo-z-datotekami)\n",
    "* [Dostop do NCBI](21 NCBI.ipynb)\n",
    "* [Izris genomskih diagramov](22 Diagrami.ipynb)\n",
    "* [Filogenetsko analizo](23 Filogenetske.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Naprej uvozimo modul Biopython."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Bio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Delo z datotekami\n",
    "\n",
    "Razred `SeqIO` omogoča delo z datotekami."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Funkcija `SeqIO.parse` poskrbi za pravilno branje zapisov v datotekah. Seveda, če ji namignemo za kateri format gre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zap1\n",
      "Seq('TGGCGTCTTTTGTTGTCGAACAAATTGATCGTGAGAGATA', SingleLetterAlphabet())\n",
      "40\n",
      "zap2\n",
      "Seq('AATGCCAACACCAAGTAAGTAAAGACAGCCTCACTAAGCTGGAAATTTCGTGGA...ATT', SingleLetterAlphabet())\n",
      "500\n"
     ]
    }
   ],
   "source": [
    "for rec in SeqIO.parse('datoteke/primer.fasta', 'fasta'):\n",
    "    print(rec.id)\n",
    "    print(repr(rec.seq))\n",
    "    print(len(rec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zgornji primer bo vrnil izjemo, ker datoteka `primer.fasta` ne obstaja.\n",
    "\n",
    "Lahko pa jo naredimo z uporabo funkcije `SeqIO.write`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Alphabet import generic_nucleotide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "zap1 = ''.join([random.choice(['A', 'T', 'C', 'G']) for _ in range(40)])\n",
    "rec1 = SeqRecord(Seq(zap1, generic_nucleotide),\n",
    "                 id='zap1',\n",
    "                 description='nakljucno',\n",
    "                )\n",
    "\n",
    "zap2 = ''.join([random.choice(['A', 'T', 'C', 'G']) for _ in range(500)])\n",
    "rec2 = SeqRecord(Seq(zap2, generic_nucleotide),\n",
    "                 id='zap2',\n",
    "                 description='tudi nakljucno',\n",
    "                )\n",
    "\n",
    "SeqIO.write([rec1, rec2], 'primer.fasta', 'fasta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zdaj lahko poskusimo prebrati datoteko `primer.fasta`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zap1\n",
      "zap1 nakljucno\n",
      "Seq('CATCGTACTGGAGATCTGTCAGGCGCTAATTAACTCTGTA', SingleLetterAlphabet())\n",
      "40\n",
      "zap2\n",
      "zap2 tudi nakljucno\n",
      "Seq('CTGACAGCTACCGCGTTGACCGGAACCCGAGTGATTCCAGTTGGCTACAGCTGA...GCT', SingleLetterAlphabet())\n",
      "500\n"
     ]
    }
   ],
   "source": [
    "for rec in SeqIO.parse('primer.fasta', 'fasta'):\n",
    "    print(rec.id)\n",
    "    print(rec.description)\n",
    "    print(repr(rec.seq))\n",
    "    print(len(rec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kot smo videli že pri ustvarjanju objektov `SeqRecord`, vsak zapis vsebuje ID, kratek opis (ki vključuje tudi ID) in zaporedje.\n",
    "\n",
    "Objekti tipa `SeqRecord` se obnašajo podobno kot nizi. Omogočajo indeksiranje, rezanje, itd:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('CTGACAGCTACCGCGTTGACCGGAACCCGAGTGATTCCAGTTGGCTACAGCTGA...GCT', SingleLetterAlphabet()), id='zap2', name='zap2', description='zap2 tudi nakljucno', dbxrefs=[])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'G'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('GACAG', SingleLetterAlphabet()), id='zap2', name='zap2', description='zap2 tudi nakljucno', dbxrefs=[])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec[2:7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('GACAG', SingleLetterAlphabet())"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec[2:7].seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'GACAG'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str(rec[2:7].seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('GACAGTA', SingleLetterAlphabet()), id='zap2', name='zap2', description='zap2 tudi nakljucno', dbxrefs=[])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rec[2:7] + rec[8:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"green\"><b>Vaja</b></font>\n",
    "\n",
    "Napiši zanko, ki bo prebrala datoteko in izpisala frekvence posameznih nukleotidov:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C: 131 (24.26%)\n",
      "A: 146 (27.04%)\n",
      "T: 122 (22.59%)\n",
      "G: 141 (26.11%)\n"
     ]
    }
   ],
   "source": [
    "cn = {}\n",
    "\n",
    "for rec in SeqIO.parse('primer.fasta', 'fasta'):\n",
    "    for n in rec.seq:\n",
    "        cn[n] = cn.get(n, 0) + 1\n",
    "\n",
    "for n, c in cn.items():\n",
    "    print('{:}: {:} ({:.2f}%)'.format(n, c, c*100/sum(cn.values())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"green\"><b>Vaja</b></font>\n",
    "\n",
    "Napiši zanko, ki bo za dano zaporedje `s` izrisala frekvenco nukleotidov v oknih dolžine `k`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def plotAT(s, k):\n",
    "    cns = []\n",
    "    for i in range(len(s)-k):\n",
    "        cns.append(s[i:i+k].count('A')+s[i:i+k].count('T'))\n",
    "    plt.plot(cns)\n",
    "    return cns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD8CAYAAAB9y7/cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXm8HGWZ739vL6c7Jzkn+0oSkpAQIAECJICyiCzCiFdU9CrjzOAMDo4L46BzHUXvdXCZ0TtXcZxxFly5oLiiInJldcQlLAkQEhIIAbJvJ+tZu0939Xv/qHqr3qp+q7q6lu6q7uf7+UBOd1dXv1X11lPP+6yMcw6CIAii88i0ewAEQRBEPJCAJwiC6FBIwBMEQXQoJOAJgiA6FBLwBEEQHQoJeIIgiA6FBDxBEESHQgKeIAiiQyEBTxAE0aHkWvljM2bM4IsWLWrlTxIEQaSe9evXH+Kcz2z2ey0V8IsWLcK6deta+ZMEQRCphzG2I8j3yERDEATRoZCAJwiC6FBIwBMEQXQoJOAJgiA6FBLwBEEQHQoJeIIgiA6FBDxBEESHQgKeMHl48wHsP15q9zAIgogIEvAEAIBzjvfdtR7feyJQPgVBEAmEBDwBAChXa9BqHKPjWruHQhBERJCAJwDoAl7+lyCI9EMCngAAlCu65l6qkAZPEJ0CCXgCAGnwBNGJkIAnAADlqmb7lyCI9EMCngAAlCqkwRNEp0ECngAgmWgqJOAJolMgAU8AkJysZKIhiI6BBDwBgDR4guhESMATAMjJShCdCAl4AgCFSRJEJ0ICngBgJTiVyERDEB0DCXgCgKzBk4mGIDoFEvAEAMu5SiYagugcSMATACzNfbxaA+e8zaMhCCIKGgp4xti3GGMHGWObFJ/9LWOMM8ZmxDM8olXImjtp8QTRGfjR4L8D4Crnm4yxBQCuALAz4jERbUCuIkmx8ATRGTQU8JzzxwAcUXx0G4CPAaD1fAdg1+DJ0UoQnUAuyJcYY28GsIdzvoExFvGQiHYga+1kokk3j2w5gG/+7lUAQC6bwS1vPAWnzOlv86ji4ZmdR/HoCwfx0Tcsb/dQEknTTlbGWC+ATwL4Xz63v5Exto4xtm5gYKDZnyNahKy1kwafbn753D6s23EUpYqGx7YO4HcvHWr3kGLjrf/2B/zLo9soMMCFIFE0JwFYDGADY2w7gPkAnmaMzVFtzDm/nXO+mnO+eubMmcFHSsSKnOBEyU7pplytYcHUCbj7xvPN153OuNb5xxiEpk00nPONAGaJ14aQX80571w1oQsgDb5zKFc1FPNZ9GQzYKw7BHy5WkMhl233MBKHnzDJuwGsBbCcMbabMXZD/MMiWo3NyUoafKrRhV0GjDEUchmzFHQnQ3NWTUMNnnN+XYPPF0U2GqJtlKs19BVyGCpXu0Lj62RKFc3UZgu5bFdcT2oWr4YyWQkA+rK+f0Le/JtIL+VqDYW8fmsXcpmuuJ7d8BALAgl4AoDuWBUCnpys6aZcqaFoaPDFfLYrzBfd8BALAgl4AoB+g0yekDP/JtJLuarZNPhuaMNIGrwaEvAEAF3r6y8KEw3dLGmmVNGdrABQyGe6QoMnG7waEvAEAF2omzb4LhAInUy52n1O1m44xiCQgCcA6BqQ0OBJG0o35WoNxW5zspJSooQEPAHOOcrVGiYVssh0SWJMJyMn/RTz2a5wmnfDQywIJOAJM827kM/qURd0s6SWqlaDVuOWDb5bNHhSSpSQgCfMm6OQyxgCgW6WtFIS19Jmoun869kN2bpBIAFPmPbLQj6rO+W6YEnfqQhBV8xLTtYOvZ5yBclueIgFgQQ8YTpVC7kMCvnuiJvuVOTVGAAUO/h6VjQS8I0gAU/Um2g6VOPrBqxraWjwHZzJaquASiYaJSTgCfNGKZKTNfXIqzHxb7mqdWRDDFlrL5EGr4QEPEFO1g5CXDvLBp9BjQPVWmcLeNLg1QTqyUqki58+sxs3/2ADPnX1qXjvRUvM96+7/XHMnVzEUzv0nurFfBY9uQx+v+0wfvvSAC5a1nwHrvfftR6/2TqAbIaBAXjw5tfhP37zMr73xE4smtGLX334YmQyye7j+6fffAJPvHIErz9lJv7zT1fXff7pn2/C+p1Hcd9NF7VhdN6U6zR4XdCfeeuDqGocEwtZ3PuhCzFtYg9WfPoBfO2Pz8bVZ8xt23ib4Vu/exWfuW8zvvSOMzGjr4C/unO9+VlcSslT24/g3d94AteePR//+LbTY/mNOCENvgvYvHcQALBpz3Hb+2tfOYx7ntmD46MV5DIMZy2cgvdeqD8AXtw/FOi3/t+m/Rgd1zBUqmKwVMWuo6PYsPsYxrUath4YxlgKNK0Nu/TxPrf7uPLzO9buwKY9gy0elT/KjjBJkdE6Oq7htHn9ODpawc4jo9h9dAwAcNvDW9sz0ABs2qtfj837BvHi/kGMVTR86PVLMaU3H5uA33pgCOPVGn79wsFY9h83JOC7ADH53W6Cao3jPa9dhEIui9ecNN1z26Z/u1KzOfnSYP5pdL6STJ2TVWpjd/mps4xtkv+QVWFdF82cUx++fBlmTirEdkzid1iyF52ukIDvAsQkdRNYeu0Sy2brtW3Tv13VUtXvVZRtANJp1y2ZcfBWNUnBZKmYHIduk0+T3DLncaWGcrWGbIYhnzVCe2OKFBJzIU3nSYYEfBcghKqbcJVT280+nhEJ4nJVvxmF2T3pIXuibEOGpTMyo16Dt25xs6FLVcN4Ko9Nn5Olas1oS2j5GWLT4I39ptVFTQK+CyhJmo8bsqYXNBZeFYqna/A1m3BJMkJA9k/IQ6txVDX386AlMDJFCCSrHrxlojHr/VdqqSxAZmnwmtlYHNBXK3EpDmI+VDzmQZIhAd8FWJqPJVydwli21RYCxsKrQvFKlZqtFHHSNXhh4vDT/CSJWrBZdkKpwYuOXbXEm8pUWCvRWstq3ov5kMYHIkACviuwbMrWJHUK42IEGrzqJhPalmn/TaBQlBHH7We8SRSS4iFuFRuzHtyTpabqSX/QqrA5WR017+PqYSD/ZhohAd8FqKJCnIJLFgR6NmvzAkB1k5WqNYzbBHyybxRx3H7Gm8SHlaXBZ2z/AnYTTRLH3gh5HpcrNdsqJa7jEeezovFEmuQaQQK+CxCCVxZWTmEsC4KgGpHqJhsqVQBY5oGkL3XFOfIz3iR2vipXa+jJZcCMuD55ZdZbyCGXYShVtUSOvRGyuaRkaywev5MVSKZJrhEk4LsAXxq800QTYDKrwgoHx6oAZJt2sgWL6WT1Md4kasG6bVq+lpJvRSoml8SxN8IZB29zssalwdvumWTPXRUk4LsA08kqCWCnMC7aBEEwjUh1kx0fExp8upysk13Ga6tBnsBjkXMaAPuDO5dhhgM9pU5WsRKt6OM3czdirJgp3zNJX32qIAHfBciJTkJAeWnwQTUiLwGfGierFCYpvxbYa5AnT0jK8eGAXYNnjKFo5Dgk/TqosNngpTDJQk6veR9HxUzS4InEI7QQzi0BVW+Dd2jwAbQVlV13UNjgiznXbZKEOG638cqhpkkUkrLgA+y+FUDXdktS+Yi0uA2rWs2M/CpXNONBZjlZ5bkdJV5mzTRAAr4LEGnd+t9WLLGMTSgE7ALk3Gc2wzDoNNEk/CaxnKzq8coPviQ+rOToEkAh4A0NXlxfr0SuJCGuQzbDFBp81tgm+utRrmjIiXuHTDRE0hC1VYRG6lZIy2a3DRoH7xB4kyfkMVgynKwdEiZZTrwGr9kiZ5ijSpZwoDeqT5Q05OsyrhmlCow5K443jmNJU4ivChLwHY6orSImqRky6RkmGY2Ttb+YM23wE/JZ9GST30xEnBc3p3DSbbJODd6JML+pHO9JxlxZGYrKYKlap8HHcSylimaV2ehEDZ4x9i3G2EHG2CbpvX9ijL3AGHuOMfZTxtiUeIdJBMXNaegspGVPdAomiJ03WP+EvBk7nJZ+r84wSaepyh6JlLxjKUvx4SqE+S1tJZGFcBWKyni1ZsXBx6jB21e/6XgYyvjR4L8D4CrHew8BWMk5PwPAVgCfiHhcRETU1VaRCjbJ2OPgg2Wy1mvweds+g9r2W0mdicZTg0+ecHQ6WZ1YGny6BLzTNwLU19uJ44ErF8pLy7mSaSjgOeePATjieO9BznnVePk4gPkxjI2QGClXMTpebbyhA2dtlRf2D+L4WAWvHhqxbeeMvNBqHNsODuPgYAnbDg7hud3HsPPwqOdvvTwwbHst2/WL+YwtOqdU0cws16gYHa/iud3H8NzuY3hx/5AZNndsdFxZDXBgqAxAP7fie7uOjCLDgEkFXWvbcXgEW/YNolbj2H5I/1tQqmg4NFxGrc0p7HuOjWHPsTFUtBpe2D9kO+9OCnnDyWo84LUaNyOdkowZ3SQJeKvmvX68L+wfjNRpXKpoGC5XzXtn++GRBt9IHlH0ZP0LAD+IYD+EBys+/QAKuQxe/NwfNfU9oXXM7i8CAD7yww3K7WRnnNDmL//ybxzbAE/cchlm9RXrvn98tIJv/347AGD57D68eGAIU3qtm7GvmLfVmb/8y7/B7qNj2P6Fq5s6Hi9uuWcjfvbsXvP1t9+zBq8/ZRZWfeYhvPH0Ofi3d59jfvboCwfwF99ZhztvOBfff3IXfrlxn/nZ1F59rMV8Bnes3YE71u7AR684GV96yN7ebseRUaz+3MP46BUn46bLlkV2HM0wNq7hgi88CgC46dKlAKyHuQrTySppo5/86Sb8y3VnxTvQkJjzWJp74jinGP9+5IcbcGCwjPdfclIkv/nZ+zYDABZM6wUA/O9fvYi/vGgJ8tn0uC5DCXjG2CcBVAF812ObGwHcCAALFy4M83NdTzCziS5Qz108DW9YMRvv+faTKFVqWDZrEm575yqUqxr6inaBoNIA33j6HNy/cT+OjVaUAv7Y2DgA4IYLF+PmK07G0ZFxTCzk8Ecr52D6pAJm9hXMLEoAZk/QKDk0PI6TZk7EjRcvwd/9ZCMODZfNAlH3b9xv23bd9qMAgGd3HsPAcBnLZ/fhY1ctBwAsnNaLTIbhZx+8AC8dGMZNdz+DFw/oPWpvvvxknL9kGv74G0+YK5qHtxxom4CXV0FiVfbRK5bbtnnylsuQy4q0fhEHr2HVgil4dtcx0xGeZMSK48oVs3HZqbNQrXGct3gaAGDVgin4/o3n4z3ffhKHh8uR/eYhY19/feky7Dk6hns37MXouIbJE7pAwDPGrgfwJgCXcY8UMs757QBuB4DVq1enJa+iYxDOqWI+g/OXTEdvTw6lyjhm9xex8oTJyu+obLgr5k3G/Rv3u9o5heA+a+EUTCrkTBPHZafOtu03zqiNclXDrL4iXr9c9B71l5JfrtYwe3LRNlYAOGVOPxZNnwjAysg9d/E0nLdkupkR2m7kh/7xsQqmTezB5F77A3tWv/VALkiZrP0T8jh38bRUtCYUxzmxkKubt4wxc25HaScvV2s4Y/5kTOjJYs3iabh3w17jmruvkJJGoEcRY+wqAH8H4M2cc2/DLNFWrA4/dodU0SvSQhFm1ygWWAj+omeIXrxhkqJGuLDJyvHebjCmO5zdHJPifRHPL9t9kxA2J1+PwVIVRQ8HK2A50EsVDUUR2ZQC56GzU5WKKFtNAvqcLrbAkRsnfsIk7wawFsByxthuxtgNAP4VQB+AhxhjzzLG/iPmcRIBETevEExWc21vQeykUS0Z8wb0DNGLr/MOYMWAW43D3WuuyPk/444CXfbtGHpyGQwZGrz8oBwpV+t31mLkh8zQWMXWok9FIZfBuMgEzWcD1/5vNc5OVSqiPhY55LQoKQ1poqGJhnN+neLtb8YwFiIG3Fq4eWpCCiHd70iUclLycwPmMjjo+D7nvC7bMiiiRrg4NtEuUIVsVGwcWpgxTTQFqYuQabuOociVX2SN9fhYBTP7Cp7bi/EPjlWMGi48FSYaZ6cqFVGbAEuVGqZNtN8vaUkME6THW0AEwqlZW82Y3S+9ysziLHXg+jueD45sXdOE8QjD2kSNcKF1+62a6KzA6KSQy1oCPmdpdEkIL5RNBscNoe2FuLaDpYqRfJYuDb6VJkB739f4kqnihAR8CghTBrW+hZsPE42HBu9qgzdNQd43oFMDivqGFL9fNBtb+HOyei/9M2YlQ8vElYmlemGzyOevWuONTTTGta1oHMV81oyLTzriOL01+GjrwpflbNkYC5rFCQn4FBBGy61zsuYba/AqYdfbY0xw1ygav04w9+qMYbFVGDQbW3jvv1rjjdP7FeV3vR4IrcQpcBpp8M5Kk8WIhWJciOPs8YhBj/phJT/44yxoFick4FOA7EhTZWR64XSy+hFQKiFhaTDq3zdt8F6mH4UTLKobknNeVyNc1A1XIUxFY+Oars02KNDl/NvrOFuJV/N0Fc6icmkoHwHo86snm0Em4+6vKeSijWySTXfm/CcbPBE1YUrUCgFn3fjMeO0tiJ00cjKV636nHhHGZmt7F5FGVK1x1LjdiewVJil+V4Q/NirQBejBMvmsOH/J0ODrGrc0ePA42/kJU5PW5nILjWi0ygJi0uBbUNAsTkjApwBZSDWrQZSrNZtgEjSylbu95+5kta8U1PvNoqJx2z6iikpw+gD01YK7k1UIgkGH81SF0O6LuawZ8ZNUDd5rJQLYj7NoOFkB1Dm/k0YjPwkQrZOVc66Hz3Z6HDzRfsJUMBR2aVXjBzdUN1IumzG66Xg7WRvZSAHYok+iuiHNFUTeocE3GO9xR3y7CpXfQj5/7dR9nQKnoZZr69yVlWzLyTY9iAgpL6KMCHI6da04+GSfJydRFBsjYsZWg7xZAV/RPE0uKtyEhIhMcRtjLsPMmifK/YqsUKn2SVQakahvL9tMSxXN1SYrzql42DRyDju3kR8I7YymadrJ6ujcJaxlSTc9lHyYaIr56OLg3fJHkn6enJAGnwLCmDTckni8wunctPBCPuvqkGuULARYWtDxMavscVQakdMHUMj70+DFw8bbZFUfWiofazu1ulKlhgyD2XO3aSdrShJ4GnWqAqLW4O1RYXF2jYoT0uBTQFgnq+rG8DKluEUqeHVk0p1g/oSLTYOPekktO1kr9jBJOWtWHIfpZG1Sg5cfCO20y8rJOGM174QtoL5GPzcMTEnXTJ29ZlWIMgxRZEc751M+y/S6RQk/T05Ig08BNht8RBp8kPnv5cTSCzP5i8GWbfCxO1krav+FeGg6SxCokPcpsGvw7RTweoE18Uz2WokAKg3eO78hKfhZIUYZ6SLmhzifjLHIM2VbAQn4FFAOY4OvK6QV3F4shKaKklG8ygtxgx6PQYO3wkEtbbtUqdlMSrIQE7Z5ER7oJy/A7qCUBHwbl+3CdCGuamNHpP0Y0uNkVa9EZYoRPqys2kr2VRvFwRORE84G33jZ7hdvDd6/eSBWE42UWl6vwWvKv/XtGzudnX1rnb/dDpzx4Y1jxe2rEMu2nGzNtDkNPrwQtmo42c9X0s+Tk66zwT+8+QAuOnlGqESVx185jJNn9wEAfrFhLy5YOgNLZ02ybfOLDXtRyGXwhhVzXPez59gYHti0HxxALsPw386ch2kTewAAu4+O4uHNB3Dlyjn4xQarnZxKmGzcfRxPbrfa5mYY8MbT52J2fxGlSs1huwxumxSRKU7WvnwYBwZLvm/AX0rdlX72zB786fknYt32I9iw+7j5/lkLp+DshVN9j63OyZrLYKSs4bndx8xt7np8Byb36uf38Mi47ft+nKzZjKzNWX+PazVsOziMBdMm4Cfr92D5nEk458RpvscehpIjfLCxlmvf1oqiiV4zXb/jCBZM7bU1HGmG4XIV63ccxalz+/DC/iHznnNDHPuh4XFsPTCMC5fNaOr37ntuLw4M6l2cdhj9V53nVpynXUdGcXys4to0Jyl0lYB/8tUjeO//XYe/vGgxPnn1aYH2odU43nX74zhtbj/eetYJ+Pz9W3DFabPx9T9bbW5zcKiEm+5+BgDwwmevchUe3/jtK2YfUwAYq2j4q9fp/SS//tgruGPtDrw8MIKHtxwwt1EJ+E/fuwlP7zxme29gqIyPXXUKylXN1qPzHavn4+EtB7BiXr/ncc7sK+CkmRPx+CtHcO3Zek/1Qj6D4bK98XdVq+G6rz8OQO925IW4WeTG1et3HMXBwRI+9pPn8MqA1dT4lDl9+NXfXOy5PxlnotWiGRMxrtWwbofemi/DgK8+uk353d6erGeZ3UUz9K5OS4x/Ab2tn8xtD23FO1bPxy0/3YhJhRw23Xql77GHQThZz1s8Db9+cQALpk7w3D6XzWD6xB4cHhnHrL4CDgyVjP1Er5le++9rMWNSD9Z96opA3//Yjzfg/o378Y5z9Pm3aHqv5/Zifr3vrnXYdWQMT33y8oblkwWHhsv40Peesb3Xk81g3mTrfMor2C89+CI27xvEgze/zvfxtIOuEvBHRvSn8/bDwZtQiSf45n2DuOI0vcWbU6sdKVuvSy5x6Pp2VczqK+Chm1+HMz/zIEYl4Tls7EP0hfzLixbj6799ValpjZQ1XHrKLNz2zlUAgAu/+ChGx/XtnAkiV66Y46vR9VOfvLzuvUIug8PDdkFQkgRDM4Wu/uKCxThlTh8+9pPnMDKuYaRcxVvPOgF//+YVuOWejdi457jHnupxlmT4k/NPxDWr5qHGgYk9WYxrNVu8OmNAXyGHwVIVeoNtd833zWfOw6WnzMLEHmubi5bNxPO3XokJ+SyuuO03GB2vmtfd+RCME+Fk/eb1a1CqaujtaXxLr/3EZahxvZqk8IfEZYM/NDzeeCMXth4YBgAMDJfBGHDzFSd7bi/m364jer/fZsyZonnLZ9+yEm8+c565P3vUkRWGOVSqYrjUuusclK4S8FGgispwpnn7DWssV2uY0JPF5N58nX3bGeUxY1Kh7vflbScVcqamLjtD652swSkonKyy06mheUAya0yekMcks8a8XlKgr6gfw6RCrmmB4wxrA2BrJu6WgCWvbrwQPWZlJhrvTe3t8Yy5j5NytYYJ+SwyGeZLuANAj8JZHHUUTTXCOv9j4xqm9vY0DH10zvNaE2W2xfyZ2pt3nRNyS0A/lUqTADlZm0QlhL1K4HrdOLJ2XS/gRZy2LuDNjkoKIeJ0QMnx6o2aWTSDyskkj7lxmnx9oStAPw+2c5FvPhzNEvCtLwJmJVW1/oYPe33jSuCJspHL6Li/Y3Ru00x9HT8tAQt5a/77bSbTbrpMwDPp/8FQaef1JXDVkRuqfTnjtp37GDSyPvsNbVStwddcl5J+og/8osoUlF83W+hKjLlU0erPRZMapZ+esHFRNCN22qPBhznmuFLwo1gRiPt0sFTxtQp1huk2c0xW3Lt3jSa7Bp/8kMmuEvBiyRamcogqZNF5oWVtyCusyt6gwq4di30IE82EnozRhk4h4B1anNw5yU+GqV/kJapznOIYPL9vC+ez0uRHxquKUr9aU52sVHHLraKQ11dMsj+iVeV3y1Wt4YPVi9gEfAT7E2fQTytCoP7aN7Mq8bMCLOatlXG5UktFmeWuEvBRPHFVhb+c2opfDd7eoEKtwcvVDlUCFhBJRg4TjZGyXa42zjD1i7phRzAnqxyD7azoWMhlUOPNFfEqVzX0KKpmtoJCTq/R4xZzHyflSjgNPpfNIOdRJTQoUZp8/Ap4p5YfRINvVGVV7LNUVSt3SaO7BHwEy0Z72QDLHmffxr+T1Vk7xdp3fSKOqmNNVatBq3GH8NQfFhWNg3PvwmLNINf6UI2zkf07K9W4Ed2EAMsMFab/pZ9ysnFh1b1RZ83GiZ866Y2II4EnyhUB5/58K87r38z88dORzLYyrqiVu6TRXQI+ihoVCi2tToP3qcnJtnOndux0EBXzag1eFT1SNJx+JR9aSTOoan3YbPBNaJKFXMY0LZgVHUP0v4wyWqhZxLULU7c/KFE40VXRUWGJen9+Vil1Ar4J4Wva4Bt2JLMrdUl3tHaZgI8uhRmQNXhnbLg/TU4uIyBrB4CiFVsuo4wucRbZ0rfVM07LEdulVYWpSk1o8LZ95TPmTess+NWo/6sKP6US4kI8eO3mu/iX7sIEF3aF5lUlNChRCz5fNnjHeWim12zZjwbvCD/W/yUTTWKIsggRYF3cca2GWk02W0gC0NMGL/V8dAmTFBRyWT1ao87eX6+lO2uhR2mikX/TOc5mBGxRcrIed7TNE+ekWSdZewV8zeEkj1+zE6GIYY9b5VsJS9QPDD9zOJwG78PJKvm2xNxMem2a7hLwxkUM4/l2E25y3K/KTq/clw8nq0BovPWJRvWah7MWelSCz2pbpj7WZpx9spPV6qoUvP+lXBe91RTyek2X4bJcRC1+zS6q6+tcPUZB5CYaH8eYN9pKWmOI2MlqXOfRcQ21GGv4RElXCXi3sMZmcBNubsvzhk5Ws+ejOpNVUBRRNA6hZ9nZnXHwlskgKsGn1uCtv5sJ15O7CdU5WQP0vxQp++0gzjLIXlgP9whMNBGPN2rN1u8cttfpD+Bk9YyiMeZrDD2F46KrBLwZ4hRi8pVd7OR2Yd/YyWraTyUNXuxPq/G6EEE981NVKqB+YoqoCGcJ3bCI33AzRTSrwWcyDD3ZjDJMEmju5nHrXNUKhIAdlFsRtmDpLuZL2DBY1bwKi7y/ZvIZZGSzp99VirxdM/d5udpMT+HWXucwdJmAj16Dn5CvdzyWq5r5vtskcy6vZQeqGN8ESTPryWbqtHz7furDJJ1NMMKi0qzlv5t1soqx1dngA6TPh83oDIOs2fX2NL/6CIr1AA+pwQcoDdEIeX/VgCbRIOY/e53+aMNszQd5DB3J4qLLBHz42FVntqkoTOQ0y0wsZJH1SCCpE/C5TF1srVz0KJNhtkQLaz/1KdYiUWi0LD6LRrMtKmzj8t9NOVml5tjihilKQh9oNoqmvU5WQA/3FCUlWuF882M39kNB4bwPS5guZKrv+TX/yfdBs05W3x3JRjvIRMMY+xZj7CBjbJP03jTG2EOMsZeMf/13ZmgjbmGNTe1DEtij47KAtws9PerFu0k1YGkFwm4uTDdAfaVDlTNMVSRJCHSnZhyWQgMna95jeVu/L0tbHzIbX4eJg3cvyxw3YtxD5arygR/b91+UAAAgAElEQVQXkTlZ85mmQgr9YL8fgu3btjoMpME3OX98ZmLbfS3p1+C/A+Aqx3sfB/AI53wZgEeM14lHTOIwyyrnpOmfkKvbZ8loo1bwCD8TgrkoafAiPV/sS+xboFpKm05WRds2p207LJYN3j1e3y89WXtIpPy3FW/fnJOsbRq8MW7OrWvWCs0uKid6HHHwtpVuUA3etjr0d4w9Nht8k/OnkQYvMq87ycnKOX8MwBHH29cAuMP4+w4Ab4l4XLEQtQYPSFUeFRq8V/iZ034qp+eLz/qLdg1ej4P352QF5PDDaJ2sUWRsCmeW03cAWDdSsxpYu5yssvnA1OBbYJu1ktwSGAdvS/Zr/lzUatwWeux3DmdChEk2+g0r81p2sqZfg1cxm3O+DwCMf2dFN6ToWPvyYXzwe09jw65jeGVgGGtfOQxA77jz/rvW44PffRobdzfbOcg+acQN/ZWHt+LJV4/ovVS3HDDivN2dV04HqCzUvvOHV237Fjg1eM45vvTgi8Z+6gWlmIhRNvwAgG/+7lW8/671eP9d6/HY1oFw+1Q8mJpxsg6VKvjbH20wOzO1A3kVIh7KP356N/767mew60jw7mGN8FPD3A9uReyCMjau4d/+62Xz9a2/2IwX9w+5bv//Nu7Dw5sP2N5z1pMPcm3v3bAXBwdLntv89qUBvP+u9Xh65zEfTtb6cFiv1YlW4/ifP9sUyX0SlNg7OjHGbgRwIwAsXLgw7p+zcc/Tu/HL5/Zh3uQiZhuNf89dNA2DpQpeHhjGSweHMX/qBJw+33/jXCFgV57QD86Ba8+Zj51HRvH0zmP46TO7ceb8KQD0ptF/2HbYt5O1mLPs249tPQQAePvq+dg2MIwLl84wts2iWuOoajXkshkMjlWx97g+gadP6jH3XXSaaCKKLpk5qYCLls3AgcESXh7Q26lNLGRx0bIZ0GocZy2c0nAfX/vjs/GbrQfN129YMRuj41UsmTHJ5nAG/GlgG/ccx4/X78aJ03txwdLmmixHhSwYRGOWTXsGsWnPIM5bMg3vPu/EWH43qhr4UTtZNxs9dycVcjhp5kT8ZusAVp84FcvnqJtmv/+7TwOArZWkGM/U3jxOmDoBqxY0nlsA8JZV81AxsrhfHhjB2lcO45pVJ7hu/6N1u/HIloNYNKMXb1gxx3PfynwHj/O268go7nx8BwB9ZXHxyTN9HUOUBBXwBxhjcznn+xhjcwEcdNuQc347gNsBYPXq1S0tniyernIhqDvfe66p8Zzx9w80vTQtVzScMqcP9910kfneBUtn4MIvPmr7nZsuXYandxx1t8Gb0S92s0SpomFcq+G6cxfgtSfNwL0futD8jiz4clnLMfb5t660OTjry/BGI+B7chncecN5ofZx9RlzcfUZc83XH7hkKT5wyVLbNs0IeHGDfeWdq3DWwvb4+mUNekJPFvksM/MYmukq1CxRZrKKchuyiSP4uPR5+Y3rV+O8xdNw0i33N32fibn9t1cub+oB+ecXLMafX7AYu4+O4sIv/rrh75YqGpbMnOirwbsz8xrwdrL6zWiPk6Az414A1xt/Xw/g59EMJ1qEfUxvCaf/3SMLwgBV9NzCqYSjSja9yMlLqv2I7eR/xVhVy26n4HNboss2eMbsx5wGGGO+zQZWqGB77O+A3QauZxxbY4kzXFLM6bAmOKFcRNVmT/YL6dcywH1mBiEEOza/BeuaKdYmh8PK33ej1cXnVPgJk7wbwFoAyxljuxljNwD4AoArGGMvAbjCeJ045MQhUYxKbggRJHrAzRkjJrEsuL0SSMqOCAink1X1G0VHopFbHLRVZ71Sd8xpwe+1iTpbNwhuvWaBeG/sqDR40zwY0cPI+dAN1mM3nPnJaibuff79OFcFRTPRyXKyevmJ2lE+2klDEw3n/DqXjy6LeCyR4+yA7tR0gkQPlCo19BXrT5tcwTFrpDwXchkcHm6QyWrWX7HKALiuEvKWli/vw3lcZlx2qdpWzTYMfldXpqbXphh4wBHqmcvYxhLnje2nAqIfLAe/BiDvvbEPnNE9qiqofvcR9Nj8mvnK1Rr6iv6OWVacAN3H4F1rSvO1XZyka+3eJGJ5LJo612m6AaIH5PoxMmISlytWizwvIeUMb7SEsrvd3IwuccTzO7eVnaztiiwJi18NPuqmJkGwRwJl7Rp8zCYaxoB8NtwKzYpaithEI/mXmk2kCltmoyebAWONNfhm8iecTtbJE/LeAl7KSG9XOGU6736fWCaamrJ3pSjK1dw+NeWyUZTylbVvr/3XOVldaqM7xwvUa/D1Dy59n8PlaltNF2Hwu7qKuqlJEIQwEeOwJdvEaKIpVWso5rKhTXCqKqFhcD50g5lCw11Xy4/TSIP3b6IR13m4rJto+oo5T8EtfruvSBp8LMgmmpIiESao80fl+BHCXG6hpqodI3CWJxVat7DvqUw0znrsbg095Akb1EnVbvzWKE+Ck1UIE8BoZCKbaGLW4KN4gKvq/IdBtTp1u8+qLo5dpwIUBK/7z/ydiv+etvJ17snpmepecfBWRjpp8LEgN8VWFaMK6vxRavCSk9WmubjGwaudrFZ/Ug8N3tEL1pnJKAuYtGrwfmuUl6s1ZCIwU4RFLnUsDyVuJ2sUK5eoNXincFZVQRXIkTv2Zu7hned+TLDNViGVr7O+MmmswTcy5cRJOu9+n5RNW7WLkzWI88fFZicEkpwy72VmKFdrNvups86FOhTTbit1c0Q5bcJpxO/qSvhE2h0p5MxIFsTtZI3i+qrq/IdBPHRzRky9lyZtL7Nd/3eY4yvkG5uGyhWtyUY11qqkkRIi5m9/kQR8LNg0eJWTVdECr+E+3ZyshjDXHySWBj9erSkbHpQd9lM/yUn2aAd3R5Qq/T9teGl9MqWIzBRhcRZKE8Qr4KNpNB6kg5YXzoeul7nN1qBeOldR9DIo5rINfSBNa/BSB7ZGfiJZg29X3fj23xkxIicEqZazzTpZRQEkVXEnsVyTbXpeRbOc9lNnEkUYJytjzHT0pVXA+02fb2cdeJlizjJHyMR5Y5cq9avSIDjnVVicD13vfBA527Ne2IeywTfQ4MX93Mz8qTPReGWyGr/dPyFnlhhpNe2/M2JCrkZXrgrnZzgnq9XFXh2jXhLO3Lqyt/UX1hmeZQl49wJhbk5W5baS0y+N+A2ta2cdeJk0a/BxOFlls0fR4z5zSwaKoplJIyer1/3sRlG6zo2UkFJVQ082Y3ZmiypTuBk6VsDLJ1OPg6/XvP2aAeT9AO4x6lqNY7Rs3XRFh0lFxlneVnjovUr8WrZSu5NVbc7Jun6WBprJZE3CMTr7yQrijJ5Qhf4GwTmvwuIMRPDSpN3aP0YR/upXw26m3LK8Om9k4hWry6hXSM3Q/jsjJsRk7TOyyNw1eP8n3SstXq7/Ul9+QGGiUTxwCrn6BtTO8cr7K0lZs27jSauTteiRJCbTzmbbMiona18xF3uxsSidrJFp8FXn6tS9JpNbA3evue2XRiZYK14/iJM10zBIQzzohLIVZ06EGx0r4MVkFeVbh0v1ST+FXAZajaPic+nkVX/brFMxVrE1lAbUmpHq5izks54lflVhkm4ajltUR1poJkwyGRp8/QO1vxivc61U0UI3+wCid7I6H7pe19JLgw9dY6eBkhBklWDPd2gURWM1/pF/r5W0/86ICWfj6pFxdRQN4F9zUTW4Nvdl7Fv/HftyXbV/OSFKUMxnMDrublfPZBh6stakUoV+WvvK2v5NG0LrU0UgyXidg1Yix3yL6JG4Mxij0uBVzdTD4FyditaVymgylwbuUVxXv2GMzShBcuZ5IaeXWdZq6jkqInSi9nE0Q+cKeBGDKvU1VZloAP920pKHBq+KiPHSjFThWU6tR4Vsm/ZysqU/ikbvUVt1uXkEUTkaw2JmOGYzpiCLoxWeTFTHnstmkM2wiE009XNZ5WT0crKGPbZGTtAgsfb2OHjDeeqara7ZNXgy0USHEMZyX1OVkxVoXoN3c7Kaf5tV9Nw1I1UPUT/x63J0iVehpKwUg5xG/K6uShE5GsNSyGWRc9iM/WRShiFK85Tf0hB+cApnr0Qqt5rpUTRRbxSJJX47kJM1l/EMotDfF07WaIu5NUPsLftaRVWr4dr/WIuqVsM9H3itedLlvqZuGvxVX3kMhbz+pP3G9atxypx+XP3V3+K9Fy3GW8+ab27/8Xs2Gt+rnxC2JanpaXd3rpQrqqgeWetRaxWFXBY/Wb8bD20+gONjFSyePlG5XS5rT6BKG+Jc/PK5vbjtoZegcY5zF0/D1/74bADAP9y/BT99Zg8OD5dx+gn+Wy7GRVFaik/t1dsn9hVzqGgc53z2Ifzgfedj6Sx7y7r33bkOK+dNxk2XLQv0m+WI4uAB/Xzf9cQO3PfcPszo68GBwTIA4C8vWowbLz6pbvuDQyVc/dXf4XNvWYnP3rcZF5w0A198+xnmuOSHrhjjlx58EZ+5ZqX5vlbj+OiPNpiv/+qup7Foei9GxzUc85jbzRzTsdEK1nz+YWQY8JlrVuJKqS1fEA3eDJPMW4J7zecfxsJpvRgsVXHOwqn4jz89B+t3HMVjWwdw7qJpdQmKraRjBPxQqYoNu44BAI6OVGxZZAKVkxXQC3ydOrUXW/YN4oV9Q1g6cxKe3zuIm3+wwSbg9xv9T89U9Idcs2gabrhwMcpVDW9ZNc+2f2UcvEKD/8AlJ+HhLQexZMZEW0VCmb+98mQ8+epR8/VFy9R9SG+6dBlOnr0fV6307jOZVMS5W/vyYRwcKmHprEn4/bZD5ue/33YIPdkM3rlmId6xer7bblrGu88/0WwZ+E/vOBP3PL0bFY3j4S0HcXhkHHf8YQc++5aVtu888PwBPPD8gUACnnNuzKFoNPi/u2o5nt11HD94aif2D5Ywq6+Aao3jyVePKAX8I1sOYmCojE/+dCMODY/jB+t2mQK+VLWn/19+6mx86meb8MI+e+NtUZVRZvvhUfT2ZHHt2fNd57Zfrj17PkbKVdQ48MN1u/DMzmMOAd98rP1/X7MAGud415qFZiPxisbx8sAIAOB3xhx9fu9xAMA71yzAnP4irjt3IWb1FUMdTxA6RsA77XeWDV7W4NVOVgD4o5VzsGXfIMpVzTUhoVqr4X0XL1FqTX3FPP7nm06z79/DyaqKErjs1Nm47NTZyt8WvPWs+baHjhsXLJ3RtibUUWD1v6yityeHS5bPwp1rd5ifl6s1rFowBf/4ttPbNUQbJ8/uw8mzdQ192sQevPeiJfj271+N7fcqGgfn1ioxLO9csxDvXAPc++wejIxrOHVuP4bL1YYmMqX50aHBz5lcxOuXz8SRkXH7dsY9+pErTsaXH9pqvj+7vxjJdV0+pw+ff6u+n/ue21unQQcpaLZi3mR87i36PrcfHqn73FkI8MqVczCpkGvbPG2/8TIiZFueXrZX2OCtZ5iq2JhAbCd/V4Zz3rRd0NliTyYp4X1JRZwb0bSkkNPtqcKBqYpCShryfOOItt98FJmeKgqOKBE3u7wIiFGaH1UhwLls3X3ljHSzto3+uqp+X4w9bN9XmYrGodV4JLV0oiDZd0gTuGnwk3v9afBC05e/K2OmNTehMXk5l5KSYp9UbHkFhrDhXL+BgOaLRLWDOG9uK+ku2jkkR4D5yUUQ18M+Nn+F/VSRbvIYokTl8A5bktjNOSvKhssVNdtFsu+QJrAnTNSU2oGbk1XeTrTdq99/gKQIFw2ec04afANsmcH5bH3D8YRksHphK0UB+41eaxD+2Yi4NER7nHdz5bQ559ZKV9GExvmwEIrPxB6HgI9B8VGVJQlbktj5PZsMMRS4dpex7hgJY9PgjeqRgD1MUlVNUjCxkDNjgd1s5kAwDd55k4xrtUjtp52I0KpkEw0gF1pLtwYftvBU2JZ2bsh1lIoetVZUcquicamAV32EmJuArW8aH4+JxnkfhjVzOb9nCviqunptO+gcJ2vFbqKR22UJGiUWiWWcurRA85Mh75JAEtfN2UnIscO6uUC81qQVULIfkF4PoLAx515lM8IgJ/JUctw1dltV+la2x6uUKecxl11WIXFcV5WJyNk2M8g+ZYSpqVTRElMjqWMEvD1hwk2DdzhZHbHrYhmpLi0QbDJ42/7aPwGSiv3aZGyJT2l5QHr5WMJmjQZJs/eDnMhT8dDgvVa5+rgaF/Zz8yPEsTJTmYjKVa0uOa0ZnM5ZIWtMDT4BK8z2jyAinN1gxMWbWJDqUnto8EUjcUF8t37/wZo7qyraxRUB0UnUr66s2vppEfDy+JxRNGHrvpjmjYi1xHzOSpDzLhKmDhww6zUpNHjN0fTC7T6Io1F8MZ+tK0kStqCZU4BPlgM1HDXx20XHaPAqJ6ssGACFk9VREElUh/N0sjb5VFY1FQm7NOwGnOVm5WxAS3tt/w3khTzfnAWpwmY1mk7WiLVE4Qwu5DOo1rx6qarDI0X4pJtWXq7WTI3Z7UHdOg2+3hnc3D69NPhktJJs/wgiwu5k1YzOSlnks8x0CDknUo+jboio7+xWWgBoXrtQe++DrQa6ifrVlRVyal2LZE9feb45lYawdUniWsWIlUbR8Em5ldMuuVRIddXK85YPRd5e/kwQm5O1LoonXC5FnZPVCMnWbfDJcLK2fwQRUWeDN06w6JQE1AvUjBSjKuy8YnnlJKjNU+29Dxd/2w0U6vwjVphkXDHgUSMfg1NpkDV4v/0I7N+Px8kqkB3bbj2FVWNy1coVWd1u2/aEaPLhhl6F1XkNwtXycY5bJEuKMMkkKHAdI2GcE0e+eGZsr4dALea9naxBNSZVRbugq4FuwmaiycuV+2qp8WHI19ftIe/82y/lAJUQ/WCaaGTHtoswrx+TtLpSOFmd33MNNoghdFwdphlOg3c6Z51hklFfmyC0fwQRISYLY9ZS0dLc7f+qyGWYp5PVSixp3snqGn+bgAmQVBhjtpr2soBIiw9Dvr5uZjpALUAbUYpbg7fVMfeu4y4WwuWqJrXBc9Pg68sD113HaKs6mL9RF6YZcay6PRs+GWG8yb5DmqBc1ZDPMvQaT2rZBlbIZdGTy3hmlQlTjixA7PsPGiZJTtagyNfPssHLdt7230BeyKYGtzmgfxZcg49aSRA2eL3eeb3dXCC/JwSb7B9xC2hw9mBtVTq/qr1e1LHqwsmq17NKRq2k9o8gIsz+h0YfRtkGJmdCeiGWcfLNKMK6gtrNyckaHCEgi3l1HHwSlsBe5LL+NPggSU9xh4r2KLKHVb8POMIDXVanqu5p4h5tRTp/IZdFtS5MM9pY9ToNPgHzM9QIGGM3M8aeZ4xtYozdzRhrfcFjA/HELBomEfkEF6VaJl6IpCSVfbTs4vFvvM/6jvJpEVDthjErJtusRVPRXLXEJFOfxRlSg6/WwFj0Dklhg2dMbTe3ft86HqftGVBUblV06HK1U8dig9d/Ry4RoQdiRDeHbOchIZmsgWcHY+wEAH8NYDXnfCWALIB3RTWwZjHbYwktXLp4fjX4Qj5TV2xMrn0i9tUMbvG3+mftnwBpwFmLJo0+DLc5oPrM3/40M0osLqxaSu5hw0B9Bqf8XWtfaidrq+4BOVHO/P2IY9XtUTTJ0ODDJjrlAExgjFUA9ALYG35I9RweLmOkrGH+1AlgTB3eJN7ryWVwdHQcw+WqeYILeZ8CPpfF2LiGwyNlab/65N57bAxA8xqTbPs7PlZBrcax//iY+RnhhWEPzmfM87732Jgvp3nSkBul92Qz2Gd0B9M/82+i4Zxj99ExHBoaj0U4Chs850ChR7Sc9K47I2qwiHkNuDtZdx8dxQ6jUcbhkXH1PRCTkxUAdh4ZxRQjXn3XkTFld7agTOjR824ODpUS42QNLOA553sYY/8HwE4AYwAe5Jw/GNnIJL7y8Eu48/Ed+B9XLkc+y/AP97+A9Z+6HNMnFcxtyhUNPbkM+oo5/PYlvW3Wa0+aDgCYMqEHY+ONb6K+Yg5D5Sp+/qz1nCpXath7bAzffWInADStMelx8Bq0GseZt9pPD4VJejNjUgGHhsfRV8yb511cBwCYVEhPIrawTy//1K/wrjUL8P2ndpmfqQSoG9/5w3bc+ovNAIATpkyIfJwr5k3G77cdxoxJBdSMtNRGGvxM4z782q9fNt9zKmB9hpZ/6y82m+MHgJUn9AMA1iyaiqe2660oF88I14tVhfj9a772e/zLdWdhz7ExjGs19EU4h3T5kzfnaJT7DkrgETDGpgK4BsBiAMcA/Igx9iec87sc290I4EYAWLhwYaDfetvZJ+DOx3dg//ESntp+BACwf7BkE/AlI+34C9eeYfZmFS3rPvWmUzHuchP94eOXomo0LbjhwsU4aeYk1DjHKwMj+Ndfb0OpquH4WAUAcNOlS5seux4HX8OY4ya57tyFrn1XCZ1/e/fZ2LxvEJcbbQxvunQp/uXRbQCAK1fMxhSjuXWSefSjr8M/P/ISHnz+gKlkCOF+9sIpeHrnsaY0+P3HS+jJZvCFa083WwRGyf+4cjmuXDEbp83rx7aDes9RNyfr606eibedfQIuXDoD4xrH3U/qgu26cxci71jpzplcxB1/cS4OD5dt7680GqZ/6z1rsOPwKEbKVZy7eFrkx3X5abPw5f9+Jj7yww04MFgy+yt/6NJgDc8Fv/u712OkrKGi6Rr7t96zBq8MDCObYbj0lFlRDD0UYR4xlwN4lXM+AACMsXsAvBaATcBzzm8HcDsArF69OtDi66yFUzF3ctHm2HE2UCgbTtaTZk7CSTMn2T6bO9ld05knaUFTenvwlrNOAAA8suWAsV/Lrnj+kulNj72Qy2K8Wqtzsl0csqFwN7Bk5iQska7l+UummwI+CTePH5bMnIQFU3tRrmp1SsZlp87WBXwTGny5WkNvIYu3nR1Po/F8NoNzTtQFbCMn69zJRVyzSr9fLl42wxTwbnP7dSfPdP3dvmLeFPZxUMhl8aYz5uEjP9xg+nFm9xcws6/Q+MsezJ/aa3u9asEUrIrQ7BOWMCrkTgDnM8Z6mb5+vgzAlmiGVY/TWSmWj4KokxbsiTXBMyfFd4ZK9g7yZH9vHmcBsrRQyGVQ48CIw0woR134pZXx1XLuQf047PebvXBfMue2qEtVNmvFpGcOBSXwleCcPwHgxwCeBrDR2NftEY2rDme4obMjTti6Ek6ssC6p9kmACSHGJMw8gm6YXFHjLCGcFtzmgBw37ZdWOu+slpNqDV6+37yqtiYFkcxYEtUeUzSHghLKC8A5/zSAT0c0Fk+cmWiq9P9YNPhKuLA8MaZ6Ad/5kytq0qAlqhBjHXQK+KLoANSMiab1GryfnsJufycNEfRQriQjjDFuUnOEzpouqi5JUWoOZmq1VF0ySNSLmOzOmzvK1Ua3IAuONEUguT3kJwfR4FsomKw4ePsDqKLxup7C9vLOyb02IrM8KWGMcZMaAV/M22u6qLokRTnx5cldCqPBG5N9sEQafFhsJoEUaV9CkDjngAjda6a7U6mqtezhJgq+OauhqgqKOVssJhVRF75U0boikzw1R+h0sqo0+CgnflGyP1qp8c2frqKpwTucrF2gPURNWp2sxbx6DohGJk1F0bTYtFDMKaqhKsr8psEGD9jLkSR5nFGRIgFvd7Kq22/FoMGHdLIW3JysXaA9RE1anaxi3CpHu6qMrRetFkwFlzrqgOOBmxL/SCGv90juFidrao6wXoO3/tZqHONaPGGSpUq4BhNu9tc02ZCTQo/NJJCe8+c6B/IZpQD1otWCSWi89jHUr2jT4h8p5rKSBp8a8ReY1BxhfRSNFDIZQ/GuHkmDL1Vq6MlmbC3+/GI6WZ02+ARrOUkl62ixmBbMKJo6P0xWKUC9aHWvT1XDGlXzG9vfCZ7boqBg1LXgk0pyr4QDEd4kUNnjo5z42QxDPsvMrLeg+zYdbA7tLY6+k91Emm5O1RxgTE+8UbWS88IZfx43qoY1qt4I8v2R5LktnKxJaakXN6k5QqHBiwRWlbCPeuIXjYbZYez7poPNkckaZDVAWCRZS3SimgNyQ/hmatG02rSgbFijCBuW53OS57bpZK3UEt+0PQpSc5eY4U2GNqFs3hvxxNcfKlqoGHszTNKhwRPhSJWJRqHBm43gA0XRtFiD76CewsV8lpysScRZ00WOPAgTp+79m/pkKIWIsXdzsBHhaEWbt6hQzQG536zfOHjOuT4XW2mDzyvi4FPcU7iQy2C4XEWNp3P8zZKaI3RmhKo1+Gg1G9tyLqgGL42bygN3JyoNXrxXzPt3spoZpG12sqa5p3AhlzGvQxrH3yypkThiSVutGU0IYnayAnokTVgnqzzubtAYiHrE6k/MXQDIZa1+s35NNGKet9LJ6swg18eRYg0+nzWvAzlZE4RzMsmTLjYnaz68xz2XYRA+p27QGIh6VIJQGJj0xBt/Gnw7BKtXT+E05SIIitK564b7MT0C3jGZth8aNbs7rX35sL5N1E7WXAZ7jo5i//FS4MmgR0pYDjWi+xA1XWSELi9swr98bh9e3D9k22ZgqIyfPbMHm/Ycx8+f3YOnXj1ifKe1TtbBsQoeeH4/ajWO7YdGsM6479LoZC2ktJ5RUNrfNNAncycXba837xvEu7/+BF747FX411/rXX7CdmdR/eYTxk31mgDdnASFfAZjRnGjaRN7kB73YPJYNmsSXjo43O5hNM3cyUXsODxqvhZdx4r5LA4Nj+OD33sac/qLePyWy8xt/vXRl3DH2h11+5rtuBfiZM7kIkbGNbzvzvX4xYcuxD89+CIe2zqAvmIOExxKVxrm9px+69zN7m/deWwXqRHwaxZNw5O3XIZKjWPGpB58+aGt+M/fvGIuF//k/IW29ntR8E/vOBMfuWI5AGDelOCTQY8XrqCQy+LJWy6OaHTdya/+5mJwHqjzY1u5/68vwuHhcczqL+Do6DimT9SVEXlVt3+wZPuOswsYoLeE82p9F68z2w4AAAoWSURBVDUfuOQkLJzWi5vufgZD5QqGShWsWTQV37h+TV3f1Selh1NSedvZJ+D8k6Yjl2Ek4JPGLPnp26f/PVTWPeJefVeDks9msHB6b+MNGyCWgoV8BrkEZ/mlAb1cQdL1xHomFnKYWNBvN3mueplbVM7XhdPCz8dmYIxh/tQJ5njKlRqmT5lg1rKXScPcZozhhIgVwSST/CvigrNLTpLt21bMc3LHSLQHrzmhCp9sxxyyuptpoXJCiNaT2ivlLMGa5LRjK+Y5uWMk2oOXsFS18muHcBW/KTR4UlTSQ2qvlJVAVLW9TiIixDLJYyTag/OhL/sX1Bp865UEubtZtzTK6BRSK3GKjlZ4SdaOrTDJ5I6RaA/Oh76cDKWywbcjOcfqbqYZ1SxTKza6jtReKWd9jyRrx2SDJ9xwPvRVJTi8tm8FYt6WKrVQZTuI1pNaiWMK+NEUCHgpioYgZJzz1t6WMllO1rGKFnnnNCJeUnulCg4TTZK1CtPJmuAxEu3BaVqUNXilk7UNwjWfZWAMGBL3GikqqSG1V8rZqT7JdsEiafCEC3U1lhpo8O3wNTHGUMxlrXuNFJXUkFqJY5ZgTZEGn+QxEu3B+dB3ayzvtn2rKOQz1r1GikpqSO2VqnOyJnjSkZOVcMPTyaoS8G1SEgq5jBTQQIpKWkitxEljFE2SQzmJ9uDmZK1qNWi1+po77ZrnhVw2FfcaYSe1VypVcfBS/02CkHFzsro1AWnXPC9KJpok32uEndRKnDRlspommgSbkYj24OZkdWsC0k4NPg33GmEn1JVijE1hjP2YMfYCY2wLY+w1UQ2sEblsBtkMS4Vd0NLgkztGoj242eDdNPhk2OBJwKeFsOWC/xnArzjnb2eM9QBoaS3TQi6D0XENjOmxukmFnKyEG85VndDc3QR8uxq3y+NMcmE/wk7g2cIY6wdwMYBvAgDnfJxzfiyqgflB2AKLuSwYS66AN8dJNwbhoL7XsNDg1Y3k21VyXY59T3LOCWEnzJVaAmAAwLcZY88wxr7BGJsY0bh8kRbbNmnwhBtOxUQI+O8/uQtAcuaMTYMnU2NqCDN7cgDOBvDvnPOzAIwA+LhzI8bYjYyxdYyxdQMDAyF+rp4Ll87ArL4CLlw6I9L9Rs2Kef0458SpZh9OgpA5d9E03PrmFQAszf1Xm/YDAP7PO87EnP4ivnrdWVg0vRfzp7a2o5PgvMXTMauvgFPn9mN2f7S9j4n4YEH7WzLG5gB4nHO+yHh9EYCPc86vdvvO6tWr+bp16wL9HkF0MpxzLP7E/fjwZctw8xUnY/XnHsIbVszBP7z19HYPjUgAjLH1nPPVzX4vsAbPOd8PYBdjbLnx1mUANgfdH0F0M4wx9OQyKBkafLlSo5ovRGjCRtHcBOC7RgTNKwD+PPyQCKI7KeQyZg34crWWeN8SkXxCCXjO+bMAml42EARRTyGXRbmqlyiguutEFNAMIoiEUMxnUK5qGDciaShahQgLCXiCSAiFXAblas2MpKF4cyIsNIMIIiEUclmUK5oZC08aPBEWEvAEkRAKeV2DF+UKyAZPhIVmEEEkhGIui3KlZmnwZKIhQkIziCASQsFwsopQSYqDJ8JCAp4gEkIhl0GpYjlZSYMnwkIziCASgh4Hr6FUIScrEQ0k4AkiITjDJMnJSoSFZhBBJIRiPmsI+Jr5miDCQAKeIBKCXotGIw2eiAyaQQSREAr5DErVmmWDJycrERKaQQSREAq5LLQax0i5ar4miDCQgCeIhCBqzwyWhICn25MIB80ggkgIQmP/v2u3G6/p9iTCEbbhB0EQEXHJ8pm4ZtU8VLQals6chFyWBDwRDhLwBJEQTpw+Ef/8rrPaPQyigyAVgSAIokMhAU8QBNGhkIAnCILoUEjAEwRBdCgk4AmCIDoUEvAEQRAdCgl4giCIDoUEPEEQRIfCOOet+zHGBgDsCPj1GQAORTictNHNx9/Nxw509/F387ED1vGfyDmf2eyXWyrgw8AYW8c5X93ucbSLbj7+bj52oLuPv5uPHQh//GSiIQiC6FBIwBMEQXQoaRLwt7d7AG2mm4+/m48d6O7j7+ZjB0Ief2ps8ARBEERzpEmDJwiCIJogFQKeMXYVY+xFxtg2xtjH2z2eqGGMfYsxdpAxtkl6bxpj7CHG2EvGv1ON9xlj7KvGuXiOMXZ2+0YeDYyxBYyxXzPGtjDGnmeMfdh4v+PPAWOsyBh7kjG2wTj2W433FzPGnjCO/QeMsR7j/YLxepvx+aJ2jj8KGGNZxtgzjLH7jNfddOzbGWMbGWPPMsbWGe9FNu8TL+AZY1kAXwPwRwBOA3AdY+y09o4qcr4D4CrHex8H8AjnfBmAR4zXgH4elhn/3Qjg31s0xjipAvgo5/xUAOcD+KBxjbvhHJQBXMo5PxPAKgBXMcbOB/BFALcZx34UwA3G9jcAOMo5XwrgNmO7tPNhAFuk19107ADwes75KikcMrp5zzlP9H8AXgPgAen1JwB8ot3jiuE4FwHYJL1+EcBc4++5AF40/v5PANeptuuU/wD8HMAV3XYOAPQCeBrAedCTW3LG++Y9AOABAK8x/s4Z27F2jz3EMc83hNilAO4DwLrl2I3j2A5ghuO9yOZ94jV4ACcA2CW93m281+nM5pzvAwDj31nG+x19Poxl91kAnkCXnAPDRPEsgIMAHgLwMoBjnPOqsYl8fOaxG58fBzC9tSOOlK8A+BiAmvF6Orrn2AGAA3iQMbaeMXaj8V5k8z4NPVmZ4r1uDv3p2PPBGJsE4CcA/oZzPsiY6lD1TRXvpfYccM41AKsYY1MA/BTAqarNjH875tgZY28CcJBzvp4xdol4W7Fpxx27xAWc872MsVkAHmKMveCxbdPHnwYNfjeABdLr+QD2tmksreQAY2wuABj/HjTe78jzwRjLQxfu3+Wc32O83VXngHN+DMB/QfdDTGGMCQVMPj7z2I3PJwM40tqRRsYFAN7MGNsO4PvQzTRfQXccOwCAc77X+Pcg9If7uYhw3qdBwD8FYJnhWe8B8C4A97Z5TK3gXgDXG39fD90uLd7/M8Ojfj6A42I5l1aYrqp/E8AWzvmXpY86/hwwxmYamjsYYxMAXA7d4fhrAG83NnMeuzgnbwfwKDcMsmmDc/4Jzvl8zvki6Pf1o5zzd6MLjh0AGGMTGWN94m8AbwCwCVHO+3Y7GXw6It4IYCt02+Qn2z2eGI7vbgD7AFSgP6VvgG5bfATAS8a/04xtGfSoopcBbASwut3jj+D4L4S+1HwOwLPGf2/shnMA4AwAzxjHvgnA/zLeXwLgSQDbAPwIQMF4v2i83mZ8vqTdxxDRebgEwH3ddOzGcW4w/nteyLYo5z1lshIEQXQoaTDREARBEAEgAU8QBNGhkIAnCILoUEjAEwRBdCgk4AmCIDoUEvAEQRAdCgl4giCIDoUEPEEQRIfy/wHGRz7weykzhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f54b0088e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotAT(rec.seq, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"green\"><b>Vaja</b></font>\n",
    "\n",
    "Napiši zanko, ki bo za dano zaporedje `s` izrisala frekvenco dinukleotidov CpG v oknih dolžine `k`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def plotCG(s, k):\n",
    "    s = str(s)\n",
    "    cns = []\n",
    "    for i in range(len(s)-k):\n",
    "        cns.append(s[i:i+k].count('CG')+s[i:i+k].count('GC'))\n",
    "    plt.plot(cns)\n",
    "    return cns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXmUHNV977+3q6u7ekbbaEVIDCMkgRDCbINsUMIOxjY2L4nj2HkhkDhPz7GNsZ/zbON3wvOSPY632PGxno2xHRu84WATAshis7FYtIEkhBYEGm1oQRpt09XrfX9U3erqmuqu6qrqrtvVv885c2a6u6br9u2qb/3qd38L45yDIAiC6H5ScQ+AIAiCiAYSdIIgiIRAgk4QBJEQSNAJgiASAgk6QRBEQiBBJwiCSAgk6ARBEAmBBJ0gCCIhkKATBEEkhHQndzZ9+nQ+NDTUyV0SBEF0PWvXrj3MOZ/htV1HBX1oaAhr1qzp5C4JgiC6HsbYLj/bkcuFIAgiIZCgEwRBJAQSdIIgiIRAgk4QBJEQSNAJgiASAgk6QRBEQiBBJwiCSAgdjUMnCKJ32bB7FE9uPYQ/uvQMnDZZi3s4odlzdAw/WbMHnHMMTuvHuy+ZG/eQSNAJgugM//LoVvx6+2GkFYYPXb0g7uGE5r7nduNrj++wHt/0ptnQVCXGEZHLhSCIDlEoVQEA+WIl5pFEw6liGROzadx102IAgF6K/3ORoBME0RHypuDlJRC+KNBLFWgZBbmMYZXL8LlI0AmC6AjCgpXBko0CvVRFTlWQM90sunkHEick6ARBdISkWej5YgWamoKmpqzHcUOCThBER0iahZ4vVZBTFWshVIYLFQk6QRAdQbgkZHBNRIFeqiBrE/QCCTpBEL2C5XKRwDURBbppoefIQicIopcoVaqoVDkAOYQvCoTLhaJcCILoKexilxQful6qGouiaYpyIQiih9CLyRP0fKmCXEaBlklZj+OGBJ0giLYjrNeMkpJC+KJAL1ag2ePQJVgbIEEnCKLtCBGf0qdK4ZqIAr1sCLpmJRaRoBME0QMIQR/oyyTCQi9VqihVOHKqAlVJIZ1iUnwuEnSCINqOsF4H+lUUy7WIl25FfB7hbsmpihR3Hp6Czhi7mzF2kDG2yfbcVMbYSsbYdvP3QHuHSRBEN2O30AGgUI7fmg2DEG+R9p9Vla6x0O8BcKPjuU8BWMU5XwhglfmYIAjClYLlQzcEvduTi4SFLvznuUxKCh+6Z4MLzvlTjLEhx9M3A7jK/Pu7AJ4A8MkIx5VoRt4Yw73Pj6DKjdvOGxafhkvOpJucXuVHz49g5+FTdc/1Z9JYfsVZsTdMaMS2Aydw/7q94DCO4XddcDrOO32y9fr+Y3n8+zO7UDZdKzsOnAQATO1XAcgR4hcGMX6RVGS4XOL/TEE7Fs3inO8HAM75fsbYzEYbMsaWA1gOAIODgwF3lyx+tGYE33jiFWTTKRQrVWw/cBJ333Zp3MMiYqBUqeKTP9sIJcWQTjEAQJVzlCocS+dNxVvOmhbzCN2557ev4YfPjiCbTqFQruLQiQK++J4Lrdcf2LAPX3/8FWTSKTDzuTlTchic2gdAjiScMFgWuplUpEnicml7CzrO+QoAKwBgeHi4u1dCImKsWMFELY2Nn3kr3vPN1RgrluMeEhETQhjufNsi/MXvngUAWD9yFL/3b7+V2i2RL1ZwxtQcfv2Ja3DDl54cN9Yx8/HWz98Ixpj1/MqXDgCQI8QvDOLzCgtdUxUpvq+gUS4HGGOzAcD8fTC6ISUfvVSxbqWNK3t3WytEcIRVl7W5VmSqDdKIfLFiRXi4WafGMZ6qE3MAUhWyCkPe4UPXVAV6Of7zOKig/wLArebftwJ4IJrh9Ab2kyGnpqTIMCPiQS8aIpCzC7oQPYmPC1GYCnC3Tu3HuB0RFdLtFrozykWW89hP2OK9AFYDOIcxtocx9n4A/wDgesbYdgDXm48Jn4iiPoC5mNLlIVxEcMR3bxc/K/NQ4uNC1AIHxDFcHfe624Ku1gUXKz+4xaHLcNfhJ8rlfQ1eujbisfQMXtYN0TuI715c4I2/5Rc9vVSxQhA1NQX9mMNCL7lb6MKdJIN7Igx6abwPXYa7DsoUjYH8OB96/AcCEQ9OS8/+d0Fi0bMLtttdpnEX2thCl8E9EYa8pFEuJOgxULAJei6joECLoj2LJQyZmvipCkOKyW6h29yGmfF3mWJR1EnSFkWtOPQMWeg9S53LJa2gWKmiXCFR70Wc8cwAwBiTxifbCFELHACy6fFjtb9uxyo1K/Fn84NYFM2mDQnV0gpKFR77eUyCHgN6qWq7spur/hLfXhPtQwiDU/xk8ck2QtQCB9zvMvVSpe4iJRACKPPFyg/OsExZzmMS9BjI225Hk2KxEMHIu/jQAXl8so0QtcABY+zOu8x8qVLnRhKkUgzZdPc3uXCGZcpyHpOgx4Ddusl2QUQD0T7colwAeXyybthrgQO22HKbdao3iEMHkrFupDuieGQ5j0nQY0Av10cIAPFf2Yl4ENEhzogQTU1JW+/ELQbb/jxgiLvboihg+JvjFr6w2CPVAHnOYxL0DiOsG23cySDnyUu0F71YAWM137IgJ3F+gjNL0i1uvlGmKGBGxXS5AeNMnJLlPCZB7zBO68Y6Gbr8ACeCoZer0NLKuJonmsQZxM5a4OK3aFrBOa+7C3WSTctROzwM9sAGQJ7zmAS9w1jWjSPKJe4DgYiHfNE9vE/mDGJnlmSt9oxxbBfKVXBeX3DMThIsdHtgAyDPeUyC3mFqccf1t6vdbrEQwWiYIi9x2KIzS9JZHdIt+9WOzJ/NL06XUjYtx3nc9nro7WDT3mM4eELHNYtmoViu4nurX8Otlw9BVaK/Pr1y6CR+8IzRXWj+jH7cctlQqPdzZpiRoPceK556BftGdQDAul1HkXVZPNTUFA6dKOAzv9iMxbMn4T2XntHpYTZkfC1wY/zf+vVOPLRxvy1yx13QNVXBtgMn8KWV23D7NQuQbsN5GyWb9h7DT9fuqXtu72ge82b0W4+tGjUk6K1z07/+BgDw2j+8A99b/Rr+5j+3AIDVICBKfrxmN+5++lWrM8sfXTqITDr4AegnQoBILm+cLODvHnoZmppCxhSyt58/e9x2Fw8O4JHNB/DD50aQYpBK0EV4ohDsM6f1Y+5ADs/sfMPaZvqELM6dPdH1/4eHBvDMzjfwlVXbccN5s+pa18nId55+Dfev34OJ2Xq5vGSw1jZSlpLHXSnodo7nSwCAE3p7uv7kixVM6VPx4asX4G/+cwv0ciWUoDutF1kOBKIziE4+n7t5Cd4z3Fik37t0EO9dOogvrdyGr6zajmqVI5ViDbfvJM7Y+ekTsvjNJ6/x/f8fvGoBzp8zGbd8+7muOO7zpTLmz5iAX/2vKxtuI8udttz3Oi3A2nSsiwSCqL4wt04nxvMUttgLFFzqnzdD3MrLVHnRy0fuh9r5JM/naoS9EFkjcpKcx10t6Jxz29/t2UfeLANaK/sZ7guzaneITNF0Mjq4EP4QkSCN/MtONAlrnzijXILQTVUXm8XUC2SpUdPVgl6q8PaZ5iZ5M00/qgOwFsNrTH0qxcysQPkPbCI8jWq3NELG/qLOKJcgyBK37QdnVqgbokZNgQQ9OHq5Ypnm7dL1QrmCnJqqVVOLyOXiTErohgObCE/NuvV36snim7Xjdgy3Sjf1Fm3UTs+JDPH13S3otgWVtrlcTAtdWCNRWejOSm3dcGAT4RHHT9andStjOzpnLfAgdFN0l7MQVyNkOI+7WtDzpQpKVUPJi20qLC8SP7SIbn2di6KAaDAr/+IQEZ5W/c+1dnTyCJ+zFngQLFeSRBeqRjRK/nKiSXAed7Wg66WqdUC068qom3WdrRMrtIU+3rrJSpzmTURLqxEiNeGT54LvZ5HQC83KrJTnczXCT5QLIEe5hq4TdHtkS75UsU6Q9gm6UTwpqkUccftmt25yakoqC4xoH15ZlE6icvVFiV+fcjNSKYZMlzS6aNSsw4kM53HXCbo9HjdftAt6G10umdS4AkSB3684vnmuW5NdIpnkHWGrXkS1GB8lfl0QXmhdUHWxUuUolqv+XS5kobdGfRH9inWFb9dECos6qkUctwUWzaXJLpFMxPHjd0FRxvC+KCx0QO6uTIJCgwYkbsjQ2LsLBb2+zZV43I6J5JxbMaiigFIUi6LO2zetCw5sIhr0UgXZdMp3Gr+MYYt+fcpeyCCAXliFyPxY6BKcx10n6PYDIF+qWejtmEhR11lTFWTTKTAWjYXuTMgwwp3kXxwiwmO48PxbtzKG97X6GRohg4vCC1GIzJegp+M/j0MJOmPsY4yxzYyxTYyxexljWlQDa0Rdm6s2L4oWbP5OxlgkbcGcnU4AI8lCdkuFiIZWI0S0iNZuoiSKKBegOxLqxPnuVuLYSS4T/3kcWNAZY3MAfATAMOd8CQAFwHujGlgj7G259FLVEvJ2TKRbIa2wbcHcFpRkSEggOoPRPNm/GCophoySkqodnV6uNOxG1AqamrKMJllpJczUsNC7VNBN0gByjLE0gD4A+8IPqTn27FC9zuUS/YHhTNM2LPQ2RLmYlgpvV7orIQ0i87gVNDUllWtCj8hC7wYfeiuJYCL1P87zOHA9dM75XsbYFwCMAMgDeJRz/mhkI3PwzSdfwa4jY9g3mree+69N+3H4RBFAOAv92FgJX31s+7j3UMxY8VpD5xTW7DqCT/98IzJKCh+4cj5+unY39h3TceXZM/DW807z3JdeHn9CZ1UFnBs++yiiBwh5eGjjfvxmx2GcMdAHvVTB5n3HMHtya57JrKrgnt++ho9ddzYm96ltGqk39zz9KrYdPIk3ThWjEfSMgpEjY/j0zzcCAPozCj52/dnoy0TTpuEXL+zDwpkTcO7sSdh24AS+t/o1qEoKH7lmIQb6M67/c+hEAV9/fAcqVY7blg25ZnY3QjPP4zvv32gtek/vz+Dac2fhR2t2445rF2LWpPZ6pQPPHGNsAMDNAOYBGAXwE8bYn3DO/92x3XIAywFgcHAw0L7GimX8/X+9jP6Mglwmjfkz+jFRU7HnaB4TtDSyaqrOcm+V1TsP49u/eRUDfSqUlGE9nyqUa3U3zC9z2YLpeGjj63hk0+t441QRU/sz+OLKbQCAjXuO+RN0F+umloVKgp40vrRyG7YfPGk9nphN47KLp7X0HvOm9+PQiQJW7zyMG5eM727UCTjn+NyDL0FTFUzKqbjkzAHvf/LgzfOm4blXj+LRzQdQqlRxLF/CVefMxLIF0yMYMfCRe9cDMDqb/WztHvz7MyMAjG5Q77zgdNf/eXLbIdzz29cAAFP6VJw/x+im5OcCduEZU3DaJA2/2nIQgJFVfqJQhpJK4YfPjmB5GzqqOQlzKbwOwKuc80MAwBi7H8DlAOoEnXO+AsAKABgeHg50LyJuNz/5tkX4U5eenv/08Mv4f7/eGeStjfc3hfv+Dy7DvOlGn8B/eXQr/vWxHQBqX+bnbl6Cz928BCf0Es7/zKM4cqo47j387Msp2vYSqZMRnwVGRI/zuPjb3z8f72ogJo34xz94E67+whOxuicK5SqqHPjQ1QvwoasXRPKet14+hFsvHwIAvLhnFO/62tNtcy05o+N8bVesjGu314xlC6bjmU9faz1+aON+fPAH63B0rOj7PcISxoc+AuAtjLE+ZuSxXwtgSzTDqsfrtkdTFZQqHKWABbqEX9x+FXYWz3LuD4D1RWXT/n2cjaJcALmSR4ho0Ev1LQuDuCqiylIOQxRdipoRdcOLSrXedswXK76ayYg7/axZlkC3SjW0LpU5h060a+7sBBZ0zvmzAH4KYB2AjeZ7rYhoXHWIBc9Ggh42VtfZdMK5L+d+VSWFdIrh6JjRz3SgL+OrhoM9USnK8RPyopeqGLD5vYMIgwy1w73OwbBEnUDlPB/1chUDfRnPfYjXpvSpRvG/EBcyEeoodELzWQM/DKH2wDn/v5zzRZzzJZzzWzjnhagGZsfLOhCZl0EjXRqVtHX72/7cUdPlMqVP9WWhF6zbt/pplzG9m4iGfKliCQkQTBhkOD5qTS3aI0pRC7rzfBTN3o2/G+tEvlSBqjBMyKah2/JcgiRSie96dKwIxoCMIrmgd4q8iwVtRwvZl1MvVcBYfX0N+4HrdmXVMop1KzW1P+MrXKnRhanWr5QEPUmUKlVUqrxO0INYuCJLOc72Zm13uUTcas/5Pnqpgv5s2rPCo7iDFnVmwrTbE5/piBkVFKZ+vF+6QtC9DqZcJrzLRUs7S9p6W+ijNpdLlZs9TpvQ6PbNcrlIlDxChEd831NtIXJBLD0rS1kCC71tLhfLKItmncD5PqIph1eFR100hTcL5uXNNRC/tXfsiIvA6FipI/5zoEsE3auGdNgFFbfaFPZMOLf9amoKJwtlAKjdynnsv5EfUsYmBkR4xB3XlDoferATO+40ed3jHAxLWklBVVhkn9Ep2iJD26vCo27bLl+qGKHEAVvtifP6ZKHcsXDk7hB0H1EuQPASuvni+HrH9seqi+/L/rqfxRb7+JyfQ8YmBkR4xAU8rA9d/F+chZ/E3WM7LU0twhIYznMpXzLKFXhdGEUmd9YstJUvBi9EZnfTRFGd0g9dIehWkawGE2v5oMvBDnijNsX4dPxm2EVZWGCegt5ggUWTsIkBER7xfduzEoMKYtwF3KzQ3giqLDYiSkF3vk+hVLX6GjS10Mv1lnyYZh72tbd2zpudrhB0r9ChXEgL3S170+sWyf668JF6nXBiUct5C0dhi8nEEnSby8VvYwsnmqrEumgeZnHQL1FUMxWMi3IxhdmrkbOotZMz6+eEaeaRUVIQrvd2zpudrhL0hlEuIeN0xVXZjtdVOedioXsdjA0tdBL0RKJbgl6z0IMsrgGmyyXGRXMrV6ONsdRRupXsxlW5UrVcKZpHmRAh4KKyqlveiF8YY9b/koVuQ/ewDsJGubhVwPM6cMU+s+kUcmra3H/zg7HRnYZIVCIferIQ32cUBbXi7jvb7rBFIFq3kr0sr16u1lwpHhdG3eaaERZ6mM9cK+xHgm7hFToUPsplfFEsb5eLWVI3o/i+oDTLtouiNC8hF4UIRTAbczecWjZ1e33o7YhDP5YvGZ3HzHO12YVRRLxpqoJCuYoxl3LXraCRoI/Hq/5y2Ey6gkvYot9FUS2t+K7F0ixaJxvzLTURPWHSxp3E3VA5X6ognWKuEV9RkcsokSVP2c9FkdGtpRXPhux5Ea9uZXmWQrlLLMOPolxqiNugRmRDJiXkS5VxC5VeB27O5hvzuygrfHduB0guE64EMCEfUUaG5CSIcml3coyX2LaC/Vy0imNlFM9GzrptURQw0vbDWNfiu6fEIhviqtkIkUkX2IfuYqF7Ib5kVaktfHhZ2LW1APe4dvKhJwuvtZ9WiDKkLwhRtZ1rhkjmiQL7uWgVx1JTno2cRQMacU6fCtBhyo747jvlcommNUib8bPSHKZNV5CFD7dSuz9+fjfW7Rpt+D+b9x2DqjCkXax/txP2a49tx6uHx7BswTT8/sVzredPFcr4p4dfxsmCsf25syfiLzpQPL9b+f7q17BswXScNWNCW/exYfexuue2HjgOIJrIkJyq4Lhexsd//ILr6+kUwwevno8zp/UH3sfIG2P4tyd2uJawWLPrSNsKcwm0BlEurxw6iW8++QoqVWDOQA4fu27huLoov91xGD9bt9d6/MKe2nn4nadfBQAzvjyFU0X3eeTgKFW4FYcuCGNdi/wWEnQbk3Mq+ECu6TZBLfRqlUMvVV2tj3decDqGG3RmuWhwCs6c1ocbl8zGxGwal8+fhl1vjOHwyTea7u/aRbNcn3cuCOmlCr7wqNENad3I0TpBf2H3KL67ehemT8iiUKrgPzbsJUFvQLlSxV8/sBkDfSrW33VD2/bzz49sRZUbx6qdK86egYySwm2XD2GSFvx0u3RoKh7atB/P7HQ/vvaO5nHWjH78zyvnB97HI5tfx33P78bpkzXXQlLXnjsz8Hv7oVFI4X++uB8/XrMHk3MqjuVLuO3yobr6OADwvdW7sOrlA5g5sdbibdmCaTh8ooiDxwtYMHMCFp02CVlVwRkDjefxzGl9uOTMAcyalMX8Gf0oVqpYOm9q4M905dkzMHJkLJIOT37oCkH/wh9e4LmNFvB2TZS0dbsK/+v7Lmr4f8NDU/Hk/77aevzD//GWlvdtJ6cqGM2XrMf2uw23uhQA8O1bh/HUtkP4l5XbUKpU27pg1a2I7GFx2922/ZSq+PPfmYdPvW2R6+ufedd5od7/usWzcN1id2OgUuWY/+mHQkfBiOPqqU9c7XoX2W4auR1FSds737YIn7p/o6vhli9VsHj2JDzw4d9puo+h6f24+hP+LkyrPn6Vr+2asfyK+Vh+RfCLbKskRgEM31jrgl6Lr413KjQ1VbfCLw5st3Kf9gSlqMuOJo1OxG6XK1UUK+1fNGyEkmLIKOEXTYVwxiHmgCHo5er4zmNW9maTYz1MAlCSSIygG2FdrVso7S4L6hendWKVXu3LjBMl8TlztsUbipBxpxMLieIuoN0+5mZoavOysH5wS7DrJI0ypgvlSv2x7vI53UKPe5HkCHrAKJFG6fidxpnwYG+FVShXUbX1SBRjztriZeNMOpGZjgi6BEZBFFEwhXLMgt6g85i40DQTdCP0mAQ9MYIe1EIR/5ON+WDIOlxGzjogBVslSSueXbXFwJPLxZVOzItXvf5OEEXIX94jga/dNCpSZ9Uyt/I9xhsvQUKPk0iCBD3YAR2mZ2CUOF1G4qAVq/luYq+ZYVgACXojOuFDL3SgVrgXYfIwBF4JfO2mUca10UUo1bQIn+5SvqMXSYyg5wKWF7Wy+STwoRfNHpRAvcsFwDj/ukjDFreZVKnRnaA18ltBHENxCkrWoyysH7wS+NpNMwtd87gb1UPWXEkKiZkBo9xl6wd0zdqNP8oFqI0n73C5OAXdquJGUS5N6YSFHmXNlqDkPMrC+iHuSJFGJTR0W8EsoHGUS9xGmQwkRtCDlheV4WS07z/vFHRXl0vVEnLLqqEoF1c6uSgaZ5RLFPXS444UyTYQbNHEXQi6s4BXqVJFucpjP4dlIDGCLnzonI9PW26GLGGLzr6oBctCH9/ezl6qwBJ0qtToSicE3Yo6inFhPYp66XFHitRcLo4oF/NC0ygOXYYoI1lIkKAbH6XQotulIMnBYFkf5QYuF9vKvm7zddYuBBS26EYnXFEyLKxHUakw7kiRRn0FxIKn1qCqqmWUUZRLcgQ9aF9OaeLQHcIsfrs1oM67WOjkQ3eno4Iecwx32FyEuCNFGka5mAueaSUFVRnf2atQkiOwQQYSJ+itnsBWhELA5r1R4byd1MsVZJQUJmTTdc8D9Rl9opIfRbm4Yxc5e3JWlEgRhx5F2GLMkSJeceiA6Vp1aQBtvJYYOQtMqBlgjE1hjP2UMfYyY2wLY+yyqAbWKkEzJvVyvPUrBOOiXIqVukzQukqM5ZollVFSYIwEvRF1aw9tWmfIS2AhRtGPM+5IEbdjXSx42lu5FRzfY74Y/x2SLIRVsa8AeJhzvgjABQC2hB9SMJyLin6Ju36FwHkwi4VPYbnbV/btLflEc484GwjLjH1e2jVHtWzjeK3bikthK7/IECmSTZvGiUsJDLuL0S2s0b5NLxO4fC5jbBKAKwDcBgCc8yKAYjTDap2gVQfDdvWOCuftZrPYW+fiVRQha0lFd9zZtGsf2SZNzDuB/TgJUkZZhkgRxphRNdX2PTkXPN1qNtVqG8V/HsdNmHroZwE4BOA7jLELAKwFcAfn/FQkI2sR4QP/wiNbMW2CERmSTSv4xI3nYNYkreH/6ZKU3RRjuPvp17DypQNYu+soJmmq9bnyxSqOjZXwtw+9hIMn9Dp/oaYqeHLbIXz4h+t87WtKn4q7bjoPmZjXDYKyae8xrHhqJ6o+QlTXj9Q619z1H5vw0evOxvlzJ7e0v4c37cfDm17H7dcuxHyXrkdxR4cANpdjsYJJmuqx9XhkiRTJZRQ8++oRfOTe9ahyXlufSIuorhTWj4ziwz9ch+sXz8LNF84hC91GGEFPA7gYwO2c82cZY18B8CkAf23fiDG2HMByABgcHAyxu+YsmDkBF8ydjAMndBw4oaNc4Rg5MobL50/DH1wyt+H/xe03FMyYmMXl86fh9eM6XtpvdBq/fvEsa2VfL1ewduQIfrxmD+YO5LBswXTrf9963ml4YttBvLT/uOd+xgoVvH5cx/uWDuK801sTNln45Qv78MsX92HedO92a1k1hRsWz8Keo3k8tvUg5k3vb1nQv/nUTqwfGcU5p03CX141XtBF4kucNIrh9ouIFIk7OEBLp/DC7lG8uGfU+n7PnT0JF5wxBQBw3bmz8PMNe/HYywex+8iYKejRNePudsII+h4Aezjnz5qPfwpD0OvgnK8AsAIAhoeH2xNmAGDahGxdt5KDx3Us/btVnq6IvC3rMk5UJdWw65FY2RcROXffdinOnjXRev2udy7GXVjsaz9PbjuEW+9+rqsXUfOlCqbkVDzWYkeZSz6/MtDCobO+zvjxVGMXk7CNTmQJ3xXn4rT+rOv3e/u1C3H7tQvxge+vxc7DJwFQlIudwDPAOX8dwG7G2DnmU9cCeCmSUUVA1uciqWFdyX0giJC0KLrIN0rO6CaCLmQ3akLshbBeGwq6BAvrjWK4/SJLpEgtPLH5OWmUy67/XuIeuwyE7Sl6O4AfMMYyAHYC+LPwQ4oGv4lGeqkyruGsbIjmBTU/Z/ALkGXJdXFUTFA3WdCa+c76Ok6MxhBxh72Gq7opw6Koff9e36+9/rss5TtkIJSgc843ABiOaCyRoioMKeZticoS5dIMsbIfxUmXhNoveqkaKKLBqDnf+ue2Rx65EXdjCCB8xrAsomiFJ3q4fuwdmvRiBYzFGzYqC4mdASs+2+MAj7tkqB+0jFHrOopby6Dx+jJhXIRbP3SD1jupWYLuxoEMC+the8vK4rawEog83Ip1gl6uQksrYCy+sFFZSKygA/7acuWL8nc60dIpy+UiGlsEfq+Qt+YyEDRMMGibNnv2bqPX4z6Gwt55iTvZ+F2UU7uNAAATlUlEQVRHZniix/ebUxWUKhzlStW4Q5IgsEEGEi3ozj6dbhQksK68EK6CKC4+tYp23bsoGjRMMMiiaLlSRaliRLk4U85r44nfKKitjQT7XmWJcqllhDaXptoFrGqW/U20lPkm0bPgx2cad9stP4h0Zz2CruxWolKXW+hBQk3tt+l+qctabGChG3cMccdvh/ShyxLlkqnVbGmGFdVTNNaWZAg9lgG5lSwkOQ+LTIb6FX7QzNR+vRheOBqVIO0mCgGbGefUVMtrB36Ke8mQWBS26qb4bHHfafiNcrG7DrshsKFTJFrQNY8TWIbGBH4wEovErWX4sQaxVGUi6F1VkJo3XsW9OOdSpP5nlJQZ1RVQ0IvxFxgDbIuiPsIWAVhrS3FfiGQh4YLefBGsW4r62BOLohCOrhf0gGGCbrW0vbBXUnS72yuUq+A8fsuWMRbo8wlEpE7ckSK1KqLNt7O7mPSAd2xJJPGC3ky49GL8daz9IBJiorLQu7ncLuc88FqCUUu72lKjCyHiA30Z12PJqoEiwTEUpuqmsbAbvxz4LRhnT5CTIVNXFuL/BtuIVxeXmt9Q7mnIqQrKVY6ThXIkiz9eawsyE8YitmrLt1BGV9zFTelTXe/2rOgQCQRFuOaCIEMsfStYPvRyta7Hbq+T6FnwSiySZWXfCyFER0+VAiXUOImiu01chEmACRLhI7ad2p9BvlQBd5Tsra3DxH8q5TJhLPTuihRxRrnIfg53iviPwjZiL+DjhkzWVTOEj//oWDGSW0uvtQWZ0UO4OIJUJBSCPdCXAedA0dERyEqZjznKBTCP9xCZojJ8Br/YazXJsCgtC8kWdI/MQGvBS3JBFwfvWEQ1Q3IZpa6lXTeRD2ERB8mS1W0uF6C27uIcjwzWrZ9SF43oNlHUHIJOPnSDsNUWpSanKiiWq1j+vTXWczMnZfG5dy1BKsWkqV/hhX18kVjoAWuaxMXqV97APb99FZwDJwtlAAFdLj7r2Pzo+RGs2nIQALB3NA8AVkXOj9y33mg3xxjOnzsZT249FHg8UaOpCjbuPVZ3vAPABC2Nz9+8BC/tP45v/Xon3Bo9bd53HEu6qOGJ3ciRIVNXFhIt6JedNQ2/mnMAI0fGAADH8iXsf0nHB69agNOn5KRJd/bi/DmTcdHgFBRKVVxx9nTvf/DAyKDtnkXR/1i/F6u2HMSCmUa3oIsGpwTqtuS3pPJ3nn4Ne4/mMWcgBwC48uwZeOt5p+GpbYdw4LgOANh64AQe3vw6MkoKlw4N1DUciYsbFs/CoRMF63gHDMEbOTKGPxo+A798cR9+teUgFs4c33XptEkabjhvVieH68q7L56LJ7YexAeunN90O3HOHsuXjMck6AASLuhvPmsaHrz9d63HD2zYizvu22ArhypHQSIvBqf14ecfXBbZ+3WbDz1fqmDuQA4Pf/SKUO9Tu033Lql89aKZ+Or7Lqp73t4R602feQTH9TLOmtGPn3zg8lDjiopbLhvCLZcN1T23dtcR/ME3Vhs1T4pVnDZJCz2P7WRyn4rvv//NntuJBKjRMaMvvezncKfoqVmwd0YHuifKJWrCLJ7FQVTVDP3WDPeTqCIsRNnv7pwp8kkRPiORKoUjY2Sh20nGt+sTp4UmS1H/ThMmASUOolr0Egup/mrke7VA81e3O27GC7rc422FnKpYFrrsF9ZO0VOC7vShFkq92elE1JIuVbrDjx5VnHE27c+H7qeao9/OOnGTsy0Ed1vykBeaquCoKehZyS+snaKnlMyejADASqWPu35Fp+m2Jhd6qRqJcNoLOjWiWuUolr1dLn6rAsaNM7xP9gtQK+RUBUdPmS6XBH2uMPSUoDu7uiTtAPeLlql3PclOVDXr7dZqI/yWkRXjyUruk66tG1SNXqwJsmTtFrrsF9ZOIffRGDHOOGS9VO3JTid+w/dkIariS36iXPwulOe6xELP2sodRFWtUxY0NYWxYnfUY+oUPTULTldD0M433Y7leuoSQS+Uo/H9KimGjNK8jo3oUOQ7ykVyQU+lGLLpFAqlill2ODmnvP3iJPv30CmS8+36wNlPU48olb7b6FULHaiVIm62L8DblSKiW7ohakQ0x46ihaFMRJ1BnQR6StCd1faSdoD7xY8vWRaM+ufRNTDwrJHvsxyEkjIW0rvBhSHKSAdtDCIrWRL0cfSUoDv7aSbtAPdL1look1/QSxWOSpVHJpw5nwXb/O6vG4REUxWcKlZQKCer5on93O2GC2sn6ClBB+ottF4t6lNzucgf5WK1CYxo8dqr6UmryWbdsBinqQqOmRmVSTre61wuPRjc4EbPzUK9oCcnFboV/MRjy0KhRYvZi6yqIN/kQiYucn7v3NSU/MePpqZs4X3yj9cv4txVFYa0kpzPFYbQs8AYUxhj6xljD0YxoHZj76eZtMw5v3RTlEvUTUhyHnVsklgOwkjASV6KvDgmZC+/0EmiuKzdAWBLBO/TEez9NJMWl+uXbopyiV7Qm9exEWKfpOMipyo4mkCXiwg5TqV6K9O7GaEEnTE2F8A7AHwrmuG0H3s/zV7tdOKsOikzYVrOuaHZ7tBc9ycyRRPkk7WXS07S8S4sc2ef114mbD30LwP4BID4q/v7RFMVbNg9ilvvfq5nF0Wz6RQYk3tR9OBxHXc9sBn7zYYSUX1POVXB7qNjuPXu51xf333UaA6RJAvdPndJcjGK76jXajE1I7CgM8ZuAnCQc76WMXZVk+2WA1gOAIODg0F3Fxk3XXA69PIejOZLuHhwCn53YfgOQN0GYwxaunm0R9ysGzmKhze/jnNmTcSyBdNwzmnR2Aw3nDcLOw+fwqjZ6cbJRE3FzRee7il8f/XWczBWrOAdb5odybjayfWLZ+GVQyfRl1Fw7uxJcQ8nMi4dGsDSoam48pwZcQ9FGljQ2xXG2N8DuAVAGYAGYBKA+znnf9Lof4aHh/maNWsavUx0kIs/vxLvOH82Pv/flsQ9FFd+vn4PPvajF/D4X12FedP74x4OQcQKY2wt53zYa7vAjkLO+Z2c87mc8yEA7wXwWDMxJ+RCSzdPgY+bfLG18EGCIHowDp0w0DwyJuPGbwo+QRA1ImkSzTl/AsATUbwX0Rm8MibjxsoQTVAiDEG0GzpbehTNFo8vI3qPtgckiDDQ2dKj5FT5XS45tffaAxJEGEjQexSvBJu46dWkL4IIAwl6j6KpqaYp8HGjl6KrgU4QvQIJeo+SU5WmRariJqrG0ATRS9AZ06N4NXqIGz3CtnME0SuQoPco0ke5RNQYmiB6CRL0HkVU4JO1Ul2+2JuljQkiDCToPYqwfgtlOa30fKmKLDUuIIiWIEHvUcSCo6zZooUebT5CEGEgQe9RcpI3uTDaA9LhSRCtQGdMjyIiSGRdGNUpsYggWiaS4lxE9yHE8o771qOviWtj3vQJ+LvfWxJJCv6Kp17BYy8f9LXtsXyJolwIokXIQu9RLjxjCq48ewY0VUGVw/Vn72ge9z43EtnC6X3P7cb2Aycb7s/+s3TeVFy9aGYk+yWIXoEs9B7ltMkavvvnS5tu852nX8Vnf/lSZO4PvVTBNYtm4p//8ILQ70UQxHjIQicaokW8cEoFtwiivZCgEw2xImEiqvmSp1BEgmgrJOhEQ6KMhOGcQy9VoVHDCoJoG3R2EQ0RyUdRuFzEwqpGFjpBtA0SdKIhOctCDy/owm1DoYgE0T5I0ImGCH93FIIummnQoihBtA8SdKIhUUa5kIVOEO2HBJ1oSJRRLuKiQBY6QbQPEnSiIVqEPnQRKUNt5QiifdDZRTSkVmI3fNiiuCiQy4Ug2gcJOtGQtvjQKWyRINoGCTrREFVJQVVYJIJOUS4E0X4CCzpj7AzG2OOMsS2Msc2MsTuiHBghB1paoTh0gugSwlRbLAP4OOd8HWNsIoC1jLGVnPOXIhobIQFaJhpB1ynKhSDaTmALnXO+n3O+zvz7BIAtAOZENTBCDnKqEknYIkW5EET7iaQeOmNsCMBFAJ6N4v0IedDUFJ7Ydgi//29PW8+984LT8WfL5rX0PhSHThDtJ7S5xBibAOBnAD7KOT/u8vpyxtgaxtiaQ4cOhd0d0WH+eOkgzp8zGf3ZNPqzabx6+BQe2LCv5ffRSxWoCoOqkIVOEO0ilIXOGFNhiPkPOOf3u23DOV8BYAUADA8P8zD7IzrPbcvm4TabNb78e2swcmSs5ffJlyrQ0mSdE0Q7CRPlwgB8G8AWzvkXoxsSITO5gIukeqlCpXMJos2Euf9dBuAWANcwxjaYP2+PaFyEpGhpJVBcul6qUsgiQbSZwC4XzvlvALAIx0J0AblMsKiXfLFCES4E0WboDCNaQlMV6OXWa7vkSxWy0AmizZCgEy2hqSkUy1VUqq2tb+ulCoUsEkSbIUEnWiJoWzoSdIJoPyToREsEbUtHLheCaD8k6ERLiFjyViNd9FKVSucSRJshQSdaQgthoVOUC0G0FzrDiJao9RltLdJFL5IPnSDaDQk60RJWW7pyiy6XMgk6QbQbEnSiJWoWun9BL1WqKFU4LYoSRJshQSdaIkifUWoQTRCdgQSdaAktQBw6NbcgiM5AZxjREkHi0Kn9HEF0BhJ0oiWE2+TLv9qOrz223XP7J7YexPu/+zwAEnSCaDck6ERLDPSpuPWyM1Gucjz44n7P7Z/adhg7D53CO940G0vnTe3ACAmid4mkpyjROzDG8Nmbl+DoWAkv7hn13D5fqmBKXwZf/+OLOzA6guhtyEInApFT/TW60EsV5DJ0mBFEJ6AzjQiEpqas6JVm6NRLlCA6Bgk6EQgt489Cz5cqVJSLIDoECToRiJyqoFiuourR6CJPNVwIomOQoBOBsBKMPGq66GVqDk0QnYIEnQiE35ouOjWHJoiOQWcaEYha1cXmC6PUqYggOgcJOhEIza+FTouiBNExSNCJQPhtFp0vVZClsEWC6Agk6EQg/FZdJAudIDoHCToRCCHSzWLRy9TYgiA6Cgk6EYiay6XxoqhYMKUoF4LoDKHONMbYjYyxrYyxHYyxT0U1KEJ+hEg3s9DFgilZ6ATRGQILOmNMAfB1AG8DsBjA+xhji6MaGCE3lg+9SZQLNbYgiM4SxkJfCmAH53wn57wI4D4AN0czLEJ2/GSKkqATRGcJUw99DoDdtsd7ALw53HCIbkG4Ub66age+v3qX6zYFy4dOgk4QnSCMoDOX58ZVamKMLQewHAAGBwdD7I6Qib6Mgr+8aj52vXGq6XaXDk3F8JkDHRoVQfQ2YQR9D4AzbI/nAtjn3IhzvgLACgAYHh5uXpqP6BoYY/jkjYviHgZBEDbC+NCfB7CQMTaPMZYB8F4Av4hmWARBEESrBLbQOedlxtiHATwCQAFwN+d8c2QjIwiCIFoiVJNozvlDAB6KaCwEQRBECCiFjyAIIiGQoBMEQSQEEnSCIIiEQIJOEASREEjQCYIgEgLjvHO5PoyxQwDc88S9mQ7gcITD6XZoPuqh+ahBc1FPEubjTM75DK+NOiroYWCMreGcD8c9Dlmg+aiH5qMGzUU9vTQf5HIhCIJICCToBEEQCaGbBH1F3AOQDJqPemg+atBc1NMz89E1PnSCIAiiOd1koRMEQRBN6ApB78Vm1IyxuxljBxljm2zPTWWMrWSMbTd/D5jPM8bYV835eZExdnF8I48extgZjLHHGWNbGGObGWN3mM/36nxojLHnGGMvmPPxWfP5eYyxZ835+JFZ1hqMsaz5eIf5+lCc428HjDGFMbaeMfag+bgn50J6Qe/hZtT3ALjR8dynAKzinC8EsMp8DBhzs9D8WQ7gGx0aY6coA/g45/xcAG8B8CHzGOjV+SgAuIZzfgGACwHcyBh7C4B/BPAlcz6OAni/uf37ARzlnC8A8CVzu6RxB4Attse9ORecc6l/AFwG4BHb4zsB3Bn3uDr02YcAbLI93gpgtvn3bABbzb+/CeB9btsl8QfAAwCup/ngANAHYB2Mfr6HAaTN563zBkbPgsvMv9PmdizusUc4B3NhXNCvAfAgjPaYPTkX0lvocG9GPSemscTNLM75fgAwf880n++ZOTJvkS8C8Cx6eD5MF8MGAAcBrATwCoBRznnZ3MT+ma35MF8/BmBaZ0fcVr4M4BMAqubjaejRuegGQffVjLrH6Yk5YoxNAPAzAB/lnB9vtqnLc4maD855hXN+IQzrdCmAc902M38ndj4YYzcBOMg5X2t/2mXTxM8F0B2C7qsZdY9wgDE2GwDM3wfN5xM/R4wxFYaY/4Bzfr/5dM/Oh4BzPgrgCRhrC1MYY6ILmf0zW/Nhvj4ZwJHOjrRtLAPwLsbYawDug+F2+TJ6cy66QtCpGXWNXwC41fz7Vhi+ZPH8n5rRHW8BcEy4IpIAY4wB+DaALZzzL9pe6tX5mMEYm2L+nQNwHYwFwccBvNvczDkfYp7eDeAxbjqRux3O+Z2c87mc8yEY2vAY5/y/owfnAoD8i6LmXL8dwDYYfsL/E/d4OvSZ7wWwH0AJhlXxfhi+vlUAtpu/p5rbMhiRQK8A2AhgOO7xRzwXvwPjtvhFABvMn7f38Hy8CcB6cz42AbjLfP4sAM8B2AHgJwCy5vOa+XiH+fpZcX+GNs3LVQAe7OW5oExRgiCIhNANLheCIAjCByToBEEQCYEEnSAIIiGQoBMEQSQEEnSCIIiEQIJOEASREEjQCYIgEgIJOkEQREL4/7diMYBERTrbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5494c905c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotCG(rec.seq, 50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Branje datotek iz spleta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seveda Python omogoča branje datotek iz spleta. Preberimo datoteko iz zbirke [ENA](http://www.ebi.ac.uk/ena) na naslovu: [ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz](ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz)\n",
    "\n",
    "Ker bomo tokrat brali datoteko FASTQ, moramo to eksplicitno navesti funkciji `SeqIO.parse`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SRR020192.1\n",
      "Seq('GATGACGGTGTCTACATTGTTCCCGACCACTCATCTCCTCTGTCATGCCCGAAA...CGT', SingleLetterAlphabet())\n",
      "[24, 23, 27, 30, 30, 30, 23, 23, 24, 23, 23, 30, 28, 27, 25, 25, 27, 27, 27, 22, 22, 24, 18, 18, 18, 30, 19, 19, 23, 23, 30, 30, 32, 32, 32, 30, 24, 23, 23, 27, 30, 32, 30, 32, 29, 28, 28, 17, 17, 17, 17, 24, 17, 17, 13, 15, 17, 25, 25, 24, 24, 23, 27, 27, 15, 15, 15, 15, 15, 17, 17, 11, 15, 15]\n",
      "74\n",
      "SRR020192.2\n",
      "Seq('GATGACGGTGTCTACATCGTTCCACCACTCATCTCCTCTGTCATGCCCGAAAGT...CCC', SingleLetterAlphabet())\n",
      "[27, 27, 27, 30, 30, 30, 23, 23, 24, 27, 27, 30, 28, 27, 27, 27, 27, 30, 30, 27, 27, 27, 27, 30, 23, 23, 23, 23, 30, 30, 32, 32, 30, 30, 27, 27, 27, 27, 30, 29, 28, 29, 29, 29, 29, 17, 15, 15, 15, 15, 15, 15, 15, 17, 26, 15, 15, 15, 15, 27, 27, 15, 15, 15, 15, 15]\n",
      "66\n",
      "SRR020192.3\n",
      "Seq('GACGACGGTGTCTACATCGTTCCACCACTCATCTCCTCTGTCATGCCCAAAGTC...CGT', SingleLetterAlphabet())\n",
      "[32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 37, 37, 37, 32, 37, 37, 34, 34, 30, 30, 30, 30, 30, 30, 31, 32, 33, 30, 30, 30, 28, 30, 29, 24, 24, 18, 18, 19, 17, 22, 24, 24, 21, 21, 23, 27, 23, 24, 22, 22, 15, 23, 23, 21, 21, 24, 24, 25, 15, 15, 15, 11, 11, 15, 15, 23, 24, 24, 23, 11, 11, 15, 15, 15, 15, 15, 15, 21, 15, 11, 11, 13, 13, 13]\n",
      "111\n"
     ]
    }
   ],
   "source": [
    "import urllib\n",
    "import gzip\n",
    "\n",
    "f = gzip.open(urllib.request.urlopen('ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz'), 'rt')\n",
    "\n",
    "for rec in list(SeqIO.parse(f, 'fastq'))[:3]:\n",
    "    print(rec.id)\n",
    "    print(repr(rec.seq))\n",
    "    print(rec.letter_annotations[\"phred_quality\"])\n",
    "    print(len(rec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "V zgornjem primeru opazite, da se poleg zaporedja zdaj izpisuje podatek o kvaliteti sekvenciranja posameznega nukleotida, t.i. ocena [PHRED](https://en.wikipedia.org/wiki/Phred_quality_score), ki je sestavni del formata [FASTQ](https://en.wikipedia.org/wiki/FASTQ_format)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"green\"><b>Vaja</b></font>\n",
    "\n",
    "Napiši funkcijo `trim_low_quality_tail(rec, phred_th)`, ki na koncu zaporedja `rec` zbriše vse poziciji z nizko kvaliteto (pod pragom `phred_th`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trim_low_quality_tail(rec, phred_th):\n",
    "    quals = rec.letter_annotations[\"phred_quality\"]\n",
    "    while quals and quals[-1] < phred_th:\n",
    "        quals = quals[:-1]\n",
    "    rec = rec[:len(quals)]\n",
    "    return rec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "111"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(rec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "trimmed = trim_low_quality_tail(rec, 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "85"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(trimmed)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}