coverage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating expected number of unique and repeated genotypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 7776  # possible number of genotypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Assuming all genotypes are equaly probable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_num_equal_prob(n_colonies):\n",
    "    \n",
    "    num_repetitions = 10000\n",
    "    duplicates = np.zeros(num_repetitions)\n",
    "    for i in range(num_repetitions):\n",
    "        sample = np.random.randint(N, size=n_colonies)    \n",
    "        duplicates[i] = n_colonies - len(set(sample))\n",
    "\n",
    "    mean_dupl = np.round(np.mean(duplicates))\n",
    "    mean_unique = n_colonies - mean_dupl\n",
    "    exp_dupl_percentage = np.round(100*mean_dupl/n_colonies)\n",
    "    print(f'Expected number of unique genotypes: {mean_unique}')\n",
    "    print(f'Expected number of duplicate genotypes: {mean_dupl} ({exp_dupl_percentage}%)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 5627.0\n",
      "Expected number of duplicate genotypes: 4373.0 (44.0%)\n"
     ]
    }
   ],
   "source": [
    "exp_num_equal_prob(n_colonies=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 466.0\n",
      "Expected number of duplicate genotypes: 14.0 (3.0%)\n"
     ]
    }
   ],
   "source": [
    "exp_num_equal_prob(n_colonies=480)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Close form solution (http://matt.might.net/articles/counting-hash-collisions/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_number_theory(n_colonies, N):\n",
    "    num_unique = np.round(N*(1-(1-1/N)**n_colonies))\n",
    "    num_dupl = n_colonies - num_unique\n",
    "    print(f'Expected number of unique genotypes: {num_unique}')\n",
    "    print(f'Expected number of duplicate genotypes: {num_dupl} ({np.round(100*num_dupl/n_colonies)}%)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 4997.0\n",
      "Expected number of duplicate genotypes: 3003.0 (38.0%)\n"
     ]
    }
   ],
   "source": [
    "exp_number_theory(n_colonies=8000, N=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 466.0\n",
      "Expected number of duplicate genotypes: 14.0 (3.0%)\n"
     ]
    }
   ],
   "source": [
    "exp_number_theory(n_colonies=480, N=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 250.0\n",
      "Expected number of duplicate genotypes: 4.0 (2.0%)\n"
     ]
    }
   ],
   "source": [
    "exp_number_theory(n_colonies=254, N=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assuming probabilities of promoters are approximately as in Fig 2C "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "promoters_prob = np.zeros((5,6))\n",
    "promoters_prob[0] = np.array([12.2, 18.1, 12.2, 15.7, 14.6, 27.2])\n",
    "promoters_prob[1] = np.array([11, 18.1, 11.4, 23.6, 15, 20.9])\n",
    "promoters_prob[2] = np.array([4.3, 1.2, 29.5, 15, 15.4, 34.6])\n",
    "promoters_prob[3] = np.array([18.9, 17.3, 21.3, 7.5, 18.1, 16.9])\n",
    "promoters_prob[4] = np.array([15.7, 20.1, 23.2, 18.5, 3.5, 18.9])\n",
    "\n",
    "for i in range(5): \n",
    "    promoters_prob[i] = promoters_prob[i]/np.sum(promoters_prob[i])\n",
    "    promoters_prob[i] = np.cumsum(promoters_prob[i])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.122     , 0.303     , 0.425     , 0.582     , 0.728     ,\n",
       "        1.        ],\n",
       "       [0.11      , 0.291     , 0.405     , 0.641     , 0.791     ,\n",
       "        1.        ],\n",
       "       [0.043     , 0.055     , 0.35      , 0.5       , 0.654     ,\n",
       "        1.        ],\n",
       "       [0.189     , 0.362     , 0.575     , 0.65      , 0.831     ,\n",
       "        1.        ],\n",
       "       [0.15715716, 0.35835836, 0.59059059, 0.77577578, 0.81081081,\n",
       "        1.        ]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "promoters_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_sample(size):\n",
    "    sample = []\n",
    "    for i in range(size):\n",
    "        r = np.random.rand(5, 1)\n",
    "        prom = []\n",
    "        for g in range(5):\n",
    "            prom.append(str(np.where(promoters_prob[g] > r[g])[0][0] + 1))\n",
    "        sample.append(''.join(prom))\n",
    "            \n",
    "    return sample\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_num_from_prob(n_colonies):\n",
    "    \n",
    "    num_repetitions = 10000\n",
    "    duplicates = np.zeros(num_repetitions)\n",
    "    for i in range(num_repetitions):\n",
    "        sample = generate_sample(size=n_colonies)    \n",
    "        duplicates[i] = n_colonies - len(set(sample))\n",
    "\n",
    "    mean_dupl = np.round(np.mean(duplicates))\n",
    "    mean_unique = n_colonies - mean_dupl\n",
    "    exp_dupl_percentage = np.round(100*mean_dupl/n_colonies)\n",
    "    print(f'Expected number of unique genotypes: {mean_unique}')\n",
    "    print(f'Expected number of duplicate genotypes: {mean_dupl} ({exp_dupl_percentage}%)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of unique genotypes: 3759.0\n",
      "Expected number of duplicate genotypes: 4441.0 (54.0%)\n"
     ]
    }
   ],
   "source": [
    "n_colonies = 10000 * 0.82\n",
    "exp_num_from_prob(n_colonies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_colonies = 8000\n",
    "exp_num_from_prob(n_colonies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_colonies = 480\n",
    "exp_num_from_prob(n_colonies)\n",
    "dupl_exp = 41\n",
    "print(f'Number of duplicates in the experiment: {dupl_exp} ({np.round(100*dupl_exp/n_colonies)}%)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_colonies = 254\n",
    "exp_num_from_prob(n_colonies)\n",
    "dupl_exp = 9\n",
    "print(f'Number of duplicates in the experiment: {dupl_exp} ({np.round(100*dupl_exp/n_colonies)}%)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ART_v3.6",
   "language": "python",
   "name": "art_v3.6"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}