Prediction of selection in euploid embryos due to lethal mutations

.gitignore

*.png
*.pdf

denovo_gl_modeling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "959f35b9-95d0-4519-9d6c-7cd2b19e435b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom, geom\n",
    "from scipy.special import logsumexp\n",
    "from arjun_plot.utils import *\n",
    "from tqdm import tqdm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f7c90ee6-0079-4bf8-90be-b2adb20bafef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-6.238324625039508, -7.735129730686877, 1.4968051056473692)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "def logprob_read(b, a1,a2, eps=1e-3):\n",
    "    l1 = np.log(1/2) + np.log((1 - eps)*(b == a1) + (eps)*(b != a1))\n",
    "    l2 = np.log(1/2) + np.log((1 - eps)*(b == a2) + (eps)*(b != a2))\n",
    "    return logsumexp([l1,l2])\n",
    "    \n",
    "def geno_likelihoods_het(m=30, q=10, seed=42):\n",
    "    \"\"\"Calculate the GL for a site.\n",
    "    \n",
    "    NOTE: we are using the GATK likelihood here\n",
    "    http://www.popgen.dk/angsd/index.php/Genotype_Likelihoods#Theory\n",
    "    \"\"\"\n",
    "    assert q > 0\n",
    "    assert seed > 0\n",
    "    np.random.seed(seed)\n",
    "    eps = 10**(-q/10)\n",
    "    # Explicitly compare the likelihood \n",
    "    b_het = binom.rvs(n=1, p=0.5, size=m)\n",
    "    loglik_het01 = 0.0\n",
    "    loglik_het10 = 0.0\n",
    "    loglik_hom00a = 0.0\n",
    "    loglik_hom00b = 0.0\n",
    "    for i in range(m):\n",
    "        loglik_het01 += logprob_read(b_het[i], a1=0, a2=1, eps=eps)\n",
    "        loglik_het10 += logprob_read(b_het[i], a1=1, a2=0, eps=eps)\n",
    "        loglik_hom00a += logprob_read(b_het[i], a1=0, a2=0, eps=eps)\n",
    "        loglik_hom00b += logprob_read(b_het[i], a1=0, a2=0, eps=eps)\n",
    "    loglik_hom = logsumexp([loglik_hom00a, loglik_hom00b])\n",
    "    loglik_het = logsumexp([loglik_het01, loglik_het10])\n",
    "    return loglik_het, loglik_hom, loglik_het-loglik_hom \n",
    "\n",
    "geno_likelihoods_het(m=10, q=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bcd88cd8-0034-4790-ae5a-f7b8598aef1b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████| 5/5 [00:24<00:00,  4.89s/it]\n",
      "100%|██████████████████████████████████████████████████████████████████████| 5/5 [00:27<00:00,  5.48s/it]\n",
      "100%|██████████████████████████████████████████████████████████████████████| 5/5 [00:30<00:00,  6.02s/it]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAEqCAYAAAA72HsuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5yUlEQVR4nO3deVxU1f8/8NewDKuICLIkIIILBoKpEOQuiWSlldJiaZpbkrl87Fv++hhp9e1Tllk5Lp9KqTS/LqX2SVMRJbOwEDVFCLFGpIgRUEFGWRzO7w8+TI6sszEXeD0fj3nkPffMue8zd0bfnXvuuTIhhAARERGRBFlZOgAiIiKixjBRISIiIsliokJERESSxUSFiIiIJIuJChEREUkWExUiIiKSLCYqREREJFlMVIiIiEiymKgYQQiBsrIycM08IiIi82CiYoRr166hc+fOuHbtmqVDISIiapeYqBAREZFkMVEhIiIiyWKiQkRERJLFRIWIiIgki4kKERERSZaNpQMgIiKyFI1Gg+rqakuH0S5ZW1vDxsYGMpnMqHaYqEiJWg04O9f+ubwccHKybDxERO1YeXk5/vjjD66FZUaOjo7w9vaGXC43uA0mKlKhVte+bt0GmKwQEZmBRqPBH3/8AUdHR3h4eBj9f/2kSwiBqqoqFBUVQalUolevXrCyMmy2CRMVAygUCigUCmg0GtM1WjeSUsfTs/a/zPSJiEyuuroaQgh4eHjAwcHB0uG0Sw4ODrC1tUVeXh6qqqpgb29vUDucTGuAhIQEZGVlIT093dKhEBGRETiSYl6GjqLotGGCOMgUyssBlervbZWqtoyIiKgD46Ufqbh9LoqTE+enEBFRh8dERUqcnDgnhYiI6Ba89ENERNTG5OfnY/r06fDx8YFcLoe/vz/mz5+PkpISg9r75Zdf8Pjjj8PX1xcODg4IDg7G+++/X69eamoq7rrrLtjZ2SEoKAhJSUlG9qR5TFSIiIjakN9//x2DBg1Cbm4utmzZgvPnz2PdunVISUlBVFQULl++rHebGRkZ6NatGzZt2oSzZ8/i5ZdfxpIlS7B69WptHaVSiXHjxmHkyJE4deoUFixYgBkzZmD//v2m7F49MsGVbgxWVlaGzp07o7S0FC4uLpYOh4iIWqiiogJKpRIBAQEG3zZrKXFxccjMzMS5c+d0bq0uLCxEYGAgpkyZgrVr1xp9nISEBGRnZ+PQoUMAgBdffBF79uxBZmamts5jjz2Gq1evYt++fQ22YYrPmSMqREREhlKrAZms9nXrop1mcvnyZezfvx9z586tt/6Ll5cXJk+ejK1bt0IIgTlz5sDZ2bnJV1NKS0vh5uam3U5LS0NMTIxOndjYWKSlpZmugw3gZFoiIiJDWGBF8dzcXAghEBwc3OD+4OBgXLlyBUVFRVi+fDkWL15s0HF+/PFHbN26FXv27NGWFRYWwrNuMdL/8vT0RFlZGW7cuGG2hfOYqBARERnCgiuKNzdrQy6Xw9XVFd26ddO77czMTIwfPx6JiYkYM2aMoSGaDC/9EBERtRFBQUGQyWTIzs5ucH92djY8PDzg6upq0KWfrKwsjB49GrNmzcI///lPnX1eXl5Q3bowKQCVSgUXFxezPoaAIypERESGKC+vvdxTN5KiUpl9oc6uXbvi3nvvxZo1a7Bw4cJ6k2k3b96MhIQEAND70s/Zs2cxatQoTJ06FW+88Ua9/VFRUdi7d69OWXJyMqKiogzsTcvwrh8j8K4fIqK2yWR3/ajVf18CKi9vlRXFc3NzER0djeDgYLz++usICAjA2bNn8cILL8DGxgbff/99sxNlb5eZmYlRo0YhNjYWK1as0JZbW1vDw8MDQO3tySEhIUhISMD06dNx6NAhPP/889izZw9iY2MbbJd3/RAREVlS3YriQrTaY0969eqF9PR09OzZE/Hx8fD390dcXBx69+6NH374Qe8kBQB27NiBoqIibNq0Cd7e3trX4MGDtXUCAgKwZ88eJCcnIywsDO+++y4+/vjjRpMUU+GIihE4okJE1Da15XVUGpKYmIiVK1ciOTkZd999t6XD0TLF58w5KkRERG3csmXL0KNHDxw7dgwRERGwsmo/F0yYqBAREbUD06ZNs3QIZtF+Ui4iIiJqdzp8ovLQQw+hS5cumDhxoqVDISIiott0+ERl/vz5+OyzzywdBhERETWgwycqI0aMQKdOnSwdBhERETVA0onKkSNH8MADD8DHxwcymQy7du2qV0ehUKBHjx6wt7dHZGQkfv7559YPlIiIiMxC0omKWq1GWFgYFApFg/u3bt2KRYsWITExESdOnEBYWBhiY2Nx6dIlbZ3w8HCEhITUexUUFLRWN4iIiMhAkr49OS4uDnFxcY3uX7lyJWbOnKm9JWvdunXYs2cPNmzYgJdeegkAcOrUKZPFU1lZicrKSu12WVmZydomIiKi+iQ9otKUqqoqZGRkICYmRltmZWWFmJgYpKWlmeWYb775Jjp37qx9+fr6muU4REREVKvNJirFxcXQaDTwrHtq5X95enqisLCwxe3ExMRg0qRJ2Lt3L7p3795kkrNkyRKUlpZqX/n5+QbHT0REZKj8/HxMnz4dPj4+kMvl8Pf3x/z581FSUmJQeyUlJRg7dix8fHxgZ2cHX19fPPfcc/WuHKSmpuKuu+6CnZ0dgoKCkJSUZILeNE3Sl35aw8GDB1tc187ODnZ2dmaMhoiIqGm///47oqKi0Lt3b2zZskXn6cnffvstjh07Bjc3N73atLKywvjx4/H666/Dw8MD58+fR0JCAi5fvowvvvgCQO3Tk8eNG4c5c+Zg8+bNSElJwYwZM+Dt7W3WBxO22UTF3d0d1tbWUKlUOuUqlQpeXl5mPbZCoYBCoYBGozHrcYiIiG6XkJAAuVyOAwcOwMHBAQDg5+eHAQMGIDAwEC+//DLWrl2rV5tdunTBs88+q9329/fH3LlzsWLFCm3ZunXrEBAQgHfffRcAEBwcjKNHj+K9994za6LSZi/9yOVyDBw4ECkpKdqympoapKSkICoqyqzHTkhIQFZWFtLT0816HCIikja1GpDJal9qtfmPd/nyZezfvx9z587VJil1vLy8MHnyZGzduhVCCMyZMwfOzs5NvhpTUFCAr776CsOHD9eWpaWl6cwLBYDY2FizzQutI+kRlfLycpw/f167rVQqcerUKbi5ucHPzw+LFi3C1KlTMWjQIERERGDVqlVQq9Xt9sFMREQkHWq1bnJS92cnJ/MdMzc3F0IIBAcHN7g/ODgYV65cQVFREZYvX47Fixfr1f7jjz+O3bt348aNG3jggQfw8ccfa/cVFhY2OC+0rKwMN27cqJc4mYqkE5Xjx49j5MiR2u1FixYBAKZOnYqkpCQ8+uijKCoqwiuvvILCwkKEh4dj37599T5IIiIiU7t9QKLunx4hzH9s0cxB5HI5XF1d0a1bN73afe+995CYmIhz585hyZIlWLRoEdasWWNMqEaTdKIyYsSIZk/Gc889h+eee66VIiIiIrKcoKAgyGQyZGdn46GHHqq3Pzs7Gx4eHnB1dcWcOXOwadOmJtsrLy/X2fby8oKXlxf69u0LNzc3DB06FEuXLoW3tze8vLwanBfq4uJittEUQOKJilRxMi0REZWX117uqRtJUanMe9kHALp27Yp7770Xa9aswcKFC3UShMLCQmzevBkJCQkAYNCln1vV1NQAgHah06ioKOzdu1enTnJystnnhcpEc0MW1KiysjJ07twZpaWlcHFxsXQ4RETUQhUVFVAqlQgICIC9vb3B7ajVf18CKi83f6IC1M5TiY6ORnBwMF5//XWd25NtbGzw/fffNzlRtiF79+6FSqXC4MGD4ezsrG3Pzc0NR48eBVA7TzQkJAQJCQmYPn06Dh06hOeffx579uxp9K4fU3zObfauHyIiIktzcqqdkyJE6yQpANCrVy+kp6ejZ8+eiI+Ph7+/P+Li4tC7d2/88MMPeicpAODg4ICPPvoIQ4YMQXBwMBYuXIgHH3wQ33zzjbZOQEAA9uzZg+TkZISFheHdd9/Fxx9/bNZbkwGOqBiFIypERG2TqUZUpCIxMRErV65EcnIy7r77bkuHo2WKz5lzVAzAOSpERCQly5YtQ48ePXDs2DFERETAyqr9XDDhiIoROKJCRNQ2tbcRFaniHBUiIiJq15ioEBERkWQxUSEiIiLJYqJCREREksVExQAKhQL9+vXD4MGDLR0KERFRu6Z3opKVlYUpU6Zg8ODBiIuLw6efftrg83g2b94Ma2trkwQpNQkJCcjKykJ6erqlQyEiImrX9EpUcnNzERkZie3bt0MIgczMTEybNg3Dhg1DYWGhuWIkIiKiDkqvROWf//wnnJ2dcebMGRw/fhz5+fn47LPPcObMGURFRSEnJ8dccRIREVEHpFeicuzYMcybNw9BQUHasieffBLHjh2DlZUVhgwZgp9//tnkQRIREdHf8vPzMX36dPj4+EAul8Pf3x/z589HSUmJ0W2XlJSge/fukMlkuHr1qs6+1NRU3HXXXbCzs0NQUBCSkpKMPl5z9EpUSkpK4OXlVa+8b9+++PHHH9G9e3eMHj0a+/fvN1mARERE9Lfff/8dgwYNQm5uLrZs2YLz589j3bp1SElJQVRUFC5fvmxU+8888wz69+9fr1ypVGLcuHEYOXIkTp06hQULFmDGjBlm/zdfr0SlR48eOH36dIP7PD098d1332HAgAF48MEHsX37dpMEKEW864eIiAAAN9XAF7La1011qxwyISEBcrkcBw4cwPDhw+Hn54e4uDgcPHgQf/75J15++WWD2167di2uXr2KxYsX19u3bt06BAQE4N1330VwcDCee+45TJw4Ee+9954x3WmWXonKiBEjsH37dty8ebPB/S4uLkhOTsbYsWPx9ddfmyRAKeJdP0REhJtq3eTk9m0zuHz5Mvbv34+5c+fCwcFBZ5+XlxcmT56MrVu3QgiBOXPmwNnZucnXrbKysrB8+XJ89tlnDT7UMC0tDTExMTplsbGxSEtLM31Hb6HX05OffvppqFQqHD9+vNHHSNvZ2WHnzp1YtGgRfvnlF5MESUREJDnbdP+hx1eetf99wnzP+s3NzYUQAsHBwQ3uDw4OxpUrV1BUVITly5c3ODLSkMrKSjz++ONYsWIF/Pz88Pvvv9erU1hYCE9PT50yT09PlJWV4caNG/USJ1PRK1EZNGhQiy7pWFlZYdWqVYbGRERERE1oaP2yW8nlcri6uqJbt24tam/JkiUIDg7Gk08+aYrwTIor0xIRERkivhx4WPX39sOq2jIzCgoKgkwmQ3Z2doP7s7Oz4eHhAVdXV70u/Rw6dAjbt2+HjY0NbGxsMHr0aACAu7s7EhMTAdReWlKpVDrHU6lUcHFxMdtoCqDniAoRERH9l41T/e3by0ysa9euuPfee7FmzRosXLhQJ0EoLCzE5s2bkZCQAAB6Xfr58ssvcePGDe12eno6pk+fju+//x6BgYEAgKioKOzdu1fnfcnJyYiKijK2W02SiebGj6hRZWVl6Ny5M0pLS+Hi4mLpcIiIqIUqKiqgVCoREBAAe3t7S4ejl9zcXERHRyM4OBivv/46AgICcPbsWbzwwguwsbHB999/X2+irL5SU1MxcuRIXLlyBa6urgBqb08OCQlBQkICpk+fjkOHDuH555/Hnj17EBsb22A7pviceemHiIioDenVqxfS09PRs2dPxMfHw9/fH3Fxcejduzd++OEHo5OUxgQEBGDPnj1ITk5GWFgY3n33XXz88ceNJimmYtCISmlpKfLy8tCrV68Gr0uVlJTgzz//RHBwMGxtbU0SqJQoFAooFApoNBqcO3eOIypERG1MWx5RaUhiYiJWrlyJ5OTkRu/KtQSLjaiUlZVh4MCBUCgUDe6fO3cu4uLi+PRkIiKiVrBs2TJ88MEHOHbsGGpqaiwdjkkZlKj4+vpi9OjR2LBhQ719V65cwddff42pU6c2uGAMERERmd60adOwYMGCdvdvr8G9mTFjBnJycvDTTz/plH/xxReoqqrCM888Y3RwRERE1LEZnKiMHz8eXbt2rTeqsnHjRgwbNkx7OxMRERGRoQxOVGxtbfHkk09i69atqKioAACcOXMGJ06c4GgKERG1CVyhw7xM8fkadSHrmWeeQVlZGXbs2AGgdjTF1dUVkyZNMjowIiIic6m72aOqqsrCkbRv169fBwCj7gA2esG3qKgoODo64sCBA/D29kZ8fDxWr15tTJNtBhd8IyJqm4QQuHjxIqqrq+Hj49PuJqBamhAC169fx6VLl+Dq6gpvb2+D2zJ6Cf1nnnkGc+bMwfvvv4+SkhLMmDHD2CY7LHWpGk57ahfqUY8rh1Nn8y7FTETUUclkMnh7e0OpVCIvL8/S4bRbrq6u8PLyMqoNo0dUysvL4e3tjaqqKoSGhuL48eNGBdSWmHJERV2qxvVrangcqX2EdtEwFRw7OTFZISIyo5qaGl7+MRNbW1uTrKdm9IiKs7MzJk2ahKSkpA4zifbWlWlNxWmPM25NSeoSFjzBiV5EROZiZWXVLlambc9M8lBClUqFX3/9FYMHD4ajo6Mp4moTTDpH5QtZw+VMVIiIqAMzekQFADw9PeHp6WmKpjos9bjyhi/9WDguIiIiSzJJokLGu30uCuenEBERMVGRFKfOTtpLPUxRiIiIjFzwjYiIiMicmKgQERGRZJk8UYmMjDR1k0RERNRBmSxRKS4uBgBkZ2c3uD86OtpUhyIiIqIOwmSTaQMCAuDk5ITKykosWLAAoaGh6N+/P0JCQlBZWYnMzExTHYqIiIg6CJONqFy7dg0//vgjrKys4OXlhZSUFEybNk27zv8jjzxiqkMRERFRB2HS25N79uyJkydPom/fvtqy6upqXLt2DW5ubqY8FBEREXUAJl9HpVevXvjoo4+QkZEBV1dX/Otf/2KSQkRERAYx+V0/s2bNwpYtWzBgwACsXbsWAHDgwAG89957pj4UERERtXMmeSjhrdzd3ZGXlwcnJyd06dIFV65cQUlJCe6++27k5uaa8lAWc+vTk8+dO2eahxISERFRPSZPVPr06YMjR47A09NTm6gAgIuLC8rKykx5KIsz6dOTiYiIqB6TX/qZPXs2Jk+ejEuXLkEmkwEAcnJy4OrqaupDERERUTtn8sm0CxcuxOXLlxEUFISKigrMmjUL+/btw7x580x9KCIiImrnTH7pp05hYSEOHjyIy5cvIzw8HMOGDTPHYSyKl36IiIjMy+QjKkqlEhs2bIAQAn379sWIESPQr18/Ux+GiIiIOgCTz1F54IEHcO7cOQDAl19+iYceegjOzs4ICwsz9aGIiIionTP5iMqff/6JM2fOaCfSAkB5eTmf9UNERER6M8mISlFRkfbPEyZMQHp6us5+Z2dn3H333aY4FBEREXUgJklUevbsCXd3dwwdOhTXr1/HxIkT8fnnn6O4uNgUzRMREVEHZbK7fpRKJTIzM3VeOTk56Nq1K/r374/+/fvjrbfeMsWhJIN3/RAREZmXyeaoBAQEICAgAA888IC2rG6J+brEhYiIiEgfZltHpSPgiAoREZF5mfz2ZCIiIiJT0StRqa6uxssvv4ygoCB4eHhg3LhxOHnyZL16OTk5WLlyJe69916TBUpEREQdj16XfpYtW4Zly5bBy8sLPj4+yMrKgo2NDY4ePYrOnTtjzZo1+PLLL6FUKiGEgIuLC65evWrG8C2Ll346HvUlNZw9nQAA5So1nLo5WTgiIqL2Ta9EpVevXujZsye++eYb2NraoqCgAOPGjYObmxsyMjJQVlaGQYMGITY2FmPGjEFUVBRsbEy+ppxkMFHpWNSX1FAXXYdniAcAQJVZBCcPRyYrRERmpFeiIpfLsWbNGsyYMUNbtn//fsTFxSEkJATbtm1D3759zRKoFDFR6VhuWWxZB6ejExGZj15zVG7evAlHR0edsrpn+CxZsqRDJSlERERkfnrf9VNcXIxbB2FsbW0BAJ6enqaLqpXk5+drn+7cv39/bN++3dIhkYSVq9RQZf79uAhVZhHKVWoLRkRE1P7pPYFk4cKFWLJkCe68806Eh4ejR48ekMlkqK6uNkd8ZmVjY4NVq1YhPDwchYWFGDhwIO677z44OXHOAdV3+1wUzk8hIjI/veaoHDhwAL/88gtOnz6NX375BTk5OdoERSaTwd/fH6GhoQgJCdH+NyQkxGzBm1pYWBi++eYb+Pr6tqg+56gQERGZl1Er01ZXVyMrK0ubuJw+fRqnT5/GpUuXahuXyaDRaAwO7siRI1ixYgUyMjLw119/YefOnZgwYYJOHYVCgRUrVqCwsBBhYWH48MMPERERofexMjIyMHXqVL2W+meiQkREZF5G3Ttsa2uLsLAwhIWF4amnntKWq1QqnDp1CmfOnDEqOLVajbCwMEyfPh0PP/xwvf1bt27FokWLsG7dOkRGRmLVqlWIjY1FTk4OunXrBgAIDw/HzZs36733wIED8PHxAQBcvnwZU6ZMwUcffWRUvERERGRaZnvWT2FhIaytreHh4WGS9mQyWb0RlcjISAwePBirV68GANTU1MDX1xfz5s3DSy+91KJ2Kysrce+992LmzJk6yVZjdSsrK7XbZWVl8PX15YgKERGRmZjtWT+ffvopvL29zdU8qqqqkJGRgZiYGG2ZlZUVYmJikJaW1qI2hBB4+umnMWrUqGaTFAB488030blzZ+2rpXNZiIiIyDBmfSihOR/MXFxcDI1GU++2aE9PTxQWFraojR9++AFbt27Frl27EB4ejvDw8CYvVy1ZsgSlpaXaV35+vlF9ICIioqa13/XtW2DIkCGoqalpcX07OzvY2dmZMSIiIiK6lV4jKtu2bcOvv/5q1pGSlnJ3d4e1tTVUKpVOuUqlgpeXl1mPrVAo0K9fPwwePNisxyEiIuro9EpUHnvsMdx5551wcnLCoEGDMH36dKxatQqHDh1CSUmJuWJskFwux8CBA5GSkqItq6mpQUpKCqKiosx67ISEBGRlZSE9Pd2sxyEiIuro9Lr085///Ee7Vsrp06fx+eefQ6PRQPbfp7V5eXmhf//+6N+/P3Jzc40Orry8HOfPn9duK5VKnDp1Cm5ubvDz88OiRYswdepUDBo0CBEREVi1ahXUajWmTZtm9LGJiIjI8oy6Pbmqqgpnz57VWfDtzJkzKCqqfR6KsQu+paamYuTIkfXKp06diqSkJADA6tWrtQu+hYeH44MPPkBkZKTBx9QHF3wjIiIyL7Oso1JYWKgddVm8eLGpm5cMJipERETmZZa7fuomsw4YMMAczVucQqGAQqEwarSIiIiImme2lWnfeustvPzyyw0uX99ecESFiIjIvNrsgm9ERETU/pk1USEiIiIyRptd8M2SuOAbERFR69BrjoqVlRVkMhns7OzQr18/7Zop/fv3R1hYGLp27aqt+9Zbb+H//b//164nnHKOChERkXlJesE3IiIi6tgkveCb1HFEhYiIyLyMWkdFLpdjwIAB9dZLuXXBNyIiIiJDmW0dlY6AIypERETmxduTDcC7foiIiFoHR1SMwBEVIiIi8+KIChEREUkWExUiIiKSLCYqREREJFlMVIiIiEiyjFpHpTElJSVQKBSQyWRYunSpOQ5hUQqFAgqFol0vZkdERCQFZrnrJycnB8HBwVyZloiIiIxilhEVPz8/HD582BxNExERUQfCdVSMwBEVIiIi8zJ6RCUzMxN79+7FhQsXAAA9evRAXFwcQkNDjW2aiIiIOjiDR1QqKysxe/ZsfP755xBCwMqq9gaimpoayGQyTJ48GR9//DHkcrlJA5YSjqgQSZRaDTg71/65vBxwcrJsPERkMINvT37xxRfx2Wef4dlnn0V2djYqKipQWVmJ7OxszJkzB5s2bcL//M//mDJWIqLmqdW1r8a2iahNMXhExd3dHePGjcOnn37a4P6nnnoK3377LYqLi40KUMo4okIkQTJZw+WcjkfUJhk8olJdXY2777670f3R0dG4efOmoc1LGp+eTERE1DoMHlGJj49HVVUVdu3a1eD+8ePHw87ODtu2bTMmPknjiAqRBNVd6vH0rN1WqWrnqHCeClGbZHCikpOTg/j4eAQGBiIhIQFBQUEAgNzcXCgUCiiVSmzduhUeHh4673NzczM+aolgokIkUZxMS9RuGJyo1N3lAwCy264J1zV5ezmAdrVSLRMVIiIi8zJ4HZVXXnmlwUSEiIiIyFS4Mq0ROKJCRERkXgbf9UNERERkbiZ5KGF5eTny8/MBAL6+vnCum8RGREREZASjRlTS09MxcuRIdOnSBSEhIQgJCUGXLl0watQoHD9+3FQxEhERUQdl8ByVn376CSNGjIBcLscTTzyB4OBgAEB2dja2bNmCqqoqpKamIiIiwqQBSwnnqBAREZmXwYlKTEwMLly4gKNHj8LLy0tnn0qlwj333IOAgAAkJyebJFApUSgUUCgU0Gg0OHfuHBMVIiIiMzH40s9PP/2E2bNn10tSAMDT0xOzZs3CsWPHjApOqhISEpCVlYX09HRLh0JERNSuGZyoWFlZNfksH41Go7MoHBEREZG+DM4koqOjoVAokJeXV2/fxYsXsWbNGtxzzz1GBUdEREQdm8FzVE6ePImhQ4dCo9HgoYceQu/evQHUPgNo9+7dsLGxwffff4+wsDCTBiwlnExLRERkXkatTJuVlYWXX34ZycnJuH79OgDA0dERY8aMweuvv45+/fqZLFApYqJCRERkXgYt+Hb9+nUMHToUM2fOxM6dO1FTU4OioiIAgIeHB+emEBERkUkYlKg4OjpCqVRqH0poZWUFT09PkwZGREREZPDQx9ixY7F//35TxkJERESkw+A5KtnZ2Zg0aRIGDBiA2bNnIyAgAA4ODvXqubm5GR2kVHGOChERkXkZnKjcOg+l7hJQQzQajSHNtwlMVIiIiMzL4Kcnv/LKK00mKERERETGMur25I6OIypERETmpfeISkVFBXbv3g2lUomuXbvi/vvvh7e3tzliIyIiog5Or0Tl0qVLiI6OhlKpRN1AjKOjI3bt2oWYmBizBChFtz49mYiIiMxHr0s/8+bNw9q1a7FgwQKMGjUK58+fx2uvvQYXFxf89ttv5oxTknjph4iIyLz0GlE5cOAApkyZgnfeeUdb5unpiSeeeAI5OTno06ePyQMkIiKijkuvBd8uXryIIUOG6JQNGTIEQgioVCqTBkZERESkV6JSWVkJe3t7nbK67Zs3b5ouKiIiIiIYcNfPhQsXcOLECe12aWkpACA3Nxeurq716t91112GR0dEREQdml6Taa2srBpc5E0IUa+8rqw93xnDybRERETmpdeIysaNG80VBxEREVE9XJnWCBxRISIiMi+9JtMSERERtSYmKkRERCRZTFSIiIhIspioEBERkWQxUSEiIiLJYqJCREREksVEhYiIiCSLiQoRERFJFhMVIiIikiwmKkRERCRZHTpRuXr1KgYNGoTw8HCEhITgo48+snRIREREdIsO/awfjUaDyspKODo6Qq1WIyQkBMePH0fXrl1b9H4+64eIiMi8OvSIirW1NRwdHQEAlZWVEEKgA+dtREREkiPpROXIkSN44IEH4OPjA5lMhl27dtWro1Ao0KNHD9jb2yMyMhI///yzXse4evUqwsLC0L17d7zwwgtwd3c3UfRERERkLEknKmq1GmFhYVAoFA3u37p1KxYtWoTExEScOHECYWFhiI2NxaVLl7R16uaf3P4qKCgAALi6uuKXX36BUqnEF198AZVK1Sp9IyIioua1mTkqMpkMO3fuxIQJE7RlkZGRGDx4MFavXg0AqKmpga+vL+bNm4eXXnpJ72PMnTsXo0aNwsSJExvcX1lZicrKSu12WVkZfH19OUeFiIjITCQ9otKUqqoqZGRkICYmRltmZWWFmJgYpKWltagNlUqFa9euAQBKS0tx5MgR9OnTp9H6b775Jjp37qx9+fr6GtcJIiIialKbTVSKi4uh0Wjg6empU+7p6YnCwsIWtZGXl4ehQ4ciLCwMQ4cOxbx58xAaGtpo/SVLlqC0tFT7ys/PN6oPRERE1DQbSwdgSRERETh16lSL69vZ2cHOzs58AREREZGONjui4u7uDmtr63qTX1UqFby8vCwUFREREZlSm01U5HI5Bg4ciJSUFG1ZTU0NUlJSEBUVZdZjKxQK9OvXD4MHDzbrcYiIiDo6SV/6KS8vx/nz57XbSqUSp06dgpubG/z8/LBo0SJMnToVgwYNQkREBFatWgW1Wo1p06aZNa6EhAQkJCRoV6YlIiIi85B0onL8+HGMHDlSu71o0SIAwNSpU5GUlIRHH30URUVFeOWVV1BYWIjw8HDs27ev3gRbIiIiapvazDoqUsRn/RAREZlXm52jYkmco0JERNQ6OKJiBI6oEBERmRdHVIiIiEiymKgQERGRZDFRMQDnqBAREbUOzlExAueoEBERmRdHVIiIiEiymKgQERGRZDFRISIiIsliokJERESSxUTFALzrh4iIqHXwrh8j8K4fIiIi8+KIChEREUkWExUifdxUA1/Ial831ZaOhhqhLv37PKlLeZ4kS60GZLLal5rniRrGRIWopW6qdZOT27dJEtSlaly/9vd5uX5NzWRFitRq3eTk9m2i/+IcFSNwjkoH84Ws4fIn+BOSFJ6ntkHWyHniP0l0G46oGIB3/RAREbUOjqgYgSMqHUzdpZ6vPGu3H1YBNk61L5KMuks/Hkdqz1PRMBUcOznBqTPPk6TUXerx/O/vSaUCnJxqX0S3sLF0AERtxu0JCZMUSbo9IWGSIlG3JyRMUqgRHFExAkdUiIiIzItzVIiIiEiymKgQERGRZDFRISIiIsliokJERESSxUTFAFxHhYiIqHXwrh8j8K4fIiIi8+I6Kkaoy/HKysosHAkREVHb1KlTJ8gae6QCmKgY5dq1awAAX19fC0dCRETUNjV3VYKXfoxQU1ODgoKCZrNBfZSVlcHX1xf5+fnt5nJSe+tTe+sPwD61FexT29De+mTu/nBExYysrKzQvXt3s7Tt4uLSLr7gt2pvfWpv/QHYp7aCfWob2lufLNUf3vVDREREksVEhYiIiCSLiYrE2NnZITExEXZ2dpYOxWTaW5/aW38A9qmtYJ/ahvbWJ0v3h5NpiYiISLI4okJERESSxUSFiIiIJIuJChEREUkWExUiIiKSLCYqrejVV1+FTCbTefXt27fJ92zfvh19+/aFvb09QkNDsXfv3laKtmV69OhRr08ymQwJCQkN1k9KSqpX197evpWj1nXkyBE88MAD8PHxgUwmw65du3T2CyHwyiuvwNvbGw4ODoiJiUFubm6z7SoUCvTo0QP29vaIjIzEzz//bKYe6GqqP9XV1XjxxRcRGhoKJycn+Pj4YMqUKSgoKGiyTUO+u6bU3Dl6+umn68U3duzYZtu11DkCmu9TQ78rmUyGFStWNNqmJc/Tm2++icGDB6NTp07o1q0bJkyYgJycHJ06FRUVSEhIQNeuXeHs7IxHHnkEKpWqyXYN/f2ZQnN9unz5MubNm4c+ffrAwcEBfn5+eP7551FaWtpku4Z+X02hJedpxIgR9eKbM2dOk+2a8zwxUWlld955J/766y/t6+jRo43W/fHHH/H444/jmWeewcmTJzFhwgRMmDABmZmZrRhx09LT03X6k5ycDACYNGlSo+9xcXHReU9eXl5rhdsgtVqNsLAwKBSKBve//fbb+OCDD7Bu3Tr89NNPcHJyQmxsLCoqKhptc+vWrVi0aBESExNx4sQJhIWFITY2FpcuXTJXN7Sa6s/169dx4sQJLF26FCdOnMBXX32FnJwcPPjgg822q89319SaO0cAMHbsWJ34tmzZ0mSbljxHQPN9urUvf/31FzZs2ACZTIZHHnmkyXYtdZ6+++47JCQk4NixY0hOTkZ1dTXGjBkDtVqtrbNw4UL85z//wfbt2/Hdd9+hoKAADz/8cJPtGvL7M5Xm+lRQUICCggK88847yMzMRFJSEvbt24dnnnmm2bb1/b6aSkvOEwDMnDlTJ7633367yXbNep4EtZrExEQRFhbW4vrx8fFi3LhxOmWRkZFi9uzZJo7MdObPny8CAwNFTU1Ng/s3btwoOnfu3LpB6QGA2Llzp3a7pqZGeHl5iRUrVmjLrl69Kuzs7MSWLVsabSciIkIkJCRotzUajfDx8RFvvvmmWeJuzO39acjPP/8sAIi8vLxG6+j73TWnhvo0depUMX78eL3akco5EqJl52n8+PFi1KhRTdaR0nm6dOmSACC+++47IUTt78bW1lZs375dWyc7O1sAEGlpaQ22Yejvz1xu71NDtm3bJuRyuaiurm60jiHfV3NpqE/Dhw8X8+fPb3Eb5j5PHFFpZbm5ufDx8UHPnj0xefJkXLx4sdG6aWlpiImJ0SmLjY1FWlqaucM0SFVVFTZt2oTp06c3+YCp8vJy+Pv7w9fXF+PHj8fZs2dbMUr9KJVKFBYW6pyHzp07IzIystHzUFVVhYyMDJ33WFlZISYmRpLnrrS0FDKZDK6urk3W0+e7awmpqano1q0b+vTpg2effRYlJSWN1m1r50ilUmHPnj0t+j91qZynussfbm5uAICMjAxUV1frfOZ9+/aFn59fo5+5Ib8/c7q9T43VcXFxgY1N04/S0+f7ak6N9Wnz5s1wd3dHSEgIlixZguvXrzfahrnPExOVVhQZGakdGly7di2USiWGDh2Ka9euNVi/sLAQnp6eOmWenp4oLCxsjXD1tmvXLly9ehVPP/10o3X69OmDDRs2YPfu3di0aRNqamoQHR2NP/74o/UC1UPdZ63PeSguLoZGo2kT566iogIvvvgiHn/88SYfNqbvd7e1jR07Fp999hlSUlLw1ltv4bvvvkNcXBw0Gk2D9dvSOQKATz/9FJ06dWr2MolUzlNNTQ0WLFiAe+65ByEhIQBqf0tyubxeQtzUZ27I789cGurT7YqLi/Haa69h1qxZTbal7/fVXBrr0xNPPIFNmzbh8OHDWLJkCT7//HM8+eSTjbZj7vPEpye3ori4OO2f+/fvj8jISPj7+2Pbtm0t+j8lqfvkk08QFxcHHx+fRutERUUhKipKux0dHY3g4GCsX78er732WmuESf9VXV2N+Ph4CCGwdu3aJutK/bv72GOPaf8cGhqK/v37IzAwEKmpqRg9erQFIzONDRs2YPLkyc1OPJfKeUpISEBmZmarzmMyt+b6VFZWhnHjxqFfv3549dVXm2xLKt/Xxvp0a6IVGhoKb29vjB49Gr/99hsCAwNbLb46HFGxIFdXV/Tu3Rvnz59vcL+Xl1e9GfEqlQpeXl6tEZ5e8vLycPDgQcyYMUOv99na2mLAgAGNfgaWVvdZ63Me3N3dYW1tLelzV5ek5OXlITk5We9Htzf33bW0nj17wt3dvdH42sI5qvP9998jJydH798WYJnz9Nxzz+Gbb77B4cOH0b17d225l5cXqqqqcPXqVZ36TX3mhvz+zKGxPtW5du0axo4di06dOmHnzp2wtbXVq/3mvq/m0FyfbhUZGQkATf5bBZjvPDFRsaDy8nL89ttv8Pb2bnB/VFQUUlJSdMqSk5N1RiSkYuPGjejWrRvGjRun1/s0Gg3OnDnT6GdgaQEBAfDy8tI5D2VlZfjpp58aPQ9yuRwDBw7UeU9NTQ1SUlIkce7qkpTc3FwcPHgQXbt21buN5r67lvbHH3+gpKSk0fikfo5u9cknn2DgwIEICwvT+72teZ6EEHjuueewc+dOHDp0CAEBATr7Bw4cCFtbW53PPCcnBxcvXmz0Mzfk92dKzfWpLp4xY8ZALpfj66+/Nmi5hea+r6bUkj7d7tSpUwDQaHxmP09GT8elFvvHP/4hUlNThVKpFD/88IOIiYkR7u7u4tKlS0IIIZ566inx0ksvaev/8MMPwsbGRrzzzjsiOztbJCYmCltbW3HmzBlLdaFBGo1G+Pn5iRdffLHevtv7tGzZMrF//37x22+/iYyMDPHYY48Je3t7cfbs2dYMWce1a9fEyZMnxcmTJwUAsXLlSnHy5EntXTD/+te/hKurq9i9e7c4ffq0GD9+vAgICBA3btzQtjFq1Cjx4Ycfarf/7//+T9jZ2YmkpCSRlZUlZs2aJVxdXUVhYaFF+1NVVSUefPBB0b17d3Hq1Cnx119/aV+VlZWN9qe5764l+3Tt2jWxePFikZaWJpRKpTh48KC46667RK9evURFRUWjfbLkOWquT3VKS0uFo6OjWLt2bYNtSOk8Pfvss6Jz584iNTVV53t1/fp1bZ05c+YIPz8/cejQIXH8+HERFRUloqKidNrp06eP+Oqrr7TbLfn9WapPpaWlIjIyUoSGhorz58/r1Ll582aDfWrp99VSfTp//rxYvny5OH78uFAqlWL37t2iZ8+eYtiwYTrttOZ5YqLSih599FHh7e0t5HK5uOOOO8Sjjz4qzp8/r90/fPhwMXXqVJ33bNu2TfTu3VvI5XJx5513ij179rRy1M3bv3+/ACBycnLq7bu9TwsWLBB+fn5CLpcLT09Pcd9994kTJ060YrT1HT58WACo96qLu6amRixdulR4enoKOzs7MXr06Hp99ff3F4mJiTplH374obavERER4tixYxbvj1KpbHAfAHH48OFG+9Pcd9eSfbp+/boYM2aM8PDwELa2tsLf31/MnDmzXsIhpXPUXJ/qrF+/Xjg4OIirV6822IaUzlNj36uNGzdq69y4cUPMnTtXdOnSRTg6OoqHHnpI/PXXX/XaufU9Lfn9WapPjZ1DAEKpVDbYp5Z+Xy3Vp4sXL4phw4YJNzc3YWdnJ4KCgsQLL7wgSktL67XTWudJ9t8DEhEREUkO56gQERGRZDFRISIiIsliokJERESSxUSFiIiIJIuJChEREUkWExUiIiKSLCYqREREJFlMVIiozUlNTYVMJkNqaqpB73/11Vchk8lMG5QeRowY0egTeIlIFxMVImpSUlISZDKZ9mVjY4M77rgDTz/9NP78809LhydZBQUFePXVV7XPSSEiw9hYOgAiahuWL1+OgIAAVFRU4NixY0hKSsLRo0eRmZlp0IPY2ruCggIsW7YMPXr0QHh4uKXDIWqzmKgQUYvExcVh0KBBAIAZM2bA3d0db731Fr7++mvEx8dbODoiaq946YeIDDJ06FAAwG+//aZT/uuvv2LixIlwc3ODvb09Bg0ahK+//lqnzuXLl7F48WKEhobC2dkZLi4uiIuLwy+//FLvOH/88QcmTJgAJycndOvWDQsXLkRlZWWL4zx69CgGDx4Me3t7BAYGYv369Y3W3bRpEwYOHAgHBwe4ubnhscceQ35+vk6duvklGRkZiI6OhoODAwICArBu3TptndTUVAwePBgAMG3aNO1ls6SkJJ22srKyMHLkSDg6OuKOO+7A22+/3eJ+EXUUHFEhIoNcuHABANClSxdt2dmzZ3HPPffgjjvuwEsvvQQnJyds27YNEyZMwJdffomHHnoIAPD7779j165dmDRpEgICAqBSqbB+/XoMHz4cWVlZ8PHxAQDcuHEDo0ePxsWLF/H888/Dx8cHn3/+OQ4dOtSiGM+cOYMxY8bAw8MDr776Km7evInExER4enrWq/vGG29g6dKliI+Px4wZM1BUVIQPP/wQw4YNw8mTJ+Hq6qqte+XKFdx3332Ij4/H448/jm3btuHZZ5+FXC7H9OnTERwcjOXLl+OVV17BrFmztElddHS0Thtjx47Fww8/jPj4eOzYsQMvvvgiQkNDERcXp9e5IGrXTPIMZiJqtzZu3CgAiIMHD4qioiKRn58vduzYITw8PISdnZ3Iz8/X1h09erQIDQ0VFRUV2rKamhoRHR0tevXqpS2rqKgQGo1G5zhKpVLY2dmJ5cuXa8tWrVolAIht27Zpy9RqtQgKChIAxOHDh5uMfcKECcLe3l7k5eVpy7KysoS1tbW49a+/CxcuCGtra/HGG2/ovP/MmTPCxsZGp3z48OECgHj33Xe1ZZWVlSI8PFx069ZNVFVVCSGESE9PFwDExo0b68VV18Znn32m04aXl5d45JFHmuwTUUfDSz9E1CIxMTHw8PCAr68vJk6cCCcnJ3z99dfo3r07gNrLOYcOHUJ8fDyuXbuG4uJiFBcXo6SkBLGxscjNzdXeJWRnZwcrq9q/fjQaDUpKSuDs7Iw+ffrgxIkT2mPu3bsX3t7emDhxorbM0dERs2bNajZejUaD/fv3Y8KECfDz89OWBwcHIzY2VqfuV199hZqaGsTHx2vjLi4uhpeXF3r16oXDhw/r1LexscHs2bO123K5HLNnz8alS5eQkZHRos/T2dkZTz75pE4bERER+P3331v0fqKOgpd+iKhFFAoFevfujdLSUmzYsAFHjhyBnZ2ddv/58+chhMDSpUuxdOnSBtu4dOkS7rjjDtTU1OD999/HmjVroFQqodFotHW6du2q/XNeXh6CgoLqrXnSp0+fZuMtKirCjRs30KtXr3r7+vTpg71792q3c3NzIYRosC4A2Nra6mz7+PjAyclJp6x3794Aai+J3X333c3G171793r96tKlC06fPt3se4k6EiYqRNQiERER2rt+JkyYgCFDhuCJJ55ATk4OnJ2dUVNTAwBYvHhxvRGLOkFBQQCA//3f/8XSpUsxffp0vPbaa3Bzc4OVlRUWLFigbac11dTUQCaT4dtvv4W1tXW9/c7OziY/ZkPHAQAhhMmPRdSWMVEhIr1ZW1vjzTffxMiRI7F69Wq89NJL6NmzJ4Da0YeYmJgm379jxw6MHDkSn3zyiU751atX4e7urt329/dHZmYmhBA6ow85OTnNxujh4QEHBwfk5ubW23f7+wMDAyGEQEBAgHZkpCkFBQVQq9U6oyrnzp0DAPTo0QMALLryLVF7wjkqRGSQESNGICIiAqtWrUJFRQW6deuGESNGYP369fjrr7/q1S8qKtL+2draut7Iwfbt2+utdHvfffehoKAAO3bs0JZdv34d//73v5uNz9raGrGxsdi1axcuXryoLc/Ozsb+/ft16j788MOwtrbGsmXL6sUlhEBJSYlO2c2bN3Vuc66qqsL69evh4eGBgQMHAoA2ibl69WqzsRJR4ziiQkQGe+GFFzBp0iQkJSVhzpw5UCgUGDJkCEJDQzFz5kz07NkTKpUKaWlp+OOPP7TrpNx///1Yvnw5pk2bhujoaJw5cwabN2/WjsrUmTlzJlavXo0pU6YgIyMD3t7e+Pzzz+Ho6Nii+JYtW4Z9+/Zh6NChmDt3Lm7evIkPP/wQd955p85ckMDAQLz++utYsmQJLly4gAkTJqBTp05QKpXYuXMnZs2ahcWLF2vr+/j44K233sKFCxfQu3dvbN26FadOncK///1v7XyWwMBAuLq6Yt26dejUqROcnJwQGRmJgIAAYz92oo7FcjccEVFbUHd7cnp6er19Go1GBAYGisDAQHHz5k0hhBC//fabmDJlivDy8hK2trbijjvuEPfff7/YsWOH9n0VFRXiH//4h/D29hYODg7innvuEWlpaWL48OFi+PDhOsfIy8sTDz74oHB0dBTu7u5i/vz5Yt++fS26PVkIIb777jsxcOBAIZfLRc+ePcW6detEYmKiaOivvy+//FIMGTJEODk5CScnJ9G3b1+RkJAgcnJytHWGDx8u7rzzTnH8+HERFRUl7O3thb+/v1i9enW99nbv3i369esnbGxsdG5VrmvjdlOnThX+/v7N9omoI5EJwZlbREQtNWLECBQXFyMzM9PSoRB1CJyjQkRERJLFRIWIiIgki4kKERERSRbnqBAREZFkcUSFiIiIJIuJChEREUkWExUiIiKSLCYqREREJFlMVIiIiEiymKgQERGRZDFRISIiIsliokJERESSxUSFiIiIJOv/A58d2eOt4YVMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nreps = 2000\n",
    "qs = [20, 30, 40]\n",
    "read_number = np.arange(5, 30, 5)\n",
    "thresh = 2\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(6,3))\n",
    "colors = [\"red\", \"blue\", \"orange\", \"green\"]\n",
    "for j,q in enumerate(qs):\n",
    "    for m in tqdm(read_number):\n",
    "        llrs = np.zeros(nreps)\n",
    "        for i in range(nreps):\n",
    "            _, _, llr = geno_likelihoods_het(m=m, q=q, seed=(i+q+m+1))\n",
    "            llrs[i] = llr\n",
    "        if m == np.min(read_number):\n",
    "            ax.scatter(m, np.sum(llrs <= thresh)/llrs.size, color=colors[j], marker=\"+\", s=20, label=f'Q={q}')\n",
    "        else:\n",
    "            ax.scatter(m, np.sum(llrs <= thresh)/llrs.size, s=20, marker=\"+\", color=colors[j])\n",
    "\n",
    "ax.legend()\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel('Prop. $LLR_{het} < 2$', fontsize=12)\n",
    "ax.set_xlabel('Read depth', fontsize=12)\n",
    "debox(ax);\n",
    "plt.savefig('denovo_fpr_vs_depth_qual.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f12823f6-ef1f-4f0b-a9cc-9b414bde0daf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

ultrasel_zygote_calcs.ipynb