GeneticDistance_Exercise8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum-Likelihood Estimation of a Genetic Distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we're going to walk through the steps of a calculating the likelihood of a genetic distance between two sequences under the Jukes-Cantor model of sequence evolution. In the case of Jukes-Cantor, the probability of change (or no change) given a distance (_t_, in expected number of changes per site) has an analytical derivation (see Yang, Chapter 1). As with most phylogenetic analyses, we will be taking advantage of the assumption that each site evolves independently.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clear()\n",
    "\n",
    "\n",
    "# Calculates the probability that one site does NOT change\n",
    "# Jukes-Cantor model\n",
    "# See Yang, Eq. X.X\n",
    "\n",
    "function probStaySame(Real t){\n",
    "    if (t >= 0.0){\n",
    "        pZero = (1/4) + ((3/4)*exp(-4*(t/3)))\n",
    "        return abs(pZero)\n",
    "    } else {\n",
    "        return abs(0.0)\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "# Calculates the probability that a site DOES change\n",
    "# Jukes-Cantor model\n",
    "# See Yang, Eq. X.X\n",
    "\n",
    "function probChange(Real t){\n",
    "    if (t >= 0.0){\n",
    "        pOne = (1/4) - ((1/4)*exp(-4*(t/3)))\n",
    "        return abs(pOne)\n",
    "    } else {\n",
    "        return abs(0.0)\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at `probStaySame()` and `probChange()` above. Verify that these are correct, based on the expressions given in the first chapter of Yang's book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Returns the entire transition probability matrix for a given branch length\n",
    "# Jukes-Cantor model\n",
    "# We won't actually use the full matrix for these calculations (since we have expressions\n",
    "# for each probability independently, but we would need to use this for more complicated models.)\n",
    "\n",
    "function jcTranMatrix(Real t){\n",
    "    return [[probStaySame(t),probChange(t),probChange(t),probChange(t)],\n",
    "           [probChange(t),probStaySame(t),probChange(t),probChange(t)],\n",
    "           [probChange(t),probChange(t),probStaySame(t),probChange(t)],\n",
    "           [probChange(t),probChange(t),probChange(t),probStaySame(t)]]\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function above uses the individual probabilities to populate the transition matrix for a given _t_. In the case of more complicated models, we do not have expressions available for each of these individual probabilities. Instead, we would need to create an instantaneous rate matrix, _Q_, and use $P(t)=e^{Qt}$ to derive the matrix of transition probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Returns the likelihood for distance t given two sequences\n",
    "# Jukes-Cantor model\n",
    "# Requires probStaySame() and probChange()\n",
    "\n",
    "function distLike(String seqOne, String seqTwo, Real t){\n",
    "    logLike = 0.0\n",
    "    numSites = seqOne.size()\n",
    "    for (i in 1:numSites){\n",
    "        isSame = (seqOne[i] == seqTwo[i])\n",
    "        if (isSame){\n",
    "            logLike += ln(probStaySame(t))\n",
    "        } else {\n",
    "            logLike += ln(probChange(t))\n",
    "        }\n",
    "    }\n",
    "    return logLike\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `distLike()`, given above, takes two strings representing DNA sequences and a _t_ as input. Copy this function into a text editor and turn it into pseudocode. How does the function work? How are we using the assumption of site independence?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# A hill-climbing algorithm to find the maximum likelihood distance t given two sequences\n",
    "# Jukes-Cantor\n",
    "# Requires distLike()\n",
    "\n",
    "function distML(String seqOne, String seqTwo, Real t, Real stepSize){\n",
    "    like = distLike(seqOne,seqTwo,t)\n",
    "    while (distLike(seqOne,seqTwo,t+stepSize) > like){\n",
    "        t = t + stepSize\n",
    "        like = distLike(seqOne,seqTwo,t)\n",
    "    }\n",
    "    while (distLike(seqOne,seqTwo,t-stepSize) > like){\n",
    "        t = t - stepSize\n",
    "        like = distLike(seqOne,seqTwo,t)\n",
    "    }\n",
    "    if(stepSize > 0.0001){\n",
    "        t = distML(seqOne,seqTwo,t,stepSize/2.0)\n",
    "    }\n",
    "    return t\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function above is a hill-climbing algorithm to attempt to find the MLE of genetic distance between two sequences (given as strings or lists). The arguments to the function include:\n",
    "\n",
    "- The first sequence\n",
    "- The second sequence\n",
    "- A starting value for genetic distance, _t_.\n",
    "- The initial step size used to explore different values of _t_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining sequences to calculate distance between\n",
    "\n",
    "seq1 = \"AAGTCCAG\"\n",
    "seq2 = \"AAGCCCCG\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've defined our sequences, let's explore the relationship between different potential genetic distances and likelihood. First, we'll loop through a set of distances between 0.05 and 1 and calculate the associated likelihoods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate likelihoods for different distances (0.05-1.0)\n",
    "for (t in seq(from=0.05,to=1.0,by=0.05)){\n",
    "    distLike(seq1,seq2,t)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do these values look like? What kinds of numbers are these? Roughly where do you think the MLE is? How confident can we be in this MLE?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now we'll use the hill-climbing function to estimate the MLE\n",
    "distML(seq1,seq2,0.05,0.1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RevBayes",
   "language": "bash",
   "name": "revbayes_kernel"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "r"
   },
   "file_extension": ".Rev",
   "help_links": [
    {
     "text": "RevBayes",
     "url": "https://revbayes.org"
    },
    {
     "text": "RevBayes Kernel",
     "url": "https://github.com/sdwfrost/revbayes_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-rsrc",
   "name": "RevBayes",
   "pygments_lexer": "R"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}