{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Sibship clustering" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tom Ellis, March 2017, updated June 2020" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created using FAPS version 2.6.6.\n" ] } ], "source": [ "import faps as fp\n", "import numpy as np\n", "import pandas as pd\n", "\n", "print(\"Created using FAPS version {}.\".format(fp.__version__))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "FAPS uses information in a `paternityArray` to generate plausible full-sibship configurations. This information is stored as a `sibshipCluster` object, and is the basis for almost all inference about biological processes in FAPS.\n", "\n", "This notebook will examine how to:\n", "\n", "1. Use a `paternityArray` to cluster offspring into plausible full-sibships.\n", "2. Compare the relative support for different partitions structures\n", "3. Make some basic inferences about the size and number of full sibships, and who is related to whom.\n", "\n", "Note that this tutorial only deals with the case where you have a single maternal family. If you have multiple families, you can apply what is here to each one, but you'll have to iterate over those families. See the specific [tutorial](https://fractional-analysis-of-paternity-and-sibships.readthedocs.io/en/latest/tutorials/07_dealing_with_multiple_half-sib_families.html) on that." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating a `sibshipCluster` object" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will begin by generating a population of 100 adults with 50 loci." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "np.random.seed(867)\n", "allele_freqs = np.random.uniform(0.3,0.5,50)\n", "adults = fp.make_parents(100, allele_freqs, family_name='a')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We take the first individal as the mother and mate her to four males, to create three full sibships of five offspring. We then generate a `paternityArray` based on the genotype data." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "progeny = fp.make_sibships(adults, 0, [1,2,3], 5, 'x')\n", "mothers = adults.subset(progeny.mothers)\n", "patlik = fp.paternity_array(progeny, mothers, adults, mu = 0.0015, missing_parents=0.01)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is straightforward to cluster offspring into full sibships. For now we'll stick with the default number of Monte Carlo draws." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sc = fp.sibship_clustering(patlik)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The default number of Monte Carlo draws is 1000, which seems to work in most cases. I have dropped it to 100 in cases where I wanted to call `sibship_clustering` many times, such as in an MCMC loop, when finding every possible candidate wasn't a priority. You could also use more draws if you really wanted to be sure you had completely sampled the space of compatible candidate fathers. Speeds are unlikely to increase linearly with number of draws:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "17.5 ms ± 1.43 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "35.5 ms ± 2.02 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "312 ms ± 2.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit fp.sibship_clustering(patlik, ndraws=100)\n", "%timeit fp.sibship_clustering(patlik, ndraws=1000)\n", "%timeit fp.sibship_clustering(patlik, ndraws=10000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We discussed this in figure 5 of the FAPS paper should you be interested in more on this." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sibling relationships" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sibship clustering calculates likelihoods that each pair of offspring are full siblings, then builds a dendrogram from this. We can visualise this dendrogram if we so wish, although the output is not pretty." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGhCAYAAADBddZJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnq0lEQVR4nO3de1zUdb7H8TfIXRwUTNAEpbxQGuUlFTMtlxWt1JLMWtvMw66ri5q6Wy3naG4clbba1XXDW+tindW8tMdKO9pRNq1TgEq1R9dWrfVC4oyVAgIyIPzOHz2Y0yimAzNfxV7Px+P3eMhvLt8PCPJy5jfz87MsyxIAAIAh/ld6AAAA8P1CfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMCogCs9wPnq6upUXFysVq1ayc/P70qPAwAALoNlWTpz5ow6dOggf//vfmzjqouP4uJixcbGXukxAABAIxQVFaljx47feZ2rLj5atWol6ZvhbTbbFZ4GAABcjrKyMsXGxrp+j3+Xqy4+6p9qsdlsxAcAAM3M5RwywQGnAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGDUVXdiuebOsiydram90mMA8JHQwBaXdeIsABdHfHiRZVl6cFmeCo+evtKjAPCRvp3aaMPkJAIEaAKedvGiszW1hAdwjdtz9DSPbgJNxCMfPrJndrLCglpc6TEAeEllda36ztt+pccArgnEh4+EBbVQWBBfXgAAzsfTLgAAwCjiAwAAGEV8AAAAo4gPAABgFPEBAACMIj4AAIBRxAcAADCKN6IArnGcb8g7KqvPNfhnNB7nyfn+Ij6AaxjnG/KNvvNyr/QI1wTOk/P9xdMuwDWM8w3hasZ5cr6/eOQD+J7gfEO4WnCeHBAfwPcE5xsCcLXgaRcAAGCUx/Fx/PhxPfroo4qKilJoaKhuueUW7dmzx3W5ZVl65pln1L59e4WGhio5OVmHDh3y6tAAAKD58ig+Tp8+rTvuuEOBgYHasmWL9u/fr9/+9rdq06aN6zrPP/+8Fi9erGXLlqmgoEAtW7ZUSkqKqqqqvD48AABofjx6Avg3v/mNYmNjlZOT49oXHx/v+rNlWVq0aJFmz56t0aNHS5JeffVVRUdH64033tDDDz/spbEBAEBz5dEjH2+99Zb69u2rsWPHql27durVq5defvll1+WHDx+W3W5XcnKya19ERIT69++vvLy8Bu/T6XSqrKzMbQMAANcuj+Ljn//8p5YuXaquXbvqnXfe0ZQpUzR9+nS98sorkiS73S5Jio6OdrtddHS067LzZWVlKSIiwrXFxsY25vMAAADNhEfxUVdXp969e2vBggXq1auXJk2apJ/+9KdatmxZowfIyMhQaWmpaysqKmr0fQEAgKufR/HRvn173XzzzW77brrpJh07dkySFBMTI0lyOBxu13E4HK7LzhccHCybzea2AQCAa5dH8XHHHXfowIEDbvsOHjyoTp06Sfrm4NOYmBjl5v7/eQ/KyspUUFCgpKQkL4wLAACaO49e7TJz5kwNHDhQCxYs0EMPPaRdu3ZpxYoVWrFihSTJz89PM2bM0Lx589S1a1fFx8drzpw56tChg+6//35fzA8AAJoZj+Lj9ttv18aNG5WRkaHMzEzFx8dr0aJFGj9+vOs6Tz31lCoqKjRp0iSVlJRo0KBB2rp1q0JCQrw+PAAAaH48PtHDfffdp/vuu++il/v5+SkzM1OZmZlNGgwAAFybOLcLAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABjlUXz8+te/lp+fn9uWkJDguryqqkrp6emKiopSeHi4UlNT5XA4vD40AABovjx+5KNHjx46ceKEa/uf//kf12UzZ87Upk2btGHDBu3cuVPFxcUaM2aMVwcGAADNW4DHNwgIUExMzAX7S0tLtXLlSq1Zs0ZDhw6VJOXk5Oimm25Sfn6+BgwY0PRpAQBAs+fxIx+HDh1Shw4ddMMNN2j8+PE6duyYJKmwsFA1NTVKTk52XTchIUFxcXHKy8vz3sQAAKBZ8+iRj/79+2vVqlXq3r27Tpw4oWeffVZ33nmn9u3bJ7vdrqCgILVu3drtNtHR0bLb7Re9T6fTKafT6fq4rKzMs88AAAA0Kx7Fx4gRI1x/TkxMVP/+/dWpUyetX79eoaGhjRogKytLzz77bKNuCwAAmp8mvdS2devW6tatmz777DPFxMSourpaJSUlbtdxOBwNHiNSLyMjQ6Wlpa6tqKioKSMBAICrXJPio7y8XJ9//rnat2+vPn36KDAwULm5ua7LDxw4oGPHjikpKemi9xEcHCybzea2AQCAa5dHT7v88pe/1MiRI9WpUycVFxdr7ty5atGihR555BFFREQoLS1Ns2bNUmRkpGw2m6ZNm6akpCRe6QIAAFw8io8vvvhCjzzyiL7++mtdd911GjRokPLz83XddddJkhYuXCh/f3+lpqbK6XQqJSVFS5Ys8cngAACgefIoPtauXfudl4eEhCg7O1vZ2dlNGgoAAFy7OLcLAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOaFB/PPfec/Pz8NGPGDNe+qqoqpaenKyoqSuHh4UpNTZXD4WjqnAAA4BrR6PjYvXu3li9frsTERLf9M2fO1KZNm7Rhwwbt3LlTxcXFGjNmTJMHBQAA14ZGxUd5ebnGjx+vl19+WW3atHHtLy0t1cqVK/W73/1OQ4cOVZ8+fZSTk6MPP/xQ+fn5XhsaAAA0X42Kj/T0dN17771KTk52219YWKiamhq3/QkJCYqLi1NeXl6D9+V0OlVWVua2AQCAa1eApzdYu3atPvroI+3evfuCy+x2u4KCgtS6dWu3/dHR0bLb7Q3eX1ZWlp599llPxwAAAM2UR498FBUV6YknntDq1asVEhLilQEyMjJUWlrq2oqKirxyvwAA4OrkUXwUFhbq5MmT6t27twICAhQQEKCdO3dq8eLFCggIUHR0tKqrq1VSUuJ2O4fDoZiYmAbvMzg4WDabzW0DAADXLo+edvnBD36gvXv3uu2bOHGiEhIS9PTTTys2NlaBgYHKzc1VamqqJOnAgQM6duyYkpKSvDc1AABotjyKj1atWqlnz55u+1q2bKmoqCjX/rS0NM2aNUuRkZGy2WyaNm2akpKSNGDAAO9NDQAAmi2PDzi9lIULF8rf31+pqalyOp1KSUnRkiVLvL0MAABoppocHzt27HD7OCQkRNnZ2crOzm7qXQMAgGsQ53YBAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKM8io+lS5cqMTFRNptNNptNSUlJ2rJli+vyqqoqpaenKyoqSuHh4UpNTZXD4fD60AAAoPnyKD46duyo5557ToWFhdqzZ4+GDh2q0aNH6+9//7skaebMmdq0aZM2bNignTt3qri4WGPGjPHJ4AAAoHkK8OTKI0eOdPt4/vz5Wrp0qfLz89WxY0etXLlSa9as0dChQyVJOTk5uummm5Sfn68BAwZ4b2oAANBsNfqYj9raWq1du1YVFRVKSkpSYWGhampqlJyc7LpOQkKC4uLilJeXd9H7cTqdKisrc9sAAMC1y+P42Lt3r8LDwxUcHKzJkydr48aNuvnmm2W32xUUFKTWrVu7XT86Olp2u/2i95eVlaWIiAjXFhsb6/EnAQAAmg+P46N79+765JNPVFBQoClTpmjChAnav39/owfIyMhQaWmpaysqKmr0fQEAgKufR8d8SFJQUJC6dOkiSerTp492796t3//+9xo3bpyqq6tVUlLi9uiHw+FQTEzMRe8vODhYwcHBnk8OAACapSa/z0ddXZ2cTqf69OmjwMBA5ebmui47cOCAjh07pqSkpKYuAwAArhEePfKRkZGhESNGKC4uTmfOnNGaNWu0Y8cOvfPOO4qIiFBaWppmzZqlyMhI2Ww2TZs2TUlJSbzSBQAAuHgUHydPntRjjz2mEydOKCIiQomJiXrnnXf0wx/+UJK0cOFC+fv7KzU1VU6nUykpKVqyZIlPBgcAAM2TR/GxcuXK77w8JCRE2dnZys7ObtJQAADg2sW5XQAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKI/iIysrS7fffrtatWqldu3a6f7779eBAwfcrlNVVaX09HRFRUUpPDxcqampcjgcXh0aAAA0Xx7Fx86dO5Wenq78/Hxt27ZNNTU1GjZsmCoqKlzXmTlzpjZt2qQNGzZo586dKi4u1pgxY7w+OAAAaJ4CPLny1q1b3T5etWqV2rVrp8LCQg0ePFilpaVauXKl1qxZo6FDh0qScnJydNNNNyk/P18DBgzw3uQAAKBZatIxH6WlpZKkyMhISVJhYaFqamqUnJzsuk5CQoLi4uKUl5fXlKUAAMA1wqNHPr6trq5OM2bM0B133KGePXtKkux2u4KCgtS6dWu360ZHR8tutzd4P06nU06n0/VxWVlZY0cCAADNQKMf+UhPT9e+ffu0du3aJg2QlZWliIgI1xYbG9uk+wMAAFe3RsXH1KlTtXnzZr377rvq2LGja39MTIyqq6tVUlLidn2Hw6GYmJgG7ysjI0OlpaWuraioqDEjAQCAZsKj+LAsS1OnTtXGjRv117/+VfHx8W6X9+nTR4GBgcrNzXXtO3DggI4dO6akpKQG7zM4OFg2m81tAwAA1y6PjvlIT0/XmjVr9Oabb6pVq1au4zgiIiIUGhqqiIgIpaWladasWYqMjJTNZtO0adOUlJTEK10AAIAkD+Nj6dKlkqS77rrLbX9OTo4ef/xxSdLChQvl7++v1NRUOZ1OpaSkaMmSJV4ZFgAANH8exYdlWZe8TkhIiLKzs5Wdnd3ooQAAwLWLc7sAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMMrj+Hjvvfc0cuRIdejQQX5+fnrjjTfcLrcsS88884zat2+v0NBQJScn69ChQ96aFwAANHMex0dFRYVuvfVWZWdnN3j5888/r8WLF2vZsmUqKChQy5YtlZKSoqqqqiYPCwAAmr8AT28wYsQIjRgxosHLLMvSokWLNHv2bI0ePVqS9Oqrryo6OlpvvPGGHn744aZNCwAAmj2vHvNx+PBh2e12JScnu/ZFRESof//+ysvLa/A2TqdTZWVlbhsAALh2eTU+7Ha7JCk6Otptf3R0tOuy82VlZSkiIsK1xcbGenMkAABwlbnir3bJyMhQaWmpaysqKrrSIwEAAB/yanzExMRIkhwOh9t+h8Phuux8wcHBstlsbhsAALh2eTU+4uPjFRMTo9zcXNe+srIyFRQUKCkpyZtLAQCAZsrjV7uUl5frs88+c318+PBhffLJJ4qMjFRcXJxmzJihefPmqWvXroqPj9ecOXPUoUMH3X///d6cGwAANFMex8eePXt09913uz6eNWuWJGnChAlatWqVnnrqKVVUVGjSpEkqKSnRoEGDtHXrVoWEhHhvagAA0Gx5HB933XWXLMu66OV+fn7KzMxUZmZmkwYDAADXpiv+ahcAAPD9QnwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIwiPgAAgFHEBwAAMIr4AAAARhEfAADAKOIDAAAYRXwAAACjiA8AAGAU8QEAAIzyWXxkZ2erc+fOCgkJUf/+/bVr1y5fLQUAAJoRn8THunXrNGvWLM2dO1cfffSRbr31VqWkpOjkyZO+WA4AADQjAb6409/97nf66U9/qokTJ0qSli1bprffflt/+tOf9Ktf/coXS16ousLMOm5r1n7rz5VStdP8DEFh7h9XV5pf83zemOFSa1zy9i2bPkMTVdYY+Lu4YM3ab/35rM7WVhufITQg1O3js+fOGl/zfN6Y4VJrXEpYYBO/p72krtL892Xdt/6trKs8q3NnzH9f+oe6//3VnfXt9+X56zWkqTNczhqXvI8wM9+XXo+P6upqFRYWKiMjw7XP399fycnJysvLu+D6TqdTTuf//5IuLS2VJJWVlTVtkKyOTbt9I1RaQapzLpUklc27Qef8zP9A4SIyvrjSE2jAmgHG17TqAlXnnC1JuvOVQfLzrzE+AxqW/6P8Kz2CJOlAn77G16xqEaS6EXMlSR8PSFLIFYhiNKx74Z5G37b+97ZlWZe+suVlx48ftyRZH374odv+J5980urXr98F1587d64liY2NjY2Nje0a2IqKii7ZCj552sUTGRkZmjVrluvjuro6nTp1SlFRUfLz87uCkwEAgMtlWZbOnDmjDh06XPK6Xo+Ptm3bqkWLFnI4HG77HQ6HYmJiLrh+cHCwgoOD3fa1bt3a22MBAAAfi4iIuKzref3VLkFBQerTp49yc3Nd++rq6pSbm6ukpCRvLwcAAJoZnzztMmvWLE2YMEF9+/ZVv379tGjRIlVUVLhe/QIAAL6/fBIf48aN05dffqlnnnlGdrtdt912m7Zu3aro6GhfLAcAAJoRP8u6nNfEAAAAeAfndgEAAEYRHwAAwCjiAwAAGEV8AAAAo5plfJSXl2vu3LkaPny4IiMj5efnp1WrVl1wvV27dunnP/+5+vTpo8DAQK++Y+rlziBJn376qYYPH67w8HBFRkbqxz/+sb788kujM9SrqanRzTffLD8/P7344otGZqirq9OqVas0atQoxcbGqmXLlurZs6fmzZunqqoqIzNI0ssvv6whQ4YoOjpawcHBio+P18SJE3XkyJEmz9CQ3bt3a+rUqerRo4datmypuLg4PfTQQzp48KBP1ruYxnyfeNPf//53jR07VjfccIPCwsLUtm1bDR48WJs2bTI2g/TNeaSefvppdejQQaGhoerfv7+2bdtmdIbzzZ8/X35+furZs6exNQ8dOqSHH35YHTt2VFhYmBISEpSZmalKgyeYKyws1PDhw2Wz2dSqVSsNGzZMn3zyibH1Jenxxx+Xn5/fRbfjx4/7fIYdO3ZcdP38fLPn/vnoo480atQoRUZGKiwsTD179tTixYt9tt4Vf3v1xvjqq6+UmZmpuLg43XrrrdqxY0eD1/uv//ov/fGPf1RiYqJuuOEGr/6jf7kzfPHFFxo8eLAiIiK0YMEClZeX68UXX9TevXu1a9cuBQUF+XyGb/vDH/6gY8eONXrNxsxQWVmpiRMnasCAAZo8ebLatWunvLw8zZ07V7m5ufrrX//apDC83K/Dxx9/rPj4eI0aNUpt2rTR4cOH9fLLL2vz5s3629/+dllvCeyJ3/zmN/rggw80duxYJSYmym6366WXXlLv3r2Vn59v7BdOY75PvOno0aM6c+aMJkyYoA4dOqiyslJ/+ctfNGrUKC1fvlyTJk0yMsfjjz+u119/XTNmzFDXrl21atUq3XPPPXr33Xc1aNAgIzN82xdffKEFCxaoZUtzZ1wuKipSv379FBERoalTpyoyMtL1s1hYWKg333zT5zN89NFHGjRokGJjYzV37lzV1dVpyZIlGjJkiHbt2qXu3bv7fAZJ+tnPfqbk5GS3fZZlafLkyercubOuv/56I3NI0vTp03X77be77evSpYux9f/7v/9bI0eOVK9evTRnzhyFh4fr888/1xdf+PCEnF45m5xhVVVV1okTJyzLsqzdu3dbkqycnJwLrme3263KykrLsiwrPT3d8uane7kzTJkyxQoNDbWOHj3q2rdt2zZLkrV8+XIjM9RzOBxWRESElZmZaUmyXnjhhSatf7kzOJ1O64MPPrjgts8++6wlydq2bZvPZ7iYPXv2WJKsrKysJs3QkA8++MByOp1u+w4ePGgFBwdb48eP9/p6F9OUr4+vnDt3zrr11lut7t27G1mvoKDggu/5s2fPWjfeeKOVlJRkZIbzjRs3zho6dKg1ZMgQq0ePHkbWnD9/viXJ2rdvn9v+xx57zJJknTp1yucz3HPPPVabNm2sr776yrWvuLjYCg8Pt8aMGePz9b/L+++/b0my5s+fb2S9d99915Jkbdiwwch6DSktLbWio6OtBx54wKqtrTW2brN82iU4OLjB88ScLzo6WqGhoVd0hr/85S+67777FBcX59qXnJysbt26af369UZmqPerX/1K3bt316OPPtqkdT2dISgoSAMHDrxg/wMPPCDpm6elfD3DxXTu3FmSVFJS0qQZGjJw4MALHtnq2rWrevTo0eTP2RNN+fr4SosWLRQbG+uTr3tDXn/9dbVo0cLtUZaQkBClpaUpLy9PRUVFRuao99577+n111/XokWLjK5bf8rz89/wsX379vL392/SI7GX6/3331dycrKioqLc1h8yZIg2b96s8vJyn89wMWvWrJGfn59+9KMfGV/7zJkzOnfunPF116xZI4fDofnz58vf318VFRWqq6vz+brNMj6ai+PHj+vkyZPq27fvBZf169dPH3/8sbFZdu3apVdeeUWLFi26as4WbLfbJX1zMkKTvv76a508eVJ79uxxveX/D37wAyNrW5Ylh8Nh/HO+GlRUVOirr77S559/roULF2rLli3Gvu4ff/yxunXrJpvN5ra/X79+kmT0eIPa2lpNmzZNP/nJT3TLLbcYW1eS7rrrLklSWlqaPvnkExUVFWndunVaunSppk+fbuQpIKfT2eB/CsPCwlRdXa19+/b5fIaG1NTUaP369Ro4cKDrPyWmTJw4UTabTSEhIbr77ru1Z88eY2tv375dNptNx48fV/fu3RUeHi6bzaYpU6Z45Zi8i2mWx3w0FydOnJD0TdWfr3379jp16pScTucFZ/X1NsuyNG3aNI0bN05JSUk+O8DSU88//7xsNptGjBhhdN3rr79eTqdTkhQVFaXFixfrhz/8oZG1V69erePHjyszM9PIeleTX/ziF1q+fLkkyd/fX2PGjNFLL71kZO0TJ05c9OdQkoqLi43MIUnLli3T0aNHtX37dmNr1hs+fLj+/d//XQsWLNBbb73l2v9v//ZvmjdvnpEZunfvrvz8fNXW1qpFixaSpOrqahUUFEiSkQM9G/LOO+/o66+/1vjx442tGRQUpNTUVN1zzz1q27at9u/frxdffFF33nmnPvzwQ/Xq1cvnMxw6dEjnzp3T6NGjlZaWpqysLO3YsUN/+MMfVFJSotdee80n6xIfPnT27FlJajAuQkJCXNfxdXysWrVKe/fu1euvv+7TdTyxYMECbd++XUuWLFHr1q2Nrr1lyxZVVVXp008/1Z///GdVVFQYWfcf//iH0tPTlZSUpAkTJhhZ82oyY8YMPfjggyouLtb69etVW1ur6upqI2tf7Ofs2z+HJnz99dd65plnNGfOHF133XVG1jxf586dNXjwYKWmpioqKkpvv/22FixYoJiYGE2dOtXn6//85z/XlClTlJaWpqeeekp1dXWaN2+e6z9rpv4uzrdmzRoFBgbqoYceMrbmwIED3Z6SHjVqlB588EElJiYqIyNDW7du9fkM5eXlqqys1OTJk12vbhkzZoyqq6u1fPlyZWZmqmvXrl5fl/jwofqHFuv/l/1t9Q9n+eqYlHplZWXKyMjQk08+qdjYWJ+udbnWrVun2bNnKy0tTVOmTDG+/t133y1JGjFihEaPHq2ePXsqPDzcp//w2u123XvvvYqIiHAdf/B9k5CQoISEBEnSY489pmHDhmnkyJEqKCjw+VOBoaGhV/TnsN7s2bMVGRmpadOmGVnvfGvXrtWkSZN08OBBdezYUdI3v2jq6ur09NNP65FHHnE7FsMXJk+erKKiIr3wwgt65ZVXJEl9+/bVU089pfnz5ys8PNyn6zekvLxcb775plJSUnz++V9Kly5dNHr0aP3nf/6n26NDvlL/vf/II4+47f/Rj36k5cuXKy8vzyfxwTEfPlT/kG590X/biRMnFBkZ6fNHPV588UVVV1dr3LhxOnLkiI4cOeJ6+dTp06d15MgRY//7lKRt27bpscce07333qtly5YZW/dibrzxRvXq1UurV6/22RqlpaUaMWKESkpKtHXrVq+/pLe5evDBB7V7924j73vSvn37i/4cSjLyd3Lo0CGtWLFC06dPV3FxsevnsaqqSjU1NTpy5IhOnTrl0xmWLFmiXr16ucKj3qhRo1RZWWnsOLT58+fL4XDo/fff1//+7/9q9+7droMcu3XrZmSGb3vjjTdUWVlp9CmX7xIbG6vq6mojj8rWf++ffxByu3btJH3ze8IXiA8fuv7663Xdddc1ePDQrl27dNttt/l8hmPHjun06dPq0aOH4uPjFR8frzvvvFPSN099xMfHa//+/T6fQ5IKCgr0wAMPqG/fvlq/fr0CAq6OB97Onj2r0tJSn9x3VVWVRo4cqYMHD2rz5s26+eabfbJOc1T/8Lqvvvbfdtttt+ngwYOuV3vUqz/OwMTP4vHjx1VXV6fp06e7fhbj4+NVUFCggwcPKj4+3ufHAjkcDtXW1l6wv6amRpKMvtqiTZs2GjRokOug2+3bt6tjx46uR8dMWr16tcLDwzVq1Cjjazfkn//8p0JCQow8CtSnTx9JFx5rU38clK+eHiQ+fCw1NVWbN292eylfbm6uDh48qLFjx/p8/enTp2vjxo1uW/1Bf48//rg2btyo+Ph4n8/x6aef6t5771Xnzp21efNmYw9z1zt37lyDBb9r1y7t3bu3wVckNVVtba3GjRunvLw8bdiwQUlJSV5fozk4efLkBftqamr06quvKjQ01EiQPfjgg6qtrdWKFStc+5xOp3JyctS/f38jT0n27Nnzgp/FjRs3qkePHoqLi9PGjRuVlpbm0xm6deumjz/++IJHm1577TX5+/srMTHRp+tfzLp167R7927NmDFD/v5mfy19+eWX2r59ux544AGFhYUZX/t8f/vb3/TWW29p2LBhRr4W9ce4rFy50m3/H//4RwUEBLheIeVtV8d/PRvhpZdeUklJiavONm3a5Ho6Ydq0aYqIiNDRo0f1H//xH5LkevSh/ojuTp066cc//rHPZ/jXf/1XbdiwQXfffbeeeOIJlZeX64UXXtAtt9ziepmnL2fo3bu3evfu7Xab+le79OjRQ/fff7/PZ/D391dKSopOnz6tJ598Um+//bbb7W+88cYm/2K+1AyWZSk2Nlbjxo1zvd353r17lZOTo4iICM2ZM6dJ6zfkF7/4hd566y2NHDlSp06d0p///Ge3y735fiuXcjnfq77ys5/9TGVlZRo8eLCuv/562e12rV69Wv/4xz/029/+1sj/7vr376+xY8cqIyNDJ0+eVJcuXfTKK6/oyJEjF/yj6ytt27Zt8Oet/r0+vPGzeClPPvmktmzZojvvvFNTp05VVFSUNm/erC1btugnP/mJkaef3nvvPWVmZmrYsGGKiopSfn6+cnJyNHz4cD3xxBM+X/9869at07lz567IUy7jxo1TaGioBg4cqHbt2mn//v1asWKFwsLC9NxzzxmZoVevXvqXf/kX/elPf9K5c+c0ZMgQ7dixQxs2bFBGRobvvieMvZ2Zl3Xq1MmS1OB2+PBhy7L+/93jGtqGDBliZAbLsqx9+/ZZw4YNs8LCwqzWrVtb48ePt+x2e5PX92SGbzt8+LDX3uH0cmaoX+9i24QJE3w+g9PptJ544gkrMTHRstlsVmBgoNWpUycrLS3tol+nphoyZMh3ft4mNeb7xFtee+01Kzk52YqOjrYCAgKsNm3aWMnJydabb77p03XPd/bsWeuXv/ylFRMTYwUHB1u33367tXXrVqMzNMTkO5xa1jfv9jpixAgrJibGCgwMtLp162bNnz/fqqmpMbL+Z599Zg0bNsxq27atFRwcbCUkJFhZWVkXvBuwKQMGDLDatWtnnTt3zvjav//9761+/fpZkZGRVkBAgNW+fXvr0UcftQ4dOmR0jurqauvXv/611alTJyswMNDq0qWLtXDhQp+u6WdZluW1kgEAALgEjvkAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYBTxAQAAjCI+AACAUcQHAAAwivgAAABGER8AAMAo4gMAABhFfAAAAKOIDwAAYNT/AUuoUWAewuBRAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.cluster.hierarchy import dendrogram\n", "import matplotlib.pyplot as plt\n", "\n", "dendrogram(sc.linkage_matrix)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Offspring individuals are labelled by their *index* in the array. Since full sibships are of size five we should expect to see clusters of {0,1,2,3,4}, {5,6,7,8,9} and {10,11,12,13,14}. This is indeed what we do see. What is difficult to see on the dendrogram are the branches connecting full siblings at the very bottom of the plot. If we bisect this dendrogram at different places on the y-axis we can infer different ways to partition the offspring into full siblings." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`sc` is an object of class `sibshipCluster` that contains various information about the array. Of primary interest are the set of partition structures inferred from the dendrogram. There are sixteen partitions - one for each individual in the array (i.e. one for each bifurcation in the dendrogram)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", " [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1],\n", " [ 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1],\n", " [ 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 1, 2, 1, 1, 1],\n", " [ 3, 3, 3, 3, 3, 4, 4, 4, 5, 4, 1, 2, 1, 1, 1],\n", " [ 3, 3, 4, 3, 3, 5, 5, 5, 6, 5, 1, 2, 1, 1, 1],\n", " [ 3, 3, 4, 3, 3, 5, 5, 5, 7, 6, 1, 2, 1, 1, 1],\n", " [ 3, 3, 4, 3, 3, 5, 5, 6, 8, 7, 1, 2, 1, 1, 1],\n", " [ 3, 4, 5, 3, 3, 6, 6, 7, 9, 8, 1, 2, 1, 1, 1],\n", " [ 3, 5, 6, 4, 3, 7, 7, 8, 10, 9, 1, 2, 1, 1, 1],\n", " [ 4, 6, 7, 5, 4, 8, 8, 9, 11, 10, 2, 3, 1, 1, 1],\n", " [ 4, 7, 8, 6, 5, 9, 10, 11, 13, 12, 2, 3, 1, 1, 1],\n", " [ 5, 8, 9, 7, 6, 10, 11, 12, 14, 13, 3, 4, 1, 1, 2],\n", " [ 6, 9, 10, 8, 7, 11, 12, 13, 15, 14, 4, 5, 1, 2, 3]],\n", " dtype=int32)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.partitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is key about partition structures is that each symbol represents a *unique but arbitrary* family identifier. For example in the third row we see the true partition structure, with individuals grouped into three groups of five individuals." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1], dtype=int32)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.partitions[2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Beyond denoting who is in a family with whom, the labels are arbitrary, with no meaningful order. This partition would be identical to `[0,0,0,0,0,1,1,1,1,1,2,2,2,2,2]` or `[10,10,10,10,10,7,7,7,7,7,22,22,22,22,22]` for example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each partition is associated with a log likelihood and equivalent log probability. We can see from both scores that the third partition is most consistent with the data. This is of course the true partition." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-4.71859462e+02 -3.07811673e+02 -1.54187774e-12 -inf\n", " -inf -inf -inf -inf\n", " -inf -inf -inf -inf\n", " -inf -inf]\n", "[1.18587580e-205 2.08491863e-134 1.00000000e+000 0.00000000e+000\n", " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n", " 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000\n", " 0.00000000e+000 0.00000000e+000]\n" ] } ], "source": [ "print(sc.lik_partitions) # log likelihood of each partition\n", "print(np.exp(sc.prob_partitions)) # probabilities of each partition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also see that the first and second partitions have non-zero, but small likelihoods. Parititons 5-8 have negative infinity log likelihood - they are incompatible with the data. These partitions split up true full siblings, and there is no way to reconcile this with the data. In real world situations such partitions might have non-zero likelihoods if they were an unrelated candidate male compatible with one or more offspring through chance alone." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases there can be rounding error when log probabilities are exponentiated and probabilities do not sum to one. This is classic machine error, and the reason it is good to work with log values wherever possible. We can check:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(sc.prob_partitions).sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can directly call the most likely partition. This is somewhat against the spirit of fractional analyses though..." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1], dtype=int32)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.mlpartition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For information about fine scale relationships, `sc.full_sib_matrix()` returns an $n*n$ matrix, where $n$ is the number of offspring. Each element describes the (log) probability that a pair of individuals are full siblings, averaged over partition structures and paternity configurations. If we plot this using a heatmap you can clearly see the five full sibships jump out as blocks of yellow (>90% probability of being full siblings) against a sea of purple (near zero probability of being full siblings)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAGiCAYAAABzmGX7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAraElEQVR4nO3df3RU9Z3/8dckIZMsJJEfS35gQrKuFeVHQCLZgN3CMWsOcqJ8e6xIKeRga9c2yI90KcQa4i+IYGWjkpJCS3HPkYLtFuoChYMpiB75mRhXThVk+ZWFTSKnJSOhJOzM/f6hzHIlwJ3cO5kZ5vk45/PH3Nw7n/ccNPPO+/259+MyDMMQAADAl2JCHQAAAAgvJAcAAMCE5AAAAJiQHAAAABOSAwAAYEJyAAAATEgOAACACckBAAAwITkAAAAmJAcAAMCE5AAAgDC1e/duFRcXKyMjQy6XS5s2bbrhNbt27dLdd98tt9utv//7v9fatWsDnpfkAACAMNXe3q7c3FzV1NRYOv/48eOaNGmSJkyYoMbGRs2dO1ff+973tH379oDmdbHxEgAA4c/lcmnjxo2aPHnyNc9ZsGCBtmzZokOHDvmPPfroozp37py2bdtmea44O4EGg8/n05kzZ5SUlCSXyxXqcAAAYcwwDH3++efKyMhQTEzwiuEXL15UZ2en7fcxDOOq7za32y232237vSVpz549KiwsNB0rKirS3LlzA3qfsEsOzpw5o8zMzFCHAQCIIE1NTbr11luD8t4XL15UzuA+am712n6vPn366Pz586ZjlZWVeuaZZ2y/tyQ1NzcrNTXVdCw1NVUej0d//etflZiYaOl9wi45SEpKkiTdqwcUp14hjgZAuNl45KNQh4Aw4jnv0+C7T/i/O4Khs7NTza1eHa8frOSk7lcnPJ/7lDP6pJqampScnOw/7lTVwElhlxxcLrfEqZfiXCQHAMySk2JDHQLCUE+0oZOTYmwlB/73SU42JQdOSktLU0tLi+lYS0uLkpOTLVcNpDBMDgAACEdewyevjSX8XsPnXDDXUFBQoK1bt5qO7dixQwUFBQG9D7cyAgBggU+G7RGo8+fPq7GxUY2NjZK+uFWxsbFRp06dkiSVl5drxowZ/vOfeOIJHTt2TD/+8Y/1ySef6Gc/+5nefPNNzZs3L6B5qRwAAGCBTz7Z+du/O1cfPHhQEyZM8L8uKyuTJJWUlGjt2rX6n//5H3+iIEk5OTnasmWL5s2bp1deeUW33nqrfvGLX6ioqCigeUkOAAAIU+PHj9f1HkfU1dMPx48frw8++MDWvCQHAABY4DUMeW08N9DOtT2N5AAAAAu6u27gyusjBQsSAQCACZUDAAAs8MmQN0oqByQHAABYQFsBAABELSoHAABYwN0KAADAxPflsHN9pKCtAAAATKgcAABggdfm3Qp2ru1pJAcAAFjgNWRzV0bnYgk2kgMAACxgzQEAAIhaVA4AALDAJ5e8ctm6PlIEXDnYvXu3iouLlZGRIZfLpU2bNl3z3CeeeEIul0vV1dU2QgQAIPR8hv0RKQJODtrb25Wbm6uamprrnrdx40bt3btXGRkZ3Q4OAAD0vIDbChMnTtTEiROve87p06f15JNPavv27Zo0adJ1z+3o6FBHR4f/tcfjCTQkAACCzmuzrWDn2p7m+IJEn8+n6dOna/78+Ro6dOgNz6+qqlJKSop/ZGZmOh0SAAC2XU4O7IxI4XhysHTpUsXFxWn27NmWzi8vL1dbW5t/NDU1OR0SAAAIgKN3K9TX1+uVV15RQ0ODXC5rGZLb7Zbb7XYyDAAAHOczXPIZNu5WsHFtT3O0cvDuu++qtbVVWVlZiouLU1xcnE6ePKkf/ehHys7OdnIqAAB6VDS1FRytHEyfPl2FhYWmY0VFRZo+fbpmzpzp5FQAACBIAk4Ozp8/r6NHj/pfHz9+XI2NjerXr5+ysrLUv39/0/m9evVSWlqa7rjjDvvRAgAQIl7FyGuj4O51MJZgCzg5OHjwoCZMmOB/XVZWJkkqKSnR2rVrHQsMAIBwYthcc2BE0JqDgJOD8ePHyzCsP+bpxIkTgU4BAEDY4TkHAAAgarHxEgAAFniNGHkNG2sOImhvBZIDAAAs8Mkln42Cu0+Rkx3QVgAAACZUDgAAsCCaFiSSHAAAYIH9NQe0FQAAQISicgAAgAVfLEi0sfESbQUAAG4uPpuPT+ZuBQAAELGoHAAAYEE0LUgkOQAAwAKfYqLmIUgkBwAAWOA1XPLa2FnRzrU9jTUHAADAhMoBAAAWeG3ereClrQAAwM3FZ8TIZ2NBoi+CFiTSVgAAACZUDgAAsIC2AgAAMPHJ3h0HPudCCTraCgAAwITKAQAAFth/CFLk/D1OcgAAgAX2H58cOclB5EQKAAB6BJUDAAAs8Mkln+wsSIycxyeTHAAAYEE0tRVIDgAAsMD+cw4iJzmInEgBAECPoHIAAIAFPsMln52HIEXQls0kBwAAWOCz2VaIpOccRE6kAACgR1A5AADAAvtbNkfO3+MkBwAAWOCVS14bzyqwc21Pi5w0BgAA9AgqBwAAWEBbAQAAmHhlrzXgdS6UoIucNAYAAPQIKgcAAFhAWwEAAJhE08ZLAUe6e/duFRcXKyMjQy6XS5s2bfL/7NKlS1qwYIGGDx+u3r17KyMjQzNmzNCZM2ecjBkAgB5nfLllc3eH0c31CjU1NcrOzlZCQoLy8/O1f//+655fXV2tO+64Q4mJicrMzNS8efN08eLFgOYMODlob29Xbm6uampqrvrZhQsX1NDQoIqKCjU0NOh3v/udDh8+rAcffDDQaQAAiHobNmxQWVmZKisr1dDQoNzcXBUVFam1tbXL89etW6eFCxeqsrJSH3/8sX75y19qw4YNeuqppwKaN+C2wsSJEzVx4sQuf5aSkqIdO3aYjq1YsUJjxozRqVOnlJWVFeh0AACEhVC0FZYvX67HH39cM2fOlCTV1tZqy5YtWrNmjRYuXHjV+e+//77GjRunb3/725Kk7OxsTZ06Vfv27Qto3qA3QNra2uRyuXTLLbd0+fOOjg55PB7TAAAg3FzeldHOkHTVd15HR0eX83V2dqq+vl6FhYX+YzExMSosLNSePXu6vGbs2LGqr6/3tx6OHTumrVu36oEHHgjoswY1Obh48aIWLFigqVOnKjk5uctzqqqqlJKS4h+ZmZnBDAkAgJDKzMw0fe9VVVV1ed7Zs2fl9XqVmppqOp6amqrm5uYur/n2t7+t5557Tvfee6969eql2267TePHjw+4rRC05ODSpUt65JFHZBiGVq5cec3zysvL1dbW5h9NTU3BCgkAgG7zfrlls50hSU1NTabvvfLycsdi3LVrl5YsWaKf/exn/rV/W7Zs0fPPPx/Q+wTlVsbLicHJkyf1xz/+8ZpVA0lyu91yu93BCAMAAMdc2Rro7vWSlJycfN3vxcsGDBig2NhYtbS0mI63tLQoLS2ty2sqKio0ffp0fe9735MkDR8+XO3t7fr+97+vn/zkJ4qJsVYTcLxycDkx+PTTT/X222+rf//+Tk8BAMBNLz4+XqNHj1ZdXZ3/mM/nU11dnQoKCrq85sKFC1clALGxsZIkwzAszx1w5eD8+fM6evSo//Xx48fV2Niofv36KT09XQ8//LAaGhq0efNmeb1ef1+kX79+io+PD3Q6AADCgk8x8tn4m7o715aVlamkpER5eXkaM2aMqqur1d7e7r97YcaMGRo0aJB/3UJxcbGWL1+uUaNGKT8/X0ePHlVFRYWKi4v9SYIVAScHBw8e1IQJE0yBS1JJSYmeeeYZvfXWW5KkkSNHmq7buXOnxo8fH+h0AACEBa/hktdGW6E7106ZMkWfffaZFi1apObmZo0cOVLbtm3zL1I8deqUqVLw9NNPy+Vy6emnn9bp06f1t3/7tyouLtbixYsDmtdlBFJn6AEej0cpKSkar4cU5+oV6nAAhJntZz4MdQgII57Pver7tWNqa2uz1Mfv1hxffi/94N1vyt2n+99LHecvaeXXfxfUWJ3C3goAAFjg1ILESEByAACABYbNXRmNCNp4ieQAAAALvHLJ283Nky5fHykiJ40BAAA9gsoBAAAW+Ax76wZ8YbX8//pIDgAAsMBnc82BnWt7WuRECgAAegSVAwAALPDJJZ+NRYV2ru1pJAcAAFgQiickhgptBQAAYELlAAAAC6JpQSLJAQAAFvhk8/HJEbTmIHLSGAAA0COoHAAAYIFh824FI4IqByQHAABYwK6MAADAJJoWJEZOpAAAoEdQOQAAwALaCgAAwCSaHp9MWwEAAJhQOQAAwALaCgAAwCSakgPaCgAAwITKAQAAFkRT5YDkAAAAC6IpOaCtAAAATKgcAABggSF7zyownAsl6EgOAACwIJraCiQHAABYEE3JAWsOAACACZUDAAAsiKbKAckBAAAWRFNyQFsBAACYUDkAAMACw3DJsPHXv51rexrJAQAAFvjksvWcAzvX9jTaCgAAwITKAQAAFkTTgkSSAwAALIimNQcBtxV2796t4uJiZWRkyOVyadOmTaafG4ahRYsWKT09XYmJiSosLNSnn37qVLwAACDIAk4O2tvblZubq5qami5/vmzZMr366quqra3Vvn371Lt3bxUVFenixYu2gwUAIFQutxXsjEgRcFth4sSJmjhxYpc/MwxD1dXVevrpp/XQQw9Jkv7t3/5Nqamp2rRpkx599FF70QIAECK0Fbrp+PHjam5uVmFhof9YSkqK8vPztWfPni6v6ejokMfjMQ0AAMKNYbNqEEnJgaMLEpubmyVJqamppuOpqan+n31VVVWVnn32WSfDwE1m+5kPQx0CwkhRRm6oQ0AY+V/jkqRjoQ7jphPy5xyUl5erra3NP5qamkIdEgAAVzEkGYaNEeoPEABHKwdpaWmSpJaWFqWnp/uPt7S0aOTIkV1e43a75Xa7nQwDAADH+eSSiyckBi4nJ0dpaWmqq6vzH/N4PNq3b58KCgqcnAoAAARJwJWD8+fP6+jRo/7Xx48fV2Njo/r166esrCzNnTtXL7zwgm6//Xbl5OSooqJCGRkZmjx5spNxAwDQo6LpboWAk4ODBw9qwoQJ/tdlZWWSpJKSEq1du1Y//vGP1d7eru9///s6d+6c7r33Xm3btk0JCQnORQ0AQA/zGS65eHxy18aPHy/DuPayCpfLpeeee07PPfecrcAAAEBosLcCAAAWXL7rwM71kYLkAAAAC6JpzUHIn3MAAADCC5UDAAAsiKbKAckBAAAWcLcCAAAwiaYFiaw5AAAAJlQOAACw4IvKgZ01Bw4GE2QkBwAAWBBNCxJpKwAAABMqBwAAWGB8OexcHylIDgAAsIC2AgAAiFpUDgAAsCKK+gpUDgAAsOLLtkJ3h7rZVqipqVF2drYSEhKUn5+v/fv3X/f8c+fOqbS0VOnp6XK73fra176mrVu3BjQnlQMAACwIxRMSN2zYoLKyMtXW1io/P1/V1dUqKirS4cOHNXDgwKvO7+zs1D/90z9p4MCB+u1vf6tBgwbp5MmTuuWWWwKal+QAAIAwtXz5cj3++OOaOXOmJKm2tlZbtmzRmjVrtHDhwqvOX7Nmjf785z/r/fffV69evSRJ2dnZAc9LWwEAAAvstBSuvNPB4/GYRkdHR5fzdXZ2qr6+XoWFhf5jMTExKiws1J49e7q85q233lJBQYFKS0uVmpqqYcOGacmSJfJ6vQF9VpIDAACsuLxuwM6QlJmZqZSUFP+oqqrqcrqzZ8/K6/UqNTXVdDw1NVXNzc1dXnPs2DH99re/ldfr1datW1VRUaGXX35ZL7zwQkAflbYCAAA9qKmpScnJyf7Xbrfbsff2+XwaOHCgVq1apdjYWI0ePVqnT5/WSy+9pMrKSsvvQ3IAAIAFTi1ITE5ONiUH1zJgwADFxsaqpaXFdLylpUVpaWldXpOenq5evXopNjbWf+zOO+9Uc3OzOjs7FR8fbylW2goAAFhhODACEB8fr9GjR6uurs5/zOfzqa6uTgUFBV1eM27cOB09elQ+n89/7MiRI0pPT7ecGEgkBwAAhK2ysjKtXr1ar7/+uj7++GP94Ac/UHt7u//uhRkzZqi8vNx//g9+8AP9+c9/1pw5c3TkyBFt2bJFS5YsUWlpaUDz0lYAAMCCUOytMGXKFH322WdatGiRmpubNXLkSG3bts2/SPHUqVOKifm/v/MzMzO1fft2zZs3TyNGjNCgQYM0Z84cLViwIKB5SQ4AALAqBI9AnjVrlmbNmtXlz3bt2nXVsYKCAu3du9fWnLQVAACACZUDAAAsiKYtm0kOAACwIop2ZSQ5AADAEteXw871kYE1BwAAwITKAQAAVtBWAAAAJlGUHNBWAAAAJlQOAACw4optl7t9fYQgOQAAwAKndmWMBLQVAACACZUDAACsiKIFiSQHAABYEUVrDhxvK3i9XlVUVCgnJ0eJiYm67bbb9Pzzz8uIpGYLAABRzPHKwdKlS7Vy5Uq9/vrrGjp0qA4ePKiZM2cqJSVFs2fPdno6AAB6hMv4Yti5PlI4nhy8//77euihhzRp0iRJUnZ2tn79619r//79Tk8FAEDPiaI1B463FcaOHau6ujodOXJEkvThhx/qvffe08SJE7s8v6OjQx6PxzQAAAg7l9cc2BkRwvHKwcKFC+XxeDRkyBDFxsbK6/Vq8eLFmjZtWpfnV1VV6dlnn3U6DAAA0E2OVw7efPNNvfHGG1q3bp0aGhr0+uuv66c//alef/31Ls8vLy9XW1ubfzQ1NTkdEgAA9hkOjAjheOVg/vz5WrhwoR599FFJ0vDhw3Xy5ElVVVWppKTkqvPdbrfcbrfTYQAA4CzWHHTfhQsXFBNjftvY2Fj5fD6npwIAAEHgeOWguLhYixcvVlZWloYOHaoPPvhAy5cv12OPPeb0VAAA9Jwoqhw4nhy89tprqqio0A9/+EO1trYqIyND//zP/6xFixY5PRUAAD0nip6Q6HhykJSUpOrqalVXVzv91gAAoAewtwIAABbwhEQAAGAWRWsOHL9bAQAARDaSAwAAYEJbAQAAC1yyuebAsUiCj+QAAAArouhWRtoKAADAhMoBAABWRNHdCiQHAABYEUXJAW0FAABgQuUAAAALeEIiAAAwo60AAACiFZUDAACsiKLKAckBAAAWRNOaA9oKAADAhMoBAABWRNHjk0kOAACwgjUHAADgSqw5AAAAUYvKAQAAVtBWAAAAJjbbCpGUHNBWAAAAJlQOAACwgrYCAAAwiaLkgLYCAAAwoXIAAIAFPOcAAABELZIDAABgQlsBAAAromhBIskBAAAWRNOaA5IDAACsiqAveDtYcwAAAEyoHAAAYAVrDgAAwJWiac0BbQUAAGBC5QAAACtoKwAAgCvRVgAAAFErKMnB6dOn9Z3vfEf9+/dXYmKihg8froMHDwZjKgAAeobhwOiGmpoaZWdnKyEhQfn5+dq/f7+l69avXy+Xy6XJkycHPKfjycFf/vIXjRs3Tr169dIf/vAH/elPf9LLL7+svn37Oj0VAAA9JwTJwYYNG1RWVqbKyko1NDQoNzdXRUVFam1tve51J06c0L/8y7/o61//euCTKghrDpYuXarMzEz96le/8h/Lycm55vkdHR3q6Ojwv/Z4PE6HBABA2Pjq95zb7Zbb7e7y3OXLl+vxxx/XzJkzJUm1tbXasmWL1qxZo4ULF3Z5jdfr1bRp0/Tss8/q3Xff1blz5wKO0fHKwVtvvaW8vDx961vf0sCBAzVq1CitXr36mudXVVUpJSXFPzIzM50OCQAA2y4vSLQzJCkzM9P0vVdVVdXlfJ2dnaqvr1dhYaH/WExMjAoLC7Vnz55rxvncc89p4MCB+u53v9vtz+p45eDYsWNauXKlysrK9NRTT+nAgQOaPXu24uPjVVJSctX55eXlKisr87/2eDwkCACA8OPQrYxNTU1KTk72H75W1eDs2bPyer1KTU01HU9NTdUnn3zS5TXvvfeefvnLX6qxsdFGoEFIDnw+n/Ly8rRkyRJJ0qhRo3To0CHV1tZ2mRxcr5wCAEDYcCg5SE5ONiUHTvn88881ffp0rV69WgMGDLD1Xo4nB+np6brrrrtMx+688079+7//u9NTAQBw0xowYIBiY2PV0tJiOt7S0qK0tLSrzv+v//ovnThxQsXFxf5jPp9PkhQXF6fDhw/rtttuszS342sOxo0bp8OHD5uOHTlyRIMHD3Z6KgAAeoxTaw6sio+P1+jRo1VXV+c/5vP5VFdXp4KCgqvOHzJkiD766CM1Njb6x4MPPqgJEyaosbExoJa945WDefPmaezYsVqyZIkeeeQR7d+/X6tWrdKqVaucngoAgJ4Tgscnl5WVqaSkRHl5eRozZoyqq6vV3t7uv3thxowZGjRokKqqqpSQkKBhw4aZrr/lllsk6arjN+J4cnDPPfdo48aNKi8v13PPPaecnBxVV1dr2rRpTk8FAMBNbcqUKfrss8+0aNEiNTc3a+TIkdq2bZt/keKpU6cUE+P88wxdhmGE1dOePR6PUlJSNF4PKc7VK9ThIAxsP/NhqENAGCnKyA11CAgj/2tc0i79Xm1tbUFZ5Cf93/fSnbOWKNad0O338XZc1McrngpqrE5h4yUAAKyIol0Z2XgJAACYUDkAAMCKKKockBwAAGCB68th5/pIQVsBAACYUDkAAMAK2goAAOBK3XnK4VevjxQkBwAAWBFFlQPWHAAAABMqBwAAWBVBf/3bQXIAAIAF0bTmgLYCAAAwoXIAAIAVUbQgkeQAAAALaCsAAICoReUAAAAraCsAAIAr0VYAAABRK2wrBxuPfKTkpNhQh4EwUJSRG+oQAIC2AgAA+AqSAwAAcCXWHAAAgKhF5QAAACtoKwAAgCu5DEMuo/vf8Hau7Wm0FQAAgAmVAwAArKCtAAAArsTdCgAAIGpROQAAwAraCgAA4Eq0FQAAQNSicgAAgBW0FQAAwJWiqa1AcgAAgBVRVDlgzQEAADChcgAAgEWR1Bqwg+QAAAArDOOLYef6CEFbAQAAmAQ9OXjxxRflcrk0d+7cYE8FAEDQXL5bwc6IFEFtKxw4cEA///nPNWLEiGBOAwBA8HG3gn3nz5/XtGnTtHr1avXt2zdY0wAAAIcFLTkoLS3VpEmTVFhYeN3zOjo65PF4TAMAgHDj8tkfkSIobYX169eroaFBBw4cuOG5VVVVevbZZ4MRBgAAzqGt0H1NTU2aM2eO3njjDSUkJNzw/PLycrW1tflHU1OT0yEBAIAAOF45qK+vV2trq+6++27/Ma/Xq927d2vFihXq6OhQbGys/2dut1tut9vpMAAAcBR7K9hw33336aOPPjIdmzlzpoYMGaIFCxaYEgMAACJGFD0EyfHkICkpScOGDTMd6927t/r373/VcQAAIkU0VQ54QiIAADDpkb0Vdu3a1RPTAAAQPFF0twIbLwEAYAFtBQAAELWoHAAAYAV3KwAAgCvRVgAAAFGLygEAAFZwtwIAALgSbQUAABC1qBwAAGCFz/hi2Lk+QpAcAABgBWsOAADAlVyyuebAsUiCjzUHAADAhMoBAABW8IREAABwJW5lBAAAYaGmpkbZ2dlKSEhQfn6+9u/ff81zV69era9//evq27ev+vbtq8LCwuuefy0kBwAAWGE4MAK0YcMGlZWVqbKyUg0NDcrNzVVRUZFaW1u7PH/Xrl2aOnWqdu7cqT179igzM1P333+/Tp8+HdC8JAcAAFjgMgzbQ5I8Ho9pdHR0XHPO5cuX6/HHH9fMmTN11113qba2Vn/zN3+jNWvWdHn+G2+8oR/+8IcaOXKkhgwZol/84hfy+Xyqq6sL6LOSHAAA0IMyMzOVkpLiH1VVVV2e19nZqfr6ehUWFvqPxcTEqLCwUHv27LE014ULF3Tp0iX169cvoBhZkAgAgBW+L4ed6yU1NTUpOTnZf9jtdnd5+tmzZ+X1epWammo6npqaqk8++cTSlAsWLFBGRoYpwbCC5AAAAAuubA1093pJSk5ONiUHwfLiiy9q/fr12rVrlxISEgK6luQAAIAwNGDAAMXGxqqlpcV0vKWlRWlpade99qc//alefPFFvf322xoxYkTAc7PmAAAAK3r4boX4+HiNHj3atJjw8uLCgoKCa163bNkyPf/889q2bZvy8vICm/RLVA4AALAiBE9ILCsrU0lJifLy8jRmzBhVV1ervb1dM2fOlCTNmDFDgwYN8i9qXLp0qRYtWqR169YpOztbzc3NkqQ+ffqoT58+luclOQAAwIJQPCFxypQp+uyzz7Ro0SI1Nzdr5MiR2rZtm3+R4qlTpxQT839NgJUrV6qzs1MPP/yw6X0qKyv1zDPPWJ6X5AAAgDA2a9YszZo1q8uf7dq1y/T6xIkTjsxJcgAAgBVsvAQAAK7k8n0x7FwfKbhbAQAAmFA5AADACtoKAADApJs7K5qujxC0FQAAgAmVAwAALHBqb4VIQHIAAIAVUbTmgLYCAAAwoXIAAIAVhiQ7zyqInMIByQEAAFaw5gAAAJgZsrnmwLFIgo41BwAAwMTx5KCqqkr33HOPkpKSNHDgQE2ePFmHDx92ehoAAHrW5bsV7IwI4Xhy8M4776i0tFR79+7Vjh07dOnSJd1///1qb293eioAAHqOz4ERIRxfc7Bt2zbT67Vr12rgwIGqr6/XP/7jPzo9HQAAcFjQFyS2tbVJkvr169flzzs6OtTR0eF/7fF4gh0SAAABi6a7FYK6INHn82nu3LkaN26chg0b1uU5VVVVSklJ8Y/MzMxghgQAQPew5sAZpaWlOnTokNavX3/Nc8rLy9XW1uYfTU1NwQwJAADcQNDaCrNmzdLmzZu1e/du3Xrrrdc8z+12y+12BysMAACcEUV7KzieHBiGoSeffFIbN27Url27lJOT4/QUAAD0PJKD7istLdW6dev0+9//XklJSWpubpYkpaSkKDEx0enpAACAwxxfc7By5Uq1tbVp/PjxSk9P948NGzY4PRUAAD2H5xx0nxFBZRMAAKyKplsZ2XgJAAAromjNARsvAQAAEyoHAABY4TMkl42//n2RUzkgOQAAwAraCgAAIFpROQAAwBK7+yNETuWA5AAAACtoKwAAgGhF5QAAACt8hmy1BrhbAQCAm4zh+2LYuT5C0FYAAAAmVA4AALAiihYkkhwAAGAFaw4AAIBJFFUOWHMAAABMqBwAAGCFIZuVA8ciCTqSAwAArKCtAAAAohWVAwAArPD5JNl4kJEvch6CRHIAAIAVtBUAAEC0onIAAIAVUVQ5IDkAAMCKKHpCIm0FAABgQuUAAAALDMMnw8a2y3au7WkkBwAAWGEY9loDrDkAAOAmY9hccxBByQFrDgAAgAmVAwAArPD5JJeNdQOsOQAA4CZDWwEAAEQrKgcAAFhg+HwybLQVuJURAICbDW0FAAAQragcAABghc+QXNFROSA5AADACsOQZOdWxshJDmgrAAAAEyoHAABYYPgMGTbaCkYEVQ5IDgAAsMLwyV5bIXJuZQxaW6GmpkbZ2dlKSEhQfn6+9u/fH6ypAAAIOsNn2B7dEej36W9+8xsNGTJECQkJGj58uLZu3RrwnEFJDjZs2KCysjJVVlaqoaFBubm5KioqUmtrazCmAwDgphTo9+n777+vqVOn6rvf/a4++OADTZ48WZMnT9ahQ4cCmtdlBKEJkp+fr3vuuUcrVqyQJPl8PmVmZurJJ5/UwoULTed2dHSoo6PD/7qtrU1ZWVk62ZCt5D6sl4T0/742PNQhAAhT/6tLek9bde7cOaWkpARlDo/Ho5SUFN2rBxSnXt1+n8uxNjU1KTk52X/c7XbL7XZ3eU0g36eSNGXKFLW3t2vz5s3+Y//wD/+gkSNHqra21nqwhsM6OjqM2NhYY+PGjabjM2bMMB588MGrzq+srLz8yCkGg8FgMLo1mpqanP468/vrX/9qpKWlORJnnz59rjpWWVnZ5byBfp8ahmFkZmYa//qv/2o6tmjRImPEiBEBfWbHFySePXtWXq9XqamppuOpqan65JNPrjq/vLxcZWVl/tfnzp3T4MGDderUqaBlgT3N4/EoMzPzqmwxkvGZIgOfKTLwmbrPMAx9/vnnysjICNocCQkJOn78uDo7O22/l2EYcrlcpmPXqhoE+n0qSc3NzV2e39zcHFCcIb9b4VrllJSUlJvmf5LLkpOT+UwRgM8UGfhMkaEnPlNP/CGZkJCghISEoM8TLhxv6g8YMECxsbFqaWkxHW9paVFaWprT0wEAcFPqzvdpWlqaI9+/jicH8fHxGj16tOrq6vzHfD6f6urqVFBQ4PR0AADclLrzfVpQUGA6X5J27NgR8PdvUNoKZWVlKikpUV5ensaMGaPq6mq1t7dr5syZN7zW7XarsrLymj2YSMRnigx8psjAZ4oMN+NnCoUbfZ/OmDFDgwYNUlVVlSRpzpw5+sY3vqGXX35ZkyZN0vr163Xw4EGtWrUqoHmDciujJK1YsUIvvfSSmpubNXLkSL366qvKz88PxlQAANy0rvd9On78eGVnZ2vt2rX+83/zm9/o6aef1okTJ3T77bdr2bJleuCBBwKaM2jJAQAAiEw8ZQgAAJiQHAAAABOSAwAAYEJyAAAATMIuObiZtnquqqrSPffco6SkJA0cOFCTJ0/W4cOHQx2WY1588UW5XC7NnTs31KHYdvr0aX3nO99R//79lZiYqOHDh+vgwYOhDqvbvF6vKioqlJOTo8TERN122216/vnnFUnrj3fv3q3i4mJlZGTI5XJp06ZNpp8bhqFFixYpPT1diYmJKiws1KeffhqaYC263me6dOmSFixYoOHDh6t3797KyMjQjBkzdObMmdAFfAM3+je60hNPPCGXy6Xq6uoeiw/dF1bJwc221fM777yj0tJS7d27Vzt27NClS5d0//33q729PdSh2XbgwAH9/Oc/14gRI0Idim1/+ctfNG7cOPXq1Ut/+MMf9Kc//Ukvv/yy+vbtG+rQum3p0qVauXKlVqxYoY8//lhLly7VsmXL9Nprr4U6NMva29uVm5urmpqaLn++bNkyvfrqq6qtrdW+ffvUu3dvFRUV6eLFiz0cqXXX+0wXLlxQQ0ODKioq1NDQoN/97nc6fPiwHnzwwRBEas2N/o0u27hxo/bu3RvU/Q/gsIC2aQqyMWPGGKWlpf7XXq/XyMjIMKqqqkIYlXNaW1sNScY777wT6lBs+fzzz43bb7/d2LFjh/GNb3zDmDNnTqhDsmXBggXGvffeG+owHDVp0iTjscceMx375je/aUybNi1EEdkjybQznc/nM9LS0oyXXnrJf+zcuXOG2+02fv3rX4cgwsB99TN1Zf/+/YYk4+TJkz0TlA3X+jz//d//bQwaNMg4dOiQMXjw4Kt2DER4CpvKQWdnp+rr61VYWOg/FhMTo8LCQu3ZsyeEkTmnra1NktSvX78QR2JPaWmpJk2aZPq3imRvvfWW8vLy9K1vfUsDBw7UqFGjtHr16lCHZcvYsWNVV1enI0eOSJI+/PBDvffee5o4cWKII3PG8ePH1dzcbPpvMCUlRfn5+TfN7wvpi98ZLpdLt9xyS6hD6Rafz6fp06dr/vz5Gjp0aKjDQQBCvivjZd3ZmjKS+Hw+zZ07V+PGjdOwYcNCHU63rV+/Xg0NDTpw4ECoQ3HMsWPHtHLlSpWVlempp57SgQMHNHv2bMXHx6ukpCTU4XXLwoUL5fF4NGTIEMXGxsrr9Wrx4sWaNm1aqENzxOXtZ53YmjZcXbx4UQsWLNDUqVMjdqfGpUuXKi4uTrNnzw51KAhQ2CQHN7vS0lIdOnRI7733XqhD6bampibNmTNHO3bsuKm2LvX5fMrLy9OSJUskSaNGjdKhQ4dUW1sbscnBm2++qTfeeEPr1q3T0KFD1djYqLlz5yojIyNiP1M0uXTpkh555BEZhqGVK1eGOpxuqa+v1yuvvKKGhga5XK5Qh4MAhU1b4Wbe6nnWrFnavHmzdu7cqVtvvTXU4XRbfX29WltbdffddysuLk5xcXF655139OqrryouLk5erzfUIXZLenq67rrrLtOxO++8U6dOnQpRRPbNnz9fCxcu1KOPPqrhw4dr+vTpmjdvnn9zlkh3+XfCzfj74nJicPLkSe3YsSNiqwbvvvuuWltblZWV5f99cfLkSf3oRz9SdnZ2qMPDDYRNcnAzbvVsGIZmzZqljRs36o9//KNycnJCHZIt9913nz766CM1Njb6R15enqZNm6bGxkbFxsaGOsRuGTdu3FW3mB45ckSDBw8OUUT2XbhwQTEx5v+9Y2Nj5fP5QhSRs3JycpSWlmb6feHxeLRv376I/X0h/V9i8Omnn+rtt99W//79Qx1St02fPl3/+Z//afp9kZGRofnz52v79u2hDg83EFZtBTtbPYej0tJSrVu3Tr///e+VlJTk74WmpKQoMTExxNEFLikp6ar1Er1791b//v0jeh3FvHnzNHbsWC1ZskSPPPKI9u/fr1WrVgW8xWk4KS4u1uLFi5WVlaWhQ4fqgw8+0PLly/XYY4+FOjTLzp8/r6NHj/pfHz9+XI2NjerXr5+ysrI0d+5cvfDCC7r99tuVk5OjiooKZWRkaPLkyaEL+gau95nS09P18MMPq6GhQZs3b5bX6/X/zujXr5/i4+NDFfY13ejf6KvJTa9evZSWlqY77rijp0NFoEJ9u8RXvfbaa0ZWVpYRHx9vjBkzxti7d2+oQ+o2SV2OX/3qV6EOzTE3w62MhmEY//Ef/2EMGzbMcLvdxpAhQ4xVq1aFOiRbPB6PMWfOHCMrK8tISEgw/u7v/s74yU9+YnR0dIQ6NMt27tzZ5f8/JSUlhmF8cTtjRUWFkZqaarjdbuO+++4zDh8+HNqgb+B6n+n48ePX/J2xc+fOUIfepRv9G30VtzJGDrZsBgAAJmGz5gAAAIQHkgMAAGBCcgAAAExIDgAAgAnJAQAAMCE5AAAAJiQHAADAhOQAAACYkBwAAAATkgMAAGBCcgAAAEz+P6rSpHY9pPfhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sibmat = sc.full_sib_matrix()\n", "plt.pcolor(np.exp(sibmat))\n", "plt.colorbar()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that real datasets seldom look this tidy!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inferring family structure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this section we will simulate a slightly more interesting family structure. This block of code creates a half-sib array of 15 offspring from five fathers, where each father contributes five, four, three, two and one offspring respectively. It then performs sibship clustering on the array. We use 1000 candidate males and 50 loci." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# Lists indexing sires and dams\n", "sires = [1]*5 + [2]*4 + [3]*3 + [4]*2 +[5]*1\n", "dam = [0] * len(sires)\n", "\n", "np.random.seed(542)\n", "allele_freqs = np.random.uniform(0.3,0.5,30)\n", "adults = fp.make_parents(1000, allele_freqs, family_name='a')\n", "progeny = fp.make_offspring(adults, dam_list=dam, sire_list=sires)\n", "mothers = adults.subset(progeny.mothers)\n", "\n", "patlik = fp.paternity_array(progeny, mothers, adults, mu= 0.0015, missing_parents=0.01)\n", "sc = fp.sibship_clustering(patlik)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Number of families" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We saw before that we could call a list of valid partitions for `sc` using `sc.partitions`. The output is not terribly enlightening on its own, however. We could instead ask how probable it is that there are *x* full sibships in the array, integrating over all partition structures. Here each number is the probability that there are 1, 2, ..., 15 families." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([9.34086609e-146, 3.27492546e-105, 1.31380725e-086, 3.81521349e-035,\n", " 7.66449138e-001, 2.09575936e-001, 2.36188753e-002, 3.35259124e-004,\n", " 2.07912635e-005, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.nfamilies()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could show the same information graphically. Its clear that almost all the probability denisty is around $x=5$ families." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa/klEQVR4nO3de7xdZX3n8c+XIEURvMDRai4kYNCiXLQh1aG1XmAmiJP4Ui6JN6iWqC9BVNQGbRkb2xmUaosttUZgcKwQI6hEiQZHAUdFSLgIJhRMwy2Rloio4CUQ+M4fa53lZuecvVeSs/Y6J3zfr9d+nb2e9ey1f/vk8tvPs9b6PbJNREQEwC5tBxAREeNHkkJERFSSFCIiopKkEBERlSSFiIioJClERESl0aQgaY6kWyWtk7RohP3TJF0h6QZJN0l6VZPxREREb2rqPgVJk4DbgCOBDcAqYIHttR19lgA32P6UpAOBFbanNxJQRET0tWuDx54NrLO9HkDSUmAesLajj4G9yudPAX7S76D77LOPp0+fPraRRkTs5K677rqf2h7q16/JpDAZuLtjewPwR119PgxcLukUYA/giH4HnT59OqtXrx6rGCMiHhck3VmnX9snmhcAF9ieArwK+JykrWKStFDSakmrN23aNPAgIyIeL5pMChuBqR3bU8q2Tm8FlgHYvhrYHdin+0C2l9ieZXvW0FDf0U9ERGynJpPCKmCmpBmSdgPmA8u7+twFvBJA0h9QJIUMBSIiWtJYUrC9BTgZWAncAiyzvUbSYklzy26nASdJ+iFwEXCiU7Y1IqI1TZ5oxvYKYEVX2xkdz9cChzcZQ0RE1Nf2ieaIiBhHkhQiIqKSpBAREZUkhYiIqDR6ojl2XtMXXbZDr7/jzKPHKJKIGEsZKURERCVJISIiKkkKERFRSVKIiIhKkkJERFSSFCIiopKkEBERlSSFiIioJClEREQlSSEiIipJChERUUlSiIiISpJCRERUGk0KkuZIulXSOkmLRtj/95JuLB+3Sfp5k/FERERvjZXOljQJOAc4EtgArJK0vFyXGQDb7+nofwrwwqbiiYiI/pocKcwG1tleb/shYCkwr0f/BcBFDcYTERF9NJkUJgN3d2xvKNu2ImlfYAbw7QbjiYiIPsbLieb5wMW2Hxlpp6SFklZLWr1p06YBhxYR8fjRZFLYCEzt2J5Sto1kPj2mjmwvsT3L9qyhoaExDDEiIjo1mRRWATMlzZC0G8V//Mu7O0l6HvA04OoGY4mIiBoaSwq2twAnAyuBW4BlttdIWixpbkfX+cBS224qloiIqKexS1IBbK8AVnS1ndG1/eEmY4iIiPrGy4nmiIgYB5IUIiKikqQQERGVJIWIiKgkKURERCVJISIiKkkKERFRSVKIiIhKkkJERFSSFCIiopKkEBERlSSFiIioJClEREQlSSEiIipJChERUUlSiIiISpJCRERUkhQiIqLSaFKQNEfSrZLWSVo0Sp/jJK2VtEbShU3GExERvTW2RrOkScA5wJHABmCVpOW213b0mQmcDhxu+35Jz2gqnoiI6K/JkcJsYJ3t9bYfApYC87r6nAScY/t+ANv3NhhPRET00WRSmAzc3bG9oWzrdABwgKTvSfqBpDkNxhMREX00Nn20De8/E3gZMAX4jqSDbP+8s5OkhcBCgGnTpg04xIiIx48mRwobgakd21PKtk4bgOW2H7Z9O3AbRZJ4DNtLbM+yPWtoaKixgCMiHu+aTAqrgJmSZkjaDZgPLO/q8xWKUQKS9qGYTlrfYEwREdFDY0nB9hbgZGAlcAuwzPYaSYslzS27rQTuk7QWuAJ4v+37moopIiJ6a/Scgu0VwIqutjM6nht4b/mIiIiW9R0pSLpO0jslPW0QAUVERHvqTB8dDzyb4uazpZL+myQ1HFdERLSgb1Kwvc72hyhOAl8InA/cKemvJT296QAjImJwap1olnQw8HHgLOAS4Fjgl8C3mwstIiIGre+JZknXAT8HzgMW2d5c7rpG0uENxhYREQNW5+qjY20/5t4BSTNs3277tQ3FFRERLagzfXRxzbaIiJjgRh0pSHoe8HzgKZI6RwR7Abs3HVhERAxer+mj5wKvBp4K/PeO9gcoSl5HRMROZtSkYPtS4FJJL7F99QBjioiIlvSaPvqA7Y8Br5e0oHu/7Xc1GllERAxcr+mjW8qfqwcRSEREtK/X9NFXy5+fHVw4ERHRpl7TR18FPNp+23NH2xcRERNTr+mjvxtYFBERMS70mj66apCBRERE+3pNHy2zfZykm3nsNJIo1sc5uPHoIiJioHpNH51a/nz1IAKJiIj2jVr7yPY95c87gc3AIcDBwOayrS9JcyTdKmmdpEUj7D9R0iZJN5aPP9++jxEREWOhznKcfw5cC7wWOAb4gaS31HjdJOAc4CjgQGCBpANH6PoF24eWj3O3KfqIiBhTdUpnvx94oe37ACTtDXyfYgW2XmYD64bLbktaCswD1m5/uBER0aQ6pbPvoyiCN+yBsq2fycDdHdsbyrZur5N0k6SLJU2tcdyIiGhIr6uP3ls+XUexytqlFFchzQNuGqP3/ypwke3Nkt4GfBZ4xQixLAQWAkybNm2M3joiIrr1GinsWT7+HfgKv7ss9VLg9hrH3gh0fvOfUrZVbN/XsbznucAfjnQg20tsz7I9a2hoqMZbR0TE9uh189pf7+CxVwEzJc2gSAbzgdd3dpD0rOGrnIC5/K4IX0REtKDviWZJQ8AHKFZhq1Zcs73VNE8n21sknQysBCYB59teI2kxsNr2cuBdkuYCW4CfASdu7weJiIgdV+fqo88DX6C4ie3twAnApjoHt70CWNHVdkbH89OB0+sGGxERzapz9dHets8DHrZ9le23MMLJ4IiImPjqjBQeLn/eI+lo4CfA05sLKSIi2lInKfyNpKcApwH/COwFvKfRqCIiohV9k4Ltr5VPfwG8vNlwIiKiTXVqH+0n6auSfirpXkmXStpvEMFFRMRg1TnRfCGwDPh94NnAF4GLmgwqIiLaUScpPMn252xvKR//Ssf9ChERsfPoVfto+Aqjr5drISylKHVxPF33HkRExM6h14nm6yiSgMrtt3XsM7npLCJip9Or9tGMQQYSERHtq1P76AnAO4CXlk1XAp+2/fCoL4qIiAmpzs1rnwKeAPxzuf2msi3rKUdE7GTqJIXDbB/Ssf1tST9sKqCIiGhPnUtSH5G0//BGeePaI82FFBERbakzUngfcIWk9RRXIu0L/FmjUUVERCt6JgVJk4BDgJnAc8vmWzuW0IyIiJ1Iz+kj248AC2xvtn1T+UhCiIjYSdWZPvqepH+iWH3tV8ONtq9vLKqIiGhFnaRwaPlzcUebqbH6mqQ5wNkUazSfa/vMUfq9DriY4kqn1TViioiIBtRZT2G71lAoz0ecAxwJbABWSVpue21Xvz2BU4Frtud9IiJi7NRZT2FvSZ+UdL2k6ySdLWnvGseeDayzvd72QxQF9eaN0O8jwEeB325T5BERMebq3KewFNgEvA44pnz+hRqvmwzc3bG9oWyrSHoRMNX2ZbWijYiIRtU5p/As2x/p2P4bScfv6BtL2gX4BHBijb4LgYUA06ZN29G3joiIUdQZKVwuab6kXcrHccDKGq/bCEzt2J5Stg3bE3gBcKWkO4AXA8slzeo+kO0ltmfZnjU0NFTjrSMiYnvUSQonUSzJubl8LAXeJukBSb/s8bpVwExJMyTtBswHlg/vtP0L2/vYnm57OvADYG6uPoqIaE+dq4/23J4D294i6WSKUcUk4HzbayQtBlbbXt77CBERMWh1zilsN9sr6Fq60/YZo/R9WZOxREREf3WmjyIi4nEiSSEiIio9k4KkSZL+bVDBREREu+pUSb1VUm4OiIh4HKhzovlpwBpJ1/LYKqlzG4sqIiJaUScp/FXjUURExLhQ5z6FqyQ9EzisbLrW9r3NhhUREW2oUyX1OOBa4FjgOOAaScc0HVhERAxenemjD1EsfnMvgKQh4P9SLIoTERE7kTr3KezSNV10X83XRUTEBFNnpPANSSuBi8rt4+kqXRERETuHOiea31+uoXx42bTE9pebDSsiItpQqyCe7UuASxqOJSIiWjZqUpD0Xdt/LOkBwJ27ANveq/HoIiJioEZNCrb/uPy5XespRETExJOCeBERUUlBvIiIqKQgXkREVBotiCdpDnA2xRrN59o+s2v/24F3Ao8ADwILba/d3veLiIgd0/fOZNtXAXcATyifrwKu7/c6SZOAc4CjgAOBBZIO7Op2oe2DbB8KfAz4xDZFHxERY6pOQbyTKOocfbpsmgx8pcaxZwPrbK+3/RCwFJjX2cH2Lzs29+Cxl75GRMSA1Zk+eifFf/DXANj+saRn1HjdZODuju0NwB91d5L0TuC9wG7AK2ocNyIiGlKnsN3m8ps+AJJ2ZQy/0ds+x/b+wF8AfzlSH0kLJa2WtHrTpk1j9dYREdGlTlK4StIHgSdKOhL4IvDVGq/bCEzt2J5Sto1mKfCakXbYXmJ7lu1ZQ0NDNd46IiK2R52ksAjYBNwMvA1YYftDNV63CpgpaYak3YD5wPLODpJmdmweDfy4VtQREdGIOucUTrF9NvCZ4QZJp5Zto7K9RdLJwEqKS1LPt71G0mJgte3lwMmSjgAeBu4HTtjeDxIRETuuTlI4geJeg04njtC2Fdsr6Fp7wfYZHc9PrfH+ERExIL2qpC4AXg/MkNQ57bMX8LOmA4uIiMHrNVL4PnAPsA/w8Y72B4CbmgwqIiLa0at09p3AneWc/29sPyrpAOB5FCedIyJiJ1Pn6qPvALtLmgxcDrwJuKDJoCIioh11koJs/xp4LfDPto8Fnt9sWBER0YZaSUHSS4A3AJeVbZOaCykiItpSJym8Gzgd+HJ5n8F+wBWNRhUREa3oe59CWS77KklPlvRk2+uBdzUfWkREDFqd0tkHSboBWAOslXSdpJxTiIjYCdWZPvo08F7b+9qeBpxGR8mLiIjYedRJCnvYrs4h2L6SYkGciIjYydSpfbRe0l8Bnyu33wisby6kiIhoS52RwluAIeBLwCUUZS/e0mRQERHRjl4F8XYH3g48h6KsxWm2Hx5UYBERMXi9RgqfBWZRJISjgLMGElFERLSm1zmFA20fBCDpPODawYQUERFt6TVSqKaKbG8ZQCwREdGyXiOFQyT9snwu4InltgDb3qvx6CIiYqBGHSnYnmR7r/Kxp+1dO57XSgiS5ki6VdI6SYtG2P9eSWsl3STpW5L23ZEPExERO6bOfQrbRdIk4BzgSGADsErScttrO7rdAMyy/WtJ7wA+BhzfVEwxfk1fdFn/Tn3ccebRYxBJxONbnfsUttdsYJ3t9bYfApYC8zo72L6iXKsB4AfAlAbjiYiIPppMCpOBuzu2N5Rto3kr8PUG44mIiD4amz7aFpLeSHFPxJ+Osn8hsBBg2rRpA4wsIuLxpcmRwkZgasf2lLLtMSQdAXwImGt780gHsr3E9izbs4aGhhoJNiIimk0Kq4CZkmZI2g2YDyzv7CDphRSluefavrfBWCIioobGkkJ5w9vJwErgFmBZuZznYklzy25nAU8GvijpRknLRzlcREQMQKPnFGyvAFZ0tZ3R8fyIJt8/IiK2TZPTRxERMcEkKURERCVJISIiKkkKERFRSVKIiIhKkkJERFSSFCIiopKkEBERlSSFiIioJClEREQlSSEiIipJChERUUlSiIiISpJCRERUkhQiIqKSpBAREZUkhYiIqCQpREREpdGkIGmOpFslrZO0aIT9L5V0vaQtko5pMpaIiOivsaQgaRJwDnAUcCCwQNKBXd3uAk4ELmwqjoiIqG/XBo89G1hnez2ApKXAPGDtcAfbd5T7Hm0wjoiIqKnJ6aPJwN0d2xvKtoiIGKcmxIlmSQslrZa0etOmTW2HExGx02oyKWwEpnZsTynbtpntJbZn2Z41NDQ0JsFFRMTWmkwKq4CZkmZI2g2YDyxv8P0iImIHNZYUbG8BTgZWArcAy2yvkbRY0lwASYdJ2gAcC3xa0pqm4omIiP6avPoI2yuAFV1tZ3Q8X0UxrRQREePAhDjRHBERg5GkEBERlSSFiIioJClEREQlSSEiIipJChERUUlSiIiISpJCRERUkhQiIqKSpBAREZUkhYiIqCQpREREJUkhIiIqSQoREVFJUoiIiEqj6ylEtGn6ost26PV3nHn0GEUSMXFkpBAREZUkhYiIqDSaFCTNkXSrpHWSFo2w//ckfaHcf42k6U3GExERvTWWFCRNAs4BjgIOBBZIOrCr21uB+20/B/h74KNNxRMREf01OVKYDayzvd72Q8BSYF5Xn3nAZ8vnFwOvlKQGY4qIiB6aTAqTgbs7tjeUbSP2sb0F+AWwd4MxRUREDxPiklRJC4GF5eaDkm4tn+8D/LSdqGpLjCPQtk8U9o1xO47ZUxMxjhMTIc7EODY6Y9y3zguaTAobgakd21PKtpH6bJC0K/AU4L7uA9leAizpbpe02vasMYu4AYlxbCTGsTMR4kyMY2N7Ymxy+mgVMFPSDEm7AfOB5V19lgMnlM+PAb5t2w3GFBERPTQ2UrC9RdLJwEpgEnC+7TWSFgOrbS8HzgM+J2kd8DOKxBERES1p9JyC7RXAiq62Mzqe/xY4dgfeYqsppXEoMY6NxDh2JkKciXFsbHOMymxNREQMS5mLiIioTMik0K98xnggaaqkKyStlbRG0qltxzQaSZMk3SDpa23HMhJJT5V0saR/k3SLpJe0HVM3Se8p/5x/JOkiSbuPg5jOl3SvpB91tD1d0jcl/bj8+bRxGONZ5Z/1TZK+LOmpLYY4HNNWcXbsO02SJe3TRmwdcYwYo6RTyt/nGkkf63ecCZcUapbPGA+2AKfZPhB4MfDOcRonwKnALW0H0cPZwDdsPw84hHEWq6TJwLuAWbZfQHFhxXi4aOICYE5X2yLgW7ZnAt8qt9t0AVvH+E3gBbYPBm4DTh90UCO4gK3jRNJU4L8Cdw06oBFcQFeMkl5OUTniENvPB/6u30EmXFKgXvmM1tm+x/b15fMHKP4j676ju3WSpgBHA+e2HctIJD0FeCnFlWrYfsj2z1sNamS7Ak8s77d5EvCTluPB9ncorurr1Fla5rPAawYZU7eRYrR9eVnhAOAHFPc4tWqU3yUUNds+ALR+cnaUGN8BnGl7c9nn3n7HmYhJoU75jHGlrP76QuCalkMZyT9Q/KV+tOU4RjMD2AT873KK61xJe7QdVCfbGym+gd0F3AP8wvbl7UY1qmfavqd8/h/AM9sMpoa3AF9vO4iRSJoHbLT9w7Zj6eEA4E/KKtRXSTqs3wsmYlKYUCQ9GbgEeLftX7YdTydJrwbutX1d27H0sCvwIuBTtl8I/Ir2pzweo5yXn0eRwJ4N7CHpje1G1V95o2jr33BHI+lDFNOwn287lm6SngR8EDijX9+W7Qo8nWIK+/3Asn5FRydiUqhTPmNckPQEioTwedtfajueERwOzJV0B8U03Csk/Wu7IW1lA7DB9vAo62KKJDGeHAHcbnuT7YeBLwH/peWYRvOfkp4FUP7sO53QBkknAq8G3jBOqxzsT/El4Iflv58pwPWSfr/VqLa2AfiSC9dSzAj0PCE+EZNCnfIZrSuz8XnALbY/0XY8I7F9uu0ptqdT/B6/bXtcfcO1/R/A3ZKeWza9EljbYkgjuQt4saQnlX/ur2ScnQzv0Fla5gTg0hZjGZGkORRTmnNt/7rteEZi+2bbz7A9vfz3swF4Ufn3dTz5CvByAEkHALvRp4jfhEsK5Qmo4fIZtwDLbK9pN6oRHQ68ieLb943l41VtBzVBnQJ8XtJNwKHA/2w3nMcqRzEXA9cDN1P8u2r9bldJFwFXA8+VtEHSW4EzgSMl/ZhihHPmOIzxn4A9gW+W/27+pc0YYdQ4x5VRYjwf2K+8THUpcEK/kVfuaI6IiMqEGylERERzkhQiIqKSpBAREZUkhYiIqCQpREREJUkhBqasJPnxju33SfrwGB37AknHjMWx+rzPsWWl1itG2HdWWYnyrB6vnz5cxVLSy7alMq2k72/n666UtNU6vZJmSfpk3ePE40OjK69FdNkMvFbS/7Ld8waaQZK0a0cBtn7eCpxk+7sj7FsIPN32I2MX3e/YHtO7pG2vBlaP5TFj4stIIQZpC8VNXe/p3tH9TV/Sg+XPl5WFvC6VtF7SmZLeIOlaSTdL2r/jMEdIWi3ptrKu0/BaEWdJWlXW539bx3H/n6TljHCHtKQF5fF/JOmjZdsZwB8D53WPBsrjPBm4TtLxo32eOiQ9v/x8N5YxzxzhGHtJukzFuiL/ImmX8rNeUMZ8s6TO3/Ox5TFvk/QnHb+Dr5XPPyzpc5KuVrHWwkll+7MkfaeM5UfDr42dV0YKMWjnADepxmIfHQ4B/oCiLPB64Fzbs1UsXHQK8O6y33SK0ur7A1dIeg7wZoqqpYdJ+j3ge5KGK5i+iKJu/+2dbybp2cBHgT8E7gcul/Qa24slvQJ4X/ktu2J7rqQHbR9aHuOobfh83d4OnG378ypKuUwaoc9sivVE7gS+AbwWuB2YXK7pgB67OM2u5e/sVcD/oLibudvBFIXT9gBukHQZsABYaftvVaxl8qQd+FwxAWSkEANVVor9PxSL0tS1qlyfYjPw78Dwf+o3UySCYctsP2r7xxTJ43kUC6C8WdKNFKXL9wZmlv2v7U4IpcOAK8sCd8NVOl+6DfHuqKuBD0r6C2Bf278Zoc+15ZoijwAXUYxg1lOUNPhHFfWDOqvyDhdkvI7H/s46XWr7N+XU3hUUiWcV8GfluZ+DyrVBYieWpBBt+AeKufnOdRG2UP59lLQLReGuYZs7nj/asf0ojx3tdtdsMSDgFNuHlo8ZHWsd/GpHPkQfvT5PT7YvBOYCvwFWlKOTrbpt/TLfTzGqupJitNG5cNLw7+wRRp8hGOmY36FIiBuBCyS9ue7niIkpSSEGzvbPgGUUiWHYHRTTNVD8h/iE7Tj0seXc+v7AfsCtFIUT36GijDmSDlD/RXquBf5U0j7llMkC4KptjOUOtvPzSNoPWG/7kxRVTA8eodtsFZWCdwGOB76rYo3gXWxfAvwl215ifJ6k3SXtDbwMWCVpX+A/bX+GIsmMt7LlMcZyTiHa8nGKarfDPgNcKumHFHPk2/Mt/i6K/9D3At5u+7eSzqWYLrlekihWcXtNr4PYvkfSIoopFAGX2d7WEtM78nmOA94k6WGK1dFGqgq7iqKa6HPKOL8MHESxQt3wl71tXdv4pvJY+wAfsf0TSScA7y9jeZDiHE3sxFIlNSIozxk8aLvvwu6xc8v0UUREVDJSiIiISkYKERFRSVKIiIhKkkJERFSSFCIiopKkEBERlSSFiIio/H9G8kikYTG/wAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "#ax.bar(np.arange(0.5, len(sc.nfamilies())+0.5), sc.nfamilies())\n", "ax.bar(np.arange(1,16), sc.nfamilies())\n", "ax.set_xlabel('Number of full sibships')\n", "ax.set_ylabel('Posterior probability')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Family size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also get the distribution of family sizes within the array, averaged over all partitions. This returns a vector of the same length as the number of offspring in the array. `family_size` returns the posterior probability of observing one or more families of size 1, 2, ... , 15. It will be clear that we are unable to distinguish a single sibship with high probability from multiple families of the same size, each with low probability; this is the price we pay for being able to integrate out uncertainty in partition structure." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2.33452339e-001, 1.95053525e-001, 1.88263368e-001, 2.29940941e-001,\n", " 1.53289828e-001, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 4.37935750e-087, 0.00000000e+000, 1.63746273e-105, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 9.34086609e-146])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.family_size()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting this shows that we are roughly equally likely to observe a family of sizes one, two, three, four and five. " ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOhUlEQVR4nO3cf6zdd13H8efL1oGAwsauqG3lVqhoRWTkUlDi/nAwiiMtf0AsiilxyWLCFAVjiiRbUhJTxCgkTmFhBYLAmAVjI8XRwNQ/dNi7wTa6OVfK3G4d7kIn/gAZZW//ON+Zs7PT3W93z+2597PnI7np+f46fd/e3uf93u/5kapCktSu75n2AJKklWXoJalxhl6SGmfoJalxhl6SGrd+2gOMOv/882t2dnbaY0jSmnLTTTd9rapmxm1bdaGfnZ1lfn5+2mNI0pqS5F9Pt81LN5LUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUuFX3ytjlmt3zqWUdf/e+SyY0iSStDp7RS1LjDL0kNa65SzdPRMu9XAVespJa5hm9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS43qFPsn2JHcmOZZkz5jtb0lye5Jbk3w2ybOHtu1Oclf3sXuSw0uSlrZk6JOsA64CXgVsBV6fZOvIbl8A5qrqBcAB4A+6Y88DrgReAmwDrkxy7uTGlyQtpc8Z/TbgWFUdr6oHgWuBncM7VNUNVfXNbvFGYGN3+5XA4ao6WVUPAIeB7ZMZXZLUR5/QbwDuHVpe6NadzqXAp8/k2CSXJZlPMr+4uNhjJElSXxN9MDbJG4A54F1nclxVXV1Vc1U1NzMzM8mRJOkJr0/oTwCbhpY3duseIcnLgbcDO6rq22dyrCRp5fQJ/RFgS5LNSc4BdgEHh3dIcgHwPgaRv39o0/XAxUnO7R6EvbhbJ0k6S9YvtUNVnUpyOYNArwP2V9XRJHuB+ao6yOBSzdOAv0gCcE9V7aiqk0neweCHBcDeqjq5Ip+JJGmsJUMPUFWHgEMj664Yuv3yxzh2P7D/8Q44bbN7PrXs+7h73yUTmESSHh9fGStJjTP0ktQ4Qy9JjTP0ktQ4Qy9Jjev1rBtN1nKfyeOzeCSdCc/oJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGre+z05JtgPvAdYB76+qfSPbLwTeDbwA2FVVB4a2fRe4rVu8p6p2TGBurbDZPZ9a1vF377tkQpNIWq4lQ59kHXAV8ApgATiS5GBV3T602z3AG4HfGXMX36qqFy5/VEnS49HnjH4bcKyqjgMkuRbYCfx/6Kvq7m7bQyswoyRpGfpco98A3Du0vNCt6+vJSeaT3JjkNeN2SHJZt8/84uLiGdy1JGkpZ+PB2GdX1Rzwy8C7kzxndIequrqq5qpqbmZm5iyMJElPHH1CfwLYNLS8sVvXS1Wd6P48DvwtcMEZzCdJWqY+oT8CbEmyOck5wC7gYJ87T3Jukid1t88HXsbQtX1J0spbMvRVdQq4HLgeuAO4rqqOJtmbZAdAkhcnWQBeB7wvydHu8J8E5pPcAtwA7Bt5to4kaYX1eh59VR0CDo2su2Lo9hEGl3RGj/sH4KeXOaMkaRl8ZawkNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNa5X6JNsT3JnkmNJ9ozZfmGSm5OcSvLakW27k9zVfeye1OCSpH6WDH2SdcBVwKuArcDrk2wd2e0e4I3AR0eOPQ+4EngJsA24Msm5yx9bktRXnzP6bcCxqjpeVQ8C1wI7h3eoqrur6lbgoZFjXwkcrqqTVfUAcBjYPoG5JUk99Qn9BuDeoeWFbl0fvY5NclmS+STzi4uLPe9aktTHqngwtqqurqq5qpqbmZmZ9jiS1JQ+oT8BbBpa3tit62M5x0qSJqBP6I8AW5JsTnIOsAs42PP+rwcuTnJu9yDsxd06SdJZsmToq+oUcDmDQN8BXFdVR5PsTbIDIMmLkywArwPel+Rod+xJ4B0MflgcAfZ26yRJZ8n6PjtV1SHg0Mi6K4ZuH2FwWWbcsfuB/cuYUZK0DKviwVhJ0sox9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY3rFfok25PcmeRYkj1jtj8pyce77Z9PMtutn03yrSRf7D7eO+H5JUlLWL/UDknWAVcBrwAWgCNJDlbV7UO7XQo8UFXPTbILeCfwS922L1fVCyc7tiSprz5n9NuAY1V1vKoeBK4Fdo7ssxP4UHf7AHBRkkxuTEnS49Un9BuAe4eWF7p1Y/epqlPAN4Bndts2J/lCkr9L8vPj/oIklyWZTzK/uLh4Rp+AJOmxrfSDsfcBP1pVFwBvAT6a5AdGd6qqq6tqrqrmZmZmVngkSXpi6RP6E8CmoeWN3bqx+yRZDzwd+HpVfbuqvg5QVTcBXwZ+fLlDS5L66xP6I8CWJJuTnAPsAg6O7HMQ2N3dfi3wuaqqJDPdg7kk+TFgC3B8MqNLkvpY8lk3VXUqyeXA9cA6YH9VHU2yF5ivqoPANcCHkxwDTjL4YQBwIbA3yXeAh4Bfr6qTK/GJSJLGWzL0AFV1CDg0su6Kodv/C7xuzHGfAD6xzBklScvgK2MlqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIa1yv0SbYnuTPJsSR7xmx/UpKPd9s/n2R2aNvbuvV3JnnlBGeXJPWwZOiTrAOuAl4FbAVen2TryG6XAg9U1XOBPwbe2R27FdgF/BSwHfjT7v4kSWdJnzP6bcCxqjpeVQ8C1wI7R/bZCXyou30AuChJuvXXVtW3q+orwLHu/iRJZ8n6HvtsAO4dWl4AXnK6farqVJJvAM/s1t84cuyG0b8gyWXAZd3ifye5s9f0j3Q+8LXHcdwjZ3nncu/hMe/TGSdjIjOuMGecjLUwI6yOOZ99ug19Qr/iqupq4Orl3EeS+aqam9BIK8IZJ8MZJ8MZJ2e1z9nn0s0JYNPQ8sZu3dh9kqwHng58veexkqQV1Cf0R4AtSTYnOYfBg6sHR/Y5COzubr8W+FxVVbd+V/esnM3AFuCfJjO6JKmPJS/ddNfcLweuB9YB+6vqaJK9wHxVHQSuAT6c5BhwksEPA7r9rgNuB04Bb6qq767Q57KsSz9niTNOhjNOhjNOzqqeM4MTb0lSq3xlrCQ1ztBLUuPWfOiXenuGaUuyKckNSW5PcjTJm6c90+kkWZfkC0n+etqznE6SZyQ5kOSfk9yR5GenPdOoJL/dfa2/lORjSZ68Cmban+T+JF8aWndeksNJ7ur+PHcVzviu7mt9a5K/TPKMKY44dsahbW9NUknOn8Zsj2VNh77n2zNM2yngrVW1FXgp8KZVOOPD3gzcMe0hlvAe4G+q6ieAn2GVzZtkA/CbwFxVPZ/BExh2TXcqAD7I4G1Ihu0BPltVW4DPdsvT9EEePeNh4PlV9QLgX4C3ne2hRnyQR89Ikk3AxcA9Z3ugPtZ06On39gxTVVX3VdXN3e3/YhCmR706eNqSbAQuAd4/7VlOJ8nTgQsZPMuLqnqwqv5jqkONtx74vu41JU8B/m3K81BVf8/gGXHDht+65EPAa87mTKPGzVhVn6mqU93ijQxeizM1p/l3hMF7fP0usCqf3bLWQz/u7RlWXUQf1r2r5wXA56c8yjjvZvAf9aEpz/FYNgOLwAe6S0zvT/LUaQ81rKpOAH/I4MzuPuAbVfWZ6U51Ws+qqvu6218FnjXNYXr4NeDT0x5iVJKdwImqumXas5zOWg/9mpHkacAngN+qqv+c9jzDkrwauL+qbpr2LEtYD7wI+LOqugD4H6Z/ueERuuvcOxn8UPoR4KlJ3jDdqZbWvcBxVZ6NAiR5O4PLoB+Z9izDkjwF+D3gimnP8ljWeujXxFssJPleBpH/SFV9ctrzjPEyYEeSuxlc/vqFJH8+3ZHGWgAWqurh34gOMAj/avJy4CtVtVhV3wE+CfzclGc6nX9P8sMA3Z/3T3mesZK8EXg18Cu1+l748xwGP9Rv6b5/NgI3J/mhqU41Yq2Hvs/bM0xV93bN1wB3VNUfTXuecarqbVW1sapmGfwbfq6qVt1ZaFV9Fbg3yfO6VRcxeNX1anIP8NIkT+m+9hexyh4wHjL81iW7gb+a4ixjJdnO4JLijqr65rTnGVVVt1XVD1bVbPf9swC8qPu/umqs6dB3D9I8/PYMdwDXVdXR6U71KC8DfpXBWfIXu49fnPZQa9hvAB9JcivwQuD3pzvOI3W/bRwAbgZuY/A9NvWXxyf5GPCPwPOSLCS5FNgHvCLJXQx+E9m3Cmf8E+D7gcPd9857V+GMq55vgSBJjVvTZ/SSpKUZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMb9H1Taz3zVY1eyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "ax.bar(np.arange(len(sires))+0.5, sc.family_size())\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Identifying fathers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the section on [paternityArray objects](https://fractional-analysis-of-paternity-and-sibships.readthedocs.io/en/latest/tutorials/03_paternity_arrays.html#paternityarray-structure) we used the method `prob_array()` to create an matrix of probabilities that each candidate is the father of each offspring, or that the father was missing from a `paternityArray` object. That was based only on comparing alleles shared between the mother, the candidate father and individual offspring.\n", "\n", "We can get a similar matrix of probabilities of paternity after accounting for information about paternity that is shared between siblings by calling `posterior_paternity_matrix` on the `sibshipCluster` object `sc`. For example, father 'a_1' is the true sire of the furst five individuals in this toy dataset. We can compare his probabilities of paternity for his real offspring to the other unrelated progeny. We see that there is an increase in support for him as the father of the first five individuals after accounting for information shared between siblings, and a decrease for the other progeny." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Before clusteringAfter clustering
08.842832e-011.000000e+00
11.000000e+001.000000e+00
28.630426e-019.999512e-01
39.696970e-011.000000e+00
47.852761e-019.498510e-01
53.906250e-724.088879e-273
61.367187e-834.088879e-273
78.750000e-604.088879e-273
88.750000e-604.088879e-273
91.070804e-832.226337e-88
107.567568e-831.814319e-84
111.296296e-712.695164e-76
121.869202e-727.826369e-132
133.926282e-607.826369e-132
141.361111e-821.361111e-82
\n", "
" ], "text/plain": [ " Before clustering After clustering\n", "0 8.842832e-01 1.000000e+00\n", "1 1.000000e+00 1.000000e+00\n", "2 8.630426e-01 9.999512e-01\n", "3 9.696970e-01 1.000000e+00\n", "4 7.852761e-01 9.498510e-01\n", "5 3.906250e-72 4.088879e-273\n", "6 1.367187e-83 4.088879e-273\n", "7 8.750000e-60 4.088879e-273\n", "8 8.750000e-60 4.088879e-273\n", "9 1.070804e-83 2.226337e-88\n", "10 7.567568e-83 1.814319e-84\n", "11 1.296296e-71 2.695164e-76\n", "12 1.869202e-72 7.826369e-132\n", "13 3.926282e-60 7.826369e-132\n", "14 1.361111e-82 1.361111e-82" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(\n", " pd.DataFrame({\n", " \"Before clustering\" : patlik.prob_array()[:, 1],\n", " \"After clustering\" : sc.posterior_paternity_matrix()[:, 1]\n", " })\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since `posterior_paternity_matrix` returns a matrix, some wrangling is required to get the output into a format that is useful for a human to read. The following sections describe two helper functions that use `posterior_paternity_matrix` to create (1) a summary of probable mating events and (2) a summary of the paternity of each offspring." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mating events" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We very frequently want to know who the fathers of the offspring were to say something about mating events. There are several levels of complexity. Firstly, you can use the `sires` method to return a list of all the males who could possibly have mated with the mother. This is essentially identifying **mating events**, but doesn't say anything about the paternity of individual offspring. For many applications, that may be all you need because it's the mating events that are the unit of interest, not the number of offspring per se.\n", "\n", "Once you have a `sibshipCluster` object, doing this is easy:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
positionlabellog_probproboffspring
01a_10.000000e+001.000000e+004.949802e+00
12a_20.000000e+001.000000e+004.000000e+00
23a_30.000000e+001.000000e+002.996754e+00
34a_4-1.110223e-161.000000e+002.000000e+00
45a_5-1.151638e-018.912201e-018.888889e-01
587a_87-2.218429e+001.087799e-011.111111e-01
6152a_152-1.078098e+012.079126e-051.163148e-06
7254a_254-2.209107e+001.097987e-012.292585e-02
8257a_257-3.677553e+002.528478e-025.731309e-03
9376a_376-2.908523e+005.455628e-021.146262e-02
10388a_388-3.677656e+002.528218e-025.731309e-03
11631a_631-5.118017e+005.987884e-031.432977e-03
12639a_639-4.750292e+008.649166e-031.295942e-03
13668a_668-7.947764e+003.534515e-043.841083e-05
14680a_680-4.542409e+001.064773e-021.295942e-03
15693a_693-4.370803e+001.264109e-022.865655e-03
16852a_852-1.286042e+012.598908e-066.001692e-07
17963a_963-5.364877e+004.678034e-036.479710e-04
181000nan-1.977494e+021.313807e-862.826406e-10
\n", "
" ], "text/plain": [ " position label log_prob prob offspring\n", "0 1 a_1 0.000000e+00 1.000000e+00 4.949802e+00\n", "1 2 a_2 0.000000e+00 1.000000e+00 4.000000e+00\n", "2 3 a_3 0.000000e+00 1.000000e+00 2.996754e+00\n", "3 4 a_4 -1.110223e-16 1.000000e+00 2.000000e+00\n", "4 5 a_5 -1.151638e-01 8.912201e-01 8.888889e-01\n", "5 87 a_87 -2.218429e+00 1.087799e-01 1.111111e-01\n", "6 152 a_152 -1.078098e+01 2.079126e-05 1.163148e-06\n", "7 254 a_254 -2.209107e+00 1.097987e-01 2.292585e-02\n", "8 257 a_257 -3.677553e+00 2.528478e-02 5.731309e-03\n", "9 376 a_376 -2.908523e+00 5.455628e-02 1.146262e-02\n", "10 388 a_388 -3.677656e+00 2.528218e-02 5.731309e-03\n", "11 631 a_631 -5.118017e+00 5.987884e-03 1.432977e-03\n", "12 639 a_639 -4.750292e+00 8.649166e-03 1.295942e-03\n", "13 668 a_668 -7.947764e+00 3.534515e-04 3.841083e-05\n", "14 680 a_680 -4.542409e+00 1.064773e-02 1.295942e-03\n", "15 693 a_693 -4.370803e+00 1.264109e-02 2.865655e-03\n", "16 852 a_852 -1.286042e+01 2.598908e-06 6.001692e-07\n", "17 963 a_963 -5.364877e+00 4.678034e-03 6.479710e-04\n", "18 1000 nan -1.977494e+02 1.313807e-86 2.826406e-10" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.sires()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The columns in the output tell you several bits of information. The most interesting of these are:\n", "\n", "1. **label** is the name of the candidate father\n", "2. **prob** is the probability that the male sired at least one offspring with the mother, as a weighted average over partition structures. For example, \n", "3. **offspring** shows the expected number of offspring sired by the male, as a weighted average over partition structures. Specifically, it's the sum over rows from `posterior_paternity_matrix`; see below.\n", "\n", "Note that if you have multiple maternal families the output will look a bit different. See the [tutorial on multiple maternal families](https://fractional-analysis-of-paternity-and-sibships.readthedocs.io/en/latest/tutorials/07_dealing_with_multiple_half-sib_families.html#clustering-multiple-families) for more.\n", "\n", "We can check this table makes sense by reviewing who the real fathers really are. This snippet gives a list of the names of the five true fathers, followed by the number of offspring sired by each." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array(['a_1', 'a_2', 'a_3', 'a_4', 'a_5'], dtype='\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
progenycandidate_1logprob_1candidate_2logprob_2candidate_3logprob_3candidate_4logprob_4
0offs_0a_1-1.229779e-01a_839-2.202419a_852-5.668155a_631-6.361302
1offs_1a_1-1.199290e-10missing-22.852463a_14-28.101025a_773-29.710463
2offs_2a_1-1.472913e-01a_668-2.226733a_839-3.613027a_852-6.385616
3offs_3a_1-3.077166e-02a_343-3.496508missing-21.443367a_605-27.661793
4offs_4a_1-2.417199e-01a_254-2.321161a_376-3.014309a_388-3.707456
5offs_5a_2-1.355254e-10missing-22.732634a_329-28.457700a_951-29.017315
6offs_6a_2-6.884449e-11missing-23.405092a_378-28.794172a_368-30.403610
7offs_7a_2-4.156320e-11missing-23.909696a_813-30.180466a_712-30.403610
8offs_8a_2-7.621637e-11missing-23.310463a_951-29.017315a_792-29.150847
9offs_9a_3-1.108144e-01a_152-2.883403a_159-3.576550a_368-4.269697
10offs_10a_3-1.451820e-01a_680-2.917771a_639-2.917771a_963-3.610918
11offs_11a_3-5.324451e-02a_254-3.518980a_711-4.212128a_950-4.905275
12offs_12a_4-6.866347e-10missing-21.107519a_39-27.071405a_382-27.407878
13offs_13a_4-6.429435e-02a_604-2.836883a_27-5.609472missing-21.792526
14offs_14a_5-1.177830e-01a_87-2.197225missing-23.714176a_852-32.377691
\n", "" ], "text/plain": [ " progeny candidate_1 logprob_1 candidate_2 logprob_2 candidate_3 \\\n", "0 offs_0 a_1 -1.229779e-01 a_839 -2.202419 a_852 \n", "1 offs_1 a_1 -1.199290e-10 missing -22.852463 a_14 \n", "2 offs_2 a_1 -1.472913e-01 a_668 -2.226733 a_839 \n", "3 offs_3 a_1 -3.077166e-02 a_343 -3.496508 missing \n", "4 offs_4 a_1 -2.417199e-01 a_254 -2.321161 a_376 \n", "5 offs_5 a_2 -1.355254e-10 missing -22.732634 a_329 \n", "6 offs_6 a_2 -6.884449e-11 missing -23.405092 a_378 \n", "7 offs_7 a_2 -4.156320e-11 missing -23.909696 a_813 \n", "8 offs_8 a_2 -7.621637e-11 missing -23.310463 a_951 \n", "9 offs_9 a_3 -1.108144e-01 a_152 -2.883403 a_159 \n", "10 offs_10 a_3 -1.451820e-01 a_680 -2.917771 a_639 \n", "11 offs_11 a_3 -5.324451e-02 a_254 -3.518980 a_711 \n", "12 offs_12 a_4 -6.866347e-10 missing -21.107519 a_39 \n", "13 offs_13 a_4 -6.429435e-02 a_604 -2.836883 a_27 \n", "14 offs_14 a_5 -1.177830e-01 a_87 -2.197225 missing \n", "\n", " logprob_3 candidate_4 logprob_4 \n", "0 -5.668155 a_631 -6.361302 \n", "1 -28.101025 a_773 -29.710463 \n", "2 -3.613027 a_852 -6.385616 \n", "3 -21.443367 a_605 -27.661793 \n", "4 -3.014309 a_388 -3.707456 \n", "5 -28.457700 a_951 -29.017315 \n", "6 -28.794172 a_368 -30.403610 \n", "7 -30.180466 a_712 -30.403610 \n", "8 -29.017315 a_792 -29.150847 \n", "9 -3.576550 a_368 -4.269697 \n", "10 -2.917771 a_963 -3.610918 \n", "11 -4.212128 a_950 -4.905275 \n", "12 -27.071405 a_382 -27.407878 \n", "13 -5.609472 missing -21.792526 \n", "14 -23.714176 a_852 -32.377691 " ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sc.paternity()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, the data are simulated, so we know who the real fathers are, and can print them with `progeny.fathers`. We see that the top candidate for all offspring is indeed the true father, with strong support (log probabilitites close to zero)." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['a_1', 'a_1', 'a_1', 'a_1', 'a_1', 'a_2', 'a_2', 'a_2', 'a_2',\n", " 'a_3', 'a_3', 'a_3', 'a_4', 'a_4', 'a_5'], dtype='