lec3-assembly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 3:  De-Novo Genome Assembly (基因组从头拼接)\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "In this lecture, we will cover the fundamental concepts as well as the popular algorithms applied in de-novo genome assembly in next-generation sequencing (NGS):\n",
    "\n",
    "- [The Genome Assembly and the related concepts](#1.-Genome-assembly:-The-concepts)\n",
    "\n",
    "- [The Overlaps and Overlap graphs of the reads](#2.-Overlaps)\n",
    "\n",
    "- [The SCS approaches](#3.-The-Simplest-approach:-Shortest-common-superstring)\n",
    "\n",
    "- [The OLC approach](#4.1-OLC-approach)\n",
    "\n",
    "- [The De-Bruijn Graph Algorithm](#4.2-De-Bruijn-Graph-Algorithm)\n",
    "\n",
    "And we need to finish the following exercises:\n",
    "- [EXERCISE 1](#$\\S$-EXERCISE-1)\n",
    "- [EXERCISE 2](#$\\S$-EXERCISE-2)\n",
    "- [EXERCISE 3](#$\\S$-EXERCISE-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Genome assembly: The concepts\n",
    "\n",
    "Whole-genome \"shotgun\" sequencing starts by amplifying and fragmenting the genome DNA. Here **shotgun** refers to the **random fragmentation** of the whole genome sequence, as if it was fired from a gun.\n",
    "\n",
    "![](../images/shotgun.png)\n",
    "\n",
    "Now the novel sequencing technologies of long read up to 1000 bases become available, which makes it easier of the  assembly task. But here we will focus on the short reads assemly which generated by comtemporary sequencing instruments, which is, the **second-generation sequencing**, while the former is known as the **third-generation sequencing**.\n",
    "\n",
    "So now we have in hands the fragmented reads and the notion is that we know that they're going to align them up like a puzzle.\n",
    "\n",
    "In next-generation sequencing, almost all genome positions are **covered** by many reads. The assembly task will be to reconstruct the genome from these short reads.  But the number of reads is perhaps 200 million or even a billion. And thus it's quite a touch task to put a huge number of pieces together if we don't really know where they are coming from.  Because we don't have the genome (reference) to guide us.\n",
    "\n",
    "That's called the **de-novo assembly**, starting from the scratch. You are given a set of reads for your organism of interest, and we're going to assemble it into contigs or even scaffolds today.\n",
    "\n",
    "This approach is quite different from the so-called **reference-guided assembly**, which is a model-based assembly. Whereas de-novo assembly is an non-model approach. For guided assembly, you can map the reads back to the reference genome, so it's pretty easier.\n",
    "\n",
    "However, as you can tell, there is a large structural variation between the reference genome and the resequencing genomes, that process might fail.\n",
    "\n",
    "So today we'll focus upon **De-Novo Genome Assembly**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Key term: Coverage\n",
    "\n",
    "Usually the **coverage** is short for the **mean coverage**, which is the average number of reads covering a position in the genome.\n",
    "\n",
    "Assume that the genome size is $G$, and the read length is $R$. If we have N reads, then the coverage  can be computed as\n",
    "$$\n",
    "\\text{coverage} = \\frac{N \\times R}{G}\n",
    "$$\n",
    "\n",
    "So if $N= 11$, $R=17$ and $G=35$, then the coverage should be $\\frac{17 \\times 11}{35} \\approx 5.3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Estimating Uncovered Bases\n",
    "\n",
    "This algorithm is proposed by Lander and Waterman.\n",
    "\n",
    "As we shown before, the coverage \n",
    "$$\n",
    "\\lambda = \\frac{NR}{G} = \\frac{\\text{length of all reads}}{\\text{length of genome}}\n",
    "$$\n",
    "\n",
    "Therefore, the probability that a base  is not covered can be approximated using a Poisson distribution:\n",
    "$$\n",
    "Poisson(0, \\lambda) \\approx e^{-\\lambda}\n",
    "$$\n",
    "\n",
    "Thus the number of uncovered bases should be\n",
    "$$\n",
    "u \\approx G e^{-\\lambda}\n",
    "$$\n",
    "\n",
    "And the nuber of gaps:\n",
    "$$\n",
    "g \\approx Ne^{-\\lambda}\n",
    "$$\n",
    "\n",
    "But unfortunately, the Poisson approximation is not a good fit for the real data, since the library is biased. Therefore, we might use the other approximation techniques like negative Binomial approximation.\n",
    "\n",
    "Coverage could also refer to the number of reads covering a particular position in the genome."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\">$\\S$ EXERCISE 1</font>\n",
    "\n",
    "Can you write a python program, to compute the per-base coverage for a set of reads and also the `*.sai` files generated by the alignment software `bwa` against the corresponding reference genome."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Sequencing error\n",
    "\n",
    "Say two reads truly originate from overlapping stretches of the genome. Why might there be differences? There are two possible reasons:\n",
    "\n",
    "- The first reason is that there is an error in either of the reads. So quality score can help us decide which is the truth.\n",
    "- The second is due to the difference between inherited copies of a chromosome, from paternal or maternal parent. We know that human are diploid. And these two copies can differ.\n",
    "\n",
    "\n",
    "We can call them allelic difference, which will pop up in terms of non-concordance of our reads. But remember for human and the other organisms, the reference genome is haploid.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2. Overlaps\n",
    "\n",
    "Overlaps are essential for genome assembly. WE CANNOT REBUILD THE GENOME UNLESS THERE ARE OVERLAPS BETWEEN THE SEQUENCED READS.\n",
    "\n",
    "Finding all overlaps is like building a **directed graph** where directed edges connect overlapping nodes (reads) where the *suffix of sources is similar to the prefix of sink*.\n",
    "\n",
    "### 2.1 Directed graph review\n",
    "\n",
    "Directed graph $G(V, E)$ consists of set of *vetices*, $V$ and set of *directed edges*, $E$.\n",
    "\n",
    "Directed edge is an *ordered pair of vertices* with first the **source** and second the **sink**.\n",
    "\n",
    "In a typical graph, the vertices (a.k.a node) are drawn as circles, while the edges (a.k.a arcs) are drawn as lines with an arrow connecting two circles.\n",
    "\n",
    "Directed graphs are also called *digraphs*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 How to compute the longest overlap between two reads?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def overlap(a, b, min_length=3):\n",
    "    \"\"\"\n",
    "    Compute the maximum overlap length of suffix of a and prefix of b.\n",
    "    @PARAMETERS\n",
    "    -a: string 1\n",
    "    -b: string 2\n",
    "    -min_length: the minimum overlap\n",
    "    @RETURN\n",
    "    The maximum overlap between a and b\n",
    "    \"\"\"\n",
    "    start = 0\n",
    "    \n",
    "    while True:\n",
    "        start = a.find(b[:min_length], start) # find the match position for b[:min_length] in a[start:]\n",
    "        if start == -1:\n",
    "            return 0\n",
    "        if b.startswith(a[start:]):                          # if perfect match\n",
    "            return len(a) - start\n",
    "        start += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function will return the length of the longest overlap between $a$ and $b$. If it is less than `minimum_length`, it will return 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6\n"
     ]
    }
   ],
   "source": [
    "seqa = \"ATCGATA\"\n",
    "seqb = \"TCGATAC\"\n",
    "overlaps = overlap(seqa, seqb)\n",
    "print overlaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 2.3 The overlap graph (重叠图)\n",
    "\n",
    "In an overlap graph, a vertex is a read; a directed edge is an overlap of at least 3 characters matching between a suffix of the source and the prefix of a sink.\n",
    "\n",
    "Here is an example:\n",
    "- Vertices (reads): $\\{a: \"CTCTAGGCC\", b:\"GCCCTCAAT\", c:\"CAATTTTT\"\\}$\n",
    "- Edges(overlaps): $\\{(a, b), (b, c)\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use the `permutations` function in the module `itertools` to construct the overlap map for a set of reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from itertools import permutations\n",
    "def naive_overlap_map(reads, k):\n",
    "    \"\"\"functions to get all the overlaps between two arbitrary reads with overlaps >= k\"\"\"\n",
    "    olaps = {}\n",
    "    for a, b in permutations(reads, 2):\n",
    "        s = overlap(a, b, k)\n",
    "        if s > 0:\n",
    "            olaps[(a, b)] = s\n",
    "    return olaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACGATTACTG          GATTACTGACG             8\n",
      "GATTACTGACG         TGACGGACT               5\n"
     ]
    }
   ],
   "source": [
    "reads = [\"ACGATTACTG\",\"GATTACTGACG\",\"CGATTTACG\",\"CCAGTAG\",\"TGACGGACT\"]\n",
    "maps = naive_overlap_map(reads, k=5)\n",
    "for key in maps:\n",
    "    print \"{:20}{:20}{:5}\".format(key[0],  key[1],  maps[key])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. The Simplest approach: Shortest common superstring\n",
    "\n",
    "One of the simplest algorithms to assemble the reads is to find the shortest common superstring for these input reads. What is the shortest common superstring (SCS)? \n",
    "> **SCS is the shortest string that contains all the input strings as the substrings**.\n",
    "\n",
    "\n",
    "### 3.1 Brute-force algorithm to find the SCS\n",
    "\n",
    "How to find the SCS for the input string? One **brute-force approach** is to loop through all the possible combinations (orderings) of the input reads and find the shortest combination:\n",
    "1. Find all the possible orderings for the input reads using the `itertools.permutations()` function;\n",
    "2. Compute the length for the superstring;\n",
    "3. Output the shortest one as the assembly result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def scs(reads):\n",
    "    \n",
    "    shortest_sup = None\n",
    "    for comb in permutations(reads):\n",
    "        sup = comb[0]\n",
    "        for i in xrange(len(reads) - 1):\n",
    "            olen = overlap(comb[i], comb[i+1], min_length=1)\n",
    "            sup +=  comb[i+1][olen:]\n",
    "            \n",
    "        if shortest_sup is None or len(sup) < len(shortest_sup):\n",
    "            shortest_sup = sup\n",
    "            \n",
    "    return shortest_sup    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACTCAGGACTTACGGAGGCTACGGACTGGA\n"
     ]
    }
   ],
   "source": [
    "print scs([\"ACTCAGGACTTACG\", \"ACTTACGGAGGCT\", \"GCTACGGACTGGA\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This brute-force algorithm can absolutely find the optimal solution. However, for a set of $n$ reads, the number of possible combinations is $n!$. When $n$ is large, the huge computational cost is unfeasible and non-practical for our task of assembling the short reads."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Greedy shortest common superstring\n",
    "\n",
    "Because brute-force approach is not pratical, let's turn to the other method, the greedy SCS, which is much faster than the former approach.\n",
    "\n",
    "> We call a method greedy, because the algorithm make a series of decison and at each decision point, it will choose the option it will reduce the most the length of the eventual superstring.\n",
    "\n",
    "This seems like a good strategy. The shorter the superstring, the closer we are to the shortest common substring.\n",
    "\n",
    "We can visualize the greedy SCS algorithm using an overlap graph. So, let's start with the overlap graph here:\n",
    "![](../images/greedy.png)\n",
    "\n",
    "Let's have a look how it works:\n",
    "\n",
    "1. In each round, find the longest overlap in the remaining overlap graph (that is, pick the edge with the greatest number as its label)\n",
    "\n",
    "2. Merge the two nodes on either terminal of the edge, and remove the other edges that connect the two nodes as well as the edges pointing to the ingoing node, and also the edges starting from the outgoing node.\n",
    "\n",
    "3. Iterate the steps (1) and (2), until no nodes can be merged.\n",
    "\n",
    "4. Stop and output the remaining contigs.\n",
    "\n",
    "Compared with the previous brute-force approach which tries all possible orderings, the greedy aprpoach greatly reduce the computational cost. And therefore we can see that it is much quicker than the previously mentioned approach.\n",
    "\n",
    "However, the speed comes at a cost of accuracy. The greedy algorithm is not guaranteed to find the correct answer.\n",
    "![](../images/greedy-drawbacks.png)\n",
    "\n",
    "Here is the practical usage of the greedy algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pick_longest_overlap(reads, k):\n",
    "    \"\"\"Find the longest overlap among a set of reads.\n",
    "    @PARAMETERS\n",
    "    -reads: a list of reads\n",
    "    -k: the threshold of overlap gap length\n",
    "    @RETURN\n",
    "    - a tuple containing the start, end and the corresponding longest overlap\n",
    "    @Note\n",
    "    -Since there might exists ties, we only pick the first occurrence.\n",
    "    \"\"\"\n",
    "    longest_overlap = 0\n",
    "    reada, readb = None, None\n",
    "    for a, b in permutations(reads, 2):\n",
    "        olen = overlap(a, b, k)\n",
    "        if olen > longest_overlap:\n",
    "            longest_overlap = olen\n",
    "            reada, readb = a, b\n",
    "    return reada, readb, longest_overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then we iteratively pick the two reads with longest overlap, and merge these two reads, until finally no reads can be merged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def greedy_scs(reads, k):\n",
    "    \"\"\"compute the shortest common superstring from the reads using greedy algorithm.\n",
    "    @INPUT\n",
    "    \n",
    "    @OUTPUT\n",
    "    \n",
    "    @NOTE\n",
    "    \"\"\"\n",
    "    read_a, read_b, olen = pick_longest_overlap(reads, k)\n",
    "    while olen > 0:\n",
    "        reads.remove(read_a)\n",
    "        reads.remove(read_b)\n",
    "        reads.append(read_a + read_b[olen:])\n",
    "        read_a, read_b, olen = pick_longest_overlap(reads, k)\n",
    "        \n",
    "    return \"\".join(reads)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here once there is no overlap left in the remaining list of reads, we will return the simple concatenation of the reads using `\"\".join(reads)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shortest common superstring for the input reads ['AAA', 'AAB', 'ABB', 'BBA', 'BBB'] is BBBAAABBA\n"
     ]
    }
   ],
   "source": [
    "myreads = [\"AAA\",\"AAB\",\"ABB\", \"BBA\", \"BBB\"]\n",
    "print \"The shortest common superstring for the input reads %s is %s\" % (myreads,  greedy_scs(myreads, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And again, let's try the **brute-force algorithm**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shortest common superstring for the input reads ['AAA', 'AAB', 'ABB', 'BBA', 'BBB'] is AAABBBA\n"
     ]
    }
   ],
   "source": [
    "myreads = [\"AAA\",\"AAB\",\"ABB\", \"BBA\", \"BBB\"]\n",
    "print \"The shortest common superstring for the input reads %s is %s\" % (myreads,  scs(myreads))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Can you draw a conclusion from the above example and scripts?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Take-home message for greedy approach\n",
    "\n",
    "1. Greedy algorithm is not guaranteed to choose overlaps yielding SCS (shortest common superstring).\n",
    "2. But greedy algorithm is a good approximation. We do know that there are some bounds of it, i.e. the superstring yielded by the greedy algorithm won't be more than around 2.5 times longer than the true SCS.\n",
    "\n",
    "That's what a compromise does in finding the shortest common superstring."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\">$\\S$ EXERCISE 2</font>\n",
    "\n",
    "In the above implementation, we need to compute the overlaps among the reads once we update the list of reads. This, will unfortunately cause the unnecessary waste of the computational resouces.\n",
    "\n",
    "Can you fix it?\n",
    "\n",
    "Hint: Maybe you can store the overlaps into a dictionary, thus what we need to do is to update some of the terms in the dictionary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 How do repeats affect the assembly?\n",
    "\n",
    "In the previous example, we could find that greedy approach might not reach at the optimal solution. However, when the genome is repetitive, the brute-force approach can also results in a wrong answer.\n",
    "\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABAAAB\n"
     ]
    }
   ],
   "source": [
    "genome = \"ABAAAAABA\"\n",
    "myreads = [\"ABA\", \"BAA\", \"AAA\", \"AAB\"]   # the genome can only generate such 3-mers\n",
    "print scs(myreads)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Over-collapsing issues\n",
    "\n",
    "We can hardly tell the actual copy number of the repeats (here 'A'), if the read length is short. This is because we can use fewer copies of 'A' to explain the above list of reads (which we call the **over-collapsing problem**), right?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But if the read length is longer, say, 7, which can span the whole repeats section:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABAAAAABA\n"
     ]
    }
   ],
   "source": [
    "myreads2 = [\"ABAAAAA\", \"BAAAAAB\", \"AAAAABA\"]   # the genome can only generate such 3-mers\n",
    "print scs(myreads2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, if the read length is long enough, the repeats will not be an issue. **But how long is long enough?** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to resulting in fewer copies of repeats, repeats can also result in other serious problems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABCBAAABBCAA\n"
     ]
    }
   ],
   "source": [
    "genome = \"AAABCBAAABBCAAA\"\n",
    "myreads = [\"AAA\", \"AAB\", \"ABC\", \"BCB\", \"CBA\", \"BAA\", \"AAB\", \"ABB\", \"BCA\", \"CAA\"]\n",
    "print scs(myreads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ABCBAABCAAABB\n"
     ]
    }
   ],
   "source": [
    "print greedy_scs(myreads, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Two state-of-the-art approaches to short read assembly\n",
    "\n",
    "Today we will talk about two general approaches for genome assembly:\n",
    "<p><img src=\"../images/assembly.png\" height=20% /></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 OLC approach\n",
    "\n",
    "- **Overlap Layout Consensus (OLC)** is a kind of string graph assemblers, which works by\n",
    "    * Construct overlap graph directly from reads;\n",
    "    * Eliminate redundant reads;\n",
    "    * Trace path for assembly;\n",
    "    * Such tools include: SGA, Fermi for longer reads.\n",
    "    \n",
    "#### Reference for SGA\n",
    "* J. Simpson, R. Durbin. Efficient de novo assembly of large genomes using compressed data structures. Genome Res, 2011."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.1 Build the overlap graph efficiently\n",
    "\n",
    "In the early phase, people are afraid of using these type of \"OLC\" assemblers because of the huge cost of computing overlaps, which was thought to be $\\mathcal{O}(n^2)$ where $n$ is the number of reads.\n",
    "\n",
    "Fortunately, a really clever algorithm was devised that used the technology we have discussed in the previous lecture, the FM-index and Burrows-Wheeler transform that allow us to index a genome and then conduct a lookup of reads in time proportional to the length of the read, $\\mathcal{O}(n)$.\n",
    "\n",
    "Here is essentially the idea. What we are going to do is we're going to take all of the reads to index them all. And we can do that roughly in $\\mathcal{O}(n \\log n)$ time. And after we index all of the reads, then we can use the index to find overlaps very very efficiently. And you can conceptualize this as simply looking at a read that you have in your hand and looking it up in the index. And you can find all the places the suffix or prefix of that read matches. And you can trace back to find all the places it matches. \n",
    "\n",
    "And those all correspond to edges in the graph.\n",
    "\n",
    "A **merged FM-index** of all reads allows us to match read prefixes and suffixes to discover overlaps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.2 Layout\n",
    "\n",
    "Now let's turn to how to process an **overlap graph** towards making **contigs**.\n",
    "\n",
    "The overlap graph is huge and messy, and the contigs will not automatically \"pop up\" to us. Here is an actual layout-graph, can we dig the contigs out immediately using our eyes? The answer is definitely no. And that comes to the question of layout, which is, how we formulate the problem of tracing a path through an overlap graph.\n",
    "\n",
    "So the first idea is the shortest common superstring (SCS), as we mentioned earlier in this lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.3 Consensus\n",
    "\n",
    "Since some of the reads contains errors, we need to use the **majorite-vote** mechanism to choose from the possible choice for a *polymorphism*  site.\n",
    "\n",
    "Furthermore, the input reads may generate **more than one contig**. With further information, we can concatenate these contigs into a scaffold or **scaffolds**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 4.2 De-Bruijn Graph Algorithm\n",
    "\n",
    "Now let's turn to another approach for solving the assembling problem which is based on the De-Bruijn graph.\n",
    "    \n",
    "- **de Bruijn graph-based assemblers**, which works by\n",
    "    * Construct k-mer graph from reads, discarding the original reads;\n",
    "    * Trace path in graph for assembly.\n",
    "    * Appropriately for shorter reads.\n",
    "    * Such tools include: velvet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.1 The De-Bruijn graph\n",
    "\n",
    "Simiar to the overlap graph, the De-Bruijn graph is also a directed graph. The difference is that it is a multigraph, too, which means that there could be multiple edges from a node, $A$, to the other node , $B$. Let's start with an example.\n",
    "\n",
    "For a genome $G$: <font color=\"blue\">AAABBBBA</font>. And the reads are the substrings of $G$ of length $k$, which we call it the $k$-mer of the genome, $G$.\n",
    "\n",
    "So the 3-mers of $G$ are: <font color=\"red\">AAA, AAB, BBB, BBB, BBA</font>.\n",
    "\n",
    "For each $k$-mer, we're going to extract the **left and right $(k-1)$-mers**.\n",
    "![](../images/debruijn.png)\n",
    "\n",
    "Now each $(k-1)$-mer will become a node in the graph. And then we add an edge from the left $(k-1)$-mer to the right. \n",
    "\n",
    "Now, we can find that\n",
    "- Each edge represents a $k$-mer, including the identical $k$-mers with more than one occurrence.\n",
    "\n",
    "- Each node is a distinct $(k-1)$-mer.\n",
    "\n",
    "Let me ask you a question:\n",
    "> Can we rebuild the genome only using the above De-Bruijn graph information?\n",
    "\n",
    "The answer is yes. But how? This is the core for the De-Bruijn graph assembler.\n",
    "\n",
    "> **The reconstruction of the original genome corresponds to a walk-through of the graph.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2.2 Eulerian walk\n",
    "\n",
    "Start at some node in this De-Bruijn graph, and then move through the graph from node to node, following the edges as we go, respecting the direction of each edge.  And then we can reconstruct the genome as we go.\n",
    "\n",
    "> ### The walk that allows us to reconstruct the genome, is the walk that crosses each edge exactly once (that is, uses each $k$-mer exactly once).\n",
    "\n",
    "![](../images/debruijn-walk.png)\n",
    "\n",
    "This, has a special name,  <font size=3 weight=\"bold\">Eulerian walk</font>.\n",
    "\n",
    "**But not every De-Bruijn graph can be written as a Eulerian walk.** But if the graph does, we call it a Eulerian graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 4.2.3 The De-Bruijn graph algorithm\n",
    "\n",
    "- Given the sequenced reads, we construct the corresponding De-Bruijn graph.\n",
    "\n",
    "- The graph is Eulerian as long as we have one read for every $k$-mer.\n",
    "\n",
    "- An Eulerian walk-through of the graph will give us a reconstruction of the original genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def DeBruijnGaph(st, k):\n",
    "    \"\"\"Build the de Bruijn graph based on `k`-mers of `st`.\n",
    "    \"\"\"\n",
    "    nodes = set()\n",
    "    edges = []\n",
    "    for i in range(len(st) - k + 1):\n",
    "        edges.append((st[i:i+k-1], st[i+1:i+k]))    # left and right k-1 mer of the k-mer\n",
    "        nodes.add(st[i:i+k-1])\n",
    "        nodes.add(st[i+1:i+k])\n",
    "    return nodes, edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set(['ACC', 'ACG', 'ATC', 'AGG', 'CAG', 'CTG', 'GGA', 'ACT', 'CAC', 'GAC', 'TGA', 'CCA', 'TCA'])\n"
     ]
    }
   ],
   "source": [
    "genome = \"ATCAGGACTGACCACG\"\n",
    "nodes, edges = DeBruijnGraph(genome,  4)\n",
    "print nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can visualize the graph like that using graphviz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ViewDeBruijnGraph(st, k):\n",
    "    \"\"\"Visualize a directed multigraph using graphviz.\n",
    "    \"\"\"\n",
    "    nodes, edges = de_bruijn(st, k)\n",
    "    dot_str = 'digraph \"DeBruijn Graph\" {\\n'\n",
    "    for node in nodes:\n",
    "        dot_str +=  '    %s [label=\"%s\"]   ;\\n' %(node, node)\n",
    "    for src, dest in edges:\n",
    "        dot_str +=  '    %s -> %s  ;\\n' %(src, dest)\n",
    "    return dot_str  +   '}\\n'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The gvmagic extension is already loaded. To reload it, use:\n",
      "  %reload_ext gvmagic\n"
     ]
    }
   ],
   "source": [
    "%load_ext gvmagic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You need to pre-install the `gvmagic` using `pip`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.30.1 (20150305.2338)\n",
       " -->\n",
       "<!-- Title: DeBruijn Graph Pages: 1 -->\n",
       "<svg width=\"180pt\" height=\"692pt\"\n",
       " viewBox=\"0.00 0.00 180.00 692.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 688)\">\n",
       "<title>DeBruijn Graph</title>\n",
       "<polygon fill=\"white\" stroke=\"white\" points=\"-4,5 -4,-688 177,-688 177,5 -4,5\"/>\n",
       "<!-- ACC -->\n",
       "<g id=\"node1\" class=\"node\"><title>ACC</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"142\" cy=\"-234\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"142\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">ACC</text>\n",
       "</g>\n",
       "<!-- CCA -->\n",
       "<g id=\"node12\" class=\"node\"><title>CCA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"142\" cy=\"-162\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"142\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">CCA</text>\n",
       "</g>\n",
       "<!-- ACC&#45;&gt;CCA -->\n",
       "<g id=\"edge11\" class=\"edge\"><title>ACC&#45;&gt;CCA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142,-215.697C142,-207.983 142,-198.712 142,-190.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"145.5,-190.104 142,-180.104 138.5,-190.104 145.5,-190.104\"/>\n",
       "</g>\n",
       "<!-- ACG -->\n",
       "<g id=\"node2\" class=\"node\"><title>ACG</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"142\" cy=\"-18\" rx=\"29.3722\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"142\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">ACG</text>\n",
       "</g>\n",
       "<!-- ATC -->\n",
       "<g id=\"node3\" class=\"node\"><title>ATC</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-666\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-662.3\" font-family=\"Times,serif\" font-size=\"14.00\">ATC</text>\n",
       "</g>\n",
       "<!-- TCA -->\n",
       "<g id=\"node13\" class=\"node\"><title>TCA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-594\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-590.3\" font-family=\"Times,serif\" font-size=\"14.00\">TCA</text>\n",
       "</g>\n",
       "<!-- ATC&#45;&gt;TCA -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>ATC&#45;&gt;TCA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85,-647.697C85,-639.983 85,-630.712 85,-622.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.5001,-622.104 85,-612.104 81.5001,-622.104 88.5001,-622.104\"/>\n",
       "</g>\n",
       "<!-- AGG -->\n",
       "<g id=\"node4\" class=\"node\"><title>AGG</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-450\" rx=\"29.3479\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-446.3\" font-family=\"Times,serif\" font-size=\"14.00\">AGG</text>\n",
       "</g>\n",
       "<!-- GGA -->\n",
       "<g id=\"node7\" class=\"node\"><title>GGA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-378\" rx=\"29.3479\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-374.3\" font-family=\"Times,serif\" font-size=\"14.00\">GGA</text>\n",
       "</g>\n",
       "<!-- AGG&#45;&gt;GGA -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>AGG&#45;&gt;GGA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85,-431.697C85,-423.983 85,-414.712 85,-406.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.5001,-406.104 85,-396.104 81.5001,-406.104 88.5001,-406.104\"/>\n",
       "</g>\n",
       "<!-- CAG -->\n",
       "<g id=\"node5\" class=\"node\"><title>CAG</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-522\" rx=\"29.3722\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-518.3\" font-family=\"Times,serif\" font-size=\"14.00\">CAG</text>\n",
       "</g>\n",
       "<!-- CAG&#45;&gt;AGG -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>CAG&#45;&gt;AGG</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85,-503.697C85,-495.983 85,-486.712 85,-478.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.5001,-478.104 85,-468.104 81.5001,-478.104 88.5001,-478.104\"/>\n",
       "</g>\n",
       "<!-- CTG -->\n",
       "<g id=\"node6\" class=\"node\"><title>CTG</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"28\" cy=\"-162\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"28\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">CTG</text>\n",
       "</g>\n",
       "<!-- TGA -->\n",
       "<g id=\"node11\" class=\"node\"><title>TGA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"48\" cy=\"-90\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"48\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">TGA</text>\n",
       "</g>\n",
       "<!-- CTG&#45;&gt;TGA -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>CTG&#45;&gt;TGA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M32.8414,-144.055C35.1016,-136.145 37.8458,-126.54 40.3749,-117.688\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"43.8212,-118.366 43.2031,-107.789 37.0905,-116.443 43.8212,-118.366\"/>\n",
       "</g>\n",
       "<!-- GAC -->\n",
       "<g id=\"node10\" class=\"node\"><title>GAC</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"85\" cy=\"-306\" rx=\"29.3722\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"85\" y=\"-302.3\" font-family=\"Times,serif\" font-size=\"14.00\">GAC</text>\n",
       "</g>\n",
       "<!-- GGA&#45;&gt;GAC -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>GGA&#45;&gt;GAC</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85,-359.697C85,-351.983 85,-342.712 85,-334.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.5001,-334.104 85,-324.104 81.5001,-334.104 88.5001,-334.104\"/>\n",
       "</g>\n",
       "<!-- ACT -->\n",
       "<g id=\"node8\" class=\"node\"><title>ACT</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"28\" cy=\"-234\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"28\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">ACT</text>\n",
       "</g>\n",
       "<!-- ACT&#45;&gt;CTG -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>ACT&#45;&gt;CTG</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M28,-215.697C28,-207.983 28,-198.712 28,-190.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"31.5001,-190.104 28,-180.104 24.5001,-190.104 31.5001,-190.104\"/>\n",
       "</g>\n",
       "<!-- CAC -->\n",
       "<g id=\"node9\" class=\"node\"><title>CAC</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"142\" cy=\"-90\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"142\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">CAC</text>\n",
       "</g>\n",
       "<!-- CAC&#45;&gt;ACG -->\n",
       "<g id=\"edge13\" class=\"edge\"><title>CAC&#45;&gt;ACG</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142,-71.6966C142,-63.9827 142,-54.7125 142,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"145.5,-46.1043 142,-36.1043 138.5,-46.1044 145.5,-46.1043\"/>\n",
       "</g>\n",
       "<!-- GAC&#45;&gt;ACT -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>GAC&#45;&gt;ACT</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M72.3499,-289.465C64.8965,-280.311 55.3225,-268.554 46.9786,-258.307\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"49.6431,-256.036 40.6148,-250.492 44.2151,-260.456 49.6431,-256.036\"/>\n",
       "</g>\n",
       "<!-- GAC&#45;&gt;ACC -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>GAC&#45;&gt;ACC</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M97.6501,-289.465C105.104,-280.311 114.678,-268.554 123.021,-258.307\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"125.785,-260.456 129.385,-250.492 120.357,-256.036 125.785,-260.456\"/>\n",
       "</g>\n",
       "<!-- TGA&#45;&gt;GAC -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>TGA&#45;&gt;GAC</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.3353,-107.964C57.9539,-118.257 62.2955,-131.732 65,-144 75.2423,-190.459 80.6755,-245.699 83.1869,-277.907\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.6981,-278.188 83.9293,-287.902 86.6788,-277.67 79.6981,-278.188\"/>\n",
       "</g>\n",
       "<!-- CCA&#45;&gt;CAC -->\n",
       "<g id=\"edge12\" class=\"edge\"><title>CCA&#45;&gt;CAC</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142,-143.697C142,-135.983 142,-126.712 142,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"145.5,-118.104 142,-108.104 138.5,-118.104 145.5,-118.104\"/>\n",
       "</g>\n",
       "<!-- TCA&#45;&gt;CAG -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>TCA&#45;&gt;CAG</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85,-575.697C85,-567.983 85,-558.712 85,-550.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"88.5001,-550.104 85,-540.104 81.5001,-550.104 88.5001,-550.104\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%dotstr ViewDeBruijnGraph(genome, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 4.2.4 When will Eulerian walk go wrong: The curse of the repetitive genome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In the previous section, we have learned how to build the De-Bruijn graph and how to walk Eulerian to rebuild the original genome. And let's see whether it has resolved the over-collapsing problem of repeats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.30.1 (20150305.2338)\n",
       " -->\n",
       "<!-- Title: DeBruijn Graph Pages: 1 -->\n",
       "<svg width=\"170pt\" height=\"620pt\"\n",
       " viewBox=\"0.00 0.00 170.47 620.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 616)\">\n",
       "<title>DeBruijn Graph</title>\n",
       "<polygon fill=\"white\" stroke=\"white\" points=\"-4,5 -4,-616 167.466,-616 167.466,5 -4,5\"/>\n",
       "<!-- ng_l -->\n",
       "<g id=\"node1\" class=\"node\"><title>ng_l</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"95\" cy=\"-594\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"95\" y=\"-590.3\" font-family=\"Times,serif\" font-size=\"14.00\">ng_l</text>\n",
       "</g>\n",
       "<!-- g_lo -->\n",
       "<g id=\"node8\" class=\"node\"><title>g_lo</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-522\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-518.3\" font-family=\"Times,serif\" font-size=\"14.00\">g_lo</text>\n",
       "</g>\n",
       "<!-- ng_l&#45;&gt;g_lo -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>ng_l&#45;&gt;g_lo</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M77.1637,-580.161C66.7237,-570.823 53.8253,-557.911 43.6145,-546.721\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"45.9808,-544.114 36.728,-538.952 40.7428,-548.757 45.9808,-544.114\"/>\n",
       "</g>\n",
       "<!-- ng_l&#45;&gt;g_lo -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>ng_l&#45;&gt;g_lo</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M85.3306,-577.116C76.6326,-566.956 63.8819,-553.696 52.4301,-542.844\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"54.499,-539.991 44.7751,-535.784 49.7533,-545.137 54.499,-539.991\"/>\n",
       "</g>\n",
       "<!-- a_lo -->\n",
       "<g id=\"node2\" class=\"node\"><title>a_lo</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-522\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-518.3\" font-family=\"Times,serif\" font-size=\"14.00\">a_lo</text>\n",
       "</g>\n",
       "<!-- _lon -->\n",
       "<g id=\"node9\" class=\"node\"><title>_lon</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-450\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-446.3\" font-family=\"Times,serif\" font-size=\"14.00\">_lon</text>\n",
       "</g>\n",
       "<!-- a_lo&#45;&gt;_lon -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>a_lo&#45;&gt;_lon</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-503.697C99,-495.983 99,-486.712 99,-478.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-478.104 99,-468.104 95.5001,-478.104 102.5,-478.104\"/>\n",
       "</g>\n",
       "<!-- long -->\n",
       "<g id=\"node3\" class=\"node\"><title>long</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-378\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-374.3\" font-family=\"Times,serif\" font-size=\"14.00\">long</text>\n",
       "</g>\n",
       "<!-- ong_ -->\n",
       "<g id=\"node4\" class=\"node\"><title>ong_</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-306\" rx=\"28.3955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-302.3\" font-family=\"Times,serif\" font-size=\"14.00\">ong_</text>\n",
       "</g>\n",
       "<!-- long&#45;&gt;ong_ -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>long&#45;&gt;ong_</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M87.5599,-361.465C85.5404,-352.946 85.0058,-342.171 85.9562,-332.458\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"89.4264,-332.923 87.5702,-322.492 82.5165,-331.804 89.4264,-332.923\"/>\n",
       "</g>\n",
       "<!-- long&#45;&gt;ong_ -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>long&#45;&gt;ong_</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-359.697C99,-351.983 99,-342.712 99,-334.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-334.104 99,-324.104 95.5001,-334.104 102.5,-334.104\"/>\n",
       "</g>\n",
       "<!-- long&#45;&gt;ong_ -->\n",
       "<g id=\"edge13\" class=\"edge\"><title>long&#45;&gt;ong_</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M110.44,-361.465C112.46,-352.946 112.994,-342.171 112.044,-332.458\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"115.484,-331.804 110.43,-322.492 108.574,-332.923 115.484,-331.804\"/>\n",
       "</g>\n",
       "<!-- ong_&#45;&gt;ng_l -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>ong_&#45;&gt;ng_l</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M106.308,-323.757C113.095,-333.555 121.881,-346.604 126,-360 149.514,-436.466 151.412,-464.144 126,-540 122.576,-550.222 116.315,-560.048 110.187,-568.404\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"107.265,-566.464 103.972,-576.534 112.826,-570.715 107.265,-566.464\"/>\n",
       "</g>\n",
       "<!-- ong_&#45;&gt;ng_l -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>ong_&#45;&gt;ng_l</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M118.337,-319.766C128.49,-329.834 139.292,-344.69 144,-360 167.514,-436.466 169.412,-464.144 144,-540 139.724,-552.765 131.025,-564.913 121.974,-574.351\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"119.52,-571.854 114.718,-581.298 124.361,-576.91 119.52,-571.854\"/>\n",
       "</g>\n",
       "<!-- ng_t -->\n",
       "<g id=\"node7\" class=\"node\"><title>ng_t</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">ng_t</text>\n",
       "</g>\n",
       "<!-- ong_&#45;&gt;ng_t -->\n",
       "<g id=\"edge14\" class=\"edge\"><title>ong_&#45;&gt;ng_t</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-287.697C99,-279.983 99,-270.712 99,-262.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-262.104 99,-252.104 95.5001,-262.104 102.5,-262.104\"/>\n",
       "</g>\n",
       "<!-- g_ti -->\n",
       "<g id=\"node5\" class=\"node\"><title>g_ti</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">g_ti</text>\n",
       "</g>\n",
       "<!-- _tim -->\n",
       "<g id=\"node10\" class=\"node\"><title>_tim</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"26.949\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">_tim</text>\n",
       "</g>\n",
       "<!-- g_ti&#45;&gt;_tim -->\n",
       "<g id=\"edge16\" class=\"edge\"><title>g_ti&#45;&gt;_tim</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-143.697C99,-135.983 99,-126.712 99,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-118.104 99,-108.104 95.5001,-118.104 102.5,-118.104\"/>\n",
       "</g>\n",
       "<!-- time -->\n",
       "<g id=\"node6\" class=\"node\"><title>time</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">time</text>\n",
       "</g>\n",
       "<!-- ng_t&#45;&gt;g_ti -->\n",
       "<g id=\"edge15\" class=\"edge\"><title>ng_t&#45;&gt;g_ti</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-215.697C99,-207.983 99,-198.712 99,-190.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-190.104 99,-180.104 95.5001,-190.104 102.5,-190.104\"/>\n",
       "</g>\n",
       "<!-- g_lo&#45;&gt;_lon -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>g_lo&#45;&gt;_lon</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M37.5771,-505.116C47.0052,-494.792 60.8009,-481.267 73.0189,-470.321\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"75.53,-472.775 80.7584,-463.56 70.9247,-467.503 75.53,-472.775\"/>\n",
       "</g>\n",
       "<!-- g_lo&#45;&gt;_lon -->\n",
       "<g id=\"edge11\" class=\"edge\"><title>g_lo&#45;&gt;_lon</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.188,-508.485C56.3394,-499.028 70.3383,-485.771 81.3679,-474.372\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"84.1483,-476.525 88.4578,-466.846 79.0531,-471.725 84.1483,-476.525\"/>\n",
       "</g>\n",
       "<!-- _lon&#45;&gt;long -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>_lon&#45;&gt;long</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M87.5599,-433.465C85.5404,-424.946 85.0058,-414.171 85.9562,-404.458\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"89.4264,-404.923 87.5702,-394.492 82.5165,-403.804 89.4264,-404.923\"/>\n",
       "</g>\n",
       "<!-- _lon&#45;&gt;long -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>_lon&#45;&gt;long</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-431.697C99,-423.983 99,-414.712 99,-406.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-406.104 99,-396.104 95.5001,-406.104 102.5,-406.104\"/>\n",
       "</g>\n",
       "<!-- _lon&#45;&gt;long -->\n",
       "<g id=\"edge12\" class=\"edge\"><title>_lon&#45;&gt;long</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M110.44,-433.465C112.46,-424.946 112.994,-414.171 112.044,-404.458\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"115.484,-403.804 110.43,-394.492 108.574,-404.923 115.484,-403.804\"/>\n",
       "</g>\n",
       "<!-- _tim&#45;&gt;time -->\n",
       "<g id=\"edge17\" class=\"edge\"><title>_tim&#45;&gt;time</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-71.6966C99,-63.9827 99,-54.7125 99,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-46.1043 99,-36.1043 95.5001,-46.1044 102.5,-46.1043\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "genome = \"a_long_long_long_time\"\n",
    "%dotstr view_de_bruijn(genome, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It did NOT over-collapse the repeats `\"long_long_long_\"`, that is the advatage over SCS-like algorithms. But let's not be too optimistic to think that we have solved all our problems.\n",
    "\n",
    "Let's see another example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "genome = \"ZABCDABEFABY\"\n",
    "k = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.30.1 (20150305.2338)\n",
       " -->\n",
       "<!-- Title: DeBruijn Graph Pages: 1 -->\n",
       "<svg width=\"189pt\" height=\"476pt\"\n",
       " viewBox=\"0.00 0.00 189.00 476.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 472)\">\n",
       "<title>DeBruijn Graph</title>\n",
       "<polygon fill=\"white\" stroke=\"white\" points=\"-4,5 -4,-472 186,-472 186,5 -4,5\"/>\n",
       "<!-- BE -->\n",
       "<g id=\"node1\" class=\"node\"><title>BE</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"71\" cy=\"-450\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"71\" y=\"-446.3\" font-family=\"Times,serif\" font-size=\"14.00\">BE</text>\n",
       "</g>\n",
       "<!-- EF -->\n",
       "<g id=\"node4\" class=\"node\"><title>EF</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"44\" cy=\"-378\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"44\" y=\"-374.3\" font-family=\"Times,serif\" font-size=\"14.00\">EF</text>\n",
       "</g>\n",
       "<!-- BE&#45;&gt;EF -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>BE&#45;&gt;EF</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M64.6014,-432.411C61.4864,-424.335 57.6663,-414.431 54.1654,-405.355\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"57.4045,-404.027 50.5403,-395.956 50.8735,-406.546 57.4045,-404.027\"/>\n",
       "</g>\n",
       "<!-- AB -->\n",
       "<g id=\"node2\" class=\"node\"><title>AB</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">AB</text>\n",
       "</g>\n",
       "<!-- BC -->\n",
       "<g id=\"node3\" class=\"node\"><title>BC</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">BC</text>\n",
       "</g>\n",
       "<!-- AB&#45;&gt;BC -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>AB&#45;&gt;BC</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M84.4297,-218.834C74.2501,-208.938 60.4761,-195.546 48.9694,-184.359\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"51.4055,-181.846 41.7957,-177.385 46.5259,-186.865 51.4055,-181.846\"/>\n",
       "</g>\n",
       "<!-- AB&#45;&gt;BE -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>AB&#45;&gt;BE</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M97.2361,-252.43C94.1769,-282.01 87.4923,-343.904 80,-396 78.7856,-404.444 77.291,-413.615 75.8711,-421.917\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"72.4093,-421.396 74.1393,-431.849 79.3052,-422.599 72.4093,-421.396\"/>\n",
       "</g>\n",
       "<!-- BY -->\n",
       "<g id=\"node9\" class=\"node\"><title>BY</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">BY</text>\n",
       "</g>\n",
       "<!-- AB&#45;&gt;BY -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>AB&#45;&gt;BY</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-215.697C99,-207.983 99,-198.712 99,-190.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-190.104 99,-180.104 95.5001,-190.104 102.5,-190.104\"/>\n",
       "</g>\n",
       "<!-- CD -->\n",
       "<g id=\"node6\" class=\"node\"><title>CD</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">CD</text>\n",
       "</g>\n",
       "<!-- BC&#45;&gt;CD -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>BC&#45;&gt;CD</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M35.3496,-144.765C39.7115,-136.283 45.1469,-125.714 50.0413,-116.197\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"53.2346,-117.641 54.6957,-107.147 47.0096,-114.439 53.2346,-117.641\"/>\n",
       "</g>\n",
       "<!-- FA -->\n",
       "<g id=\"node7\" class=\"node\"><title>FA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"44\" cy=\"-306\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"44\" y=\"-302.3\" font-family=\"Times,serif\" font-size=\"14.00\">FA</text>\n",
       "</g>\n",
       "<!-- EF&#45;&gt;FA -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>EF&#45;&gt;FA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M44,-359.697C44,-351.983 44,-342.712 44,-334.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"47.5001,-334.104 44,-324.104 40.5001,-334.104 47.5001,-334.104\"/>\n",
       "</g>\n",
       "<!-- ZA -->\n",
       "<g id=\"node5\" class=\"node\"><title>ZA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"154\" cy=\"-306\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"154\" y=\"-302.3\" font-family=\"Times,serif\" font-size=\"14.00\">ZA</text>\n",
       "</g>\n",
       "<!-- ZA&#45;&gt;AB -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>ZA&#45;&gt;AB</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M142.066,-289.811C134.786,-280.546 125.337,-268.52 117.146,-258.094\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"119.841,-255.86 110.911,-250.159 114.337,-260.185 119.841,-255.86\"/>\n",
       "</g>\n",
       "<!-- DA -->\n",
       "<g id=\"node8\" class=\"node\"><title>DA</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"108\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"108\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">DA</text>\n",
       "</g>\n",
       "<!-- CD&#45;&gt;DA -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>CD&#45;&gt;DA</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M73.2112,-73.1159C78.8497,-64.345 85.9734,-53.2637 92.3038,-43.4162\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"95.3657,-45.1257 97.8292,-34.8212 89.4775,-41.3404 95.3657,-45.1257\"/>\n",
       "</g>\n",
       "<!-- FA&#45;&gt;AB -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>FA&#45;&gt;AB</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M55.9342,-289.811C63.2136,-280.546 72.6626,-268.52 80.8544,-258.094\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"83.663,-260.185 87.0891,-250.159 78.1587,-255.86 83.663,-260.185\"/>\n",
       "</g>\n",
       "<!-- DA&#45;&gt;AB -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>DA&#45;&gt;AB</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M115.615,-35.3225C128.288,-64.7748 150.828,-128.529 135,-180 131.714,-190.685 125.46,-201.149 119.075,-209.919\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"116.116,-208.023 112.75,-218.069 121.646,-212.314 116.116,-208.023\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%dotstr view_de_bruijn(genome, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above De-Bruijn graph, the Eulerian walk is not unique. Thus two possible assembly can be generated like:\n",
    "- ZAB<font color=\"blue\">CD</font>AB<font color=\"magenta\">EF</font>ABY \n",
    "\n",
    "- ZAB<font color=\"magenta\">EF</font>AB<font color=\"blue\">CD</font>ABY\n",
    "\n",
    "Here the repeats \"AB\" is not tandem repeats, instead they are anologous the **transposon elements** in the genome. Therefore, you cannot determine whether it is *A* or *B* in the above example.\n",
    "\n",
    "For all the possible Eulerian paths, only one of them is the correct assembly. All the rests correspond to incorrect reshufflings that portions of the genome that occur between repetitive elements are going to be shuffled around and put in the wrong order.\n",
    "\n",
    "So, again, the repetitive elements make the assembly difficult, specifically the transposon-like repetitive elements.\n",
    "\n",
    "This will also introduce ambiguity, as will be done by the SCS-based algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ZABCDABEFABY'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reads = [\"ZAB\", \"ABC\", \"BCD\", \"CDA\", \"DAB\", \"ABE\", \"BEF\", \"EFA\", \"FAB\", \"ABY\"]\n",
    "scs(reads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ZABEFABCDABY'"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "reads = [\"ZAB\", \"ABC\", \"BCD\", \"CDA\", \"DAB\", \"ABE\", \"BEF\", \"EFA\", \"FAB\", \"ABY\"]\n",
    "reads2 = list(np.random.permutation(reads))\n",
    "scs(reads2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\">$\\S$ EXERCISE 3</red>\n",
    "\n",
    "Write a python function `eulerian_walk(nodes, edges)` to execute the Eulerian walk on the De-Bruijn graph and test it on the above genome sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class DeBruijnGraph():\n",
    "    \n",
    "    def __init__(self, reads, k):\n",
    "        \"\"\"build de bruijn graph from a set of reads using k-mer \n",
    "        \"\"\"\n",
    "        pass\n",
    "    \n",
    "    def eulerian_walk(self):\n",
    "        \"\"\"infer the eulerian path from the de bruijn graph, get the assembly result\n",
    "        \"\"\"\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many of the above implementations are based on the assumption that\n",
    "> None of the $k$-mer of the genome is missed by the sequencing.\n",
    "\n",
    "But in practical, this is seldomly true, let alone the prevalent sequencing errors.\n",
    "\n",
    "And the built graph are never Eulerian, in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Summary\n",
    "\n",
    "Now we know that the two categories of computational methods to implement de-novo genome are flawed. But in the long term, they are stilll useful in doing such tasks. \n",
    "\n",
    "The OLC approach are based upon the **overlap graph** which are built upon directly from the reads, while the DBG methods upon the  **De-Bruijn graph** from the $k$-mers of the reads.\n",
    "\n",
    "Once we have built the graph, the problem is still far from resolved, since the graph is big and messy.\n",
    "\n",
    "To reach at the optimal assembly, the former methods are relying on the **shortest-common-superstring (SCS)** assumption, while the latter dependent upon the **Eulerian path**.\n",
    "\n",
    "- The slow (brute-force) SCS-finding algorithm can surely find the shortest superstring to contain all the reads, but this is at the cost of huge computational burdens. Moreover, it cannot avoid the **over-collapsing problem of repeats**.\n",
    "\n",
    "- The greedy SCS-finding approach is faster, but this heuristic algorithm is not guaranteed to finding the optimal solution.\n",
    "\n",
    "- The De-Bruijn graph (DBG)-based algorithm can somewhat to some extent avoid the problems caused by repeats, but not completely.\n",
    "\n",
    "- The OLC-based approach is appropriate for the assembly tasks of longer reads, while the De-Bruijn Graph algorithm, for shorter reads. \n",
    "\n",
    "\n",
    "So, if the genome sequence is repetitive, no matter which agorithms we pick, we can hardly avoid the mistake. This can only be resolved by the sequencing technologies of longer reads because they can span repetitive portions of the genome and the surrounding non-repetitive parts.\n",
    "\n",
    "Instead of generating a whole genome sequence, the de-novo assembler tends to report a set of contigs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Stop and think\n",
    "\n",
    "And finally, can you summarize the possible reasons that the graph, either overlap graph or De-Bruijn graph, is messy?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}