Skip to content

Instantly share code, notes, and snippets.

@hardingnj

hardingnj/poisson.ipynb

Created February 27, 2018 10:01
Show Gist options
  • Save hardingnj/41b59011b511871cef855310dbe5d1c4 to your computer and use it in GitHub Desktop.
Save hardingnj/41b59011b511871cef855310dbe5d1c4 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
poisson.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import allel\n",
"import pomegranate"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"('0.9.0', '1.1.9')"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pomegranate.__version__, allel.__version__"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create test genotypes"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"gv = np.random.random((100000, 2)) < 0.05\n",
"gv = gv.astype(\"int8\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Positions\n",
"pos = np.random.choice(range(0, 1000000), size=100000)\n",
"pos = allel.SortedIndex(np.sort(pos))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### add known ROH:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"roh = [[0, 10000], [100000, 150000], [600000, 800000]]\n",
"\n",
"for start, stop in roh:\n",
" loc = pos.locate_range(start, stop)\n",
" size = loc.stop - loc.start\n",
" \n",
" new_gv = np.random.random((size, 2)) < 0.001\n",
" gv[loc] = new_gv"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"is_access = np.ones(1000000).astype(\"bool\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'/home/njh/miniconda3/envs/poisson/lib/python3.6/site-packages/allel/__init__.py'"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"allel.__file__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Poisson model"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"hetp [ 0.001 0.0025 0.01 ]\n",
"start [ 0.33333333 0.33333333 0.33333333]\n",
"tx [[ 9.99000000e-01 5.00000000e-04 5.00000000e-04]\n",
" [ 5.00000000e-04 9.99000000e-01 5.00000000e-04]\n",
" [ 5.00000000e-04 5.00000000e-04 9.99000000e-01]]\n",
"CPU times: user 28 ms, sys: 1.05 ms, total: 29.1 ms\n",
"Wall time: 28.7 ms\n"
]
}
],
"source": [
"%%time\n",
"pdf, pfroh = allel.roh.roh_poissionhmm(\n",
" gv, pos, transition=0.001, is_accessible=is_access, window_size=1000)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"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>start</th>\n",
" <th>stop</th>\n",
" <th>length</th>\n",
" <th>is_marginal</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>10000</td>\n",
" <td>9999</td>\n",
" <td>True</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>99001</td>\n",
" <td>150000</td>\n",
" <td>50999</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>600001</td>\n",
" <td>800000</td>\n",
" <td>199999</td>\n",
" <td>False</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop length is_marginal\n",
"0 1 10000 9999 True\n",
"1 99001 150000 50999 False\n",
"2 600001 800000 199999 False"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pdf"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.26125825825825827"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pfroh"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multinomial model"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/njh/miniconda3/envs/poisson/lib/python3.6/site-packages/hmmlearn/hmm.py:377: RuntimeWarning: divide by zero encountered in log\n",
" return np.log(self.emissionprob_)[:, np.concatenate(X)].T\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 8.57 s, sys: 43.1 ms, total: 8.61 s\n",
"Wall time: 8.6 s\n"
]
}
],
"source": [
"%%time\n",
"mdf, mfroh = allel.stats.roh_mhmm(gv, pos, is_accessible=is_access)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"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>start</th>\n",
" <th>stop</th>\n",
" <th>length</th>\n",
" <th>is_marginal</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>10171</td>\n",
" <td>10171</td>\n",
" <td>True</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>98890</td>\n",
" <td>150047</td>\n",
" <td>51158</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>599932</td>\n",
" <td>800807</td>\n",
" <td>200876</td>\n",
" <td>False</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop length is_marginal\n",
"0 1 10171 10171 True\n",
"1 98890 150047 51158 False\n",
"2 599932 800807 200876 False"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mdf"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.26220500000000002"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mfroh"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## plot"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'Multinomial model')"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAEyCAYAAADZQ4ufAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYHFd5L/7v2zOakWaRRrJGsqzVYNlmsbEWb4FfLjeEJZCLkwcHZAMxxmCSG9Z7fzfXDjchhktibggXEwhYXohjvBFCErOFxWF58sM20sj7BopsjUaS7dEyo9GMNEv3+/uj6vScPn2qurqnq7pn5vt5nnmk7qmuOn3qVKl03nrfElUFERERUTPLNboBRERERJXwgoWIiIiaHi9YiIiIqOnxgoWIiIiaHi9YiIiIqOnxgoWIiIiaHi9YiIiIqOnxgoWIiIiaHi9YiIiIqOm1NroB1Vq+fLlu2LCh0c0gIiKiOujr6zukqr2Vlpt1FywbNmzAzp07G90MIiIiqgMR2ZtkOYaEiGjeKRQUgyPj4LPUiGaPWTfDQkQ0E4WC4rKbHkDf3qPYsn4p7nr/RcjlpNHNIqIKOMNCRPPK4dEJ9O09iqmCom/vURwenWh0k4goAV6wENG8sryrDVvWL0VrTrBl/VIs72prdJOIKAGGhIhoXhER3PX+i3B4dALLu9ogwnAQ0WyQyQyLiKwVkR+LyFMi8oSIfCR8f5mI/FBEfhX+uTSL9hDR/JbLCXq723mxQjSLZBUSmgLw31X1ZQAuAvBHIvJyANcAuE9VNwK4L3xNlApmhhARzV6ZXLCo6kFV3RX+fQTAUwBWA7gEwG3hYrcB+J0s2kPzj8kMufgv78O27Q+gUOBFCxHRbJL5TbcisgHAJgAPAlipqgeB4KIGwIqIz1wtIjtFZOfg4GBWTaU5hJkhRESzW6YXLCLSBeAfAXxUVY8l/ZyqblfVraq6tbe3YvVeojLMDCEimt0yyxISkQUILlbuUNVvhm+/ICKrVPWgiKwC8GJW7aH5hZkhRJSmQkF5fklZVllCAuAWAE+p6uesX90L4Irw71cA+Jcs2kPzEzNDiCgNvEcuG1mFhF4N4N0AfkNEHg5/3gzgegCvF5FfAXh9+JqIiGjW4D1y2cgkJKSq/w4g6r+1r8uiDUQ0O3GqnZqduUfOPJ9qNt8j18zHGyvdElHT4oMKaTaYK/fINfvxxmcJEVHT4lQ7zRZz4R65Zj/eeMFCs049K9YmWZdZJp8vsFKuZSb7IelnmY5O8419bGRdnbvZjzeZbSffrVu36s6dOxvdDGqQek5ZJlmXWWbnc0fQ0d6KsYk8tjbhVGnWZrIfqv1sM8fUierJPjY2r+sBINjVn214phHHm4j0qerWSstxhoVmlXpOWSZZl1kmr8DIySnkm3SqNGsz2Q/VfnYuTLUTJeEeG3392Ydnmvl44wVLQllN083mB/Rl0fZ6Tlm661rWsaBsH5/SuQBb1i9FiwDdC1vR0qRTpVmbyX5o5LTzbD6+qNxc25/usVGv42Su9BNDQglkNU3X7Hdox8my7fWcsjTrWtaxAJff/KB3H99x1YU4emISyzoW4MjYJEMToZnsh0ZMO8/m44vKzdX9aR8bqpjxcTIb+okhoTrKapqu2e/QjpNl2+s5ZWnWdWRsMnIfHz0xid7udrS05Jp2qrQRZrIfGjHtPJuPLyo3V/enfWzU4ziZS/3EC5YE0pqmq7Sd2RR2mM1tB7Lbx9Q4s32MUinuz2TmUj8xJJSQPU2Xzyt2Dx7HmSu7kMtVd81XaSq8WTIiamlH2m3Pcv31mIptFm6/NcsYa4Sk39233Hzut2bVjOepqalC5L8PaW27nv+uNHOWECvdJmSm5goFxTtveTC1dE6znUaqNeaZZtuziMPa7RdBw/dDPbj9dsdVF9Y8fueCJGPUN9YANP19APNRteectM8jU1MFbPrfP8TIySl0L2zFQ//r9WhtzaW67Xr+u9Ls97swJFSlLNM5G6UZ29mMbZoN3H7bPXic/ViBb6xx/M0Nae/H3YPHMXJyCkBQBmH34PHUt511qYdG4gVLlZZ1LMA5a5YgJ8C5a5bglM4FiT/riyVWSjeL+n2aaWq1xjx9bSoUFC8cO4kXj51M3FbfekybWmro96TbiHs/6TqqqZxrlqmUMj+TfW3vy83rerC0YwE2R6RyV/puSdXy2XqP57j1VdqWOcbtFPa5dB/AfJZ0P8Ydo77Xxpkru9C9MAhcdC9sxZkru4rLq2rksVfr+Pet13ynWtbpK/VQ7fk7TbyHpQpmumzHs0fQ0d6CE5OFqqueuvdJxE2/RU3PZTFtV20cM2oafdv2+/GL544CAC7YsAx3X125qmnUd5uaKuD3tt+PRweGZ1Rtth79WmvYIC5E40uZT7LOJN938Pg4PnTXQ9gVbudvLtuM5V1txVTueo2vWj5b7/Ect75K27IrG79qbQ/+4QMXo6Vlekqf97DMfknu94gLo1YKq7r3sLhlMdxjr9ZSGb71rljcXrxPbSbHsCn1cNlND1R1/q4V05pTYKbLCgCOj+drqnpqp6lVmn6L+n0W03bVptPFTaMbff2V2xr33Y6emMRjA8MzrjZbj36tNWxQKUTjpszXY1/ncoKcCHaF69nVP4RcTspSuesxvmr5bL3Hc9z6kh5zeQUeHRjGkbHJ4u+auQIoJVdpP1Y6RiuFVVtbczh71eLiDbf2+qKOvVpKZfjWa77TTI4pt9SDkeT8nTZesMRwp9TssIRb9bSa0IdbRTVqejIqDJLl9PTUVAFPHzyGQqEQu5yvTeY9Iy4EYfim4+O2UY1K/V7NdLFvGnZ5Vxs2r1uKlpxgc8Tno0I0Zh9vWddT/J2qlrW1Uv9FqSbMUe8qtnHHRqGgyBcKOGf1dNt6FraWjDlfyCzqQZT2vvGFD4N91IMWATav6ymbPo8aGzMJMVH6ZroP7M+7592NKzpLxsSZK7u85+WoNvhuI3DHYS3HeNxxGnce9X3nuDC80QyhUIaEIsSFDcx0mal6qpo89OGbbjx6Irp6alQYJIvp6bg73n2iUkEHj49DAJzS6Q9B2J+Pmo6P20YSSfu9muli3zTstu33F7dx99UXR44BO0Szae0STBWAR/cPY8u6HtywbRM+cs/D2OW01a7GW8vUcTVhjnpVsY07Nkx/md9tWtuDu993Ibb+5X3FMdf3J7+Jd3/1FyVT5317/Q+itPeN3afuMtu2P4C+cPr9bk/YzR0bMwkxUfpmug98ny8UtOS8644J97wcFSaKuo3gjqsuDI7lcBzeWeMxHnXOTXIeTVK93T5/pzm7yJDQDEVNqZnpMrvqaTWhD3e9popq1ECICoNkMT0dd8e7j69NuZxg5eKFWLF4YWQIwoibjo/bRhJJ+72a6WLfNOyu/iHkFdjVPxQ5Bnwhmkf3B/t4V/8Qhk9MFn9nt7VS/1X67tWEOepVxTbu2HB/9+j+YezaN1Qy5vr6j5ZNnUc9iNLdN6ZP3WV29R8N+joi7OaOjZmEmCh9M90Hvs+75113TLi/jwoTRd1GsHvweMk4rPUY9x2nSc6j7neOCknZ5+9mCIXygiVCNdPi1UydVTvd3sjshKg73mtV6buk+V3rte649dQyZnyVdc2Uc71CNY0cQ3HHhu93F5y+tGTMXXD60rJ+inoQZZJqxXGhy1rGJbOHGm+m+6AeY6LSMeuO2Xof49X2R5JjpRkxJBQjybS4HSI6NDoBVUVOpq96k1TMTPq6Z2Er/uPQaE0Vdmvlq9o4k8qp7hSjW1HWF3KLq9xo+t2esqwUmoqbOXDDGe42TulsK9ue7/u54yDJttx+8b32ja+o72XvO0AS77OpqQJ++eIITulsS/Q/q7h12f0BAAVVHB2dwPKudqxYvBD5vBa3tbwr+B+mO86j9om97Bm9nTg8NlncTk4Ep3S24cXj4xgancBZp3aXjF87THlkrPKDLePGvD2G51KF5NkkyXk27nO+/e/+zl0m7vj1hVfjPm9v84WRkyVjttLx5ftdVLXdqPNm1HnYvJ6YyKOv/yguOH0pWlpa6rHLSiQNCfGCZQZmmv6WZB3mM9XeT9Ko75wk5lopdlpN6umithYcH88DCO6PuPN95e0BkqUG+9rWt7d0G93trRibzEemVbv3ZSRNBayU7qwK7Nhbuk4g+v6QWituTk0VcN6nfmD1afS9OJXW5fZnoaDY2T9U/Oz563sgksOu/tpSO+3vmANgbgs3+2jz2iV4+vnjGBkv7QN7/PjuhamkkeUGKJkk+6KaZeLGSr32u+/c4TufVdpupfcrjXv387ddcT5ecd33kVegRYCnrnsT2trqe9HCe1gyMNP0tyTrMJ+p9n6StCRtb9LP+2Kn1aSemn9YgeD+CF97ksa4fW1ztzEyXn7vhG8ddptqSe31pTu764zbVq0VN3cPHi/t073R9+JUWld5fw6VfLZv71DJ/q82tdP+jnYO2/Q+GsLIeHkf2OPHdy9MJVHfmfezNI8k+6KaZeLGSr32u+94jju/VjsOk4579/M//uUg8uG8Rl5Rci7KGi9YIkSlT9rv5wsFvPK0xciFqWluXDIqTmlzU97cFDqTFrq0Y0HZ/SSNSKc06Xg5AOesXhzZ3qjUUxM7zQF42apubAlT+85ZvRiFQsGbVrh0UWtkenlX+/SVflSfm3TjnADnOGmudh/aqYCb1/XgnNMWQwB0Wf+b6F7YWpKe6H5XkxZbbNO6Hm+1SDdd3I0pn9HbWZKW6LsPxL0HxE6xPKO3E51h33QvbMUZvZ3Fdsal/Z65squkTzevW4J8oVCWkmy+99JFrSXjd+mi1uJ3dVOEt6zvKRlL565ZXEzjduPobmqnLzXavsfK/v/s9D5ajO726WPmjN5OvHDsJAqFQvH7R+1Pd7vPD58obtsef75xOxvuBZgrfPvr4PAJDB47UUwbjqqMbc4LbhmCqBRne6zEjfOkqfDuePYdzy9d3oFXrp7+N8Y+Lnypy1ElF+z15xCcN815TlVLUpvN9zF995svW46W8ABrEeCC05e6XyczDAl5RE2dmWl6NxQBTE+dA/GxQN92fClvJoXOroa7ae0SXHfJK3H2qd0ApCHTz4WC4h033l8SnrjzfeXtjZt2nJjIY9Onf4jR8Tw623LYuKIbDw8MF9d399XTaYWP7BtCZ0QKa9J7WII2/xw79g6VbAOANwR07polaM1Jcfnz1/fgC5dtRktOsHTRArz9pgfw6ECQguxLs739ygvw9pvuL1nG7q+vvfcCbPmLH5WFa+zvdPnND5akJYqIN0ZuYt7/9Wu7iim8ZpzuePYIXrF6Mf7xAxfj3V/dUbJPNkek/QLT97As7ViAD9/9MHY8Vx6Ksvfx8ZNT6GxvwdhEHh2eEJ09Nl4YOYlDIyfxZ//yZDGN26SG2/en2KmdQQqov+LmxEQeb7vx53jiwDG8as0SfOVdW7Cso624jzaHx8yZK7pw+c0PFtdx/vql+OLlm7Gso3x/mrCUu113fF5648/x8L7hsr7hPSzZ8IVRL7vp/uJx27kghzNPXewd4+bzbhkCoDx8DKB4D+Hbb3oAj+wbih3nSVLh41L9zfH8yEDpdnwh1L69peeIqJILhknHfmTfEF61ZgkWtOSwq3+oLCR7+5UX4B03P1BM277tivPx0MBQw+9hyexpzSJyK4DfBvCiqr4yfG8ZgHsAbADwHIC3q2rj5ptC7tQZgJJpejdMAAB9YRprb3d7yVMx456S6aa8me2YFDcAOHR8vDg999C+YfR2L0Qul8PgyHjZtF8WTxc2aaFGX7+/vW7f2e3bc3gUo+H3HZ0o4NHwYsWsz0xRPjYwjELEeux+Xbl4YUkb3T4P2jwdjrC3YfchRIqpgPZl/K59w2gN09gHR8ZL0hnNZ+w27jk8isf3H0MhTG+219XXH4Q+3HBNUBlTittw0xJ7u9vLvqf5rq25HB7bX55iWQDw1MER7Dk8VrZPdvUPBW33jJ/W1hxeftoSDI6MY1dEyMldnxm/bojOHhsiwKoli0raa6eGmydku2PbfB+3Hb3d7Rgen8JTB0dQUOCx/ceQy+UwdHKquI/MMXP0xFTJOnbtC7ZrL1vcnxHbdfvgMc+4dY9/So8/jDp9nI9OBueWvPrPkea84JYh8J1X7WO/4AlF2+Pc1zZ721GhXHNeM8dH+XaGINb4tM9XJnXZbNMtuWCYdOxC+Dm46wv/vufwaMlxMTKZx8VnLK/HbpuRLENCfwfgTc571wC4T1U3ArgvfN1wlVLR3FAEUFrF1Q6F+Kb+DDsc5E5Nu8tUSuFMY/rZF9pxpy3tKUW7XQKgo62l2G47LLJxRWdxKr+rvaUkVGCWNdVPc+JPYa3U5snJPJ48MIwXhk8UwxbG2ad2YSo/hUKhgE1WeMuEJzatXYKXndpdXN5XZdiuTGuPh7NP7cIZvR1ly9jfb+v6nmK4prMth4Iq7ErC7r6t9IBCd0raDSfZ49YeZ3ZVXRP2mZzMF0NVUfv6lM4F2LyuBwKgs60lCJu1tyDnHBeb1y4pGRtR3y8uZbRYDdjpQ1/aqF1F1Peebx3u8WVXHzX9aHM/52sTZcOfWjy9j825JSpUY84LcdWf7XO6qhbHfUfb9D+d9rpNqHdZR6s3bGiHbHyft79XDsH5wTAhdDucXk0qtlm3r8Lu5nU9OHd1+Tmj2cKbmYaERGQDgG9bMyzPAHitqh4UkVUAfqKqZ8WtI8tKt3GpaG4as6niak+7b4nI7nDv2H7V2h7cddWF2HbLgyUVbQGULONWLJxJNdIk3z8qtANMp9h++O7piqxm2rMk5NPegl1htVJ7CvRr770Aew6PFlNu7TRT+4Fbm9b24OtXX4Shk1MVv6cdYgNKs0ZGJ6Zwzmnd+NWLoxibnL446GxrweiEmXJdis9vOw9v/PzPiv+z6Wxrwcmp0odc2lUuTWXaD93Zh539w8Xt7bj2dbjs1gdLlmnJTY+THc8ewctWdaH/yMmyLBbzXXzhkahMALcarjtm7NR4OwRyw7ZN+PDdDxXDPi0S3Fhn2pPLSdm+Nhcuxan3thacmMzjvHAcv+PmB/DIvmF0xWRUVRq7JvX4Q3c9hL4wBDs6nsd5a5fgG3/wayXHga8atK9S6eU3P4id4RT6N6wpdNNX97z/Irzzll+gL5xyFytE+Lfv3IyWXK6YAmo+d24YhmqWwlrzjS/8a6cF+9L57fPa6PgUzouo/mwfS2bZRQtyGJ0Izh/nrVmCG989ve/dzLwd17wO2255oBg2PH/9UohIMaRjzgm+0gd26Oac07qxe3AMoxN5dCwQiOQwOpFHV3sLfvSx/4SVS6bHXpLjyq70fOdVF+Lw2ETxOPOdM7IIbzZdSCjCSlU9CADhRcsK30IicjWAqwFg3bp1mTTMDitEhXjsKXp7Gt8OD7jhADP1Z4edHh0YxrNHxrwVbX2hAV9b6i0qLGbav3LxwmLIwJ32LAn5jOexa99Q2RTo8PgUzl61uPie6UvTj8aj+4cxdHIq0fe0Q2w2kyny2P6Rst+ZixUgCBP0Hx4rmYY1v7e/n13lsliZtn86PDAyPoWHBoZKlrFDStPhmuPF9tihIQBl4SG3j90pZzMlHRVOstdX1narv002gN0e3762x7XpIzOOH99/DGr1u286vtLYtasB2yHYx/YfKzsOoqpBu5VId/UfDUJHzhS66av/ODRaUn3UTLk/tv8YWltavGFPE4bixUpjuOMolxOsWrIIq5YsKr7nC9XY5zX33OoLzZplzcUKADx2oHTfu5l55hxg9PUfLYZ07HOCjx26sc9bY5MKYDr8OnRiEqf2TH/XSseVCenbFXbt48x3zmgmsyJLSFW3q+pWVd3a29vb6OZ4RYWRfNkd9vJxGUWNzDqI+j6+qUu3fW6FXFOt1NcPUdtNsmzUZ3MoHdjd7ab9PWWhPDfL6ILTl5a8N/3ZSvtueiraV6HV/9meipWEo/q42vej1mfGnWGyAdz2uJ+zPxMVNq0mlOeTZAxGfedaji93mah+ZDbQ7JV0TPmWdUOe7ueSnPeSjpuo80RnW64k863a6uNJjpVmHs8MCcWoVOnQXdatnmlnRhw5Po5lnW3FaWUzhbh78DheckoH9hwewxm9nTh6Yjr0EbXOrP43Z7dv96FRLO1YABEphsCi+sRMyx45Pg5T0dT0g1t11NeHqoqCKobHJovTuoPHx5EvFEre8+0buwLqMy+OAGHfmX43lVXNdzHs6rF2pVdTfTUuNGh+b6qqblzRhaMnpsoqtpr+NPt5WccCPD9yEs8OHsfGFV1Y3rUwclu+qrd2f/iqWfqqX5rv3pLLlexDk2m1dNGC2Cqz9li2x4IbHjWvff1bbWVkNzTrqyYbVw03rgq1m3kVVcE2ye+yeEgclYsbI77zg121Omll5Z6FrfjV4PGSc6BvzOTzWlJd2q1O7atgnWS85fOKp184hpwIzljeWfLvRVQF9GqqMTf635qmrHTruWD5KwCHVfV6EbkGwDJV/eO4dWT9tGY7fulLr7WXTZLCBpRWMNzx7BGIc9+AXY3TTrmtpgpoPb+/m8INRFd8rZQSHpWG7faTWb97HxAQ1EV5+WmL0bf3aOS+AeKfoF1NdcpqKkoC/rTu268sTWXu+5PfxLtufbBkXJh7bXzfxR4Lbn8kqabrVrC1UySjvr/7/W6/8gJs/vQPy6rgun1tYvVuleBKlTuTSFoFuVJKtN2nSdpRa5oqpct3vnHLDZjzg+88UShoyX0n9tPBtyQ4l9RabdY9r1c6x7v3x5h21uvfj5lWf56ppruHRUTuAvBaAMtFZADAJwBcD+DrInIVgH4Av5dVeyrxxTqj0nSrSWEDSisYFgCYGwLs+wbcdUaloGbx/d2LFSD6/oRKKeFRadhuP9nrdy+pj0/k0bd3KHbfmL8b9v1D9vaSpIVHLRtVUdL3/d1UZvMUYl+fRn0X3/0jvu/m46tgK7n4MeV+v77+o5FVcN2+lvD+jySViKsZy5X2mz/VNTo1O2k7ak1TpXRFnW/ccgNR54nDx8djnw5e6VxSzbnB935UWrE7ftz7Y0w76/XvR6V7FptFZvewqOplqrpKVReo6hpVvUVVD6vq61R1Y/jnkazaU4kvfmnuZ9jsxPmiUo/t9djOW7M4SNUM77cw9w10trfgpcs7MDgyjiXtLTj71K5ihUNfilxW39+97wMI+kQQpNf1LGwpS3u2qym69zb4qoou61iAc1YvKVl/Lvy+51rvA9PpimYbJracE2DT2qAyq6lmaszkPgRfmmyhoMXUa19apBsjj3oKsa0z/A5ubN1eZ9Afi0s+56bCm2qfTx0YLqZLuxVst6zvKe6nl63qxrKO0v+72N/PpEGfv6GnrAqu3T67r31jx71PZNPaJTg0crIkpds2NVUopqa7KaG+lG93n2z23KPjpj7HpY7bJQlMKrdJB/Ud//b3bOb7AOaSIE13aTHV2D1fmvNUV/t0+r19bNn3nXS1t2D9KYumx4+Vxu+7ry7qHGDa5RtfZdWf1/WgJSfYZFXWPmf14rLzo93OzvYWnL+hp5hi7d535jtO3PMX4K/o6563m20cs9JtDDeF+YN37iqm8JoHwrnpyW7qsVnPwWMn8IbP/QyjE/lggIkUq3z2LGzFpdvvxxMHjqEzrBwqmE7LPX/DUtzx3gvx9pvuL6usmeZ0nS+FGwjuR1jS3orNf/EjjI5Pfx83pfSRfUMlKYNuOrA9HWqHOc5bswT3vP8ivOPm6ZTA89YswZfeuQnHTkzhrFO7USigpGJja0suCBGVVIZcii9cFp06WM29FG7b7fZuWtuDb/yBPy3SjgW7T1A1ceN8oYA/uuMhPDIQ9Nc97y9P456aKpRUVgWAc0/rRktLS0klTwAlVX2721vx0J8G08TuU5gnJwvF9HN3OtmeArdTu2+74nxcetP9eOrgSHAcWFPdvth8VCXiF0ZO4g3/92eRD/OMC2GZKp7Lu9qKKd++UNn565fing8EISu3GrBdwdSXOg74wlxB8S73ezPFuXHcKtZAEKq8830X4eiJyWIa/8P9Q8VqzG4aszkuPvEvjxfXc96aJWhrDarAmtDi4bHSpxvb48M9B5i2+caXXYIgGDvhecvKWLRToM2YnJoq4G03/hxPHhwp1nbq6x/CprVL8ElPBXS72m0+X57m74Zmo87bWWi6kNBsZKd15UTwkFMV0Zee7KZcmvWMnJgqpn/mFYBOVyM8NpEvVus003EllVb3DhUrDxpZTDvHVZN9+uCxYuqy+T5uSmnB6RM3/dSeprS/72MHjhXTvO332lpb8bLTOgEAh0fHyyo2FpwQxK598amD1aTtuW232/vofn9aJFCaUtnamitJ5c7lpJgybCpbPjrgT+M227c9fnCkpPLl9PT19Ml7ZHx6mthUsDXs9HPfdLJhp3Y/d3QMTx8cKaZF2hU67TFiqtYC/krEw2OT3mq/RlwIyxw3R8YmY0Nlu/ZNH6duNWC3gmmSMEAxJdUzxc8U58Yw1Wptff1DJWn+j4WVq814cs/Tra059HYvLFnPo/uHS46toycmS8axSW23l48qO+GOL7ttu/qHgkrnE6Vhd3u82WPyqfDYs8+fD0dUQLer3R49MVGW5u8b877zdjOZFWnNzcAXQoh60JTvgWn2lF5OUFLV1g0jCEp3zJb15VU33QqzWbO/T4ugGB5a3JbDoZGT2GxVUVXVkofuuWmjpsqi4Uu3tacno6Y8BcCiBeUVKH0PSPNN/8c9SDIurdcXpqu0vrh12+EJt3Kmza5Waz4XTFFP/8Pvq6RruNPMJhzpToF3t0+H/0wVX18oza3qbPeB+zs3BdSkZ5rleha1ekNY9nf1jQFb3BR9XHjQrrRsr8sco688rRuDYSgrLhxM6XPDcYC/MrWpGmtCRss6FpSETd2w5ua1S0oqv0al0hsmjAOUh2TdkLJ5cKL9EE5f5fS48e5WpjWVxE31bhO67FnYiicPDCOfz5ccK9EPiQ0+u8mqUl3NuSxtDAlVwU0T8029uVOFwHT4xkzpPb7/GDqthx2609Mm1e5Xg8eLqbWm4qKpuvmRex4pqzCbNRPi2LB0Ebb85X3F/60DwQH4g4+S5pZeAAAgAElEQVT+Oj769UfQ59zBb/eVXc3UnVKPSjv19bsd3uhoy+G+j/0nnNqzqGQ/melYdyoUSJYx4kvx/cPb+8oe3Jh0fb5129PHm9YuwdPPHy9Wwu37k9/E0ZOTxRCS/fDA5V2lFYJftXoxPvW7r8Q7b/qFt5KuYT88sLPd/wDHTWsW4+kXRjE6kS+2Y3i8NGTlhpHcqp6+is+FgpaFyOx1bF3Xgz+/5BXotVLj3T5yjz03jdSe7YgLAdr9b/ejXc00n1e87Sv/Hx4ZOBaMcStbLSocTOlxw3FfvHwTPnjnw2UPO3TDqeevD8bjzv7SsCkAXHrjz4tVmkcnyqvg2qamCrj0Kz8vOf7vfF/pAxjtdZtwixu2/uLlwfHrpt5Hjfcbtm3CR+55uFiZNqjQPJ1x2NXeihOT+eD8cXCkOHtz/voefPHyLcUHIvqqBG/b/gB27j1SPBe4DwRN698ahoRSYE/12w8ltKfe3KlCoDQz4amDIyXTk+70NDA9lW6m7+2Ki4/tP4bhE1PeCrNZMyEOOzxkHB/PY+/hsWIFRfuhe/bD7uwqi+6UuhtmAEqzNex12eGNsYkChk9OYZW1P9ysEd/0f6X+tPe/CIphBqPWDBR73aXTx0MllXD3HB7F2asWF0NIvvCI8fjBEYyezJdkHrlhFwAlDw+MeoDjQ/uOedvh2y92X0hEyM8OZ9rrcdfx0MAwVi5eVNLnbh+548n3gEi3j5P0v2FXMw2m1Y8Vf2dnqzXrFPpc5objRk7mvWE/N5za1z8E+//pJmx6Sld7WZXmuP169MRk2fHvPoDRXrfvwYnmIZwtLbmysRs13k11artCsz1mj4dt32WdPwBgV/9wyQMR/Q+JPVp2Lsg6QzUO/ztQgS+8A5Q/oM1Mn/mmKN0pvagqi76HDfq2tWRRq/fO9Kz6wn2Q4xm9ncWH+Rl2Joz7fe2Qh7nL35d95W7bd5e9WZf9QEU7xOBO2dtToaUP8/O3IS6cVJYdE2aPRI2NJFOr9gMxz12z2FvV0hfC8D1Qzc1Msqtiuv0Z9aBPt9JmR1sOZ/R2RK7HcKes48J75viayudLwjC+7Cd3Kt1dphq+fauq2BTzoMXN66bbZ7LVzLiJekglpSOuUnHcOXnL+qUl+7F7YSvO6O0shgHtbCL7waBRD/G01+uren3myq7IrMuoDLVK39M+X7nh82Lbncre565ZXJYhlM8XSjLhfOeCqFBqIzAkFCOq6JtdNMiENMxD4UzBIvuBaeaBd/YU9vKutrJqpnGFe8y2PnjnruKD6nx3pmfVF274Y8ezR/CK1YvxD++/CM8dPVFWNdKuNOneMX/5zdPhE5N95W477kF/5v3br5x+oKLZti+Dy7ffVFGc9nQzwOLCSYWC4tIbf168+97+rG9sJCnetOPZI+gIQ4ab1y7BdWEGgK+KpW+62M6KcjOTKvWnr5rv8q42jI/nsfkvfoSxiemMolxOytZzeGyi+KBE+2FqIuIN77nHl+E+dBJAWbgoryib/q9mPPv27c7njqCjLQgHmIck2tlf27Y/gB3PBWP9mx+4GO+69RfFvq9UiI/qzxfWiDsnf+VdW7C8qx2X3fRg2X60s+JOTOaD7MOcFMM7vszMqLC1r6q3m3Vpqsr6Hm4a9z3NbQf2OdOchx7ZN4zOtpZiUc+vvfdC/HJwBH/6z08UjxV7rPsKO7rnAl9l3HpjSKgOooq+2ZkRJqRRUqDHeWCabwq7xcleqVS4x96W4bszPau+cMMfwcP8RjAyWSjLhLFDXe4d88FD6cqzr3zbjnrQn3nffaBiVAZX3H7zZYBVCieZ0In92chtxEytmu0VMD1t/JCVAWDzhZB8D1RzM5Pi+tPsI9++e+7oGMYmSjOKTulqL1uP/aBEt9+jwns+7kMnzd+NXfuGyrIoZlKAzuzbvFoPy3TCAWbKXAE8fXAEew6PlYzdZpo6ny98Dz+MOyfncjkcGZss24++rLhHw+wiw5eZ6Qtb+x7A6LbVfdhrNeHoQ8fHy86ZAIrhLHPPyq7+IQyPT2Hl4kUloTJ7rLuFHX3nAjvjr9EYEorhK/pWKcPAN32WpEhZkody+aYgs5qii9p2tQ/OipvGrdQ3M33Q30z3W7UPqEyyDV97q3lwYLX9X+tnfFk9te4Xtx0+UcXzjJlOVUft22oe9pn0YYmUrUrHXdR+NJI8vDatttZ6/MadO+LGetzDHJsRQ0IVmCk/+8Fu7oOh7Gl396F8vun1qAfARRUc87Un64esuf1gbztpATZfeMj3ULq4z7jbiNp20r60p1nj1pNkv/m24Rsbcf3kmzaO28/2eDAPHXQ/k6TNSfZhVHgpSb/FvWfGlRF1nLlj3zxozn1oaFJRx3bSMeM+BDPJ/qL0JTlnR4WSzHFkh0PiHpBZ6RxQ63GV5DvGHUtRD/l0x7j90NakD1FNQ1M+/LAeGp3WHPWwu2rSZH3raeZ4dz3a67tnwPdQunr0U7P070zakeSzboq3L3UYqJxinUV/1Xsble75SrqOejy0sFnGG6X/EL+k/wbEPQCx0vpqbWeltvnuo5tJ2+sp6QULQ0JV8D3QKul7ldbTzOrR3qh7Bgz7npiZbrdZ+ncm7UjyWXeZvv7y/qxlPWn0V7234d7z5VYErWYdhm8MVrOeRo83qs+4SLL+JOf7rI+9JG3r669f2xuBFyxViIsdllcMXBo8vM+TdllL3NKWdeXBerTXrUoaV8k26XaTpIHXMy5bbb/bcWXfAx+TfDbpvT2+asFJ7zFKq7/S3EY19/tEpeO798VUWz06KtW+2e8DmMtquQ+slvUnuY/NPT59Yyvqs7Wc4yu1zS01kLTtzTSmGRKqUtL7AewHckWFO2qJDzZq+rke7bWrkkZVsk263SRp4PWMv9ba71EPfEyyjiTfwU7j9FULrmY9acer09gnSe5Tigv72HH/D9/9cOLq0VGp4VnE+yleknFRj/UnvYfFTrOOukXA/fdjJqHkqHtbzDnCrco8k/tv6oUhoZSY9DJ7B7rvBemP09UOfVPNvvUk0aipunq0165Kata5cvHC2KfbRm230tRvre1N8j2q6Xf3oYnulGycJN/BTuOMegBf0vWkfbNovbdh1mdKBPjWWynsY8ZgSy5XVj06TlRqOC9WGi/JuKjH+pMcY74060r/FszkHB/VNvscYUoNJGl7s41pXrCkwJeCWa+wT5JpxrT42lRp6jKtqcW0p36jtlft9yhLsVzXU7Gqb5QkY2ImfVDLNHTW4clqt5n0WIzrQ9/2mnXKnPyiwsdZqHaspDG2fOtsxLE7UwwJpWQm6ceVpgSTTDPWWzV3x/s+m8bUYtpTv1HbqyUsZtpp95evqm/cOiqNiZn0cS3T0I0IT9baziTHYtR0etT2spwyp9qlnTmUtA3VjJU0xpZbxqGZMtsYEmqwJOGOKJWmBJNMM9ZbLdlQdnvTnJpNa+o3anvVbsd8LqiwWV6hMokkY2ImfVDLNHQjwpO1bDPpsejrw7jtNeOUOZVLO3MoiWrHShpjy15ns2YBVcILlibULJkdlbbHafHqzKS/0u7rWtbfiP3fDOOeZpesw8ezwWwd1wwJOZplmrdZMjsqba8Z2jAb1CN8lfZ3r2X9jdgfzT7msg5Vzif1CMtynwSa6d8YVrqtAStWNrfZun9ma7upes1wv8RcxeMoW1n2N+9hqcFsjevNF7N1/8zWdlP1muF+ibmKx1G2mrG/ecFiWdaxAOesCSoBZp0yPFNR1TwrLWentjV7mluzVLCt1myNF1P1qq3AG5Vq645J+7idnMzjyQPDODg0hueHT+D54RNlx33c56OWs5fJ5wslfzftdNfhO3/Yyw+OjGNyMo/H9w/h8f1DeGH4RMXjzF6P3Qa3WjaPo3Q143mLIaGQPZV77polWNCSw67+oVkx9Zj0IW7ucuevXwqRoPKq78FYzfidm6WCbS3bmY333lD1klbgjQodASiroms/KDQHoODZbtQDL6MeNGov5z48s6utBccn8gCA7vZWjE5MoaOtBcfH88V13Pm+C/HOWx4sfh4Q9O0Nvs/o+BQ621tx/OQUAEBL2hmd0m/3yyJre93trRibzGOLUy2b0tVs97C0ptaChETkTQBuANAC4GZVvb4R7bCnch8dGAZESqZze7vbG9GsRKKqebpt9i0n1vecDd/ZpObVi2/aM43vXe92U/Oy93XUPndDRwBKpt3jHhTqu1gBSqv4xn0+ajn7osJcrADAyHjQPnPxYNZh1jtlnz+s72P+LGvn3qHI48zuF3t7pg1utWxKV7Odtxp6wSIiLQC+BOD1AAYA7BCRe1X1yazb0rOwFWef2oUnDozgFau6UEAOTz0/UlNF0ixNTRVwYGgUL+3txC9fHAUAnLmiA/uHjuPJg0NY2rEAQHBw5wsFnLGiE8+8ML1cS0sLnjgwgpecshDtbQvwxIERnLmyEweHRzE4chJnn9qNXG7uRg7NtKf532gz72tqfoWC4oWRkzg0chIigqUdC3BkdAIKoKCK4ROTWNqxAKrAWeH5ZmGr4OSU4swVHTg4PIp8QfHS3k4888JxnLmiA+NTkzjr1G48eXAEQHA0++bF7c+fsaITzzx/HGeu7Cz7vFkuaEM3njxwDGeu6EAu14InwmVMmwBg0YIcTkwW0N4CmGsI064zV3biyYPHsbF3EVpbF+DxA8eKn120QHBisrylZ65YhIPDozg4PFbSJ4BAVXHWyi48cXCk5PMdYRvOWtmJialJPDYwBgXQmsvN+XMUTWtoSEhELgbw56r6xvD1tQCgqn8Z9Zk0QkJTUwWc96kflFzRG+ev78E9H/i1pgyPTE0VcO51/4oxz0mhXrrbW/HQn74era1z94TAcA3Vg/vQU8rGfDhHzXWzJUtoNYB91uuB8L0SInK1iOwUkZ2Dg4N1b8TuwePeixUA2NU/3BR3R/vsHjye6sUKEEzF7h48nuo2Go0VS6kegnAGL1ayNh/OURRo9AWL71+Isn+BVXW7qm5V1a29vb11b8SZK7vQ1d7i/d2W9T1NGyY4c2UXOhak+49sd3srzlzZleo2iOaCILzYk8m2Gn3irlaa7eU5av5o9E23AwDWWq/XADiQdSNyuRwe/tM34JcvjgTxZQBHRyewvKu9pmcBZSWXy+HRT7wJTz4/jOETkzijtxMKweHj48g7sWEgiKGPnJwqWU5EsGRRK549NIplnW1Y2tFWjLkzPkyUnIjgng/8WqJ7WIDgeWBLFrVi7+ExnL68A0dGJ4vLjZycwkus9wAgF67v2IkpbFzRhcNjk8gXCjg6OoEeZztxnzfLuXIiOKWzDXlVPHcouC9OITg6OoHFYTtfap073HW4649qr/1Zt0987VGg+Bl3GzxHzS+NvoelFcAvAbwOwH4AOwBcrqpPRH1mtjytmYiIiCqbNaX5ReTNAD6PIK35VlX9dIXlBwHsTak5ywEcSmndVIp9nS32d3bY19lif2cnrb5er6oV7/do+AVLMxGRnUmu8mjm2NfZYn9nh32dLfZ3dhrd1wz8ERERUdPjBQsRERE1PV6wlNre6AbMI+zrbLG/s8O+zhb7OzsN7Wvew0JERERNjzMsRERE1PR4wUJERERNjxcsRERE1PR4wUJERERNjxcsRERE1PR4wUJERERNjxcsRERE1PR4wUJERERNjxcsRERE1PR4wUJERERNrzWtFYvIrQB+G8CLqvpKz+8FwA0A3gxgDMB7VHVXpfUuX75cN2zYUOfWEhERUSP09fUdUtXeSsuldsEC4O8AfBHA30f8/rcAbAx/LgTw5fDPWBs2bMDOnTvr1EQiIiJqJBHZm2S51EJCqvozAEdiFrkEwN9r4AEAPSKyKq32EBER0ezVyHtYVgPYZ70eCN8jIiIiKpFmSKgS8byn3gVFrgZwNQCsW7cutQZtuOY7Ja+fu/4tqW2rXtw219Ns+P5EzSbNY5LK8Tw1fzRyhmUAwFrr9RoAB3wLqup2Vd2qqlt7eyvel0NERERzTCMvWO4F8PsSuAjAsKoebGB7iIiIqElFhoREZATTIRoTvtHw76qqi+NWLCJ3AXgtgOUiMgDgEwAWIPjwVwB8F0FK824Eac1X1vwtiIiIaE6LvGBR1e6ZrFhVL6vwewXwRzPZBhEREc0PiUJCIvIaEbky/PtyETk93WYRERERTat4wSIinwDwPwFcG77VBuBraTaKiIiIyJZkhuV3AbwVwCgAqOoBADMKFxERERFVI8kFy0R4v4kCgIh0ptskIiIiolJJLli+LiI3Iiid/34APwJwU7rNIiIiIppWsdKtqn5WRF4P4BiAswD8mar+MPWWEREREYUSleYPL1B4kUJEREQNkbRwXJlKheOIiIiI6qVi4TgR+SSA5wHcjqDK7TvBLCEiIiLKUJKbbt+oqn+rqiOqekxVvwzgbWk3jIiIiMhIcsGSF5F3ikiLiORE5J0A8mk3jIiIiMhIcsFyOYC3A3gBwIsAfi98j4iIiCgTSdKanwNwSfpNISIiIvJL8iyhNSLyTyLyooi8ICL/KCJrsmgcEREREZAsJPRVAPcCOA3AagDfCt8jIiIiykSSC5ZeVf2qqk6FP38HoDfldhEREREVJblgOSQi7wqzhFpE5F0ADqfdMCIiIiIjyQXLexFkCT0P4CCAS8P3iIiIiDKRJEuoH8BbM2gLERERkVfFCxYROR3AhwBssJdXVV7EEBERUSaSPK35nwHcgiA7qJBuc4iIiIjKJblgOamqX0i9JUREREQRklyw3CAinwDwAwDj5k1V3ZVaq4iIiIgsSS5YzgHwbgC/gemQkIaviYiIiFKX5ILldwG8RFUn0m4MERERkU+SOiyPAOhJuyFEREREUZLMsKwE8LSI7EDpPSxMayYiIqJMJLlg+UTqrSAiIiKKkaTS7U9rXbmIvAnADQBaANysqtc7v38PgL8CsD9864uqenOt2yMiIqK5KckMS01EpAXAlwC8HsAAgB0icq+qPukseo+qfjCtdhAREdHsl+Sm21pdAGC3qu4JM4zuBnBJitsjIiKiOSrNC5bVAPZZrwfC91xvE5FHReQbIrLWtyIRuVpEdorIzsHBwTTaSkRERE2s4gWLiLxaRH4oIr8UkT0i8qyI7EmwbvG8p87rbwHYoKrnAvgRgNt8K1LV7aq6VVW39vb2Jtg0ERERzSVJ7mG5BcDHAPQByFex7gEA9ozJGgAH7AVU9bD18iYAn6li/URERDRPJLlgGVbV79Ww7h0ANorI6QiygLYBuNxeQERWqerB8OVbATxVw3aIiIhojktywfJjEfkrAN9EFQ8/VNUpEfkggO8jSGu+VVWfEJFPAtipqvcC+LCIvBXAFIAjAN5T29cgIiKiuSzJBcuF4Z9brfcSPfxQVb8L4LvOe39m/f1aANcmaAMRERHNY0kKx/3nLBpCREREFCVJltASEfmcSSsWkb8WkSVZNI6IiIgISFaH5VYAIwDeHv4cA/DVNBtFREREZEtyD8tLVfVt1uvrROThtBpERERE5Eoyw3JCRF5jXojIqwGcSK9JRERERKWSzLD8IYDbwvtWBEw/JiIioowlyRJ6GMCrRGRx+PpY6q0iIiIisiTJEvpIeLEyAuBzIrJLRN6QftOIiIiIAknuYXlvOKvyBgArAFwJ4PpUW0VERERkSXLBYp66/GYAX1XVR+B/EjMRERFRKpJcsPSJyA8QXLB8X0S6ARTSbRYRERHRtCRZQlcBOA/AHlUdE5FlCMJCRERERJlIMsNyMYBnVHVIRN4F4H8BGE63WURERETTklywfBnAmIi8CsAfA9gL4O9TbRURERGRJckFy5SqKoBLANygqjcA6E63WURERETTktzDMiIi1wJ4F4BfF5EWAAvSbRYRERHRtCQzLO8AMA7gKlV9HsBqAH+VaquIiIiILElK8z8P4HPW637wHhYiIiLKUOQMi4j8e/jniIgcc//MrolEREQ030XOsKjqa8I/eYMtERERNVSSm24RpjT/P+HLn6nqo+k1iYiIiKhUoqc1A7gDwYMPVwC4Q0Q+lHbDiIiIiIykpfkvVNVRABCRzwC4H8DfpNkwIiIiIiPp05rz1us8+LRmIiIiylCSGZavAnhQRP4pfP07AG5Jr0lEREREpZLUYfmciPwEwGsQzKxcqaoPpd0wIiIiIiPygkVEFgL4AwBnAHgMwN+q6lRWDSMiIiIy4u5huQ3AVgQXK78F4LPVrlxE3iQiz4jIbhG5xvP7dhG5J/z9gyKyodptEBER0dwXFxJ6uaqeAwAicguAX1Sz4vAhiV8C8HoAAwB2iMi9qvqktdhVAI6q6hkisg3AZxA8u4iIiIioKG6GZdL8pcZQ0AUAdqvqHlWdAHA3gEucZS5BMJMDAN8A8DoRYQYSERERlYibYXmV9cwgAbAofC0AVFUXV1j3agD7rNcDAC6MWkZVp0RkGMApAA7ZC4nI1QCuBoB169ZV2Gztnrv+LamtOy2zsc1EcxmPSaJ0xD1LqGWG6/bNlGgNy0BVtwPYDgBbt24t+z0RERHNbUkKx9VqAMBa6/UaAAeilhGRVgBLABxJsU1EREQ0C6V5wbIDwEYROV1E2gBsA3Cvs8y9AK4I/34pgH9TVc6gEBERUQlJ8/pARN4M4PMAWgDcqqqfFpFPAtipqveGtV5uB7AJwczKNlXdU2GdgwD2ptTk5XDun6HUsK+zxf7ODvs6W+zv7KTV1+tVtbfSQqlesMw2IrJTVbc2uh3zAfs6W+zv7LCvs8X+zk6j+zrNkBARERFRXfCChYiIiJoeL1hKbW90A+YR9nW22N/ZYV9ni/2dnYb2Ne9hISIioqbHGRYiIiJqerxgISIioqbHCxYiIiJqerxgISIioqbHCxYiIiJqerxgISIioqbXGvULEfkWgMicZ1V9ayotIiIiInJEXrAA+GxmrSAiIiKKkahwnIgsArBOVZ9Jv0lEREREpSrewyIi/wXAwwD+NXx9nojcm3bDiIiIiIy4kJDx5wAuAPATAFDVh0VkQ6UPicitAH4bwIuq+krP7wXADQDeDGAMwHtUdVel9S5fvlw3bKi4eSIiIpoF+vr6Dqlqb6XlklywTKnqcHB9UZW/A/BFAH8f8fvfArAx/LkQwJfDP2Nt2LABO3furLYtRERE1IREZG+S5ZKkNT8uIpcDaBGRjSLyNwB+XulDqvozAEdiFrkEwN9r4AEAPSKyKkmjiYiIaH5JMsPyIQAfBzAO4C4A3wfwqTpsezWAfdbrgfC9g+6CInI1gKsBYN26dXXYtN+Ga75T9t5z178lte3Vg6/N9dLs352oGaV5TFI5nqfmj4oXLKo6huCC5eN13rYvxuRNWVLV7QC2A8DWrVsrpzURERHRnNLIwnEDANZar9cAODDDdRIREdEcFHcPy2cB/DWAZwGcAHBT+HMcwON12Pa9AH5fAhcBGFbVsnAQERERUeQMi6r+FABE5FOq+uvWr74lIj+rtGIRuQvAawEsF5EBAJ8AsCBc91cAfBdBSvNuBGnNV9b4HYiIiGiOS3LTba+IvERV9wCAiJwOoGK+tKpeVuH3CuCPErWSiIiI5rUkFywfA/ATEdkTvt6AMGOHiIiIKAtJsoT+VUQ2Ajg7fOtpVR1Pt1lERERE0ypesIjIAgAfAGDuY/mJiNyoqpOptoyIiIgolCQk9GUEN8v+bfj63eF770urUURERES2JBcs56vqq6zX/yYij6TVICIiIiJXkmcJ5UXkpeaFiLwEQD69JhERERGVSjLD8j8A/DjMEhIA68GaKURERJShJFlC94VZQmchuGBhlhARERFlKkmWUAuANyKov9IK4HUiAlX9XMptIyIiIgKQLCT0LQAnATwGoJBuc4iIiIjKJblgWaOq56beEiIiIqIISbKEvicib0i9JUREREQRksywPADgn0QkB2ASwY23qqqLU20ZERERUSjJBctfA7gYwGPhE5aJiIiIMpUkJPQrAI/zYoWIiIgaJckMy0EEDzz8HoBi/RWmNRMREVFWklywPBv+tIU/RERERJlKUun2uiwaQkRERBQlyT0sRERERA3FCxYiIiJqerxgISIioqYXeQ+LiPwNgMhUZlX9cCotIiIiInLE3XS7M7NWEBEREcWIvGBR1duybAgRERFRlIppzSLSC+B/Ang5gIXmfVX9jRTbRURERFSU5KbbOwA8BeB0ANcBeA7AjhTbRERERFQiyQXLKap6C4BJVf2pqr4XwEUpt4uIiIioKElp/snwz4Mi8hYABwCsSa9JRERERKWSzLD8bxFZAuC/A/h/AdwM4GNJVi4ibxKRZ0Rkt4hc4/n9e0RkUEQeDn/eV1XriYiIaF5I8iyhb4d/HQbwn5OuWERaAHwJwOsBDADYISL3quqTzqL3qOoHk66XiIiI5p+4wnF/rKr/J6qAXILCcRcA2K2qe8L13Q3gEgDuBQsRERFRrLgZlqfCP2stILcawD7r9QCACz3LvU1Efh3ALwF8TFX3uQuIyNUArgaAdevW1dgcIiIimq3iCsd9K/yz1gJy4lut8/pbAO5S1XER+QMAtwEoq++iqtsBbAeArVu3Rj4ugIiIiOamJIXjtgL4OID19vKqem6Fjw4AWGu9XoMgw6hIVQ9bL28C8JlK7SEiIqL5J0la8x0A/geAxwAUqlj3DgAbReR0APsBbANwub2AiKxS1YPhy7diOgxFREREVJTkgmVQVe+tdsWqOiUiHwTwfQAtAG5V1SdE5JMAdobr/LCIvBXAFIAjAN5T7XaIiIho7ktywfIJEbkZwH0Axs2bqvrNSh9U1e8C+K7z3p9Zf78WwLWJW0tERETzUpILlisBnA1gAaZDQgqg4gULERERUT0kuWB5laqek3pLiIiIiCIkKc3/gIi8PPWWEBEREUVIMsPyGgBXiMizCO5hEQCaIK2ZiIiIqC6SXLC8KfVWEBEREcWoGBJS1b0AelkPDu8AAAkUSURBVAD8l/CnJ3yPiIiIKBMVL1hE5CMIisetCH++JiIfSrthREREREaSkNBVAC5U1VEAEJHPALgfwN+k2TAiIiIiI0mWkADIW6/z8D/YkIiIiCgVSWZYvgrgQRH5p/D17wC4Jb0mEREREZWqeMGiqp8TkZ8CeDWCmZUrVfWh1FtGREREFEoywwIADwM4aJYXkXWq2p9aq4iIiIgsFS9YwoygTwB4AdP3rygAFo4jIiKiTCSZYfkIgLNU9XDajSEiIiLySZIltA/AcNoNISIiIoqSZIZlD4CfiMh3EDxLCEBwM25qrSIiIiKyJLlg6Q9/2sIfIiIiokwlSWu+LouGEBEREUWJvGARkc+r6kdF5FsIsoJKqOpbU20ZERERUShuhuX28M/PZtEQIiIioiiRFyyq2hf+9TxVvcH+XfgE55+m2TAiIiIiI0la8xWe995T53YQERERRYq7h+UyAJcDOF1E7rV+1Q2AReSIiIgoM3H3sPwcwfODlgP4a+v9EQCPptkoIiIiIlvcPSx7AewFcHF2zSEiIiIqFxcSGoEnnRnhww9VdXFqrSIiIiKyxM2wdGfZECIiIqIoFbOERGSd7yfJykXkTSLyjIjsFpFrPL9vF5F7wt8/KCIbqv8KRERENNcleZbQd6y/LwRwOoBnALwi7kMi0gLgSwBeD2AAwA4RuVdVn7QWuwrAUVU9Q0S2AfgMgHdU0X4iIiKaByrOsKjqOdbPRgAXAPj3BOu+AMBuVd2jqhMA7gZwibPMJQBuC//+DQCvExFJ3nwiIiKaD5LMsJRQ1V0icn6CRVcD2Ge9HgBwYdQyqjolIsMATgFwyF5IRK4GcDUArFuXKBpVk+euf0tq607LbGwz0VzGY5IoHRUvWETkv1kvcwA2AxhMsG7fTImbdZRkGajqdgDbAWDr1q2+zCUiIiKaw5KU5u+2ftoR3NPihnZ8BgCstV6vAXAgahkRaQWwBMCRBOsmIiKieaTiDIuqXlfjuncA2CgipwPYD2AbglL/tnsRPKvofgCXAvg3VeUMChEREZWQqOsD5/lBZVT1rRVXLvJmAJ8H0ALgVlX9tIh8EsBOVb1XRBYCuB3AJgQzK9tUdU+FdQ4iqMCbhuVw7p+h1LCvs8X+zg77Olvs7+yk1dfrVbW30kJxFyyDCG6IvQvAg3DuN1HVn9ahkU1FRHaq6tZGt2M+YF9ni/2dHfZ1ttjf2Wl0X8eFhE5FUEPFPLX5OwDuUtUnsmgYERERkRF5062q5lX1X1X1CgAXAdgN4Cci8qHMWkdERESECjfdikg7gLcgmGXZAOALAL6ZfrMaZnujGzCPsK+zxf7ODvs6W+zv7DS0r+PuYbkNwCsBfA/A3ar6eJYNIyIiIjLiLlgKAEbDl/ZCAkBVdXHKbSMiIiICEHPBQkRERNQsklS6nfNE5E0i8oyI7BaRaxrdnmYmImtF5Mci8pSIPCEiHwnfXyYiPxSRX4V/Lg3fFxH5Qti3j4rIZmtdV4TL/0pErrDe3yIij4Wf+YJ5IGbUNuY6EWkRkYdE5Nvh69NF5MGwH+4Rkbbw/fbw9e7w9xusdVwbvv+MiLzRet879qO2MdeJSI+IfENEng7H+MUc2+kQkY+F55DHReQuEVnIsV0/InKriLwoIo9b7zVsLMdtIzFVndc/CIra/QeAlwBoA/AIgJc3ul3N+gNgFYDN4d+7AfwSwMsB/B8A14TvXwPgM+Hf34zgPihBkG32YPj+MgB7wj+Xhn9fGv7uFwAuDj/zPQC/Fb7v3cZc/wHw3wDcCeDb4euvIyiyCABfAfCH4d//K4CvhH/fBuCe8O8vD8d1O4DTw/HeEjf2o7Yx138QPD3+feHf2wD0cGyn0s+rATwLYJE13t7DsV3XPv51BM/+e9x6r2FjOWobVX2nRndqo3/CDv++9fpaANc2ul2z5QfAvyCo1/MMgFXhe6sAPBP+/UYAl1nLPxP+/jIAN1rv3xi+twrA09b7xeWitjGXfxA8g+s+AL8B4NvhwX4IQGv4++L4BfB9ABeHf28NlxN3TJvlosZ+3Dbm8g+AxQj+ERXnfY7t+vf1agSFSZeFY/XbAN7IsV33ft6A0guWho3lqG1U830YEpo+cIyB8D2qIJyW3YSgEvJKVT0IAOGfK8LFovo37v0Bz/uI2cZc9nkAfwygEL4+BcCQqk6Fr+3+KfZp+PvhcPlq90HcNuaylyB4Ev1XJQjB3SwineDYrjtV3Q/gswD6ARxEMFb7wLGdtkaO5Rn/W8sLFueRAyHeiVyBiHQB+EcAH1XVY3GLet7TGt6fd0TktwG8qKp99tueRbXC77gPkmlFMIX+ZVXdhCBLMu6eNvZrjcL7Gi5BEMY5DUAngN/yLMqxnY0s+nHGfc8LluAqb631eg2AAw1qy6wgIgsQXKzcoaqmkOALIrIq/P0qAC+G70f1b9z7azzvx21jrno1gLeKyHMA7kYQFvo8gB4RMUUf7f4p9mn4+yUIHipa7T44FLONuWwAwICqPhi+/gaCCxiO7fr7TQDPquqgqk4iKEj6a+DYTlsjx/KM/63lBQuwA8DG8M7xNgQ3dMU+qXo+C+8EvwXAU6r6OetX9wIwd5BfgeDeFvP+74d3iF8EYDicJvw+gDeIyNLwf1tvQBBLPghgREQuCrf1+866fNuYk1T1WlVdo6obEIzLf1PVdwL4MYBLw8Xcvjb9c2m4vIbvbwszLU4HsBHBDXPesR9+Jmobc5aqPg9gn4icFb71OgBPgmM7Df0ALhKRjrAvTF9zbKerkWM5ahvJNfqmoGb4QXD38i8R3FX+8Ua3p5l/ALwGwTTeowAeDn/ejCA2fB+AX4V/LguXFwBfCvv2MQBbrXW9F8EzqnYDuNJ6fyuAx8PPfBHT9YK825gPPwBei+ksoZcgOCnvBvAPANrD9xeGr3eHv3+J9fmPh/35DMK7+cP3vWM/ahtz/QfAeQB2huP7nxFkRnBsp9PX1wF4OuyP2xFk+nBs169/70Jwf9AkgtmNqxo5luO2kfSHheOIiIio6TEkRERERE2PFyxERETU9HjBQkRERE2PFyxERETU9HjBQkRERE2PFyxERETU9HjBQkRERE3v/wfce0upBfgw5QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f08418ccdd8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(3, figsize=(9, 5), sharex=True)\n",
"\n",
"ax1, ax2, ax3 = axes\n",
"\n",
"vals, windows, counts = allel.windowed_statistic(\n",
" pos, allel.GenotypeVector(gv).is_het(), statistic=np.sum, size=1000)\n",
" \n",
"ax1.scatter(windows.T[0], vals, marker=\"o\", s=5)\n",
"\n",
"ax2.broken_barh(tuple(zip(pdf.start.values, pdf.length.values)), \n",
" (0, 1))\n",
"ax2.set_ylabel(\"Poisson model\")\n",
"\n",
"ax3.broken_barh(tuple(zip(mdf.start.values, mdf.length.values)),\n",
" (0, 1))\n",
"ax3.set_ylabel(\"Multinomial model\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Test on some real data..."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"import sys, os"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"sys.path.insert(0, 'selection_paper/ag1000g-phase1-selection-paper/agam-report-base/src/python')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"ag1k_dir = '/kwiat/vector/ag1000g/release'\n",
"from ag1k import phase1_ar3"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"phase1_ar3.init(os.path.join(ag1k_dir, \"phase1.AR3\"))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# get a kenyan sample"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"pos = allel.SortedIndex(phase1_ar3.callset_pass[\"2L\"][\"variants/POS\"][:])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"kes_df = phase1_ar3.df_samples.query(\"population == 'KES'\")"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"304"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kes_ix = kes_df.index[5]\n",
"kes_ix"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"gv = allel.GenotypeVector(\n",
" phase1_ar3.callset_pass[\"2L\"][\"calldata/genotype\"][:, kes_ix])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((10377280, 2), (10377280,))"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gv.shape, pos.shape"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"real_accessibility = phase1_ar3.accessibility[\"2L\"][\"is_accessible\"][:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Compare to phase1 AR3.1 ROH calls"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"roh_data = \"/kwiat/vector/ag1000g/release/phase1.AR3.1/roh/main/runs_of_homozygosity.h5\""
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"import h5py"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"roh_fh = h5py.File(roh_data, \"r\")"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'AK0070-C'"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kes_sample = phase1_ar3.df_samples.loc[kes_ix]\n",
"kes_sample.ox_code"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"xx = roh_fh[\"2L\"][kes_sample.ox_code][:]\n",
"\n",
"phase1df = pd.DataFrame(xx, columns=[\"start\", \"stop\"])\n",
"phase1df[\"length\"] = phase1df.stop - phase1df.start"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"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>start</th>\n",
" <th>stop</th>\n",
" <th>length</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>35831727</td>\n",
" <td>35831726</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>37842374</td>\n",
" <td>37871903</td>\n",
" <td>29529</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>37905959</td>\n",
" <td>37917042</td>\n",
" <td>11083</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>38468052</td>\n",
" <td>49364325</td>\n",
" <td>10896273</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop length\n",
"0 1 35831727 35831726\n",
"1 37842374 37871903 29529\n",
"2 37905959 37917042 11083\n",
"3 38468052 49364325 10896273"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phase1df"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"callset = phase1_ar3.callset_pass\n",
"accessibility = phase1_ar3.accessibility"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"hetp [ 0.001 0.0025 0.01 ]\n",
"start [ 0.33333333 0.33333333 0.33333333]\n",
"tx [[ 9.99000000e-01 5.00000000e-04 5.00000000e-04]\n",
" [ 5.00000000e-04 9.99000000e-01 5.00000000e-04]\n",
" [ 5.00000000e-04 5.00000000e-04 9.99000000e-01]]\n"
]
}
],
"source": [
"df, froh = allel.roh_poissionhmm(\n",
" gv, pos, is_accessible=real_accessibility, window_size=1000)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"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>start</th>\n",
" <th>stop</th>\n",
" <th>length</th>\n",
" <th>is_marginal</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>25044</td>\n",
" <td>35832833</td>\n",
" <td>35807789</td>\n",
" <td>True</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>37842854</td>\n",
" <td>37869174</td>\n",
" <td>26320</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>37906143</td>\n",
" <td>37915810</td>\n",
" <td>9667</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>38468037</td>\n",
" <td>49355530</td>\n",
" <td>10887493</td>\n",
" <td>True</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" start stop length is_marginal\n",
"0 25044 35832833 35807789 True\n",
"1 37842854 37869174 26320 False\n",
"2 37906143 37915810 9667 False\n",
"3 38468037 49355530 10887493 True"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'Multinomial model')"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAFMCAYAAADoa+cuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4HOWV7/HvacmWsSVjGwvjC9gCAjhA2CxsCElmwhYSEkgmk4khJIQAvnMnZJu5d4asBJLMkJmETMg2cVjCFjuzZDET1iwkMIAtCbOEzTG2ZBtvsiXbWrBkqc/9o6vb3a1uqSSrqtvS7/M8/UjVXf3W6epW69Sp933L3B0RERGRcpYodQAiIiIiQ1HCIiIiImVPCYuIiIiUPSUsIiIiUvaUsIiIiEjZU8IiIiIiZU8Ji4iIiJQ9JSwigpk9ambtZlaVd/+PzeyrWcsnmtlmM/u7YHmGmf3czLrMrMXMLst7/mXB/V1m9gszm5H1WGferd/MvpP1+Llm9rKZdZvZ78xs7hCvYbaZ3RbE1xE89wYzmxLi9R9nZr80s1YzazOzh8zs+KzHv2xm9wzVjohERwmLyDhnZnXAWwEHLh5kvVOB3wFfc/dvBnd/D+gFZgEfAn5gZicG658I/BD4cPB4N/D9dHvuXp2+BY+/DvxH8NyZwM+ALwIzgEbgp4PENgN4EjgIOMvda4DzgWnAMSF2wzRgOXB8EMtK4JchniciMTHNdCsyvpnZl4B3ACuA49z93VmP/RjYSOqf+f3Ade5+a/DYFKAdOMndVwf33Q285u7Xmdk/AnXuflnw2DHAS8Ah7t6RF8MVwPXAMe7uZrYY+Ki7vzlrW9uB09z95QKv4avAJcAp7p4chX0yA9gBzHT3HWb2ZeAN7n75/rYtIiOjCouIfAS4N7i9w8xm5T2+AHgQ+Ew6WQkcB/Snk5XAs8CJwe8nBssAuPurpKoxxxWI4QrgLt93BJX/3C7g1ay2850H/GywZMXM/tvMriv2eJ63AVvcfUfI9UUkYpWlDkBESsfM3gLMBf7d3beb2avAZcC3slY7k1S14YG8p1cDu/Lu2wXUhHw8HcMc4M+Aq/Labh3quVkOATYXeQyA7MrRYMzsCFKnuv42zPoiEg9VWETGtyuAh919e7D8k+C+bN8DGoBHzGx61v2dwNS8dacCHSEfT/sI8Li7rwvTtpm9Nauj7gvBYzuA2YVe4HCYWS3wMPB9d1+6v+2JyOhRwiIyTpnZQcBfAX9mZlvMbAvwGeAUMzsla9V+Uh1q1wMPmVk6kVgNVJrZsVnrngKkk4gXguX09o4GqoLnZfsIcGfeffnPnUKq8+wL7v5YVofd9CmiXwPvM7MRf6cFydjDwHJ3/9pI2xGRaChhERm/3ksqGTkBODW4vRF4jFQSkeHue4EPkOr4er+ZTQn6lfwMuNHMppjZ2aQ6vt4dPO1e4D1BRWQKcCOpfiaZCouZvRk4nGB0UJafAyeZ2fvNbBLwJeC5Qh1uAzeTqsDcmR7+bGaHm9nNZnbyUDsiSMIeAv7H3Yv1c0mY2aSsW1WR9UQkAkpYRMavK4A73H29u29J34DvAh8ys5w+bu7eC/wFsAe4L6jQ/A2pocTbgKXA/3H3F4L1XwD+mlTiso1U/5O/KRBDThITPLcVeD/wNVIjkRYCi4q9EHdvA94M7AVWmFkH8BtS/V7WAJjZA2b2uSJNvA84A7gyb26YOVnrXEpq6HX69mqxeERk9GlYs4iIiJQ9VVhERESk7ClhERERkbKnhEVERETKnhIWERERKXtKWERERKTsHXBT88+cOdPr6upKHYaIiIiMgqampu3uXjvUegdcwlJXV0djY2OpwxAREZFRYGYtYdbTKSEREYldMum0dvSgucAkrAOuwiIiIge2ZNK59EdP0dTSzvy501l6zZkkElbqsKTMqcIiIiKRy66o7Ojqpamlnb6k09TSzo6u3lKHJwcAJSwiIhKpdEXlrH/6DYuWPMWMyROYP3c6lQlj/tzpzKyeWOoQ5QCgU0IiIhKp/IpKW/dell5zJju6eplZPREznQ6SoanCIiIikZpZPXFARSWRMGprqpSsSGiqsIiISKTMbEBFJZl0VVhkWJSwiIhI5NIVFdAoIRkZnRISEZFYaZSQjIQSFhERiVWhPi0iQ9EpIRERiVWhPi0iQ1HCIiIiscvu0yIShk4JiYhIyenaQjIUVVhERKSkNGpIwlCFRUREhmWoashwqyUaNSRhqMIiIiKhDVUNGUm1JD1qKP0cjRqSQpSwiIhIaIWqIdMPmsCa1k6Om1XNjq69NDa30e/Q2NzGjq7enAnj0iOD3MkZJaRRQzIUJSwiIhJafjVk2qRKTvvqI3Ts6aNmUiWNnz2XyVWVdOzpY3JVJTMmTwByKy+nz5kGGE+vz63CaNSQDEYJi4iIhJZfDXllSwcde/oA6NjTx9MbdtLd2w9Ad28/bd17qa2pGlCZwYz+rCqNkhUZijrdiojIsGRfafm4WdXUTEod+9ZMqmTBUdOpD2axrc/qj5I/u61mupXhsgNtzHt9fb03NjaWOgwREQn09SUzfVgSiUSmr8qMyRNo695b8Of2rl4MqK2pGtCfRcYXM2ty9/qh1tMpIRER2S+VlQnmzZ6aWU4kjEOmTOTSHz1FY3Mbk6sq6e7tp37udO69aiGX3boi0wfm3qsW8qHbVmgOFhmSTgmJiMh+Syadrbv3sG33Hvr7k7y8ZTeNLe30e6pvS7q/yprWzpy+LPnLmoNFilGFRURE9ksy6Sxa8iQrm9sBqJ5YQWfQ8RagpqqS7r39zJ87neNmVeeMMspfVn8WKUYJi4iI7Jf0CKC07GQlYfDI376NikQiM//KLZeelum/kh51tLVjDzu7enF39WORgpSwiIjIfkmPAEpXWLLNnzudWVMnYWYFZ8E1S1VoLvjWHzJzuaz6wvlUVqrHguTSJ0JERIalry/Jy5t3k0wmSSad7Z29LL3mTFZ87lzu/FjuYI9PnXNM5vdi1wxa09qZM5fLn7Z1xPdi5IChCouIiITW15fMzGxbXVXBCbOn8vT6ncyfO507rziDK27PnXbi8tsbWVA3nWWLzyp6zaA31E4hASSD51y//AWWLT5Lo4UkhxIWEREJLbsa0tnTT0NzOw6sXNfGL599reBzGpvbWb21g+MPqynYX2VH195MsgLQ2LJvffVnkTRNHCciIqElk0lOuTFVYTEgzH+Q9Cih+rnTufvKBcz/x19n+qs0fe48Lrv1SRrX78qsX11VyevB+pqXZewLO3FcLH1YzOxIM/udmb1kZi+Y2aeC+2eY2SNm9qfg5/Q44hERkZFJJBI0fe48vn/ZaYQpfnz1kjfS3dtPf9JpbGmnYX1bTn+VpvXtPJ2VrAC8HqyveVkkW1ydbvuAv3P3NwJnAh83sxOA64DfuPuxwG+CZRERKVPJpHP57Sv4m5+sIhmivPKP979Cf1DJP6gywbcefiXzWE1V6tpDpx05Nec5p8+ZpusMyQCx9GFx983A5uD3DjN7CTgcuAT482C1O4FHgX+IIyYRERm+HV29A4YvZ3eYzde9d98jnb39rNq4GwADHvnMW6moqOAr7zuZd93yeGa9Gy45kdqaSbq2kOSIfVizmdUBpwErgFlBMpNOag6NOx4REQlv6sSKnOX5Rx7MQRPD/yvxoCzjwCd/+izJpDPvsJqcKz7PO6wmM6mcSFqsCYuZVQP/BXza3XcP43mLzazRzBpbW1ujC1BERAb19IadOctXv/Vountz6yuDpRlJUrPfAjSua+OVLbsxM1Z94Xwe/NRbWfWF89ja0cNLm3aRTBar28h4FNuwZjObQCpZudfdfxbcvdXMZrv7ZjObDWwr9Fx3XwIsgdQooVgCFhGRAU49PLe/ydlHTx8wUmioL+l035ck8M5bHmdB3QyWLT6T42bV8MEfPkFDSyopqqmqZNUXNeutpMQ1SsiA24CX3P3mrIeWA1cEv18B/DKOeEREZGSaNuaO6FnWWHjuleFoaG7jhU07eXzNNhpb9lVwOnr6WNPaud/ty9gQyzwsZvYW4DHgefb1zfocqX4s/w7MAdYDH3D3tsHa0jwsIiKl0deX5KTrH2BP/9Drjoaaqkqevf58EglVWMaysPOwxDVK6HGKn9Y8N44YRERk/6xp7YwtWQH49d++TcmKZOiTICIioRw3q5qDYur5WF1Vyczqqng2JgcEJSwiIhJKIpHgP/732bFsq7Onj7buvbFsSw4MSlhERCSUZNL5/M9Wxba9/DlfZHxTwiIiIqHs6Orl2c3dsW0vf84XGd+UsIiISCgzJk+IdXvzDp0c6/akvClhERGRUOLuU3LnUxti3Z6UNyUsIiISyozJEwaddn+0XXXWkTFuTcqdEhYREQmlrXvvkNPuj6bnt3TFuDUpd0pYREQklLgrLN98ZDXJpC4fJylKWEREJJS4KyyrNuxkR1dvjFuUcqaERUREQol7XpTTjpzGzOqJsW5TypcSFhERCSXueVE+9uY6zOI8CSXlTAmLiIiEUj93Wqzbu+OJZvVhkQwlLCIiEsrOOC/VDDStL+8+LMmk09rRg7uSqjgoYRERkVBmVk/kpEPju4KyGUybFNPloYcpmXQu/dFTnPVPv2HRkqdUCYqBEhYRkRI7UI7UzYyL5x8V2/aSDq9uL8+5WHZ09dLU0k5f0mlqaS/rStBYoYRFRKSEDqQj9WTS+WXj2li3efQh5Xk9oZnVE5k/dzqVCWP+3OkazRSD8qy1iYiME4WO1Gtr4jvtMhw7unp5oTXeSsLaHd3Mmz011m2GYWYsveZMdnT1MrN6okYzxUAVFhGREjqQjtTj7k+SMHhD7ZRYtzkciYRRW1OlZCUmqrCIlJFk0sf9EVvU+2C029/f9g6kI/W4+5MkHdpf76O2JtoJ64bzHupvtHSUsIiUiXRfhqaWdubPnc7Sa84kkRhfX4hR74PRbn+02ksfqZe7uPuTnH7kwZFXnIbzHupvtLR0SkikTJT7qIM4RrJEvQ+KtT/Uayv2eLm/Z6Nt7Y7uWLf3zhMPi7yKMZz3cLy93+VGCYtImSjnvgxxjWSJeh8Uan+o1zbY4+X8nkUh7v4kl9YfHvk2hvMejrf3u9zolJCMa319Sda0dnLcrGoSidLm73H0ZRjp+fe4RrJEvQ8Ktb+9s2fAaztkysTMOoO99pHGG3U/mqj6WbS/3jdqbYXx2NqdXPimSYO+htHqQ9Ta2cNQzw77fqufSzSUsMi41deX5LSvPkLHnj5qJlWy6gvnU1lZ2qQlyr4M+3P+PX1kmX5ulEeWUffnyG8//7XNmDwhZz/95OqFg7724cYbdT+ae69ayIduWxFJP4uZ1RM5/Ygant7YMSrtDeXaZauof7L4axjNffnJpatCtTPU+61+LtFRwiJDivtoYSTbG8lz1rR20rEndcTYsaePNa2dBed7SLc9Y/IE2rr3Rj66JIr93deXpKG5jYbmNpIOjc1tQ1ZJ8uMIexQ6VDvZ943GPk0mPRPXIVMmDmgve1vbu3oxyBmKmn/UvL0zt6LS1r236FH1SN6r7IpNY3Mbq7d2cPxhNUXb7e/3TBUQbMD28itAa1o7h2x/KMUqj2bGdy+r583//LvQbe2P/uA1vLxlNzOrqwYMId7R1Utjcxv9RT7TYd6fZNJZvbWDxpZ2+oN92NrZQ8JsRJ/L/Pdja8cednXvLYsq7oFOCYsMKu6jhZFsb6QxHjermppJlZkKS+ofQuG2G5vbmFxVSXdvP/URji6JYn9nV5LSJldVMmPyhGHFB+GPQodqZ7T2aTLpLFryJCub2wGoqaqke+++9rK3ddDECjp7UhfvW1A3g2WL920v+6i5UDXJbOBR9Ujfq3T76dd/0Xcez3n92e2eduTBvLylk46ePqqrKjhh9lSeXr8zZ3v58R43q3rQ9ocyWOWxry/Jed96NNybM0oOmlDBu255HBj4vs2YPIHJVam/4fzPdJj3p9Df9+lzpvGJpat4eoR/g9nvx2lHHswF3/pDWVVxD2RKWMpAHBWFsEfz+cvZRwsr16WOdN44eyruFN3+YLFlP1aojeH2lcg/OmpY18aKtTtYePSMIY9mEokEq75w/qB9WNLx9DuZf/jZR3J9fUlWb+vgkCkTOXTqpMw+TB/x19ZUDXidhY7o0v0mXt6ye99ryTs6Hmk/hexKUlp3bz/bu3oLHkUW2qdPvrqdaZMnFIy5tTM1eiYR/FPPfr07unozz2kMjlzb0/dl7dMVwTbecGg1h06dlGkjvwKTfs3TJlXy6vYupk+eQFNLeyb2jp5Ue+mj5B2dPZltpZMVgMaW1L499tBq2rr3Zqove/v7Wb+jm3s+dgY79/QzbVIlr2zp4A21U2h/vS8nntaOnqIVq/z3pq8vyctbd5MwY95hNSy95kxWb+3gXbc8lnl++kg8/Zr6kk5jy07SXXw7e/ozn8WVefsru/qVrhit3trBRd95PPMevrJlN8fNqmF7V2/m/TpkykS2d/XSn0zS3tXLzOoqdnT15lQen1jbyvTJqdex6/W9dO+N99IB3Xv3vW8NzW1s7+xhZnUVWzv2sK61k+7e1OPZn+kZkyfwp22dg1ZNcj7nDt09fdz3ibfgwHu+8/iglci+viR/3LyTjW2vc/qcg0kkKjKf/+xq5Oqtu/nwbQ1Aal8+8uIW3nHSYZnvmuF+j6crX+nP42DfscUqimG3W479cKzcL7aVr76+3hsbG0sdxqiJo6IQ9mi+0LlvM/irf3uChpadmfbOmDsNswRPrx+4/cFiy37s9DnTABvQhruzaMm+5y9bfOagf1DZR0ede/pIGPQ7o3Y0k46nYV0bltX2M188H3c49SsPZx21T+cnV5/JZbc+lTniP2PudMz2vc70Ps6vLvzk6oVc+qN9zwOoSP33ob7AezOcfgrJZJKTb3g45x92AphfN51VeUfrhfapAcm8NhfUzWDpNQNjzn+993xsAad/7dd07NlXIWhsbs/sy0LSn6+mltx9lH7N2e9FdVUF82ZV07h+V84+mz9nGg40ZMWWrWZSJV09fUypSv2cnFV9Sbfz/BcvYOE//5aOPX1UGDgwJYhn/pxpJJNO4/rU30VNVSXPfOl8KioSA/4G7r5yAad99WG6epOZdVd98XzM4NSvpCoZUyYmSFgip5LS1NLOxAoYqp9r+nNX6G/3gz98Muf9qZ5YQWdvf9FlgCkTEnTtzX/HS+OMuan3sTHv+8edzL5Py95vk6v2vb/p9wuzTNWk0N9hoc9M9vua1teX5E1ffqDg+5KuAAF88Ie535vZba764vkkEjas7/Hsylf2d0Oh79jBKorDqTzFVVk3syZ3rx9qPdWm8hSab2Ek808Ua2fr7j1s270nc392RaFhXRsvbd7Fpp3dvLRpF8lksmBb6aPWvuCoNX8ugL6+JC9v3p15fv76Wzv28OKmXZmj+XQ1Y/W2jpzlHV299PT0c868Q3Pab2jZSdP63PXSMWaPuGhsbmNbxx5aO3ro709mjmZSj7XT2NKWWS/dxvbOXu766Bnc/bEF/OTqBTn9EDa0d/HAc5vo7d3L5l2v88Srraxc15Y5Uv+n952U+SeY7pOSjqu/P1n0/diy63U2tnfxxJpW+vv7c/Y5wL1XLeTuqxekrnUPdPX08crWDp5q3pHzT25lczvPb2rP+QfR0NKesz8amttyqgs/+shpXHHWXFo7Xs+pFEDqH3Lm6Hjrvvdq5bo2Hnlxa2a5ITjHn/25yX5t2zp6uOuqM3I/n8DTLTsz7bV27KGvL8mKtTtoaN63T7/23hMHJCsAV5w1hz9uah+QEDS2tGc+GyvWtfHAi5vpCo7Wu3r6aWxpJ0nxZCXVRurzlY6hP+/zmf38zp5+LnzTrPRbQ7/D/z3naM6bN5PGIsnK3VedQVdPH8mg/WRe9SXdzr899mqm0tDvZNbvTzorm9tz/mF29PTx38+9xmN/2srLW3Zlfc7beGJdayZZSa/7xNpWHl+zja6gItTVm8xUhzp7+jnnuEO448r57AkxKGdlczv3Pb+Rhua2zPv5ypbdANyy6LScPkf5yUn+MjCsZOW+a98c2T+Rz77jWL696BS+fPGJOa+haf3OAckK7Pt89We9rx17+vjF35zFteccQ1Pzvr+fFc07cr47/vn9J3LuvNoBf4NdvX2Z76bXdnbz+J+28cTabUWTyJXNbfzPmlb+8KetBZMVSL3/T63bzkt537+tnT0F/8+kv4v+tK0j5/OYXelM27Z7T+Z15VcU0+sNNpdMoe/xcppvRhWWLIOdbx+N8/bZ59rTGW+ho6C0mqpKmj5/Hh++Y2VOW+6eOTJLH+2njwAKnX/OP5IzLPNFlT7nn390saBuBndeUc8bb3i44GvM7ivwk6sXctmtKzKVE/fUP2oY/Ei2Iqti8fTnz+Py21cOqGSs+kLqSOQDP3icpg27Q7+XqaOY8/jQbSsL9pXIfz/SKgxeuP4dXHFnQ8HXA4WPSsPG1NlTuGIBMP/Ig2nasKvgcw0I+5daPbGCN86uKfqFWciUCQkSiUTmn+ZIVU+sYN5h+yoeI3VQpfF6375XXD9nGljukfZITaqAPcN/+4ZlysQKukbwGRlNZ8ydTl9fH6tei2ZET+o7bCGLljxV8PsrKpMSsKfAH1D+d9toOaNuOsn+ZNG/zVHZxtzpWMIG9Jsp1p8pu8KSrkL39SU5+YYHC56yy64UFatiD1YBH6zSPRrCVlhK3ofFzC4Evg1UALe6+02liqVQ7/0ZUyYW7IWePr83dWIFjevbecOhUzh06kH09/uA0RjZ59LTGltSR8Uzpkzk+otP5KKgU1m21FHbBhrWtZEkdZ5/084u1u3ozjnH/OAfN1NfN51Dpx40YOTLPSvWUpGwzH3ZR3rpbVz7Z3XU103nqrueydy/YO7BfOfRl4vuq46ePm754JtYePRMdnT1snJdG07qaO8rl8yjKTgSzo4FcrP+7GrI9x9dzYp1bak7su6//bGXea29Z1jJCsAn3z6Hmx96MdNmdl+Jx9ZsozKRKHi6oN/hy8ufzTxvZXP7gFExI/1H1NnTh1M88ThjTk3RL8XhHFZ09vYP+x976qh6/08DdPb28843HbrfCUt2sgJQndjLo82jcx2bqJMVGPlnZCiffOsRfOexjaE+Dw0t0SYR33z/CSQSCZYtPouXt+zOdIwt5J6rzuCKOxrpT+6LPAHc+N4T+MIvXhzWdgslKzDwu220FDutOJree9psvrT8JfqDyuRPVjRz7GGT+Z/V7Zn/JQ0tO/nFx8+irauHzj19HDVzClMPmsCDz2/m7cfN5LertxftX9TR08eyhhbOnXcoh0yZxBff/UYAaqur6OzuZVnDBhYePS3ne3z5tW9m9+t7OaZ2Cq0dPWVxkceSVljMrAJYDZwPbAQagEvdvegnOMoKSzrzzOkxnpXVpqsZZqlzj+k3N23+kQezeltXzlFqdVUFbzxs4NFugtH495DrjLnTueMj8znpK78e5ZYHd+r/msIzm+K9KJqIlN5zXziPqdVV9PcnOfmGhwomDQuOmsGyaxZy6Y8G9hnp3buXZ17rzKw7nEqijMxI93H+CK3RdKBUWBYAa9x9LYCZLQMuAYaXco+C3t5+Hnt1G2fPncTW9gk079wLkJNodOzp47jPP8Cx0+HlAkl3oaPjzp7+gqX5KI4FGlraWdq4MYKWB6dkRWR8uv2JZj59wfG0de/l9by+Lwbc/8m3MG/21Jy5btIjrdydBf/4m9IEPo6NNCFc2dzG92+/n2uvvmhU4xmOUicshwMbspY3AgvjDqK3t5/jvvRgqHX7KZyslIuvPVD8NI6IyGg65rDUqa+Z1ROpnzs9d8TYUTMyyQrkznVTW1PF1vbOAe3NnzOtYIdaKQ/fWAMb7/kVN11emqSl1AlLodrSgATQzBYDiwHmzJkz6kE0rS/jDEREpEx9+ifreM/JJ2BmLFt81oA5eQbr83DDrwYeXN343pMG7Qsjpbfsj1CqjqalHta8ETgya/kIYFP+Su6+xN3r3b2+trZ21INYcNT0UW+zFCZVlDoCERlPHr5mXub3RMKYNXUShx18UGYSxcH8y/vemLNcXVXBvMNq+OWlR0QSq4yOr59Wum2XOmFpAI41s6PMbCKwCFgedxAVFRWsvvFCbrvidP7u7Yczr7aKv3/7LI6ZnuCth+1b7221MAE46RDjhLzJVw9LwJxqOG7SvvtmxRJ9yi3nTOaFG9/Jis+dy82XHxvjluGmCw6KdXsSjcMrU0P1ytm0UgcwTJMr4PPvmFzqMCLxm2vmccwxx4z4+VOmTGHldX/OhSfU8vPF9Tx3/QUkEglOOeUUbvxLJS3l6OunwQc/WLo+LCWfh8XM3gX8K6nvytvd/WuDrT/WZroVEREZz8KOEip5wjJcZtYKtETU/Exge0Rty0Da3/HTPo+X9nf8tM/jNRr7e667D9nf44BLWKJkZo1hsjwZHdrf8dM+j5f2d/y0z+MV5/4udR8WERERkSEpYREREZGyp4Ql15JSBzDOaH/HT/s8Xtrf8dM+j1ds+1t9WERERKTsqcIiIiIiZU8Ji4iIiJQ9JSwiIiJS9pSwiIiISNlTwiIiIiJlTwmLiIiIlD0lLCIiIlL2lLCIiIhI2VPCIiIiImVPCYuIiIiUvcqoGjaz24F3A9vc/aQCjxvwbeBdQDfwUXd/eqh2Z86c6XV1daMcrYiIiJRCU1PTdnevHWq9yBIW4MfAd4G7ijz+TuDY4LYQ+EHwc1B1dXU0NjaOUogiIiJSSmbWEma9yE4JufsfgLZBVrkEuMtTngKmmdnsqOIRERGRA1cp+7AcDmzIWt4Y3CciIiKSI8pTQkOxAvd5wRXNFgOLAebMmRNZQHXX/SqytkVExormmy4atbbqrvtVwfb0fVxeRvM9H6lSVlg2AkdmLR8BbCq0orsvcfd6d6+vrR2yX46IiIiMMUUTFjPrMLPdwa0ja7nDzHaPwraXAx+xlDOBXe6+eRTaFRERkTGm6Ckhd6/Zn4bNbCnw58BMM9sIXA9MCNr+N+B+UkOa15Aa1nzl/mxPRERExq5QfVjM7C3Ase5+h5nNBGrcfd1gz3H3S4d43IGciu9yAAAVT0lEQVSPh45URERExq0h+7CY2fXAPwCfDe6aCNwTZVAiIiIi2cJ0un0fcDHQBeDum4D9Ol0kIiIiMhxhEpbe4PSNA5jZlGhDEhEREckVJmH5dzP7IamZaK8Bfg38KNqwRERERPYZstOtu3/DzM4HdgPHA19y90cij0xEREQkEGqUUJCgKEkRERGRkiiasJhZB0Wmygdw96mRRCQiIiKSZ8iJ48zsRmALcDep6/98CI0SEhERkRiF6XT7Dnf/vrt3uPtud/8B8P6oAxMRERFJC5Ow9JvZh8yswswSZvYhoD/qwERERETSwiQslwF/BWwFtgEfCO4TERERiUWYYc3NwCXRhyIiIiJSWJhrCR1hZj83s21mttXM/svMjogjOBEREREId0roDmA58L+Aw4H7gvtEREREYhEmYal19zvcvS+4/RiojTguERERkYwwCct2M7s8GCVUYWaXAzuiDkxEREQkLUzC8jFSo4S2AJuBvwzuExEREYlFmFFC64GLY4hFREREpKAhExYzOwr4BFCXvb67K4kRERGRWIS5WvMvgNtIjQ5KRhuOiIiIyEBhEpY97n5L5JGIiIiIFBEmYfm2mV0PPAz0pO9096cji0pEREQkS5iE5U3Ah4Fz2HdKyINlERERkciFSVjeBxzt7r1RByMiIiJSSJh5WJ4FpkUdiIiIiEgxYSoss4CXzayB3D4sGtYsIiIisQiTsFwfeRQiIiIigwgz0+3vR9q4mV0IfBuoAG5195vyHv8o8C/Aa8Fd33X3W0e6PRERERmbwlRYRsTMKoDvAecDG4EGM1vu7i/mrfpTd782qjhERETkwBem0+1ILQDWuPvaYITRMuCSCLcnIiIiY1SUCcvhwIas5Y3Bffneb2bPmdl/mtmREcYjIiIiB6ghExYzO9vMHjGz1Wa21szWmdnaEG1bgfs8b/k+oM7dTwZ+DdxZJIbFZtZoZo2tra0hNi0iIiJjSZg+LLcBnwGagP5htL0RyK6YHAFsyl7B3XdkLf4I+Hqhhtx9CbAEoL6+Pj/pERERkTEuTMKyy90fGEHbDcCxZnYUqVFAi4DLslcws9nuvjlYvBh4aQTbERERkTEuTMLyOzP7F+BnDOPih+7eZ2bXAg+RGtZ8u7u/YGY3Ao3uvhz4pJldDPQBbcBHR/YyREREZCwLk7AsDH7WZ90X6uKH7n4/cH/efV/K+v2zwGdDxCAiIiLjWJiJ494eRyAiIiIixYQZJXSwmd2cHqVjZt80s4PjCE5EREQEws3DcjvQAfxVcNsN3BFlUCIiIiLZwvRhOcbd35+1fIOZPRNVQCIiIiL5wlRYXjezt6QXzOxs4PXoQhIRERHJFabC8n+AO4N+K4aGH4uIiEjMwowSegY4xcymBsu7I49KREREJEuYUUKfCpKVDuBmM3vazC6IPjQRERGRlDB9WD4WVFUuAA4FrgRuijQqERERkSxhEpb0VZffBdzh7s9S+ErMIiIiIpEIk7A0mdnDpBKWh8ysBkhGG5aIiIjIPmFGCV0FnAqsdfduM5tB6rSQiIiISCzCVFjOAl5x951mdjnwBWBXtGGJiIiI7BMmYfkB0G1mpwB/D7QAd0UalYiIiEiWMAlLn7s7cAnwbXf/NlATbVgiIiIi+4Tpw9JhZp8FLgfeZmYVwIRowxIRERHZJ0yF5YNAD3CVu28BDgf+JdKoRERERLKEmZp/C3Bz1vJ61IdFREREYlS0wmJmjwc/O8xsd/7P+EIUERGR8a5ohcXd3xL8VAdbERERKakwnW4JhjS/NVj8g7s/F11IIiIiIrlCXa0ZuJfUhQ8PBe41s09EHZiIiIhIWtip+Re6exeAmX0deBL4TpSBiYiIiKSFvVpzf9ZyP7pas4iIiMQoTIXlDmCFmf08WH4vcFt0IYmIiIjkCjMPy81m9ijwFlKVlSvdfVXUgYmIiIikFU1YzGwS8NfAG4Dnge+7e19cgYmIiIikDdaH5U6gnlSy8k7gG8Nt3MwuNLNXzGyNmV1X4PEqM/tp8PgKM6sb7jZERERk7BvslNAJ7v4mADO7DVg5nIaDiyR+Dzgf2Ag0mNlyd38xa7WrgHZ3f4OZLQK+TuraRSIiIiIZg1VY9qZ/GeGpoAXAGndf6+69wDLgkrx1LiFVyQH4T+BcM9MIJBEREckxWIXllKxrBhlwULBsgLv71CHaPhzYkLW8EVhYbB137zOzXcAhwPaQ8YuIiMg4MNi1hCr2s+1ClRIfwTqY2WJgMcCcOXP2M6zimm+6KLK2RURkoGLfu/o+lnxhJo4bqY3AkVnLRwCbiq1jZpXAwUBbfkPuvsTd6929vra2NqJwRUREpFxFmbA0AMea2VFmNhFYBCzPW2c5cEXw+18Cv3X3ARUWERERGd9CXa15JII+KdcCDwEVwO3u/oKZ3Qg0uvtyUjPm3m1ma0hVVhZFFY+IiIgcuOxAK2iYWSvQElHzM1GH3zhpf8dP+zxe2t/x0z6P12js77nuPmR/jwMuYYmSmTW6e32p4xgvtL/jp30eL+3v+GmfxyvO/R1lHxYRERGRUaGERURERMqeEpZcS0odwDij/R0/7fN4aX/HT/s8XrHtb/VhERERkbKnCouIiIiUPSUsIiIiUvaUsIiIiEjZU8IiIiIiZU8Ji4iIiJS9otcSMrP7gKJDiNz94kgiEhEREckz2MUPvxFbFCIiIiKDCDUPi5kdBMxx91eiD0lEREQk15B9WMzsPcAzwIPB8qlmtjzqwERERETSwnS6/TKwANgJ4O7PAHXRhSQiIiKSa7A+LGl97r7LzIbVsJndDrwb2ObuJxV43IBvA+8CuoGPuvvTQ7U7c+ZMr6urG1YsIiIiUp6ampq2u3vtUOuFSVj+aGaXARVmdizwSeCJEM/7MfBd4K4ij78TODa4LQR+EPwcVF1dHY2NjSE2LyIiIuXOzFrCrBfmlNAngBOBHmApsBv49FBPcvc/AG2DrHIJcJenPAVMM7PZIeIRERGRcWbICou7dwOfD26j6XBgQ9byxuC+zfkrmtliYDHAnDlzRjmMfequ+1VkbYuIjBXNN100am3VXfergu3p+7i8jOZ7PlKlnDiuUKeYgttz9yXAEoD6+vqhx2GLiIjImBJm4ri/AA4D7gmWLwWaR2HbG4Ejs5aPADaNQrsiIiIyxhRNWNz99wBm9hV3f1vWQ/eZ2R9GYdvLgWvNbBmpzra73H3A6SARERGRMKOEas3saHdfC2BmRwFDDj8ys6XAnwMzzWwjcD0wAcDd/w24n9SQ5jWkhjVfOZIXICIiImNfmITlM8CjZrY2WK4j6AA7GHe/dIjHHfh4iO2LiIjIOBdmlNCDwfwr84K7Xnb3nmjDEhEREdlnyITFzCYA/xtI92N51Mx+6O57I41MREREJBDmlNAPSPU9+X6w/OHgvqujCkpEREQkW5iE5Qx3PyVr+bdm9mxUAYmIiIjkCzM1f7+ZHZNeMLOjgf7oQhIRERHJFabC8v+A3wWjhAyYi4Ygi4iISIzCjBL6TTBK6HhSCYtGCYmIiEiswowSqgDeQWr+lUrgXDPD3W+OODYRERERINwpofuAPcDzQDLacEREREQGCpOwHOHuJ0ceiYiIiEgRYUYJPWBmF0QeiYiIiEgRYSosTwE/N7MEsJdUx1t396mRRiYiIiISCJOwfBM4C3g+uGChiIiISKzCnBL6E/BHJSsiIiJSKmEqLJtJXfDwASAz/4qGNYuIiEhcwiQs64LbxOAmIiIiEqswM93eEEcgIiIiIsWE6cMiIiIiUlJKWERERKTsKWERERGRsle0D4uZfQcoOpTZ3T8ZSUQiIiIieQbrdNsYWxQiIiIigyiasLj7nXEGIiIiIlLMkMOazawW+AfgBGBS+n53PyfCuEREREQywnS6vRd4CTgKuAFoBhoijElEREQkR5iE5RB3vw3Y6+6/d/ePAWdGHJeIiIhIRpip+fcGPzeb2UXAJuCI6EISERERyRWmwvJVMzsY+Dvg/wK3Ap8J07iZXWhmr5jZGjO7rsDjHzWzVjN7JrhdPazoRUREZFwIcy2h/w5+3QW8PWzDZlYBfA84H9gINJjZcnd/MW/Vn7r7tWHbFRERkfFnsInj/t7d/7nYBHIhJo5bAKxx97VBe8uAS4D8hEVERERkUINVWF4Kfo50ArnDgQ1ZyxuBhQXWe7+ZvQ1YDXzG3TcUWEdERETGscEmjrsv+DnSCeSsULN5y/cBS929x8z+GrgTGDC/i5ktBhYDzJkzZ4ThiIiIyIFqyE63ZlZvZj83s6fN7Ln0LUTbG4Ejs5aPIDXCKMPdd7h7T7D4I2B+oYbcfYm717t7fW1tbYhNi4iIyFgSZljzvcD/A54HksNouwE41syOAl4DFgGXZa9gZrPdfXOweDH7TkOJiIiIZIRJWFrdfflwG3b3PjO7FngIqABud/cXzOxGoDFo85NmdjHQB7QBHx3udkRERGTsC5OwXG9mtwK/AdKnb3D3nw31RHe/H7g/774vZf3+WeCzoaMVERGRcSlMwnIlMA+YwL5TQg4MmbCIiIiIjIYwCcsp7v6myCMRERERKSLM1PxPmdkJkUciIiIiUkSYCstbgCvMbB2pPiwGuLufHGlkIiIiIoEwCcuFkUchIiIiMoghTwm5ewswDXhPcJsW3CciIiISizAz3X6K1ORxhwa3e8zsE1EHJiIiIpIW5pTQVcBCd+8CMLOvA08C34kyMBEREZG0MKOEDOjPWu6n8IUNRURERCIRpsJyB7DCzH4eLL8XuC26kERERERyDZmwuPvNZvZ74GxSlZUr3X1V5JGJiIiIBMJUWACeATan1zezOe6+PrKoRERERLIMmbAEI4KuB7ayr/+KA5o4TkRERGIRpsLyKeB4d98RdTAiIiIihYQZJbQB2BV1ICIiIiLFhKmwrAUeNbNfkbqWEJDqjBtZVCIiIiJZwiQs64PbxOAmIiIiEqsww5pviCMQERERkWKKJixm9q/u/mkzu4/UqKAc7n5xpJGJiIiIBAarsNwd/PxGHIGIiIiIFFM0YXH3puDXU93929mPBVdw/n2UgYmIiIikhRnWfEWB+z46ynGIiIiIFDVYH5ZLgcuAo8xsedZDNYAmkRMREZHYDNaH5QlS1w+aCXwz6/4O4LkogxIRERHJNlgflhagBTgrvnBEREREBhrslFAHBYYzE1z80N2nRhaViIiISJbBKiw1cQYiIiIiUsyQo4TMbE6hW5jGzexCM3vFzNaY2XUFHq8ys58Gj68ws7rhvwQREREZ68JcS+hXWb9PAo4CXgFOHOxJZlYBfA84H9gINJjZcnd/MWu1q4B2d3+DmS0Cvg58cBjxi4iIyDgwZIXF3d+UdTsWWAA8HqLtBcAad1/r7r3AMuCSvHUuAe4Mfv9P4Fwzs/Dhi4iIyHgQZuK4HO7+NHBGiFUPBzZkLW8M7iu4jrv3AbuAQ4Ybk4iIiIxtQ54SMrO/zVpMAKcDrSHaLlQpyR91FGYdzGwxsBhgzpxQ3WdGpPmmiyJrW0REBir2vavvY8kXpsJSk3WrItWnJf/UTiEbgSOzlo8ANhVbx8wqgYOBtvyG3H2Ju9e7e31tbW2ITYuIiMhYMmSFxd1vGGHbDcCxZnYU8BqwiNRU/9mWk7pW0ZPAXwK/dfdCc7+IiIjIODbYxHHLiz0G4O4XD/F4n5ldCzwEVAC3u/sLZnYj0Ojuy4HbgLvNbA2pysqi4b4AERERGfusWEHDzFpJdYhdCqwgr7+Ju/8+8uiKx9USUfMzge0RtS0DaX/HT/s8Xtrf8dM+j9do7O+57j5kf4/BEpYKUnOoXAqcTKrvylJ3f2E/AytbZtbo7vWljmO80P6On/Z5vLS/46d9Hq8493fRTrfu3u/uD7r7FcCZwBrgUTP7RByBiYiIiKQN2unWzKqAi0hVWeqAW4CfRR+WiIiIyD6Ddbq9EzgJeAC4wd3/GFtUpbOk1AGMM9rf8dM+j5f2d/y0z+MV2/4erA9LEugKFrNXMsDdfWrEsYmIiIgAgyQsIiIiIuVi2NcSGovM7EIze8XM1pjZdaWOZ6wzs9vNbJuZjYfTjCVnZkea2e/M7CUze8HMPlXqmMY6M5tkZivN7Nlgn490Ak4ZBjOrMLNVZvbfpY5lPDCzZjN73syeMbPGyLc33isswfDt1aSGcG8kNUPvpe7+YkkDG8PM7G1AJ3CXu59U6njGOjObDcx296fNrAZoAt6rz3h0gqvOT3H3TjObQOoK959y96dKHNqYFlz7rh6Y6u7vLnU8Y52ZNQP17h7LvDeqsMACYI27r3X3XmAZ4a6VJCPk7n+gwDWjJBruvjm4yjru3gG8xMArp8so8pTOYHFCcBvfR4cRM7MjSI1qvbXUsUg0lLCkvrg3ZC1vRF/mMkaZWR1wGqnZqyVCwemJZ4BtwCPurn0erX8F/h5IljqQccSBh82sycwWR70xJSx5lxwI6EhIxhwzqwb+C/i0u+8udTxjXTD55qmkrlS/wMx0+jMiZvZuYJu7N5U6lnHmbHc/HXgn8PHgdH9klLCkKipHZi0fAWwqUSwikQj6UfwXcK+7a/LHGLn7TuBR4MIShzKWnQ1cHPSpWAacY2b3lDaksc/dNwU/twE/J9XFIjJKWFKdbI81s6PMbCKpK0YPeqVqkQNJ0AH0NuAld7+51PGMB2ZWa2bTgt8PAs4DXi5tVGOXu3/W3Y9w9zpS3+G/dffLSxzWmGZmU4JO/JjZFOACINKRn+M+YXH3PuBa4CFSnRH/fSxf4LEcmNlS4EngeDPbaGZXlTqmMe5s4MOkjjqfCW7vKnVQY9xs4Hdm9hypg6JH3F1DbWUsmQU8bmbPAiuBX7n7g1FucNwPaxYREZHyN+4rLCIiIlL+lLCIiIhI2VPCIiIiImVPCYuIiIiUPSUsIiIiMiLDuZitmX0ra6TiajPbOaxtaZSQiIiIjMRIL2ZrZp8ATnP3j4V9jiosIiIiMiKFLmZrZseY2YPBNYYeM7N5BZ56KbB0ONuq3I84RURERPItAf7a3f9kZguB7wPnpB80s7nAUcBvh9OoEhYREREZFcFFVt8M/EfqqiAAVOWttgj4T3fvH07bSlhERERktCSAncGVyotZBHx8JA2LiIiI7Dd33w2sM7MPQOriq2Z2SvpxMzsemE7qenLDooRFRERERqTIxWw/BFwVXBjxBeCSrKdcCizzEQxR1rBmERERKXuqsIiIiEjZU8IiIiIiZU8Ji4iIiJQ9JSwiIiJS9pSwiIiISNlTwiIiIiJlTwmLiIiIlD0lLCIiIlL2/j9GvOBxJ1Q02QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f08216fdc88>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(3, figsize=(9, 5), sharex=True)\n",
"\n",
"ax1, ax2, ax3 = axes\n",
"\n",
"eqw = allel.equally_accessible_windows(real_accessibility, size=1000)\n",
"\n",
"vals, windows, counts = allel.windowed_statistic(\n",
" pos, \n",
" allel.GenotypeVector(gv).is_het(), \n",
" statistic=np.sum, \n",
" windows=eqw)\n",
" \n",
"ax1.scatter(windows.T[0], vals, marker=\"o\", s=5)\n",
"ax1.set_title(kes_sample.ox_code + \": 2L\")\n",
"\n",
"ax2.broken_barh(tuple(zip(phase1df.start.values, phase1df.length.values)), \n",
" (0, 1))\n",
"ax2.set_ylabel(\"Poisson model\")\n",
"\n",
"ax3.broken_barh(tuple(zip(df.start.values, df.length.values)),\n",
" (0, 1))\n",
"ax3.set_ylabel(\"Multinomial model\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@hardingnj
Copy link
Author

In the last plot, I have the labels the wrong way around.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment