Week10_Molecular_Evolution_Solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular evolution modeling\n",
    "*Matthew Wakefield*\n",
    "\n",
    "In this tutorial we are going to construct and assess a model of molecular evolution on a tree. To simplify the analysis we will use a specific simple tree shown below. \n",
    "\n",
    "\n",
    "```\n",
    "         Sequence1            Sequence2\n",
    "             \\                   /\n",
    "              \\                 /\n",
    "               \\               /\n",
    " length = t1    \\             / length = t2\n",
    "                 \\           /\n",
    "                  \\         /\n",
    "                   \\       /\n",
    "                   Ancestral\n",
    "```\n",
    "\n",
    "This tree is used as it has the fundamental features that would be iterated over when implementing a pruning algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence1 = 'GATCGATCGATCGATCGATCGATCGATCGATCGATCGATC'\n",
    "sequence2 = 'GATCGATCGATCGATCGATCGATCGATCGATCGATCGATC'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to construct a Q matrix for the modeling. We do this using numpy, the numerical library for python, and use the array data structure.  This is by default floats.\n",
    "As we will not be varying the mutation rate we dont need to keep remaking this matrix and it will be a constant in our calculations. Similarly because the JC69 assumes motif probabilities (base probabilities) are constant and equal we dont need to worry about varying them either. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu=1\n",
    "jc69Q = numpy.array([[-3.0*mu/4.0,     mu/4.0,     mu/4.0,     mu/4.0],\n",
    "                     [     mu/4.0,-3.0*mu/4.0,     mu/4.0,     mu/4.0],\n",
    "                     [     mu/4.0,     mu/4.0,-3.0*mu/4.0,     mu/4.0],\n",
    "                     [     mu/4.0,     mu/4.0,     mu/4.0,-3.0*mu/4.0],\n",
    "                    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.75,  0.25,  0.25,  0.25],\n",
       "       [ 0.25, -0.75,  0.25,  0.25],\n",
       "       [ 0.25,  0.25, -0.75,  0.25],\n",
       "       [ 0.25,  0.25,  0.25, -0.75]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jc69Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be multiplying the Q matrix by time (t) which is a scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.25,  0.75,  0.75,  0.75],\n",
       "       [ 0.75, -2.25,  0.75,  0.75],\n",
       "       [ 0.75,  0.75, -2.25,  0.75],\n",
       "       [ 0.75,  0.75,  0.75, -2.25]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jc69Q*3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When calculating the probability we will need to get a single column of the matrix as this will be the ancestral node probabilities. To make this easy we can create a dictionary corresponding to the matrix index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.75,  0.25,  0.25,  0.25])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jc69Q[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_index = {'A':0,'G':1,'C':2,'T':3,}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.75,  0.25,  0.25,  0.25])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jc69Q[base_index['A']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.75"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# you can do multi dimensional lookups as well\n",
    "jc69Q[base_index['A']][base_index['A']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the probability across our tree we need to calculate five probabilities, the probability at each of the three nodes and the transition probability at each of the two branches.\n",
    "First we will calculate the probabilities of the tip nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Exercise 1:**\n",
    "\n",
    "Calculate the log-probability of each of the two tip nodes.\n",
    "\n",
    "The tip node probability is the product of the motif frequencies or $\\pi$ for\n",
    "each motif (i.e. each base in our sequence).  Under the Jukes-Cantor model this is the same for all\n",
    "motifs A, C, G, and T, with pi = 1/4 for each.  When accumulating probabilities across sites we want to multiply probabilities. However it is better to \n",
    "log the probabilities to avoid underflow, and sum the\n",
    "logs. We can then just keep the log of the probability.\n",
    "\n",
    "*Hints: under this model the probabilities of seq1 and seq2 will be the same*\n",
    "\n",
    "You will need numpy.log()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_pr_seq(seq):\n",
    "    return numpy.log(0.25)*len(seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-55.4517744448\n",
      "-55.4517744448\n",
      "8.27180612553e-25\n",
      "8.27180612553e-25\n"
     ]
    }
   ],
   "source": [
    "print(log_pr_seq(sequence1))\n",
    "print(log_pr_seq(sequence2))\n",
    "print(numpy.exp(log_pr_seq(sequence1)))\n",
    "print(numpy.exp(log_pr_seq(sequence2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to calculate the transition probability from each of the nodes, i.e. the matrix $P$. This is given by \n",
    "\n",
    "$$P = e^{tQ}$$\n",
    "\n",
    "We need to calculate the product of tQ and the exponential of the resulting matrix. To calculate the exponential we will use the scipy.linalg.expm() function that uses the Padé approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import expm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_transition_matrix(t, Q):\n",
    "    return expm(t*Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that after a very small time has passed, the chance of a base substitution is small (what do you expect to get if $t=0$?):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.99253738,  0.00248754,  0.00248754,  0.00248754],\n",
       "       [ 0.00248754,  0.99253738,  0.00248754,  0.00248754],\n",
       "       [ 0.00248754,  0.00248754,  0.99253738,  0.00248754],\n",
       "       [ 0.00248754,  0.00248754,  0.00248754,  0.99253738]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_transition_matrix(0.01, jc69Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a larger time has passed, the chance of a base substitution increases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.35150146,  0.21616618,  0.21616618,  0.21616618],\n",
       "       [ 0.21616618,  0.35150146,  0.21616618,  0.21616618],\n",
       "       [ 0.21616618,  0.21616618,  0.35150146,  0.21616618],\n",
       "       [ 0.21616618,  0.21616618,  0.21616618,  0.35150146]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_transition_matrix(2, jc69Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** What happens as you increase time `t` to larger and larger numbers? What does this mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to apply the probability matrix across the sequence. To do this we will need a function that calculates the transition probability matrix, and then processes the probabilities for each site in the sequence.\n",
    "\n",
    "We will have two probabilites at the ancestral site, the probability of the transition to sequence1 and the probability of the transition to sequence2.\n",
    "As we are interested in the probability for both data points we will multiply these probabilities to get the combined probability of sequence1 and sequence2 at that site.\n",
    "\n",
    "If we know the motif at this site in the ancestral sequence, we can just multiply the probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0249735250559\n"
     ]
    }
   ],
   "source": [
    "# Leaf motifs (nucleotides)\n",
    "motif1 = 'A'\n",
    "motif2 = 'T'\n",
    "\n",
    "# Ancestral motif\n",
    "motif_ancestral = 'G'\n",
    "\n",
    "# Branch lengths\n",
    "t1 = 1\n",
    "t2 = 1\n",
    "\n",
    "# substitution rate matrix\n",
    "Q=jc69Q\n",
    "\n",
    "# Then we can calculate the probability of this data:\n",
    "\n",
    "P1 = make_transition_matrix(t1,Q)\n",
    "P2 = make_transition_matrix(t2,Q)\n",
    "\n",
    "prob_substitution1 = P1[base_index[motif1], base_index[motif_ancestral]]\n",
    "prob_substitution2 = P1[base_index[motif2], base_index[motif_ancestral]]\n",
    "probability = prob_substitution1 * prob_substitution2\n",
    "\n",
    "print(probability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't actually know the ancestral sequence. We need to generate the probability for each possible alternative ancestral motif (A, C, G or T), and sum these probabilities to get the total probability of this model.\n",
    "\n",
    "We can vectorise this efficiently with numpy's element-wise multiplication, by pulling the relevant rows out of `P1` and `P2` and multiplying them together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.52590958,  0.15803014,  0.15803014,  0.15803014])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1[base_index[motif1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.15803014,  0.15803014,  0.15803014,  0.52590958])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P2[base_index[motif2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.08310956,  0.02497353,  0.02497353,  0.08310956])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1[base_index[motif1]] * P2[base_index[motif2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.21616617919084685"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is the sum over all possible ancestral motifs\n",
    "sum(P1[base_index[motif1]] * P2[base_index[motif2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a function that puts that all together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_site_probability(motif1, motif2, t1, t2, Q=jc69Q):\n",
    "    base_index = {'A':0,'G':1,'C':2,'T':3,}\n",
    "    P1 = make_transition_matrix(t1,Q)\n",
    "    P2 = make_transition_matrix(t2,Q)\n",
    "    return sum(P1[base_index[motif1]]*P2[base_index[motif2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.35150146242745961"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_site_probability('A','A',1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.21616617919084685"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_site_probability('A','T',1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above gives the probability of the model for one motif (base) in the sequence.\n",
    "\n",
    "To get the likelihood of the entire sequence we need to iterate over each site in the sequence, take the log and accumulate the sum of the logs. This is equivalent to just multiplying the probabilities, but more numerically stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sequence_loglikelihood(seq1, seq2,\n",
    "                               t1, t2,\n",
    "                               Q=jc69Q):\n",
    "    likelihood = 0\n",
    "    assert len(seq1)==len(seq2)\n",
    "    for i in range(len(seq1)):\n",
    "        likelihood += numpy.log(\n",
    "            get_site_probability(seq1[i],seq2[i],\n",
    "                                 t1,t2,Q))\n",
    "    # Technically we should multiply by tip probabilities too:\n",
    "    #likelihood += log_pr_seq(seq1) + log_pr_seq(seq2)\n",
    "    return likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#This is a special case, identical sequences and zero branch length\n",
    "get_sequence_loglikelihood('AAAA','AAAA',0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.5324468055704896"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_sequence_loglikelihood('AAAA','AATA',0.1,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.480386086500874"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_sequence_loglikelihood('AAAA','ACTA',0.1,0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above code works as expected, however we end up doing lots of calls to the expensive expm function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Optional exercise:** Refactor the functions to reduce the number of calls to expm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:**\n",
    "\n",
    "Try some different sequences. Does it matter what the base changes are? If you have the same proportion of changes do you get the same results? Same proportion of conserved sites? What happens as the sequence gets longer? Does swapping t1 and t2 give the same likelihood result?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:**\n",
    "\n",
    "Plot the likelihood against t using a single value of t for both branch lengths. Do this for identical sequences and sequences with 25% difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1250dd3c8>"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAE9lJREFUeJzt3X+s3XV9x/HXi1IYdG5guCreFktq18S2pJfcIMZlWRa0BcR7GNmsg5nNxc4MEpEFY4Uoi2VqyNSZTU11/uHaDFjE2ihQIXMxEu/glpb+AKoXhbZXFq5h1Sk30B/v/XG/lbP23PM9537P93zP9/t9PpIT7vl+zz3fz0nL63z6/vz4OiIEAKi+M4puAACgPwh8AKgJAh8AaoLAB4CaIPABoCYIfACoCQIfAGqCwAeAmiDwAaAmziy6Ac0uuOCCWLp0adHNAIBS2blz588jYijtdQMV+EuXLtXExETRzQCAUrH9XCevo6QDADVB4ANATRD4AFATBD4A1ASBDwA1kfssHdvrJP2jpAWSvhoRn+71NbbtmtJdOw7oZ0dm9MbzztGta1eoMTLc68sAQKnlGvi2F0j6Z0nvkHRY0mO2t0fEk726xrZdU9p4317NHD0uSZo6MqON9+2VJEIfAJrkXdK5TNJkRPwkIl6RdLeksV5e4K4dB34T9ifNHD2uO7bv7+VlAKD08g78YUmHmp4fTo79hu0NtidsT0xPT3d9gZ8dmWl5/MjMUW3bNdX1+wFAVRU+aBsRmyNiNCJGh4ZSVwaf5o3nnTPnubt2HMjSNAColLwDf0rSkqbni5NjPXPr2hVznpur9w8AdZR34D8mabnti22fJWm9pO29vEBjZFjnLmz9MX73nIW9vBQAlFqus3Qi4pjtmyTt0Oy0zK9FRM9HU89euEAvHT1x2vFXjh1v8WoAqKfca/gRcX9E/F5ELIuIO/O4xpGXjrY8/tLREwzcAkCi8EHbXmDgFgDSVSLw2w3cTjFwCwCSKhL4jZFhneHW5+Y4DAC1U4nAl6QT0fp4SNTxAUAVCvzhNnV8tlkAgAoFfrs6PtssAECFAr8xMqzzz517oRWzdQDUXWUCX5I+cc3KOc8xWwdA3VUq8JmtAwBzq1TgS8zWAYC5VC7w283WoY4PoM4qF/isugWA1ioX+NTxAaC1ygW+RB0fAFqpZOCz6hYATlfJwGfVLQCcrpKBz6pbADhdJQNfYtUtAJyqsoHPbB0A+P8qG/gSs3UAoFmlA59VtwDwqkoHPqtuAeBVlQ586vgA8KpKB75EHR8ATqp84LPqFgBmVT7wWXULALMyBb7tP7G93/YJ26OnnNtoe9L2AdtrszVz/lh1CwCzsvbw90n6Y0nfbz5o+y2S1ktaKWmdpC/aXpDxWvPGqlsAyBj4EfFURLTqIo9JujsiXo6In0qalHRZlmtlwWwdAMivhj8s6VDT88PJscIwWwdA3aUGvu2Hbe9r8RjrRQNsb7A9YXtienq6F2/ZErN1ANRdauBHxBURsarF41ttfm1K0pKm54uTY63ef3NEjEbE6NDQUHet7wKzdQDUXV4lne2S1ts+2/bFkpZLejSna3WE2ToA6i7rtMxrbR+W9DZJ37G9Q5IiYr+keyU9KelBSTdGxPGsjc2K2ToA6izrLJ1vRsTiiDg7Il4fEWubzt0ZEcsiYkVEPJC9qdkxWwdAnVV+pe2pmK0DoK5qF/jskQ+grmoX+OyRD6Cuahf47er4EmUdANVVu8CX5q7jSyzCAlBdtQz8dnV8FmEBqKpaBv6ta1e0nYbJ4C2AKqpl4DdGhnX95RfNeZ7BWwBVVMvAl6RNjdUswgJQK7UNfIlFWADqpdaBzyIsAHVS68BnERaAOql14LOZGoA6qXXgS+3r+Ldv29vXtgBAnmof+O3q+FvHDzJ4C6Ayah/47er4IQZvAVRH7QM/7daHDN4CqIraB77U/taHDN4CqAoCX7O9/LmwCAtAVRD4CRZhAag6Aj/BIiwAVUfgJ7gTFoCqI/CbcCcsAFVG4DfhTlgAqozAb8KdsABUGYHfhDthAagyAv8U3AkLQFVlCnzbd9l+2vYe29+0fV7TuY22J20fsL02e1P7hzthAaiirD38hyStiohLJP1I0kZJsv0WSeslrZS0TtIXbS/IeK2+aTd4y2wdAGWVKfAj4rsRcSx5Oi5pcfLzmKS7I+LliPippElJl2W5Vj+1W4R1ZOZoH1sCAL3Tyxr++yU9kPw8LOlQ07nDybFSaLe3jkRZB0A5pQa+7Ydt72vxGGt6zW2Sjkna2m0DbG+wPWF7Ynp6uttfz027LZMp6wAoo9TAj4grImJVi8e3JMn2X0h6l6TrI+LkcOeUpCVNb7M4Odbq/TdHxGhEjA4NDWX6ML3UbstkFmEBKKOss3TWSfqIpHdHxEtNp7ZLWm/7bNsXS1ou6dEs1+q3tBuj0MsHUDZZa/j/JOk1kh6yvdv2lyUpIvZLulfSk5IelHRjRBzPeK2+S+vlA0CZZJ2l8+aIWBIRa5LHB5vO3RkRyyJiRUQ80O59BhWDtwCqhJW2KSjrAKgKAj8FZR0AVUHgp6CsA6AqCPwOUNYBUAUEfgeYkw+gCgj8DjAnH0AVEPgdYvAWQNkR+B1KG7y9fdvePrUEAOaHwO9Cu7LO1vGD1PIBDDQCvwvtyjohbnIOYLAR+F1IG7zlJucABhmB36V2vXyJhVgABheB36W0wVvKOgAGFYE/D+1uck5ZB8CgIvDnod1NziXKOgAGE4E/D2llHVbeAhhEBP48tSvrsL8OgEFE4M/TrWtXyG3O08sHMGgI/HlqjAzr+ssvmvM8++sAGDQEfgabGqvbnqesA2CQEPgZsW0ygLIg8DNi22QAZUHgZ8S2yQDKgsDvgXZlnS1smwxgQBD4PZC2oRq1fACDgMDvgbRtk6nlAxgEBH6PsG0ygEGXKfBtf9L2Htu7bX/X9huT47b9BduTyflLe9PcwdUYGdaisxbMeX7jfXv62BoAOF3WHv5dEXFJRKyR9G1JH0+OXylpefLYIOlLGa9TCndeO/dCrJmjJ+jlAyhUpsCPiF82PV2k2Vu7StKYpK/HrHFJ59m+MMu1yoBdNAEMssw1fNt32j4k6Xq92sMflnSo6WWHk2Otfn+D7QnbE9PT01mbU7i0wVt6+QCKkhr4th+2va/FY0ySIuK2iFgiaaukm7ptQERsjojRiBgdGhrq/hMMmLTBW2r5AIqSGvgRcUVErGrx+NYpL90q6brk5ylJS5rOLU6OVV5jZFg3tNlFk1o+gKJknaWzvOnpmKSnk5+3S3pfMlvnckm/iIjns1yrTNJ20eRG5wCKcGbG3/+07RWSTkh6TtIHk+P3S7pK0qSklyT9ZcbrlM755y7U/7zUesEVNzoHUISss3SuS8o7l0TENRExlRyPiLgxIpZFxOqImOhNc8sjrZbPpmoA+o2VtjlJm6K5lU3VAPQZgZ+jdjc6D1HLB9BfBH6O0m50Ti0fQD8R+DlKu9G5RC0fQP8Q+DlLm6LJDVIA9AuB3wftavkSq28B9AeB3wdptXxW3wLoBwK/Dzqp5TNjB0DeCPw+2dRY3fYGKczYAZA3Ar+P2t0gRWLGDoB8Efh9lLb6lhk7APJE4PcZM3YAFIXA7zNm7AAoCoHfZ53M2KGXDyAPBH4B0mbszBw9wQAugJ4j8AuSNmOHAVwAvUbgF6QxMty2ly9Jd2zf36fWAKgDAr9Aab38IzOtb5EIAPNB4BeoMTKsG1IGcCnrAOgVAr9gaQO4zNgB0CsE/gBoV9phxg6AXiHwBwBbLgDoBwJ/QJx/7sK25yntAMiKwB8Qn7hmZdvzbLkAICsCf0B0MmPnw/fsJvQBzBuBP0DSZuyEpFvuJfQBzE9PAt/239oO2xckz237C7Ynbe+xfWkvrlMHaYuxTgT1fADzkznwbS+R9E5JB5sOXylpefLYIOlLWa9TF51suUA9H8B89KKH/zlJH9FsxeGkMUlfj1njks6zfWEPrlULab18iV4+gO5lCnzbY5KmIuKJU04NSzrU9Pxwcgwd6GQAlwVZALp1ZtoLbD8s6Q0tTt0m6WOaLefMm+0Nmi376KKL2odcnWxqzPbyt4wfnPM1W8YPavRNr01duAUAUgc9/Ii4IiJWnfqQ9BNJF0t6wvazkhZLetz2GyRNSVrS9DaLk2Ot3n9zRIxGxOjQ0FDWz1MpabN2JEo7ADo375JOROyNiNdFxNKIWKrZss2lEfHfkrZLel8yW+dySb+IiOd70+R6Savnzxw9oeu/8sM+tQZAmeU1D/9+zf4LYFLSVyT9TU7XqbxO6vmPPPMioQ8gVc8CP+np/zz5OSLixohYFhGrI2KiV9epo05KO4888yJTNQG0xUrbkmCqJoCsCPySYKomgKwI/BLZ1Fitty97bdvXbBk/SOgDaInAL5mtH3ibzj6z/R8bN0wB0AqBX0Kfue6S1NdQzwdwKgK/hDqt5zNVE0AzAr+kOp2qSegDOInAL7FOpmo+8syLDOICkETgl1onpR2JmTsAZhH4JbepsZrQB9ARAr8COpmfLzFdE6g7Ar8iOpmfL0kfvoeboAN1ReBXyGeuu0RnuP1rQtIt9xL6QB0R+BXSGBnWZ/90jVIyXyeCnj5QRwR+xTRGhvW596xJfV1Iuvme3QzkAjVC4FdQp9M1JQZygToh8Cuq0+maEuUdoC4I/ArrNPQp7wD1QOBXXDc9fRZnAdVG4NdAt6HPhmtANRH4NdFN6D/yzIt6x2f/M98GAeg7Ar9Gugn9H7/wa624/QEGc4EKIfBrppvQf/nYCQZzgQoh8Guom9CXZuv6Kz/+IL19oOQI/Jra1Fitz79njRZ2+Dfg168cp7cPlByBX2ONkWH9+O+v1vLXLer4d5i6CZQXgQ89dMsfdrSf/klbxg8yoAuUUKbAt32H7Snbu5PHVU3nNtqetH3A9trsTUWetn7gbV3V9U8O6BL8QHn0oof/uYhYkzzulyTbb5G0XtJKSeskfdH2gh5cCznqdjBXejX4GdQFBl9eJZ0xSXdHxMsR8VNJk5Iuy+la6KGTg7nndDqamzg5qMsqXWBw9SLwb7K9x/bXbJ+fHBuWdKjpNYeTY6exvcH2hO2J6enpHjQHWTVGhvXUJ6/surcvza7SXfrR7zCwCwyg1MC3/bDtfS0eY5K+JGmZpDWSnpf0D902ICI2R8RoRIwODQ11/QGQn26nbjbbMn5Qb/7Y/ZR5gAHiiOjNG9lLJX07IlbZ3ihJEfGp5NwOSXdERNt/74+OjsbExERP2oPeun3bXm0ZPzjv37/h8ou0qbG6hy0CcJLtnRExmva6rLN0Lmx6eq2kfcnP2yWtt3227YslLZf0aJZroVibGqv17Kev7mr6ZrMt4wcp9QAFy9TDt/2vmi3nhKRnJf11RDyfnLtN0vslHZN0c0Q8kPZ+9PDLYduuKW28b49mjp6Y93ucf+5CfeKalWqMtBzaAdCFTnv4PSvp9AKBXy7bdk3p1n/frQy5L0ladNYC3XntasIfmCcCH32zbdeUbrlntzLmviTCH5gPAh99l3Vg91SEP9AZAh+F6XXwS5IlXc9MH6AlAh+F62Wp51RnWPqzt/IFAEgEPgZIHj3+VigBoa4IfAycbbumdMf2/Toyc7Rv1+RLAHVA4GOgFRH+p2ItAKqCwEdpDEL4n4pBYpQJgY9SGsTwT0PZCEUj8FEJvdjGoQz40kAWBD4q6fZte7V1/KAG529t9THWMfgIfNRGv6Z9AnnKsraEwEdtlXEcADhpPveOIPCBFuoyJoDyWmDrmU9d1dXvdBr4Z867VUAJNUaG29ai+UJA0Y7n2Akn8IEmaV8IJ1E2Ql4W2Lm9N4EPzEOnXwzt8KWBVt771iW5vTeBDxSkF18aeeHLqP/6sQMsgQ/gNIP8ZYT5O6PoBgAA+oPAB4CaIPABoCYIfACoCQIfAGpioLZWsD0t6bkMb3GBpJ/3qDllwWeuvrp9XonP3K03RcRQ2osGKvCzsj3RyX4SVcJnrr66fV6Jz5wXSjoAUBMEPgDURNUCf3PRDSgAn7n66vZ5JT5zLipVwwcAzK1qPXwAwBwqEfi219k+YHvS9keLbk8/2P6a7Rds7yu6Lf1ge4nt79l+0vZ+2x8quk15s/1bth+1/UTymf+u6Db1g+0FtnfZ/nbRbekH28/a3mt7t+1cb/lX+pKO7QWSfiTpHZIOS3pM0nsj4slCG5Yz238g6VeSvh4Rq4puT95sXyjpwoh43PZrJO2U1Kjyn7NtS1oUEb+yvVDSDyR9KCLGC25armzfImlU0u9ExLuKbk/ebD8raTQicl93UIUe/mWSJiPiJxHxiqS7JY0V3KbcRcT3Jb1YdDv6JSKej4jHk5//V9JTkiq9f2/M+lXydGHyKHcPLYXtxZKulvTVottSRVUI/GFJh5qeH1bFg6DubC+VNCLpv4ptSf6S8sZuSS9Ieigiqv6ZPy/pI5LqdFPhkPRd2zttb8jzQlUIfNSI7d+W9A1JN0fEL4tuT94i4nhErJG0WNJltitbvrP9LkkvRMTOotvSZ78fEZdKulLSjUm5NhdVCPwpSc03gVycHEPFJHXsb0jaGhH3Fd2efoqII5K+J2ld0W3J0dslvTupad8t6Y9sbym2SfmLiKnkvy9I+qZmy9S5qELgPyZpue2LbZ8lab2k7QW3CT2WDGD+i6SnIuKzRbenH2wP2T4v+fkczU5MeLrYVuUnIjZGxOKIWKrZ/4//IyJuKLhZubK9KJmEINuLJL1TUm4z70of+BFxTNJNknZodiDv3ojYX2yr8mf73yT9UNIK24dt/1XRbcrZ2yX9uWZ7fbuTx1VFNypnF0r6nu09mu3YPBQRtZiqWCOvl/QD209IelTSdyLiwbwuVvppmQCAzpS+hw8A6AyBDwA1QeADQE0Q+ABQEwQ+ANQEgQ8ANUHgA0BNEPgAUBP/B7d3O22xnLXVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib\n",
    "\n",
    "lnLs = []\n",
    "ts = []\n",
    "\n",
    "for t in range(0,500):\n",
    "    t= t/100.0\n",
    "    ts.append(t)\n",
    "    lnLs.append(get_sequence_loglikelihood(sequence1,\n",
    "                                           sequence2,\n",
    "                                           t,t))\n",
    "\n",
    "matplotlib.pyplot.scatter(x=ts,y=lnLs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:9: RuntimeWarning: divide by zero encountered in log\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x12501af28>"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFsdJREFUeJzt3X+Q3HV9x/HX23DQI7YTHM4Cl8SjDBOncJrDHUwmnQ46SMQfzQJ2igN2Ou2QdqozIjYO0YwghEKbVh3r1Blo/cMm1eCQnihggAGnSg164S65RIz8KAQWLIEQKHjVM3n3j9uDvXP3+9293e+vz/f5mNnh9rvf++5nB3jt596fH19zdwEAwveGrBsAAEgHgQ8AJUHgA0BJEPgAUBIEPgCUBIEPACVB4ANASRD4AFASBD4AlMRxWTeg0cknn+xDQ0NZNwMACmX37t3Pu/tA3Hm5CvyhoSGNjY1l3QwAKBQze7Kd8yjpAEBJEPgAUBIEPgCUBIEPACVB4ANASeRqlk5aNo1O6usPPqWj7lpkpg+/c5k2V4ezbhYAJKpUgT86XtOGb05o+tjrx466a+uug3rw8Rd0z1XnZdY2AEhaaUo6m0YndeX2uWHf6JHnXtWKTXdpdLyWbsMAICWlCPzR8Zq27joYe94vf31MV906QegDCFIpAn/LzgNtn3vMpY079ibYGgDIRikCv3ZkqqPzp6aPadPoZEKtAYBsBB/4Cy3PbN11kNIOgKAEH/if+/b+Bf8upR0AIQk68EfHa3rxF9MtX++L+fRT08fo5QMIRtCBHzVYa5K2/PFKXb5qeeQ16OUDCEXQgf9MxGDtZauWqzoyqM3VYS0+flHL8+jlAwhF0IG/5MS+psdP7HvDnK0UbrgoeluFa29f+DgAAORF0IHv3vz48cfN7dFXRwYje/lHpqbp5QMovK4C38yuN7O9ZjZhZneb2Wn14+eZ2Uv14xNm9tneNLczR6aaD9i+1OQ4vXwAoeu2h7/F3d/m7islfUdSY7B/391X1h/Xdfk+HRsdr8lavHbakv7fOFYdGYwcwG315QEARdFV4Lv7yw1PF0tqUURJ35adB5o2xiRtWLui6e/EbZFMWQdAkXVdwzezG8zsKUmXaW4Pf7WZ7TGzu8zsrG7fp1OtZui4ZnrzrZzUYqBXoqwDoNhiA9/M7jWzfU0e6yTJ3T/j7sskbZP0sfqvPSTpLe7+dkn/JGk04vrrzWzMzMYOHTrU/Seqa1a2kaTBFsdnXfPB1t9NlHUAFFls4Lv7+e5+dpPHt+aduk3SJfXfedndX6n/fKekPjM7ucX1b3b3irtXBgYGuvw4r3vXW5tfq9XxWVG9f0lsqgagsLqdpXNmw9N1kn5aP36KmVn953Pr7/NCN+/Vqft/2vyvhVbHG0WVdbaxqRqAguq2hn9TvbyzV9IFkj5eP/4hSfvMbI+kL0m61L3VrPhktKrhR62+nRVV1nF1tr8+AORFV/e0dfdLWhz/sqQvd3Ptbp22pL/pPvitavuNqiOD+ty397fceK3T/fUBIA+CXWn7rrcO/MY8/P6+RS2nZM4X1cuXmKIJoHiCDPzR8Zpu212bMw/fJF3yjsHYQdlZcedR1gFQNEEG/padBzQ1fXTOMVd7A7aNoqZwUtYBUDRBBn43A7aNoso/rbZtAIC8CjLwWw3MtjNg2yiqrONiTj6AYgky8DesXaH+vrnbHXcyYNsoqqzDnHwARRJk4EvSCce9/tFOOrFPN1483PaAbaOoLwnm5AMokuACf3S8po07Jufse/N/08cWfL3qyGDkylsGbwEURXCB32yGztT00a564lFz8hm8BVAUwQV+r2boNIobvKWOD6AIggv8Xs3QmS9q8JZ98gEUQXCB38sZOvOv2wr75AMoguACvzoyqBsvHtbgkn6ZZnrmC52hM/+6USjrAMi74AI/Sdz+EECRBRf4s9Mya0em5JqZNrlxx2RPeuBxtz+klw8gz4IL/CSmZc6Km5PPIiwAeRZc4CcxLbNRVC+fRVgA8iy4wE9qWuas6sig3tBitRWLsADkWXCBn9S0zEbHWtydlx00AeRZcIGf1LTMRuygCaCIzL1FdzUDlUrFx8bGsm5GrNHxmq7cPtHy9cEl/Xrg6nen2CIAZWZmu929EndecD38NLCDJoAiCi7wR8drWnPTfTr96ju05qb7EiuvsIMmgKIJKvCTXHQ1HztoAiiaoAI/yUVXzUQN3rIIC0DeBBX4SS+6mi9qqid1fAB5E1TgJ73oaj4WYQEokqACP41FV/OxCAtAUQQV+GksupqPRVgAioKFV11iERaArLHwKiUswgJQFAR+D7AIC0AREPg9wCIsAEUQXOCntbXCfCzCApB3QQV+mlsrzMciLAB5F1Tgp721QqOoRVgSZR0A2etZ4JvZJ83Mzezk+nMzsy+Z2aNmttfMzunVe7WS9tYK87VahCVJ196+P5U2AEArPQl8M1sm6QJJBxsOXyjpzPpjvaSv9OK9oqS9tcJ8UXX8I1PT9PIBZKpXPfwvSPqUZialzFon6Ws+Y5ekJWZ2ao/er6kstlaY//5R0zAZvAWQpa4D38zWSaq5+555Lw1Keqrh+dP1Y4nJYmuF+e9/2arlLV9n8BZAlo5r5yQzu1fSKU1e+oykT2umnLMgZrZeMyUfLV/eOizbVR0ZTC3gm9lcHda/P3iwaT2fRVgAstRWD9/dz3f3s+c/JD0u6XRJe8zsCUlLJT1kZqdIqkla1nCZpfVj8699s7tX3L0yMDDQ7efJhagdNKnjA8hKVyUdd5909ze7+5C7D2mmbHOOu/9c0u2S/rQ+W2eVpJfc/dnum5x/LMICkEdJzsO/UzN/ATwq6RZJf53ge+UKi7AA5FFPA7/e03++/rO7+0fd/Qx3H3b3Yu173AUWYQHIo6BW2uYJi7AA5A2BnxAWYQHIGwI/ISzCApA3BH5CWIQFIG+CCvys9sJvZXN1uOXgLYuwAKQtmMDPci/8KCzCApAXwQR+lnvhR4kavGW2DoA0BRP4We+F30rUIqwjU9MptgRA2QUT+Fnvhd9K3EZulHUApCWYwM96L/woJ53Y1/I1yjoA0hJM4Ge9F36Uaz54VsvXWIQFIC3mHrEHQMoqlYqPjYW55c7IdXfrxV80r9kv6e/TxDULvqUAgJIzs93uXok7L5geft7F9fIBIGkEfkoYvAWQNQI/RQzeAsgSgZ8iBm8BZInAT1F1ZJBePoDMEPgpY/AWQFYI/JQxeAsgKwR+BqLKOllv9gYgXAR+BqLKOtwYBUBSCPwMVEcGW94YRaKsAyAZBH5GWt0YRWK2DoBkEPgZiboxCnPyASSBwM/IhrUrIu9rSy8fQK8dl3UDemV0vKYtOw/omSNTOm1JvzasXZGLrZFbqY4MauzJw9q662DT15mTD6DXgujh5/UG5nE2V4cjX897+wEUSxCBn9cbmLeDOfkA0hJE4Of1BubtYE4+gLQEEfh5vYF5O5iTDyAtQQR+nm9g3g7m5ANIQxCBn+cbmLeDOfkA0sBNzHNgdLymT2yfUKt/E9zkHEAUbmJeINWRQV22annL15mTD6AXCPycYE4+gKQR+DnC7Q8BJKkngW9mnzQzN7OT68/PM7OXzGyi/vhsL94ndNz+EECSug58M1sm6QJJ8zeF+b67r6w/ruv2fcogblbRptHJlFoCIES96OF/QdKnpJaTTNCBqLLO1l0HqeUDWLCuAt/M1kmqufueJi+vNrM9ZnaXmbWuVWCOqLKORC0fwMLFBr6Z3Wtm+5o81kn6tKRm9fmHJL3F3d8u6Z8kjUZcf72ZjZnZ2KFDhxb6OYJRHRmM7OVTywewULGB7+7nu/vZ8x+SHpd0uqQ9ZvaEpKWSHjKzU9z9ZXd/pf77d0rqmx3QbXL9m9294u6VgYGBnn2wIovr5VPLB7AQCy7puPuku7/Z3YfcfUjS05LOcfefm9kpZmaSZGbn1t/nhZ60uASqI4NafPyilq9vo5YPYAGSmof/IUn7zGyPpC9JutTztIdDAdxwUeuFWC72ygfQuZ4Ffr2n/3z95y+7+1nu/nZ3X+Xu/9Wr9ymLuFo+e+UD6BQrbXMsrpZPWQdAJwj8HItbiMUUTQCdIPBzLm6vfABoF4Gfc3F37WKKJoB2Efg5FzdFk+0WALSLwC+AqCmaErV8AO0h8AuA7RYA9AKBXxBM0QTQLQK/IOJq+Rt37E2xNQCKiMAvkKha/tT0MXr5ACIR+AXCQiwA3Qgm8EfHa1pz0306/eo7tOam+4Lt7cYN3ob6uQF0L4jAHx2vaeOOSdWOTMk1s7HYxh2TQYZf3OAttXwArQQR+Ft2HtDU9NE5x6amjwa5hXB1ZFCXr1re8nVq+QBaCSLwn2mxVXCr40W3uRq9EItePoBmggj801psMNbqeAiiavn08gE0E0Tgb1i7Qv19c+eo9/ctit14rMjiavnM2AEwXxCBXx0Z1I0XD2twSb9MM1sK33jxcOw0xiKLW4jFjB0A81mebjVbqVR8bGws62YUxuh4TVdun2j5en/fG/Tw9Rem2CIAWTCz3e5eiTsviB5+WTFjB0AnCPyCY8YOgHYR+AGIm7HDXbEASAR+EOJm7HBXLAASgR+EuBk7EqUdAAR+MOJug0hpBwCBH4i4GTsSpR2g7Aj8gGyuDlPaAdASgR+Ydko79PKBciLwA9NOaYdePlBOBH6A4ko7DOAC5UTgByqutMMALlA+BH6gmJsPYD4CP2AM4AJoROAHjAFcAI0I/MAxgAtgFoFfAgzgApC6DHwzu9bMamY2UX+8r+G1jWb2qJkdMLO13TcVC8UALgCpNz38L7j7yvrjTkkys9+XdKmksyS9V9I/m1l04iBRDOACSKqks07SN9z9l+7+35IelXRuQu+FNjCAC6AXgf8xM9trZl81s5PqxwYlPdVwztP1Y8hQOwO4l93ywxRbBCBNsYFvZvea2b4mj3WSviLpDEkrJT0r6R87bYCZrTezMTMbO3ToUMcfAJ2JK+088NhhQh8IVGzgu/v57n52k8e33P1/3P2oux+TdIteL9vUJC1ruMzS+rFm17/Z3SvuXhkYGOj28yBGOwO4Dzx2mHo+EKBuZ+mc2vD0Ikn76j/fLulSMzvBzE6XdKakH3XzXuiduF6+RD0fCFG3Nfy/N7NJM9sr6V2SPiFJ7r5f0q2SfiLpu5I+6u5Hu3wv9Eg7A7jU84HwmLtn3YbXVCoVHxsby7oZpXHZLT/UA48djjxnzRlv0rYrVqfUIgALYWa73b0Sdx4rbUts2xWrdcJx0f8JPPDYYbZeAAJB4Jfc313ytthztu46SOgDASDwS66der7EfjtACAh8aHN1WGvOeFPseZ/YPkHoAwVG4EPSTD0/LvRd0lW3EvpAURH4eM22K1bHLso65tKGb06k1CIAvUTgY452FmVNHxNz9IECIvAxR7uDuOy5AxQPgY/fsLk63Hbov+fz30u+QQB6gsBHU+2G/iPPvUpPHygIAh8ttTtdk/IOUAwEPiK1M11TorwDFAGBj1jthv4jz72qFZvuYp4+kFMEPtqy7YrVOvPNi2PP++Wvj+nK7RPsvQPkEIGPtt1z1Xlt9fSlmb13qOsD+ULgoyPtlnck6vpA3hD46Fi75R2Juj6QJwQ+FqST8g51fSAfCHws2LYrVre1OGsWN1IBskXgoyubq8P64p+sVF+b/yVt3XWQEg+QEQIfXauODOqRv31/23V9SjxANgh89EwndX2J3j6QNgIfPdVpXX+2t3/WZ79L8AMJI/DRc+3utNno1V8dJfiBhAUR+KPjNa256T6dfvUdWnPTfQRGDnQ6mDtrNvhZpQv0XuEDf3S8po07JlU7MiWXVDsypY07Jgn9HJgdzO2krj/rgccOa+jqOxjYBXqo8IG/ZecBTU0fnXNsavqotuw8kFGLMN+2K1YvqLcvzQzsDl19h0auu5svcaBLhQ/8Z45MdXQc2eimty9JL/5imlIP0KXCB/5pS/o7Oo5szfb2+xfS3dfrpR4Gd4HOFT7wN6xdof6+RXOO9fct0oa1KzJqEeJURwb18PUXdhX8s4O7hD/QPnP3rNvwmkql4mNjYx3/3uh4TVt2HtAzR6Z02pJ+bVi7QtWRwQRaiCSMjte04ZsTmj7W/bUuX7Vcm6vD3V8IKBAz2+3uldjzQgh8hGF0vKartk+oB7kvSVp8/CLdcNEwX/4IHoGPwto0Oqmtuw729JqEP0JG4KPwelnqmc8kXUb5B4Eg8BGMJIO/0Ukn9umaD57FXwEoHAIfwRkdr+na2/fryNR0au9JKQhFkErgm9m1kq6QdKh+6NPufqeZDUl6WNLsctdd7v5Xcdcj8NGJJGr9neDLAHmRZuC/4u7/MO/4kKTvuPvZnVyPwMdCzOyntFdTSdd8OkSJCGlpN/CPS6MxQJKqI4NzQjXrnv+s2e0grtw+EXsuXw5IQy96+H8m6WVJY5I+6e4v1nv4+yX9rP7aJnf/ftz16OGj1zaNTmrbroPKz0hVcvjSKK+elXTM7F5JpzR56TOSdkl6XpJLul7Sqe7+52Z2gqQ3uvsLZvYOSaOSznL3l5tcf72k9ZK0fPnydzz55JNxbQYWLK/lnzJhSmzvpT5LJ6pub2bfk/Q37h7ZfaeHj7TxBYC86OaLMJUavpmd6u7P1p9eJGlf/fiApMPuftTMfk/SmZIe7+a9gCTMr//PKlMpCPng0mtjT0n99dNtDf/fJK3UTFufkPSX7v6smV0i6TpJ05KOSbrG3b8ddz16+Mi7LNYCoFwWmemxG9/X0e+k0sN394+0OH6bpNu6uTaQR63+IpjFFwK6dTTBxbBMywR6KO4LoRHjB2hmkVli1ybwgYx08uXQCn9RhOfD71yW2LUJfKDAevGlkRb+oomWxnRVAh9AKor05RSqwt/TFgDQHgIfAEqCwAeAkiDwAaAkCHwAKIlc3eLQzA5JWuh2mSdrZufOMuEzlwOfuRy6+cxvcfeBuJNyFfjdMLOxdvaSCAmfuRz4zOWQxmempAMAJUHgA0BJhBT4N2fdgAzwmcuBz1wOiX/mYGr4AIBoIfXwAQARggh8M3uvmR0ws0fN7Oqs25M0M/uqmT1nZvuybktazGyZmd1vZj8xs/1m9vGs25Q0M/stM/uRme2pf+bPZd2mNJjZIjMbN7PvZN2WtJjZE2Y2aWYTZpbYbf8KX9Ixs0WSfibpPZKelvRjSR92959k2rAEmdkfSnpF0tea3TQ+RGZ2qqRT3f0hM/ttSbslVQP/92ySFrv7K2bWJ+kHkj7u7rsyblqizOwqSRVJv+PuH8i6PWkwsyckVdw90bUHIfTwz5X0qLs/7u6/kvQNSesyblOi3P0/JR3Ouh1pcvdn3f2h+s//K+lhSUHvteszXqk/7as/it1Di2FmSyW9X9K/ZN2WEIUQ+IOSnmp4/rQCD4KyM7MhSSOSHsy2JcmrlzcmJD0n6R53D/0zf1HSpySV7S4pLuluM9ttZuuTepMQAh8lYmZvlHSbpCvd/eWs25M0dz/q7islLZV0rpkFW8Izsw9Ies7dd2fdlgz8gbufI+lCSR+tl217LoTAr0lqvAnk0voxBKZex75N0jZ335F1e9Lk7kck3S/pvVm3JUFrJP1RvZ79DUnvNrOt2TYpHe5eq//zOUn/oZlSdc+FEPg/lnSmmZ1uZsdLulTS7Rm3CT1WH8D8V0kPu/vns25PGsxswMyW1H/u18zEhJ9m26rkuPtGd1/q7kOa+f/4Pne/PONmJc7MFtcnIsjMFku6QFIiM/AKH/ju/mtJH5O0UzMDebe6+/5sW5UsM/u6pB9KWmFmT5vZX2TdphSskfQRzfT6JuqP92XdqISdKul+M9urmY7NPe5emqmKJfK7kn5gZnsk/UjSHe7+3STeqPDTMgEA7Sl8Dx8A0B4CHwBKgsAHgJIg8AGgJAh8ACgJAh8ASoLAB4CSIPABoCT+H8z35VhuQcYAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lnLs = []\n",
    "ts = []\n",
    "\n",
    "sequence3 = sequence2.replace('C','T')\n",
    "\n",
    "for t in range(0,500):\n",
    "    t= t/100.0\n",
    "    ts.append(t)\n",
    "    lnLs.append(get_sequence_loglikelihood(sequence1,\n",
    "                                           sequence3,\n",
    "                                           t,t))\n",
    "\n",
    "matplotlib.pyplot.scatter(x=ts,y=lnLs)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}