kmer_counting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory Python II: k-mer counting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we'll work out how to do some basic tasks required for de Bruijn graph-based genome assembly. These are:\n",
    "\n",
    "- split a read (or set of reads) into k-mers\n",
    "- find the unique k-mers, and count how many times they occur (calculate k-mer coverage)\n",
    "- find the suffix or prefix of a k-mer (this is needed to find overlaps)\n",
    "\n",
    "In order to do these things, we'll learn about string slicing, and about two new Python data types: sets and dictionaries. Like lists, sets and dictionaries are used for building collections of things, but they have some different features.\n",
    "\n",
    "*Note: This workshop has been written for Python 3. If you are using Python 2, you should be able to make everything work by first running the Python command:*\n",
    "\n",
    "    from __future__ import print_function, division"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## k-mers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The k-mers in a string are all the substrings of length `k`. So, given the string `\"ELEPHANT\"`, and k-mer length `k=4`, the k-mers are:\n",
    "    \n",
    "    ELEP\n",
    "    LEPH\n",
    "    EPHA\n",
    "    PHAN\n",
    "    HANT\n",
    "\n",
    "When we count k-mers, we want to find out which k-mers occur in our string, and how often they occur. So, for the string `\"BANANA\"` and `k=2`, we would find counts like\n",
    "\n",
    "    BA  1\n",
    "    AN  2\n",
    "    NA  2\n",
    "\n",
    "In the Assembly lectures we will see how breaking data into k-mers helps us to efficiently assemble genomes. In applying these algorithms, we will be dealing with large amounts of DNA sequence, and values of `k` will of course be much larger (e.g. in the range 20-80, and sometimes larger; values depend strongly on data). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding k-mers: string slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [first workshop](hammingdist.html) you were introduced to string indexing, using square brackets `[]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:05.425314Z",
     "start_time": "2017-07-23T11:58:05.413331Z"
    }
   },
   "outputs": [],
   "source": [
    "dna_string = \"GATTACA\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:07.529188Z",
     "start_time": "2017-07-23T11:58:07.521933Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G T A\n"
     ]
    }
   ],
   "source": [
    "# First, third and last\n",
    "print(dna_string[0], dna_string[2], dna_string[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:08.094378Z",
     "start_time": "2017-07-23T11:58:08.085837Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C\n"
     ]
    }
   ],
   "source": [
    "# Second-last\n",
    "print(dna_string[-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python also allows us to get multiple letters from a string at once using *slicing*. For instance, `s[2:4]` will give a *slice* of the string, starting at index `2` and going up to, but not including, index `4`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:09.181458Z",
     "start_time": "2017-07-23T11:58:09.175777Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TT\n"
     ]
    }
   ],
   "source": [
    "print(dna_string[2:4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also leave one (or both) of the numbers out, to go all the way from the start, or all the way until the end of the string.\n",
    "\n",
    "String slicing will be useful for finding the k-mers in a read, and for finding the overlapping suffixes and prefixes of length (k-1).\n",
    "\n",
    "Here are some slicing examples: play with this yourself until you get a feel for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:10.316242Z",
     "start_time": "2017-07-23T11:58:10.308614Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AT\n"
     ]
    }
   ],
   "source": [
    "print(dna_string[1:3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:10.914955Z",
     "start_time": "2017-07-23T11:58:10.908553Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATT\n"
     ]
    }
   ],
   "source": [
    "print(dna_string[1:4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:11.258570Z",
     "start_time": "2017-07-23T11:58:11.251537Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GATT\n"
     ]
    }
   ],
   "source": [
    "# From the start up until just before index 4\n",
    "print(dna_string[:4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:11.662554Z",
     "start_time": "2017-07-23T11:58:11.655042Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GATT\n"
     ]
    }
   ],
   "source": [
    "# From index zero (first letter) until index 4 - same result as last example\n",
    "print(dna_string[0:4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:12.037053Z",
     "start_time": "2017-07-23T11:58:12.030550Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACA\n"
     ]
    }
   ],
   "source": [
    "# From index 4 until the end\n",
    "print(dna_string[4:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:12.461198Z",
     "start_time": "2017-07-23T11:58:12.453537Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CA\n"
     ]
    }
   ],
   "source": [
    "# From the second-last character until the end\n",
    "print(dna_string[-2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:12.868409Z",
     "start_time": "2017-07-23T11:58:12.858735Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GATTACA\n",
      "GATTACA\n"
     ]
    }
   ],
   "source": [
    "# From index 0 to index 6\n",
    "# Since the length of the string is 7, this returns the entire string!\n",
    "# It's the same as dna_string[:], which means 'from the start till the end'\n",
    "print(dna_string[0:7])\n",
    "print(dna_string[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that because `dna_string[:4]` gives us the string until just *before* index `4`, and `dna_string[4:]` gives us the string starting *at* index `4`, this actually splits the string into two parts with no overlaps or missing pieces. That is:\n",
    "\n",
    "    dna_string == \"GATTACA\"\n",
    "    dna_string[:4] == \"GATT\"\n",
    "    dna_string[4:] == \"ACA\"\n",
    "    dna_string[:4] + dna_string[4:] == \"GATT\" + \"ACA\" == dna_string\n",
    "\n",
    "Similarly, `dna_string[:2]` and `dna_string[2:]`, for instance, split the string neatly into `\"GA\"` and `\"TTACA\"`.\n",
    "\n",
    "This can be a useful way to remember and understand why a slice starts on the first index, but only goes up to just *before* the second index, and does not include it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge\n",
    "\n",
    "1.  Write a function which takes as input a string and a kmer size `k`, and prints out all the kmers. It's ok to print a kmer twice. Some output should look like:\n",
    "\n",
    "    ~~~python\n",
    "    print_kmers(\"GATTACA\",4)\n",
    "    GATT\n",
    "    ATTA\n",
    "    TTAC\n",
    "    TACA\n",
    "\n",
    "    print_kmers(\"MISSISSIPPI\",5)\n",
    "    MISSI\n",
    "    ISSIS\n",
    "    SSISS\n",
    "    SISSI\n",
    "    ISSIP\n",
    "    SSIPP\n",
    "    SIPPI\n",
    "    ~~~\n",
    "\n",
    "    *Hints: You will need to use string slicing for this! You will probably want to use a `for` loop, and the `len()` function. *\n",
    "    \n",
    "2.  Modify your function so that intead of printing out each kmer, it prints out, for each k-mer, the length (k-1) prefix (all but the last letter) and the length (k-1) suffix (all but the first letter). Some sample output should look like:\n",
    "\n",
    "    ~~~python\n",
    "    print_prefixes_and_suffixes(\"GATTACA\",4)\n",
    "    GAT ATT\n",
    "    ATT TTA\n",
    "    TTA TAC\n",
    "    TAC ACA\n",
    "    \n",
    "    print_prefixes_and_suffixes(\"MISSISSIPPI\",5)\n",
    "    MISS ISSI\n",
    "    ISSI SSIS\n",
    "    SSIS SISS\n",
    "    SISS ISSI\n",
    "    ISSI SSIP\n",
    "    SSIP SIPP\n",
    "    SIPP IPPI\n",
    "    ~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Collections of things: Lists and Sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this workshop we're going to want to find all the kmers in our reads, and count kmers. To do this we'll have to be able to build collections of strings. \n",
    "\n",
    "Python has a few different data types for collecting several values together, which will be useful for us here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You've already seen lists: they are ordered collections of values. We'll review them and then introduce one new useful function: `append()`.\n",
    "\n",
    "We define lists with square brackets and commas, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:18.784859Z",
     "start_time": "2017-07-23T11:58:18.779710Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 3, 5, 7]\n"
     ]
    }
   ],
   "source": [
    "odd_numbers = [1,3,5,7]\n",
    "print(odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values in a list don't have to be numbers, but can be of any type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:19.949835Z",
     "start_time": "2017-07-23T11:58:19.944475Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['one', 'three', 'five', 'seven']\n"
     ]
    }
   ],
   "source": [
    "odd_number_strings = [\"one\",\"three\",\"five\",\"seven\"]\n",
    "print(odd_number_strings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can refer to elements of the list by indexing with square brackets, `[]`. Remember that we count from zero when doing this. We can also slice lists, like strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:20.888626Z",
     "start_time": "2017-07-23T11:58:20.882828Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "# Second element of list (index 1)\n",
    "print(odd_numbers[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:21.500559Z",
     "start_time": "2017-07-23T11:58:21.495106Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 5, 7]\n"
     ]
    }
   ],
   "source": [
    "print(odd_numbers[1:4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several built-in functions that can be applied to lists. One that we've already seen is `len()`, to find the length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:22.329286Z",
     "start_time": "2017-07-23T11:58:22.308759Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also special functions in Python which apply only to certain data types. These functions are called *methods* and are used by writing them after the variable name, with a dot, like `variable.function()`. For instance, lists have a method called `append()` which adds an item to the end of the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:23.368606Z",
     "start_time": "2017-07-23T11:58:23.358369Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before append: [1, 3, 5, 7]\n",
      "After append: [1, 3, 5, 7, 9]\n"
     ]
    }
   ],
   "source": [
    "odd_numbers = [1,3,5,7]\n",
    "print(\"Before append:\",odd_numbers)\n",
    "odd_numbers.append(9)\n",
    "print(\"After append:\",odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `append()` function doesn't return a new list. It modifes the list it is attached to. \n",
    "\n",
    "`append()` is useful when building up a list piece by piece, for example, in a `for` loop.\n",
    "\n",
    "`append()` is not the only method that works on lists. For instance, try `reverse()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:24.573536Z",
     "start_time": "2017-07-23T11:58:24.565441Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7, 5, 3, 1]\n"
     ]
    }
   ],
   "source": [
    "odd_numbers = [1,3,5,7]\n",
    "odd_numbers.reverse()\n",
    "print(odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's no need to learn all these functions at once. You will learn about them as you need them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge\n",
    "\n",
    "Write a function which takes a string and a kmer length k as input. This time, instead of printing out the kmers directly inside the function, build and return a *list* of all the kmers in the string. As before, it's ok to include a kmer more than once. Some example output should look like:\n",
    "\n",
    "~~~python\n",
    "print(find_kmers(\"GATTACA\",4))\n",
    "['GATT', 'ATTA', 'TTAC', 'TACA']\n",
    "\n",
    "print(find_kmers(\"MISSISSIPPI\",5))\n",
    "['MISSI', 'ISSIS', 'SSISS', 'SISSI', 'ISSIP', 'SSIPP', 'SIPPI']\n",
    "\n",
    "print(find_kmers(\"MISSISSIPPI\",3))\n",
    "['MIS', 'ISS', 'SSI', 'SIS', 'ISS', 'SSI', 'SIP', 'IPP', 'PPI']\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sets are a Python data type which, like lists, hold collections of values. The main differences between sets and lists are:\n",
    "\n",
    "- a set can only have **one** copy of each value at most: values must be unique\n",
    "- the values in a set are not in order - Python may return them in a random order each time!\n",
    "\n",
    "If you're familiar with sets as used in mathematics, Python's sets are actually quite similar, and can be used for set calculations by calling set methods like `.union()` and `.intersection()`. However we won't use those features in this workshop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can define a set by passing in a list to `set()`, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:27.967693Z",
     "start_time": "2017-07-23T11:58:27.961458Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1, 3, 5, 7}\n"
     ]
    }
   ],
   "source": [
    "odd_number_set = set([1,3,5,7])\n",
    "print(odd_number_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, this looks a lot like a list! But what if there are repeated numbers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:28.842130Z",
     "start_time": "2017-07-23T11:58:28.833498Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1, 3, 5, 7}\n"
     ]
    }
   ],
   "source": [
    "odd_number_set = set([1,3,3,5,7,1])\n",
    "print(odd_number_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The set will contain only unique values - i.e. only one of each value that's been added to it. This can be useful when you're trying to detect all the unique items in a large collection of things:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:29.786462Z",
     "start_time": "2017-07-23T11:58:29.765439Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'7', '3', '6', '1', '5'}\n"
     ]
    }
   ],
   "source": [
    "big_list = ['3', '7', '3', '3', '5', '3', '5', '1', '1', '7', '3', '5', '5',\n",
    "       '3', '5', '5', '5', '1', '7', '7', '1', '1', '7', '7', '3', '7',\n",
    "       '3', '5', '1', '5', '5', '3', '3', '5', '1', '1', '1', '1', '7',\n",
    "       '3', '7', '1', '5', '1', '7', '3', '5', '5', '6', '5', '7', '1',\n",
    "       '5', '1', '5', '1', '3', '7', '7', '5', '3', '1', '1', '5', '5',\n",
    "       '5', '1', '5', '1', '7', '3', '7', '1', '7', '5', '7', '1', '5',\n",
    "       '3', '5', '3', '5', '3', '5', '5', '1', '3', '5', '3', '7', '1',\n",
    "       '5', '7', '3', '1', '7', '5', '7', '7', '1']\n",
    "unique_values = set(big_list)\n",
    "print(unique_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One function which might be useful for building a collection of k-mers is the `add()` method for sets, called like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:30.795502Z",
     "start_time": "2017-07-23T11:58:30.780910Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before add: {1, 3, 5, 7}\n",
      "After add: {1, 3, 5, 7, 11}\n"
     ]
    }
   ],
   "source": [
    "odd_numbers = set([1,3,5,7])\n",
    "print(\"Before add:\",odd_numbers)\n",
    "odd_numbers.add(11)\n",
    "print(\"After add:\",odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the values are not in any particular order! Sometimes you may notice Python apparently printing out sets in a certain order, but you should never rely on this - it could be random another time.\n",
    "\n",
    "Adding an item to a set that is already there has no effect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:31.835868Z",
     "start_time": "2017-07-23T11:58:31.827507Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before add: {1, 3, 5, 7}\n",
      "After add: {1, 3, 5, 7}\n"
     ]
    }
   ],
   "source": [
    "odd_numbers = set([1,3,5,7])\n",
    "print(\"Before add:\",odd_numbers)\n",
    "odd_numbers.add(3)\n",
    "print(\"After add:\",odd_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's also possible to use a `for` loop with sets, just like lists. The difference is that we can't guarantee what order we will go through the items in, since sets have no order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:32.717486Z",
     "start_time": "2017-07-23T11:58:32.709582Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "3\n",
      "9\n",
      "7\n"
     ]
    }
   ],
   "source": [
    "odd_number_set = set([1,3,1,7,9,3])\n",
    "for num in odd_number_set:\n",
    "    print(num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge\n",
    "\n",
    "Write a function which takes a string and a kmer length k, and returns a set of all the *unique* kmers in that string. Some example output should look like this:\n",
    "\n",
    "~~~python\n",
    "print(find_unique_kmers(\"GATTACA\",4))\n",
    "set(['GATT', 'TACA', 'TTAC', 'ATTA'])\n",
    "\n",
    "print(find_unique_kmers(\"MISSISSIPPI\",5))\n",
    "set(['SISSI', 'MISSI', 'SSIPP', 'ISSIP', 'ISSIS', 'SSISS', 'SIPPI'])\n",
    "\n",
    "print(find_unique_kmers(\"MISSISSIPPI\",3))\n",
    "set(['SIS', 'SIP', 'ISS', 'PPI', 'IPP', 'SSI', 'MIS'])\n",
    "~~~\n",
    "\n",
    "Since sets don't preserve the order of the items inside them, your sets may display the kmers in a different order - this is still correct.\n",
    "\n",
    "If you compare your output to the output from the previous challenge, you should see that if the input string has repeats, there are less *unique* kmers than total observed kmers:\n",
    "\n",
    "~~~python\n",
    "print(find_kmers(\"MISSISSIPPI\",3))\n",
    "['MIS', 'ISS', 'SSI', 'SIS', 'ISS', 'SSI', 'SIP', 'IPP', 'PPI']\n",
    "\n",
    "print(len(find_kmers(\"MISSISSIPPI\",3)))\n",
    "9\n",
    "\n",
    "print(find_unique_kmers(\"MISSISSIPPI\",3))\n",
    "set(['SIS', 'SIP', 'ISS', 'PPI', 'IPP', 'SSI', 'MIS'])\n",
    "\n",
    "print(len(find_unique_kmers(\"MISSISSIPPI\",3)))\n",
    "7\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Counting k-mers: Dictionaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to look at one more Python data type for storing collections of values: the **dictionary**.\n",
    "\n",
    "This is another built-in Python data type. Dictionaries are a bit like sets - they store unordered collections of items. But dictionaries store key-value pairs. \"Keys\" are used for indexing and must be unique. \"Values\" do not have to be unique.\n",
    "\n",
    "So, just like in sets, all keys in a dictionary must be unique - only one of each can be stored. Also like a set, the keys in a dictionary are stored in no particular order. \n",
    "\n",
    "Here's an example dictionary, which we've defined using curly braces `{}`, and by putting a colon `:` between the each key and its corresponding value. This dictionary stores some people's heights in centimetres."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:35.519963Z",
     "start_time": "2017-07-23T11:58:35.507135Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Sam': 201, 'Fiona': 167, 'Quentin': 167}\n"
     ]
    }
   ],
   "source": [
    "heights = {\"Sam\":201, \"Fiona\":167, \"Quentin\":167}\n",
    "print(heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:36.105973Z",
     "start_time": "2017-07-23T11:58:36.097712Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'dict'>\n"
     ]
    }
   ],
   "source": [
    "print(type(heights))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can retrieve a value using its key. We use square brackets `[]`, just like getting an item from a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:38.301612Z",
     "start_time": "2017-07-23T11:58:38.296886Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "201\n"
     ]
    }
   ],
   "source": [
    "print(heights[\"Sam\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying to look up a key that isn't in the dictionary will give an error (a `KeyError`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:39.600013Z",
     "start_time": "2017-07-23T11:58:39.278287Z"
    }
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'Anne'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-28-96f1ce2d9468>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mheights\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Anne\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m: 'Anne'"
     ]
    }
   ],
   "source": [
    "print(heights[\"Anne\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also add a key-value pair to the dictionary this way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:42.364266Z",
     "start_time": "2017-07-23T11:58:42.359231Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Sam': 201, 'Fiona': 167, 'Quentin': 167, 'Mary': 180}\n"
     ]
    }
   ],
   "source": [
    "heights[\"Mary\"] = 180\n",
    "print(heights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can assign a new value to an existing key. If we do this, the old value is lost. This means keys will always be unique - trying to add a new key that is the same simply overwrites the old one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:44.035592Z",
     "start_time": "2017-07-23T11:58:44.029443Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Sam': 201, 'Fiona': 167, 'Quentin': 170, 'Mary': 180}\n"
     ]
    }
   ],
   "source": [
    "# apparently Quentin is still growing.. update his height to the new value\n",
    "heights[\"Quentin\"] = 170\n",
    "print(heights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can avoid errors like the one we got with `height[\"Anne\"]` by checking to see if a key is in the dictionary or not. We use the Python keyword `in`, which returns `True` or `False`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:45.019295Z",
     "start_time": "2017-07-23T11:58:45.013947Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n"
     ]
    }
   ],
   "source": [
    "print(\"Anne\" in heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:53.919004Z",
     "start_time": "2017-07-23T11:58:53.913088Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(\"Fiona\" in heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:58:54.771754Z",
     "start_time": "2017-07-23T11:58:54.764328Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fiona 167\n"
     ]
    }
   ],
   "source": [
    "if \"Fiona\" in heights:\n",
    "    print(\"Fiona\",heights[\"Fiona\"])\n",
    "if \"Anne\" in heights:\n",
    "    print(\"Anne\",heights[\"Anne\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dictionaries have methods called `keys()` and `values()` which return lists of keys and values (watch out - these can be in random order!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:00.050188Z",
     "start_time": "2017-07-23T11:59:00.044514Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['Sam', 'Fiona', 'Quentin', 'Mary'])\n"
     ]
    }
   ],
   "source": [
    "print(heights.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:00.832022Z",
     "start_time": "2017-07-23T11:59:00.827302Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_values([201, 167, 170, 180])\n"
     ]
    }
   ],
   "source": [
    "print(heights.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, our keys were strings and our values were numbers, but that doesn't have to be the case. Our keys could be numbers, for instance. Our values could be strings, or could even be lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:02.386722Z",
     "start_time": "2017-07-23T11:59:02.380839Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{12: [3, 4], 100: [2, 2, 5, 5], 9: [3, 3]}\n"
     ]
    }
   ],
   "source": [
    "prime_factors = {12: [3,4],\n",
    "                 100: [2,2,5,5],\n",
    "                 9: [3,3]}\n",
    "print(prime_factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:03.369656Z",
     "start_time": "2017-07-23T11:59:03.364841Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 4]\n"
     ]
    }
   ],
   "source": [
    "print(prime_factors[12])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Don't be fooled - even though writing `prime_factors[12]` with a number between the square brackets looks a bit like list indexing, it's not. The difference is that our dictionary only contains the keys `9`, `12`, and `100`, because those are the only ones we've added. A list that contains the index `9` must also contain `0,1,2,3,4,5,6,7` and `8`, because the indices just give the order of the items in the list."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One common use of dictionaries, which is what we'd like to do, is to count things. We can make the keys the things we want to count, and the values the number of times we've seen each one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:28.970646Z",
     "start_time": "2017-07-23T11:59:28.947278Z"
    }
   },
   "outputs": [],
   "source": [
    "big_list = ['3', '7', '3', '3', '5', '3', '5', '1', '1', '7', '3', '5', '5',\n",
    "       '3', '5', '5', '5', '1', '7', '7', '1', '1', '7', '7', '3', '7',\n",
    "       '3', '5', '1', '5', '5', '3', '3', '5', '1', '1', '1', '1', '7',\n",
    "       '3', '7', '1', '5', '1', '7', '3', '5', '5', '6', '5', '7', '1',\n",
    "       '5', '1', '5', '1', '3', '7', '7', '5', '3', '1', '1', '5', '5',\n",
    "       '5', '1', '5', '1', '7', '3', '7', '1', '7', '5', '7', '1', '5',\n",
    "       '3', '5', '3', '5', '3', '5', '5', '1', '3', '5', '3', '7', '1',\n",
    "       '5', '7', '3', '1', '7', '5', '7', '7', '1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many times does each number appear in `big_list`? We could write a `for` loop and store the counts in variables. But we'd need a new variable for the number of `'1'`s, the number of `'3'`s, the number of `'5'`s, etc, and we might not know in advance how many variables we need.\n",
    "\n",
    "This is where a dictionary is handy: each time we see a new value that needs counting, we can add it to the dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-07-23T11:59:37.119902Z",
     "start_time": "2017-07-23T11:59:37.100823Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'3': 21, '7': 22, '5': 31, '1': 25, '6': 1}\n"
     ]
    }
   ],
   "source": [
    "# Start by making counts an empty dictionary, using {}\n",
    "counts = {}\n",
    "\n",
    "# Now loop through the values in big_list and count them\n",
    "for num in big_list:\n",
    "    \n",
    "    # Check to see if this key is already in the dictionary\n",
    "    # If not, add it with an initial count of zero\n",
    "    if not (num in counts):\n",
    "        counts[num] = 0\n",
    "    \n",
    "    # Now that we are sure the key is in the dictionary, we can increment the count\n",
    "    counts[num] += 1\n",
    "\n",
    "# After the loop is finished, counts should contain the right counts for each number seen\n",
    "print(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we could get the unique items in `big_list` using `counts.keys()`. This will give us the same unique items that we found before using `set()`. The keys of a dictionary are a lot like sets, and sets behave a lot like dictionaries with no values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Side note: a more advanced solution is to use the [defaultdict](https://docs.python.org/3/library/collections.html#collections.defaultdict) type. This doesn't require us to initialise each count to zero. Instead we make `counts` a defaultdict of type int, and any key not present will automatically be added when needed, with starting value zero:\n",
    "\n",
    "```\n",
    "from collections import defaultdict\n",
    "\n",
    "counts = defaultdict(int)\n",
    "for num in big_list:\n",
    "    counts[num] += 1\n",
    "\n",
    "print(counts)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dictionaries: algorithmic considerations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, dictionaries (and sets) are an implementation of **hash tables**, and work by carrying out a hash function on the keys to decide where to store them in the computer's memory. This means that dictionaries are very fast at retrieving information. Even if you have a very large dictionary, finding an item in it using a key takes a roughly constant amount of time that does not grow with the size of the dictionary - the performance of the algorithm is O(1). \n",
    "\n",
    "\"Exact string matching\" is extremely fast and simple since we just need to test for equality, not carry out an alignment. If we store our known k-mers as keys of a dictionary, and we are given a new k-mer which we would like to match against them, then we also can search for matches across a huge dataset of known k-mers in O(1) time simply by looking up the key.  \n",
    "\n",
    "Performing *inexact* string matching can't be done using a hash table, since the hash function only gives the same result if the key is exactly the same. For approximate string matching, we need some kind of alignment algorithm. Not only is this slower for each string, we may need to carry it out on a large number of strings to find a match. This is why, when there are many reads, a hash-based algorithm like De Bruijn assembly is faster than an alignment-based assembly algorithm. It's also why some alignment algorithms use *seeds*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge\n",
    "\n",
    "Write a function which takes a string and a kmer length k, and returns a dictionary. The keys of the dictionary should be the kmers found in the string, and the values should be the kmer counts.\n",
    "\n",
    "~~~python\n",
    "print(count_kmers(\"GATTACA\",4))\n",
    "{'ATTA': 1, 'GATT': 1, 'TACA': 1, 'TTAC': 1}\n",
    "\n",
    "print(count_kmers(\"MISSISSIPI\",3))\n",
    "{'IPI': 1, 'ISS': 2, 'MIS': 1, 'SIP': 1, 'SIS': 1, 'SSI': 2}\n",
    "\n",
    "print(count_kmers(\"CTCGGGGGCTATTAATACCTAAGTGCTCGGGGGCTATTAATACCTAAGTGGGGCTATT\",3))\n",
    "{'AAG': 2,\n",
    " 'AAT': 2,\n",
    " 'ACC': 2,\n",
    " 'AGT': 2,\n",
    " 'ATA': 2,\n",
    " 'ATT': 3,\n",
    " 'CCT': 2,\n",
    " 'CGG': 2,\n",
    " 'CTA': 5,\n",
    " 'CTC': 2,\n",
    " 'GCT': 4,\n",
    " 'GGC': 3,\n",
    " 'GGG': 8,\n",
    " 'GTG': 2,\n",
    " 'TAA': 4,\n",
    " 'TAC': 2,\n",
    " 'TAT': 3,\n",
    " 'TCG': 2,\n",
    " 'TGC': 1,\n",
    " 'TGG': 1,\n",
    " 'TTA': 2}\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "228px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}