Week6_debruijn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 6 lab - de novo assembly with De Bruijn graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this week's lab you'll be writing programs to count k-mers and to find the edges of a de Bruijn graph.\n",
    "\n",
    "As a reminder, we have a relevant tutorial document intended for those new to programming or learning Python: [kmer-counting tutorial](kmer_counting.ipynb). The first couple of challenges below are covered (more gently) in this tutorial. By now you have probably already learned the programming concepts in the tutorial, but you might still find it worth looking at the _Dictionaries: algorithmic considerations_ section. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll define some toy \"reads\" from English strings, which you can use to see your functions in action."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ters_no',\n",
       " 'ers_not',\n",
       " 'matners',\n",
       " 'size_ma',\n",
       " 'tters_',\n",
       " 'atters_n',\n",
       " 'size_m',\n",
       " 'e_matter']"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yoda_reads = '''\n",
    "ters_no\n",
    "ers_not\n",
    "matners\n",
    "size_ma\n",
    "tters_\n",
    "atters_n\n",
    "size_m\n",
    "e_matter'''.strip().split(\"\\n\")\n",
    "yoda_reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "caesar_reads = '''\n",
    "_your_e\n",
    "l,nd_me\n",
    "nd_me_y\n",
    "d_me_y\n",
    "omans,_\n",
    "s,_coun\n",
    "me_you\n",
    "eids,_\n",
    "mans,_c\n",
    "our_ea\n",
    "'''.strip().split(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a version with no read errors\n",
    "caesar_reads_perfect = '''\n",
    "me_your_e\n",
    "ymen,_len\n",
    "_your_e\n",
    "ds,_Roma\n",
    "nds,_Roma\n",
    "trymen,_\n",
    "lend_me_\n",
    "omans,_co\n",
    "n,_lend_\n",
    "riends,_R\n",
    "ntrymen\n",
    ",_country\n",
    "e_your_ears\n",
    "untryme\n",
    "riends,_\n",
    "'''.strip().split(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And some error-free \"reads\" from a very small tRNA gene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt_te_reads = '''TGTAGTTGAAATACAACGAT\n",
    "GGTCGTGGTTGTAGTCCGTGC\n",
    "TGGTCGTGGTTGTAGTCCG\n",
    "TATCATTGGTCGTGGTTGTAG\n",
    "TTGTAGTTGAAATACAACGAT\n",
    "ATCATTGGTCGTGGTTGTAG\n",
    "ATGGTTTTTCATATCATTGG\n",
    "GGTTTTTCATATCATTGGTCGT\n",
    "GGTCGTGGTTGTAGTCCGT\n",
    "TCATTGGTCGTGGTTGTAGTCC\n",
    "GGTCGTGGTTGTAGTCCGTGCG\n",
    "GATGGTTTTTCATATCATT\n",
    "TGGTCGTGGTTGTAGTCC\n",
    "AACGATGGTTTTTCATATCA\n",
    "GGTTTTTCATATCATTGGTCG\n",
    "GTTTTTCATATCATTGGTCGTG\n",
    "ATTGGTCGTGGTTGTAGTCCGT\n",
    "AACGATGGTTTTTCATAT\n",
    "GTTTTTCATATCATTGGT\n",
    "AACGATGGTTTTTCATAT'''.strip().split()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Counting k-mers "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1:**\n",
    "\n",
    "Complete the function below to return all the k-mers from a given string. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_string_kmers(read, k):\n",
    "    '''\n",
    "    Return a set of all k-mers of length k from string read.\n",
    "    '''\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return {'_yo', 'our', 'r_e', 'ur_', 'you'}\n",
    "get_string_kmers(caesar_reads[0], 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:**\n",
    "\n",
    "Write a function which applies `get_string_kmers()` to a list of strings and returns all kmers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kmers(reads, k):\n",
    "    '''\n",
    "    Given a list of strings representing reads, and a value k, \n",
    "    return a set of all k-mers of length k.\n",
    "    '''\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return {'tters_','ize_ma','atners','ers_no','matner','rs_not','size_m',\n",
    "#                '_matte','e_matt','atters','ters_n','matter'}\n",
    "# (i.e. 12 unique kmers)\n",
    "get_kmers(yoda_reads, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a graph with networkx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could represent a graph using Python data structures. For instance, here is one way to represent a graph where nodes B and C are linked from A:\n",
    "    \n",
    "    B <- A -> C\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = ['A','B','C']\n",
    "edges = [('A','B'), ('A','C')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that a tree is a kind of graph, and we have already built trees using data structures in past labs.\n",
    "\n",
    "This time we'll use a library `networkx` intended specifically for graph manipulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how to build that same graph in networkx. Note that we use `DiGraph` instead of `Graph`, which gives a directed graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = nx.DiGraph()\n",
    "g.add_edge('A','B')\n",
    "g.add_edge('A','C')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have explicitly added nodes with code like `g.add_node('A')`, but since every node is connected to at least one edge in this case, networkx automatically adds the nodes for us when the edges are added. \n",
    "\n",
    "Remember you can look at networkx method documentation with `?` or `help()`. Have a look at the documentation for the `add_edge` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NodeView(('A', 'B', 'C'))"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.nodes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OutEdgeView([('A', 'B'), ('A', 'C')])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.edges()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a way to draw the graph with a spring layout, where networkx will try to place nodes so they are not too close together. This is a simple layout algorithm - it's easy to see the result for this tiny graph, but may be difficult see for graphs of any size. The layout is somewhat random and for larger graphs, will be different every time you draw it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFCCAYAAADGwmVOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAD9xJREFUeJzt3VGIneldx/HfyZrdmQEZghuV2nTaSrHiZiI0U7W4mrawNXuVFLVguyWLvZkVRZogZL0oFEJ6MeneCGPLgqNgwUJgm5QNC73ZKiidlsREWitLqWGh21klDZZMOtnkePEmm8nMnMmZmXPO8573fD4Qwpzzvi/P3Z/vOc9531a73W4HABi4XaUXAACjyhAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAoxBAGgEIMYQAo5OdKLwAAtmRpKVlYSC5fTq5fTyYnk+np5Nlnk717S69uS1rtdrtdehEA8FCLi8np08mFC9XfN2/ef298PGm3k8OHk5Mnk5mZMmvcIkMYgPqbn09OnEiWl6th20mrVQ3kublkdnZw69smH0cDUG/3BvCNGw8/tt2ujjtxovq75oNYCQNQX4uLyaFDDwzgdyf5cZJHkuxO8qEkf5tk39pzJyaSV19NDh4cyFK3w+5oAOrr9OnqI+g1zif5aZIfJfmlJH++0bnLy9X5NWYIA1BPS0vVJqxNPrAdS/KHSb670ZvtdvLyy8mbb/ZpgTtnCANQTwsLDz3kRpJ/SvLbnQ5otbq6Tik2ZgFQT5cvP/gzpFWOpBpgP03yi0le6XSN5eXkypW+LK8XlDAA9XT9ese3XkrykyQ/S/I3SX4/yRudDr52rdcr6xlDGIB6mpx86CGPJPn43f//pdNBe/b0bk09ZggDUE/T08nY2KaHtJN8Lcm1JL++0QHj48n+/b1fW4/4nTAA9bS0lExNrfte+N25/zvhVpKpJCeTfHKja4yNJVev1vae0jZmAVArb731Vl5//fV861vfypGPfjSPvvzyAz9T+mG3F2q1kqefru0ATnwcDUANrKys5Mknn8zjjz+esbGxvO9978snPvGJ/NuHP1x9pLwd4+PVwxxqzBAGoLjdu3fnzp07uXbtWm7fvp3bt2/nIx/5SH7v+PHqYQwTE1u74MREdV6Nb1mZGMIA1MCtW7cyMzOTO3fuJEkee+yxvPjii9Wbs7P3B3GrtfmFWq37A7jmD29IDGEACrt06VI++MEP5rXXXssXvvCFJMlzzz2X97znPfcPmp2tHsZw9Gi12WrtR9Tj49XrR49Wxw3BAE5szAKgkJWVlZw6dSrz8/OZm5vLM888kyR59NFH85nPfGb9CQcPJmfPVveCXlio7oR17Vr1O+D9+5Njx2q9CWsjfqIEwMBdunQpx44dyzvf+c58+ctfzjve8Y7SSyrCx9EADMzKyko+97nP5amnnspnP/vZnD9/fmQHcOLjaAAGZHX9Xrp0aaSH7z1KGIC+Ur+dKWEA+kb9bk4JA9Bz6rc7ShiAnlK/3VPCAPSE+t06JQzAjqnf7VHCAGyb+t0ZJQzAtqjfnVPCAGyJ+u0dJQxA19RvbylhAB5K/faHEgZgU+q3f5QwABtSv/2nhAFYR/0OhhIG4G3qd7CUMABJ1G8JShhgxKnfcpQwwAhTv2UpYYARpH7rQQkDjBj1Wx9KGGBEqN/6UcIAI0D91pMSBmgw9VtvShigodRv/SlhgIZRv8NDCQM0iPodLkoYoAHU73BSwgBDTv0OLyUMMKTU7/BTwgBDSP02gxIGGCLqt1mUMMCQUL/No4QBak79NpcSBqgx9dtsShightTvaFDCADWjfkeHEgaoCfU7epQwQA2o39GkhAEKUr+jTQkDFKJ+UcIAA6Z+uUcJAwyQ+mU1JQwwAOqXjShhgD5Tv3SihAH6RP3yMEoYoA/UL91QwgA9pH7ZCiUM0CPql61SwgA7pH7ZLiUMsAMXL17MsWPHsm/fPvXLlilhgG24V78f+9jHcvz4cfXLtihhgC1Sv/SKEgbokvql15QwQBfUL/2ghAE2oX7pJyUM0IH6pd+UMMAa6pdBUcIAq6hfBkkJA0T9UoYSBkae+qUUJQyMLPVLaUoYGEnqlzpQwsBIUb/UiRIGRob6pW6UMNB46pe6UsJAo6lf6kwJA42kfhkGShhoHPXLsFDCQGOoX4aNEgYaQf0yjJQwMNTUL8NMCQNDS/0y7JQwMHTUL02hhIGhon5pEiUMDAX1SxMpYaD21C9NpYSB2lK/NJ0SBmrpXv2+613vUr80lhIGamV1/Z44cSLnzp0zgGksJQzUhvpl1ChhoDj1y6hSwkBR6pdRpoSBItQvKGGgAPULFSUMDIz6hQcpYWAg1C+sp4SBvlK/0JkSBvpG/cLmlDDQc+oXuqOEgZ5Sv9A9JQz0hPqFrVPCwI6pX9geJQxsm/qFnVHCwLaoX9g5JQxsifqF3lHCQNfUL/SWEgYeSv1CfyhhGBVLS8nCQnL5cnL9ejI5mUxPJ88+m+zd2/E09Qv902q32+3SiwD6aHExOX06uXCh+vvmzfvvjY8n7XZy+HBy8mQyM/P2WysrKzl16lTm5+dz5syZfOpTn0qr1Rrw4qHZDGFosvn55MSJZHm5GradtFrVQJ6bS2ZnH6jfL33pS+oX+sQQhqa6N4Bv3Oj6lPbERL5+6FD+dHFR/cIAGMLQRIuLyaFDGw7gQ0n+PckbSR7b4NSbjzyS/zt/PnsPH+7rEgG7o6GZTp+uPoJe44dJ/jlJK8m5Dqc+dudO9r74Yv/WBrxNCUPTLC0lU1MPbsC66/NJXknyW0n+K8nXO11jbCy5enXTXdPAzilhaJqFhY5v/UOST97990qSH3c6sNXa9DpAbxjC0DSXL29Ywf+S5L+T/HGSDyT51SRf6XSN5eXkypV+rRC4yxCGprl+fcOX/z7JU0kev/v3n9x9raNr13q6LGA9d8yCppmcXPfScpKvJrmd5JfvvvazJD9JtVP6wEbX2bOnP+sD3qaEoWmmp6uNVau8lOSRJN9Ncunuv+8leTLV98TrjI8n+/f3d52A3dHQOBvsjv6DJL+R5MyaQ7+a5C+SvJ41H4vZHQ0DYQhDw6ysrOS16en82ve/n0e2c4FWKzl6NDl7ttdLA9bwcTQ0yMWLFzMzM5MX9+5Na3x8excZH68e5gD0nSEMDbD2eb9nvvnN7DpzJpmY2NqFJiaqhzgcPNifhQIPsDsahlzH5/3Ozlb/b+MpSsBgKGEYUmvr99y5c+sfOTg7m7z6avUd79hYNWhXGx+vXj96tDrOAIaBsjELhtC2nvf75pvVrSivXKluxLFnT/UzpGPH7IKGQgxhGCIrKys5depU5ufnPe8XGsB3wjAkOn73Cwwt3wlDzXX13S8wlJQw1Jj6hWZTwlBD6hdGgxKGmlG/MDqUMNSE+oXRo4ShBtQvjCYlDAWpXxhtShgKUb+AEoYBU7/APUoYBkj9AqspYRgA9QtsRAlDn6lfoBMlDH2ifoGHUcLQB+oX6IYShh5Sv8BWKGHoEfULbJUShh1Sv8B2KWHYAfUL7IQShm1Qv0AvKGHYIvUL9IoShi6pX6DXlDB0Qf0C/aCEYRPqF+gnJQwdqF+g35QwrKF+gUFRwrCK+gUGSQlD1C9QhhJm5KlfoBQlzMhSv0BpSpiRpH6BOlDCjBT1C9SJEmZkqF+gbpQwjad+gbpSwjSa+gXqTAnTSOoXGAZKmMZRv8CwUMI0hvoFho0SphHULzCMlDBDTf0Cw0wJM7TULzDslDBDR/0CTaGEGSrqF2gSJcxQUL9AEylhak/9Ak2lhKkt9Qs0nRKmltQvMAqUMLWyun6PHz+ufoFGU8LUxr363bdvn/oFRoISpri19Xv+/HkDGBgJSpii1C8wypQwRahfACVMAeoXoKKEGRj1C/AgJcxAqF+A9ZQwfaV+ATpTwvSN+gXYnBKm59QvQHeUMD2lfgG6p4TpCfULsHVKmB1TvwDbo4TZNvULsDNKmG1RvwA7p4TZEvUL0DtKmK6pX4DeUsI8lPoF6A8lzKbUL0D/KGE2pH4B+k8Js476BRgMJczb1C/AYClhkqhfgBKU8IhTvwDlKOERpn4BylLCI0j9AtSDEh4x6hegPpTwiFC/APWjhEeA+gWoJyXcYOoXoN6UcEOpX4D6U8INo34BhocSbhD1CzBclHADqF+A4aSEh5z6BRheSnhIqV+A4aeEh5D6BWgGJTxE1C9AsyjhIaF+AZpHCdec+gVoLiVcY+oXoNmUcA2pX4DRoIRrRv0CjA4lXBPqF2D0KOEaUL8Ao0kJF6R+AUabEi5E/QKghAdM/QJwjxIeIPULwGpKeADULwAbMYT74JVXXskTTzyRGzdu5OLFi5mZmcl3vvOdXLp0KZ/+9KfTarVKLxGAGmi12+126UXU3tJSsrCQXL6cXL+eTE4m09PJs88me/c+cOitW7fy3ve+N2+88UYOHDiQq1evZm5uLs8884zhC8ADDOHNLC4mp08nFy5Uf9+8ef+98fGk3U4OH05OnkxmZpIkL7zwQp5//vncvHkzu3btytmzZ3PkyJECiweg7gzhTubnkxMnkuXlath20mpVA3luLj86ciRTU1O5detWkmTXrl2ZmprKD37wgwEtGoBhYnf0Ru4N4Bs3Hn5su10dd+JE/uMb38ju3btz4MCBPPHEE3n/+9+fAwcO9H+9AAwlJbzW4mJy6NC6AfyVJF9M8p9Jfj7Jbyb56yS/u/qgiYnk1VeTgwcHs1YAhprd0WudPl19BL3KF5P8ZZLnk/w4ydUkzyX52tpzl5er8wGgC0p4taWlZGrqgQ1Y15P8SpK/S/JH3VxjbCy5enXdrmkAWEsJr7awsO6lf01yM8nRbq/Ram14HQBYyxBe7fLlB3+GlOR/kzyeLexgW15Orlzp8cIAaCJDeLXr19e99AtJ/ifJW1u5zrVrPVoQAE1mCK82Obnupd9JMpbkpa1cZ8+eHi0IgCYzhFebnq42Vq0ymeTzSf4s1SC+keRWkgtJ/mqja4yPJ/v393edADSC3dGrbbA7+p5/TPJCku+l+p3wB1L9TvhDaw+0OxqALhnCa33848lLL21+q8pOWq3k6NHk7NnerwuAxjGE1+pwx6yuuGMWAFvgO+G1ZmaSublqoG7FxER1ngEMQJc8wGEjs7PV/1t8itLb5wFAF3wcvZlvf7u6F/TLL1fDdvU9pe89T/jpp6vnCStgALbIEO7Gm29Wt6K8cqW6EceePdXPkI4dswsagG0zhAGgEBuzAKAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACjGEAaAQQxgACvl/K+q0DivDn9IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx.draw_spring(g, with_labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFCCAYAAADGwmVOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEhxJREFUeJzt3W+sZHV9x/Hvmbmzs8vdZSWI2GWXdmM1IiSYgKI+sfpgSVuIVKtGBEPVR7aNZSEGbSJJ4z4wSn1gddPE2KARaNxE/gRIY9KkVuPGNGFLSVGjCWVtK43KunD37uydO6cPDgPLcvfuvfPn/M6f1yu5uSrMzPHR577POTOT5XmeBwBQuk7qAwCAtjLCAJCIEQaARIwwACRihAEgESMMAIkYYQBIxAgDQCJGGAASMcIAkIgRBoBEjDAAJGKEASARIwwAiRhhAEjECANAIgupDwAANiPPI1ZXi995HpFlxU+3W/yuEyMMQKXlecRwGDEYFL/z/Oz/bpZFLCxE9PvF76qPcpbn6/3fAYA08rwY3sEgYjTa/OM7nWKM+/3qjrERBqByhsOIpaXJxvdMnU7E4mJRxlVjhAGojDyPWF4u6nfW+v2IbduqVcVGGIBKyPOifldW5vcavV5RxVUZYm9RAiC5MgY4onj+paX1b+4qkxEGILnl5fkP8NjKSvF6VWCEAUhq/PajMo3f7pSaEQYgmfFp6BSqcFraCAOQzKTvAZ6F0aj8Aj+TEQYgifGHcaQ0GKStYSMMQBLDYboKHhuN0l4bNsIAJDFNBR86dE+8+91Xx5492+Oyy34n3v/+P4zDh79f+nFMywgDkMSkBfqVr/xtfOYzfxW33vqZ+PGPn4nHH386PvaxT8Qjjzww0fOtrk52HLPgE7MAKF2eRxw7tvnHHT/+27j88kviy1/+h7jhhvfP7Hhe9ao0n6KlhAEo3aT1+aMf/TBOnjwZ1133J5U4nmkZYQBKN+k52Gef/XVceOGrY2HGX4mU6pywEQagdJOO3gUXXBi//vWvYjjjW5qNMACtMen117e+9e2xdevWePjh+ytxPNMywgCUbtLRO//8nXHHHX8Tn/rUn8fDD98fJ06ciJWVlfjudx+NO+/8VOnHMy13RwNQuknvjh779re/FQcPfil++tMnY/v2HXHllVfF/v1/Hddc846Jni/V3dFGGIAkjh1L/wUKERGdTsTOnYleO83LAtB2M77BeWLdbrrXNsIAJNHvpz6CQsrjMMIAJLGwUJwKTqnTSVvkRhiAJLIsfQ33++nujI6IqMgZeQBaYTSKeOaZiKeeivjZz6L/L9+LwRf/PkZ5+U3Y6aT/I8Dd0QCU4/HHI97ylpfOAS8tRSwsxPDZ4/Hcqa2lH86OHelvDnM6GoByvOlNEa97XcTJkxHPPx+xZUvE178eC4tbSy/Sfj/9AEcYYQDK8sQTxeno8d1Yr399xIc/HBER27ZF9HrlHEavV7xeFRhhAObr1KmIO++M2Lcv4tOfLv5zRMTXvvbiXVFZFrG4OP8h7vWK10l5M9bpKhDjADTWkSMRt9wSsXt3xGOPRVxyScRwGPHOd0Zcc83L/tXxEC8vRwwGsz+Ufr8o4KoMcIQbswCYh1OnIg4ciDh4MOILX4j4yEc2tX7DYXHf1mg0/aF0OsW4V+Ea8JkqeEgA1Npa9btJCwsR559fFPFgMNkYj9+ClPq9wOtRwgDMxpT1ezZ5XpTxYBCxurr+IHc6xWdBj+9+rur4jilhAKY3g/o9mywrbqga37SV58UY53nxk2XFT7db/dE9k7ujAZjc6Xc+33prxEMPzXSA15JlReX2esVbjXu9elTvWpQwAJOZY/22hRIGYHMS1G9TKWEANk79zpQSBuDc1O9cKGEA1qd+50YJA7A29Tt3ShiAV1K/pVDCALxE/ZZKCQNQUL+lU8IAbad+k1HCAG2mfpNSwgBtpH4rQQkDtI36rQwlDNAW6rdylDBAG6jfSlLCAE2mfitNCQM0lfqtPCUM0DTqtzaUMECTqN9aUcIATaB+a0kJA9Sd+q0tJQxQV+q39pQwQB2p30ZQwgB1on4bRQkD1IX6bRwlDFB16rexlDBAlanfRlPCAFWkfltBCQNUjfptDSUMUBXqt3WUMEAVqN9WUsIAKanfVlPCAKmo39ZTwgBlU7+8QAkDlEn9cholDFAG9csalDDAvKlfzkIJA8yL+uUclDDAPKhfNkAJA8yS+mUTlDDArKhfNkkJA0xL/TIhJQwwDfXLFJQwwCTULzOghAE2S/0yI0oYYKNOr9/9+9UvU1PCABtxev0eORKxa1fqI6IBlDDAetaqXwPMjChhgLNRv8yZEgY4k/qlJEoY4HTqlxIpYYAI9UsSShhA/ZKIEgbaS/2SmBIG2kn9UgFKGGgX9UuFKGGgPdQvFaOEgeZTv1SUEgaaTf1SYUoYaCb1Sw0oYaB51C81oYSB5lC/1IwSBprhsceK+t2zR/1SG0oYqLdx/V57bcRtt6lfakUJA/Wlfqk5JQzUj/qlIZQwUC/qlwZRwkA9qF8aSAkD1ad+aSglDFSX+qXhlDBQTeqXFlDCQLWoX1pECQPVoX5pGSUMpKd+aSklDKSlfmkxJQykoX5BCQMJqF+ICCUMlEn9wssoYaAc6hdeQQkD86V+4ayUMDA/6hfWpYShJfI8YjiMWFkp4nRlpfjveT6HFzu9fm+/Xf3CWShhaKjx6A4G5x7bLItYWIjo94vfWTbFC4/r99JL1S+cQ5bnc/k7GEgkz4vhHQwiRqPNP77TKca439/kGJ86FXHgQMTBgxF33RVx001Trjk0nxKGBhkOI5aWJhvfsdEoYnm5GPHFxaKMz0n9wkRcE4YGyPOIEycinntuugE+3WhUPN+JE+ucyj7z2u+DDxpg2AQlDDWX50X9rqzM5/nHp7UXF884u6x+YWpKGGps3gM8trJSvE6eh/qFGVLCUGPLy/Mf4LGVlYjlp56J827Yp35hRtwdDTU1HBbXbEu1vBw7/u2fY+H6P3LnM8yAEYYayvOI48dndxPWZl64083i/PNtMMyCa8JQQ5O+B3hqWRajUfH6wPSMMNTM+MM4UhoM5vRxl9AyRhhqZjhMVMGnGY2K4wCmY4ShZqat4Ouv/4PYu/eCGEz5RKlrHJrACEPNTFOgTz/9VPzwh/8aWZbFo48+ONVxrK5O9XAgjDDUSp5Pdy32vvu+EVdf/bb40Iduifvuu3uqYxmNXBeGaXmLEtTItO8Nvuqq349PfGJ/XHXVNbFv39viiSd+Ea95zcUTP9+OHRv8ggdgTUoYamSaP5kPH/5+HD36X3HDDR+IN7/5qti793Vx6NA9yY4HMMJQK9OM3r333h3vete+uPDCV0dExPved+PUp6SNMEzH6WiokZWViOef3/zjlpeX441vfG2MRquxuLg9IiIGg0H89rfH4nvfOxJXXHHlRMezfXtErzfRQ4HwBQ5QK5N+VOQjj9wf3W43fvCD/4heb8uL//tHP/qBuO++b8TnPndXqccDFJyOhhrpdid73L333h033vhnsXv3pXHxxa998efjH/+LOHToWzGc8H1Pkx4PUHA6Gmrm2LFqXIvtdCJ27kx9FFBvShhqpipvCVLBMD0jDDXT76c+gkJVjgPqzAhDzSwsFKeCU+p0qlPkUGdGGGomy8YVmu7CcL/vzmiYBSMMNdR/8kh0/ue/k3ynYafjVDTMihGGOjl1KuLOOyO7dl8s/vzfIzrl5+jiogqGWXFVB+risccibrkl4tJLI44ciYVdu6J/otzv9e33XQuGWVLCUHUv1G9ce23E7bdHPPhgxK5dERGxbVt5HxvZ6xWvB8yOv2mhys6o3/H4jmVZcXp4aan4XOl56fWchoZ5UMJQRevU75nGQzyvm6X6fQMM86KEoWrOUb9rybKI886L2LKlqOJZ3DTd6RTj6xowzI/PjoaqOHUq4sCBiIMHI+66K+KmmybKzzwvbtYaDCYb4/FbkLwXGObP37hQBRPU79lkWcTWrcWIDofFGK+urj/InU7xWdDju5+NL5RDCUNKM6rfjcjzYozzvPjJsuKn2zW6kIoShlRmWL8bkWWu70LVuDsayraJO5+BZvN3MZSp5PoFqk0JQxnUL7AGJQzzpn6Bs1DCMC/qFzgHJQzzoH6BDVDCMEvqF9gEJQyzon6BTVLCMK3T6/e229QvsGFKGKYxrt89e9QvsGlKGCZxZv0+9JABBjZNCcNmqV9gRpQwbJT6BWZMCcNGqF9gDpQwrEf9AnOkhOFs1C8wZ0oYzqR+gZIoYTid+gVKpIQhQv0CSShhUL9AIkqY9lK/QGJKmHZSv0AFKGHaRf0CFaKEaQ/1C1SMEqb51C9QUUqYZlO/QIUpYZpJ/QI1oIRpHvUL1IQSpjnUL1AzSphmUL9ADSlh6k39AjWmhKkv9QvUnBKmftQv0BBKmHpRv0CDKGHqQf0CDaSEqb4jR4r63b1b/QKNooSprnH97tsXsX+/+gUaRwlTTeoXaAElTLWoX6BFlDDVoX6BllHCpKd+gZZSwqSlfoEWU8KkoX4BlDAJqF+AiFDClEn9AryMEqYc6hfgFZQw86V+Ac5KCTM/6hdgXUqY2VO/ABuihJkt9QuwYUqY2VC/AJumhJme+gWYiBJmcuoXYCpKmMmoX4CpKWE2R/0CzIwSZuPUL8BMKWHOTf0CzIUSZn3qF2BulDBrU78Ac6eEeSX1C1AKJcxL1C9AqZQwBfULUDol3HbqFyAZJdxm6hcgKSXcRuoXoBKUcNuoX4DKUMJtoX4BKkcJt4H6BagkJdxk6heg0pRwU6lfgMpTwk2jfgFqQwk3ifoFqBUl3ATqF6CWlHDdqV+A2lLCdaV+AWpPCdeR+gVoBCVcJ+oXoFGUcF2oX4DGUcJVp34BGksJV5n6BWg0JVxF6hegFZTwBuR5xOpq8TvPI7Ks+Ol2i98zpX4BWsMIryHPI4bDiMGg+J3nZ/93syxiYSGi3y9+TzzKp05FHDgQcfBgxBe/GHHzzXNYeACqxAifJs+L4R0MIkajjT9mZaX46XSKMe73N7mf6heglVwTfsFwGHH8eMTy8sYH+EyjUfH448eL51vT8eMR731vxG9+49ovQMu1voTzvBjOwWB2zzkaRTz3XFHE27adUcWf/WzEAw8UY/yrX6lfgBbL8ny9K57NlucRS0vFqeR56fUiFhdfGOKf/zziiisiTp4s/uEnPxnxpS+59gvQUq0d4TIGeOzFIX7H2yMOH37pH+zcGXH0aMSOHfM/CAAqp7Wno5eXyxngiOJ1lp/+vzjv8OGILVsiLrkkYu/eiMsvL+cAAKikVpbwcFhcsy3bjt7JWNi+tfwXBqCSWnd39Pg0dApLq1vXfc8xAO3SuhHezHuAZ200mu1d2ADUW6tGePxhHCkNBut/AhcA7dGqER4O01Xw2Gi0zgd5ANAqrRrhSSv4yit/L3bt2hZ79myPvXsviA9+8I/jF784WvpxANAsrRrhaQr0nnseiqNHn48nn/zfuOiii+OOO/5y4udaXZ38OABojtaM8PhrCKe1devWeM97/jR+8pP/nPg5RiPXhQFo0QjPqj5PnDgR3/nOP8bVV7+tEscDQH215hOzpi3Pm2++IbrdhVhaej4uuug1cejQPyU9HgDqrzUlPO3offOb98dTTx2LX/5yEJ///N/Fdde9M5555pfJjgeA+mvNCM/qi4q63W5cf/17o9vtxuHD309+PADUV2tOR89q9PI8j0cffTCOHXs23vCGy5IfDwD11ZoR7nane/yNN14f3W43siyL3bt/N7761bvjsssm/xakaY8HgPpr1bcoHTtWjWuxnU7xVcIAtFtrrglHRCxUpPtVMAARLRvhfj/1ERSqchwApNWqEV5YKE4Fp9TpVKfIAUirVSOcZekrtN93ZzQAhVaNcEQxgqlquNNJ/0cAANXRuhHOsojFxTSvvbioggF4SetGOKK4Jlt2kfb7rgUD8HKtHOGIiG3bInq9cl6r1yteDwBO19oRHp+WnvcQ93pOQwOwtlZ9YtZa8jxieTliMJj9c/f7RQEbYADW0voRHhsOI5aWIkaj6Z+r0ynq1zVgANZjhE+T50URDwaTjfH4LUjeCwzARhjhNeR5UcaDQcTq6vqD3OkUnwU9vvvZ+AKwUUZ4A/K8GOM8L36yrPjpdo0uAJMzwgCQSGvfogQAqRlhAEjECANAIkYYABIxwgCQiBEGgESMMAAkYoQBIBEjDACJGGEASMQIA0AiRhgAEjHCAJCIEQaARIwwACRihAEgESMMAIkYYQBIxAgDQCJGGAASMcIAkIgRBoBEjDAAJGKEASARIwwAiRhhAEjECANAIkYYABIxwgCQyP8DSECNJXZooBAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx.draw_spring(g, with_labels=True, node_size=1200, node_color='#eeeeff', edge_color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:**\n",
    "\n",
    "Write functions `get_suffix()` and `get_prefix()` which just return the `k-1`-length strings which are, respectively the suffix and prefix of the supplied k-mer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write both functions here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'GGTA'"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should return 'GGTA'\n",
    "get_suffix('AGGTA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'AGGT'"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should return 'AGGT'\n",
    "get_prefix('AGGTA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a set of reads and a value k, build a directed graph, using networkx, where nodes are the prefixes/suffixes of kmers and an edge exists for every kmer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_graph(reads, k):\n",
    "    '''\n",
    "    Given a set of reads and a value k, return the networkx de Bruijn graph object.\n",
    "    '''\n",
    "    kmers = get_kmers(reads, k)\n",
    "    graph = nx.DiGraph()\n",
    "    # Complete this function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "yoda_graph = build_graph(yoda_reads, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.draw_spring(yoda_graph, with_labels=True, node_size=1200, node_color='#eeeeff', edge_color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Challenge:** If you want a challenge, you can try to extract contigs from your graph as described in lectures. To extract a contig, we need to find an unbalanced node, then walk along the directed graph only so long as our path is unambiguous.\n",
    "\n",
    "You can find the in-degree and out-degree of nodes (number of edges leading in and out), and the edges themselves, with networkx methods like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = nx.DiGraph()\n",
    "g.add_edge('A','B')\n",
    "g.add_edge('A','C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.in_degree('A')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.out_degree('A')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OutEdgeDataView([('A', 'B'), ('A', 'C')])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.out_edges('A')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}