lectures/lecture08/lecture08.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 08: Practical analyses in Python\n",
    "\n",
    "Last week Phil gave a great intro to Python!\n",
    "This week we are going to introduce a few more concepts that will help you make your code more sophisticated (and flexible)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining new functions\n",
    "We use the `def` keyword to define a new function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hello!\n"
     ]
    }
   ],
   "source": [
    "def say_hello():\n",
    "    print('Hello!')\n",
    "\n",
    "s = say_hello()\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nice to meet you, phil\n",
      "Nice to meet you, charlie\n"
     ]
    }
   ],
   "source": [
    "# functions can take inputs, called \"arguments\" or \"parameters\"\n",
    "def greet(name):\n",
    "    print('Nice to meet you,', name)\n",
    "\n",
    "greet('phil')\n",
    "greet('charlie')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions can return results to the caller, using the `return` statement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reverse_complement(seq):\n",
    "    \"\"\"returns the reverse complement of a nucleic acid sequence\"\"\"\n",
    "    base_partner = {'A':'T', 'T':'A', 'C':'G', 'G':'C'}\n",
    "    bwd = ''\n",
    "    # iterate through all bases in the sequence\n",
    "    for base in seq:\n",
    "        # look up the complementary base in the dictionary\n",
    "        pair = base_partner[base]\n",
    "        # add the complentary base to the beginning of the string (reverse comp)\n",
    "        bwd = pair + bwd\n",
    "    return bwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fwd= ACGGTAATGATCCTCAG rev= CTGAGGATCATTACCGT\n"
     ]
    }
   ],
   "source": [
    "fwd = 'ACGGTAATGATCCTCAG'\n",
    "rev = reverse_complement( fwd )\n",
    "\n",
    "print('fwd=',fwd,'rev=',rev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# an assert statement checks to see if something is true, and will interrupt the program if it's not.\n",
    "assert reverse_complement(reverse_complement(fwd)) == fwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions can have OPTIONAL ARGUMENTS whose DEFAULT VALUES are pre-specified in the function definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reverse_complement(seq, unk_partner='N'):\n",
    "    \"\"\"Returns the reverse complement of a nucleic acid sequence\n",
    "    \n",
    "    Uses unk_partner as the partner of unrecognized letters\n",
    "    \"\"\"\n",
    "    base_partner = {'A':'T', 'T':'A', 'C':'G', 'G':'C'}\n",
    "    rseq = ''\n",
    "    for a in seq:\n",
    "        if a in base_partner:\n",
    "            # look up the complementary base in the dictionary\n",
    "            pair = base_partner[a]\n",
    "            rseq = pair + rseq\n",
    "        else:\n",
    "            rseq = unk_partner + rseq\n",
    "    return rseq\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GTCGNANTCNGCTACAGT\n",
      "GTCG-A-TC-GCTACAGT\n"
     ]
    }
   ],
   "source": [
    "fwd = 'ACTGTAGCxGAcTNCGAC'\n",
    "\n",
    "print(reverse_complement(fwd))\n",
    "print(reverse_complement(fwd, unk_partner='-'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions can call other functions, even themselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 factorial = 1\n",
      "1 factorial = 1\n",
      "2 factorial = 2\n",
      "3 factorial = 6\n",
      "4 factorial = 24\n",
      "5 factorial = 120\n",
      "6 factorial = 720\n",
      "7 factorial = 5040\n",
      "8 factorial = 40320\n",
      "9 factorial = 362880\n"
     ]
    }
   ],
   "source": [
    "def factorial(n):\n",
    "    \"\"\"Calculate the factorial of a number recursively. Bad things will happen \n",
    "    if the number is negative or not an integer \"\"\"\n",
    "    if n==0:\n",
    "        return 1\n",
    "    else:\n",
    "        return n * factorial(n-1) # this is called \"recursion\"\n",
    "\n",
    "for i in range(10):\n",
    "    print(i,'factorial =',factorial(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try using `help(factorial)`, `factorial?`, and `factorial??` to see the docstring and source code of our new function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function factorial in module __main__:\n",
      "\n",
      "factorial(n)\n",
      "    Calculate the factorial of a number recursively. Bad things will happen \n",
      "    if the number is negative or not an integer\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(factorial)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0;31mSignature:\u001b[0m \u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mDocstring:\u001b[0m\n",
      "Calculate the factorial of a number recursively. Bad things will happen \n",
      "if the number is negative or not an integer \n",
      "\u001b[0;31mFile:\u001b[0m      /tmp/ipykernel_509/2221558712.py\n",
      "\u001b[0;31mType:\u001b[0m      function\n"
     ]
    }
   ],
   "source": [
    "factorial?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0;31mSignature:\u001b[0m \u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mSource:\u001b[0m   \n",
      "\u001b[0;32mdef\u001b[0m \u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0;34m\"\"\"Calculate the factorial of a number recursively. Bad things will happen \u001b[0m\n",
      "\u001b[0;34m    if the number is negative or not an integer \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0;32mif\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m        \u001b[0;32mreturn\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m        \u001b[0;32mreturn\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# this is called \"recursion\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFile:\u001b[0m      /tmp/ipykernel_509/2221558712.py\n",
      "\u001b[0;31mType:\u001b[0m      function\n"
     ]
    }
   ],
   "source": [
    "factorial??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Practice time:\n",
    "\n",
    "Take a few minutes to write a function that uses a dictionary to count the number of times each base occurs.\n",
    "Be sure to account for sequences that contain uppercase and lowercase bases.\n",
    "\n",
    "I've provided an empty function here (with some docstrings already written). Use the provided test cases check your work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_bases(seq):\n",
    "    \"\"\"Count the number of times each base occurs in the sequence.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    seq : string\n",
    "        DNA sequence.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Keyed by each nucleotide, value is number of times the nucleotide\n",
    "        is observed.\n",
    "    \n",
    "    \"\"\"\n",
    "    base_counts = {}\n",
    "    \n",
    "    seq = seq.upper()\n",
    "    \n",
    "    for base in seq:\n",
    "        if base in base_counts:\n",
    "            base_counts[base] += 1\n",
    "        else:\n",
    "            base_counts[base] = 1\n",
    "    \n",
    "    return base_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': 2, 'T': 2, 'C': 2, 'G': 2}"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_bases('AATCGGCT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': 2, 'T': 3, 'G': 2, 'C': 1}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_bases('aatTGGcT')"
   ]
  },
  {
   "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",
    "A regular expression is a sequence of characters that forms a search pattern.\n",
    "They can be used to check if a string contains the specified search pattern.\n",
    "\n",
    "The `re` package offers a set of functions that allows us to search a string for a match. \n",
    "The main functions we will be using are called `search` which allows us to find string matches at any position in the string.\n",
    "We can search for a pattern in a string like this:\n",
    "\n",
    "```\n",
    "import re\n",
    "string = \"This is an example string\"\n",
    "\n",
    "# compile the search pattern\n",
    "search_pattern = re.compile(\"search pattern here\")\n",
    "\n",
    "# search for the search pattern in the string\n",
    "search_pattern.search(string)\n",
    "```\n",
    "\n",
    "Here are some common elements to have in your search pattern:\n",
    "* letter characters which returns a match where the string contains the specified letter (e.g. `A`, `B`, `C`, ...)\n",
    "* special characters which returns a match where the string contains the specified special character; these must be preceded by a `\\`\n",
    "* `\\d` which returns a match where the string contains digits (numbers from `0`-`9`)\n",
    "\n",
    "You may also want to add the following customizations:\n",
    "* `[]` specifies a set of characters to search for (e.g. `[a-n]`)\n",
    "* `()` capture and group everything contained inside, and search the string for everything together\n",
    "* `?P<name>` indicates a search pattern group with name `name`\n",
    "* `+` specifies one or more occurrences of a certain pattern element\n",
    "* `{}` specifies exactly the specified number of occurrences of a certain pattern element\n",
    "* `$` specifies the end of the string\n",
    "* `^` specifies the beginning of the string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Here is a common example dealing with influenza.\n",
    "You download some strains from the database, and they have names that look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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 subtype.\n",
    "Let's build a regular expression that gets the subtype out of `strain2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<re.Match object; span=(22, 29), match='(H12N1)'>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import re\n",
    "\n",
    "# compile a re for the subtype\n",
    "strainmatch = re.compile(\n",
    "        '\\(H\\d+N\\d+\\)')                        \n",
    "\n",
    "# search for the search pattern in the string\n",
    "strainmatch.search(strain2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's extend that a bit and build a regular expression that gets only the subtype out and then use a dictionary to count how many sequences there are of each subtype:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'H3N2': 2, 'H12N1': 1}\n"
     ]
    }
   ],
   "source": [
    "# compile a re for the subtype (with a named search pattern)\n",
    "strainmatch = re.compile(\n",
    "        '\\((?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",
    "    # search for re in each strain\n",
    "    m = strainmatch.search(strain)\n",
    "    # isolate named pattern\n",
    "    subtype = m.group('subtype')\n",
    "    # add pattern/count to dictionary\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.\n",
    "\n",
    "I like to test my regular expression calls using [this website](https://regex101.com)"
   ]
  },
  {
   "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": null,
   "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. The `match` function finds a match at the beginning of the string (whereas the `search` function we were using before searches for matches throughout the whole string)\n",
    "    \n",
    "Below I've written the function documentation, try to implement it.\n",
    "__Take a few minutes in groups to work through this.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "    \n",
    "    # your code here ...\n",
    "    \n",
    "    return None\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_barcodes(seqs)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "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.9.12"
  },
  "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
  },
  "vscode": {
   "interpreter": {
    "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}