Day 03 - DNA Fragmentation/Starter/notebook.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "# Day 3 notebook\n",
                "\n",
                "The objectives of this notebook are to practice\n",
                "* how data are generated from universal DNA arrays\n",
                "* how data are generated from a DNA shotgun sequencing experiment\n",
                "* plotting with matplotlib\n",
                "* raising and handling exceptions"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## The Ebola virus genome\n",
                "\n",
                "In 2014, an outbreak of the Ebola virus devastated human populations in West Africa, killing more than [11,000 people](https:\/\/www.cdc.gov\/vhf\/ebola\/history\/2014-2016-outbreak\/index.html).  During the outbreak, a group of researchers obtained virus samples from 78 patients, sequenced and assembled the Ebola virus genomes from these samples, and used the genome sequences to identify the origin of the outbreak ([Gire et al. 2014](http:\/\/science.sciencemag.org\/content\/345\/6202\/1369.long)).  In this notebook, we will work with one of the genome assemblies produced by this study.\n",
                "\n",
                "The Ebola virus is a negative-sense single-stranded RNA virus.  This means that its genome is actually an RNA, rather than a DNA, molecule.  It has a genome that is in the \"negative-sense\" orientation, meaning that the mRNAs encoded by the genes in its genome are in the reverse complement orientation.  For this class of viruses, transcription is actually an RNA to RNA process!\n",
                "\n",
                "Although the virus's genome is an RNA molecule, we will work with it as though it were DNA, as one typically needs to convert from RNA to DNA (via reverse transcription) before using the latest sequencing technologies."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Problem 1: Reading from a FASTA file (1 POINT)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "The genome assembly that we will work with in this notebook is saved in the file `ebola_virus_isolate_G3686.fasta` and was downloaded from its [Genbank record](https:\/\/www.ncbi.nlm.nih.gov\/nuccore\/KM034562).  This file is in [FASTA format](https:\/\/en.wikipedia.org\/wiki\/FASTA_format), one of the most common formats for storing biological sequence data.  To save you some time, I have provided you with a `fasta` module that you can use to read in the viral sequence that is in this file.  In the code block below, use one of the functions from the `fasta` module to read in the name and sequence of the ebola virus genome assembly."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "import fasta\n",
                "help(fasta.read_sequences_from_fasta_file)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Here is an example of the output of the `read_sequences_from_fasta_file` on a small file:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "fasta.read_sequences_from_fasta_file(\"small_fasta_example.fasta\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "In the code block below, use the `read_sequences_from_fasta_file` on the Ebola genome FASTA file and assign the name of the viral sequence to the variable `ebola_sequence_name`\n",
                "and assign the sequence to the variable `ebola_sequence`."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "ebola_sequence_filename = \"ebola_virus_isolate_G3686.fasta\"\n",
                "\n",
                "###\n",
                "### YOUR CODE HERE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "reading_ebola_fasta",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for reading in the ebola sequence from a FASTA file\n",
                "assert isinstance(ebola_sequence_name, str) and len(ebola_sequence_name) == 101, \"Failed to assign ebola_sequence_name\"\n",
                "assert isinstance(ebola_sequence, str) and len(ebola_sequence) == 18957, \"Failed to assign ebola_sequence\"\n",
                "print(\"SUCCESS: reading in the ebola sequence from a FASTA file passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Simulating universal DNA array data\n",
                "\n",
                "We will start by generating the sort of data that one would hope to get from a universal DNA array.  Recall that a universal DNA array gives you information about the presence or absence of each possible DNA \"$k$-mer\", which is a DNA substring of length $k$, within your sample of DNA.  Here, we will assume that we have universal DNA arrays that give us perfect error-free data.  In practice, such arrays are far from perfect, and are no longer used for sequencing purposes.  However, the type of data they generate is useful to think about from a computational perspective."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Problem 2: $k$-mer spectrum of a DNA sequence (1 POINT)\n",
                "\n",
                "Write a function `kmer_spectrum` below that takes as input a DNA sequence and a length $k$ and returns a `set` object containing all *unique* $k$-mers within the input sequence, which we will refer to as its \"$k$-mer spectrum\".  This is an idealized form of the data that a universal DNA array would give us.  For simplicitly, in this activity we will treat a DNA sequence as a single-stranded molecule.  In practice, one would have to handle the fact that DNA is double-stranded.  You may want to get familiar with how to work with Python [`set`](https:\/\/docs.python.org\/3\/library\/stdtypes.html#set) objects."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "locked": false,
                    "solution": false
                }
            },
            "outputs": [],
            "source": [
                "def kmer_spectrum(dna_sequence, k):\n",
                "    \"\"\"Gives the k-mer spectrum of a DNA sequence.\n",
                "    \n",
                "    Args:\n",
                "        dna_sequence: a DNA sequence\n",
                "        k: an integer length\n",
                "    Returns:\n",
                "        A set object containing all unique k-mers within dna_sequence\n",
                "    \"\"\"\n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "kmer_spectrum",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for kmer_spectrum\n",
                "assert kmer_spectrum(\"ACGTCGA\", 3) == {'ACG', 'CGA', 'CGT', 'GTC', 'TCG'}\n",
                "assert kmer_spectrum(\"TCTATCT\", 3) == {'ATC', 'CTA', 'TAT', 'TCT'}\n",
                "assert kmer_spectrum(\"ACGTCGA\", 2) == {'AC', 'CG', 'GA', 'GT', 'TC'}\n",
                "assert kmer_spectrum(\"AAAAAAA\", 2) == {'AA'}\n",
                "assert kmer_spectrum(\"T\", 3) == set()\n",
                "print(\"SUCCESS: kmer_spectrum passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Problem 3: $k$-mer spectrum of the Ebola genome sequence (1 POINT)\n",
                "\n",
                "Use your `kmer_spectrum` function to compute the spectrum for the ebola genome sequence for $k = 1, 2, \\ldots, 25$.  Assign to the variable `ebola_spectrum_sizes` a *list* of the *sizes* (number of elements) of these spectra."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "k_values = range(1, 26)\n",
                "###\n",
                "### YOUR CODE HERE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "ebola_spectrum_sizes",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for ebola_spectrum_sizes\n",
                "assert len(ebola_spectrum_sizes) == 25, \"Failed required size\"\n",
                "assert ebola_spectrum_sizes[0] == 4, \"Failed on value of first entry\"\n",
                "###\n",
                "### AUTOGRADER TEST - DO NOT REMOVE\n",
                "###\n",
                "print(\"SUCCESS: ebola_spectrum_sizes passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Next let us make a plot of the size of the ebola spectra vs. $k$.  To do this, we will use the Python package [`matplotlib`](https:\/\/matplotlib.org\/).  For simplicity, we will use the Pyplot interface of this package.  I encourage you to read through the [Pyplot tutorial](https:\/\/matplotlib.org\/tutorials\/introductory\/pyplot.html)."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "import matplotlib.pyplot as plt\n",
                "plt.figure(figsize = (8, 5))\n",
                "plt.plot(k_values, ebola_spectrum_sizes, 'o')\n",
                "plt.xlabel('k')\n",
                "plt.ylabel('# of distinct k-mers in ebola genome')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Examine the plot above.  \n",
                "* Why does it have this shape?\n",
                "* For which value of $k$ are all $k$-mers in the Ebola genome unique?  \n",
                "* Is this value of $k$ surprising?\n",
                "Talk with your neighbors and write your thoughts about these questions below.\n",
                "\n",
                "*WRITE YOUR THOUGHTS HERE*"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Simulating reads from a DNA shotgun sequencing experiment\n",
                "\n",
                "Recall that a DNA shotgun sequencing experiment involves three steps to generate read data:\n",
                "1. Fragmentation of many identical copies of the DNA of interest.\n",
                "2. Random sampling of a subset of the resulting DNA fragments.\n",
                "3. Sequencing a read of limited length from one end from each of the sampled fragments.\n",
                "\n",
                "We will write functions to model each aspect of this proecss and then put them all together to simulate a set of reads from a DNA shotgun sequencing experiment."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### PROBLEM 4 - Fragmentation simulation (1 POINT)\n",
                "\n",
                "We will first model the fragmentation of the DNA sequence, for which we have a very large number of identical copies.  Each copy of the DNA will be fragmented randomly, likely resulting in a unique fragmentation pattern for each DNA molecule.  Fragmentation is a complex process to model, however, for the purposes of this activity, we will simulate it in a very simple manner, by choosing `k` distinct positions along a given DNA molecule at which to cut the molecule.  Each cut position will be chosen randomly and uniformly across the entire DNA molecule.  We will specify the cut position by the index (0-based) of the base to the *right* of the cut.\n",
                "\n",
                "Write a function that takes as input the length of the DNA sequence that is to be fragmented and the number of positions, `k`, at which to cut the DNA, and returns a *sorted* list of the intervals of the resulting DNA fragments.  Each interval should be specified by a tuple, of the form `(start, end)`, which will represent the fragment of DNA beginning at `start` (0-based coordinates) and spanning up to, **but not including**, `end` (another way to think about it is that these are the indices that you would use to slice the Python string to obtain the fragment's sequence).\n",
                "\n",
                "To generate random cut points, you will use the pseuodorandom number generating functionality provided by Python's [random module](https:\/\/docs.python.org\/3\/library\/random.html).  For the purposes of being to able to test your code, you will need to use **a single call** of the [`random.sample`](https:\/\/docs.python.org\/3\/library\/random.html#random.sample) function to randomly generate all of the cut points simultaneously within your function."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "import random\n",
                "def fragment_dna(length, k):\n",
                "    \"\"\"Randomly cuts a sequence of a given length into k fragments.\n",
                "    \n",
                "    Args:\n",
                "        length: an integer specifying the length of the sequence\n",
                "        k: the number of fragments that should result from random cuts\n",
                "    Returns:\n",
                "        A sorted list of fragment intervals, where each interval\n",
                "        is specified as a tuple (start, end).  For each interval, \n",
                "        the coordinates are 0-based and the end coordinate  \n",
                "        *not* included in the fragment.\n",
                "    \"\"\"\n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "fragment_dna",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for fragment_dna\n",
                "# We will use the seed function to reset the pseudorandom number generator to a specific state\n",
                "# so that we can test inputs with known outputs\n",
                "def fragment_dna_with_seed(length, k, seed):\n",
                "    random.seed(seed)\n",
                "    return fragment_dna(length, k)\n",
                "assert fragment_dna_with_seed(100, 3, 42) == [(0, 4), (4, 15), (15, 82), (82, 100)], \"Failed on 100, 3 (seed 42)\"\n",
                "assert fragment_dna_with_seed(100, 2, 43) == [(0, 5), (5, 37), (37, 100)], \"Failed on 100, 3 (seed 43)\"\n",
                "assert fragment_dna_with_seed(100, 1, 44) == [(0, 53), (53, 100)], \"Failed on 100, 3 (seed 44)\"\n",
                "assert fragment_dna_with_seed(100, 0, 45) == [(0, 100)], \"Failed on 100, 3 (seed 45)\"\n",
                "assert fragment_dna_with_seed(5, 4, 46) == [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5)], \"Failed on 5, 4 (seed 46)\"\n",
                "print(\"SUCCESS: fragment_dna passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "For your convenience, below is a function that will fragment multiple copies of a genome of a given length and return a list with all fragment intervals combined and sorted."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def fragment_multiple_dna_copies(length, k, num_copies):\n",
                "    \"\"\"Randomly cuts multiple sequences of a given length, each sequence\n",
                "    cut into k fragments.\n",
                "    \n",
                "    Args:\n",
                "        length: an integer specifying the length of each sequence\n",
                "        k: the number of fragments that should result from each sequence\n",
                "        num_copies: the number of sequences\n",
                "    Returns:\n",
                "        A sorted list of fragment intervals from all sequences, where intervals\n",
                "        are represented in the same way as for the fragment_dna function.\n",
                "    \"\"\"\n",
                "     \n",
                "    intervals = []\n",
                "    for i in range(num_copies):\n",
                "        intervals.extend(fragment_dna(length, k))\n",
                "    return sorted(intervals)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### PROBLEM 5 - Filtering fragments (1 POINT)\n",
                "\n",
                "We will next model the random sampling of a set of fragments to be sequenced by our sequencing technology of chioce.  Typically, only fragments longer than a certain minimum length will be sequenceable, so we will first filter the set of possible fragments by length before selecting a subset.\n",
                "\n",
                "Let us first write a function that filters fragments by length.  This function should take a list of fragment intervals (as would be produced by `fragment_multiple_dna_copies`) and a minimum length and returns a *new* list that contains only those intervals whose lengths are at least the minimum, in the same order as they appeared in the input list.\n",
                "\n",
                "You might consider using a list comprehension with an `if` clause for this purpose.  Alternatively, you might consider using the [`filter`](https:\/\/docs.python.org\/3\/library\/functions.html#filter) higher-order function with a [`lambda`](https:\/\/docs.python.org\/3\/tutorial\/controlflow.html#lambda-expressions) expression."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def filter_for_long_fragments(fragment_list, min_length):\n",
                "    \"\"\"Returns a list of all fragments that are at least as long as a\n",
                "    given minimum length.\n",
                "\n",
                "    Args:\n",
                "        fragment_list: a list of fragments, in the format returned by fragment_dna\n",
                "        min_length: the minimum length of a fragment to return in the output\n",
                "    Returns:\n",
                "        A list of fragments.\n",
                "    \"\"\"    \n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "filter_for_long_fragments",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for filter_for_long_fragments\n",
                "test_fragment_list = [(0, 4), (4, 15), (15, 82), (82, 100)]\n",
                "assert filter_for_long_fragments(test_fragment_list, 20) == [(15, 82)], \"Failed for min_length = 20\"\n",
                "assert filter_for_long_fragments(test_fragment_list, 5) == [(4, 15), (15, 82), (82, 100)], \"Failed for min_length = 5\"\n",
                "assert filter_for_long_fragments(test_fragment_list, 0) == test_fragment_list, \"Failed for min_length = 0\"\n",
                "assert filter_for_long_fragments(test_fragment_list, 100) == [], \"Failed for min_length = 100\"\n",
                "print(\"SUCCESS: filter_for_long_fragments passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### PROBLEM 6 - Sampling fragments (1 POINT)\n",
                "Next, we will write a function to sample fragments.  We will again use the `random.sample` function from the `random` module.  Note that this function raises an exception if you ask it to sample more elements than are present in the sequence input (since it samples *without* replacement).  For example,"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "try:\n",
                "    random.sample([1, 2, 3, 4], 6)\n",
                "except Exception as err:\n",
                "    print(err)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "This error is not as user-friendly as we might like it to be for someone who is using our functions to randomly sample DNA fragments from a large pool of DNA fragments.  Let us remedy this by explicitly checking that the number of fragments requested is smaller than the size of the pool and if it is not, raising a `ValueError` exception with a more informative message.  Let the format of the error message be `\"The number of fragments to be sampled (INSERT_NUM_FRAGMENTS_HERE) is larger than the fragment pool (INSERT_FRAGMENT_LIST_LENGTH_HERE)`\"."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def sample_fragments(fragment_list, num_fragments):\n",
                "    \"\"\"Randomly samples a specified number of fragments from a \n",
                "    a list of fragments without replacement.\n",
                "    \n",
                "    Args:\n",
                "        fragment_list: a list of fragments, in the format returned by fragment_dna\n",
                "        num_fragments: the number of fragments to randomly sample from the list\n",
                "    Returns:\n",
                "        A sorted list of fragments.\n",
                "    \"\"\"\n",
                "    \n",
                "    # Check for num_fragments being valid\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "    return sorted(random.sample(fragment_list, num_fragments))"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "sample_fragments",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for sample_fragments\n",
                "test_fragment_list = [(0, 4), (4, 15), (15, 82), (82, 100)]\n",
                "random.seed(42)\n",
                "assert sample_fragments(test_fragment_list, 2) == [(0, 4), (82, 100)]\n",
                "random.seed(42)\n",
                "assert sample_fragments(test_fragment_list, 4) == [(0, 4), (4, 15), (15, 82), (82, 100)]\n",
                "\n",
                "# The tests below check that you are raising an exception properly\n",
                "random.seed(42)\n",
                "try:\n",
                "    sample_fragments(test_fragment_list, 5)\n",
                "except ValueError as err:\n",
                "    assert str(err) == \"The number of fragments to be sampled (5) is larger than the fragment pool (4)\", \"Failed to raise the correct error message\"\n",
                "except Exception as err:\n",
                "    raise AssertionError(\"Raised something other than a ValueError\")\n",
                "else:\n",
                "    raise AssertionError(\"Failed to raise an exception with num_fragments = 5\")\n",
                "print(\"SUCCESS: sample_fragments passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### PROBLEM 7 - Sequencing fragments (1 POINT)\n",
                "Finally, a sequencer will be applied to the sampled fragments to read one or both of its ends. We will assume that we are using a sequencing technology (such as Illumina) that reads a fixed length from one strand of the fragment starting at the 5' end of that strand and reading in the 5' to 3' direction.  Thus, the end of the fragment that is sequenced depends on which strand is chosen.  We will assume that for each fragment independently, the strand that is sequenced is chosen at random."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Write a function `read_fragment` that takes as input a fragment interval and a maximum read length and returns an \"oriented\" read interval `(start, end, strand)` representing the interval of the genome that is read by the sequencer and the orientation (strand) of that read.  For example, with a read length of 100 the fragment interval (500, 800) would result in an oriented read interval of either `(500, 600, +)` or `(700, 800, -)`.  We will use `+` to represent the forward strand (the strand whose sequence is represented by the genome string) and `-` to represent the reverse strand."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def read_fragment(fragment, read_length):\n",
                "    \"\"\"Returns the interval sequenced from one end of a given fragment, \n",
                "    given a maximum read length and a randomly sampled orientation.\n",
                "    \n",
                "    Args:\n",
                "        fragment: a fragment, represented as a tuple (start, end)\n",
                "        read_length: the maximum length of a read\n",
                "    Returns:\n",
                "        The interval of the read, represented as a tuple (start, end, orientation),\n",
                "        where orientation is \"+\" for a forward strand read and \"-\" for a reverse strand\n",
                "        read. If the fragment is shorter than the maximum read length, the entire fragment\n",
                "        is read.\n",
                "    \"\"\"\n",
                "    \n",
                "    strand = random.choice(\"+-\")\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "read_fragment",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for read_fragment\n",
                "random.seed(42)\n",
                "assert read_fragment((500, 800), 100) == (500, 600, '+'), \"Failed on (500, 800), 100 (seed 42)\"\n",
                "random.seed(44)\n",
                "assert read_fragment((500, 800), 100) == (700, 800, '-'), \"Failed on (500, 800), 100 (seed 44)\"\n",
                "random.seed(42)\n",
                "assert read_fragment((500, 800), 400) == (500, 800, '+'), \"Failed on (500, 800), 400 (seed 42)\"\n",
                "print(\"SUCCESS: read_fragment passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Problem 8 - Putting it all together: simulating shotgun sequencing oriented read intervals (1 POINT)\n",
                "\n",
                "We will now put all of the functions you have written for simulating fragmentation, filtering fragments by length, random selection, and fragment sequencing together.  Write a function below to do this using the documentation string\n",
                "as the specificiation."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def simulate_shotgun_sequencing_read_intervals(genome_length, \n",
                "                                               num_copies, \n",
                "                                               k,\n",
                "                                               min_fragment_length,\n",
                "                                               num_fragments, \n",
                "                                               read_length):\n",
                "    \"\"\"Simulates shotgun sequence read intervals.\n",
                "    \n",
                "    Args:\n",
                "        genome_length: the integer length of the DNA sequence\n",
                "        num_copies: the number of identical copies of the DNA in the sample\n",
                "        k: the number of cuts in each DNA during fragmentation\n",
                "        min_fragment_length: the minimum length of a fragment to be sequenced\n",
                "        num_fragments: the number of fragments to sequence\n",
                "        read_length: the length of the reads\n",
                "    Returns:\n",
                "        A list of oriented read intervals, each represented as a tuple (start, end, orientation)\n",
                "    \"\"\"    \n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "simulate_shotgun_sequencing_read_intervals",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for simulate_shotgun_sequencing_read_intervals\n",
                "random.seed(42)\n",
                "assert simulate_shotgun_sequencing_read_intervals(genome_length = 100,\n",
                "                                                  num_copies = 10,\n",
                "                                                  k = 5,\n",
                "                                                  min_fragment_length = 10,\n",
                "                                                  num_fragments = 3,\n",
                "                                                  read_length = 10) == [(34, 44, '-'), (36, 46, '+'), (77, 87, '-')]\n",
                "print(\"SUCCESS: simulate_shotgun_sequencing_read_intervals passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Simulating read intervals from the Ebola genome\n",
                "Let us use our function from the last section to simulate some read intervals along the Ebola genome."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "simulated_ebola_read_intervals = simulate_shotgun_sequencing_read_intervals(\n",
                "    genome_length = len(ebola_sequence),\n",
                "    num_copies = 10000,\n",
                "    k = 10,\n",
                "    min_fragment_length = 300,\n",
                "    num_fragments = 1000,\n",
                "    read_length = 100)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "A common feature to examine in shotgun sequencing data is the \"read coverage\" distribution, which is how many reads cover (overlap) each position along the genome.  Below is a function that you can use to quickly compute read coverage given a list of oriented read intervals."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def read_coverage(read_intervals, genome_length):\n",
                "    \"\"\"Computes the number of reads covering (overlapping) each position\n",
                "    along the genome.\n",
                "    \n",
                "    Args:\n",
                "        read_intervals: a list of oriented read intervals\n",
                "        genome_length: the integer length of the genome sequence\n",
                "    Returns:\n",
                "        A list of coverage values with the ith entry in the list giving\n",
                "        the read coverage at position i in the genome.\n",
                "    \"\"\"    \n",
                "    num_reads_starting = [0] * genome_length\n",
                "    num_reads_ending = [0] * genome_length\n",
                "    for start, end, strand in read_intervals:\n",
                "        num_reads_starting[start] += 1\n",
                "        num_reads_ending[end - 1] += 1\n",
                "    coverage = [0] * genome_length\n",
                "    coverage[0] = num_reads_starting[0]\n",
                "    for i in range(1, genome_length):\n",
                "        coverage[i] = coverage[i - 1] - num_reads_ending[i - 1] + num_reads_starting[i] \n",
                "    return coverage"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "And here is a function to plot the read coverage."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def plot_read_coverage(coverage):\n",
                "    import matplotlib.pyplot as plt\n",
                "    plt.figure(figsize=(16,4))\n",
                "    plt.plot(range(len(coverage)), coverage)\n",
                "    plt.xlabel('position along Ebola genome')\n",
                "    plt.ylabel('# of reads covering position')\n",
                "    plt.show()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Now let us compute and plot the read coverage across the genome."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "simulated_ebola_read_coverage = read_coverage(simulated_ebola_read_intervals, len(ebola_sequence))\n",
                "plot_read_coverage(simulated_ebola_read_coverage)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Examine the plot above and discuss the following questions with your neighbors:  \n",
                "* Why does it have this shape?\n",
                "* Try changing the parameters of the simulation slightly.  Which parameters have the greatest effect on the read coverage distribution?\n",
                "* What might be important about this distribution for the purposes of assembling the genome?\n",
                "Write your thoughts about these questions below.\n",
                "\n",
                "*WRITE YOUR THOUGHTS HERE*"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Problem 9 - Read intervals to read sequences (1 POINT)\n",
                "\n",
                "Although it has been convenient and efficent to work with read intervals, we will ultimately want to simulated actual read sequences, which will be the input to a genome assembly method.  For simplicity, we assume that the sequencer is error free.  Later, we will consider how sequencing errors might occur and complicate the assembly task."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Write a function that takes as input the genome sequence and a single oriented read interval and outputs the DNA  sequence of the read for the given interval and orientation.  You will need to use your `reverse_complement` function from the Day 1 notebook."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def sequence_for_read_interval(read_interval, genome_sequence):\n",
                "    \"\"\"Returns the sequence of a read given its oriented interval along\n",
                "    a genome sequence.\n",
                "    \n",
                "    Args:\n",
                "        read_interval: an oriented read interval, represented as a tuple (start, end, orientation)\n",
                "        genome_sequence: the genome sequence as a string\n",
                "    Returns:\n",
                "        The sequence of the read as a string.\n",
                "    \"\"\"    \n",
                "    \n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n",
                "    \n",
                "# You will likely want to use your reverse_complement function from day 1\n",
                "# which you can paste below\n",
                "###\n",
                "### YOUR CODE HERE\n",
                "###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "nbgrader": {
                    "grade": true,
                    "grade_id": "sequence_for_read_interval",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "scrolled": false,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for sequence_for_read_interval\n",
                "assert sequence_for_read_interval((0, 5, '+'), \"CGATGCTTAA\") == 'CGATG', \"Failed on (0, 5, '+')\"\n",
                "assert sequence_for_read_interval((0, 5, '-'), \"CGATGCTTAA\") == 'CATCG', \"Failed on (0, 5, '-')\"\n",
                "assert sequence_for_read_interval((4, 10, '+'), \"CGATGCTTAA\") == 'GCTTAA', \"Failed on (4, 10, '+')\"\n",
                "assert sequence_for_read_interval((4, 10, '-'), \"CGATGCTTAA\") == 'TTAAGC', \"Failed on (4, 10, '-')\"\n",
                "print(\"SUCCESS: sequence_for_read_interval passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "### Putting it all together: simulating shotgun sequencing read sequences\n",
                "\n",
                "For completeness, Below you are provided with a function that puts all of your work above together and simulates a set of read sequences for shotgun sequencing given a genome sequence."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "simulate_shotgun_sequencing_reads",
                    "locked": true,
                    "points": "1",
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def simulate_shotgun_sequencing_reads(genome_sequence, \n",
                "                                      num_copies, \n",
                "                                      k,\n",
                "                                      min_fragment_length,\n",
                "                                      num_fragments, \n",
                "                                      read_length):\n",
                "    \"\"\"Simulates shotgun sequence reads from a genome.\n",
                "    \n",
                "    Args:\n",
                "        genome_sequence: the genome sequence as a string\n",
                "        num_copies: the number of identical copies of the DNA in the sample\n",
                "        k: the number of cuts in each DNA during fragmentation\n",
                "        min_fragment_length: the minimum length of a fragment to be sequenced\n",
                "        num_fragments: the number of fragments to sequence\n",
                "        read_length: the length of the reads\n",
                "    Returns:\n",
                "        A list of read sequences.\n",
                "    \"\"\"    \n",
                "    \n",
                "    read_intervals = simulate_shotgun_sequencing_read_intervals(len(genome_sequence),\n",
                "                                                                num_copies,\n",
                "                                                                k,\n",
                "                                                                min_fragment_length,\n",
                "                                                                num_fragments,\n",
                "                                                                read_length)\n",
                "    reads = [sequence_for_read_interval(interval, genome_sequence)\n",
                "             for interval in read_intervals]\n",
                "    return reads"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "locked": true,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "Let's try simulating from the Ebola genome to see what sequencing data look like.  For example,"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "simulate_shotgun_sequencing_reads(\n",
                "    genome_sequence = ebola_sequence,\n",
                "    num_copies = 10000,\n",
                "    k = 10,\n",
                "    min_fragment_length = 300,\n",
                "    num_fragments = 10,\n",
                "    read_length = 30)"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3 [3.6]",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text\/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.6.4"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}