HW 3 - Neighbor Joining, Parsimony, Branch & Bound/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 - HW3\n",
                "The objectives of this homework are to practice phylogenetic tree reconstruction algorithms.  Specifically, you will gain practice with the following techniques:\n",
                "\n",
                "* Neighbor-joining\n",
                "* branch and bound\n",
                "* unweighted and weighted parsimony\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: Neighbor-joining (20 points)\n",
                "\n",
                "Construct a tree from these distances using the neighbor-joining algorithm. Show your updated distance matrix after each merge and give branch lengths for the final tree.\n",
                "\n",
                "|       | A | B | C | D | E |\n",
                "|-------|---|---|---|---|---|\n",
                "| **A** |   | 7 | 13|  9| 10| \n",
                "| **B** |   |   |  8|  4|  5|\n",
                "| **C** |   |   |   |  6|  9|\n",
                "| **D** |   |   |   |   |  5|\n",
                "| **E** |   |   |   |   |  &nbsp; |\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "###\n",
                "### Your solution to Problem 1 here\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-1f115a62bb2dd441",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Trees for Problem 2 and 3\n",
                "![parsimony_trees](parsimony_trees.png)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-d388d0f01da01423",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 2: Unweighted parsimony (20 POINTS)\n",
                "\n",
                "Suppose we are given five DNA sequences (1, 2, 3, 4, and 5), each of which is one base long. The figure above gives two possible trees relating these five sequences.\n",
                "\n",
                "**(a)** For each of the two trees, find the minimal cost of the tree using Fitch\u2019s algorithm (unweighted parsimony). Show the intermediate computations of the algorithm.\n",
                "\n",
                "**(b)** For each of the two trees, determine an assignment of ancestral bases that achieves the minimal cost that you found in (a).\n",
                "\n",
                "**(c)** Which tree would be preferred when using unweighted parsimony?"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "###\n",
                "### Your solution to Problem 2 here\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-47a34506ea0189d7",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 3: Weighted parsimony (20 POINTS)\n",
                "\n",
                "(25points) Suppose we are given the same five sequences as in Problem 2 and the same two possible trees. Given the weighted parsimony costs given in the matrix below:\n",
                "\n",
                "![weighted_parsimony_weights](weighted_parsimony_weights.png)\n",
                "\n",
                "**(a)** For each of the two trees, find the minimal cost of the tree using the weighted parsimony algorithm. Show the intermediate computations of the algorithm.\n",
                "\n",
                "**(b)** For each of the two trees, determine an assignment of ancestral bases that achieves the minimal cost that you found in (a).\n",
                "\n",
                "**(c)** Which tree would be preferred when using weighted parsimony with the given costs?\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "###\n",
                "### Your solution to Problem 3 here\n",
                "###\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-8cc69ce21bfbf312",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 4: Branch and bound (40 POINTS)\n",
                "\n",
                "Implement the basic branch and bound algorithm (page 7 of the \"Searching through tree space\" lecture slides) for finding an optimal tree with an unweighted parsimony objective function.  Your implementation should be broken down into three functions, `best_tree_branch_and_bound`, `branch`, and `bound`, which you are to implement below.\n",
                "\n",
                "Some implementation specifications for each function:\n",
                "### `best_tree_branch_and_bound`\n",
                "\n",
                "* You should use the [Heap queue algorithm](https:\/\/docs.python.org\/3\/library\/heapq.html) module functions to efficiently maintain the queue (Q) \n",
                "* Q should be sorted by lower bound of the trees, with ties broken by lexicographical ordering of the trees' newick strings.\n",
                "* The algorithm should begin with an unrooted tree consisting of the first three names in `sequence_names`, and add leaves in the order in which they appear in the sequence_names list.\n",
                "\n",
                "### `branch`\n",
                "\n",
                "* This function should call the `add_leaf` function (provided) to grow an unrooted tree with the next leaf in all possible ways\n",
                "\n",
                "### `bound`\n",
                "\n",
                "* You will likely want to take advantage of your work in the day 16 notebook for this function\n",
                "* To convert an unrooted tree to a rooted tree (such that you can use your parsimony scoring code), it is recommended that you call the provided `root` function"
            ]
        },
        {
            "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 4\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",
                "import heapq\n",
                "import toytree\n",
                "\n",
                "def best_tree_branch_and_bound(alignment, sequence_names, branch, bound):\n",
                "    \"\"\"Computes an optimal (lowest scoring) tree using a branch and bound algorithm.\n",
                "    \n",
                "    Args:\n",
                "        alignment: a list of strings corresponding to the rows of a multiple alignment.\n",
                "        sequence_names: a list of the names of the sequences in the same order as the\n",
                "                        rows of the multiple alignment.\n",
                "        branch: a function that grows a partial tree in multiple ways\n",
                "        bound: a function that computes the lower bound of a partial tree\n",
                "    Returns:\n",
                "        A tuple (score, newick_string) where newick_string is a Newick formatted string\n",
                "        representing the optimal tree (unrooted) and score is its score.\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "    \n",
                "def branch(newick_tree, next_leaf_name):\n",
                "    \"\"\"Grows a partial unrooted tree by adding the next leaf in all possible ways.\n",
                "    \n",
                "    Args:\n",
                "        newick_tree: a partial unrooted tree as a Newick string\n",
                "        next_leaf_name: the name of the next leaf to add to the tree\n",
                "    Returns:\n",
                "        A list of Newick strings representing all possible ways in which to add the next leaf.\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "\n",
                "def bound(newick_tree, alignment, sequence_names):\n",
                "    \"\"\"Computes a lower bound for the unweighted parsimony score of a full tree that can be \n",
                "    grown from a given partial tree.\n",
                "    \n",
                "    Args:\n",
                "        newick_tree: a partial unrooted tree as a Newick string\n",
                "        alignment: a list of strings corresponding to the rows of a multiple alignment.\n",
                "        sequence_names: a list of the names of the sequences in the same order as the\n",
                "                        rows of the multiple alignment.\n",
                "    Returns:\n",
                "        The lower bound as an integer.\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "\n",
                "def add_leaf(newick_tree, leaf_name, edge_index):\n",
                "    \"\"\"Adds a new leaf on to the specified edge in an unrooted tree.\n",
                "    \n",
                "    Args:\n",
                "        newick_tree: a partial unrooted tree as a Newick string\n",
                "        leaf_name: the name of the next leaf to add to the tree\n",
                "        edge_index: the index (from 1 to the number of edges) of the edge on\n",
                "                    which to add the leaf.  Edges are ordered by the order\n",
                "                    in which their child node is encountered in a preorder\n",
                "                    traversal of the tree.\n",
                "    Returns:\n",
                "        A newick string representing the tree with the added leaf.\n",
                "    \"\"\"\n",
                "    new_tree = toytree.tree(newick_tree)\n",
                "    if len(new_tree.treenode.children) != 3:\n",
                "        raise ValueError(\"Tree does not look unrooted: \" +  newick_tree)\n",
                "    for i, node in enumerate(new_tree.treenode.traverse(\"preorder\")):\n",
                "        if i == edge_index:\n",
                "            break\n",
                "    parent = node.up\n",
                "    node.detach()\n",
                "    new_internal_node = parent.add_child()\n",
                "    new_internal_node.add_child(node)\n",
                "    new_internal_node.add_child(name=leaf_name)\n",
                "    return new_tree.treenode.write(format=9)\n",
                "\n",
                "def root(newick_tree):\n",
                "    \"\"\"Converts an unrooted tree into a rooted tree.\n",
                "    (useful for scoring an unrooted tree with parsimony).\n",
                "    \n",
                "    Args:\n",
                "        newick_tree: an unrooted tree as a Newick string\n",
                "    Returns:\n",
                "        A rooted version of the tree as a newick string.\n",
                "    \"\"\"\n",
                "    unrooted_tree = toytree.tree(newick_tree)\n",
                "    if len(unrooted_tree.treenode.children) != 3:\n",
                "        raise ValueError(\"Tree does not look unrooted: \" +  newick_tree)\n",
                "    unrooted_tree_root = unrooted_tree.treenode\n",
                "    first_child = unrooted_tree.treenode.children[0]\n",
                "    first_child.detach()\n",
                "    rooted_tree_root = toytree.TreeNode.TreeNode()\n",
                "    rooted_tree_root.add_child(first_child)\n",
                "    rooted_tree_root.add_child(unrooted_tree_root)\n",
                "    return rooted_tree_root.write(format=9)\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 4"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-34048ad0320bda93",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Data sets for testing"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-34048",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "import fasta\n",
                "def read_names_and_alignments_from_fasta(filename):\n",
                "    return zip(*fasta.read_sequences_from_fasta_file(filename))\n",
                "\n",
                "v3_sequence_names, v3_alignment = read_names_and_alignments_from_fasta(\"v3_alignment.fasta\")\n",
                "\n",
                "v3_big_sequence_names, v3_big_alignment = read_names_and_alignments_from_fasta(\"v3_big_alignment.fasta\")\n",
                "\n",
                "medium_num_seqs = 7\n",
                "v3_medium_alignment = v3_big_alignment[:medium_num_seqs]\n",
                "v3_medium_sequence_names = v3_big_sequence_names[:medium_num_seqs]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Tests"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch3",
                    "locked": true,
                    "points": 4,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch3\n",
                "assert sorted(branch('(D,PA,PB);', 'C1')) == ['(D,PA,(PB,C1));', \n",
                "                                              '(D,PB,(PA,C1));', \n",
                "                                              '(PA,PB,(D,C1));']\n",
                "print(\"SUCCESS: branch passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch4",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch4\n",
                "assert sorted(branch('(D,PA,(PB,C1));', 'C2')) == ['(D,(PB,C1),(PA,C2));',\n",
                "                                                   '(D,PA,((PB,C1),C2));',\n",
                "                                                   '(D,PA,(C1,(PB,C2)));',\n",
                "                                                   '(D,PA,(PB,(C1,C2)));',\n",
                "                                                   '(PA,(PB,C1),(D,C2));']\n",
                "print(\"SUCCESS: branch4 passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch5",
                    "locked": true,
                    "points": 3,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch5\n",
                "assert sorted(branch('(D,(PB,C1),(PA,C2));', 'C3')) == ['((PB,C1),(PA,C2),(D,C3));',\n",
                "                                                        '(D,(C1,(PB,C3)),(PA,C2));',\n",
                "                                                        '(D,(PA,C2),((PB,C1),C3));',\n",
                "                                                        '(D,(PB,(C1,C3)),(PA,C2));',\n",
                "                                                        '(D,(PB,C1),((PA,C2),C3));',\n",
                "                                                        '(D,(PB,C1),(C2,(PA,C3)));',\n",
                "                                                        '(D,(PB,C1),(PA,(C2,C3)));']\n",
                "print(\"SUCCESS: branch5 passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "bound_one_column",
                    "locked": true,
                    "points": 5,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for bound_one_column\n",
                "v3_alignment_column0 = [s[0] for s in v3_alignment]\n",
                "v3_alignment_column3 = [s[3] for s in v3_alignment]\n",
                "v3_alignment_column29 = [s[29] for s in v3_alignment]\n",
                "assert bound('(D,PA,PB);', v3_alignment_column0, v3_sequence_names) == 0\n",
                "assert bound('(D,PA,PB);', v3_alignment_column3, v3_sequence_names) == 1\n",
                "assert bound('((C1,C2),(PA,PB),D);', v3_alignment_column29, v3_sequence_names) == 1\n",
                "print(\"SUCCESS: bound_one_column passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "bound_alignment",
                    "locked": true,
                    "points": 5,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for bound_alignment\n",
                "assert bound('(D,PA,PB);', v3_alignment, v3_sequence_names) == 25\n",
                "assert bound('((D, C1),PA,PB);', v3_alignment, v3_sequence_names) == 48\n",
                "assert bound('((C1,PA),(C2,PB),D);', v3_alignment, v3_sequence_names) == 64\n",
                "print(\"SUCCESS: bound_alignment passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-0b40486d279ff93d",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "class CallCounter:\n",
                "    def __init__(self, func):\n",
                "        self._func = func\n",
                "        self._num_calls = 0\n",
                "    def __call__(self, *args, **kwds):\n",
                "        self._num_calls += 1\n",
                "        return self._func(*args, **kwds)\n",
                "    def num_calls(self):\n",
                "        return self._num_calls\n",
                "    def reset(self):\n",
                "        self._num_calls = 0"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch_and_bound_v3",
                    "locked": true,
                    "points": 8,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch_and_bound_v3\n",
                "counting_bound = CallCounter(bound)\n",
                "assert best_tree_branch_and_bound(v3_alignment, v3_sequence_names, branch, counting_bound) == (58, '(D,PA,(PB,(C1,C2)));')\n",
                "assert counting_bound.num_calls() == 19\n",
                "print(\"SUCCESS: branch_and_bound_v3 passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch_and_bound_v3_medium",
                    "locked": true,
                    "points": 7,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch_and_bound_v3_medium\n",
                "counting_bound = CallCounter(bound)\n",
                "assert best_tree_branch_and_bound(v3_medium_alignment, v3_medium_sequence_names, branch, counting_bound) == (75, '(C09,(C35,(PD1,PD3)),(PA5,(PB6,PB8)));')\n",
                "assert counting_bound.num_calls() == 133\n",
                "print(\"SUCCESS: branch_and_bound_v3 passed all tests\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "branch_and_bound_v3_big",
                    "locked": true,
                    "points": 5,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for branch_and_bound_v3_big\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
}