postmutloc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import numpy.random as npr\n",
    "import scipy.optimize as opt\n",
    "import sys\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import sys\n",
    "import glob\n",
    "\n",
    "from scipy.special import logsumexp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "This notebook interprets and analyzes the posterior probabilities of de novo mutation produced by `run_de_novo_probabilities.sh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>logprob</th>\n",
       "      <th>family</th>\n",
       "      <th>position</th>\n",
       "      <th>nodename</th>\n",
       "      <th>rep</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4000000</td>\n",
       "      <td>-572.620340</td>\n",
       "      <td>F196</td>\n",
       "      <td>6791</td>\n",
       "      <td>m1c1_fblo</td>\n",
       "      <td>2953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4000001</td>\n",
       "      <td>-inf</td>\n",
       "      <td>F196</td>\n",
       "      <td>14573</td>\n",
       "      <td>m1c1_fblo</td>\n",
       "      <td>2953</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     index     logprob family  position   nodename   rep\n",
       "0  4000000 -572.620340   F196      6791  m1c1_fblo  2953\n",
       "1  4000001        -inf   F196     14573  m1c1_fblo  2953"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load in the data.\n",
    "\n",
    "# results/posterior_mut_probs.tsv.gz is the posterior probabilities\n",
    "# from the MCMC samples provided by the authors.\n",
    "dat = pd.read_csv('results/posterior_mut_probs_subsample.tsv.gz', sep='\\t', comment='#', compression='gzip')\n",
    "dat['position'] = dat['position'].astype(np.int)\n",
    "dat['rep'] = dat['rep'].astype(np.int)\n",
    "dat.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>pprob</th>\n",
       "      <th>tissue</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>family</th>\n",
       "      <th>position</th>\n",
       "      <th>nodename</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">F098</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">2377</th>\n",
       "      <th>m1_bl</th>\n",
       "      <td>0.0</td>\n",
       "      <td>bl</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m1_ch</th>\n",
       "      <td>0.0</td>\n",
       "      <td>ch</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          pprob tissue\n",
       "family position nodename              \n",
       "F098   2377     m1_bl       0.0     bl\n",
       "                m1_ch       0.0     ch"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the posterior means across posterior samples.\n",
    "pprobs = dat.groupby(['family', 'position', 'nodename']).apply(lambda x: np.exp(x['logprob']).mean())\n",
    "pprobs_frame = pprobs.to_frame(name='pprob')\n",
    "pprobs_frame['tissue'] = pprobs.index.to_frame()['nodename'].str.split('_', expand=True, n=1).iloc[:,1]\n",
    "pprobs_frame.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# For each heteroplasmic variant, calcualte the maximum posterior probability, across tissues.\n",
    "pprobs_max = pprobs.groupby(level=[0,1]).max()\n",
    "pprobs_max_frame = pprobs_frame.groupby(level=[0,1]).apply(lambda x: x.loc[x['pprob'].idxmax(),:])\n",
    "tissues_max, pprobs_max_by_tissue = zip(*list(pprobs_max_frame.groupby('tissue')['pprob']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Preparing various things for the plot.\n",
    "tissues_max_translated = {\n",
    "    'bl':'adult blood',\n",
    "    'ch':'adult cheek',\n",
    "    'fblo':'early blood',\n",
    "    'fbuc':'early cheek',\n",
    "    'loo':'late oogenesis, meiotic arrest',\n",
    "    'som':'pre-gastrulation somatic',\n",
    "    'eoo_post': 'early oogenesis, post-oogonial divergence',\n",
    "    'eoo_pre': 'early oogenesis, pre-oogonial divergence',}\n",
    "labels = [tissues_max_translated[tis] for tis in tissues_max]\n",
    "\n",
    "oplcolors = '''fblo #1f77b4   fixed\n",
    "blo  #1f77b4   rate\n",
    "fbuc #ff7f0e   fixed\n",
    "buc  #ff7f0e   rate\n",
    "som  #2ca02c   fixed\n",
    "loo  #9467bd   rate\n",
    "eoo_post  #d62728   fixed\n",
    "eoo_pre  #ff8080   fixed'''\n",
    "import cStringIO\n",
    "oplcolorsio = cStringIO.StringIO(oplcolors)\n",
    "\n",
    "colordat = pd.read_csv(oplcolorsio, sep='\\s+', header=None, names=['tissue', 'color', 'type'])\n",
    "color_dict = {row.tissue:row.color for row in colordat.itertuples(index=False)}\n",
    "color_dict['bl'] = color_dict['blo']\n",
    "color_dict['ch'] = color_dict['buc']\n",
    "tissues_max_colors = [color_dict[tis] for tis in tissues_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEQCAYAAACAxhKnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VGX2wPHvGwiQnpBiSAKEFhCChQ6KIkWMsChVBCXo\nImBZ+SEgRaUpIAqiYkFYF3ClCbEgRVyULqCgK7oERCSEHoiUhEAKOb8/JhlTJpM+SeB8nmeeJPd9\n773nvTOZM7fMPUZEUEoppRzFqawDUEopdWPRxKOUUsqhNPEopZRyKE08SimlHEoTj1JKKYfSxKOU\nUsqhNPEopZRyKE08SimlHEoTj1JKKYfSxKOUUsqhKpd1AGXFz89PQkNDyzoMpZSqsPbu3XtORPwL\nO98Nm3hCQ0PZs2dPWYehlFIVljHmaFHm00NtSimlHEoTj1JKKYfSxKOUUsqhbthzPBXdpUuXiIuL\nIzU1taxDUUpdx5ydnQkICMDT07PElqmJpwK6dOkSZ86cITg4GBcXF4wxZR2SUuo6JCJcuXKFEydO\nAJRY8tFDbRVQXFwcwcHBuLq6atJRSpUaYwyurq4EBwcTFxdXYsvVxFMBpaam4uLiUtZhKKVuEC4u\nLiV6WF8TTwWlezpKKUcp6fcbTTxF0HRxU5sPVXD33XcfixYtyrff5s2badSoUYGXGxoayq5du2y2\nTZ48meHDhxd4WQWxa9cu9A4YShWOJh5VoRhjOH36dFmHoZQqBr2q7ToROm5tqS075tVupbZspdSN\nR/d4VIl5+eWXqVWrFp6entx5553s37/f2vbbb7/RunVrPDw8iIyM5Nq1a9a2wYMH8+qrr1r/fvXV\nVxk8eHCu5UdERABQr1493N3d+eGHH2zGsXnzZurVq4e/vz8vv/xynvG+9dZb1KlTh4CAAIYNG8aV\nK1cK1DZlyhT8/f2pX78+W7ZsyX/DKKWy0cSjSkzjxo3Zu3cv8fHx3HHHHTz++OPWtgEDBtCtWzfi\n4+Np3749mzZtKvTy169fD8Dhw4dJTEykZcuWNvt9+umn7Nq1i927dzN//nzrfFlt2LCBWbNm8fXX\nX3Po0CEOHz5sTVL22r788kv+9a9/sXv3bnbu3MnKlSsLPQ6lbnSaeFSJ6d27N/7+/jg7OzNhwgR+\n+OEHUlJSiImJ4cCBA4wbN44qVaowZMiQUj0h/3//93/4+/tTt25dhg0bZjM5rFixgmHDhtGgQQO8\nvLyYOHEiy5cvz7ctKiqKYcOGUbduXfz9/fm///u/UhuHUtcrTTyqxCxYsIDGjRvj5eVF7dq1SU9P\n5/z585w+fZqAgACqVKli7VuzZs1SiyMkJCTbek6dOpWrz8mTJ7PFULt2bWs/e22nT5/OtXylVOFo\n4lElIiYmhueee46PP/6YCxcucPSopUyHiBAYGJjrvnLHjh2z/u7m5pbtHEpxvyF9/PjxbOsJDAzM\n1ScoKChbDLGxsdSoUSPftsDAwFzLV0oVjiYeVSISExNxcnLC39+flJQUJk+ebG0LDQ0lLCyMmTNn\nkpqaysKFC4mJibG2N23alDVr1pCYmEh0dDQrVqzIcz0BAQHZ5rXl7bff5ty5cxw5coT58+fTt2/f\nXH369u3L/Pnz+f3337l06RIvv/wy/fr1y7etV69eLFiwgJiYGM6ePcvbb79d8I2klAI08agSEh4e\nzpAhQwgPD6du3bo0btw4W/vSpUtZvXo11atXZ/Pmzdxzzz3WtkGDBlGzZk2CgoJ48sknrW/ytrzw\nwgs88MADeHt751lB9sEHH6R169a0bNmSxx9/nPvvvz9Xn4iICEaOHEmnTp2oX78+tWvXZuLEifm2\n9ejRg0cffZQWLVrQtm1bevfuXehtpdSNzohIWcdQJlq0aCFFLX2d110Kfon8pTghFVh0dDQ333yz\nQ9allFJg+33HGLNXRFoUdlm6x6OUUsqhNPEopZRyKE08SimlHEoTj1JKKYfSxKOUUsqhNPEopZRy\nKE08SimlHEoTj1JKKYfSxKPKJXslrK9n06dP1zteF1NBt2FMTAzVqlUr0DJzljh3d3cv9j0Fb2QO\nTTzGmL7GmNXGmBPGmERjzF5jzMM2+j1hjDlkjLma0aeTjT7BxpjPjDEJxphzxph3jDGujhmJUqVj\nwoQJvPnmm2UdRjYdOnSwloWoCByxDRMTEwkICCjVdVzPHF36+jngCDASOAfcDyw1xviJyFyAjEQ0\nD5gMbAceA9YYY1qKyK8ZfZyBDUAK0B/wBt7I+PmIIwdUXkQ3Kr1b6Nx8ILrUlp1TWloalStrRXZ1\nY7nRXveOPtT2NxEZICKfiMi3IjIaWIYlIWWaDCwWkZdFZBMwGPgdGJelTx/gZqC3iKwVkSXAP4AB\nxpgGjhiIyi0mJoaIiAiqV69OkyZN+Oabb6xtCxYsoEGDBnh4eHDbbbexfft2a1uHDh2YOHEizZo1\nw8vLK9syd+zYQYMG2Z/S0aNHM2bMGJsxFLWc9eTJk63lrF977bVsh1V+/vln2rdvj4+PD61ateK/\n//2vtc0Yw/vvv0+tWrUICAjggw8+sLYlJSXx1FNPERQURM2aNXnjjTesbatXr7Zuj/r16/P1119b\n4xg+fDgABw8epF27dnh6elKjRo1s5cHzknn46M0338TPz4+6deuyZs0aa/vRo0fp2rUr3t7e3Hrr\nrdkqwU6dOpXAwEC8vLxo1qwZp06dYvr06Wzbto3Bgwfj7u6ebQxZbdq0iVtvvRVvb2+6du1qLYuR\nX9vq1aupU6cOfn5+vP3221SrVs1693F7sYaGhjJ79mzCwsKoXr269SauObdhfHw8Xbt2xdfXl4CA\nAJ5++uls5TnssVfi3BjD6dOn+fe//829996bra1v37688847QP6vnXfffZfQ0FDrjWw/+OADgoKC\nCAkJYe7cuRhjrP3t/X/Z2x5paWlMnDiR0NBQvLy86NChQ7bnplmzZnh7e9OpU6d87/xeYkSkTB/A\nGCA54/e6gAD35egzEYjP8vdHwK4cfaoAV4F/FGS9zZs3l6IKXxRu8+Eo+/fvzz2tYaNSexREWlqa\nNGnSRN555x1JTU2VrVu3ip+fn8THx4uIyJo1ayQ2NlbS0tJk7ty5EhISIqmpqSIicvfdd0uDBg3k\njz/+kKSkJBERqV27tuzcuVPS09MlNDRUfvjhBxERSU9Pl5CQEPnxxx9zxfDVV19JSEiI/Pbbb3Lh\nwgXp1KmTjB8/Pt+21atXS2hoqBw5ckTOnj0rrVu3ltq1a4uIyKVLl6RGjRry2WefSVpamqxYsULq\n1KkjKSkpIiICyIABA+Ty5cuydetWcXFxkfPnz4uIyPDhw+WRRx6RxMREOXbsmDRo0ED+85//iIjI\nTTfdJN99952IiBw7dkx+//13ERGZNGmSDBs2TEREHnroIZkxY4akp6fLxYsXZe/evfk+D0eOHBFj\njAwdOlSuXLki69evFw8PDzl79qyIiLRp00ZefPFFSU5OltWrV4uXl5ecO3dOoqOjpWbNmnLmzBm5\ndu2a/Pzzz3LhwgXr87Ns2bI813n27Fnx8vKSNWvWSHJysowfP17atm2bb9uZM2fEw8NDNm7cKFev\nXpUnn3xSjDFy5MgRu7Fmvj46duwof/75pxw8eFBcXV0lOjo61zaMi4uTL7/8Uq5evSrHjx+X8PBw\nef/9963bqmrVqjbHtHr1aqlVq5YcPnxY4uLipHnz5tbXRObzfurUKbl48aK4ublZt29iYqK4u7vL\nqVOnCvTa6dmzp1y6dEmSkpJk37594u3tLT/99JMkJibKAw88IJa36Pz/v+xtj1deeUVat25t/f/b\nsmWLiIjExMSIn5+fbN++XdLS0uT111+XO+64I8/n2db7DrBHivC+Xx4uLmgL/Jbxe6OMnwdy9IkG\nqhtj/LP0y9ZHRFKAw1mWoRzo+++/Jz09naeffprKlSvTvn17Wrdubf0k361bN2rWrEmlSpV45pln\nuHz5crZPV0888QR16tTBxcUl23KNMQwYMMB6jmH79u24ublx++2354qhOOWsn3jiCUJDQ/Hz8+PZ\nZ5+1LnPNmjXcdtttPPjgg1SqVIl+/frh6urK3r17rX3Gjh2Lq6sr7du3x8/Pj0OHDiEiLF68mNmz\nZ+Pm5kZISAhDhw5l1apVADg7OxMdHU1SUhIhISHUq1cv13icnZ2JjY0lLi4OT09PmjVrVqDnQkSY\nNGkS1apV47777qNVq1asW7eO2NhY9u3bx0svvUSVKlX429/+RrNmzVi7di2VK1cmOTmZ6OhoRIRb\nbrkl195nXtatW0eLFi3o1q0bVapUYeLEifz3v/8lNjY237Y2bdrQqVMnqlatyosvvpj5IdJurJlG\njBiBj48PYWFhNG/enF9+yX13eH9/f7p3707VqlUJDg5myJAh2fa281LQEueenp506tSJTz/9FLC8\nXpo3b05gYGCBXjvjx4/Hw8MDFxcXoqKi6Nu3L7fddhtubm48//zz1n75/X/Z2x6LFi1ixowZ1v+/\nu+66C7CUKunbty933HEHlSpVYtSoUezbt4/Tp0/nu32Kq0wTT8ZFAw8CszMm+WT8vJCj6/kc7T42\n+mT287ExPXN9Q40xe4wxe86ePVu0oJVNx44d49ChQ3h7e1sfmzZtspaM/vzzz6279N7e3ly4cIH4\n+Hjr/FnLSec0cOBAVqxYgYiwbNkyHn441/UoQMmVsw4ODs42rm+++SbbuP744w9Onjxp7ZP1JLOr\nqyuXL1/m7NmzXLlyhbCwMOt8U6ZM4cyZMwCsXLmSlStXEhQUxAMPPGCzkunMmTO5fPky4eHhtG7d\nOtfhnrw4OTkRFBRk/Tuz/PfJkydzlSDP3A6ZhxjHjBlDYGAgzz77LMnJyTaXHxERgbu7O+7u7mzb\nti3Xtq1WrRo33XSTdZ15teXc7jVq1MDJyfKWZC/WTLa2e04JCQkMGjSI4OBgPD09GT9+fLbXXV4K\nU+L8oYce4pNPPgHgk08+4aGHHgIK9trJuo78Xof2/r/sbY/jx49Tp06dXHEfO3aMhQsXWpfn4+ND\namoqJ06cyHf7FFeZJR5jTCiwFPhCRBY5Yp0iMl9EWohIC39///xnUAUWHBxM06ZNuXDhgvVx+fJl\nRo4cSXJyMg8//DAzZ84kPj6eCxcu4OXlZf10C2Q7lp1T48aN8fPzY8uWLURFReWZeIpTzjrrP1vW\n34ODg+nWrVu2cSUlJdGrVy+728PPz49q1apx9OhR63wJCQl89tlnALRp04b169db32DHjRuXaxlB\nQUEsXryYM2fOMHToUAYMGGB3nZnS09Ozvblllv8OCgoiLi6OlJQUm9shMjKS77//nl9++YUdO3bw\n0UcfAbmfm/Xr15OYmEhiYiLt27fPtW2Tk5M5c+YMNWrUsNuWc7ufOnWK9PR069jtxVpQb7zxBn/+\n+Sf79u3j0qVLzJgxI9vrLi+FKXHeo0cPdu/eTUxMDF9//bW1OGBBXjtZt21+r8O8/r/yU7NmTZvn\nboKDg3nyySezLfPKlSs0b94832UWV5kkHmNMdWA9cBQYmKUpc88m5z6+T4728zb6ZPY7b2O6KmWt\nW7cmNTWVBQsWkJqaSnJyMlu2bOHkyZMkJyeTkpKCv78/IsJrr73GpUuXCrX8gQMH8swzzxASEkJY\nWJjNPsUtZ3306FHi4+OtJ4YBunfvzp49e/jyyy+5du0aSUlJrF271uan66ycnJwYNGgQo0eP5uLF\ni6Snp7N//3727t1LSkoKy5cvJzExkapVq+Lq6mr9pJ9VVFQUp0+fxsnJCU9Pz2x9jDF5ngg2xvDy\nyy+TnJzM119/zffff8/9999PrVq1aNq0KdOmTSM1NZV169axZ88e7r//fg4ePMjWrVtJSUnBzc2N\nKlWqWNeXX7nxiIgI9uzZw1dffUVqaiovv/wyt956K7Vq1cq3befOnWzatImUlBSmT59ufSO2F2th\nJCQk4OrqipeXF0eOHMl28Yc9hSlx7u7uTpcuXfj73/9Oq1atrHsehX3t9OzZk5UrV7Jv3z4uX77M\nrFmzrG32/r/yM3jwYF544QWOHz/OtWvX2Lp1KwADBgxgxYoV7NixAxHh0qVL1kOGpc3hiSfjuzZr\nsFwM0F1EkrI0Z563yXmephHwp4iczdIvWx9jTBUsFyfkPD+kHKBy5cqsWbOGtWvXUqNGDUJCQpg5\ncybp6el4enry+uuv06lTJ2rUqMHVq1ftHrqw5eGHHyY6Otrup/7ilLN+5JFHaN68Oa1bt6Zbt25U\nrVoVAC8vL9asWcObb76Jv78/devWZeHChQWKec6cObi5uREeHk716tV57LHHrAn3X//6FyEhIfj6\n+rJ//36mTZuWa/5du3Zx++234+7uziuvvMLixYsByyfhoKCgbIdisqpSpQoNGzYkODiYYcOG8fHH\nH5O5h79s2TK+++47/P39GTduHJ9++il+fn4kJyfz3HPP4evrS506dQgPD2fQoEEAPP3008yfPx9v\nb2+b34/x9/cnKiqK559/Hn9/f3bv3s3SpUvzbbvppptYvHgxjz32GEFBQTRo0IBKlSpZt31esRbG\niBEjOHbsGN7e3vTr14+ePXsWaL7Cljh/6KGH+Pbbb7OVbS/sa+fWW29l2rRpdO3alYYNG3LnnXda\nDzXa+//Kz5gxY+jQoQNt27alevXqTJ48GYA6deqwdOlSRo0ahY+PD40bN2bdunUF2DrF59DS18aY\nysAXQCugnYgcstHnILBNRIZk/O0E/BfYJyKPZEx7GPg3UE9EjmZM6wN8AjS0tdyctPR1xXL16lUC\nAgLYv3+/3fNBJeGjjz5i4cKF2S7fLU8++eQTDh8+zPjx43O1xcTE0KhRI65evVoGkRVPbGwsdevW\nJTk5mUqVKpV1OGVu69atREZGcuTIkbIOBSjZ0teO/sbSe1i+NDoC8DXG+GZp+0lEkrF8j+djY0wM\nsAOIBBoAWT/qrgJeAD41xryE5bDbHGBpQZKOqng+/PBD2rVrV2pJ54svvuD+++/n1KlTzJo1i8ce\ne6xU1lMSsn6qrug2bNjAHXfcAVjuONCtW7cbOumsXbuWzp07k5iYyJQpU+jRo0dZh1QqHJ14Mr9p\n9ZaNtjpAjIgsM8a4A2OBl4D/YTkk92tmRxFJNcbcB7yDZS8nGViO5TtB6jrTokULzpw5w5dffllq\n63jjjTd45JFHcHV1pW/fvjz11FOlti71l02bNtG/f39EhHbt2rFgwYKyDqlMLVu2jAEDBlC5cmW6\ndu3KlClTyjqkUuHQQ23liR5qU0qpgivJQ23l4QukSimlbiCaeJRSSjmUJh6llFIOpYlHKaWUQ2ni\nUUop5VCaeFS5pKWvVX62bdtW4PuKFeb1lLUmUEREBFFRUUUNUeXhxil5d70rzev9J00qvWWrbCZM\nmFDWIVQY7du3z1ZioDSsX7++VJd/o9I9HlWupKWllXUINwzd1o6n29xCE48qMVr6uvyXvg4NDeX1\n11+nYcOGNG3aNN/x5WSvHLW9tt9++41WrVrh4eHB4MGD6dKlC4sWLQLgypUrDBs2jICAAOrUqcNb\nb/11Y5PBgwfzf//3f9x99914enrSs2dPkpIs9xXevHkzjRr9da/gp556iho1auDt7U1ERESB68qs\nXr2aunXrWstvZ9WhQweWL19ObGwsXl5e2e6B9+6779KnTx8Azp07R//+/fH396devXosW7Ys2zJy\nvr537dpFkyZN8Pb2Zty4cTRq1IjNmzcD9l839rYHWAoehoeH4+HhQdOmTTl48CBg/3+zLGjiUSXi\n2rVrdO/ene7duxMXF8e8efPo378/f/75J2Cpr/Ltt99y4cIFhgwZwsMPP5zt09/y5cuJiori3Llz\n2Zbbrl070tLSyLzLhIiwYsUKm3ep3rBhA7NmzeLrr7/m0KFDHD58mJdffjnfti+//JLFixfzww8/\nsGvXrmy3hk9ISCAiIoJRo0Zx7tw5Ro8eTa9evUhNTbX22b59OwcOHCAqKoqRI0dy4YKlRuGoUaNI\nSEjg0KFD7Ny5k3nz5rFx40YAhg4dykcffURCQgKbN2+2WYF00qRJ9OjRg4sXL3Lw4EHuvffeXH1s\nSUlJITo6muPHj/Pee+8xYMCAbNv1s88+Y/v27fz4448FGl9W/fv3p1WrVsTFxfHKK6/Qs2dPa2E1\ne20DBgygR48e/Pnnn3Ts2NH6JgswdepUjhw5wu+//86GDRt47bXXslXWXLlyJe+//z4nTpwgJiaG\nf//73zZju/POOzl48CCnTp2ievXqjBo1Kt9tFRcXx8CBA5k3bx7Hjx9n37592WoAZapVqxaNGzfm\nq6++sk775JNPrPfNe/TRRwkNDeXEiROsXr2akSNHWt/0IfvrOzk5md69ezN27FjOnj2Lu7s7hw8f\ntva197qxtz127NjBP/7xDxYsWMClS5dYuXIlnp6e+f5vlgVNPKpEaOnr8l/6OtPIkSPx9/fHxcWl\nQOPLZK8ctb22mJgYDhw4wNixY3F2dmbQoEHZKmKuWLGCiRMn4unpSVhYGMOGDbM+N2ApOdC4cWM8\nPDz429/+xr59+2yOe8CAAXh6euLi4sLzzz9foBLX69ato23bttx7771Uq1aNyZMn51korl+/ftZK\no6dOneKnn36ie/funDp1ih07djBt2jSqVKlCkyZN6NOnj7XoH2R/fe/cuRM3NzcGDRqEs7MzY8eO\ntZaCyO91Y297LFq0iCeffJK2bdtijKFRo0bUqFGjQGWzHU0TjyoRWvraojyXvs6Udaz2xjd9+nRr\nievp06fbLUdtr+306dMEBATg7Oxscxvbe27y2r62vPLKK9SvXx9PT0/at29fpBLXwcHBNovygaWY\n4Nq1a7l69SqrVq0iIiICV1dXjh07xuXLl/H19bVuw0WLFnH69GnrvPZKXDs7O1vrJeX3urG3PeyV\nuM6vbLaj6VVtqkRklub98ccfc7Vllr5evXo1HTt2pFKlSvj4+BS59PW2bdts9ivN0teFrcyYtfR1\nzvNW8Ffp66SkJEaMGMG4ceNYsmRJrvEsXryY9PR0Fi5cyIABAwp03iKz9HVm8jl27BgdOnSwtmfd\n1vmNL+tVdpl7XykpKdYEExsbyz333JOtVHXOtsDAQOLi4khNTbUmn6yJO/O5qV27tnW+wpa43rJl\nC/Pnz+fbb7+lXr167Nu3j9atW+c7X2BgIN9++6317xMnTuRZXC0kJISmTZuyfv16Vq5cyYgRIwDL\nNvT29rab6OyVuE5NTeXsWUuNy/xeN/bYK3Gd1/9mWdE9HlUitPR1duWx9LUthRmfvXLU9tpCQ0MJ\nCwvj9ddfJzU1lY8//pg//vgj2/P2yiuvcOnSJQ4dOsQHH3xQ6JpDCQkJODs74+vrS2JiItOnTy/Q\nfJnltzdu3EhycjJTp061+yGoX79+vPXWW/z888/WbRocHEzLli2ZNGkSV65csZ6TzHqOJ6u2bduS\nkJDAxx9/TGpqKq+//jrJycmA/ddNfiIjI3n//ffZvXs3ImI931WcstmlRkRuyEfz5s2lqMIXhdt8\nOMr+/fsdtq7CiImJkQceeEB8fX3Fz89PIiIi5NixYyIiMnv2bKlevbr4+fnJlClTpHbt2rJz504R\nEbn77rtl2bJl2ZaVtV1E5Pjx4+Lk5CSzZs2yG8Mbb7whtWrVEn9/fxkyZIhcvny5QG0vvfSS+Pr6\nSr169eSVV16RsLAwa9vPP/8sHTt2FB8fH7npppukd+/ekpiYKCIigJw6dcrat2HDhrJp0yYREbl8\n+bKMGDFCQkJCxMvLS1q1aiXffvutJCcnS5cuXcTLy0u8vLykc+fOcvToURERmTRpkgwbNkxEREaP\nHi2BgYHi5uYmt9xyi3zzzTfWbREUFCQpKSm5xn/kyBGpWrWqzJkzR3x9fSU0NFS++OKLPLdrfuPL\n6Y8//pDOnTuLl5eXNG3aVDZu3FigtujoaGnRooW4u7tLZGSk3H333bJ06VLrdhoyZIj4+/tLrVq1\n5I033rDOFxkZKTNmzLD+PWPGDImMjBQRkU2bNknDhg1FRCQ1NVUGDBgg7u7uUqdOHXn77belatWq\ndsed6bPPPpPQ0FDx9fWVN998U6pWrSpHjhwRkdyvzZMnT4qTk5P0798/2zLOnTsngwYNkptuukl8\nfHzk7rvvlp9//tnmMkREtm/fLjfffLN4enrK2LFjpV69evLdd99Zt4et101+20NEZMmSJdKoUSNx\nd3eXpk2byoEDB0TE/v9mQdl63wH2SBHef7UeTxFoPR7H09LXf7keSl/XrVuXhQsXcvfdd5d1KGUu\nOTkZHx8foqOjrYcbyyOtx6NuOI4ofZ2amkpsbCyzZs0q1yWH+/XrZzPplGe7d+/m+PHjpKWl8c47\n75CUlFSgczDXq82bNxMfH09ycjITJ04kLCysXCedkqYXF6hyT0tfV3zHjh2jZ8+eXLp0ibCwMFat\nWkW1atXKOqwy88svv9C3b1+Sk5O59dZbWbp0aVmH5FB6qK0I9FCbUupGo4falFJKVViaeJRSSjmU\nJh6llFIOpYlHKaWUQ2niUeVSUSqQLlq0iPvuuy/PdmNMtvtnlYT77rvPent/pVTBaOJRSinlUPo9\nnutE6Li1pbbsmFe7ldqyc0pLS6NyZX1ZKnU90z0eVWLKQwXStLQ0Bg0ahIeHB23atOH333+32a84\nlTRbt26Nh4cHkZGRXLt2reAbSCkFaOJRJaQ8VCAFy61I7rnnHuLj4+natSsDBw602a84lTS7detG\nfHw87du3L7f3c1OqPNPEo0pEeahACpabTz722GNUqVKFCRMm8OuvvxIbG5utT3EraY4bN44qVaow\nZMgQQkNDS24jKnWD0MSjSkR5qECaczlVq1YlICAgV6XF4lbSzNqWtXKmUqpg9CyuKhHloQIpWMr/\nZkpJSSEuLo7AwMBsfexVyyxsJU1bJauVUvY5fI/HGFPfGPOBMWafMeaaMWazjT4xxhjJ8cj1BQxj\nTGNjzDfGmCRjzEljzFRjTCWHDERlUx4qkAL88ccfLF68mNTUVGbMmEGTJk1y3W6+uJU0Z86cSWpq\nKgsXLszUI2egAAAgAElEQVSzCqhSKm9lcaitCXA/cBD4zU6/pUDbLI9stXuNMT7ARkCAB4CpwChg\nSsmHrPJTuXJl1qxZw9q1a6lRowYhISHMnDmT9PR0PD09ef311+nUqRM1atTg6tWrhT5E9fDDDxMd\nHZ3nRQWZOnTowMaNG6levTpr167l448/ttlv2bJlfPfdd/j7+zNu3Dg+/fRT/Pz88m1bunQpq1ev\npnr16tYLGZRShePwsgjGGCcRSc/4fRXgJyIdcvSJAVaJyGg7yxkPPA/UFpFLGdOeByYDgZnT8qJl\nESoWR1YgVUrlVpJlERx+jicz6ZSACGBDjgSzHJgJ3A2UWtWw4Tvfst0QWVprVKVdgVQp5Tjl+eKC\nvxtjngWuAP8BRonI0SztjYBvs84gIrHGmKSMttIrV6kcyhEVSJVSjlNeE88XwC7gOHAzMAnYZoxp\nKiIXM/r4ABdszHs+o01dJ4p6SFQpVT6Vy8QjIiOy/LnNGPMd8F/gMeDNoi7XGDMUGAqWK5uUUko5\nXoX4AqmI/IrlKrhmWSafB7xsdPfJaLO1nPki0kJEWvj7+5d8oEoppfJVIRJPBsl4ZDqA5VyOlTGm\nJuCa0aaUUqocqhCJxxgTjiXJ7M0yeT3Q1RjjkWXaQ1guRtjiwPCUUkoVgsPP8RhjXPnry6DBgKcx\npk/G3+uAe4BHgDXASSwJ50UgFliUZVHzgGeBT40xM4G6WL7D80Z+3+FRSilVdspijycAWJnxaAM0\nzvJ3AHAs4+ebwNdYrmj7D3Bn1oQiIueBTkAlLJdOTwHmZPRXFVxRSl/bUxplrydPnszw4cNLdJlK\n3QjK4gukMUDed4S06FTAZe0HOhY3JqWUUo5T4MRjjKkFnBKRVBttlYEgEYnNPadyiMm2LvArqWVf\nzL9PCdHS10pd/wpzqO0IYLv6Ftya0a5uYOWh9PWKFSsIDw/Hw8ODpk2bcvDgQWvbZ599Rq1atQgI\nCOCDDz6wTk9KSuKpp54iKCiImjVr8sYbb1jb0tLSeOmll6hduzaBgYGMHTvWZrnrixcv0rZtW6ZO\nnVrAraXUjaswicfe4bFqQHIxY1EVWHkofb1jxw7+8Y9/sGDBAi5dusTKlSvx9PS0tm/fvp0DBw4Q\nFRXFyJEjuXDBcuOLUaNGkZCQwKFDh9i5cyfz5s1j48aNAMyePZtdu3bx448/Eh0dzfbt21m4cGG2\n9Z4/f54uXbrQo0cPJk6cWAJbU6nrm93EY4y5xRgzyBgzKGPS/Zl/Z3kMxXJi316JA3WdKw+lrxct\nWsSTTz5J27ZtMcbQqFEjatSoYW0fO3Ysrq6utG/fHj8/Pw4dOoSIsHjxYmbPno2bmxshISEMHTqU\nVatWAfCvf/2LadOm4evri4+PDyNHjrS2AcTHx9O5c2f69evH+PHjS2x7KnU9y+9gek/+ukpMgLw+\nzh0BhpVUUKriyVr6OlNqaiqdOlmuE/n888+ZOnUqf/zxBwCXLl0iPj6e+vXrA/mXvu7atSuvv/66\n3dLXx48f54477shzOQEBAdbfXV1duXz5MmfPnuXKlSvZistdu3aNzp07W8fVpUsXa4XU9PR0a8wA\nW7ZswcXFhaFDh+a5XqVUdvklnunALCyH2S5huYLshxx9UmxdcKBuLOWh9HXNmjULXRHUz8+PatWq\ncfTo0VznlzLHFRUVxS233GJz/l69euHu7k737t3ZsGFDrj02pVRudg+1iUiqiFwWkUQRcRKRzRl/\nZ31o0lHlovR1ZGQk77//Prt370ZEOHjwIKdOnbK7XCcnJwYNGsTo0aO5ePEi6enp7N+/n717LTfJ\nePzxx3nhhRc4c+YMIsIff/yRK/HNmjWLBg0a0Lt3b1JT9d9BqfwU+gukxpgwY0xHY8z9OR+lEaCq\nGMpD6es77riDOXPmMHjwYDw9Penbt2+BEtycOXNwc3MjPDyc6tWr89hjj1nnGzNmDC1btqRt27Z4\neXnRs2dPm8ls/vz5uLq68sgjj5CeXlK1DpW6PhW49LUxpjGWCp9NsH2Fm4hIpRKMrVQVp/T1u8O/\ntTn96XmO+S6rlr7WKqRKOVpZlb7+AKgK9AL2AymFXZlSRaWlr5W6fhQm8dwO9BeRNaUVjFK2aOlr\npa4vhUk8h7F8UVQph9LS10pdXwpzccEoYIIxpm5pBaOUUur6V5g9nhlY6uccMMbEABdydhCRViUU\nl1JKqetUYRLPrxkPpZRSqsgKnHhE5LHSDEQppdSNoSwqkCqllLqBFTjxGGM+ye9RmoGq8q+ky1Vf\nr5YsWULPnj3LOgwAhg8fnq3+UGFEREQQFRVVwhGpG0FhzvH425jmAzQC4oGDNtqVymXy5MmcPn2a\nefPmlXUoZWLgwIEMHDiwrMMAKPBzsGjRIpYvX85XX31lnbZ+/frSCqvEdejQgeHDh9O/f/+yDkVR\nuHM899iaboypCXwGzCmpoFTh5XUbn5LgqFsBKZXVtWvXqFTpr7twaVn060exz/GIyDEsl1q/Vvxw\n1PXiiy++sJagbtiwIZ9//jkA27ZtY/r06Xz44Ye4u7tbDznZK5ud0yeffEJYWBi+vr7069fPWuU0\nv7YPPviAoKAgQkJCmDt3brZSDPbWHxoayuzZswkLC6N69erZqozaK4393Xffceutt+Lh4UGtWrVY\nvHgxYNl7uO+++wA4c+YM9957L15eXvj7+/Pss88WaPsaY3j33XepWbMm/v7+rFy5kjVr1lC3bl0C\nAgJYsWKFte+5c+fo378//v7+1KtXj2XLllnbBg8ezKuvvgpY3uhfeOEFgoODCQ4O5oUXXuDatWvE\nxsYyfPhwNm7ciLu7O82bNwcsexGZBfrS0tKYOHEioaGheHl50aFDB5tx79y5kxYtWuDp6UloaGi2\nEuSTJ09m4MCBPPjgg7i5ubFz506bZdHtPVdTp04lMDAQLy8vmjVrxqlTp5g+fTrbtm1j8ODBuLu7\nF/nQoio5JXVxwTVAb6KlrLy8vIiKiuLixYu89tprDBo0iPj4eNq3b8+ECRP4+9//TmJiIp999lm+\nZbOz2r9/P8OGDePjjz/m+PHjVK1alWeeeSbftl9++YVx48axbt06Dh48mO3NqiDrX7duHbt372bX\nrl3Mnj2bAwcOAPZLYz/33HOMGzeOhIQEfvzxR1q0yH0vxTlz5tCgQQPi4+OJjY0t1CG4nTt3cujQ\nId577z2GDx/OJ598wq+//srixYt55plnrAnw0UcfJTQ0lBMnTrB69WpGjhzJwYO5j4z/85//ZO3a\ntezdu5c9e/bw5Zdf8uGHH1KrVi3mzZtH586dSUxMtJaMyGrmzJl8/fXXbNu2jT///JOpU6fajNnZ\n2Zn58+dz4cIFli9fzujRo9m/f7+1/dNPP+Uf//gHCQkJtGzZEsheFt3ec3XgwAH++c9/sm/fPs6f\nP8+iRYtwdXVlwoQJtG/fnkWLFpGYmMhzzz1X4G2sSkdhLi5obONxmzHmYSzF4nIWiFM3sA4dOtCw\nYUOcnJx44IEHqFevHj/99JPNvvmVzc5q1apV9O7dm1atWuHi4sK0adNYtWoV165ds9sWFRVF3759\nue2223Bzc+P5558v1PpHjBiBj48PYWFhNG/enF9++QWwXxrb2dmZw4cPc+HCBfz8/GjSpEmu8Tg7\nO3Pq1ClOnDiBi4sLrVu3LvA2fv7556lWrRq9evXiwoULPPPMM7i6uhIREUFSUhInT57k1KlT7Nix\ng2nTplGlShWaNGlCnz59+Oyzz3Itb8WKFYwePZrAwEBq1KjB6NGjrXs0+Vm0aBEzZsywlje/6667\nbPZr0aIFzZo1w8nJiTZt2tC5c2e+++47a3vHjh3p1KkTTk5OVK1aFcheFt3ec1W5cmWSk5OJjo5G\nRLjllltsFvdTZa8wezy/Ar/keOwFlgB/AkNKPDpVYW3fvp127dpRvXp1vL29+eWXX4iPj7fZN2vZ\n7MzHpk2bbNa9OXnyZLZaPjVr1uTatWucO3fObtvp06ez3dk6ODi4UOu3VTY7c94uXbpY53v88ceJ\ni4sDYMGCBfz000/UqVOHjh07WpNVVmPGjCE4OJh27doRHh5eqKvEMmOqVKkSzs7O2WJ0cXHh8uXL\nHDt2jMuXL+Pr62uNcdGiRZw+fTrfbVu7du18C+llOn78OHXq1Mm33//+9z+6dOmCv78/Xl5erF27\nNtvrwtbdx7NOs/dc1a9fn9dee40xY8YQGBjIs88+S3JycoHiV45VmDN1ti4uuAocF5ETJRSPuk48\n+uijvPTSSzz66KM4Oztz++23W0td5yxzba9sdk5BQUHExsZa/z5+/DiVKlXCz8/PbltgYCAnTvz1\nMs36e2HWn5O90tiNGjUiKiqK1NRUpk2bxvDhw9mxY0e2Pp6ensydO5e5c+fyn//8h+7duxMfH4+7\nu3uhY8krPm9v7zyTflZBQUEcO3bM+ndsbCw1atQA7Jcmh7/KjoeGhtrt98wzz9ChQwdWr16Ni4sL\nPXv2zLcEetZp+T1XkZGRREZGcvr0abp168ZHH33EE088kW/8yrEKvMcjIltsPHZr0lG2JCQk4Ovr\nS6VKlVi2bBm//vrX3ZYCAgI4evSo9W97ZbNz6t27N1FRUezZs4erV6/y0ksv0atXLypVqmS3rWfP\nnqxcuZJ9+/Zx+fJlZs2aVaT152SvNPbSpUs5f/48lStXxt3dHSen3P9u69atIyYmBrCcFzPGWK/k\nCg0NZfPmzQXa3nkJDg6mZcuWTJo0iStXrpCWlsaePXtsnuPp27cvs2fP5syZM5w+fZpZs2bRr18/\nwPKcHT9+3HreKKfBgwfzwgsvWPts3brVZr+EhAR8fHyoVq0amzZtsnk41R57z9XBgwfZunUrKSkp\nuLm5UaVKFes2DwgIsG5nVfYKdXGBMaayMeYhY8xcY8ySjJ/9jDE31DWOTwf2tPlQf5k7dy5PPvkk\n1atXZ+vWrbRt29ba1qtXL/788098fHzo06eP3bLZOTVp0oT33nuPhx9+mKCgIBITE3nnnXfybbv1\n1luZNm0aXbt2pWHDhtx5551UqVIFsF+2Oz/2SmOvXbuWBg0aWC+0yIwlqwMHDnDXXXfh7u7OwIED\nWbx4MS4uLqSmppKYmEh4eHjhN34OS5YsISYmhjp16hAQEMDo0aNtHoIaOnQo9913H7fffjvNmjXj\n/vvv54knngDgnnvuISgoCD8/P1q1yn0v4DFjxtChQwfatm1L9erVmTx5ss1YZs6cyVtvvYWnpyfv\nvfce3bp1K9RY7D1XycnJPPfcc/j6+lKnTh3Cw8MZNGgQAE8//TTz58/H29ubN998s1DrVCWvMKWv\nA4CvgVuAGOAMcBMQCvwM3CsiZ0slylJQnNLXTM7jhOXki0UPqBBuxNLXJW3r1q1ERkZy5MiRsg7F\npu+//54PP/ww2+XGSpWlkix9XZg9njcAX6CNiNQVkbYiUhdonTFdL45X5dratWtJTk4mPj6eKVOm\n0KNHj7IOKU+tWrXSpKOuW4VJPPcDY0Xk+6wTReQHYDxQuH1mpRxs2bJlBAQEEBYWxk033cSUKVPK\nOiSlbkiFOTdTFUjIoy0BqFL8cJQqPR9//HFZh6CUonB7PLuAscYYt6wTM/4em9GeL2NMfWPMB8aY\nfcaYa8aYzTb6GGPMBGPMMWPMFWPMVmPMbTb6NTbGfGOMSTLGnDTGTDXGVMrZTymlVPlRmD2eUcBm\n4Jgx5mssFxcEAF0BA3Qo4HKaYDlstwtwzqPPOOAlYAxwAHgO2GiMCReR0wDGGB9gI7AfeACoB8zG\nkkxfLMS4lFJKOVBhvsfzX6A+MB9LiYQuWBLPPKCBiPxcwEV9KSI1RaQv8L+cjcaYalgSzwwReUdE\nNgJ9AQGeydJ1OOAC9BKR/4jIPGAK8JwxxrOg41JKKeVYhblX261AKxEZJyKdRKRxxs8JQCtjTO6v\nbtsgIvl9OaId4AlYC8uJyGXgSyAiS78IYIOIXMoybTmWZHR3QWJRSinleIU5xzMHy6XTtrSk5Orx\nNMJyt+tDOaZHZ7Rl7XcgawcRiQWScvRTSilVjhQm8TQDduTRthO4vfjhAJaqpokikvPeHOcBV2NM\nlSz9LtiY/3xGm1KlKmstm8IqTsnpiuJGGKMqmsJcXFAJcMujzY0KcDm1MWYoMBSgVq1aZRxNyWq6\nuGmpLfuXyNx3Va7oNm/ezPDhw611dRy9ruut7PeNMEZVcgqzx/MDGW/aNgwFinj/mVzOA+42Lov2\nAZJEJCVLP1v3rvHJaMtFROaLSAsRaeHv719C4arCSktLK+sQ8lURYlSqoipM4pkMdDLG7DbGPGWM\n6WWMedoYsxtLyYSXSiimA1j2rurnmJ7znM4BcpzLMcbUBFxz9FMOYozhrbfespZO/uc//2ltK2wJ\n45zOnDlDly5d8PT0pHv37jzyyCPWG1EeOnSI9u3b4+3tTXBwMJMmTbLOd/DgQdq1a4enpyc1atSw\nHhqLiIjgt99+w93dHT8/vzxjNMZkq13TqFEjm3eMtheDrXUVpOQ0WMpBP/roozz44IN4eHjQoUMH\na72fnFavXk2DBg3w8PCgfv361js/X7lyhWHDhhEQEECdOnV46623rPMMHjyYESNG0L59e9zd3Xn0\n0Uc5ffo0nTt3xsvLiwEDBlhjKc4YAd5++21rfK1atSpQqQZ1fSrM5dRbgXuBdGAusAp4C0gDuojI\nthKK6TvgEpZLqAEwxrgCfwPWZ+m3HuhqjPHIMu0h4AqwpYRiUYW0ceNGDh48yOrVqxk9ejT/+99f\nV8wXtISxLU899RR169bl7Nmz2ap8Zpo2bRrnzp3jm2++YcGCBXz11VcATJo0iR49enDx4kUOHjzI\nvffeC8D69esJCwsjMTGRc+fO2YyxsPKKIa91Zcqr5HSmzz//nDFjxhAfH4+rqytz5ti+jmfo0KF8\n9NFHJCQksHnzZurVqwfA1KlTOXLkCL///jsbNmzgtddey1aOYNWqVSxYsICjR4+ydetWunXrxptv\nvklsbCx79+5lzZo1xR7jkiVLmDt3Ll988QUXL15k3rx51ruDqxtPocoiiMhmEWkLeAA1AU8RuaMw\nSccY42qM6WOM6QMEA/6ZfxtjXEXkKvAqMCFjj6oTsDIj1rlZFjUPSAY+NcZ0zjh/Mxl4I8cl1sqB\nxo8fj7u7O82bN6dPnz7ZKmoWtIRxTmlpaXzxxRdMmjSJqlWr0qlTp2yllRs0aMBdd91F5cqVadSo\nEQ899BDbt28HLKWlY2NjiYuLw9PTk2bNmtmNP2uMhWEvhvzkV3K6S5cu3HHHHVSpUoU+ffqwb98+\nm8txdnYmOjqapKQkQkJCrIlnxYoVTJw4EU9PT8LCwhg2bFi25T/00EM0atQIX19fOnToQOvWrQkP\nD8fLy4vOnTtb11ecMS5atIgJEybQuHFjnJycaNasGR4eHvnPqK5LhUo8mUQkSUROiEhSEWYPwJJI\nVgJtgMZZ/s6s3fsqMA3LzUfXYPleTxcROZMlhvNAJyyH5b7E8uXROcBf+//K4bKWKa5Zs2a20skF\nLWG8ZMkS3N3dcXd3Z/jw4Zw7dw4RISgoyDp/1tLVJ06c4MEHHyQwMBAvLy/ef/9962GcmTNncvny\nZcLDw2ndujVbttjfGbZVerkg7MWQn/xKTudVdjunlStXsnLlSoKCgnjggQes1UQLs3wXFxebJbSL\nO8aClsZWNwaHF3ATkRgst9ix10ewJJ5p+fTbD3QsseBUsR0/ftx6xeCxY8eyveEVpoTxwIEDrb+n\npaVhjOHUqVPWUswnTpygdu3aALz44osEBATw22+/4enpyciRI7ly5QpgKee8ePFi0tPTWbhwIQMG\nDODEiRN5lkLOOd3V1dW6LBHh7FnbJafsxZBf2WV7JacLo02bNqxfv56kpCRGjBjBuHHjWLJkiXX5\nmdurqMsvzhgzS2MrBUXc41EqLzNnziQxMZGffvqJVatW0bt3b5v9ClNuunLlyvTo0YOpU6eSkpLC\npk2brOWlwVJO2cPDAw8PD37++WeWLVtmbYuKiuL06dM4OTnh6emZrRRyXFyc9Y0zL02bNmX58uVc\nu3aNt956i4sXbRf7sxdDfuuyV3K6oFJSUli+fDmJiYlUrVoVV1dX61j79u3LK6+8wqVLlzh06BAf\nfPBBoZdf3DEOHjyYGTNmcODAAUSEH3/8kYSEvG52r653N1TJ6utZefmuTceOHWnYsCEiwsyZM/Ms\n3ZxZwnjEiBGMHz8eYwwtW7Zk/vz5Nvu/++67PPLII/j6+nLXXXfRo0cPqlatCsDEiRMZOHAgHh4e\ntGnThl69elnLVu/atYtnnnmGhIQE6tWrx+LFiwG4+eabiYiIIDg4mCpVqmS7ci2rOXPmMGjQIF57\n7TWeffZZ6tfPebEl+caQ37qGDh1KbGwst99u+Q52ZGSkteR0YfzrX/9i+PDhALRs2dJ6gcKkSZMY\nMWIE9evXx8XFhVGjRnHfffcVevnFGeOAAQOIi4vj/vvv5+zZszRp0oR169YVOgZ1fShw6evrjZa+\nLnmZh8MCAwNLfV0dO3YkMjKSyMjIUl+XUqrsSl8rVWb+97//ceDAAdLT0/n888/ZtWuX9dJopVTF\noofaVIVw4cIFBg4cSFxcHLVr12bJkiVFOkGulCp7mnhUiSnNw7Z33HGHXhWl1HVCD7UppZRyKE08\nFVTm1URKKVXaSvr9RhNPBeTm5saJEydISUkp1cNbSqkbm4iQkpLCiRMncHPLqypO4ek5ngooJCSE\nc+fOcfToUb19v1KqVFWuXBkvLy/rXcdLZJkltqQbSOjVpTanxzho/U5OTgQEBGS7p5ZSSlUUeqhN\nKaWUQ2niUUop5VCaeJRSSjmUJh6llFIOpRcXFEGMS143F+3m0DiUUqoi0j0epZRSDqWJRymllENp\n4lFKKeVQmniUUko5lCYepZRSDqWJRymllENp4lFKKeVQmniUUko5lCYepZRSDqWJRymllENp4lFK\nKeVQmniUUko5lCYepZRSDqWJRymllENp4lFKKeVQmniUUko5VLlMPMaYwcYYsfEYnqWPMcZMMMYc\nM8ZcMcZsNcbcVpZxK6WUyl95r0DaEbiS5e8/svw+DngJGAMcAJ4DNhpjwkXktONCVEopVRjlPfH8\nICKJOScaY6phSTwzROSdjGk7gRjgGeBFRwaplFKq4MrlobYCaAd4Ap9kThCRy8CXQERZBaWUUip/\n5T3xHDbGpBljDhpjhmWZ3gi4BhzK0T86o00ppVQ5VV4PtZ3Ccv7me6AS0B+YZ4xxFZE5gA+QKCLX\ncsx3HnA1xlQRkZScCzXGDAWGAtSqVas041dKKZWHcpl4RGQDsCHLpPUZ53VeNMa8VYzlzgfmA7Ro\n0UKKF6VSSqmiKO+H2rJaBVQHQrHs2bgbYyrl6OMDJNna21FKKVU+VKTEI1l+HsByCK5+jj6NMtqU\nUkqVUxUp8fQBzgFHge+AS0DfzEZjjCvwN2B9mUSnlFKqQMrlOR5jTBSWCwv2YdmzeSjj8ayIpANX\njTGvAi8ZY87z1xdInYC5pR1f9LLlNqffPGlSaa9aKaUqvHKZeICDwONATcAA+4FBIvLvLH1exZJo\nxgO+wB6gi4iccXCsSimlCqFcJh4RmQBMyKePANMyHkoppSqIcpl4yruIB2fZnB7j2DCUUqpCqkgX\nFyillLoOaOJRSinlUJp4lFJKOZQmHqWUUg6liUcppZRDaeJRSinlUJp4lFJKOZQmHqWUUg6liUcp\npZRD6Z0LlFLqOvLu8G9tTn96XkcHR5I33eNRSinlUJp4lFJKOZQmHqWUUg6liUcppZRD6cUFSilV\nweR1AUFFoXs8SimlHEoTj1JKKYfSxKOUUsqhNPEopZRyKE08SimlHEoTj1JKKYfSxKOUUsqhNPEo\npZRyKE08SimlHEoTj1JKKYfSxKOUUsqh9F5tSilVTlX0e7LlRfd4lFJKOZQmHqWUUg6liUcppZRD\nVehzPMaYxsBcoC1wAfgnMEVErpVpYEopVc7kdb7o6XkdHRxJBU48xhgfYCOwH3gAqAfMxrIX92KZ\nBDXZy07bRcfFoZQql67XiwUKq8ImHmA44AL0EpFLwH+MMZ7AZGPMaxnTHCr06tI822IKuazy9OlE\nqYpI/4fKr4p8jicC2JAjwSzHkozuLpuQlFJK5aci7/E0ArJ9pBGRWGNMUkbbl2USVTlTHj/1lWRM\nhT10oZ92syuPr4/ySF9nJasiJx4fLBcU5HQ+o83h1n8+Os+2dy+42JxeUi9QRxw71jepv9jb3qW9\nPa7350Hf5K9/RkTKOoYiMcakAmNE5M0c048DH4nIBBvzDAWGZvzZEDhYxNX7AeeKOO/14EYe/408\ndtDx38jjtzX22iLiX9gFVeQ9nvOArcvIfDLachGR+cD84q7YGLNHRFoUdzkV1Y08/ht57KDjv5HH\nX5Jjr8gXFxzAci7HyhhTE3DNaFNKKVUOVeTEsx7oaozxyDLtIeAKsKVsQlJKKZWfipx45gHJwKfG\nmM4Z528mA2844Ds8xT5cV8HdyOO/kccOOv4befwlNvYKe3EBWG+Z8w7Zb5kzWW+Zo5RS5VeFTjxK\nKaUqnop8qK3EGWMaG2O+McYkGWNOGmOmGmMqFWA+L2PMQmPMeWPMRWPMEmOMryNiLklFGb8xpmXG\n2H/PmO+gMWaSMaaao+IuKUV9/rPM72SM2WOMEWNM99KMtaQVZ+zGmF7GmB+MMVeMMfHGmK+MMW6l\nHXNJKsb/fgtjzNfGmD8zHhuNMa0dEXNJMcbUN8Z8YIzZZ4y5ZozZXMD5ivy+V5Evpy5Rxbzp6CdA\nGDAESAdmAp8D7Usr3pJWjPE/lNF3JnAIuAV4OeNn71IMuUSV0E1nhwAhpRJgKSrO2I0xQ7Ac7n4N\nGIPl6wwdqUDvLUUdf8ZVtBuBH4FHMyaPwXLfyKYicrQ04y5BTYD7gV2AcyHmK/r7nojow3K4cTyW\n70H9KYIAAA/9SURBVP94Zpn2PJCUdZqN+doCAtyVZVqrjGmdy3pcDhi/n41pQzPGX7usx1Xa48/S\n1wc4C/w9Y+zdy3pMjnjugQTgibIeQxmNfzhwDfDK8Tq4BjxZ1uMqxPidsvy+CthcgHmK9b6nh9r+\nUtSbjkYAZ0Rka+YEEfkeOJLRVlEUafwiYutb3D9l/AwqufBKXXFvOvsysAP4phRiK21FHXu/jJ+L\nSyswBynq+J2BNOBylmmJGdNMSQdZWkQkvQizFet9TxPPXxqR44unIhKL5VNPI5tz5DFfhuh85itv\nijp+W9pi2fU+XDKhOUSRx2+MuQV4HMj7Zn3lW1HH3hrLbaf+bow5boxJNcbsNsa0K71QS0VRxx+V\n0We2MSbAGBMAzMGy97SylGItL4r1vqeJ5y9FveloubtZaRGVyDiMMYFYjov/W0TiSig2RyjO+OcC\n74jI7yUelWMUdeyBWO55+CIwFvgblk//XxljbirpIEtRkcYvIieBe7CcyzyT8egFdBWRs6UQZ3lS\nrPcLTTyqxBhjqmA54ZgIjCzjcBzCGNMfy5vvK2UdSxkwgDvwdxFZIiJfAQ9iOcfxTJlG5gDGmBpY\n9mz2Yjm8FJHx+1pjTK2yjK2808Tzl0LfdLSY85U3xRqHMcYAH5FxhYyIVKSxQxHGb4xxBl7HcjWP\nkzHGG/DMaHbLcTun8qw4r30BNmdOyDhPshdoXILxlbaijn8MlvM8fUTkq4zE2xtL4q2oh10Lqljv\nF5p4/lLUm47mmi9DXsdAy6vi3nT1TSyXoj4gIhVp3JmKMn43LJdPv4Hln+088HNG23L+usiivCvq\ncx+NZa8n54l0g+UcX0VR1PE3Av4nIqmZE0QkBfgflkuyr2fFet/TxPOXot50dD0QaIy5M3OCMaYF\nUDejraIo8k1XjTHjsRxaeUREtpdeiKWqKONPxHKMP+vj4Yy2CcDA0gm1xBX1uV+T8fOezAnGGC+g\nOX8l4IqgqOM/CoRnHGIGwBhTFQgHYkohzvKkeO97ZX0NeXl5YNlFPAX8B+iM5bsoicArOfr9DnyY\nY9oG4A8sJxYfxHKlz7ayHpMjxg8MwHK4ZSHQJsfDv6zH5YjnP0d7KBXvezzFee1/njFvJNANyxv1\nWcCnrMdV2uPHkmBTgbUZY++e8aabCtxa1uMqxPhdgT4Zj51Y9tgy/3a189wX+X2vzAddnh5Yjkt/\ni+WTziks382olKNPDLAoxzTvjDfeC8AlYCk2vlhZ3h9FGT+wKOON1tZjcFmPyRHPf472Cpd4ijN2\nLBcXvA/EZ8y7EWha1uNx4Pg7AVuBPzMeW4AOZT2eQo498zVr6xFqZ+xFft/Tm4QqpZRyKD3Ho5RS\nyqE08SillHIoTTxKKaUcShOPUkoph9LEo5RSyqE08SillHIoTTzXEWNMjDFmVlnHUZKMMasKWopX\nFZwxpooxZrIx5rYSXu4iY8yeklxmWcgoX17sG52a/2/vzIOsKq4w/vvUACpihbAEtyiiFZWkXCHG\nDYlRIqUIKomgESsJaDCJRqvUuIBbGbcsZYxBUVwikeASo0NACh2NJiBoVNSIMYIKjiCLgDA4Rez8\ncc7FnjvvvXkz82YGqf6qXs3cc0+fPre7by+nz+0j7V5OKHSvi+XR9QBP17eYTpJGSzqppTq2Bz43\n4WkTysJQ7EO+hITG0AEYh30Y+FIF5V6NBVBLaBomAo81wnMoFmgtw2jgVez0iM8V0sCzBSGE8Hk5\nlDJhC4OkbUMItSGEigT/y+RVQpbLE9AxhLChUjIriRDCYmBxIzyz20idVkcytVUImYlB0mBJr0ta\nL6lKUldJfSQ9JWmd83w9l/YCSXMlrZa0VNJjkvpE9/eQtFbSVRHtC5JekPSMpK2cVs/U1lydipkH\n8maUzDwgqb/LqJX0rOvbQ9JfJH0s6d+SBpZRhrtKmuZyFkn6YRG+vv4ca/031QPQlZKd6XqApNle\nFv+SdESOb2vnfVfSJ5JekzQiuj9KUp0sBEKcbj8vs2Mi2rmS/uNy3pLUaIwiSdVuXhztZVDrz7pz\njq+bpHskrfBnqZYd0hjznOhtZJ2kVbLooFko57X+d5LrHSTt7uk6SbpB0nuu+8uSjs/JXiTpZkmX\nS1qMHZlS0NQmaX9Js1zPVZLuVxQoLmpvIyXdK+kjisz+I94Rku7z+l8maVyOL6vvwyXNBTYAp/q9\nPbxtrvH09d63CB0k/VbSSkkfSbpF9Q8E7SXpLklvez29KemamCdCl3L0LfTMEc8mU5vM/HwQcGZU\nf6O83t6WpFzarN12L5VHm6G9zwnaUn7YmWXLsFgkw4DTsWPyHwTmAWOwQFEvAa+DHVfkaX+NHbI4\nADgRmOaydox4xmCHDx7s11diBxn2jngWATe1VCeKnDfm8uZF1+Ox0L8vYycxnwS8CzwLzMJikhyL\nnd+1Aj9wsEj5CXjR049wfecDS4DqiK8PsNrlD8Hin7wOzI3LtID8TNdXgLP8uWdjB1puF/Fd6+V8\nGXAccLuXxWnhs/OpPgHOysm/CvgAP98L+JGnu9nL4DosVMDFjbSjan/m+V4GI4D3gLk5vmc9v7Ow\nyJ/PYINJH7+/J1CHxQsaCBwPXA4M9ftHu35X89mhrh393uPebs5x3ScCG4H9c22txuv2RGBYkTbS\nHTvL65/ePk7HZvavAB1y7a0GuBX4NjCwSPlkvEuACV5H13rZji1Q3//F2vnRWMC+jtjBlguwE6hP\nxsxVS4CuUfosj4e8rVzo9X5jxPM14CZ/rqO8zpcAE5qp7/LoeoCn65vT6Vz/f18sLEVVVH/dsbAE\ngdx5cdgZcg+1dz+5SZ/2VmBL+fkLtxHYM6Ld4I3g+xHteKftU0TO1piNfG2czu9Nx06O/SbWOZ6T\nu7+IhgNPk3WiaQNPAI6KaD922hURbV+nfadE+WU69I9oX3H9qyPafVin0SGi7YUF3xpcQn6m68CI\ntr/TBvl1Vyx087hc2mnAguj6UWB6jmcBFv4azJKwBJiU4/k9Nmh2KqFntdftbhHtsJyegwqU+/bY\nIDrBr08BVpTIpzMFDnLFDr2sJ9vpzwBTc22tJv8sBdrIL7GBp0tE60/9wTxrb4+U8Z5lvE/k6Hd4\nmW+Vq+8hOb6zvU3FE7ZdsEH6kogWsLgyW0W0S7HBrGsR3bbBJgobaDiolqNv2QOPX8+jwIG12KTk\nnui6NzbQbTYH1yZTW2WxKNS3cb/lf58sQNtkOpH0DUkzJa3AXor1WMewd07+D4CdXN5TIYTbWkun\nJqAO+HsF5PcDloYQ5mSEEMI72GotxjHAI8CnkraRtA224boIOJjSqCOKlomtlMA6HrA4Ktth4Yxj\nTAH2jswUU4BvSfoSmCkJq6spkbydisjpgs2US+HFEMK72UUI4TlsBdLPSf2AZSGEpyOeddhKJYuP\nMh/Y0c1xx0ravpE8MxyDraSey8rXy3gWDct3Vmh8z6Qf1umuiXSdg9XX4TneqjJ1BGsDMR7GynyX\niBZoGBumH1a+b0f6LAaeK6DPoyGEOKDdw9iksC/YvpGk82Rm7FpswnA/tqrKh74uR99K4U7gZEmd\n/XoUsBSbuG4WSANPZfFR7rquAD2jdQKQxWZ/AjM1jcFmt4dgHU2nWFgIYQnWyXfEjqJvFZ2aiLW5\nl7OB/GBRGRuT/2XsmfPI07oBF2EvefzrDezaFF0L6NXL/y7Npcuuu/rfv3qeJ/v1dzHzURYEr1w5\nxVCsHDK5vYrwLM1khxAWYKbI3tiKbbmkyWXY+LthdZEv3/E0LN/88xVCryJ8m3RtorwM+efPrntF\ntFVRHTdHn8byOA8ztT2ClXU/YKzfy7f1cvStFP6MrXCG+17PmcC9IYSNrZBXs5C82tofg7BZ9hCf\nteIzzAadk6ShWLCpl4EbJM0IIaxvBZ2yWWx+k/SLrZBXhg+AHgXoPbAYKRlWYi/6xAK8JTdny0BN\nlGfslp5thK8ECCF8LKkKG3BuB4ZjZqhQQE6MenJKoFg5ZHJrivD0jGWHEKqAKllU0MFYePJbgO+V\nyHslZgIq5/uQ0DhLSV3zq9ly5GXIy8yuayJaIXk1wH5F9MnXS2N5nAo8GEK4NGOQtG8L9K0IQgjr\nJD2ArXTewVZfkyqdT0uQVjztj22x2Uk8GxlOblIgqRvwB2yfYBA2MF3fSjotw2a5+0T5d8b2lloL\nc4GekvpHee4GHJjjm4V1HC+EEOblfotaqMOrmJnz1Bx9OPBmCOHDiPYAcJSkE7BVxQPRvcXA+0Xk\nrMHMYKVwoD87AJIOwzqq5500B+gh6ciIZztscGkQejyEsDqEMBkbsLOOsdgqdBa24vm4QPk258PQ\nOeTCSks6BNv7aEmY9KG562FYJ17SJdn1OUjSHpE+O2NtO6/PELnHaJRHLdZOwN7dT3JpioU7b66+\npVBHcSvCncAR2Ep1dgjhjRbkU3GkFU/740nMoWCSpDuxTvVCGprIbsMcDi7yGc1YYLKkh0MIT1VS\noRDCp5IeBc6X9I7rcgH1Vx6VxjRsJTdV0kXYC30lDU0U47EOuErSXdgqZ2fME+ruEEJ1cxUIIayU\n9BvgMkkbsc3bYZjjw2kF9F2PeSotDCE8H8n5VNJ4YILv283EvJ7OAX5Rxr7Ih/5847CO5XpsX2K6\ny58h6R/AFEkXY6uzC7GO8EYASWOwDw6nY4PgXthAeK/LqJO0EDPHvIqtcl9xXWcAMyVdjzmzdMEc\nMTqFEC4pqzA/w6/8uWe4vM6Yw8F8zGOsudhP0gSXcSS2//mznNm3EO7GTLV/k3QF5pQyDmtHE3K8\nO2Dt8Q7svbwcuDWEkK2MZgI/lTQH854biXldVlLfUngDG9SPw9rAwhDCCrB9NEmvYftWY1qQR+ug\nvb0btpQfOW8ep43ClvudI9ru5DzGgDOwhluLufj2J/JQwzq9/wGH5+RPdb4d/HpTmgro1BPz3lqD\nLddH5+WR88Rx2gBy3jhOr+eRU6QMd8M6ylrPcwzm+l2d4/uq01c671tYp7FLCdkNdC2kFzYJuBJz\nYa7DHBBGFpH5R09/XZH7P3Hd6jAX3vPLaEfV/mxnY67ltdgG+a45vu7YILLKeZ4GDonuH4pt1r+P\nDSoLsQGsY8RzLDbYbKB+mOOOXgaZ7h94vQyO0tZra420uQOwCdZ6bBIzGehZqv2VKJ+MdyTwJ2wy\n9qHrG3+iULC+/V5v7Gv/tdgnCY8DexVoFz8HfudlvBpz9Y7LrzNmwsrCXk/ETOGb2n9z9aU8r7be\nmDv7agp7KF7jZd6lsXJt618KfZ2QsBnBPwxcHkI4pb112Rwh+8h1IXBCCOHx9tVm84ak57HPAM5o\nb13ySKa2hISEhC0IshMsBmLesWMbYW8XpIEnISEhYcvCXMykeUkIYW57K1MIydSWkJCQkNCmSO7U\nCQkJCQltijTwJCQkJCS0KdLAk5CQkJDQpkgDT0JCQkJCmyINPAkJCQkJbYo08CQkJCQktCn+D3wb\nf6hWlAfzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0599b5f950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Making the plot of posterior probabilities in the paper.\n",
    "with mpl.rc_context({'font.size':15}):\n",
    "    f, ax = plt.subplots(figsize=(6,4))\n",
    "    ax.hist(pprobs_max_by_tissue, 50, stacked=True, label=labels, color=tissues_max_colors)\n",
    "    with mpl.rc_context({'font.size':11.5}):\n",
    "        ax.legend();\n",
    "    ax.set_xlabel('maximum de novo posterior probability')\n",
    "    ax.set_ylabel('count')\n",
    "    f.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "87"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the number of sites that are de novo according to OPL.\n",
    "# This will vary a little with MCMC sampling noise.\n",
    "is_opl_denovo = pprobs_frame.groupby(level=[0,1])['pprob'].apply(lambda x: np.any(x >= 0.8))\n",
    "is_opl_denovo.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See `germline_de_novo_regressions.ipynb` for regressions of number of *de novo* mutations against age."
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [anaconda2]",
   "language": "python",
   "name": "Python [anaconda2]"
  },
  "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
}