Week2_Alignment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sequence alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This week we'll look at some of the alignment algorithms discussed in lectures.\n",
    "\n",
    "If you are new to programming, you can try this [introductory Python tutorial](hammingdist.ipynb) first, as it will give you a more guided introduction to the first exercise on Hamming distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequence data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll read in some real data to play with. Nucleotide (cDNA) sequences for an insulin gene in mice and the insulin gene in humans have been provided in the `data/` directory. \n",
    "\n",
    "We'll use nucleotide sequence rather than protein sequence, because the substitution matrix is very important when aligning protein sequences, and we won't implement substitution matrices today."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">ENST00000397262.5 INS-203 cdna:protein_coding\n",
      "TCCAGGACAGGCTGCATCAGAAGAGGCCATCAAGCAGGTCTGTTCCAAGGGCCTTTGCGT\n",
      "CAGGTGGGCTCAGGATTCCAGGGTGGCTGGACCCCAGGCCCCAGCTCTGCAGCAGGGAGG\n",
      "ACGTGGCTGGGCTCGTGAAGCATGTGGGGGTGAGCCCAGGGGCCCCAAGGCAGGGCACCT\n",
      "GGCCTTCAGCCTGCCTCAGCCCTGCCTGTCTCCCAGATCACTGTCCTTCTGCCATGGCCC\n",
      "TGTGGATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCAGCCG\n",
      "CAGCCTTTGTGAACCAACACCTGTGCGGCTCACACCTGGTGGAAGCTCTCTACCTAGTGT\n",
      "GCGGGGAACGAGGCTTCTTCTACACACCCAAGACCCGCCGGGAGGCAGAGGACCTGCAGG\n",
      "TGGGGCAGGTGGAGCTGGGCGGGGGCCCTGGTGCAGGCAGCCTGCAGCCCTTGGCCCTGG\n",
      "AGGGGTCCCTGCAGAAGCGTGGCATTGTGGAACAATGCTGTACCAGCATCTGCTCCCTCT\n",
      "ACCAGCTGGAGAACTACTGCAACTAGACGCAGCCCGCAGGCAGCCCCACACCCGCCGCCT\n",
      "CCTGCACCGAGAGAGATGGAATAAAGCCCTTGAACCAGC\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Look at the this file with the linux cat command\n",
    "!cat data/Homo_sapiens_INS_203_sequence.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">ENSMUST00000105933.7 Ins2-205 cdna:protein_coding\n",
      "GATCCGCTACAATCAAAAACCATCAGCAAGCAGGAAGGTACTCTTCTCAGTGGGCCTGGC\n",
      "TCCCCAGCTAAGACCTCAGGGACTTGAGGTAGGATATAGCCTCCTCTCTTACGTGAAACT\n",
      "TTTGCTATCCTCAACCCAGCCTATCTTCCAGGTTATTGTTTCAACATGGCCCTGTGGATG\n",
      "CGCTTCCTGCCCCTGCTGGCCCTGCTCTTCCTCTGGGAGTCCCACCCCACCCAGGCTTTT\n",
      "GTCAAGCAGCACCTTTGTGGTTCCCACCTGGTGGAGGCTCTCTACCTGGTGTGTGGGGAG\n",
      "CGTGGCTTCTTCTACACACCCATGTCCCGCCGTGAAGTGGAGGACCCACAAGTGGCACAA\n",
      "CTGGAGCTGGGTGGAGGCCCGGGAGCAGGTGACCTTCAGACCTTGGCACTGGAGGTGGCC\n",
      "CAGCAGAAGCGTGGCATTGTAGATCAGTGCTGCACCAGCATCTGCTCCCTCTACCAGCTG\n",
      "GAGAACTACTGCAACTAGACCCACCACTACCCAGCCTACCCCTCTGCAATGAATAAAACC\n",
      "TTTGAATGAGCA\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!cat data/Mus_musculus_Ins2_205_sequence.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could read the FASTA file using just Python code (note that this code assumes there is only ONE sequence in each FASTA file!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('data/Homo_sapiens_INS_203_sequence.fa') as f:\n",
    "    sequence = \"\"\n",
    "    for row in f.readlines():\n",
    "        if not row.startswith('>'):\n",
    "            sequence += row.strip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TCCAGGACAGGCTGCATCAGAAGAGGCCATCAAGCAGGTCTGTTCCAAGGGCCTTTGCGTCAGGTGGGCTCAGGATTCCAGGGTGGCTGGACCCCAGGCCCCAGCTCTGCAGCAGGGAGGACGTGGCTGGGCTCGTGAAGCATGTGGGGGTGAGCCCAGGGGCCCCAAGGCAGGGCACCTGGCCTTCAGCCTGCCTCAGCCCTGCCTGTCTCCCAGATCACTGTCCTTCTGCCATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCAGCCGCAGCCTTTGTGAACCAACACCTGTGCGGCTCACACCTGGTGGAAGCTCTCTACCTAGTGTGCGGGGAACGAGGCTTCTTCTACACACCCAAGACCCGCCGGGAGGCAGAGGACCTGCAGGTGGGGCAGGTGGAGCTGGGCGGGGGCCCTGGTGCAGGCAGCCTGCAGCCCTTGGCCCTGGAGGGGTCCCTGCAGAAGCGTGGCATTGTGGAACAATGCTGTACCAGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAACTAGACGCAGCCCGCAGGCAGCCCCACACCCGCCGCCTCCTGCACCGAGAGAGATGGAATAAAGCCCTTGAACCAGC\n"
     ]
    }
   ],
   "source": [
    "print(sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, a better practice usually, we could use a standard library such as scikit-bio to parse the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skbio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequences = skbio.io.read('data/Homo_sapiens_INS_203_sequence.fa', format='fasta')\n",
    "# This gives us a skbio Sequence object\n",
    "human_ins_object = list(sequences)[0]\n",
    "# This gives us the actual sequence as a string\n",
    "human_ins = str(human_ins_object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat for mouse INS2 gene\n",
    "mouse_ins = str(list(skbio.io.read('data/Mus_musculus_Ins2_205_sequence.fa', format='fasta'))[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Edit distance "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1: Hamming distance**\n",
    "\n",
    "Edit the Hamming distance function below so that it returns the correct Hamming distance for two strings `a` and `b`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamming(a,b):\n",
    "    \"\"\"\n",
    "    Calculate the Hamming distance between strings a and b.\n",
    "    The strings must be the same length.\n",
    "    \"\"\"\n",
    "    # 0 is the wrong answer. Edit this function to give the right answer.\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think also: what will your function do if the strings are of different length? What *should* it do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 2\n",
    "hamming(\"GATTACA\",\"GACTATA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 6\n",
    "hamming(\"tuesday\",\"sundays\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These strings are of different length!\n",
    "hamming(\"happiness\",\"applying\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2: Levenshtein distance** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Edit the `lev` function below to calculate Levenshtein distance recursively. You can use the costs \n",
    "* 1 for an indel\n",
    "* 1 for a mismatch\n",
    "* 0 for a match\n",
    "\n",
    "This is the same function as shown during lectures, but try to implement it without looking back at the slides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lev(a,b):\n",
    "    \"\"\"\n",
    "    Recursively calculate Levenshtein distance between strings a and b.\n",
    "    \"\"\"\n",
    "    if len(a)==0:\n",
    "        return len(b)\n",
    "    if len(b)==0:\n",
    "        return len(a)\n",
    "    # Add code below to calculate the distance in terms of \n",
    "    # lev(a,b[1:]), lev(a[1:],b), and lev(a[1:],b[1:])\n",
    "    # We're interested in the incremental cost of aligning the current characters,\n",
    "    # a[0] and b[0]\n",
    "    return 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 2\n",
    "lev(\"GATTACA\",\"GACTATA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 4\n",
    "lev(\"tuesday\",\"sundays\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 6\n",
    "lev(\"happiness\",\"applying\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment scores, global and local alignment "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a function which implements a recursive alignment function like `lev()`, but returns an alignment score rather than an edit distance. Notice that it uses `max()` rather than `min()`, as we're trying to find the maximum possible score, not the minimum possible edit distance.\n",
    "\n",
    "Avoid looking at the function below until you've done Exercise 2 above as it may give away the answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def align_recursive(a, b, indel_score=-1, match_score=2, mismatch_score=-1):\n",
    "    \"\"\"\n",
    "    Recursively calculate alignment score between strings a and b,\n",
    "    using supplied scores for matches, mismatches and indels.\n",
    "    \"\"\"\n",
    "    if len(a)==0:\n",
    "        return len(b)\n",
    "    if len(b)==0:\n",
    "        return len(a)\n",
    "    if a[0]==b[0]:\n",
    "        match_mismatch_score = match_score\n",
    "    else:\n",
    "        match_mismatch_score = mismatch_score\n",
    "    return max(align_recursive(a[1:],b) + indel_score,\n",
    "               align_recursive(a,b[1:]) + indel_score,\n",
    "               align_recursive(a[1:],b[1:]) + match_mismatch_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that you can now change the scoring system for indels, matches, and mismatches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "align_recursive(\"GATTACA\",\"GACTATA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "align_recursive(\"GATTACA\",\"GACTATA\",match_score=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, this function is recursive and will be slow for large strings. `align_recursive(human_ins, mouse_ins)` is not practical to run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global alignment: Needleman-Wunsch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do global alignment with the Needleman-Wunsch algorithm, we need two steps:\n",
    "\n",
    "1. Fill out the grid of alignment scores. This is enough to give the final alignment score.\n",
    "2. Trace-back from the bottom-right corner of the grid to get the actual alignment of the strings.\n",
    "\n",
    "Here, we've provided a function to do the traceback, and given an incomplete function to calculate the alignment score grid. Complete in the `calculate_scoregrid()` function to correctly fill out the grid of scores.\n",
    "\n",
    "For traceback, we have two options:\n",
    "* Keep track of which cell(s) was/were the origin of the best score(s) for each given cell, and use this information for traceback. This increases storage requirements by a constant factor (i.e. they are still O(N^2)).\n",
    "* Or, during traceback, calculate which cells(s) could have been the origin of the best score(s) for each cell. This increases the computational cost of traceback by a constant factor (i.e. it is still O(N)).\n",
    "\n",
    "In this case, the provided traceback function will work out which path to follow, so you don't need to keep track of the path as you calculate the scores."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** Complete the `calculate_scoregrid()` function to calculate the scores needed for global alignment via the Needleman-Wunsch algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_scoregrid(a, b,\n",
    "                        indel_score=-1, match_score=2, mismatch_score=-1):\n",
    "    \"\"\"\n",
    "    Given two strings a and b, calculate the maximum score grid, using\n",
    "    specified scores for indels, matches and mismatches. Return the grid.\n",
    "    Grid row and column 0 correspond to \"before\" the start of each string,\n",
    "    so grid indexes are offset by 1 from string indexes. That is,\n",
    "    grid position [1,1] represents the result of matching a[0] to b[0].\n",
    "    \"\"\"\n",
    "    # The grid needs to be 1 bigger in each direction than the string lengths\n",
    "    X = len(a)+1\n",
    "    Y = len(b)+1\n",
    "    # Initialise the grid with zeroes\n",
    "    scoregrid = np.zeros((X,Y), np.int)\n",
    "    #### YOUR CODE HERE #### \n",
    "    # You need to:\n",
    "    # * initialise the top edge of grid, i.e. scoregrid[x,0] for all x, with indel scores\n",
    "    # * initialise the left edge of grid, i.e. scoregrid[0,y] for all x, with indel scores\n",
    "    # * loop over x and y, filling out each cell of the grid by looking for the\n",
    "    #   maximum possible score from each of the three earlier cells\n",
    "    \n",
    "    return scoregrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pre-defined functions to get the traceback given a correct scoregrid\n",
    "# Use help(traceback) or help(get_alignment) to see how to call them\n",
    "from alignment_functions import traceback, get_alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `calculate_scoregrid()` works correctly, the below will work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = \"GATTACA\"\n",
    "b = \"GACTATA\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [0, 0, 0, 0, 0, 0, 0, 0]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Once you've implemented calculate_scoregrid, this should show the correct\n",
    "# values instead of all zeroes\n",
    "scoregrid = calculate_scoregrid(a,b)\n",
    "scoregrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alignment score: 0\n"
     ]
    }
   ],
   "source": [
    "print(\"Alignment score:\",scoregrid[-1,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the score grid isn't correct and consistent with the scoring system\n",
    "# and the strings, traceback won't be able to find a path and will give an error\n",
    "trace = traceback(a,b,scoregrid)\n",
    "aligned_string_a, aligned_string_b = get_alignment(trace)\n",
    "print(aligned_string_a)\n",
    "print(aligned_string_b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try aligning the cDNA strings `human_ins` and `mouse_ins`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Challenge exercise 4:** Change your `calculate_scoregrid()` function to perform local instead of global alignment. You can import the `traceback_local()` function to help you test your result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from alignment_functions import traceback_local"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}