Week5/DeterminingGenotypes/DeterminingGenotypes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Determining genotypes\n",
    "\n",
    "This excersize aims at applying Bayesian statistics to calculate the probability of a specific genotype. In particular we will:\n",
    "1. Combine the use of Markov changes with the probability theorem of Bayes.\n",
    "2. Compare the computed results with simple simulations.\n",
    "\n",
    "The exercise is made to play around with a simulated sample of genes and DNA. The goal is to calculate the probability of an original state, based on data and prior probability and therefore a direct application of Bayes theorem.\n",
    "\n",
    "### References:\n",
    "- Lecture on Bayesian statistics\n",
    "- Lecture on Markov Chains\n",
    "\n",
    "### Authors: \n",
    "- Mathias L. Heltberg (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 29-12-2021 (latest update)\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always we will include a collection of libraries (Standard)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "from numpy.linalg import matrix_power                  # This is to raise matrices to powers\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also (big surprise...) need a random generator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.random # Random generator\n",
    "r.seed(42)    # Set a random seed (but a fixed one)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# And now for the actual problem: \n",
    "\n",
    "You are given a sample of a particular gene (for simplicity say eyecolor). This gene is found in two versions: A and B. Given that we all have DNA from our fathers and mothers, there are 3 possibilities of combinations for genotypes: AA, AB and BB.\n",
    "\n",
    "So say you measure everything optimally, and from the laboratory you extract genes from some cells. If you have 8 versions of the gene, coming from AB it might look like: \n",
    "- [A,B,A,B,A,A,A,B]\n",
    "\n",
    "However, you don't have the best equipment and there is a possibility that you measure a gene wrong. This is given by a matrix, XX.\n",
    "\n",
    "**Assume you get the sequence: [A,A,A,A,B]. We want to compute the probability for each of the three genotypes given sequences of data.**\n",
    "\n",
    "Start by writing on paper, thinking about how you would calculate this probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information 1 (Transformations matrix): \n",
    "\n",
    "The matrix of transformations called XX (meaning the probability that your machine measures A if it was really a B), takes the form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XX = np.zeros((2,2)); XX[0,0] = 0.85; XX[0,1] = 0.15; XX[1,0] = 0.05; XX[1,1] = 0.95;\n",
    "print(\"XX (i.e. the matrix of transformations error) = \")\n",
    "print(XX)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information 2 (Prior knowledge): \n",
    "\n",
    "You know that in the population, there are the probabilities for each combination (this is your prior!):\n",
    "- AA = 0.45\n",
    "- BB = 0.40\n",
    "- AB = 0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pri = np.zeros(3); Pri[0] = 0.45; Pri[1] = 0.15; Pri[2] = 0.40;\n",
    "print(\"Prior = \")\n",
    "print(Pri)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that we do measure a gene (say B), then the probability that it comes from BB is larger than it coming from AB (make sure you understand this). This we put in the matrix S:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = np.zeros((2,3)); \n",
    "S[0,0] = 2./3; S[0,1] = 1./3; S[1,1] = 1./3; S[1,2] = 2./3 ## Matrix to end up in a genotype\n",
    "print(\"S = \")\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now construct a full matrix of the seven states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.zeros((7,7)); M[0:2,2:4] = np.transpose(XX); M[2:4,4:] = S; M[6,6] = 1; M[5,5] = 1; M[4,4] = 1;\n",
    "print(\"M = \")\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what is the probability a measured base B comes from the genotype BB?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = np.zeros([1,2]); V1[0,1] = 1;\n",
    "XX2 = np.transpose(XX) \n",
    "P = np.matmul(V1,np.matmul(XX,S)); \n",
    "print(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think about why do we use the transposed matrix? \n",
    "\n",
    "So now the probability that the sequence [A,A,A,A,B] comes from the 3 genotypes must be the product of these probabilities (think about why this has to be the case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Seq = [0,0,0,0,1]\n",
    "LLout = np.ones(3)\n",
    "LLoutFull = np.ones(3)\n",
    "for igeno in range(3):    \n",
    "    V2 = np.zeros(3); V2[igeno]= 1;\n",
    "    V2Full = np.zeros(7); V2Full[4+igeno]= 1;\n",
    "    for iseq in range(len(Seq)):\n",
    "        V1Full = np.zeros([1,7]); V1Full[0,Seq[iseq]] = 1;\n",
    "        LLoutFull[igeno] *= np.matmul(np.matmul(V1Full,np.linalg.matrix_power(M,2)),V2Full)\n",
    "        V1 = np.zeros([1,2]); V1[0,Seq[iseq]] = 1;        \n",
    "        LLout[igeno] *= np.matmul(np.matmul(V1,np.matmul(XX2,S)),V2)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we could use the full matrix and square it, and we could use the combination of our XX matrix and S. These give the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The likelihood using the squared matrix is:\")\n",
    "print(LLoutFull)\n",
    "print(\"The likelihood using the two submatrices is:\")\n",
    "print(LLout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we should produce a probability based on this likelihood. We can get a probability by including the prior or not - but these results differ of course!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Probabilities with the included prior \")\n",
    "for i in range(3):\n",
    "    print(LLout[i]/np.dot(LLout,Pri)*Pri[i])\n",
    "print(\" \")\n",
    "print(\"Probabilities without the included prior \")\n",
    "for i in range(3):\n",
    "    print(LLout[i]/np.sum(LLout))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we have obtained the predicted probabilites based on our calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets make a simulation to test if this is true!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 0 ## counter \n",
    "nsave = 100 ## 100 datapoints to save - only for plotting!\n",
    "lseq = 5 ## length of the sequence\n",
    "Ntest = 100000 ## Number of points to generate\n",
    "PlAr = np.zeros((nsave,3))\n",
    "ntrue = np.zeros(Ntest) ## array where we save the actual genotype that gave rise to a sequence\n",
    "ns = np.zeros(5) ## Here we put the test-sequence for all runs\n",
    "click = 1\n",
    "\n",
    "while(c<Ntest):\n",
    "    A = np.random.random() ## This is a random number to determine genotype    \n",
    "    if (A<Pri[0]): ## For genotype AA\n",
    "        ntrue[c] = 0\n",
    "        for i in range(lseq):\n",
    "            B = np.random.random()\n",
    "            if (B<XX[0,0]):\n",
    "                ns[i] = 0\n",
    "            else:\n",
    "                ns[i] = 1\n",
    "    elif (A<Pri[0]+Pri[1]): ## For genotype AB\n",
    "        ntrue[c] = 1\n",
    "        for i in range(lseq):\n",
    "            if (np.random.random()<0.5):\n",
    "                B = np.random.random()\n",
    "                if (B<XX[0,0]):\n",
    "                    ns[i] = 0\n",
    "                else:\n",
    "                    ns[i] = 1\n",
    "            else: \n",
    "                B = np.random.random()\n",
    "                if (B<XX[1,0]):\n",
    "                    ns[i] = 0\n",
    "                else:\n",
    "                    ns[i] = 1\n",
    "    else: ## For genotype BB\n",
    "        ntrue[c] = 2\n",
    "        for i in range(lseq):\n",
    "            B = np.random.random()\n",
    "            if (B<XX[1,0]):\n",
    "                ns[i] = 0\n",
    "            else:\n",
    "                ns[i] = 1\n",
    "                \n",
    "    if (np.sum(ns) == 1): ## If our simulated sequence matches the desired sequence\n",
    "        c+= 1\n",
    "        \n",
    "        if (c > click*Ntest/nsave): \n",
    "            l1 = np.ones(c)\n",
    "            PlAr[click,0] = np.sum(l1[ntrue[0:c]==0])/len(ntrue[0:c])\n",
    "            PlAr[click,1] = np.sum(l1[ntrue[0:c]==1])/len(ntrue[0:c])\n",
    "            PlAr[click,2] = np.sum(l1[ntrue[0:c]==2])/len(ntrue[0:c])            \n",
    "            click+=1\n",
    "\n",
    "fig, ax = plt.subplots(3)\n",
    "ax[0].plot(np.linspace(1*Ntest/nsave,Ntest,nsave-1),PlAr[1:,0],'-g')\n",
    "ax[0].set_xlabel('# of tested sequence')\n",
    "ax[0].set_ylabel('Probability')\n",
    "ax[0].set_title('Probability of state AA')\n",
    "\n",
    "ax[1].plot(np.linspace(1*Ntest/nsave,Ntest,nsave-1),PlAr[1:,1],'-b')\n",
    "ax[1].set_xlabel('# of tested sequence')\n",
    "ax[1].set_ylabel('Probability')\n",
    "ax[1].set_title('Probability of state AB')\n",
    "\n",
    "ax[2].plot(np.linspace(1*Ntest/nsave,Ntest,nsave-1),PlAr[1:,2],'-r')\n",
    "ax[2].set_xlabel('# of tested sequence')\n",
    "ax[2].set_ylabel('Probability')\n",
    "ax[2].set_title('Probability of state BB')\n",
    "fig.tight_layout()\n",
    "fig.set_size_inches(18.5, 10.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = np.ones(len(ntrue))\n",
    "print(\"When we measure a sequence [A,A,A,A,B] the fractions of genotypes are:\")\n",
    "print(np.sum(l1[ntrue==0])/len(ntrue))\n",
    "print(np.sum(l1[ntrue==1])/len(ntrue))\n",
    "print(np.sum(l1[ntrue==2])/len(ntrue))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly this agrees well with the calculated result when we use the prior!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculations by updating bayesian prior\n",
    "Below here, this code can be used to calculate the probabilities by updating the prior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######### Try to make an updated prior ##########\n",
    "Seq = [0,0,0,0,1]\n",
    "LLout = np.ones(3)\n",
    "XX2 = np.transpose(XX) \n",
    "Pri0 = np.zeros(3); Pri0[0] = 0.45; Pri0[1] = 0.15; Pri0[2] = 0.40;\n",
    "for iseq in range(len(Seq)): \n",
    "    V1 = np.zeros([1,2]); V1[0,Seq[iseq]] = 1;        \n",
    "    LLout = np.matmul(V1,np.matmul(XX2,S))\n",
    "    Pout = 0 ####### You write here but use Pri0\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulations of specific sequences\n",
    "This piece of code can be used to generate sequences of length Nseq in order to calculate the probabilities for the questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######### Generate probabilities of arbitrary length ##########\n",
    "Nseq = 8\n",
    "ProbAA = np.zeros((Nseq,Nseq+1))\n",
    "ProbAB = np.zeros((Nseq,Nseq+1))\n",
    "ProbBB = np.zeros((Nseq,Nseq+1))\n",
    "MaxProb = np.zeros((Nseq,Nseq+1))\n",
    "for tseq in range(Nseq): \n",
    "    Seq = np.zeros(tseq, dtype=int)\n",
    "    for vseq in range(tseq+1):\n",
    "        if (vseq > 0):\n",
    "            Seq[vseq-1] = 1\n",
    "        LLout = np.ones(3)    \n",
    "        for iseq in range(tseq): \n",
    "            V1 = np.zeros([1,2]); V1[0,Seq[iseq]] = 1;   \n",
    "            L1 = np.matmul(V1,np.matmul(XX2,S))\n",
    "            LLout *= L1.flatten()\n",
    "        ProbAA[tseq,vseq] = LLout[0]/np.dot(LLout,Pri)*Pri[0]\n",
    "        ProbAB[tseq,vseq] = LLout[1]/np.dot(LLout,Pri)*Pri[1]\n",
    "        ProbBB[tseq,vseq] = LLout[2]/np.dot(LLout,Pri)*Pri[2]\n",
    "        MaxProb[tseq,vseq] = np.max(LLout/np.dot(LLout,Pri)*Pri)\n",
    "        \n",
    "fig, ax = plt.subplots(3)\n",
    "ax[0].set_title('Probability for AA')\n",
    "ax[0].imshow(ProbAA)\n",
    "ax[0].set_aspect('equal')\n",
    "ax[0].set_xlabel('Number of B in sequence')\n",
    "ax[0].set_ylabel('Sequence length')\n",
    "\n",
    "ax[1].set_title('Probability for AB')\n",
    "ax[1].imshow(ProbAB)\n",
    "ax[1].set_aspect('equal')\n",
    "ax[1].set_xlabel('Number of B in sequence')\n",
    "ax[1].set_ylabel('Sequence length')\n",
    "\n",
    "ax[2].set_title('Probability for BB')\n",
    "ax[2].imshow(ProbBB)\n",
    "ax[2].set_aspect('equal')\n",
    "ax[2].set_xlabel('Number of B in sequence')\n",
    "ax[2].set_ylabel('Sequence length')\n",
    "\n",
    "#fig.tight_layout()\n",
    "fig.set_size_inches(28.5, 20.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the original values - why is the P(BB) band more narrow than the P(AA)?\n",
    "\n",
    "Here we can plot some relevant probabilities based on the simulation above. For instance the maximal probability of AB and the minimal value of the most probable sequence of bases for each length of sequences (make sure you understand this)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MinProb = np.zeros(Nseq)\n",
    "MaxProbAB = np.zeros(Nseq)\n",
    "for i in range(Nseq):\n",
    "    MinProb[i] = np.min(MaxProb[i,0:i+1])\n",
    "    MaxProbAB[i] = np.max(ProbAB[i,0:i+1])\n",
    "xseq = np.linspace(1,Nseq,Nseq)\n",
    "plt.plot(xseq,MinProb,'--xr')\n",
    "plt.plot(xseq,MaxProbAB,'--+g')\n",
    "plt.title('Probabilities for specific values at different sequence lengths')\n",
    "plt.ylabel('Probability')\n",
    "plt.xlabel('Sequence')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions:\n",
    "\n",
    "This excersize should make you think of how to combine two different types of uncertainty with the concept of a bayesian prior.\n",
    "\n",
    "1) Try to change the sequence of length 5 so no genotype probability is larger than 0.66? Can you make a sequence of length 10 that does not have any probability larger than 0.66?\n",
    "\n",
    "2) For the sequence of length 5: What happens if you update the prior for each base you are reading? Does it give the same result?\n",
    "\n",
    "3) Assume for this question that the priors P(AA) = P(BB) = 0.49. How long a sequence do you need to make the call (p>0.95) of a the genotype AB if you have equal number of A's and B's?\n",
    "\n",
    "4) Assume again that the priors P(AA) = P(BB) = 0.49 but that the matrix XX now has the values \n",
    "- XX[0,0] = 0.99. \n",
    "- XX[1,0] = 0.01. \n",
    "- XX[0,1] = 0.01. \n",
    "- XX[1,1] = 0.99. \n",
    "\n",
    "For the sequence: [A,A,A,A,B] - does the prior change the result in an important way? \n",
    "\n",
    "5) Lets say you are in charge of deciding the strategy of the lab department that should determine genotypes. Use the original prior. \n",
    "- How much do you need to improve the measurement accuracy, if you want p>0.95 for all sequences of length 10?\n",
    "\n",
    "### Advanced question:\n",
    "\n",
    "6) Assume that you had a 10% error in one of the measurement estimates (i.e. the entries in the matrix XX). Which one would afffect the probabilites the most?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points\n",
    "\n",
    "In this excersize we want to compute probabilities in a setting where we have meaningful priors: the probability of observing specific genotypes.\n",
    "\n",
    "Here the observation of each genotype is intependent, and you do not have to use the Markov transfer matrix - but it may simplify the calculations. However, you should appreciate the use of prior probabilities and how they can be of extreme importance in the calculation of resulting probabilities.\n",
    "\n",
    "From this exercise, you should have learned:\n",
    "\n",
    "1) How to apply a Bayesian prior and calculate outcome probabilities. \n",
    "\n",
    "2) How to use the structure of the Markov Transfer matrix to obtain probabilities effectively when many outcomes are possible\n",
    "\n",
    "3) Understand how the change in prior probabilities affect the probability estimates \n",
    "\n",
    "4) Understand how increasing the length of a sequence, and thereby obtaining more data, determines the outcome probability compared to a change in prior probability."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}