Day 12 - Multiple Alignment/Starter/notebook.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-9889a0a8d7a46373",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "# Day 12 notebook\n",
                "\n",
                "The objectives of this notebook are to practice\n",
                "\n",
                "* computing scores of multiple alignment columns\n",
                "* progressive multiple alignment\n",
                "* iterative refinement"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-6c93ada6842a576c",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Modules used in this notebook"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-12d427fcb5be2ff7",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "import matplotlib.pyplot as plt # for plotting\n",
                "import fasta                    # For reading in an alignment\n",
                "import submatrix                # For reading in a substitution matrix"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-f7634a8cb5079638",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Multiple alignment of the CFTR protein"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-99a2a388e0eb6b4b",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "In this activity we will analyze a multiple alignment of the CFTR (Cystic fibrosis transmembrane conductance regulator) protein from a variety of vertebrate species.  We will implement the sum-of-pairs column scoring approaches and apply it to analyze the columns of this multiple alignment."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-42d6bb8c1577c3ec",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "We will first read in the CFTR multiple alignment."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-92d698980e71d063",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "cftr_aligned_sequences = fasta.read_sequences_from_fasta_file(\"cftr_alignment.fasta\")\n",
                "cftr_names, cftr_alignment = zip(*cftr_aligned_sequences)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-4c29c6fe856f2534",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Here are the species of the CFTR proteins in this alignment:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-349c7cfc2573cf9c",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "print(*cftr_names, sep=\"\\n\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-1d0dbf8c0db2c551",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "And here is a slice of this alignment in relatively interesting region:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-9e19dce890aa159f",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def slice_alignment(alignment, start, length):\n",
                "    return [s[start: start + length] for s in alignment]\n",
                "\n",
                "slice_alignment(cftr_alignment, 660, 100)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-71c59ee3f00b8d3f",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Test sequences for multiple alignment algorithms"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-23402d036daca5fd",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Below are three sets of homologous amino acid sequences, taken from the CFTR alignment. Each sequence is a short substring of a CFTR protein sequence.  These sets will be used in testing the multiple sequence alignment algorithms that you implement below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-bdf0376e7c7f73a5",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "tiny_set = {\n",
                "    'Zebrafish': 'SETPKPD',\n",
                "    'Human': 'TEGKPTKS',\n",
                "    'Mouse': 'TEESMYTQI'\n",
                "}\n",
                "tiny_seqs = sorted(tiny_set.values())\n",
                "\n",
                "small_set = {\n",
                "    'Guinea pig': 'EEGAPVKS',\n",
                "    'Zebrafish': 'SETPKPD',\n",
                "    'Human': 'TEGKPTKS',\n",
                "    'Mouse': 'TEESMYTQI',\n",
                "    'Spiny dogfish': 'PEGSETK'\n",
                "}\n",
                "small_seqs = sorted(small_set.values())\n",
                "\n",
                "large_set = {\n",
                "    'Guinea pig': 'FSIVPKTPLQISGIEEDSDDPVERRL',\n",
                "    'Zebrafish': 'ASFIQVPEEEVRRTLPDRKF',\n",
                "    'Human': 'FSIVQKTPLQMNGIEEDSDEPLERRL',\n",
                "    'Mouse': 'ISIVQKTPLCIDGESDDLQEKRL',\n",
                "    'Small-eared galago': 'FSVVQKTPLPMNGIDEEDSEEPVERRL',\n",
                "    'Rat': 'ISIVQKTPLSIEGESDDLQERRL',\n",
                "    'Spiny dogfish': 'FSLVQTAMSYPQTNGMEDATSEPGERHF'\n",
                "}\n",
                "large_seqs = sorted(large_set.values())"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-f86769537a07fcf9",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Sum of pairs scoring"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-a79f064cc5ccc879",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### PROBLEM 1: Computing the sum of pairs score of a column (1 POINT)\n",
                "\n",
                "Write a function `sum_of_pairs_column_score` below, which takes as input a column from a multiple alignment (represented as a string) and a substitution matrix, and computes the sum of pairs score of that column.  Recall that the sum of pairs score for a column is defined as:\n",
                "\n",
                "$S_{sum\\_of\\_pairs}(m_i) = \\sum_{k < l} s(m_i^k, m_i^l)$\n",
                "\n",
                "where $s(m_i^k, m_i^l)$ is the substitution score of the character in the $k$th row with the character in the $l$th row.\n",
                "\n",
                "For simplicity, we will simply treat a space ('-') like any other character in a column, and will add the space character to the substitution matrix, with the score of aligning a character to a space being the space score of a linear gap penalty function."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-3c51c13c3843e4b8",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "First we will set up the substitution matrix, which will be the BLOSUM62 matrix with a space score of -2 added to it."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-59a701d7a78ad043",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def submatrix_with_spaces(matrix, space_score):\n",
                "    \"\"\"Returns an extended substitution matrix representing a linear gap penalty function.\"\"\"\n",
                "    extended_matrix = matrix.copy()\n",
                "    chars = {a for (a, b) in matrix}\n",
                "    for a in chars:\n",
                "        extended_matrix[a, '-'] = space_score\n",
                "        extended_matrix['-', a] = space_score\n",
                "    extended_matrix['-', '-'] = 0\n",
                "    return extended_matrix\n",
                "\n",
                "blosum62 = submatrix.read_substitution_matrix(\"BLOSUM62.txt\")\n",
                "space_score = -2\n",
                "blosum62_with_spaces = submatrix_with_spaces(blosum62, space_score)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-42a5b19afdc9c536",
                    "locked": false,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "outputs": [],
            "source": [
                "def sum_of_pairs_column_score(column, matrix):\n",
                "    \"\"\"Returns the sum of pairs score of an alignment column.\n",
                "    \n",
                "    Args:\n",
                "        column: a column of a multiple alignment, represented as a string with the ith character\n",
                "                specifying the character from the ith sequence in the alignment.\n",
                "        matrix: a substitution matrix, represented as a dictionary\n",
                "    Returns:\n",
                "        The sum of pairs score of the column.\"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "sum_of_pairs_column_score",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for sum_of_pairs_column_score\n",
                "assert sum_of_pairs_column_score('SQSQ', blosum62_with_spaces) == 9\n",
                "assert sum_of_pairs_column_score('WWWWWW', blosum62_with_spaces) == 165\n",
                "assert sum_of_pairs_column_score('ACDEFG', blosum62_with_spaces) == -27\n",
                "assert sum_of_pairs_column_score('GEGEKRK', blosum62_with_spaces) == 4\n",
                "assert sum_of_pairs_column_score('--S-', blosum62_with_spaces) == -6\n",
                "print(\"SUCCESS: sum_of_pairs_column_score passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-b094675700e01d28",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Sum of pairs scores of columns in the CFTR alignment\n",
                "We will now use your `sum_of_pairs_column_score` to compute and examine the sum of pairs scores of the columns within the CFTR alignment."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-671ed1d3abd6f627",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def transpose_alignment(alignment):\n",
                "    \"\"\"Returns a column-based alignment from a row-based alignment or vice versa\"\"\"\n",
                "    return list(map(''.join, zip(*alignment)))\n",
                "\n",
                "cftr_alignment_columns = list(zip(*cftr_alignment))\n",
                "sum_of_pairs_scores = [sum_of_pairs_column_score(column, blosum62_with_spaces) for column in cftr_alignment_columns]"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-ec4800f42e2b600f",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "plt.figure(figsize=(15, 3), dpi=150)\n",
                "plt.plot(sum_of_pairs_scores, \".\")\n",
                "plt.xlabel(\"alignment column\")\n",
                "plt.ylabel(\"sum of pairs score\")\n",
                "plt.title(\"Sum of pairs scores of columns in the CFTR alignment\")\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-5bc8b478c2f679f0",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "What does this plot tell you regarding which regions of the CFTR proteins are more or less conserved?"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "###\n",
                "### *Your thoughts here*\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Multiple alignment algorithms"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-2f9553f4b1341f75",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### A core global pairwise alignment function\n",
                "Below is a function, `align_global`, that implements standard global pairwise alignment, but that generically supports aligning two alignments to each other.  You will use this function to implement both progressive alignment and iterative refinement."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-0c8b21160a5de420",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def align_global(x, y, substitution_matrix):\n",
                "    \"\"\"Computes a global pairwise alignment of sequences\/alignments x and y\n",
                "    \n",
                "    Uses a linear gap scoring function with the space score encoded in the given substitution matrix.\n",
                "    In the case of multiple optimal alignments, this function returns a \"middle-road\" \n",
                "    alignment (prefers match\/mismatch first).\n",
                "    \n",
                "    Args:\n",
                "        x: an alignment\/sequence represented as a list of strings. A single \n",
                "           unaligned sequence should be represented by a list with a single element (the sequence).\n",
                "        y: an alignment\/sequence represented as a list of strings. A single \n",
                "           unaligned sequence should be represented by a list with a single element (the sequence).\n",
                "        substitution_matrix: the substitution matrix, represented as a dictionary \n",
                "    Returns:\n",
                "        An alignment, represented as a list of aligned strings of the same length.\n",
                "    \"\"\"\n",
                "    # create lists of the columns of the two alignments\n",
                "    xt, yt = transpose_alignment(x), transpose_alignment(y)\n",
                "    num_rows, num_cols = len(xt) + 1, len(yt) + 1\n",
                "    x_all_spaces, y_all_spaces = '-' * len(x), '-' * len(y)\n",
                "    m = matrix(num_rows, num_cols) # score matrix\n",
                "    t = matrix(num_rows, num_cols) # traceback pointer matrix\n",
                "\n",
                "    # constants representing the traceback pointers, with higher\n",
                "    # values representing higher priority for traceback in case of ties\n",
                "    M, IX, IY = 3, 2, 1\n",
                "    \n",
                "    # A function that scores the alignment of two columns\n",
                "    def score_column_pair(x_column, y_column):\n",
                "        return sum(substitution_matrix[a, b] for a in x_column for b in y_column)\n",
                "\n",
                "    # initialization\n",
                "    m[0][0] = 0\n",
                "    for i in range(1, num_rows):\n",
                "        m[i][0] = m[i - 1][0] + score_column_pair(xt[i - 1], y_all_spaces)\n",
                "        t[i][0] = IX\n",
                "    for j in range(1, num_cols):\n",
                "        m[0][j] = m[0][j - 1] + score_column_pair(x_all_spaces, yt[j - 1])\n",
                "        t[0][j] = IY\n",
                "  \n",
                "    # main fill\n",
                "    for i in range(1, num_rows):\n",
                "        for j in range(1, num_cols):\n",
                "            m[i][j], t[i][j] = max((m[i - 1][j - 1] + score_column_pair(xt[i - 1], yt[j - 1]),     M),\n",
                "                                   (m[i - 1][j]     + score_column_pair(xt[i - 1], y_all_spaces), IX),\n",
                "                                   (m[i][j - 1]     + score_column_pair(x_all_spaces, yt[j - 1]), IY))\n",
                "\n",
                "    # traceback\n",
                "    columns = []\n",
                "    i, j = num_rows - 1, num_cols - 1\n",
                "    while i > 0 or j > 0:\n",
                "        traceback = t[i][j]\n",
                "        if traceback == M:\n",
                "            i -= 1\n",
                "            j -= 1            \n",
                "            columns.append(xt[i] + yt[j])\n",
                "        elif traceback == IX:\n",
                "            i -= 1\n",
                "            columns.append(xt[i] + y_all_spaces)\n",
                "        elif traceback == IY:\n",
                "            j -= 1\n",
                "            columns.append(x_all_spaces + yt[j])\n",
                "        else:\n",
                "            assert False\n",
                "    \n",
                "    return transpose_alignment(reversed(columns))\n",
                "\n",
                "def matrix(num_rows, num_cols, initial_value=None):\n",
                "    \"\"\"Returns a matrix (a list of rows, each of which is a list) \n",
                "    with num_rows and num_cols and with initial_value in each entry\"\"\"\n",
                "    return [[initial_value] * num_cols for i in range(num_rows)]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-28bc551441691e68",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 2: Progressive alignment (1 POINT)\n",
                "Write a function that computes a multiple alignment of a set of sequences using the progressive alignment strategy.  Your function should take in an \"alignment order\" which specifies the ordering of the pairwise alignments to perform.  For example, if there are four sequences, the alignment order \n",
                "\n",
                "`[((1,),   (3,)),\n",
                " ((0,),   (2,)),\n",
                " ((1, 3), (0, 2))]`\n",
                " \n",
                " would specify the following sequence of pairwise alignments:\n",
                " \n",
                " 1. align sequence 1 to sequence 3\n",
                " 2. align sequence 0 to sequence 2\n",
                " 3. align the partial alignment of sequences 1 and 3 to the partial alignment of sequences 0 and 2\n",
                " \n",
                " "
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-394b085441e5b9e1",
                    "locked": false,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "outputs": [],
            "source": [
                "def align_progressive(sequences, alignment_order, substitution_matrix):\n",
                "    \"\"\"Computes a multiple alignment with the progressive alignment strategy\n",
                "    and a given ordering of pairwise alignments to perform.\n",
                "    \n",
                "    Args:\n",
                "        sequences: a list of strings\n",
                "        alignment_order: a list of two-element tuples, with each element\n",
                "                         of a tuple being a tuple of integer indices, \n",
                "                         specifying a partial alignment of a subset of\n",
                "                         the sequences (indicated by the indices). Each element\n",
                "                         of the alignment_order list specifies the next pair of\n",
                "                         partial alignments to align with each other.\n",
                "        substitution_matrix: a substitution matrix, represented as a dictionary\n",
                "    Returns:\n",
                "        A multiple alignment, represented by a list of aligned strings.\n",
                "    \"\"\"\n",
                "    # we will store the partial alignments that we construct in a dictionary with tuples\n",
                "    # of (sorted) indices of the sequences in each partial alignment as keys\n",
                "    # we start with each sequence as its own partial alignment\n",
                "    partial_alignments = {(i,): [seq] for i, seq in enumerate(sequences)}\n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "    \n",
                "    # by this point there should be an entry in partial_alignments for the tuple\n",
                "    # containing the indices of all sequences\n",
                "    full_set = tuple(range(len(sequences)))\n",
                "    full_alignment = partial_alignments[full_set]\n",
                "    # return the full alignment with the sequences in the same order as in the input\n",
                "    return sorted(full_alignment, key=lambda s: sequences.index(s.replace('-', '')))"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "align_progressive",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for align_progressive\n",
                "tiny_order1 = [((0,), (1,)), \n",
                "               ((0, 1), (2,))]\n",
                "tiny_result1 = ['S-E-TPKPD-', \n",
                "                'TEE-SMYTQI', \n",
                "                'T-EGKP-TKS']\n",
                "assert align_progressive(tiny_seqs, tiny_order1, blosum62_with_spaces) == tiny_result1\n",
                "\n",
                "tiny_order2 = [((0,), (2,)), \n",
                "               ((0, 2), (1,))]\n",
                "tiny_result2 = ['SE-TP--KPD', \n",
                "                'TEESMYTQ-I', \n",
                "                'TEGKP-TK-S']\n",
                "assert align_progressive(tiny_seqs, tiny_order2, blosum62_with_spaces) == tiny_result2\n",
                "\n",
                "tiny_order3 = [((1,), (2,)), \n",
                "               ((1, 2), (0,))]\n",
                "tiny_result3 = ['S-ET-P-KPD', \n",
                "                'TEESMYTQ-I', \n",
                "                'T-EGKPTK-S']\n",
                "assert align_progressive(tiny_seqs, tiny_order3, blosum62_with_spaces) == tiny_result3\n",
                "\n",
                "small_order1 = [((0,), (1,)),\n",
                "                ((2,), (3,)),\n",
                "                ((2, 3), (4,)),\n",
                "                ((0, 1), (2, 3, 4))]\n",
                "    \n",
                "small_result1 = ['E-EGAP-VKS', \n",
                "                 'P-EGSE-TK-', \n",
                "                 'S-E-TPKPD-', \n",
                "                 'TEE-SMYTQI', \n",
                "                 'T-EGKP-TKS']\n",
                "assert align_progressive(small_seqs, small_order1, blosum62_with_spaces) == small_result1\n",
                "print(\"SUCCESS: align_progressive passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-9030879f44e08c85",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 3: Iterative refinement (1 POINT)\n",
                "We will next write a function for performing iterative refinement on a multiple alignment.  The main function, `refine_alignment` is written for you below.  You are to write `realign_first_sequence` which takes as input an alignment and realigns the first sequence in the alignment back to the remainder of the alignment.  You will need to make sure that you\n",
                "\n",
                "1. Remove space characters from the first sequence before aligning it back\n",
                "2. Remove columns with all space characters in the partial alignment that results after removing the first sequence."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-3ccfb3f70087db18",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def refine_alignment(alignment, substitution_matrix):\n",
                "    \"\"\"Performs iterative refinement on a multiple alignment.\n",
                "    \n",
                "    Args:\n",
                "        alignment: a multiple alignment, represented by a list of aligned strings.\n",
                "        substitution_matrix: a substitution matrix, represented by a dictionary.\n",
                "    Returns:\n",
                "        A multiple alignment, represented by a list of aligned strings.\"\"\"\n",
                "    previous_alignment = None\n",
                "    current_alignment = alignment\n",
                "    while current_alignment != previous_alignment:\n",
                "        previous_alignment = current_alignment\n",
                "        for i in range(len(current_alignment)):\n",
                "            # realign the first sequence\n",
                "            current_alignment = realign_first_sequence(current_alignment, substitution_matrix)\n",
                "            # move first sequence to the bottom of the alignment\n",
                "            current_alignment = current_alignment[1:] + current_alignment[:1]\n",
                "    return current_alignment"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-a3c70d0b3b12f20a",
                    "locked": false,
                    "schema_version": 1,
                    "solution": true
                }
            },
            "outputs": [],
            "source": [
                "def realign_first_sequence(alignment, substitution_matrix):\n",
                "    \"\"\"Returns a new alignment that results from realigning the first\n",
                "    sequence of the given alignment back to the remaining alignment.\n",
                "    \n",
                "    Args:\n",
                "        alignment: a multiple alignment, represented by a list of aligned strings.\n",
                "        substitution_matrix: a substitution matrix, represented by a dictionary.\n",
                "    Returns:\n",
                "        A multiple alignment, represented by a list of aligned strings.\"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "realign_first_sequence",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for realign_first_sequence\n",
                "tiny_result1 = ['S-E-TPKPD-', \n",
                "                'TEE-SMYTQI', \n",
                "                'T-EGKP-TKS']\n",
                "tiny_refine1 = ['S-E-TP--KPD', \n",
                "                'TEE-SMYTQ-I', \n",
                "                'T-EGKP-TK-S']\n",
                "assert realign_first_sequence(tiny_result1, blosum62_with_spaces) == tiny_refine1\n",
                "\n",
                "tiny_result2 = ['SE-TP--KPD', \n",
                "                'TEESMYTQ-I', \n",
                "                'TEGKP-TK-S']\n",
                "tiny_refine2 = ['SE-TP--KPD', \n",
                "                'TEESMYTQ-I', \n",
                "                'TEGKP-TK-S']\n",
                "assert realign_first_sequence(tiny_result2, blosum62_with_spaces) == tiny_refine2\n",
                "\n",
                "# for the last test we will run refine_alignment rather than realign_first_sequence\n",
                "small_result1 = ['E-EGAP-VKS', \n",
                "                 'P-EGSE-TK-', \n",
                "                 'S-E-TPKPD-', \n",
                "                 'TEE-SMYTQI', \n",
                "                 'T-EGKP-TKS']\n",
                "\n",
                "small_refine1 = ['EEGAP-VK-S', \n",
                "                 'PEGSE-TK--', \n",
                "                 'SE-TP--KPD', \n",
                "                 'TEESMYTQ-I', \n",
                "                 'TEGKP-TK-S']\n",
                "assert refine_alignment(small_result1, blosum62_with_spaces) == small_refine1\n",
                "print(\"SUCCESS: realign_first_sequence passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-7992089c7e7c05da",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## CHALLENGE: Find the best multiple alignment of the \"large\" test set\n",
                "Using your progressive alignment and iterative refinment functions above, try your hand at finding the best scoring alignment of the \"large\" test set.  Use the function below to compare the scores of your alignments."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-37de1582f617061e",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def score_alignment_sum_of_pairs(alignment, substitution_matrix):\n",
                "    \"\"\"Returns the sum of pairs score for the given alignment using the given\n",
                "    substitution matrix, which must be extended for a linear gap penalty function.\"\"\"\n",
                "    return sum(sum_of_pairs_column_score(column, substitution_matrix) \n",
                "               for column in transpose_alignment(alignment))"
            ]
        }
    ],
    "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
}