{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Introduction to Bioinformatics**<br>\n",
    "A masters course by Blaž Zupan and Tomaž Curk<br>\n",
    "University of Ljubljana (C) 2016-2018\n",
    "\n",
    "Disclaimer: this work is a first draft of our working notes. The images were obtained from various web sites, but mainly from the wikipedia entries with explanations of different biological entities. Material is intended for our bioinformatics class and is not meant for distribution.\n",
    "\n",
    "## Lecture Notes Part 5\n",
    "\n",
    "# Sequence Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Molecular Roots of Evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We, living creatures on earth, evolved. In fact, we are still\n",
    "evolving. The heritable characteristics of our population change over\n",
    "successive generations. When cells divide, the DNA of a cell is\n",
    "replicated: two identical replicates of the DNA are made from one\n",
    "original DNA molecule. Well, not totally identical. When replicating\n",
    "human DNA, for instance, there are at least a few errors left after\n",
    "the replications. Like, on average, there are about seven nucleotides\n",
    "that do not match those from the original DNA. The cell includes\n",
    "machinery to detect and repair such mistakes, but even that machinery\n",
    "sometimes overlooks the mistakes and this could be passed to the\n",
    "organisms of the next generation. These replication mistakes are called\n",
    "mutations. Mutations can be harmless (changes in \"junk\" part\n",
    "of the DNA) or harmful (the organism with a mutated DNA does not\n",
    "survive), but some could be beneficial and make the organism better\n",
    "fit than the old generation. Organisms with better fitness would\n",
    "proliferate throughout the population. This process, called evolution,\n",
    "happens all the time, for all organisms on Earth.\n",
    "\n",
    "Our description of molecular background of evolution is very\n",
    "simplistic. But will suffice for now. In the previous lecture we\n",
    "learned about genes, units of heredity and segments of the DNA that code\n",
    "for proteins. Particularly interesting are mutations that change the\n",
    "coding sequence of the gene, as the resulting protein may\n",
    "have a different sequence of amino acids. It is \"may\", as different\n",
    "sequence of nucleotides may result in the same amino acid chain (think\n",
    "about the number of different codons for each amino acid and\n",
    "redundancy in the DNA-protein coding tables). \n",
    "\n",
    "Two different species of organisms with a common ancestor may\n",
    "therefore have some differences in the genome sequence, and therefore\n",
    "may have (slightly) different genes that encode similar proteins. If\n",
    "these genes are coming from the same ancestor gene, two\n",
    "gene variants are called *orthologous genes*. The difference between\n",
    "genes arises due to nucleotide substitution, insertion of the \n",
    "sequence or deletion of the DNA sequence. The later two types of\n",
    "changes in the DNA are called indels. Another kind of mistake in the DNA\n",
    "replication appears where a part of the DNA is duplicated, that is,\n",
    "mistakenly appears several times in the new genome. As there could be\n",
    "now several genes encoding the same protein, all but one of these\n",
    "genes can be changed with no harm. Multiple copies of the same genes\n",
    "that are then, within a single species, changed accross the evolution,\n",
    "give raise to gene families, a set of similar genes that are all stem from the\n",
    "same predecessor. Such genes are called *parlogous genes*. Orthologs\n",
    "and paralogs are jointly referred as homologous genes.\n",
    "\n",
    "There is a great value of finding homologous genes within the same\n",
    "species or accross the species. To find them, we need to define the\n",
    "distance between sequences and means to express these distances\n",
    "numerically. Sequences of DNA that we would like to compare need first\n",
    "to be aligned to expose the similarities and differences. This\n",
    "procedure is referred to as *sequence alignment*. Sequence alignment\n",
    "and scoring of similarity of sequences give rise to many different\n",
    "procedures of bioinformatics, and help us in tasks such as:\n",
    "\n",
    "* gene function prediction: knowing the function of a gene in one\n",
    "  species, and finding a similar gene in the other species, we can now\n",
    "  make an inference,\n",
    "* gene finding: knowing genome structure in one organism can help us\n",
    "  annotate the gene of another organism\n",
    "* phylogeny: based on comparison of gene sequences of organisms we can\n",
    "  infer which organisms had common predecessors, and we reconstruct\n",
    "  the tree of life,\n",
    "* sequence assembly: knowing the sequence of one genome we can\n",
    "  assemble the sequence of another gene from small fragments of the\n",
    "  DNA that have been sequence using modern high-throughput sequencing\n",
    "  methods.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So What Is Sequence Alignment?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider two sequences $s$ and $t$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "s, t = (\"ATACGTA\", \"TATGATA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and consider their possible alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A1 = (\"ATACGTA-\", \"-TATGATA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us pretty-print the alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pp_alignment(s, t):\n",
    "    print(s)\n",
    "    print(\"\".join(\"|\" if si==ti else \" \" for si, ti in zip(s, t)))\n",
    "    print(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATACGTA-\n",
      " || |   \n",
      "-TATGATA\n"
     ]
    }
   ],
   "source": [
    "pp_alignment(*A1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aligned sequences have the same length. We actually talk about the\n",
    "length of alignment. In our case, $|A_1(s,t)|=c=7$. The symbol \"-\"\n",
    "denotes an indel, a skip in the sequence due to deleted nucleotide in\n",
    "one sequence or added nucleotide in the other one.\n",
    "\n",
    "Consider a different alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATACG-TA\n",
      " || | ||\n",
      "-TATGATA\n"
     ]
    }
   ],
   "source": [
    "A2 = (\"ATACG-TA\", \"-TATGATA\")\n",
    "pp_alignment(*A2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and another one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-AT-ACGTA\n",
      " || |  ||\n",
      "TATGA--TA\n"
     ]
    }
   ],
   "source": [
    "A3 = (\"-AT-ACGTA\", \"TATGA--TA\")\n",
    "pp_alignment(*A3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which alignment is better? We need some means to quantitatively assess the alignments and compute the alignment score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Scoring Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $x$ and $y$ be aligned sequences, that is, $x$ a version of $s$\n",
    "after the alignment, and $y$ a version of $t$ after the\n",
    "alignment. Let $x_i$ be an $i$-th symbol of aligned sequence $x$, and $y_i$ be an\n",
    "$i$-th symbol of aligned sequence $y$. We define a scoring function\n",
    "for position $i$ to be $\\sigma(x_i,y_i)$. The alignment score is then\n",
    "defined as:\n",
    "\n",
    "$$M(A(s,t))=\\sum_{i=1}^c \\sigma(x_i, y_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All we need now is a concrete scoring function $\\sigma$. Let us define one, and score our alignments. We will first use a very simple scoring function, which returns -1 if the symbols are different, and 1 if the symbols at some position of the alignment are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simple_sigma(a, b):\n",
    "    return -1 if a != b else 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def M(A, sigma):\n",
    "    return sum(sigma(xi, yi) for xi, yi in zip(*A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATACGTA-\n",
      " || |   \n",
      "-TATGATA\n",
      "Score: -2\n",
      "\n",
      "ATACG-TA\n",
      " || | ||\n",
      "-TATGATA\n",
      "Score: 2\n",
      "\n",
      "-AT-ACGTA\n",
      " || |  ||\n",
      "TATGA--TA\n",
      "Score: 1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "alignments = [A1, A2, A3]\n",
    "for A in alignments:\n",
    "    pp_alignment(*A)\n",
    "    print(\"Score: %d\\n\" % M(A, simple_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best alignment, according to out scoring function `simple_sigma` is alignment $A_2$. What if change the scoring function? Say, assuming that indels really hurt, because they can change the whole reading frame, we should punish them more than we punish a mismatch of the nucleotides. Let us try is out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def punished_indels(a, b):\n",
    "    return -2 if (a == \"-\" or b == \"-\") else (-1 if a != b else 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, a check if this works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, -2, 2)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "punished_indels(\"A\", \"C\"), punished_indels(\"-\", \"T\"), punished_indels(\"G\", \"G\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fine, looks like it works. Now we use the scoring function `punished_indels` to score the alignments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATACGTA-\n",
      " || |   \n",
      "-TATGATA\n",
      "Score: -1\n",
      "\n",
      "ATACG-TA\n",
      " || | ||\n",
      "-TATGATA\n",
      "Score: 5\n",
      "\n",
      "-AT-ACGTA\n",
      " || |  ||\n",
      "TATGA--TA\n",
      "Score: 2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for A in alignments:\n",
    "    pp_alignment(*A)\n",
    "    print(\"Score: %d\\n\" % M(A, punished_indels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, alignment $A_2$ looks best, but now the margin is bigger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where Are We Now?\n",
    "\n",
    "We can score the alignments. Given the scoring functions. The scoring\n",
    "functions can be more complicated, but for now, the ones we have\n",
    "define will do fine. One problem though. How do we get an alignment?\n",
    "We would actually like to find an alignment $A^*$ that, given a scoring\n",
    "function, maximizes the alignment score $M$, such that there is no\n",
    "other alignment $A$ where $M(A^*)<M(A)$.\n",
    "\n",
    "One think that we could do is check out all possible alignments. There\n",
    "are quite many. Say, for alignment sequence of length $2n$ there are\n",
    "$n$ possible insertion sites. The number of possible alignments is\n",
    "thus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\binom{2n}{n}={(2n)! \\over (n!)^2} \\approx {2^{2n} \\over \\sqrt{\\Pi n}}$$\n",
    "\n",
    "We used the [Stirling's approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation) to factorial. Let's write this in code and see what comes out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import sqrt, pi\n",
    "def foo(n):\n",
    "    return round(2**(2*n) / sqrt(pi * n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(187079, 138710677319)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "foo(10), foo(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for a short sequences of only 20 nucleotides, there are too many possible alignments. No way to use brute force to find optimal alginment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Representation Of A Space of All Possible Alignments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us start with with two very short sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s, t = \"ATGA\", \"TCA\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now construct a table with symbols from s in columns and symbols\n",
    "from t in rows. We also add an indel symbol \"-\" at the first column\n",
    "and row of the table. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![DTP](figs/04-dynamic-programming-table.jpg \"A dynamic programming table.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We claim that any walk from upper left\n",
    "corner to lower right corner of the table, where we can move just one\n",
    "place right, down or diagonally, represents an alignment. Horizontal\n",
    "or vertical move inserts the indel in one of the sequences and\n",
    "consumes the symbol of the other, while a diagonal move consumes a\n",
    "symbol from both sequences. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks complicated, but let us consider an example. We will use the walk that is marked with solid squares. In the first horizontal move we consume \"T\" from $t$ and introduce an indel in $s$. The next move is vertical: we consume a symbol \"A\" from $s$ but nothing from $t$, essentially introducing an indel to aligned sequence $t$. At this stage, our alignment is (-A, T-). Next move is diagonal, we consume T from $s$ and C from $t$. Next is a horizontal move, consuming A from $t$ and introducing indel to $s$. At this point, the alignment is (-AT-, T-CA). The final two moves are vertical, consuming the remaining two letters from $s$ and introducing two indels to $t$. Our final alignment is therefore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-AT-GA\n",
      "      \n",
      "T-CA--\n",
      "Score: -11\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A1 = \"-AT-GA\", \"T-CA--\"\n",
    "pp_alignment(*A1)\n",
    "print(\"Score: %d\\n\" % M(A1, punished_indels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hm, probably not the best alignment, but still a valid one. The path in our table that starts with the same two moves as before, but then takes the dotted path, represents the following alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-ATGA--\n",
      "       \n",
      "T----CA\n",
      "Score: -14\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A2 = \"-ATGA--\", \"T----CA\"\n",
    "pp_alignment(*A2)\n",
    "print(\"Score: %d\\n\" % M(A2, punished_indels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to our scoring function this is an even worse alignment. But again, it is a valid one. And then again we are just trying to convince you, the reader, the walks in the table we have constructed represent all possible alignments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also something else we would like to point out. Both our alignments, A1 and A2, start with the same two moves in the table, a horizontal one and a vertical one, reaching a shadowed square. The two alignment then branch from the shadowed square. Our aligned sequence at that point is (-A, T-), and the score of this partial alignment is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_partial = \"-A\", \"T-\"\n",
    "M(A_partial, punished_indels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this score does not change with the rest of the alignment. That is, whatever follows, the score for this partial alignment with which we started with stays the same. That is simply a result from the fact that our scoring function refers only to a pair of symbols at a certain position, and does not take into consideration the aligned sequence before (or after) this position. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the total alignment score, we can thus break the aligned sequence to subsequences, compute the alignment score for aligned subsequences and sum this up. For alignment A1, we thus can compute the score until the shaded square, and the score for the remaining part of the alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-11"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M((\"-A\", \"T-\"), punished_indels) + M((\"T-GA\", \"CA--\"), punished_indels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same for A2. Again, we have obtained the same result as before, where we have considered the complete aligned sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-14"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M((\"-A\", \"T-\"), punished_indels) + M((\"TGA--\", \"---CA\"), punished_indels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Needleman-Wunsch Algorithm for Optimal Global Alignment\n",
    "\n",
    "Which means, if know the alignment score at the point of the shaded square, this can be used to compute the score for all alignments that start with this particular alignemnt sequence, and we just need to add the score that we get for the alignemnt after the shaded square.\n",
    "\n",
    "Great. But that does not tell us anything about the optimal alignment. Yet, if we would know what the optimal alignment until the shaded square was, any alignment subsequence that would follow the shaded square would precede the optimal alignment until the shaded square.\n",
    "\n",
    "Right! We can now start with empty alignment. This is an upper left corner of our table. And we can now visit all the squares of the table to the right and down. Each square can be reached from three previous positions through a vertical, horizontal or diagonal move. But now that we know that each previous position denotes a starting subsequence of alignment with the score that will just add to our score of the remaining alignment sequence, we can compute the score at each cell by using the score from any of the previous three cells plus the score for the particular move.Horizontal and vertical moves introduce indels, so we add -2 to the score of the corresponding previous cell (-2 due to our `punished_indels` function). Diagonal move consumes symbols from $s$ and $t$, so we add either 2 if the symbols match or -1 if they do not match.\n",
    "\n",
    "Above, we have defined a recursive procedure that can compute the scores for all of the cells of the alignment table. We will refer to the procedure as *dynamic programming*, and to the table as *dynamic programming table*. Let $M_{i,j}$ be the $i$-th row and $j$-th column entry of the table. Then,\n",
    "\n",
    "$$M_{i,j}= \\max\n",
    "\\begin{cases}\n",
    "M_{i-1,j}+\\sigma(s_i,-)\\\\\n",
    "M_{i,j-1}+\\sigma(-,t_i)\\\\\n",
    "M_{i-1,j-1}+\\sigma(s_i,t_i)\n",
    "\\end{cases}$$\n",
    "\n",
    "In this equation, the first case refers to the vertical, the second to the horizontal and the third to the diagonal move.\n",
    "\n",
    "The scores in each of the cells are the scores of optimal alignments at that particular cell, and the path to that particular cell (where the path was chosen to obtain the maximal score) is the path that defines the optimal alignment that would lead us to particular cell. The global alignment, that is, the one that consumes all the symbols from $s$ and $t$ is the alignment where we have reached the lower right cell in the dynamic programming table. If we follow the path that lead us to the score computed in this cell, we obtain the optimal alignment. So simple! For the sequences of length $|s|$ and $|t|$, our computation complexity is only in the order of $|s|\\times|t|$.\n",
    "\n",
    "Above, we forgot about one thing. Intialization. But this is simple. With empy aligment at the start, we are in the upper left cell. The score is 0 (we have not done anything yet). We cal also fill out all the entries in the first row, going from left to right and adding -2, the punishment for introducing the indel. Similar for the first column.\n",
    "\n",
    "The algorithm we have described above is called [Needleman-Wunsch](https://en.wikipedia.org/wiki/Needleman–Wunsch_algorithm) algorithm. It belongs to the [dynamic programming](https://en.wikipedia.org/wiki/Dynamic_programming) type of algorithms. It consist of three stages:\n",
    "\n",
    "* initialization\n",
    "* computation of entries of dynamic programming table\n",
    "* trace-back to reconstruct the optimal alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the alignment of our sequences $s$ and $t$, this is how our dynamic programming table looks like just after the initialization:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|-|T|C|T|A\n",
    "-|-|-|-|-|-|\n",
    "-|0|-2|-4|-6|-8\n",
    "A|-2||||\n",
    "T|-4||||\n",
    "G|-6||||\n",
    "A|-8||||"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this is how it looks like after all its entries have been computed:\n",
    "\n",
    "|-|T|C|T|A\n",
    "-|-|-|-|-|-|\n",
    "-|**0**|-2|-4|-6|-8\n",
    "A|**-2**|-1|-3|-5|-4\n",
    "T|-4|**0**|-2|-1|-3\n",
    "G|-6|-2|**-1**|**-3**|-2\n",
    "A|-8|-4|-3|-2|**-1**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal alignment has a score of -1. We have bolded the cells that are on the path of the contruction of the optimal alignment. Remember, to find this path, we have to start at lower right, and trace back the path, each time taking the previous cell from which we have yielded the maximal score. The optimal alignment corresponding to the path marked in our dynamic programming table is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = \"ATG-A\", \"-TCTA\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us check out if its alignment score really matches the one from our dynamic programming table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M(A, punished_indels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does. Huh. We are fine."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}