HW 2 - Sequence Alignment/Starter/notebook.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-fd91a97c0a7166df",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "# BMI\/CS 576 Fall 2019 - HW2\n",
                "The objectives of this homework are to practice\n",
                "\n",
                "* implementing sequence alignment algorithms\n",
                "* devising dynamic programming algorithms\n",
                "* working with probabilistic models\n",
                "\n",
                "## HW policies\n",
                "Before starting this homework, please read over the [homework policies](https:\/\/canvas.wisc.edu\/courses\/167969\/pages\/hw-policies) for this course.  In particular, note that homeworks are to be completed *individually*.\n",
                "\n",
                "You are welcome to use any code from the weekly notebooks in your solutions to the HW."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-e632c3d68ef6e040",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 1: Affine gap global alignment with homopolymer indels (60 points)\n",
                "\n",
                "### Homopolymer indels\n",
                "Some cutting-edge DNA sequencing technologies, such as Oxford Nanopore, perform poorly in regions of sequences that are *homopolymers*, which are runs of consecutive and identical bases.  For example, in the DNA sequence `GCTAGCCCCCTATC`, the substring `CCCCC` is a homopolymer of the base `C`.  Some sequencers have difficulty in accurately determining the length of such homopolymer regions.  If the example sequence just given were the true sequence, a sequencer could make a mistake and report a sequence of `GCTAGCCCTATC` (truncating the homopolymer region to three `C`s) or `GCTAGCCCCCCTATC` (expanding the homopolyer region to six `C`s), just to give a few examples of possible errors.\n",
                "\n",
                "When aligning pairs of DNA sequences (say for the purpose of determining overlapping reads in genome assembly), a homopolymer error in a sequence results in an insertion or deletion in that sequence with respect to the other.  Given the knowledge that homopolymer errors are common in sequences from certain sequencing technologies, we may wish to use an alignment scoring scheme that penalizes gaps resulting from homopolymer errors *to a lesser degree* than other, non-homopolymer gaps.\n",
                "\n",
                "### An algorithm for pairwise alignment with homopolymer indels\n",
                "In this problem you will implement the dynamic programming algorithm for global alignment with affine gap penalties that distinguishes between homopolymer and non-homopolymer insertions and uses different penalties for the two cases.  There are many ways in which this could be approached.  We will use one specific approach that considers a gap to be a *homopolymer gap* if the inserted sequence is a homopolymer and the base preceding the inserted sequence is the same as the base within the homopolymer.  Note that in pairwise alignment, a deletion in one sequence is indistinguishable from an insertion in the second sequence, and therefore we will consider all gaps to be insertions.  We will score a homopolymer gap of length $k$ using the  function:\n",
                "$$w_h(k) = h + tk$$\n",
                "and non-homopolymer gaps with the function\n",
                "$$w(k) = g + sk,$$\n",
                "and we will generally have that $h > g$ and $t > s$ such that homopolymer gaps penalized less than non-homopolymer gaps (note that all parameters are negative).\n",
                "\n",
                "It turns out that we can find an optimal alignment with this scoring scheme using a modification of the standard affine gap global alignment dynamic programming algorithm.  Specifically, we will introduce two additional matrices, $H_x$ and $H_y$, which will keep track of alignments that end with a homopolymer insertion in $x$ and $y$, respectively.  To be precise, $H_x[i, j]$ will be defined as the optimal score of an alignment of the first $i$ characters of $x$ and the first $j$ characters of $y$ that ends in a homopolymer insertion in $x$.  The $H_x$ and $H_y$ matrices are similar to the $I_x$ and $I_y$ matrices in meaning, but are only for alignments ending in homopolymer gaps.  The dynamic programming recurrences thus become:\n",
                "\n",
                "$M(i, j) = \\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "M(i - 1, j - 1) + S(x_i, y_j) \\\\\n",
                "I_x(i - 1, j - 1) + S(x_i, y_j) \\\\\n",
                "I_y(i - 1, j - 1) + S(x_i, y_j) \\\\\n",
                "H_x(i - 1, j - 1) + S(x_i, y_j) \\\\\n",
                "H_y(i - 1, j - 1) + S(x_i, y_j) \\\\\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "$I_x(i, j) = \\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "M(i - 1, j) + g + s \\\\\n",
                "I_x(i - 1, j) + s \\\\\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "$I_y(i, j) = \\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "M(i, j - 1) + g + s \\\\\n",
                "I_y(i, j - 1) + s \\\\\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "$H_x(i, j) = \\delta(x[i] \\neq x[i - 1]) + \n",
                "\\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "M(i - 1, j) + h + t \\\\\n",
                "H_x(i - 1, j) + t \\\\\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "$H_y(i, j) = \\delta(y[j] \\neq y[j - 1]) +\n",
                "\\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "M(i, j - 1) + h + t \\\\\n",
                "H_y(i, j - 1) + t \\\\\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "In the recurrences for $H_x$ and $H_y$, the function $\\delta$, is defined as \n",
                "$$\\delta(condition) = \\left\\{\n",
                "\\begin{array}{ll}\n",
                "-\\infty, &\\textrm{if}\\ condition \\\\\n",
                "0, & \\textrm{otherwise}\n",
                "\\end{array}\n",
                "\\right.$$\n",
                "which serves to restrict the scores in the $H$ matrices to those of alignments that end in homopolymer indels.\n",
                "\n",
                "Note that the definition of a homopolymer gap provided above does *not* include insertions of a homopolymer at the beginning of either sequence, since there would be no base preceding such an insertion.  Therefore, for initialization, the entries of the first row and column of both the $H_x$ and $H_y$ matrices are set to $-\\infty$.\n",
                "\n",
                "### Your task\n",
                "Implement the algorithm described above as a function `align_global_affine_hp_gaps` below, that takes as input two sequences, `x` and `y`, a substitution matrix, the gap scoring parameters ($g$, $h$, $s$, and $t$).  The substitution matrix will be represented as a dictionary with two-element tuples, `(a, b)`, as keys and scores as values.  Your function should output a tuple of two elements, the first being the score of an optimal alignment, and the second being a single alignment that obtains that score. Your alignment should be represented as a list of two strings.  See the \"Tests for PROBLEM 1\" section at the bottom of this notebook for examples of the inputs and outputs.  \n",
                "\n",
                "In the case that there are multiple optimal alignments, during the traceback, if there are ties for which matrix to jump back to at each step, the order of preference for which matrix to jump to should be ($I_x$, $H_x$, $M$, $I_y$, $H_y$).  For example, in a $H_y$ cell, if the $M$ and $H_y$ cases in the recurrence both tie for the maximum, then during the traceback, you should prefer to jump back to the $M$ cell.\n",
                "\n",
                "Your implementation must use an efficient (polynomial-time) dynamic programming algorithm (i.e., either top-down or bottom-up)."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-0db645b3e786815e",
                    "locked": false,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "outputs": [],
            "source": [
                "# Code for PROBLEM 1\n",
                "# You are welcome to develop your code as a separate Python module\n",
                "# and import it here if that is more convenient for you.\n",
                "\n",
                "# You may find the following constant value for negative infinity of use\n",
                "NEGATIVE_INFINITY = float(\"-inf\")\n",
                "\n",
                "\n",
                "def align_global_affine_hp_gaps(x, y, submatrix, g, h, s, t):\n",
                "    \"\"\"Computes an optimal global pairwise alignment \n",
                "    with an affine gap homopolymer-aware scoring function.\n",
                "        \n",
                "    Args:\n",
                "        x: a string representing the first sequence\n",
                "        y: a string representing the second sequence\n",
                "        submatrix: a substitution matrix\n",
                "        g: the gap existence score for non-homopolymer gaps\n",
                "        h: the gap existence score for homopolymer gaps\n",
                "        s: the space score for non-homopolymer gaps\n",
                "        t: the space score for homopolymer gaps\n",
                "    Returns:\n",
                "        A tuple, (score, alignment), where score is a numeric value giving the score of the\n",
                "        alignment and alignment is a list of two strings\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-b927d1764be1b281",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 2: A strange optimizing ribosome (20 POINTS)\n",
                "\n",
                "Recall that the ribosomes are the molecular machines (which are actually complexes of proteins and RNAs) that perform the translation of RNA to protein.  Usually, the ribosome begins translation at a start codon and then translates each immediately following codon in turn until a stop codon is reached.  Suppose that we encounter a very strange ribosome that occassionally skips over some number of bases after it translates a codon, resulting in one or more codons being skipped and\/or the reading frame changing.  For example in the scenario below, the strange ribosome skips the first two bases, then translates a codon, then skips one base, and then translates two consecutive codons.\n",
                "\n",
                "![strange ribosome](hw2_strange_ribosome.png)\n",
                "\n",
                "It appears that this strange ribosome is optimizing an objective function that rewards the translation of specific amino acids and penalizes for the number of bases that are skipped. Letting $\\mathbf{x}$ denote an RNA sequence and $\\mathbf{C}$ denote the set of start positions of the (non-overlapping) codons that are translated by the ribosome, the ribosome *maximizes* the objective function:\n",
                "\n",
                "$f(\\mathbf{C}) = s(|\\mathbf{x}| - 3|\\mathbf{C}|) + \\sum_{i \\in C} r\\left(translation\\left(x_i x_{i+1} x_{i+2}\\right)\\right)$\n",
                "\n",
                "where $r(a)$ denotes the reward that the ribosome gets for translating amino acid $a$, $s$ denotes the (negative) penalty that the ribosome incurs for skipping a base (which includes bases skipped at the beginning of the RNA, at the end of the RNA, and in between translated codons), and $translation$ denotes the function that translates a codon sequence to its corresponding amino acid.\n",
                "\n",
                "For example, with $s = -1$ and \n",
                "\n",
                "$r(a) = \\left\\{\\begin{array}{ll}\n",
                "6 & \\textrm{if } a < \\textrm{'M'} \\\\\n",
                "1 & \\textrm{if } a \\geq \\textrm{'M'} \\\\\n",
                "-\\infty & \\textrm{if } a \\textrm{ is the stop codon symbol}\n",
                "\\end{array}\\right.$\n",
                "\n",
                "where $a < \\textrm{'M'}$ indicates that the character for the amino acid $a$ is lexigraphically smaller than 'M', \n",
                "the value of the objective function for the set of codons translated by strange ribosome in the figure above is \n",
                "\n",
                "$\\begin{eqnarray}\n",
                "f(\\{3, 7, 10\\}) & = & s(12 - 3 \\cdot 3) + r(\\textrm{G}) + r(\\textrm{H}) + r(\\textrm{R}) \\\\\n",
                "& = & -3 + 6 + 6 + 1 \\\\\n",
                "& = & 10\n",
                "\\end{eqnarray}\n",
                "$\n",
                "\n",
                "In this problem, we will figure out how to compute the set of codons that is translated by the strange ribosome for any sequence, $\\mathbf{x}$, by using a dynamic programming algorithm to maximize the objective function.  Let a subproblem $F[i]$ denote the score of the best set of codon start positions for the prefix of length $i$ of $\\mathbf{x}$.\n",
                "\n",
                "**(A)** Give initialization equations for $F[0]$, $F[1]$, and $F[2]$.\n",
                "\n",
                "**(B)** To devise a recurrence equation for $F[i]$, for $i > 2$, we consider two possible cases for the best set of codon start positions for the $i$th prefix: either (i) the last three bases of the prefix form a codon or (ii) the last three bases do not form a codon (which implies that the last base of the prefix is not in a codon).  Fill in the blanks in the recurrence equation below.\n",
                "\n",
                "$F[i] = \\max\\left\\{\n",
                "\\begin{array}{l}\n",
                "F[i - 3] + \\underline{\\hspace{2in}} \\\\\n",
                "F[i - 1] + \\underline{\\hspace{2in}}\n",
                "\\end{array}\n",
                "\\right.$\n",
                "\n",
                "**(C)** What amino acid sequence would this strange ribosome produce for the RNA sequence `CGCCCAAGGUUG` with the example scoring scheme given above?  Show your work."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "problem2",
                    "locked": false,
                    "points": 20,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "source": [
                "###\n",
                "### Your written solutions to PROBLEM 2\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-f2373176973d3b7b",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 3: Classifying bacterial DNA fragments (20 POINTS)\n",
                "\n",
                "Suppose that we obtain some microbial DNA from a sample of water taken from Lake Mendota and sequence a number of short DNA fragments within this sample.  A common task in such experiments is to identify the microbial species to which each DNA fragment belongs.  Suppose that one (very short) sequence obtained in this experiment is:\n",
                "\n",
                "`GTGATCTTCA`\n",
                "\n",
                "which has 2 `A`s, 2 `C`s, 2 `G`s, and 4 `T`s.  \n",
                "\n",
                "Further, suppose that there are only three microbial species, A, B, and C, to which the sequence could belong.  Unfortunately, we do not have the genome sequences for these species, but we do know the overall frequencies of bases within the genomes of the three species, which are given in the table below:\n",
                "\n",
                "Species | % A | % C | % G | % T\n",
                "--- | --- | ---\n",
                "A | 25 | 25 | 25 | 25\n",
                "B | 30 | 20 | 20 | 30\n",
                "C | 10 | 40 | 40 | 10\n",
                "\n",
                "In addition, we know that the relative abundances of the three species within Lake Mendota are:\n",
                "\n",
                "Species | % abundance\n",
                "--- | --- \n",
                "A | 50\n",
                "B | 25\n",
                "C | 25\n",
                "\n",
                "Let $X$ be a random variable representing the vector of *counts* of the bases within a sequence fragment, and let $Y$ be a random variable representing the species to which that sequence belongs.  Let $x = [2, 2, 2, 4]$ be the vector of counts of `A`, `C`, `G`, and `T`, respectively, in the sequence above.  Compute the following and *show your work*:\n",
                "\n",
                "**(A)** $P(X = x\\ |\\ Y = A)$\n",
                "\n",
                "**(B)** $P(X = x\\ |\\ Y = B)$\n",
                "\n",
                "**(C)** $P(X = x\\ |\\ Y = C)$\n",
                "\n",
                "**(D)** $P(X = x)$\n",
                "\n",
                "**(E)** $P(Y = A\\ |\\ X = x)$"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "problem3",
                    "locked": false,
                    "points": 20,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "source": [
                "###\n",
                "### Your written solutions to PROBLEM 3\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-74ce6b4272ffdd9b",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Tests for PROBLEM 1"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-34048ad0320bda93",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Substitution matrices for testing"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-83725a31e0f0f936",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "DNA = \"ACGT\"\n",
                "BASE_TYPE = {'A': 'purine',\n",
                "             'G': 'purine',\n",
                "             'C': 'pyrimidine',\n",
                "             'T': 'pyrimidine'}\n",
                "\n",
                "# A simple match=+10 and mismatch=-10 substitution matrix\n",
                "basic_submatrix = {(a, b): 10 if a == b else -10 for a in DNA for b in DNA}\n",
                "\n",
                "# A matrix that penalizes transitions less than transversions\n",
                "tt_submatrix = {(a, b): 10 if a == b else (-10 if BASE_TYPE[a] != BASE_TYPE[b] else -5)\n",
                "                for a in DNA for b in DNA}\n",
                "\n",
                "# A utility function for displaying a substitution matrix\n",
                "def print_dict_matrix(m, width=5):\n",
                "    \"\"\"Prints the dict-based matrix m with the specified width for each column.\"\"\"\n",
                "    def print_row(fields):\n",
                "        print(\"\".join([\"{:>{width}}\".format(field, width=width) for field in fields]))\n",
                "    labels = sorted({c for pair in m for c in pair})\n",
                "    print_row([\"\"] + labels)\n",
                "    for a in labels:\n",
                "        print_row([a] + [m.get((a, b), \"\") for b in labels])"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-17c18dc8fb0fdcfe",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Functions for scoring an a global alignment with affine and context-sensitive gap scores"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-cca52f4490120d70",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def score_alignment(alignment, submatrix, g, h, s, t):\n",
                "    \"\"\"Returns the score of the alignment given a substitution matrix and gap scores.\"\"\"\n",
                "    aligned_pair_scores = [submatrix[pair] for pair in zip(*alignment) if '-' not in pair]\n",
                "    insert_scores = list(insertion_scores(alignment, g, h, s, t))\n",
                "    return sum(aligned_pair_scores) + sum(insert_scores)\n",
                "    \n",
                "def insertion_scores(alignment, g, h, s, t):\n",
                "    \"\"\"Returns a list of the scores of the insertions within the alignment.\"\"\"\n",
                "    for i, start, end in insertions(alignment):\n",
                "        hp_base = alignment[i][start - 1: start]\n",
                "        length = end - start\n",
                "        if alignment[i][start: end] == hp_base * length:\n",
                "            yield h + length * t\n",
                "        else:\n",
                "            yield g + length * s\n",
                "\n",
                "def insertions(alignment):\n",
                "    \"\"\"Returns an iterator over the insertions within the alignment.\n",
                "    Each element of the iterator has the form (seq_index, start, end)\"\"\"\n",
                "    for i, seq in enumerate(alignment):\n",
                "        for gap in gaps(seq):\n",
                "            yield (1 - i, gap[0], gap[1])\n",
                "\n",
                "import re\n",
                "def gaps(s):\n",
                "    \"\"\"Returns an iterator over the gaps of the string s.\n",
                "    Each element of the iterator has the form (start, end)\"\"\"\n",
                "    gap_pattern = re.compile('-+')\n",
                "    for match in gap_pattern.finditer(s):\n",
                "        yield (match.start(), match.end())"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-4f8eeeadd962abe9",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "###  Test cases and testing functions"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-22b0304ca5404964",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "test_case_inputs = {\n",
                "    'small_1':  (\"AGTA\", \"AGA\",  basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_2':  (\"AGT\",  \"AGT\",  basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_3':  (\"AGT\",  \"G\",    basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_4':  (\"G\",    \"AGT\",  basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_5':  (\"AGT\",  \"\",     basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_6':  (\"\",     \"AGT\",  basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_7':  (\"A\",    \"\",     basic_submatrix, -20, -5, -10, -1),\n",
                "    'small_8':  (\"CCAG\", \"AGTT\", tt_submatrix,    -20, -5, -10, -1),\n",
                "    'small_9':  (\"CTC\",  \"TGCT\", basic_submatrix,   0,  0, -10, -1),\n",
                "    'small_10': ('TTTA', 'TAAA', basic_submatrix, -20, -5, -10, -1),\n",
                "    'large_1':  (\"GAATCGGTCCATTTCACAGTTCACACGCGAGGATTTTTGACCTCTGCACAAGCGCGCACATCG\",\n",
                "                 \"TTTTTTATACCCTCTCTACGCCCGGAATCCACAGAACGACCTACCATCTGGGTTACATCGCGCAGATTTAA\",\n",
                "                 basic_submatrix, -20, -5, -10, -1)\n",
                "}\n",
                "\n",
                "test_case_correct_outputs = {\n",
                "    'small_1': (0,   ['AGTA', \n",
                "                      'AG-A']),\n",
                "    'small_2': (30,  ['AGT', \n",
                "                      'AGT']),\n",
                "    'small_3': (-50, ['AGT', \n",
                "                      'G--']),\n",
                "    'small_4': (-50, ['--G', \n",
                "                      'AGT']),\n",
                "    'small_5': (-50, ['AGT', \n",
                "                      '---']),\n",
                "    'small_6': (-50, ['---', \n",
                "                      'AGT']),\n",
                "    'small_7': (-30, ['A', \n",
                "                      '-']),\n",
                "    'small_8': (-37, ['CCAG-', \n",
                "                      'A-GTT']),\n",
                "    'small_9': (-10, ['--CTC', \n",
                "                      'TGCT-']),\n",
                "    'small_10': (6,  ['TTTA--',\n",
                "                      'T--AAA']),\n",
                "    'large_1': (18, ['G-----AATC--GGTCCATTTCACAGT-TC-ACACG-CGAGGATTTTTGACCTCTG--C-ACAAGCGCGCACAT--CG',\n",
                "                     'TTTTTTATACCCTCTCTACG-CCCGGAATCCACAGAACGAC-CT----ACCATCTGGGTTACA-TCGCGCAGATTTAA'])\n",
                "}\n",
                "\n",
                "import numbers\n",
                "def check_valid_alignment_result(result, x, y):\n",
                "    \"\"\"Checks that the alignment result is valid for sequences x and y.\"\"\"\n",
                "    assert isinstance(result, tuple), \"Output is not a tuple\"\n",
                "    assert len(result) == 2, \"Output does not have exactly two elements\"\n",
                "    score, alignment = result\n",
                "    assert isinstance(alignment, list), \"Alignment is not a list\"\n",
                "    assert isinstance(score, numbers.Number), \"Score is not a number\"\n",
                "    assert len(alignment) == 2, \"Alignment does not have exactly two elements\"\n",
                "    assert all(isinstance(element, str) for element in alignment), \"Alignment elements are not strings\"\n",
                "    assert len(alignment[0]) == len(alignment[1]), \"Alignment strings do not have the same length\"\n",
                "    assert alignment[0].replace('-', '') == x, \"First string of alignment is not x\"\n",
                "    assert alignment[1].replace('-', '') == y, \"Second string of alignment is not y\"\n",
                "\n",
                "def check_valid_alignment_score(result, submatrix, g, h, s, t):\n",
                "    \"\"\"Checks that the computed score of the alignment is equal to the score given in the result.\"\"\"\n",
                "    score, alignment = result\n",
                "    computed_score = score_alignment(alignment, submatrix, g, h, s, t)\n",
                "    assert computed_score == score, \"Computed score ({}) does not equal the returned score ({})\".format(computed_score,\n",
                "                                                                                                        score)\n",
                "def check_test_case(case_name, test_name=None, valid_result=True, valid_score=True, correct_alignment=True):\n",
                "    inputs = test_case_inputs[case_name]\n",
                "    correct_output = test_case_correct_outputs[case_name]\n",
                "    result = align_global_affine_hp_gaps(*inputs)\n",
                "    if valid_result:\n",
                "        check_valid_alignment_result(result, *inputs[:2])\n",
                "    if valid_score:\n",
                "        check_valid_alignment_score(result, *inputs[2:])\n",
                "    if correct_alignment:\n",
                "        error_template = (\"Incorrect output\\n\"\n",
                "                          \"returned (score = {0[0]}):\\n\"\n",
                "                          \"{0[1][0]}\\n\"\n",
                "                          \"{0[1][1]}\\n\"\n",
                "                          \"correct (score = {1[0]}):\\n\"\n",
                "                          \"{1[1][0]}\\n\"\n",
                "                          \"{1[1][1]}\")\n",
                "        assert result == correct_output, error_template.format(result,\n",
                "                                                               correct_output)\n",
                "    print(\"SUCCESS:\", test_name if test_name else case_name, \"passed!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-d0b3479a99a9317a",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Visible tests"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_1_valid",
                    "locked": true,
                    "points": 10,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_1 valid output\n",
                "check_test_case(\"small_1\", test_name=\"small_1 valid output\", valid_score=False, correct_alignment=False)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_1_valid_score",
                    "locked": true,
                    "points": 10,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_1 valid score\n",
                "check_test_case(\"small_1\", test_name=\"small_1 valid score\", correct_alignment=False)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_2",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_2\n",
                "check_test_case(\"small_2\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_3",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_3\n",
                "check_test_case(\"small_3\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_4",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_4\n",
                "check_test_case(\"small_4\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_5",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_5\n",
                "check_test_case(\"small_5\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_6",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_6\n",
                "check_test_case(\"small_6\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_7",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_7\n",
                "check_test_case(\"small_7\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_8",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_8\n",
                "check_test_case(\"small_8\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_9",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_9\n",
                "check_test_case(\"small_9\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "small_10",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: small_10\n",
                "check_test_case(\"small_10\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_1",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_1\n",
                "check_test_case(\"large_1\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-9d707db41eb05083",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Hidden tests"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_2",
                    "locked": true,
                    "points": 2,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_2\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_3",
                    "locked": true,
                    "points": 2,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_3\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_4",
                    "locked": true,
                    "points": 2,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_4\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_5",
                    "locked": true,
                    "points": 2,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_5\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "large_6",
                    "locked": true,
                    "points": 2,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# TEST: large_6\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n"
            ]
        }
    ],
    "metadata": {
        "celltoolbar": "Create Assignment",
        "kernelspec": {
            "display_name": "Python 3 [3.6]",
            "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.4"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}