Assignment2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COMPSCI 369 – Assignment 2 – Axton Pitt (2671457), Friday May 13th\n",
    "\n",
    "## Problem 1 – Jukes-Cantor Model of DNA sequence evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1a) \n",
    "Below is a function for simulating the Jukes-Cantor Model of DNA sequence evolution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTAAGCGTGACGCGCGCGGTAGTACGCGGTGGTCCTTTAGAGTACATCGC\n",
      "CTAAGCGTGACGGGCCCGGTAGTACGCGGTGGTCCTTCATAGTACATCGA\n",
      "CTAAGCGTGATGCGCCCAGTAGTACGAGGTGGTCCTTTAGAGTACATAGC\n",
      "Number of sites at which ancestral and decendent 1 differ:\n",
      "5\n",
      "Number of sites at which ancestral and decendent 2 differ:\n",
      "5\n",
      "Number of sites at which decendent 1 and decendent 2 differ:\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "def jukesCantorSimulation(sequenceLength, time, mutationRate):\n",
    "\n",
    "    # generate random sequence of length perscribed\n",
    "    from random import randint\n",
    "    from random import random\n",
    "    from math import log\n",
    "\n",
    "    ancestoralSequence = []\n",
    "    decendentSequence1 = []\n",
    "    decendentSequence2 = []\n",
    "\n",
    "    # create DNA seq from 0,1,2,3 as\n",
    "    for index in range(0, sequenceLength):\n",
    "        randomBase = randint(0,3)\n",
    "        ancestoralSequence.append(randomBase)\n",
    "        # print(ancestoralSequence)\n",
    "\n",
    "    decendentSequence1 = ancestoralSequence.copy()\n",
    "    decendentSequence2 = ancestoralSequence.copy()\n",
    "\n",
    "    time1 = time\n",
    "    time2 = time\n",
    "\n",
    "    # create decendentSequence1\n",
    "    while time1 > 0:\n",
    "\n",
    "        # get a uniform random number representing U(0, 1)\n",
    "        lengthOfTime = random()\n",
    "\n",
    "        # sum probability of mutation\n",
    "        sumProbability = sequenceLength * mutationRate * 0.75\n",
    "\n",
    "        # turn this into time to mutation\n",
    "        timeMinus = -log(lengthOfTime) / sumProbability\n",
    "\n",
    "        if time1 <= 0:\n",
    "            break\n",
    "        else:\n",
    "            bases = [0, 1, 2, 3]\n",
    "\n",
    "            # index of base to mutate\n",
    "            baseIndex = randint(0, sequenceLength-1)\n",
    "            candidateIndex = randint(0, 3)\n",
    "            candidateBase = bases[candidateIndex]\n",
    "\n",
    "            # get current base\n",
    "            currentBase = decendentSequence1[baseIndex]\n",
    "\n",
    "            while currentBase == candidateBase:\n",
    "\n",
    "                candidateIndex = randint(0, 3)\n",
    "                candidateBase = bases[candidateIndex]\n",
    "\n",
    "            decendentSequence1[baseIndex] = candidateBase\n",
    "\n",
    "        time1 -= timeMinus\n",
    "\n",
    "    # create decendentSequence2\n",
    "    while time2 > 0:\n",
    "\n",
    "        # get a uniform random number representing U(0, 1)\n",
    "        lengthOfTime = random()\n",
    "\n",
    "        # sum probability of mutation\n",
    "        sumProbability = sequenceLength * mutationRate * 0.75\n",
    "\n",
    "        # turn this into time to mutation\n",
    "        timeMinus = -log(lengthOfTime) / sumProbability\n",
    "\n",
    "        if time2 <= 0:\n",
    "            break\n",
    "        else:\n",
    "            bases = [0, 1, 2, 3]\n",
    "\n",
    "            # index of base to mutate\n",
    "            baseIndex = randint(0, sequenceLength-1)\n",
    "            candidateIndex = randint(0, 3)\n",
    "            candidateBase = bases[candidateIndex]\n",
    "\n",
    "            # get current base\n",
    "            currentBase = decendentSequence2[baseIndex]\n",
    "\n",
    "            while currentBase == candidateBase:\n",
    "                candidateIndex = randint(0, 3)\n",
    "                candidateBase = bases[candidateIndex]\n",
    "\n",
    "            decendentSequence2[baseIndex] = candidateBase\n",
    "\n",
    "        time2 -= timeMinus\n",
    "\n",
    "    outputSequences = []\n",
    "\n",
    "    # convert to bases\n",
    "    # return ancestral sequence\n",
    "    letterArrayAncestoralSequence = numberToBaseConversion(ancestoralSequence)\n",
    "\n",
    "    # return descendant sequences\n",
    "    letterArrayDecendentSequence1 = numberToBaseConversion(decendentSequence1)\n",
    "    letterArrayDecendentSequence2 = numberToBaseConversion(decendentSequence2)\n",
    "\n",
    "    outputSequences.append(letterArrayAncestoralSequence)\n",
    "    outputSequences.append(letterArrayDecendentSequence1)\n",
    "    outputSequences.append(letterArrayDecendentSequence2)\n",
    "\n",
    "    return outputSequences\n",
    "\n",
    "def numberToBaseConversion(array):\n",
    "    newLetterArray = []\n",
    "\n",
    "    for item in array:\n",
    "        if item == 0:\n",
    "            newLetterArray.append(\"A\")\n",
    "        elif item == 1:\n",
    "            newLetterArray.append(\"C\")\n",
    "        elif item == 2:\n",
    "            newLetterArray.append(\"G\")\n",
    "        elif item == 3:\n",
    "            newLetterArray.append(\"T\")\n",
    "\n",
    "    return newLetterArray\n",
    "\n",
    "def baseToNumberConversion(array):\n",
    "    newArray = []\n",
    "\n",
    "    for item in array:\n",
    "        if item == \"A\":\n",
    "            newArray.append(0)\n",
    "        elif item == \"C\":\n",
    "            newArray.append(1)\n",
    "        elif item == \"G\":\n",
    "            newArray.append(2)\n",
    "        elif item == \"T\":\n",
    "            newArray.append(3)\n",
    "\n",
    "    return newArray\n",
    "\n",
    "def printSequence(array):\n",
    "    outputString = \"\"\n",
    "    \n",
    "    for item in array:\n",
    "        outputString += item\n",
    "        \n",
    "    print(outputString)\n",
    "\n",
    "def numberOfBaseDifferences(sequence1, sequence2):\n",
    "    minimumSequence = min(sequence1, sequence2)\n",
    "    minimumLength = len(minimumSequence)\n",
    "    diff = 0\n",
    "    # return number of sites at which sequence differs\n",
    "    for index in range(0, minimumLength):\n",
    "        if sequence1[index] != sequence2[index]:\n",
    "            diff += 1\n",
    "    return diff\n",
    "\n",
    "# run with sequence length of 50, 0.01 mutation rate, and time = 10\n",
    "sequences = jukesCantorSimulation(50, 0.01, 10)\n",
    "\n",
    "# print sequences\n",
    "for item in sequences:\n",
    "    printSequence(item)\n",
    "\n",
    "# print differences\n",
    "print(\"Number of sites at which ancestral and decendent 1 differ:\")\n",
    "print(numberOfBaseDifferences(sequences[0], sequences[1]))\n",
    "\n",
    "print(\"Number of sites at which ancestral and decendent 2 differ:\")\n",
    "print(numberOfBaseDifferences(sequences[0], sequences[2]))\n",
    "\n",
    "print(\"Number of sites at which decendent 1 and decendent 2 differ:\")\n",
    "print(numberOfBaseDifferences(sequences[1], sequences[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1b) \n",
    "You would expect the average number of mutations that occur to be Poisson distributed equal to parameter \\[2tL\\frac{3}{4}\\mu\\]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "296.399\n",
      "206.23379899999418\n"
     ]
    }
   ],
   "source": [
    "def distributionSimulation():\n",
    "\n",
    "    # repeat 1000 times\n",
    "    # run with sequence length of 1000, 0.01 mutation rate, and time = 25\n",
    "    count = 0\n",
    "    totalCount = 0\n",
    "    totalSquaredSum = 0\n",
    "\n",
    "    for index in range(0, 1000):\n",
    "        simulatedSequences = jukesCantorSimulation(1000, 0.01, 25)\n",
    "        count = numberOfBaseDifferences(simulatedSequences[1], simulatedSequences[2])\n",
    "        # print(count)\n",
    "        totalCount += count\n",
    "        totalSquaredSum += count**2\n",
    "\n",
    "\n",
    "    # determine mean and variance\n",
    "    mean = totalCount / 1000\n",
    "    variance = (totalSquaredSum / 1000) - mean**2\n",
    "\n",
    "    print(mean)\n",
    "    print(variance)\n",
    "\n",
    "distributionSimulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the distribution to be Poisson distributed, the mean and variance should be equal to \\(2tL\\frac{3}{4}\\mu\\). As \\(2tL\\frac{3}{4}\\mu\\) = \\(2(25)(1000)\\frac{3}{4}(0.01)\\) = 375. As the mean is 296.399 and variance is 206.234 it looks like the number of mutations that occur is not Poisson distributed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1c) \n",
    "Simulating a pair of sibling sequences, B and C, of length 10000 with time = 10 units and mutation rate 0.03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability Pab a = b:\n",
      "0.6656262192742879\n",
      "Probability Pab a =/= b:\n",
      "0.11145792690857069\n"
     ]
    }
   ],
   "source": [
    "def numberOfSpecificBaseDifferences(sequence1, sequence2):\n",
    "    minimumSequence = min(sequence1, sequence2)\n",
    "    minimumLength = len(minimumSequence)\n",
    "\n",
    "    results = []\n",
    "    AtoA = 0\n",
    "    AtoT = 0\n",
    "    AtoC = 0\n",
    "    AtoG = 0\n",
    "    totalAto = 0\n",
    "\n",
    "    # return number of sites at which sequence differs\n",
    "    for index in range(0, minimumLength):\n",
    "        if sequence1[index] == \"A\" and sequence2[index] == \"A\":\n",
    "            AtoA += 1\n",
    "        elif sequence1[index] == \"A\" and sequence2[index] == \"T\":\n",
    "            AtoT += 1\n",
    "        elif sequence1[index] == \"A\" and sequence2[index] == \"C\":\n",
    "            AtoC += 1\n",
    "        elif sequence1[index] == \"A\" and sequence2[index] == \"G\":\n",
    "            AtoG += 1\n",
    "\n",
    "    totalAto = AtoA + AtoT + AtoC + AtoG\n",
    "\n",
    "    results.append(AtoA)\n",
    "    results.append(AtoT)\n",
    "    results.append(AtoC)\n",
    "    results.append(AtoG)\n",
    "    results.append(totalAto)\n",
    "\n",
    "    return results\n",
    "\n",
    "\n",
    "def calculateProbabilityAB():\n",
    "    import numpy\n",
    "\n",
    "    oneCSequences = jukesCantorSimulation(10000, 10, 0.03)\n",
    "    b = oneCSequences[1]\n",
    "    c = oneCSequences[2]\n",
    "\n",
    "    values = numberOfSpecificBaseDifferences(b, c)\n",
    "    estimatePABSame = values[0] / values[4] # theoretical value is 0.6616\n",
    "    estimatePABTotalDiff = (values[1] + values[2] + values[3]) / values[4] # theoretical value is 0.112797\n",
    "    estimatePABDiff = estimatePABTotalDiff / 3\n",
    "    print(\"Probability Pab a = b:\")\n",
    "    print(estimatePABSame)\n",
    "    print(\"Probability Pab a =/= b:\")\n",
    "    print(estimatePABDiff)\n",
    "\n",
    "calculateProbabilityAB()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimates above a similar to the theoretical values of 0.6616 and 0.112797"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2 – Jukes-Cantor Model of DNA sequence evolution with insertions and deletions\n",
    "\n",
    "### 2a)\n",
    "Simulation extended to include insertions and deletions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A---TAC---TAATTACAATGTCATAG---AAGGCAGATCAGGGTAGCTA\n",
      "T---TGC---TATTTACGATATTAGAG---AAGGAAGATCAAGATAGTTG\n"
     ]
    }
   ],
   "source": [
    "def simulationWithInsertionsAndDeletions(sequence1, sequence2, sequenceLength, time, mutationRate):\n",
    "    import random\n",
    "\n",
    "    decendentSequence1 = sequence1\n",
    "    decendentSequence2 = sequence2\n",
    "\n",
    "    # draw a Poisson random variate hi with rate Lt\\mu/10 - for number of insertions\n",
    "    rate = (sequenceLength * time * mutationRate) / 10\n",
    "    noOfInsertions = poissonDistribution(rate, 1) # T = 1\n",
    "\n",
    "    # draw a Poisson random variate hd with rate Lt\\mu/10 - for number of deletions\n",
    "    noOfDeletions = poissonDistribution(rate, 1) # T = 1\n",
    "\n",
    "    # print(noOfInsertions)\n",
    "    # print(noOfDeletions)\n",
    "\n",
    "    for count in range(0, noOfInsertions):\n",
    "        inserts = insertionSequence()\n",
    "        # add insertions to sequence\n",
    "        indexToInsert = random.randint(0, sequenceLength)\n",
    "        indexToInsertP1 = indexToInsert + 1\n",
    "        indexToInsertP2 = indexToInsertP1 + 1\n",
    "        decendentSequence1.insert(indexToInsert, inserts[0])\n",
    "        decendentSequence1.insert(indexToInsertP1, inserts[1])\n",
    "        decendentSequence1.insert(indexToInsertP2, inserts[2])\n",
    "        decendentSequence2.insert(indexToInsert, inserts[0])\n",
    "        decendentSequence2.insert(indexToInsertP1, inserts[1])\n",
    "        decendentSequence2.insert(indexToInsertP2, inserts[2])\n",
    "\n",
    "    for count in range(0, noOfDeletions):\n",
    "        deletion = \"-\"\n",
    "        # add deletions to sequence\n",
    "        indexToInsert = random.randint(0, sequenceLength)\n",
    "        indexToInsertP1 = indexToInsert + 1\n",
    "        indexToInsertP2 = indexToInsertP1 + 1\n",
    "        decendentSequence1[indexToInsert] = deletion\n",
    "        decendentSequence1[indexToInsertP1] = deletion\n",
    "        decendentSequence1[indexToInsertP2] = deletion\n",
    "        decendentSequence2[indexToInsert] = deletion\n",
    "        decendentSequence2[indexToInsertP1] = deletion\n",
    "        decendentSequence2[indexToInsertP2] = deletion\n",
    "\n",
    "    printSequence(decendentSequence1)\n",
    "    printSequence(decendentSequence2)\n",
    "\n",
    "    return [decendentSequence1, decendentSequence2]\n",
    "\n",
    "def insertionSequence():\n",
    "    import random\n",
    "    insertionSeq = []\n",
    "    bases = [\"A\", \"T\", \"C\", \"G\"]\n",
    "    for count in range(0, 3):\n",
    "        baseIndex = random.randint(0,3)\n",
    "        base = bases[baseIndex]\n",
    "        insertionSeq.append(base)\n",
    "\n",
    "    return insertionSeq\n",
    "\n",
    "\n",
    "\n",
    "def poissonDistribution(lbd, T):\n",
    "    times = []\n",
    "    t = randexp(lbd)\n",
    "    while t < T: # draw exp(lambda) variables until we are past T\n",
    "        times.append(t)\n",
    "        t += randexp(lbd)\n",
    "    return len(times)\n",
    "\n",
    "def randexp(lbd): # lambda is a keyword in Python\n",
    "    from math import log\n",
    "    from random import random\n",
    "    return - log(random()) / lbd\n",
    "\n",
    "newSequences = jukesCantorSimulation(50, 20, 0.01)\n",
    "sequences = simulationWithInsertionsAndDeletions(newSequences[1], newSequences[2], 50, 20, 0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b) \n",
    "If this model was used, it would alter the sequence length. This would not result in a Poisson distributed number of insertions and deletions with parameter 2Ltμ/10, as the total sequence length has the potential to change and alter the total probability of mutation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3 – Overlap alignment algorithm\n",
    "\n",
    "### 3a)\n",
    "Implementation of overlap alignment algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def overlapAlignment(sequence1, sequence2, scoreMatrix, gapPenalty):\n",
    "    import numpy\n",
    "\n",
    "    fMatrix = numpy.zeros((len(sequence1), len(sequence2)))\n",
    "\n",
    "    # construct the matrix\n",
    "    for i in range(0, len(sequence1)):\n",
    "        fMatrix[i][0] = 0*i # as this is overlap do not need penalty here\n",
    "    for j in range(0, len(sequence2)):\n",
    "        fMatrix[j][0] = 0*j\n",
    "    for i in range(1, len(sequence1)):\n",
    "        for j in range(1, len(sequence2)):\n",
    "            matchValue = fMatrix[i-1][j-1] + scoreMatrix[sequence1[i]][sequence2[j]]\n",
    "            deleteValue = fMatrix[i-1][j] + gapPenalty\n",
    "            insertValue = fMatrix[i][j-1] + gapPenalty\n",
    "            fMatrix[i][j] = max(matchValue, deleteValue, insertValue)\n",
    "\n",
    "    # backtracking and form alignment - need to get max from last column and row\n",
    "    maxX = 0\n",
    "    maxY = 0\n",
    "\n",
    "    for i in range(0, len(sequence1)):\n",
    "        if fMatrix[i][0] > maxX:\n",
    "            maxX = fMatrix[i][0]\n",
    "\n",
    "    for j in range(0, len(sequence2)):\n",
    "        if fMatrix[0][j] > maxY:\n",
    "            maxY = fMatrix[0][j]\n",
    "\n",
    "    alignment1 = \"\"\n",
    "    alignment2 = \"\"\n",
    "    i = maxX\n",
    "    j = maxY\n",
    "\n",
    "    while i > 0 or j > 0:\n",
    "        if i > 0 and j > 0 and fMatrix[i][j] == fMatrix[i-1][j-1] + scoreMatrix[sequence1[i]][sequence2[j]]:\n",
    "            alignment1 = sequence1[i] + alignment1\n",
    "            alignment2 = sequence2[j] + alignment2\n",
    "            i -= 1\n",
    "            j -= 1\n",
    "        elif i > 0 and fMatrix[i][j] == fMatrix[i-1][j] + gapPenalty:\n",
    "            alignment1 = sequence1[i] + alignment1\n",
    "            alignment2 = \"-\" + alignment2\n",
    "            i -= 1\n",
    "        elif j > 0 and fMatrix[i][j] == fMatrix[i][j-1] + gapPenalty:\n",
    "            alignment1 = \"-\" + alignment1\n",
    "            alignment2 = sequence2[j] + alignment2\n",
    "            j -= 1\n",
    "\n",
    "    print(alignment1)\n",
    "    print(alignment2)\n",
    "\n",
    "import numpy\n",
    "\n",
    "newSequences = jukesCantorSimulation(50, 20, 0.01)\n",
    "sequences = simulationWithInsertionsAndDeletions(newSequences[1], newSequences[2], 50, 20, 0.01)\n",
    "scoreMatrix = numpy.zeros((len(sequences[0]), len(sequences[1])))\n",
    "for i in range(1, len(sequences[0])):\n",
    "    for j in range(1, len(sequences[1])):\n",
    "        scoreMatrix[i][j] = 2\n",
    "\n",
    "overlapAlignment(sequences[0], sequences[1], scoreMatrix, -3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}