IAB-notebooks/2/3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2.3]() Generalized dynamic programming for multiple sequence alignment <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L2' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Table of Contents**\n",
    "0. [Progressive alignment](#1)\n",
    "    0. [Building the guide tree](#1.1)\n",
    "    0. [Generalization of Needleman-Wunsch (with affine gap scoring) for progressive multiple sequence alignment](#1.2)\n",
    "    0. [Putting it all together: progressive multiple sequence alignment](#1.3)\n",
    "0. [Progressive alignment versus iterative alignment](#2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Until now we worked with alignments between two sequences, but it is likely that you will want to align many sequences at the same time. For example, if you are trying to gain insight on the evolutionary relationships between all of the 16S bacterial genes in a given sample, it would be time consuming and very inefficient to compare them two at a time. It would be more efficient and useful to compare all of the 16S sequences from the bacteria in the same alignment.\n",
    "In the pairwise sequence alignment chapter, we went over dynamic programming algorithms. It's possible to generalize Smith-Waterman and Needleman-Wunsch, the dynamic programming algorithms that we explored for pairwise sequence alignment, to identify the optimal alignment of more than two sequences. Remember that our scoring scheme for pairwise alignment with Smith-Waterman looked like the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "& F(0, 0) = 0\\\\\n",
    "& F(i, 0) = F(i-1, 0) - d\\\\\n",
    "& F(0, j) = F(0, j-1) - d\\\\\n",
    "\\\\\n",
    "& F(i, j) = max \\begin{pmatrix}\n",
    "& F(i-1, j-1) + s(c_i, c_j)\\\\\n",
    "& F(i-1, j) - d\\\\\n",
    "& F(i, j-1) - d)\\\\\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generalize this to three sequences, we could create $3 \\times 3$ scoring, dynamic programming, and traceback matrices. Our scoring scheme would then look like the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "& F(0, 0, 0) = 0\\\\\n",
    "& F(i, 0, 0) = F(i-1, 0, 0) - d\\\\\n",
    "& F(0, j, 0) = F(0, j-1, 0) - d\\\\\n",
    "& F(0, 0, k) = F(0, 0, k-1) - d\\\\\n",
    "\\\\\n",
    "& F(i, j, k) = max \\begin{pmatrix}\n",
    "F(i-1, j-1, k-1) + s(c_i, c_j) + s(c_i, c_k) + s(c_j, c_k)\\\\\n",
    "F(i, j-1, k-1) + s(c_j, c_k) - d\\\\\n",
    "F(i-1, j, k-1) + s(c_i, c_k) - d\\\\\n",
    "F(i-1, j-1, k) + s(c_i, c_j) - d\\\\\n",
    "F(i, j, k-1) - 2d\\\\\n",
    "F(i, j-1, k) - 2d\\\\\n",
    "F(i-1, j, k) - 2d\\\\\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However the complexity of this algorithm is much worse than for pairwise alignment. For pairwise alignment, remember that if aligning two sequences of lengths $m$ and $n$, the runtime of the algorithm will be proportional to $m \\times n$. If $n$ is longer than or as long as $m$, we simplify the statement to say that the runtime of the algorithm will be be proportional to $n^2$. This curve has a pretty scary trajectory: runtime for pairwise alignment with dynamic programming is said to scale quadratically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "from functools import partial\n",
    "from IPython.core import page\n",
    "page.page = print"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "seq_lengths = range(25)\n",
    "s2_times = [t ** 2 for t in range(25)]\n",
    "\n",
    "plt.plot(range(25), s2_times)\n",
    "plt.xlabel('Sequence Length')\n",
    "plt.ylabel('Runtime (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exponent in the $n^2$ term comes from the fact that, in pairwise alignment, if we assume our sequences are both of length $n$, there are $n \\times n$ cells to fill in in the dynamic programming matrix. If we were to generalize either Smith-Waterman or Needleman-Wunsch to three sequences, we would need to create a 3 dimensional array to score and trace back the alignment. For sequences of length $n$, we would therefore have $n \\times n \\times n$ cells to fill in, and our runtime versus sequence length curve would look like the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "s3_times = [t ** 3 for t in range(25)]\n",
    "\n",
    "plt.plot(range(25), s3_times)\n",
    "plt.xlabel('Sequence Length')\n",
    "plt.ylabel('Runtime (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That curve looks steeper than the curve for pairwise alignment, and the values on the y-axis are bigger, but it's not really clear how much of a problem this is until we plot runtime for three sequences in the context of the run times for pairwise alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(range(25), s2_times)\n",
    "plt.plot(range(25), s3_times)\n",
    "plt.xlabel('Sequence Length')\n",
    "plt.ylabel('Runtime (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for four sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "s4_times = [t ** 4 for t in range(25)]\n",
    "\n",
    "plt.plot(range(25), s2_times)\n",
    "plt.plot(range(25), s3_times)\n",
    "plt.plot(range(25), s4_times)\n",
    "plt.xlabel('Sequence Length')\n",
    "plt.ylabel('Runtime (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We clearly have a problem here, and that is that the runtime for multiple sequence alignment using full dynamic programming algorithms grows exponentially with the number of sequences to be aligned. If $n$ is our sequence length, and $s$ is the number of sequences, that means that runtime is proportional to $n^s$. In pairwise alignment, $s$ is always equal to 2, so the problem is more manageable. However, for the general case of $s$ sequences, we really can't even consider Smith-Waterman or Needleman-Wunsch for more than just a few sequences. The pattern in the plots above should illustrate why."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we explored with database searching, we need to figure out how to align fewer sequences. This is where *progressive alignment* comes in."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [2.3.1](#1) Progressive alignment<a name='1'></a> <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L97' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Table of Contents**\n",
    "0. [Building the guide tree](#1.1)\n",
    "0. [Generalization of Needleman-Wunsch (with affine gap scoring) for progressive multiple sequence alignment](#1.2)\n",
    "0. [Putting it all together: progressive multiple sequence alignment](#1.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In progressive alignment, the problem of exponential growth of runtime and space is managed by selectively aligning pairs of sequences, and aligning alignments of sequences. What we typically do is identify a pair of closely related sequences, and align those. Then, we identify the next most closely related sequence to that initial pair, and align that sequence to the alignment. This concept of aligning a sequence to an alignment is new, and we'll come back to it in just a few minutes. The other concept of identifying the most closely related sequences, and then the next most closely related sequence, and so on should sound familiar. It effectively means that we're traversing a tree. And herein lies our problem: we need a tree to efficiently align multiple sequences, but we need an alignment to build a good tree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You probably have two burning questions in your mind right now:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. How do we build a tree to guide the alignment process, if we need an alignment to build a good tree?\n",
    "2. How do we align a sequence to an alignment, or an alignment to an alignment?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll explore both of those through-out the rest of this notebook. First, let's cover the process of progressive multiple sequence alignment, just assuming for a moment that we know how to do both of those things."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process of progressive multiple sequence alignment could look like the following. First, we start with some sequences and a tree representing the relationship between those sequences. We'll call this our guide tree, because it's going to guide us through the process of multiple sequence alignment. In progressive multiple sequence alignment, we build a multiple sequence alignment for each internal node of the tree, where the alignment at a given internal node contains all of the sequences in the clade defined by that node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/gregcaporaso/An-Introduction-To-Applied-Bioinformatics/master/book/fundamentals/images/msa-tree-input.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the root node, descend the bottom branch of the tree until you get to the an internal node. If an alignment hasn't been constructed for that node yet, continue descending the tree until to get to a pair of nodes. In this case, we follow the two branches to the tips. We then align the sequences at that pair of tips (usually with Needleman-Wunsch, for multiple sequence alignment), and assign that alignment to the node connecting those tips."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/gregcaporaso/An-Introduction-To-Applied-Bioinformatics/master/book/fundamentals/images/msa-tree-a1.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we want to find what to align the resulting alignment to, so start from the root node and descend the top branch of the tree. When you get to the next node, determine if an alignment has already been created for that node. If not, our job is to build that alignment so we have something to align against. In this case, that means that we need to align `s1`, `s2`, and `s3`. We can achieve this by aligning `s1` and `s3` first, to get the alignment at the internal node connecting them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/gregcaporaso/An-Introduction-To-Applied-Bioinformatics/master/book/fundamentals/images/msa-tree-a2.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can next align the alignment of `s1` and `s3` with `s2`, to get the alignment at the internal node connecting those clades."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/gregcaporaso/An-Introduction-To-Applied-Bioinformatics/master/book/fundamentals/images/msa-tree-a3.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, we can compute the alignment at the root node of the tree, by aligning the alignment of `s1`, `s2`, and `s3` with the alignment of `s4` and `s5`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/gregcaporaso/An-Introduction-To-Applied-Bioinformatics/master/book/fundamentals/images/msa-tree-final.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alignment at the root node is our multiple sequence alignment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [2.3.1.1](#1.1) Building the guide tree<a name='1.1'></a> <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L130' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Let's address the first of our outstanding questions. I mentioned above that *we need an alignment to build a good tree*. The key word here is *good*. We can build a very rough tree - one that we would never want to present as representing the actual relationships between the sequences in question - without first aligning the sequences. Remember that building a UPGMA tree requires only a distance matrix, so if we can find a non-alignment-dependent way to compute distances between the sequences, we can build a rough UPGMA tree from them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute distances between the sequences based on their *word* composition. We'll define a *word* here as `k` adjacent characters in the sequence. We can then define a function that will return all of the words in a sequence as follows. These words can be defined as being overlapping, or non-overlapping. We'll go with overlapping for this example, as the more words we have, the better our guide tree should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio import DNA\n",
    "%psource DNA.iter_kmers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in DNA(\"ACCGGTGACCAGTTGACCAGTA\").iter_kmers(3):\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in DNA(\"ACCGGTGACCAGTTGACCAGTA\").iter_kmers(7):\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in DNA(\"ACCGGTGACCAGTTGACCAGTA\").iter_kmers(3, overlap=False):\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we then have two sequences, we can compute the word counts for each and define a distance between the sequences as the fraction of words that are unique to either sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import kmer_distance\n",
    "%psource kmer_distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use this as a distance function..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = DNA(\"ACCGGTGACCAGTTGACCAGT\")\n",
    "s2 = DNA(\"ATCGGTACCGGTAGAAGT\")\n",
    "s3 = DNA(\"GGTACCAAATAGAA\")\n",
    "\n",
    "print(s1.distance(s2, kmer_distance))\n",
    "print(s1.distance(s3, kmer_distance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to override the default to create (for example) a 5-mer distance function, we could use ``functools.partial``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "fivemer_distance = partial(kmer_distance, k=5)\n",
    "\n",
    "s1 = DNA(\"ACCGGTGACCAGTTGACCAGT\")\n",
    "s2 = DNA(\"ATCGGTACCGGTAGAAGT\")\n",
    "s3 = DNA(\"GGTACCAAATAGAA\")\n",
    "\n",
    "print(s1.distance(s2, fivemer_distance))\n",
    "print(s1.distance(s3, fivemer_distance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now apply one of these functions to build a distance matrix for a set of sequences that we want to align."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_sequences = [DNA(\"ACCGGTGACCAGTTGACCAGT\", {\"id\": \"s1\"}),\n",
    "                   DNA(\"ATCGGTACCGGTAGAAGT\", {\"id\": \"s2\"}),\n",
    "                   DNA(\"GGTACCAAATAGAA\", {\"id\": \"s3\"}),\n",
    "                   DNA(\"GGCACCAAACAGAA\", {\"id\": \"s4\"}),\n",
    "                   DNA(\"GGCCCACTGAT\", {\"id\": \"s5\"})]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio import DistanceMatrix\n",
    "\n",
    "guide_dm = DistanceMatrix.from_iterable(query_sequences, metric=kmer_distance, key='id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scikit-bio also has some basic visualization functionality for these objects. For example, we can easily visualize this object as a heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = guide_dm.plot(cmap='Greens')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can next use some functionality from SciPy to cluster the sequences with UPGMA, and print out a dendrogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster.hierarchy import average, dendrogram, to_tree\n",
    "\n",
    "for q in query_sequences:\n",
    "    print(q)\n",
    "\n",
    "guide_lm = average(guide_dm.condensed_form())\n",
    "guide_d = dendrogram(guide_lm, labels=guide_dm.ids, orientation='right',\n",
    "                     link_color_func=lambda x: 'black')\n",
    "guide_tree = to_tree(guide_lm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import guide_tree_from_sequences\n",
    "%psource guide_tree_from_sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = guide_tree_from_sequences(query_sequences, display_tree=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a guide tree, so we can move on to the next step of progressive alignment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [2.3.1.2](#1.2) Generalization of Needleman-Wunsch (with affine gap scoring) for progressive multiple sequence alignment<a name='1.2'></a> <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L234' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Next, we'll address our second burning question: aligning alignments. As illustrated above, there are basically three different types of pairwise alignment we need to support for progressive multiple sequence alignment with Needleman-Wunsch. These are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Alignment of a pair of sequences.\n",
    "2. Alignment of a sequence and an alignment.\n",
    "3. Alignment of a pair of alignments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard Needleman-Wunsch supports the first, and it is very easy to generalize it to support the latter two. The only change that is necessary is in how the alignment of two non-gap characters is scored. Recall that we previously scored an alignment of two characters by looking up the score of substitution from one to the other in a substitution matrix. To adapt this for aligning a sequence to an alignment, or for aligning an alignment to an alignment, we compute this substitution as the average score of aligning the pairs of characters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, if we want to align the alignment column from $aln1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "A\n",
    "C\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to the alignment column from $aln2$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "T\n",
    "G\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we could compute the substitution score using the matrix $m$ as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "s = \\frac{m[A][T] + m[A][G] + m[C][T] + m[C][G]}{aln1_{length} \\times aln2_{length}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code adapts our implementation of Needleman-Wunsch to support aligning a sequence to an alignment, or aligning an alignment to an alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import format_dynamic_programming_matrix, format_traceback_matrix\n",
    "from skbio.alignment._pairwise import _compute_score_and_traceback_matrices\n",
    "\n",
    "%psource _compute_score_and_traceback_matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio.alignment._pairwise import _traceback\n",
    "%psource _traceback"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio.alignment import global_pairwise_align_nucleotide\n",
    "%psource global_pairwise_align_nucleotide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the sake of the examples below, I'm going to override one of the ``global_pairwise_align_nucleotide`` defaults to penalize terminal gaps. This effectively tells the algorithm that we know we have a collection of sequences that are homologous from beginning to end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "global_pairwise_align_nucleotide = partial(global_pairwise_align_nucleotide, penalize_terminal_gaps=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, we can still use this code to align pairs of sequences (but note that we now need to pass those sequences in as a pair of one-item lists):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "aln1, _, _ = global_pairwise_align_nucleotide(query_sequences[0], query_sequences[1])\n",
    "print(aln1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can align that alignment to one of our other sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "aln1, _, _ = global_pairwise_align_nucleotide(aln1, query_sequences[2])\n",
    "print(aln1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can align another pair of sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "aln2, _, _ = global_pairwise_align_nucleotide(query_sequences[2], query_sequences[3])\n",
    "print(aln2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then align that alignment against our previous alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "aln3, _, _  = global_pairwise_align_nucleotide(aln1, aln2)\n",
    "print(aln3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [2.3.1.3](#1.3) Putting it all together: progressive multiple sequence alignment<a name='1.3'></a> <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L317' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can now combine all of these steps to take a set of query sequences, build a guide tree, perform progressive multiple sequence alignment, and return the guide tree (as a SciPy linkage matrix) and the alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio import TreeNode\n",
    "guide_tree = TreeNode.from_linkage_matrix(guide_lm, guide_dm.ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can view the guide tree in [Newick format](http://scikit-bio.org/docs/latest/generated/skbio.io.newick.html) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(guide_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import progressive_msa\n",
    "%psource progressive_msa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa = progressive_msa(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, guide_tree=guide_tree)\n",
    "print(msa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now build a (hopefully) improved tree from our multiple sequence alignment. First we'll look at our original distance matrix again, and then the distance matrix generated from the progressive multiple sequence alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = guide_dm.plot(cmap='Greens')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa_dm = DistanceMatrix.from_iterable(msa, metric=kmer_distance)\n",
    "fig = msa_dm.plot(cmap='Greens')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The UPGMA trees that result from these alignments are very different. First we'll look at the guide tree, and then the tree resulting from the progressive multiple sequence alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = dendrogram(guide_lm, labels=guide_dm.ids, orientation='right',\n",
    "               link_color_func=lambda x: 'black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa_lm = average(msa_dm.condensed_form())\n",
    "d = dendrogram(msa_lm, labels=msa_dm.ids, orientation='right',\n",
    "               link_color_func=lambda x: 'black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can wrap this all up in a single convenience function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import progressive_msa_and_tree\n",
    "%psource progressive_msa_and_tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa = progressive_msa(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, guide_tree=guide_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa, tree = progressive_msa_and_tree(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide,\n",
    "                                     display_tree=True, display_aln=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [2.3.2](#2) Progressive alignment versus iterative alignment<a name='2'></a> <a class='iab-edit' href='https://github.com/caporaso-lab/An-Introduction-to-Applied-Bioinformatics/edit/master/book/fundamentals/multiple-sequence-alignment.md#L382' target='_blank'>[edit]</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In an iterative alignment, the output tree from the above progressive alignment is used as a guide tree, and the full process repeated. This is performed to reduce errors that result from a low-quality guide tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iab.algorithms import iterative_msa_and_tree\n",
    "%psource iterative_msa_and_tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa, tree = iterative_msa_and_tree(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, num_iterations=1, display_aln=True, display_tree=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa, tree = iterative_msa_and_tree(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, num_iterations=2, display_aln=True, display_tree=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa, tree = iterative_msa_and_tree(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, num_iterations=3, display_aln=True, display_tree=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa, tree = iterative_msa_and_tree(query_sequences, pairwise_aligner=global_pairwise_align_nucleotide, num_iterations=5, display_aln=True, display_tree=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some references that I used in assembling these notes include [1](http://statweb.stanford.edu/~nzhang/345_web/sequence_slides3.pdf), [2](http://math.mit.edu/classes/18.417/Slides/alignment.pdf), [3](http://www.sciencedirect.com/science/article/pii/0378111988903307), [4](http://bioinformatics.oxfordjournals.org/content/23/21/2947.full), and [5](http://nar.oxfordjournals.org/content/32/5/1792.full)."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 4
}