lectures/lecture13/lecture13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 13: Practical analyses in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dictionaries\n",
    "Python defines a powerful type called a dictionary or `dict`.\n",
    "\n",
    "### Mutable versus immutable types\n",
    "Before using dictionaries, we need to understand the difference between [mutable and immutable types](https://medium.com/@meghamohan/mutable-and-immutable-side-of-python-c2145cf72747).\n",
    "\n",
    "Let's review the types we've learned so far:\n",
    "#### Integers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "int"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_int = 1\n",
    "\n",
    "type(x_int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An integer is immutable, we can't change it's value by reassigment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_int = 2, x = 1\n"
     ]
    }
   ],
   "source": [
    "y_int = x_int\n",
    "y_int += 1\n",
    "print(f\"y_int = {y_int}, x = {x_int}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Floats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "float"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_float = 3.7\n",
    "\n",
    "type(x_float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Floats are also immutable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_float = 4.7, x = 3.7\n"
     ]
    }
   ],
   "source": [
    "y_float = x_float\n",
    "y_float += 1\n",
    "print(f\"y_float = {y_float}, x = {x_float}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "list"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_list = [1, 2, 3]\n",
    "\n",
    "type(x_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But lists are mutable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_list = [1, 2, 3, 4], x_list = [1, 2, 3, 4]\n"
     ]
    }
   ],
   "source": [
    "y_list = x_list\n",
    "y_list.append(4)\n",
    "\n",
    "print(f\"y_list = {y_list}, x_list = {x_list}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the value of `x_list` has changed too! Lists are mutable, meaning they point to the mutable object, not to its value.\n",
    "Be careful with mutable objects!\n",
    "\n",
    "There is an immutable object that has some similar properties to a list called a `tuple`, but we aren't going to go into those right now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "str"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_str = 'hello'\n",
    "\n",
    "type(x_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Strings are immutable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_str = hello friend, x_str = hello\n"
     ]
    }
   ],
   "source": [
    "y_str = x_str\n",
    "y_str += ' friend'\n",
    "\n",
    "print(f\"y_str = {y_str}, x_str = {x_str}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is a dictionary?\n",
    "Dictionaries are \"look up tables\" that can be used to map keys to values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'best_teacher': 'Jesse', 'worst_teacher': 'Rasi'}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict = {'best_teacher': 'Jesse',\n",
    "                'worst_teacher': 'Rasi'}\n",
    "\n",
    "teacher_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dictionary maps the best and worst teachers to their names.\n",
    "We can look up the value for a key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Jesse'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict['best_teacher']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we get an error if they key doesn't exist:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'second_best_teacher'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-c5a65f8e63d2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mteacher_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'second_best_teacher'\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: 'second_best_teacher'"
     ]
    }
   ],
   "source": [
    "teacher_dict['second_best_teacher']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get all of the keys:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['best_teacher', 'worst_teacher'])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_values(['Jesse', 'Rasi'])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict.values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the items, which are tuples with keys and values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_items([('best_teacher', 'Jesse'), ('worst_teacher', 'Rasi')])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict.items()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add entries to a dictionary like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'best_teacher': 'Jesse',\n",
       " 'worst_teacher': 'Rasi',\n",
       " 'most_fashionable_teacher': 'Trevor'}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict['most_fashionable_teacher'] = 'Trevor'\n",
    "\n",
    "teacher_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the above cell shows that dictionaries are mutable, as we changed the dictionary by adding to it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An advantage of a dictionary is that it is very fast to look up key / value pairs even if there are lots of entries in the dictionary.\n",
    "\n",
    "So why did I make a big deal about mutable versus immutable?\n",
    "Dictionary keys can **only** be immutable. \n",
    "So keying a dictionary with a mutable type (such as a `list`) is not allowed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unhashable type: 'list'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-f0904dbf7205>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mteacher_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'previous_teacher'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'current_teacher'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Erick'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Jesse'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'list'"
     ]
    }
   ],
   "source": [
    "teacher_dict[['previous_teacher', 'current_teacher']] = [['Erick', 'Jesse']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But dictionary values can be mutable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'best_teacher': 'Jesse',\n",
       " 'worst_teacher': 'Rasi',\n",
       " 'most_fashionable_teacher': 'Trevor',\n",
       " 'previous_and_current_teachers': [['Erick', 'Jesse']]}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teacher_dict['previous_and_current_teachers'] = [['Erick', 'Jesse']]\n",
    "\n",
    "teacher_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a dictionary in a function\n",
    "Now we will use a dictionary to implement a function that gets the reverse complement of a sequence.\n",
    "\n",
    "Here is the desired signature of the function. Can you write this function?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reverse_complement(seq):\n",
    "    \"\"\"Get reverse complement of a DNA sequence.\n",
    "    \n",
    "    Parameters\n",
    "    -----------\n",
    "    seq : str\n",
    "        Uppercase DNA sequence.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    str\n",
    "        Reverse complement of the sequence in upper case.\n",
    "        \n",
    "    Example\n",
    "    --------\n",
    "    >>> reverse_complement('ATGCAC')\n",
    "    'GTGCAT'\n",
    "    \n",
    "    \"\"\"\n",
    "    rc_dict = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G', 'N': 'N'}\n",
    "    revcomplement = []\n",
    "    for nt in reversed(seq.upper()):\n",
    "        revcomplement.append(rc_dict[nt])\n",
    "    return ''.join(revcomplement)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the function is written, you can use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'GTGCAT'"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reverse_complement('ATGCAC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'NTTACCTGTATCCTGAACTTTTAGCAT'"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reverse_complement('aTGCTAAAAGTTCAGGATACAGGTAAN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regular expressions\n",
    "The `re` module is for \"regular expressions\". \n",
    "These are very useful for parsing strings.\n",
    "\n",
    "Here is a common example from my work with influenza.\n",
    "You download some strains from the database, and they have names that look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "strain1 = 'A/New York/3/1994 (H3N2)'\n",
    "strain2 = 'A/California/3/X/2003 (H12N1)'\n",
    "strain3 = 'A/Perth/2009 (H3N2)'\n",
    "\n",
    "strains = [strain1, strain2, strain3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You want to get some information out of these, like the year.\n",
    "\n",
    "Regular expressions allow you to define patterns in a string.\n",
    "Here we build a regular expression that gets the subtype out and then use a dictionary to count how many sequences there are of each subtype:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'H3N2': 2, 'H12N1': 1}\n"
     ]
    }
   ],
   "source": [
    "import re\n",
    "\n",
    "# regular expression that gets year and subtype\n",
    "strainmatch = re.compile(\n",
    "        '(?P<year>\\d{4}) \\((?P<subtype>H\\d+N\\d+)\\)$')\n",
    "\n",
    "subtype_counter = {}  # dict to store the results\n",
    "\n",
    "for strain in strains:  # loop over all strains\n",
    "    m = strainmatch.search(strain)\n",
    "    year = m.group('year')\n",
    "    subtype = m.group('subtype')\n",
    "    if subtype in subtype_counter:\n",
    "        subtype_counter[subtype] += 1\n",
    "    else:\n",
    "        subtype_counter[subtype] = 1\n",
    "        \n",
    "print(subtype_counter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are lots of handy special codes in the Python regular expression module (see [here](https://docs.python.org/3.7/library/re.html)), and you can use them to do almost any type of string matching."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using regular expressions to parse barcodes\n",
    "Now we will use regular expressions to parse barcodes from nucleotide sequences.\n",
    "For instance, you might have to do this in a single-cell RNA-seq experiment where there is a barcode at the end of each read telling you the cell that the read came from.\n",
    "\n",
    "Imagine that our valid molecules should have sequences like this:\n",
    "\n",
    "`CTAGCNNNNNNGATCA`\n",
    "\n",
    "See how there is a 6-nucleotide barcode in the center of the sequence.\n",
    "We have a list of sequences, and want to parse through them to figure out which ones meet the expected pattern--and get the barcode from those that do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "seqs = ['CTAGCatcgatGATCA',  # has barcode ATCGAT\n",
    "        'CCAGCatagcaGATCA',  # does not have expected 5' sequence\n",
    "        'CTAGCtacagGATCA',   # barcode too short\n",
    "        'CTAGCgaccatGATCA',  # has barcode GACCAT\n",
    "        'CTAGCatcgatGATCA',  # has barcode ATCGAT\n",
    "        'CTAGCatcgatGGTCA',  # does not have expected 3' sequence\n",
    "        ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function that parses these barcoded sequences and gets the ones with valid barcodes.\n",
    "In doing this, note that:\n",
    "\n",
    "  1. If you have a string `s`, `s.upper()` makes it all uppercase.\n",
    "  2. Consider using the following regular expression options (all described [here](https://docs.python.org/3.7/library/re.html):\n",
    "    - `[]` indicates a set of characters.\n",
    "    - `{n}` indicates number of occurrences.\n",
    "    - `(?P<name>)` indicates a group with name `name`.\n",
    "    - The `match` method finds a match at the beginning of the string.\n",
    "    \n",
    "Below I've written the function documentation, try to implement it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_barcodes(seqs, bclen=6, upstream='CTAGC', downstream='GATCA'):\n",
    "    \"\"\"Parse and count barcodes.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    seqs : list\n",
    "        DNA sequences.\n",
    "    bclen : int\n",
    "        Length of barcode\n",
    "    upstream : str\n",
    "        Sequence upstream of barcode.\n",
    "    downstream : str\n",
    "        Sequence downstream of barcode.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Keyed by each valid barcode, value is number of times the barcode\n",
    "        is observed.\n",
    "        \n",
    "    Note\n",
    "    ----\n",
    "    The function is **not** case-sensitive, and all barcodes are reported\n",
    "    in upper-case.\n",
    "    \n",
    "    \"\"\"\n",
    "    regex = re.compile(upstream +\n",
    "                       f\"(?P<bc>[ACGT]{{{bclen}}})\" +\n",
    "                       downstream\n",
    "                       )\n",
    "    counts = {}\n",
    "    for seq in seqs:\n",
    "        seq = seq.upper()\n",
    "        m = regex.match(seq)\n",
    "        if m:\n",
    "            bc = m.group('bc')\n",
    "            if bc in counts:\n",
    "                counts[bc] += 1\n",
    "            else:\n",
    "                counts[bc] = 1\n",
    "    return counts\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the function once you've implemented it. Does it give the right result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ATCGAT': 2, 'GACCAT': 1}"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_barcodes(seqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Biopython\n",
    "[Biopython](https://biopython.org/) is a package that has lots of useful functions for computational biology.\n",
    "\n",
    "It is very handy for things like reading in sequences in many different formats: [Bio.SeqIO](https://biopython.org/wiki/SeqIO) is your friend!\n",
    "\n",
    "(Do note that if you are analyzing truly large datasets, `Biopython` is not very fast and you may want to use something like [pysam](https://pysam.readthedocs.io/en/latest/api.html); but `Biopython` is a good starting point).\n",
    "\n",
    "### Reading in a file\n",
    "I have included the file [barcodes_R1.fastq](barcodes_R1.fastq), which has some FASTQ sequences in it.\n",
    "\n",
    "First, let's just see what the beginning of that file looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! head -n 8 barcodes_R1.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use `Biopython` to read the FASTQ entries.\n",
    "\n",
    "First, import `Biopython.SeqIO`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Bio.SeqIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now read in the sequencing reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seqreads = list(Bio.SeqIO.parse('barcodes_R1.fastq', format='fastq'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many reads were there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Found {len(seqreads)} sequencing reads.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the first read:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seqreads[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that it has a lot of information, including the header, quality scores, etc.\n",
    "\n",
    "For our purposes, we will just convert the sequence part into a string for each sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seqreads_str = [str(s.seq) for s in seqreads]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure we still have the same number of sequencing reads, and look at the first one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(seqreads_str) == len(seqreads)\n",
    "\n",
    "seqreads_str[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A real biological analysis: parsing barcodes\n",
    "The reads that we just read as `seqreads_str` come from a real sequencing run of influenza virus HA and NA genes.\n",
    "\n",
    "The sequences are as follows:\n",
    "\n",
    "    5'-[end of HA]-AGGCGGCCGC-[16 X N]-3'\n",
    "    \n",
    "or \n",
    "\n",
    "    5'-[end of NA]-AGGCGGCCGC-[16 X N]-3'\n",
    "    \n",
    "The end of NA is:\n",
    "\n",
    "    ...CACGATAGATAAATAATAGTGCACCAT\n",
    "    \n",
    "The end of HA is:\n",
    "\n",
    "    ...CCGGATTTGCATATAATGATGCACCAT\n",
    "    \n",
    "The sequencing run reads from the reverse end of the molecules, so the first thing in the sequencing reads is the barcode followed by the constant sequence and the end of HA or NA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to determine which reads have valid sequences, get the barcodes out of strings, figure out if the barcode matches to HA or NA, and count the barcodes.\n",
    "So this requires setting up an analysis that does the following:\n",
    "\n",
    " 1. Get the reverse complement of each read.\n",
    " 2. Determine if it matches the expected pattern for HA and NA, and if so which one.\n",
    " 3. If it matches, extract the barcode and add it to a dictionary to keep track of counts.\n",
    " \n",
    "Can you write code that does this?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "623.8px",
    "left": "0px",
    "right": "829.6px",
    "top": "33.2px",
    "width": "322.4px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}