Fig4c_biasing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_FFC(base_run, pcr_run):\n",
    "    a = base_run/np.sum(base_run)\n",
    "    b = pcr_run/np.sum(pcr_run)\n",
    "    old_settings = np.seterr(divide='ignore',invalid='ignore')\n",
    "    ampratio = b/a\n",
    "    return ampratio\n",
    "\n",
    "def func(x,a,b):\n",
    "    return a*x+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "ready-to-sequence pool\n",
    "run36: biased chip; no PCR\n",
    "run42: biased chip; PCR\n",
    "\"\"\"\n",
    "counts36 = np.load(\"./data/run36.npy\")\n",
    "counts42 = np.load(\"./data/run42.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seqs = []\n",
    "with open(\"./seqs/P9-20161114-150-62-1.6M-ID37-AdapterSequencingPrimer-IndexAdapter.txt\") as fr:\n",
    "    for line in fr:\n",
    "        seqs.append(line.rstrip(\"\\n\").split(\"-\")[1])\n",
    "        \n",
    "GCs = []\n",
    "for seq in seqs:\n",
    "    GCCounts = seq.count('G') + seq.count('C')\n",
    "    GCs.append(float(GCCounts)/len(seq))\n",
    "GC_array = np.array(GCs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope: 5.233144119255713e-05\n",
      "p_value: 0.7619775547678641\n",
      "r-squared: 5.972257451638687e-08\n"
     ]
    }
   ],
   "source": [
    "ampratio = np.log2(compute_FFC(counts36, counts42))\n",
    "ffc = ampratio[np.isfinite(ampratio)]\n",
    "gc = GC_array[np.isfinite(ampratio)]\n",
    "xdata, ydata = 100*gc, ffc\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(xdata, ydata)\n",
    "y_fit = func(xdata, slope, intercept)\n",
    "print \"slope:\", slope\n",
    "print \"p_value:\", p_value\n",
    "print \"r-squared:\", r_value**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAEoCAYAAAAg8EqiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYHFXVuN+TySQxMGxhDUsgIIQYDHsgIAkg4sfyKYiI\nyCKIH4rrh4orsiggCvqhICJIXMDlJ7JGRLaEBLJAAIEQAmELJMO+DltIZs7vj1vNVFe6btftqZ7u\n6Trv89Rz+3adunW7pqZO3XvPIqqKYRiGYRSRQY3ugGEYhmE0ClOChmEYRmExJWgYhmEUFlOChmEY\nRmExJWgYhmEUFlOChmEYRmEplBIUkWEicpKI3Ckir4jImyKySET+KiK71tDe4SIyU0ReFZEuEblL\nRE4QEalH/w3DMIx8GdzoDvQXIrIZcCMwGngGuBVYAYwCDgLuB2YHtHcB8EXgbeAWYDmwN3A+sBdw\nSI7dNwzDMOqAFMFZXkSG45TcpsB3gHM19sNFZE1ghKo+mrG9TwB/BzqBPVT18ej7dYDpwBjg66r6\nqxx/hmEYhpEzRVGCZwHfBn6pql/Pob15wHbAUap6eWLfHjhF+IyqbtjXcxmGYRj1o+WVoIi040Zs\nawGbqepTfWxvQ+BpYBmwhqouqyDzNDAS2E1V5/TlfIZhGEb9KMKa4A7ACGCJqj4lItsDHwfWBZ4D\nblTVOwLa2y4qH6ykACPuAj4WyZoSNAzDaFKKoAS3icqlIvIz4BtAfPh7sohcDXxGVd/O0N5mUbnY\nI1MabW7mkTEMwzAaTBFcJNaKyu2BE4GfA1sAa+JGa0ui8tcZ21s1Kt/0yLwBCNAR2lnDMAyj/yiC\nEiz9xsHAn1T1m6r6hKq+rqpTce4RAEdGbhSGYRhGQSiCEuyKfb4kuVNV7wbuxo3cJmVo742oXMUj\nsypuyrXLI2MYhmE0mCKsCT6R8jkpswOwfob2nozKUR6ZjROyKyEirW2WaxiGUQdUNdeIXEUYCd4b\n+zwiRWbtqHwjZX+l9j4gIkNTZHaqcO6VUFXbcthOOeWUhvehlTa7nnY9m3WrBy2vBFW1E5gbVfdO\n7heRNXBGMwDzMrS3BLgHGAJ8skJ7k4CNcM7ymcOwGYZhGP1PyyvBiDNwa37fE5EdSl9GI7nfAKsD\n8zTm2C4iZ4rIQyJyRoX2zoraO1tENo8dsy7OylSBn9TllxiGYRi5UYQ1QVR1qoicg/MRnCUic4CX\ngJ1xkV2eBg5PHLYBsFVUJtv7h4j8GhdA+wERuZneANodwFXABXX6OUaCyZMnN7oLLYVdz3yx69nc\ntHzYtDgi8nHgy7hILsNxTu3XAGer6ksJ2SnAUcDvVfVzKe0dBnwJ55DfBiwEfqeqv8nQFy3StTcM\nw+grIoLmbBhTKCXYTJgSNAzDCKMeSrAoa4KGYRiGsRKmBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQ\nMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzC\nYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3D\nMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCym\nBA3DMIzCUgglKCJfEZG/icgCEXlRRN4VkedF5CYR+Uwf2h0mIieJyJ0i8oqIvCkii0TkryKya56/\nwTAMw8gfUdVG96HuiMjTwDrAfGAp8CYwCpiAexG4WlUPDmxzM+BGYDTwDDAXWBG1ux1wmqqe6Tle\ni3DtDcMw8kJEUFXJtc0iPIhFZCJwr6q+nfh+a+BWYF3gWFX9Q8b2hgP3A5sC3wHOjWs0EVkTGKGq\nj3raMCVoGIYRQD2UYCGmQ1V1VlIBRt8/BFwACLBPQJMn40aA56vqOUltpqqv+BSgYRiG0RwUQglW\nYUVULssiLCLtwHGAAj+vV6cMwzCM+jO40R1oJNG63hdxCu2ajIftAIwAlqjqUyKyPfBx3JTqc8CN\nqnpHPfprGIZh5EuhlKCIfBaYBLQDGwETcVOhZ6jqtRmb2SYql4rIz4Bv4JRoiZNF5GrgM5WmYA3D\nMIzmoWjTobsBRwGfBj6EU14nA2cEtLFWVG4PnIibEt0CWBP4GLAkKn+dT5cNwzCMelEoJaiqn1fV\nNmA48AHg/4DTgDkisn7GZkrXbDDwJ1X9pqo+oaqvq+pU4KBo/5HRdKth5E5XF8ye7UrDMGqnUNOh\nJVR1GbAQ+LaIPAecA5wPHJLh8Phj55IKbd8tInfj1g4nAU+kNXTqqae+93ny5MlMnjw5w+mNotPV\nBR/6EDz4IHzgAzBzJnR0NLpXhpE/06dPZ/r06XU9RyH8BH2IyFrAi8ByYLiqdleR3x+4DjeVuomq\nLq0g8zecQv2+qv4kpR3zEzRqYvZs2GMPWLEC2tthxgzYZZdG98ow6o/5CdaHV3BuEoPpXe/zcW/s\n84gUmbWj8o0+9MswKjJunBsBtrfD2LHus2EYtWFK0E1ZDgZexY0IvahqJy5EGsDeyf0isgbOaAZg\nXk59NIz36OhwU6AzZthUqGH0lZZXgiKym4jsLyJtlfYBv8NNbV6SCH12pog8JCKVLEfPwLlWfE9E\ndogdMxT4DbA6ME9V5+T8cwwDcIpvl11MARpGXymCYcwWwBTgVRG5B3gW6AA2B8biFOBU4IeJ4zYA\ntorKMlR1qoicg/MRnCUic4CXgJ2BkcDTwOF1+TWGYRhGbhRBCd4GnI7zC9wC2BU3insW+Dtwmape\nl3KsUu4I37tD9SQRmQV8GZc1YjjwFM7S9GxVfSnPH2EYhmHkT+GtQxuFWYcahmGEYdahhmEYhpEj\npgQNwygMFmnHSGJK0DCMQlCKtLPHHq40RWiAKUHDMArC/Pku1NyKFbBggftsGKYEDcMoBBZpx6iE\nWYc2CLMONYz+p6urN/C4BRoYeNTDOtSUYIMwJWgYhhGGuUgYhmEYRo6YEjSMFqezE377W1cahlFO\nEcKmGUZh6eyEzTeHd96BYcPgscdg5MhG98owmgcbCRotiTlFO6ZOdQoQXHn99Y3ph/09jGbFlKDR\ncphTdC8HHOBGgODK/fbzy4coq6yyXV0wcaL7e0ycWOy/h9F8mBI0Wg5ziu5l5Eg3BXrxxdWnQkNe\nHkJk5851f5MVK1x55521/x7DyBtTgkbLYU7R5YwcCccdV30tMOTlIUT2rbf8dcNoJKYEjZajowNm\nzoQZM1xZdKforNOWIS8PIbLDh/vrzYqtYxYDc5ZvEOYsb/QHpWnLUpSUai8FIRFVssqW1gQXLoQx\nY2DWrOZ/MQm9bkb/YM7yhtFg6jU6qFe7oeujHR2wyy7ZH/hZ3uM6OpzimzkzmwJshhGYrSsXB1OC\nhpGRelmd1tOatV7ro7X0OYvCbBbLXltXLg6mBA0jI/UaHdRz1FGv9dH583stPh980N/nri6YMMEp\ntQkT/IqtlmtRj5GjrSsXB1OChpGReo0O6j3qCJnizKpQRo2CwVG8qbY22GSTdNlp0+Chh6C725XT\np6fLhl6Leo4cQ6eGjYGJhU0zjIyURgd5p+KpV7uhhBiDLFgAy5a5z8uWOeWW5oLx5JP+epzQa1Fp\n5LjLLv5jDCOOjQSNlqRexhX1Gh3Uc9SR9VrUa1r2kENg6FD3eehQ+MQn/PIh1yJkRAr1iYhjDGxM\nCRotR7MYVzQDXV1Ooey+uyt91yJkKnLs2HLFtvXW6bIjRzqr0COOcGWeAbwXL4bly93nFSvgqafS\nZesVEccY2NSkBEVkhIh8QkTOEJFLROSKqDwj+n5E3h01jKw0i3l7PVMYZR2lTJvmrkFPjyt963Eh\nxiDz5pVPh959d7psZyfsthtcdpkr87we48a5rb3dKW2f4g4x5mmWe8ioP5nXBEVkMPBJ4ARgV0Ci\nLYkCKiKzgF8DV6jqihz6ahiZKI1oFizI39Ckq8s9IMeN8yuJ0BRGWdstyWZdu1u4sLz+8MNw4IH+\n9rO4MiTX9RYvTpetlMniuOPS5UOuRcga4qhRTll2d7spVN/UaT3vIaPJUNWqG3Ak8DTQDfQAzwBX\nAmcAJwLHReUZwFXAs5FcN/AUcESW8xRpc5feqBevv646e7Yr82xz3DjVwYNd6Wv7ootUnTpx28UX\n+9sdP961O3589T7PmqXa1ubaHTzY/c40rrmmvB/XXuvvxxZbOLkttvD3Y+lS1fZ2J9ve7uppPPyw\nqoiTFXF1Xx/GjXO/r9o1DmXWLHe9Sn32XbdSX/K+h4y+ET03830WVxWAuZEy6wR+CnwgU8MwDjgn\nOq4bmJN35wfyZkpw4HHTTeUK5eab02WXLlUdOtTJDR3qVxKzZqkOGuRk29qqP5xD2n744fI++xTQ\npZeWy06Z4u/DkCFObsgQfx9CrluIbLwvF13k70NJLut1C+X1193f0RRmfamHEsyyJrgx8FVglKqe\npKqZZsdVdb6qfhMYBXwNqGK3VT9EZEsR+aqI/ElEHhKRbhHpEZGDa2irXUT2F5ELReQuEXlWRJaJ\nyBIR+buITKrHbzAaT0g2hI6O3um2TTbxT9ONGNE7BdnTA2ut5e9HJfeENG64obx+443pslOn+utx\nrrgC3n3XfX73XbjyynTZRx7x1+MsWeKvJ+nshNGj4fjjXelbb5wxo/y6zZzpbzskX6IZ0QxcsijB\nzVX1AlVdXssJVHW5qp4PbF7L8TnxReD/gMOBLaPvao1evQdwHfA/wEjgTtzU8EvAwcA0ETm1L501\n+pesD7uQbAjTpsGiRe7zokV+g5QbbuhVgqrw73/7+/Hii/56nPb28vpgjxXAjjv663FWW6287lPy\nV19dXr/mmnTZO+7w15NccUW5YvvHP9Jl//Uvfz1OiGILNaIx14vmoqoSVNW38zhRXu3UyAO4qdxD\ngS2AGX1oqwe4AviQqm6oqv+tqp9W1fHAYbip35NtRDgwCHnYTZjgDCYGD3blzjuny1YySEkjqUCq\nGYMkLTF9lplz5vjrcZKjrqVL02Wffrq87nNNeOGF8vrzz6fLvv56eb2aolixwl+PM2RIeb3k4lGJ\n+fPhgQd6EwH7FNu4cS47RlsbbLWV34imlFFjjz1caYqw8dTsJyiOESKynoi0Vz+icajqpar6HVX9\nh6o+0ce2pqnqoao6q8K+vwO/j6pH9OU8Rv8QYjbf0eFGGt/4hit9yippebjxxumyjz5aXn/sMX+f\n584tr/sytb/2Wnk9qWTijBtXXvc9zJNTmqVRbyWSytSnXJNTlNWmLG++2V+Pc9RR5fUjPP+hI0a4\nqWlw1qTVpqgBJEOCn7lze++3+fP9fzujfwhSgiKyloicLCJ3AcuA53GGL2+JyAMi8mMR2aAeHR1A\n3ItzHdmo0R1pNeoxjVQym4fqZvOdnbDNNnD22a70rT+tvba/HueJxGvZ44/7+5wcST33XLrsyy/7\n63Guu85fj5Ocsk2uPcZZvtxf9+0rTXWmkVTqPiWfXDv1jc7/9jd/PY75FA5sMitBETkIeBQ4FdgB\n52NY8hVsAz4AfBd4RESOSRwrIrJdTn1udkprjs80tBctRqjxQVZH9cWLe6fQurv903qXXVZuDHL5\n5emyyZHDmmumy4YoKoCXXiqv+9YE304sQrz5Zrpsco3s+uvTZZMjzGQ9zquvltdfeSVdNnQ6NKlw\nfApo1Kjyum90HjKdvWJFuWGTT8mHRNoBWz/sDzIpQRH5JPB3YA1gPvAtYDKwdbRNBk4CHgRWAS4R\nkc9Hx7YDfwOquOgOfERkPeBonNHNFQ3uTksR8rZdclQ//nhX+hRhyHpOyKjjkkvK65demi6bnP6s\nNh0aYhiTfHj7LDNDKL0MpNVrJbmm51Molfb75H/wg/L6ySenyyZfFnwvD7/4RXn9vPPSZRcvdi9b\n4BRmtTBvEye6lz5bP6wfVZWgiKwD/C6qfk1Vx6vquao6Q1UfjrYZqnqOqn4Q+F+cEjhPRMYA1wCf\noHZrzAGBiLQBlwOrAzerquc92gglJK5lpQglPkoPpVKZRnIElqzHSa57+Uz9kxabPgvOUJIPTp/i\nbgZK63Bp9SRJ5fTGG+my8+aV1++6K102ZHS+7rrl9XXWSZcdNQoGRU9dEf/0e2n9sLvb1g/rSZaR\n4FeAVYHvqeqvqgmr6nnA94FhwDzgo7hpVM+7cEtwEbAXsBgXYcfIkZC4lgcc0LvO194O++2XLjt3\nbnmuO9+D5tZb/fU4N93kr8dJTr35/P6M/iFp0Zqsx0layibrcebNK59S91n2hrRr1E4WJfhfOB+4\ncwPaPRd4GRiOmyLdQ1U9NmEDGxE5DzgWZyS0t6p6jMCNWsmaYueNN3qnxZYv948OQpyzQ9wCkmtx\nPsd6o/lIxkZNGi/F+c9//PU4IZa9yZFrsm7kQ5aJl9HAHapaZbKoF1VdEQXQ3h+YpKpVlvoHLiJy\nLm60/BxOAVax7evl1FNPfe/z5MmTmTx5ct7d63dCgh/XizPOKK+fdRZMmVJZdtq0letHH11ZNmmp\nmNc6mNF8lKbT0+pxkta5PmvdEMve5JR/liDezfD/lyfTp09nui/SRA5kUYKrALUsyXYBK1pcAf4U\ntwb6AvBhVfXYkK1MXAm2AiHZDepJ8q3d9xb/7LPldd9DKbk+VW0N0SgGyfvA57CfXC/0uc6EWAFD\nryHNwoXO2GvWrIGvCJODg9NOOy33c2SZDn0R2LSGtkfhlENLIiI/Ab6JmyreJ2tM1Vam3v5SWd0e\nkiM031t80lrSZwpvGH3l978vr//xj+mySbeTatah9XTEf+QR+Pa387MubiayKMG7gZ1FJHMAbBEZ\nBUyIjh2QiMiZUbDtMyrs+zHOJeQV4COqen+/d7AJCbHgLJHVD6qzEzbbzLk9bLaZXxGG+JqFvm0b\nRl94JuE97LuP588vrz/wQP79ycIjj7iR5U9/6spWU4RZlODfcM7wl4rIkGrCkcylUdueOAv9h4hs\nJyJzRGS2iMwGtsc5+Z9V+i76Ps4GwFZRGW/rQOB7OJePR4GvisiUCtu36//LmosQC04Ii6N4wQXl\nVnUXXpgumwxf5QtnlTyn+WIZ/Yl6HMeSLhw+lw5YORZqMlZqkqyju9/9rjzAe9r6+oClWq4lnLKY\nhwsMPRvY3iO7AzAnkr0bkLxzP9WyAZOiPvm2FYljpkTf/y7x/dEZ2uoGbq3SJx0o1CtXWkjuuAkT\nymUnTEiX3WqrctkxY9Jl43KlrZllm6UfzSDbLP2ol2wp92FpGzo0XVZVdc89y+X32itdNiTP5N13\nl8vefbe/H/WEOuQTrGoYo6oqIh8HZuKmOO8SkQW4ZLslM4L1on1jI6X5NPCxqNMNR1Vvw41mQ445\nBjimwvd/AP6QU9eanlBjlxDrtJD8fIsXl9d97gm2zme0AklDrGpP05DIQ5Wi56TFR00GCvCFvRuI\nZIpNoapLRGR74NfAJ3FxQj9AeRQYwaUZ+jvwJVX1xNMwBgqVjF122aWybKjCDPHRSz4AfJFEkrLN\n8SpmGGFkyUoRJySEXNKX8d57w87VSmQOoK2qr6jqp3EBor8O/An4d7RdjnMV2EpVDzMF2DqEGLuE\nWodekYiu6kuImvSv8uWkM4xWIDQ+a3ImxRfvdOzY8noyjVackDyaMPCCfgdHKVTnDP7LOvTFaEJK\nxi6l0Z1vZDdqlAtEvWKFi4/oi4sIK0+z+GJx2ujOMPyEGHlNmgTXXNNb32OPdNmODpc665//hP33\nr74c0gy+wiHUnFTXKBZZlM6CBb1RVZYtqx4DMzSjumEY6YQEHw8JydbV5eLvnnCCK33KdSDmVjQl\nWFCyTlmEpHMJDfibNHZJ1g3DqA9JYzGf8dj8+XDffU6x3XefX7HV4ivcaLKkUjpTRFbvy0lEZHUR\nObMvbRj5EZKgNiSdy/33++tJNtywvD5ypF/eMIx8SFpQ+3wFk6EFkw7/cTo6XOqyCy5w5UCY3cky\nEvw28LiInBISNQZARDYRkVOBx3ERVowmoF5TFgcfXF4/6CC/fPIfb9GifPphGIafkDyTp59eXv/x\nj/3t7ruvmzrdd99sM02NNqLJogR3Ax4DTsEpw5tF5LsiMllE1hORwQAiMjiq7yki3xORW3HK74fA\nImBivX6EEUbIlMUmm/SaaovAxhuny/773+X1G2/09yPEMMYwjMaQNJr50IfSZUPil4bMSNWTLM7y\nc3CxQw/HuUbsBewZlxGRZUA8aE/Jw2UOcJ6qNkX4NMMRYvE5fXp5yKQZM2DLLSvLJm94y4RtGAMD\nn0/i0UfDL2P+AEcdlS4bEgAjxAe5noT4Cf5ZVXcGdgbOBGYBb+EU3rCofAu4HTgdF15toinA5iRr\ngtoDDoBhw9znYcP8WdqTzrmWb88wBgY+6++k/+5VV+VzzmYxoqnFT3AeLpYoACIyHFgDeFVVLX92\nizFypAu/dP31TgH6jFdsjc8wWo/Zs/31OMOH++txQmak6kmfXSRU9S1V7TQF2HhCFplDF6R9Pkcl\nSiPGEsmo9oZhDDySSYJ94dgmTHDLJYMGubJadBlofOCL4JGg0ZyERGro6nI35yOPuBv1zjvTZTs7\nYfRo5/w+dCg8/nj6aDAkhZFhGAODkJikXV3wxBPupfmJJ1zd9xyaONEF1dh6a5g1qzGjQXOWbxFC\n3B6uvx4WLnQ36sKF8K9/pctecUV5FJgrr0yXtQgwhtF6JCM/+SJB/eQnvUpy+XI4++x02RAf5HqS\nixIUkcNF5DYReVBELhSRtWP7PiciZ4vIaXmcy6hMyCJzck5/zpx02dVWK6/7FNujj5bXfalcDMMY\nGIwY4a/HmT7dX29GpK8p/0TkOOBXwJ3AC8B4YDVg/8iIBhH5L2Cqqgbl9GtlRCT3dIvv/OT/GPbd\n/821TcMwjHqzzTjNNB0qIqhqrgsteawJHgGMVdUnSl+IyH8DU0TkKFW9F5dn0Kgzw446FHinqtyy\nZXDZL1/icy+fU/9OGYZheOg8/JvM+k0Dl0/6mpoeODvl+7WAv+BGhvsC3X09Vytt7tI3hlmzVAcP\nVgXV9nbV2bPTZW+4wcmVthtuSJcdP75cdttt/f2Iy5Y2k/XLNks/mkG2WfrRDLL1bHu99crl1lsv\nXfYvfymX/etf/X0OJXpuVn2+hmx5rAkui0KmbSAiO8WU68vAZ4CDgAk5nMfIiZD1w1tvLa/75vj3\n3NNfNwxj4PHGG+V1X7Le0OTXAyV2aDV+DZwLXBtt76GqPap6KvAMsGLlQ41GEBLp/XOfK68fc0y6\nbNLAxmdwA86XyFc3DKPxJJ8Pq66aLrv77uX13XZLl22W2KF5OMs/q6pfAz4N7JAiczGwTV/PZeRD\nSJLMVVd1I0Zwpe8fIGk1ttZa/n6EJAE1DKMxbLRRed0XRP/228vrd9yRLjt/fm+w7QcfbFwC3tze\nvVX1UVXt9Oz3ZKwy+pMQn8LLLiv3+7n88nRZbXDkB8MwstGWsNMfnFPYlE039dfjjBrV+4I9eLDL\nWNMIbAKqgISsCZYc5Uu84zE+Tc7/v/BC7X00DKN+dHeX15Oh0eIkZ4p8M0d77umivwwa5MrJk9Nl\nFy/uPW93Nzz1lLfLdcOUYAEpBa6dMcMfXg1g/Pjy+rbbpst+8IPl9W2qTICX3gLT6oZhNJ7DDvPX\n43R0uEgwd9zhSt+zZcBmkRCR7upSACwHXsRlnPi9ql4dei6jfpRSKVWj9GZXijPqe7NLvl0m60ls\nTdAwmp8Q61DI/mwZyFkkJOM2BBgJ/DfwDxGZkkeHjf6lowNuvhkuvNCVvht19dXL68mQa0mSI78h\nQ2rro2EYYSQzvCQzwMR5/HF/vS9kzWtaT4KVoKoOAs4B3gB+CmwLrAmsjnOMPxvoAn4ObAIcjQun\ndlSUnd4YQIRYkibXFaqNBEOVpmEY+bD22v56nP9NRGL82tfy708jCVaCInIM8HXgI6r6HVW9X1Vf\nU9UuVX1AVb8LfAT4KrCvqv4J+DhudPjZHPtu9AMhlqRJn8Jjj/W3nUy1ZH6ChtE/JF2dfK5PcctR\nkepr983gAB9CLY+dLwEzVTXVFVpV5wIzgROi+hzgXmC7WjrZF0RkeJTl4hcicruIvCEiPSJybfWj\nU9v8gIicLCI3ichjIvK2iLwqIrNE5Gsi0jITeyGL10uWlNeXLq1v3wzDqI2QyC6jRvVOnw4Z4ndl\naBYH+BBqUYJjgGczyD0LbBWrP47LLtHfvB+4DDcy3RV4H9BXj7Z/A6cCE4ElwJU4A6DxwC+A2SKy\nRh/P0RSERJd58sny+uLF/ra3395fNwyjPqy/vr8eZ8GC8pyivnyCITNHzUItSnAZbh2wGttGsiWG\n4NYK+5su4He4UekE4Iu4qdm+8BBwLLCOqk5S1c+o6oeBrYEHcb/9F308R1MQsiZ4yCG9xi1DhsDB\nB/vb/uQn/fU4SeVr64eGUU5yOSHpEB/nuOP89VoZNw7GjHFTqFtt1Ti3hxBqUYK3A2NE5IdpAiLy\nA5xCmBn7ejNcDNF+RVUfV9XPq+pF6vIbLqt6UPU291HVP6jqW4nvnwK+gFOyh4pITnEYGkfIm11H\nB4we7f4ZR4+ubvGVDKvmC7OW9FdM+iQaRtFJuhj5DNO22MJfjzN2bPnL7dZb+/vR3e36Us0wrlmo\nRQn+EKdIThGRh0TkVBE5RkQ+KyKniMiDwGm4xHanAojIJsA44Lac+t3M3BuVwwBPDuaBQcia4Ny5\nsHCh+wdYuBDuvNPf9ltv+etxXnnFXzeMovO+95XXhw9Pl03u88kuWADvvus+v/uufzp07ly3v6fH\nldWeAc1A8EhFVe8TkQNw62xbAScnRAR4DjhSVf8Tffc2sA+wsA99HSi8PyrfBV5uZEfyoJ4Orckw\nSb6wScl9Tz+dXz8MoxUICVYxYYIb0T38sJu23Hnn+vatmanJKF1VbwU2x/kATgFujLbfA8cAm6vq\nzTH5F1T1FlUtgr3gd6PyOlVd3tCe5ERWh9YJE9zIsa3NldX+scaM8dfjHHSQv24YRWe99fz1JCUX\npaSrUpIy+EwuAAAgAElEQVQJE9wsUFubK33/16VnwODB2Z4BzUDNa1aq+jbwp2gzABH5LHAo8Cbw\n/cb2Jj86O2HqVDjgABg5Ml2uowNmzco+aiyFZCu9jfpCsiXx/eMOH14+tbrKKtnbNYxmoa2tfDTn\nM3QBuPhi+OhHe+u//W267Ny5bpoTXHnnnbD33unyJaObar68oc+AZqDPhhsisgGwYVRdqqq5Gb+I\nyE+BA2s4dK88+5EFEdkb+A3QA/yPqi7qz/PXi85O2Hxzlz1i2DB47LHqijBL3MCS7Ny52f5hkiHV\nfA67q61WrgTNktQYiIwZU26IVs0gZeJEl+vv6addOXFiuuyLL/rrcebPd2v83d3uhfXBB/3/4yHP\ngGag5hgdIvJ5EXkY5yc3N9qWiMhCEcnJ4JYNgC0Dt/cD/ZqPQER2B66JzvtVVf1Lf56/nkyd2ps+\n6Z13nK9gnmSdav3iF8vrX/hCumwye4Uvm4VFrTGalQ99qLyezNqeZNq03rXyp5+G6dPTZUPW40tu\nD21tA8ftIYSaRoIi8nvgSJwRjAKlZLojcYroIhHZTVWP6UvnVPXI6DxNi4hMBK7HOeF/S1V/nfXY\nU0899b3PkydPZnLIfGA/ccABbgRYGgnut19j+vFywsTIZx26cKG/HqetrTzmabUpJ8PoL5Kjs+T/\nQJJK9/2BKfNoIevxJaqtHdaD6dOnM92nzfNAVYM24NO4Kb9ngeOBobF9Q6PvngG6gcNC26/3hjPm\n6QGuzaGtXYDXo9/6ncBjdaCwdKnqxRe7slH85S+qLne92/7613TZXXctl504MV22ra1ctq0tXTYu\nV9r6W7ZZ+tEMss3Sj3rJTpxYLrfbbumyqmH/I0uXqg4Z4uSGDPH/b8+apTp4sJNtb1edPdvfj3oS\nPTf79NxObrVM/nweZ/6/lzoH9Pecz1V1mapeBOwNrAD+p4b2BwQisjNwA7AKcIqq/qTBXaobI0e6\niBK+tcB6ExLrcMIEfz1Ocq3R0jkZzcLbb/vrSUIyQyxe3Gt009NTfTq0GZLf1otalOC2wHRVXZAm\nEO2bRrbwak2JiOwUrW+u9DtFZEdc/NAO4Eeq+uN+72DBOOSQ3iC+Q4fCJz6RLnvLLf56nKTRTDK9\nU5zkeqGtHxqwspGWz2hr443L675g1ElrTZ/1JoS5J4wb57b2dqfUfIqt5Cs8Y4YrB4LFZwi1rAkO\nJ5sT+Mu4dbKGIyJXAaUQsetE5e4iMjsmdrqq/itWH45b36wUbPtGXP7EV4BRnoTB31DVAe8wH0JX\nl7MmGzcu33+WUki2hx+uHpJt4kR44IHyehpJ94kQd4pGrJEYzUd7OyxfXl5PY/To8kAPo0eny37+\n83DOOb31avE9Q9wTQoNgDDSLzxBqUYJLgZ1FRKI52pUQEQF2otdgptFsi0vwW0JxSiz+rrQOK6NU\nVoKrR9+vARyVck4FTqEFosZkpZRGpfSPleWtMavSnD8fFi1yUzePPuo3095vP7joovJ6GltsUZ4p\ne/PN02Xb2srjM9pI0IDyfHvgV4LJzCq+TCsvveRetFTdvVbNMAbCXZRaVbGFUMu/8b9xwbB/JiIr\n2dKJyCBcdvnRuDWzhqOqm6lqW5Xtj4ljbou+X+lFIUNbbao6WF1A7cIQmkYlJPdYyLpE0tJtq60q\nywHcf//KvyGNNdf0141ikhzN+UZ3SXyzCaNG9a5Rt7f7p05LdHY6J/nOBg4/ipBU9yfAa8D/Ao+I\nyI9E5NgoiPaPgEeAbwCvRrJGQQhdQA/NUJF1XSJpUT1jRrrsPvuU133rLq++6q8bxWTddf31OGPH\n+utxQvL4QW9gi+OPd2UjFGEhkupGo5v9cNOimwHfAy4GLok+j8Y50O+nqhbmuECELqDXYnVWeQK+\nnMmTy+Mi7rFHuuwhh5TXfTkNk1NHu+6aLps0uLGoNQOLECOoN98sr/uyoWy3XXl92xxNB+sd2CIL\nRUmqi6rOwUVmOQq4FDdF+u/o81HAlqo6N69OGgOHrBFgSrJZlWbIG+ZTT/UqS1V/xokdd+xdw2lv\nhx12SJdNuk+UrFUrEbJOZDQfIa4zJyfy6HzfEzV43Dh/Pc7YseUW0dXCppUCW0DjAlsMRHeKvgTQ\nXoZLp3RZft0xikbWxflKb5h5LOovWNBr2bd8uZtySvOHTE5zrVPJlCpizBhnqVfC9wAbOrR32gt6\nH2RG41hzTXjmmfJ6Gs8kohQ/+2y6bEgev8WLe6MZdXe7lzufr+7IkXDffTBlChxzTGP8euuZeq1e\nmH2bMSAIecMMeYMOSeybPKevD8nwa75wbOuvX16vlgLHqI3kA9k3Rb3hhv56nGuvLa9fd11Yv9IY\nNap3BmHw4OqGMV1dcOihzq3i0EMbtx4XMhvUDFQdCUZZ4WumaBaSRn0IecOsZFCQ9lYc8mY+Z46/\nHieZ0DQen7RaH/JM/dTqaaXWXLM8juxaa6XLbrhheXxN30gpJEv77rvDVVf11nfbLV02dCRYmqVY\nsaL6SLBesyWtTpbp0Cep7CuXBc14DsOoSj38msaOdes9777rSt+occsty+s+14uQSDShRjSDBpX7\nK/pGmSHRTOpFydetRJ7+lUmjlGQ9Tr2mLUP6EJLRvRTVZcGCbOtrpdmSrPKGI4uCeoralaBh9Dul\n8FEPPeQeOL4HzYIFTgGCK32jxqQRg+8hc8ABcEPMS9ZnpFA6f1o9ydCh5XEkfUYb664Lr71WXk8j\nNIlrM5Dso0/BrrFG+bXwrfO98IK/Huf11/31JFkzutcS1WWgrcc1A1XfyVR108jZvKatP36EYcQp\nhY+6/XZX+h4GIWuCyYebb83lkkvK65demi6bHKFU8+9KKl9fvsTkOpJvXalesVGTlrLJel/YaKPy\nuu/3ffWr5fWvfCVdNjnK96UaClGYpYzuPT29Gd19hK6vDbT1uGbADGOMAUNIJIqsD4OQaa+PfrS8\n/pGPpMsmp219mSySodq22CJdFlZer/JZk4ZkEE8afyQVjO+cvj4kR9YbbJAuu8Ya5fVqUXmOSWQs\n/exn02WTLw8+R/VvfrO8fuKJ6bLJ6+S7bkbzYUrQGBDUKxJFaeq0ra165P1kuhmf/+GkSf56nM99\nzl9PkjSySRrhxEn6MvqUVYj1a/Kc8TXKJElrV58SDFFU4OLI+upxfvSj8voZZ6TLhiRxTr60+F5i\nQn3/jPpjStAYENQrEkW9pk5Dcrvtu2+5w75vhAkrKzKf036Iq0ZyFJwcccZJGvr4DH+Sa1++qD/J\nKDy+qDwAt93mr8dJKt+ka0qcek2TV/L9MxqLKUFjQFDPSBT1mDodO7Zcsfne+Bcv7h1JqVZ/MCbX\nAH1rgkl8CihpkONz6whJQVUK5VUiHhggSTL4tC+rB8CBB/rrcZI5KH05KUNI5ro8+OB02ZA8fkb/\nYErQGBA0Q2LPkKmsefPKI9HcfXe6bOiDMTmqHDEiXTY5avSN7pJTfr7UPaNGldc33TRdNmStMZnF\nw5fVA8Isa5MuLu9/f7psyAvPyJEuG8l3vuNKny9fM9zHRjmmBI0BQ6Mt3xYv7l0L6+nxj9iefNJf\nj9PR4YIdX3CBK6v9vuTDO/lwj7P99v56nOS0qm+aNTk96PON68vUaTU3gp128tfjnHdeef2Xv0yX\n3WSTcleGZEb4OF1dblT5s5+5stp6daPvY6McU4KGkZGQKdnkFJlv6q2rCz78YfjCF1xZ7SGaHLH5\njDZ23LG87gsQfvjh/rqv3WQ9TnJk5BspxSOvVKonCbHMTLqiJOtxpk8vD8LuS8c1d64bsXZ3u7Ka\n24PRXJgSNIyMhExljRzpDG2OOMKVvgf/tGnOSb+nx5XJfIhJ/vUvfz1OyLTeRz5Svo6ZzLUYJzkF\n6wtXlhz5+a5b0g0lWU9SivgD1SP+JC1/faPGyZPL6750XCFGNEbzYUrQMALIOpXV2eliSF52mSt9\nDvAhU6cQFr4thJB1zJCp06TC9FnKJv38jj46XRYqR/xJ46ij/PU48Rij4MKcpRHyomE0H6YEjZYk\nxLG+HoQkOP3oR8vXn/bd1992yBRnCCHKePz48rovOey99/rrcf761/L63/6WLhtKiJ9nyLUo+ZoO\nHlzd19RoPkwJGi1HvRzrQwhJcPrSS73+e4MH+60yIcwHMWSqLkQZ77mnm4psa3NlcvowTtKoxBfa\nLOlY7/Plg/o5n4es6ZZ8TWfOrO5rajQfpgSNlqNejvUhjBwJjz0GF1/sSt+aYKgPZMiD//nn/fU4\nL73UGy+0rc2vjDs64Kab4MILXel78B92WHn90EPTZY87rnyN79hj02XBWeyWpkOXL/db7IZEBwpx\newCz+BzImBI0Wo56OtaHMHKke6hneYCG+I6FRB0JGdGEJHHt6nKj2xNOcKVvtB2yvtbR4Rzk29pc\nWe1axNcbVf0GOiHRgULdHoyBiylBo+UYiA7JISOJEOf6kSPh8cfdiPTxx/0KuVIS1zRCRttJJZis\nJ9tdtMgp90cfrT6KD3FlgOzX2dweioMpQaMlqdf0VGcn/Pa31dMd1bPdUOf6jo5s+eVClOu4cc5p\nf9AgFzDaJ5tMQ+RLSxQ6ig9ZezWMSoj6ggkadUNE1K79wKKz08W2XLbMTS1WG1mFtLv55s6KdNiw\n6muIJcOfUvJU32g3RLYknyUpa8i1eOSRcjeOhx/2R7l55BGXf/HYY/1y8b5cf71TgHn8PcBdh4kT\n3ah1zBgzeGkWRARVrRJHKAwbCRpGRq64ojf487JlcOWV+bQb4k4BYVORtRgJZXk3mzq1/Fr4+nz1\n1eX1a69Nl+3sdO4XZ5/tyiwj46xrryGYxWdxMCVoGBlJBolOBpFOktVXMXRKL2TKMES2q8tZUO6+\nuyt9/e7LNKRPyYa+ENQTs/gsBqYEDSMjIb5xIb6KIe4UEGb4EyIbEr5t5Ei47z7nQnDfff4+f/jD\n5fW9906XnTy53FfRF67MMPKg5ZWgiOwqIheKyBwRWSoi74hIl4jcJyJniYjH1ThT+/uKyFUi0hm1\n3Ski00TkxLx+g9EcdHTAnDnOxH7OHL9CCZ2GDJ3SCxmlZJUNiZLS1eX8/c45x5U+JT9vXnn9nnvS\nZV96qVcJDhpUPXCAYfSVljeMEZEfAd8DngQeA14A1gJ2isrngUmq6vFeqtiuABcBxwHvALOBZ4D1\ngW2A11Q1dVnfDGNam9JIcMECN2rM21Wjq8sp2nHj8mu3sxM228w5nw8ZAk88ka6UZ892o7QVK9xU\n64wZTtGmtZvV8Kfe180Y2NTDMKYISnBL4C1VXZL4/n3ApcCngNtUdc/Adn8CnATcAhyhqs/F9gmw\ng6rO8xxvSrDFyWppWUu7oRafWRRmV5dTZCWLSN9oN9SiNcSCs17XzRj4mBLMGRHZCHgK6AHep6rL\nMx63NXA/buS3tap6UoqmtmFK0KiJkFFYiMIMaXf2bNdud7eLLjNzZrqsYeSFuUjkz4pY2R1w3Am4\na3dJLQrQMPpCiMVnyNpkqNVpVsd6w2hmCjsSFJF24ELgWOBqVT044NiFwPuBfYD/AIcDY3Frg/cA\nV6jqO1XasJGgUTNZpwxD19hCpiJt2tLob2w6tA+IyBbADwAF1sYZxqwD3AkcpKrPZmxnCPB2VP0S\n8COcgc17IkBn1OZdnnZMCRr9gikro1UwJdgHRGQ3IBle9xbgf1T1yYB21sOtBSpuGvV+4CtRuSlw\nBvAxnBXqB1T1xZR2TAkahmEEUDglKCI/BQ6s4dC9VPUZT7sbApNwo7j1gSNVNVMQLBHZAFiKU4Kv\nAaNV9dXYfgHuxblJnK6qp6W0Y0rQMAwjgHoowcF5NlYHNgAyhNAtQ4F2r4DqUuDPIjILt6Y3RURm\nZZwSjbsFXxlXgFHbKiK/Bc4H9gQqKkGAU0899b3PkydPZrIvBIlhGEbBmD59OtN9oYtyoKlHgv2B\niFwL7A8cq6p/yHhMyeH+h6p6RoX9/wX8E3hEVSsmjrGRoGEYRhjmIlEfXojKdQOOKQV+GpGyvxSK\n7Y2aemQYhmH0C4VWgiLSBpRC9C4KOPRKnBVoWpSZfaIy1TrUMAzDaDwtrwRF5NsistKITUTWAaYA\nm+OixtyQ2L+TiCwUkQUVmv0DsAT4oIickjjuEJzf4AqcH6JhGIbRpLT8mqCI9OCiwdyPC6DdDWwE\nbA8MA54FDlDVexPHTQKm4Wxd2iq0OwG4EVgVN4p8AOcisUN0jq+qaqoStDVBwzCMMArnIpEHIvJF\n3JTntsB6wCo414aFwHXARar6eoXjJgG34pRgRStaERmFc8DfF7em+BpwO3Cuqs6q0i9TgoZhGAGY\nEmwhTAkaRmXqkSbKaA3MOtQwjJamFOt0jz1c6UvWaxh5YErQMIymISTrhWHkgSlBwzCahpB0ToaR\nB7Ym2CBsTdAwKmNZL4w0zDCmhTAlaBiGEYYZxhiGYRhGjpgSNAzDMAqLKUHDMAyjsJgSNAzDMAqL\nKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzD\nMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgS\nNAzDMAqLKUHDMAyjsJgSNAzDMApLoZSgiOwrIleJSKeIvBOV00TkxBrb21JELhORpVF7T4rIr0Vk\n/bz7bhiGYeSPqGqj+1B3RESAi4DjgHeA2cAzwPrANsBrqrplYJuTgOuBYcA9wCJgPLA18AKwm6o+\n6jlei3DtDcMw8kJEUFXJs83BeTbWxJyFU4C3AEeo6nOlHZGC3CGkMREZDvwVpwC/rKoXxvb9FPgm\n8Bdgp7533TAMw6gXLT8SFJGtgftxI7+tVfXNHNr8MvBL4BZV3SexbxDwMDAa2F9Vb0hpw0aChmEY\nAdRjJFiENcETcL/zkjwUYMTHAAX+nNyhqj24UaIAH8/pfIZhGEYdKMJ0aGmkdruIrAUcDozFrQ3e\nA1yhqu8EtrldVN6Vsv+uhJxhGIbRhLT0SFBEhgDvj6pb4qYpzwOOB74G/BF4VEQyr92JSAewVlRd\nnCL2VFRuFtpnI5zp06c3ugsthV3PfLHr2dy0tBIE1sRNS4JTfk8CuwEdOKvQa4CRwFQRWTtjm6vG\nPqdNr74RlR0hnTVqwx4y+WLXM1/sejY3TT0dGllaHljDoXup6jOUK/k3gX1U9dWovkBEDgbuxSnE\nLwGn9aW/hmEYxsCiqZUgsAFuGjMEBdqjz12x76+MKUAnqKoi8lvgfGBPsinBN2KfV0mco0RptFhp\nn2EYhtEkFMFF4gXcGt4PVfWMCvv/C/gn8IiqjsnY5kvAGsB4VZ1fYf+BuKnWeaq6c0obrX3hDcMw\n6oA5y4dzD/BhYETK/tJa4Bsp+9Pa3AvnDL+SEgRKiu/etAby/kMahmEY4bS6YQzAlTjjmD1T9pdc\nKNLcHSpxTdTmZ5I7Imf5w3DTslcFtGkYhmH0M0VQgn8AlgAfFJFT4jtE5BCc3+AK4MLEvp1EZKGI\nLKjQ5hTgWWBPETkhse9sYHPgnrRoMYZhGEZz0PJKMHKE/yRuuvOHkWL7u4jcBfw/3Ijta6p6f+LQ\n4TijnK0qtPkmbrT3NnC+iNwlIn8WkQeBbwDPA58WkTNFpCfaUjNViMjhIjJTRF4Vka6ovROiuKYG\nUO1aisiU2P5KW6WXmcLQl+tj92c5tVxLuz+rIyLDROQkEblTRF4RkTdFZJGI/FVEdk05ps/3ZhHW\nBFHVuSLyQeAHwL44t4vXcNOV56rqrLRDo61SmzNEZDvgh8DewDjgOdyI8nRgE+BbQA+9voorISIX\nAF/EKdRbgOVRe+fj1h0PCfmtrUgUzKDqtcT9re4AKmXv6KxD1wYawdfH7s9UarnX7P5MQUQ2A27E\nxVx+BrgVN0M3CjgIF/95duKYfO5NVbUt5w0YAiwAngb+AXQDJ1aQ+wTuwb4EGB37fh3gwei4rzT6\n9wyQazkl2ndUo/vcjFst18fuz1yvpd2f6ddmOO7FYAUuA48k9q8JbJH4Lrd7s+WnQxvEj3DTqF8A\nXvfIfRf3dvhtVX289KWqvoB7wxHgO3Xs50Ag67U08sfuT6M/OBk3AjxfVc/RSJuVUNVXdOXcrLnd\nm6YEc0ZEJgAnAper6j89chsC2wPvAlck96vqDGApsL6I7FKn7jY1Wa+lkT92fxr9gYi043K9KvDz\njMfkem8WYk2wvxCRoThr1JdwAbp9lDJMPKiqy1Jk7sKlbdoOmJNLJwcIgdfyvcOAvURkPC5qz3PA\n7cBNybfLghJyfez+9FPLvWb358rsgPPhXqKqT4nI9rgUdOvirs+NqnpH4phc701TgvlyJi5rxadU\n9ZUqsqUME2mZKKDY2ShCrmUJBY5MfCe4OLGHaYXoPgUj5PrY/emnlnvN7s+V2SYql4rIz3DW9fEX\ngpNF5GrgM6r6dvRdrvemTYfmhIhMxI1YrlLVlYboFSjFF/Ul+n0D909SqGwUNVxLcNF5vorLFbkq\nLjvIAcB/ou9uEpEN6tDdgULo9bH7M51a7jW7PytTSku3PW7p4+fAFjhjmI/hDF8+Bvw6dkyu96Yp\nwRwQkWHA74FXcZnsjRqp9Vqq6i9V9QJVfVhV31bV51T1X7gQdnNw0yvfrUefBwJ2ffKjlmtp1z+V\nkg4aDPxJVb+pqk+o6uuqOhXnHgFwZORGUbcOGH3jLFyUmG+o6vMZjynFKl3FI7MqbmqgSNkoarmW\nqajqiqhNAfbra3uthuf62P0ZSC33mt2fZffOJcmdqno3cDfu+kyKvs713rQ1wXz4OM5n5bMi8tnE\nvlJmihPEZZd4VFU/j0vwC84ZNI2No/JJj0yrUcu1rMbCqNwwny62HJWuz5NRafdnGLXca0W+P59I\n+ZyU2QFYP6o/GZW53JumBPNjELCHZ/9m0bZaVC9lmPiAiAxNsXLaKSFbFEKvZTVKGURCMoUUiUrX\nx+7P2qjlXivy/Rm/d0bgXBuSJDP95Hpv2nRoDqjqZqraVmkD/hiJfTP6bofomCW4lExDcLFNyxCR\nScBGwDOqOju5v1Wp5Vpm4FNReWf+PW4JVro+dn/WTC33WmHvT1XtBOZG1b2T+0VkDZzRDMC86Jhc\n701Tgv1HpZiXpbWAs0Vk8/cERdbFWUMp8JP+6d6Aouxaish4EdlfXBqr+PdtIvIN4Cu4a/mLfuxj\n09CH62P3Z4JarqXdn1U5A3effU9E3nuxjXyFfwOsjktQHvf3y+3etOnQ/qFi0GdV/YeI/BoX5ucB\nEbmZ3iCwHbgA3xf0Wy8HBpWu5aa4a/WyiNyDy+IxAueDNBIXR/AkVb25vzrZZGxKDdfH7s+KbEr4\ntazlmMKgqlNF5Bycj+AsEZmDC5KxM+76PI1LeRc/Jr97M0uAUdv6FBy2FDh3paDPMZnDgJk4t4Au\nXLSDLzS67822pV1L3EPm57joG08Db+F8iB4GLga2a3TfG3zd+nR97P7s27W0+zPztf04cHOkAN+O\nrs9PgRGeY/p8b0rUkGEYhmEUDlsTNAzDMAqLKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNAzDMAqL\nKUHDMAyjsJgSNAzDMAqLKUHDMAyjsJgSNFoSEdlHRKaIyEIReVVElonI8yJyu4j8VER2ytDGfiLy\nJxFZJCJdIvKOiDwtIlNF5HgRWbVaG0ZjEZFvi0iPiHwk8f2qInJRdE+8JSLTRGR7Tztfj9rZxyNz\ntYi8LiLr5fkbjPpiEWOMliIKoPs3XAJOBR4DHsSlYRkBbIfL4i24TNZHV2hjHeDvuHROCjyEy/n2\nLi46/U7AUOBFYAdVfbq+v6r+iMiTwCbApqr6VIO7g4hMx13/yao6o8Y21gcewQVf3iux7yrgY8B/\ngCXAR3GhurZV1ScSshvj7qFrVPVIz/m2idr7vap+rpY+G/2PBdA2WgYRWROYjYvVOBP4sqo+UEFu\nV+DbwNYV9q0OzAJGA3fg4hA+mJBZBfgC8H1gTVw8yIGORluzkEd/TsdlHz8t/qWIbItTgDeo6n7R\nd0cBvwdOxGV1iPMr3AvQid4Oqz4QKdejReQXqjq/j/03+gGbDjVaiQtxyXZvB/aupAABVHW2qn4c\nOKHC7guAzXE5zvZOKsDo+DdV9VxgR+C5vDpv5IeIrAUcATymqrcldm+PU7B/iH33J+AdYNdEOwcB\n/w18S1VfyHDqS3HP1aQiNZoUU4JGSyAiWwCH4B5uX1TVFdWOUdV5iTZG46LSK24E+G6V4x9X1SAl\nKCL7isiVIrI0Wqd8JlqnPCnKn5aUP1JEpovIyyLytog8KiLni8hGKe33iEh39PlTIjIrWs98XURu\nFpHdEvJHi0gPbipUgCejNnpEpFtENknIjxGR34nI41F/XhaRm0TkwJT+lNrbJFqnvSVao31TRGYn\njxORSVF/JkX9mR7rT4+I7JHxUn8OGEZvIuY4pUzur5S+ULcu9Fp0TKkvqwLnAbep6pSM570B92J0\nuIislvEYo4GYEjRahf1x9/N9qrqgxjYOiNp4QFXvz61nESJyIfAv3FTcEuAK3BrSRrgkoesl5C/H\njVZ2wWUdvyradQLwn3gC0grnOg24HFgGTMVN2e4F3CwiE2Kij+KmAd/CKf8ronppeyPW5mFRfz8b\nfX8dcB+wO3CNiJxaoSulac3jcApiOPBP3DrrzsBVInJwTP7Z6LzPRcfdEOvLlGh/Fj4WHX9LhX1P\n4hTsmNhvWxNYB4ivB54ZfXd8xnOiqj3AdNzv/HDW44wG0ugcUrbZlseGe+PvBn7bhzb+APQAF9eh\nf1+P2l4K7FRh/2SgI1Y/IZLvBMbEvhfg/6J9jwPtiXZ6ouvwAs7II77vomj/vyuc/4nouE1S+r8N\nbrrwNeAjiX1bA4uj4ydVaLcHZ3SyT2Lf96J9D1c437SovT1quNbvi/r6DjCkwv61gddxynBbYC3c\ndAH2Ni8AAAW7SURBVGg3cFwksyOwAvhhDef/avS7ftXf/we2hW82EjRahbWjsuK6TcxlIrnFp/vW\nwY0ens+zYyLShjOiUeCzqnpXUkZVp6tqV+yrb0TyP1DVhTE5Bb4FPAWMwk0BV+KHqvqfxHcnR+WH\noj6F8AOgHbc2dmOi7w/hjEYE+HKFYxX4parelPj+ZzilukXa9G6NjAWGAE9ohSltVX0Rp4A3Bu7B\nWfkejltL/p2IDAJ+i7MsPat0nDiGJdurwENRmepyYTQPpgSNojAWOCq2HR2Va/sOyokdcetQSyoo\ngpUQkQ1xBj49wGXJ/aq6HDfVKbgRZCX+WeG453HrYEPpXReriogIsC9Omf0jRazkxrBryv5K/VmO\nG80CjMzanwyUppVfShNQ1fNx07jnAr8BjgX2il4yvg6MB45X1eUiso6I/D9cNvi3ROQxETncc/7S\nec1fcABgLhJGq/BiVK5TaaeqnoczcgBARJ7AGYPEKY0i1825b6Oi8uGM8htG5TOVRjIRjydky9B0\nX7/XgTWIGYBkYASwGk4JvuB0YioVrz9u5JrWHwL7U401Em1XRFVn41xq3iOaGTgNuFRVZ0ajwn8D\nW+FGw4/ipjsvE5EuVb2uQtOl865RYZ/RZJgSNFqFe3Am8Tv2oY27gSNxzvB5Uqu/W7P47ZWmTrup\nMDJNkNbnnvy6U5VXo7IW68zzcSO+k6L6Prh1wx+p6s8BROQ2nGHTd3DGQUlK5321wj6jyTAlaLQK\n/8RNbY0XkbFam4XoP4GfA9uIyHhVvS+nvpVGQVtllF8alSNFpD2aNkwyOiFbT17EGbYMwwUgeKsf\nztkXSm4rmad8AUTkEzgL4c+oasl9YjxOsc8pyanqayLyULSvEqXzmg/pAMDWBI2WQFUX4cz7BfiN\niAS/4KnqY7iQawJcKCLtPnkR2TwKzVWNu3GKZCNf7MlYP5bipjsH4Ua3yfMOBj6DezhPz3D+LJSm\nXVe6bqraDdwcVdMMcfImtT8ZeDA6fjMRGZLlABHpwE2X36iqf4ntejMqhycOWYX0Ue/YqLwnW3eN\nRmJK0GglTsCZve8G3CoiFd/UoxiPaVNlX8YpoAnANBEZV+H4VUXkWzjlVnX9UJ3j/lk45TpFKgTv\nFpHJ0YO4xM8j+R+JyFYxuUE4q8pNcL/1imrnz0hpRLlSKLmI03EuA78UkU9VEhCRnbIo+Zz6k4qq\nvoMbubWTfXr8LFwIvC8kvr8f93c4OjIQQkR2wfkYpim5XXEKclpYz42G0GgfDdtsy3PDWeTdhlu/\n6sGZuV+F8yO8Eme+XvKluwkYkaGN+Thl8xecFeQ79Pr8bRTQt4uiNrtxYdkuxznPl3zsNknIXxZ9\n/w7OafzPOMOMHpwRz/YVztEDdHv6UNEfEBfmqwfnsvB34OJoWzMmcyjOSb4nauefOP+6G3BO7N3A\nmVnOF9tf0R8QOJBe/8JrYv15f8ZrXfLLrOrnh3PaXwGclLL/lqiP/4nuoa5Ift8Ksm3RtegCVmv0\n/4NtGe6VRnfANtvqseFM+qfgLDJfixTJc7ig2D8Dds7Qxn7RQ/7R6KH2Nm5971qcSf37aujXftHx\nz0Z9Whop3BOp7Nj9GdyU58vR+R/DBXTeMKX9HmCF5/xPRA/wpBIUnO/cg7joMd0pcpsBv8C9GHRF\n2yKcMv8SsH6W88X2T4v2r+QUj4vUck+keEsvD5mc53GWmW8Ci6rItUXn+A/Q5mnr4uj+eRu4Czgg\nRfaA6G/wm0b/D9iWbbNUSoZhtCRRmLr/wQVCn95P57wSF3B7W7UsEgMCU4KGYbQk4pLbPgLco6p7\n9sP5PgjcC0xR1ePqfT4jH8wwxjCMlkRdho8fA3tIIrN8nTgdNz38/X44l5ETNhI0DMMwCouNBA3D\nMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzCYkrQMAzDKCymBA3DMIzC8v8BQT2fmwYx\n2mQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12cdbc250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matplotlib.rcParams.update({'font.size':22})\n",
    "fig1, (ax1) = plt.subplots(nrows=1, ncols=1, figsize=(6,4))\n",
    "ax1.plot(xdata, ydata, '.', color = 'b', label=\"data\")\n",
    "ax1.plot(xdata, y_fit, 'r-', label='fit')\n",
    "ax1.set_xlabel('GC content (%)')\n",
    "ax1.set_ylabel('log$_2$(Q)')\n",
    "ax1.set_ylim([-6, 6])\n",
    "#ax1.set_xlim([25, 75])\n",
    "#ax1.set_xticks(np.linspace(25,75,6))\n",
    "ax1.set_yticks(np.linspace(-6,6,6))\n",
    "plt.savefig('./figs/fig_GC_VS_ampLog_NHMO.png', dpi=300, bbox_inches='tight', format='png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'df' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-ddc52c5d68f7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mseaborn\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrcParams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m'font.size'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m22\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mJointGrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"GC content (%)\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"log$_2$(Q)\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfont_scale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_joint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"b\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0medgecolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"white\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'df' is not defined"
     ]
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "plt.matplotlib.rcParams.update({'font.size':22})\n",
    "g = sns.JointGrid(x=\"GC content (%)\", y=\"log$_2$(Q)\", data=df)\n",
    "sns.set(font_scale=1.5)\n",
    "g = g.plot_joint(plt.scatter, color=\"b\", edgecolor=\"white\")\n",
    "_ = g.ax_marg_x.hist(df[\"GC content (%)\"], color=\"b\", alpha=.6, bins=28)\n",
    "_ = g.ax_marg_y.hist(df[\"log$_2$(Q)\"], color=\"b\", alpha=.6,orientation=\"horizontal\", bins=40)\n",
    "plt.xlim([20,80])\n",
    "plt.savefig('./figs/fig_GC_VS_ampLog_NHMO_dist2.png', dpi=300, bbox_inches='tight', format='png')"
   ]
  },
  {
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}