automation_python_instructor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Outline of Defensive Programming and Automation\n",
    "\n",
    "We've gone over the basics of Python and Bash, and we've done some data exploration using the Jupyter Notebook. This is a good entrypoint to discuss some programming *practices*: how to program in ways that are understandable/transparent, will save you time, and will prevent headaches down the road.\n",
    "\n",
    "Outline:\n",
    "1. **Defensive Programming** - defending yourself and others against \"bad\" code\n",
    "2. **Working with files** - allowing our Jupyter/Python functions to interact with files outside this notebook\n",
    "3. **Scripts** - saving our code for later\n",
    "4. **Automation** - letting the computer run everything for you\n",
    "\n",
    "Code is read much more often than it is written. Ensuring your original code when first written is documented and efficient means less headscratching when you inevitably revisit it to re-do analyses, add more functions, use other datasets, and many other real scenarios that you'll probably experience in lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Defensive Programming\n",
    "\n",
    "It is very easy to start coding, and then subsequently lose track of why and how you wrote your code. You may also have collaborators with whom you share your code, you may have auditors or reviewers that need to check your work, or you may even step away for your computer for the day and forget what you were doing! \"Defensive programming\" is a set of tactics that protects you and others from code that is difficult to read, understand, or troubleshoot. There are many components to defensive programming, some of which we'll cover in today's module:\n",
    "* Style guides\n",
    "* Naming conventions\n",
    "* Comments\n",
    "* Docstrings\n",
    "* Don't Repeat Yourself - DRY\n",
    "* Error handling\n",
    "* Outlining your code - Pseudocode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Style Guides\n",
    "\n",
    "Programming languages have style guides that ensure your code is readable, usable, and debugable by others. Most Python developers choose to follow the PEP8 style guide: https://www.python.org/dev/peps/pep-0008/. Other languages, such as R, have their own style guides (e.g. Tidyverse: https://style.tidyverse.org/). Many automated tools rely on these standards to parse your code for debugging and generating documentation. We won't go through all of PEP8, but we'll start with a subset of conventions that are common across different languages.\n",
    "\n",
    "Adopting a style guide provides **consistency** to your code, making it easier for you and others to read and edit. That being said, while the above guides are commonly used, they aren't required. If you have a reason to format something differently, as long as others can understand it and you've documented the reason, go for it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Naming Conventions\n",
    "\n",
    "How do you decide which names to give to your variables, functions, and scripts? In general, they should be **descriptive**: the name you give a variable should tell you about what information it holds. For example:\n",
    "* A variable named `path_to_fasta` would tell you that it stores the location of a [FASTA file](https://en.wikipedia.org/wiki/FASTA_format) on your computer, probably as a `String` object.\n",
    "* A variable named `phred_score` would tell you that it stores a [PHRED score](https://en.wikipedia.org/wiki/Phred_quality_score), probably as an `int`.\n",
    "* A variable named `temp2` probably doesn't tell you anything about what data it holds.\n",
    "\n",
    "There are many common conventions as to how variable names should be formatted. The ones below are not specific to Python; you'll see them across other languages:\n",
    "\n",
    "1. **snake_case**: All words are lower case with underscores between them\n",
    "2. **CamelCase**: Words start with capital letters and are not separated\n",
    "3. **mixedCase**: Like CamelCase but the first word is lowercase\n",
    "4. **UPPERCASE_WITH_UNDERSCORES**: All letters are uppercase, separated by underscores \n",
    "\n",
    "Python's style guide outlines rules for using naming conventions:\n",
    "\n",
    "1. Variables: **snake_case**\n",
    "    - dna_sequence, file_path\n",
    "2. Functions: **snake_case**\n",
    "    - combine_replicates(), align_sequences()\n",
    "3. Errors: **CamelCase**\n",
    "    - ValueError, SyntaxError, FileNotFoundError\n",
    "\n",
    "There are more guidelines, but these are common ones that are encountered early on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Using naming conventions\n",
    "\n",
    "Edit the code block below to conform to PEP8 naming conventions. Post in the collaborative document your answers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Velocity(TOTALDISTANCE, time):\n",
    "    \"This calculates the distance over time\"\n",
    "    Velocity_Result = TOTALDISTANCE / time\n",
    "    return(Velocity_Result)\n",
    "Velocity(10, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dangers in variable naming\n",
    "\n",
    "Aside from making names look professional, there are some rules about naming to help prevent errors and bugs. Specifically, you should never name something the same as a function included in python. \n",
    "\n",
    "Let's use the function sum() to see what happens of we use sum as a variable name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_of_two_numbers = sum([5, 4])\n",
    "print(\"Sum data type\", type(sum))\n",
    "sum = 10 + 5\n",
    "print(\"Sum data type\", type(sum))\n",
    "print(sum_of_two_numbers)\n",
    "print(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try calling the sum() function again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum([5, 4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments\n",
    "\n",
    "Comments are human-readable annotations that document various lines/functions in your code document. In Python, they are denoted using `#` - everything written after it on the same line is not run. This means we can use it to help both ourselves and others understand what the code is trying to do. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_function(x):\n",
    "    # This is a comment.\n",
    "    # print(x + x)\n",
    "    # The print function above will not run due to the '#'\n",
    "    print(x)\n",
    "\n",
    "my_function(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should write comments for your code often, specifically when you are doing a task that is specialized like using a formula or applying a custom function to do a task. Some good advice from PEP8:\n",
    "\n",
    ">Comments that contradict the code are worse than no comments. Always make a priority of keeping the comments up-to-date when the code changes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Docstrings\n",
    "\n",
    "We previously covered docstrings within functions as a string just after the def statement. However, what if you want to write **more than just a single line of documentation**? Fortunatly there is a way to do so by using **triple-quotes**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "def my_function(x):\n",
    "    \"\"\"Docstring for my_function\n",
    "    \n",
    "    Prints the input `x` to the console.\n",
    "    \"\"\"\n",
    "    # This is a comment.\n",
    "    # print(x + x)\n",
    "    # The print function above will not run due to the '#'\n",
    "    print(x)\n",
    "\n",
    "my_function(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PEP guidelines on docstrings\n",
    "\n",
    "Python PEP guidelines suggest the following format:\n",
    "\n",
    "\"\"\"One line description\n",
    "\n",
    "More details about your function, in triple-quotes\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "**or**\n",
    "\n",
    "\"\"\"Only a single line description in triple-quotes\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... are supplemented by community formats\n",
    "\n",
    "Docstrings are pretty flexible, even using PEP standards. There are a couple common format guidelines for docstrings that you can choose from. Why start with these?\n",
    "\n",
    "1. It gives a good overview of what information people expect in your docstrings\n",
    "2. The formats here can be parsed by common tools\n",
    "\n",
    "While the formats for documentation may differ slightly depending on the language choice, the information expected from them is fundamental. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Google format\n",
    "\"\"\"Takes a string and returns a list of letters\n",
    "\n",
    "Args:\n",
    "    string (list): A string to parse for letters\n",
    "    upper (bool): The letters are returned uppercase \n",
    "        (default is False)\n",
    "\n",
    "Returns:\n",
    "    list: A list of each letter in the string\n",
    "\"\"\"\n",
    "\n",
    "#Numpy format\n",
    "\"\"\"Takes a string and returns a list of letters\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "string : str\n",
    "    A string to parse for letters\n",
    "upper : bool, optional\n",
    "    The letters are returned uppercase  (default is False)\n",
    "\n",
    "Returns\n",
    "-------\n",
    "list\n",
    "    A list of each letter in the string\n",
    "\"\"\"\n",
    "\n",
    "#reStrucured text \n",
    "\"\"\"Takes a string and returns a list of letters\n",
    "\n",
    ":param string: A string to parse for letters\n",
    ":type string: str\n",
    ":param upper: A string used to join each string (default is False)\n",
    ":type upper: bool\n",
    ":returns: A list of each letter in the string\n",
    ":rtype: list\n",
    "\"\"\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Importance of naming and documentation\n",
    "\n",
    "Given the following function with poor naming and no documentation, determine:\n",
    "1. What are the 2 inputs\n",
    "2. What does it return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FUNCTION(number, words):\n",
    "    Smallest = 0\n",
    "    for dictionary in number:\n",
    "        LETTER = (dictionary / words) * 100\n",
    "        if LETTER > Smallest:\n",
    "            Smallest = LETTER\n",
    "    return Smallest\n",
    "\n",
    "# Answers:\n",
    "# What are the two inputs:\n",
    "#     * number = list of int\n",
    "#     * word = int\n",
    "# What does the function return:\n",
    "#     * Maximum percent of elements in list “number”, from the total “words”"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Refactoring a Function\n",
    "\n",
    "**Refactoring** is a term that means re-writing code without changing the task it performs. \n",
    "\n",
    "Refactor the above function with poor naming and no documentation. You will want to:\n",
    "\n",
    "1. Rename variables to an appropriate name\n",
    "2. Write a docstring explaining what the function does using one of the example formats (Google, Numpy, reStructured)\n",
    "\n",
    "Test your function by running it after you refactor to see if it still produces the same output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your answer here!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An advanced answer: simplifies the function using list comprehension (probably shouldn't teach it at this point)\n",
    "def greatestPC(num_list, num2):\n",
    "    return max([x / num2 * 100 for x in num_list])\n",
    "\n",
    "greatestPC([1, 2, 6, 4, 9, 1, 35, 51], 159)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Don't Repeat Yourself - DRY\n",
    "\n",
    "DRY is a concept in programming to avoid writing redundant code. A good sign your code is redundant would be if you copy and paste parts of it and then edit the new copies.\n",
    "\n",
    "Suppose we had three lists of proteins and wanted to check if there are any matches to a list of proteins of interest and wrote the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_data1 = [\"CREG1\", \"ELK1\", \"SF1\", \"GATA1\", \"GATA3\", \"CREB1\"]\n",
    "protein_data2 = [\"ATF1\", \"GATA1\", \"STAT3\", \"P53\", \"CREG1\"]\n",
    "protein_data3 = [\"RELA\", \"MYC\", \"SF1\", \"CREG1\", \"GATA3\", \"ELK1\"]\n",
    "proteins_of_interest = [\"ELK1\", \"MITF\", \"KAL1\", \"CREG1\"]\n",
    "\n",
    "# Are there any matches in the first list?\n",
    "match_list1 = []\n",
    "for protein in protein_data1:\n",
    "    if protein in proteins_of_interest:\n",
    "        match_list1.append(protein)\n",
    "\n",
    "# Are there any matches in the second?\n",
    "match_list2 = []\n",
    "for protein in protein_data2:\n",
    "    if protein in proteins_of_interest:\n",
    "        match_list2.append(protein)\n",
    "\n",
    "# Are there any matches in the third?\n",
    "match_list3 = []\n",
    "for protein in protein_data3:\n",
    "    if protein in proteins_of_interest:\n",
    "        match_list3.append(protein)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to do an analysis on multiple datasets. But there are a few problems with this approach:\n",
    "\n",
    "1. You have to copy/paste this code for each list you want to compare.\n",
    "2. If you make a change to the analysis, you need to edit every copy.\n",
    "    - Very hard to remain consistent\n",
    "3. If you want to compare to another list, you need to either:\n",
    "    - Edit every copy.\n",
    "    - Change the variable `proteins_of_interest`.\n",
    "        - May affect your analysis somewhere else\n",
    "\n",
    "To prevent repeating code, we can write a function to do our repeated task. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Refactor an Analysis Into a Function\n",
    "\n",
    "Write a function that does the repeated task and run it on the three lists.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your function here!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One potential answer\n",
    "def match_protein(protein):\n",
    "    \"\"\" Takes a list of str and return a list of str that contains proteins that match with proteins of interest\n",
    "    Args:\n",
    "        protein (list): list of str containing protein data\n",
    "        \n",
    "    Returns:\n",
    "        list: list of str containing proteins that matche proteins of interest\n",
    "    \"\"\"\n",
    "    match_list = []\n",
    "    for p in protein:\n",
    "        if p in proteins_of_interest:\n",
    "            match_list.append(p)\n",
    "    return match_list\n",
    "\n",
    "# Alternate function using list comprehension\n",
    "def match_protein(protein):\n",
    "    return [p for p in protein if p in proteins_of_interest]\n",
    "\n",
    "match_list1 = match_protein(protein_data1)\n",
    "match_list2 = match_protein(protein_data2)\n",
    "match_list3 = match_protein(protein_data3)\n",
    "print(match_list1)\n",
    "print(match_list2)\n",
    "print(match_list3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Helpful Way to Format Strings\n",
    "\n",
    "In the next sections, we will learn how to make custom error messages. Having a good way to format strings for these messages will make our job easier. Previously, we concatenated strings using the `+` operator. \n",
    "\n",
    "Suppose we have an integer, and we want to add it to a string:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Adding an int to a string\n",
    "my_int = 5\n",
    "message = \"This number does not work:\"\n",
    "message + my_int"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we directly add a string and an integer, we get a `TypeError`. One way to neatly get around this is to use what is called an **f string**. This is a python-specific method that takes care of formating the string for us.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example with f strings\n",
    "m = f\"This number does not work: {my_int}. The number is lower than 10.\"\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error Handling in Python\n",
    "\n",
    "One fundamental concept in programming is troubleshooting errors. A wide variety of errors can be reported when something goes wrong in Python. For example, look at the error message when the following code is run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(hello)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in the command `print(hello)` when `hello` is not defined, Python returns something called a `NameError`. Python has many different types of errors built-in including:\n",
    "\n",
    "- NameError\n",
    "- ZeroDivisionError\n",
    "- TypeError\n",
    "- IndexError\n",
    "- ... and more!\n",
    "\n",
    "These errors are also referred to as **exceptions**, meaning \"exceptional events\". More information on the other exceptions built into Python are here: https://docs.python.org/3/library/exceptions.html.\n",
    "\n",
    "We can also call an error on purpose using the `raise` keyword with a custom message. This message is part of the **traceback**, which is a report Python gives on what happened.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Raising a custom error\n",
    "raise NameError(\"My custom error message\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For errors we get without directly calling `raise`, both the class of `Error` and the message are predetermined. This is useful in some situations, but less so in others. Suppose you wrote a reverse complement function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# live demo\n",
    "def reverse_complement(dna_sequence):\n",
    "    \"\"\"Reverses the complement of a dna sequence\"\"\"\n",
    "    complements = {\"T\":\"A\", \"A\":\"T\", \"C\":\"G\", \"G\":\"C\"}\n",
    "    reverse = dna_sequence[::-1]\n",
    "    result = \"\"\n",
    "    for letter in reverse:\n",
    "        result = result + complements[letter]\n",
    "    return(result)\n",
    "\n",
    "help(reverse_complement)\n",
    "print(reverse_complement(\"CAAT\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and suppose someone imports this function from a script we wrote (which we will do later today). Maybe they will run our script with a sequence that contains lowercase letters (indicating masked genomic regions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reverse_complement(\"CACGtgcatggTGAAA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a user, this is a really confusing error. Instead of returning the reverse complement, the function raises an error saying `KeyError: 'g'`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Try and Except Keywords\n",
    "\n",
    "One way we can pre-empt errors is though `try` and `except`. A `try-except` code block runs the code under `try`, and if there is a particular error we think would happen, we can write a custom response under `except`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Try-Except demonstration\n",
    "try:\n",
    "    print(hello)\n",
    "except NameError:\n",
    "    print(\"A custom error message\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `try-except` block works as follows:\n",
    "1. The code under `try` is run.\n",
    "2. If no exception occurs after running all of the code under `try`, the code under `except` is skipped.\n",
    "3. If an exception occurs, *and* it matches the one in the `except` line, then any remaining code under `try` is skipped and the code under `except` is now run.\n",
    "4. If an exception occurs, but it *doesn't* match the exception we gave it, this exception is **unhandled** and the program stops.\n",
    "\n",
    "Note that, even though we **caught** the `NameError`, the program actually keeps going. This is because the `try-except` block only **controls the flow** of the program by specifying what code should be run in the case an error is encountered; it doesn't `raise` an exception by itself. We can visualize this by adding a `print()` statement after the whole `try-except` block:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example with print() before and after - continue from previous block\n",
    "print(\"Before\")\n",
    "try:\n",
    "    print(hello)\n",
    "except NameError:\n",
    "    print(\"A custom error message\")\n",
    "print(\"After\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want the program to stop, we can use `raise` within the code run under `except`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example with raise added within the except code - continue from previous block\n",
    "print(\"Before\")\n",
    "try:\n",
    "    print(hello)\n",
    "except NameError:\n",
    "    print(\"A custom error message\")\n",
    "    raise\n",
    "print(\"After\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, why go through the trouble of **catching exceptions**?\n",
    "\n",
    "1. If our program receives input that causes an error, we want the program to **fail fast**: don't do anything else, and immediately notify the user.\n",
    "2. We want to avoid returning **incorrect or unexpected** results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Printing Error Messages\n",
    "\n",
    "Before we use our new keywords, let's go over the output for errors. When we print something in Jupyter or in the command line, by default the text goes to a destination called **standard output**, or `stdout`. Simultaneously, errors and warnings go to a destination called **standard error**, or `stderr`. It is dangerous to print error messages in `stdout` because some workflows utilize it for data and the error messages can get mixed in (eg, when using pipes `|` on the command line).\n",
    "\n",
    "We can set a destination for out `print()` command using the `file` argument. The default destination is `sys.stdout`. To send to `sys.stderr`, we will need to import the `sys` module.\n",
    "\n",
    "### A brief digression on modules\n",
    "* `sys` is from the **standard library**, a collection of modules included in Python but not available by default.\n",
    "* Libraries are incredibly useful: there are libraries for working with numeric and scientific data, generating plots, fetching data from the web, working with image and document files, databases, and more.\n",
    "* The full documentation for `sys` is part of the general Python documentation: https://docs.python.org/3/library/sys.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Help for print\n",
    "help(print)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Import sys, compare stdout vs stderr\n",
    "import sys\n",
    "\n",
    "print(\"Hello this is going to stdout\", file = sys.stdout)\n",
    "print(\"Hello this is going to stderr\", file = sys.stderr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Jupyter, `stderr` has a red background and `stdout` has a no background color. \n",
    "\n",
    "**Note**: The **traceback** messages shown when an error occurs are going to `stderr`. In Jupyter, they are formatted differently than other output to `stderr`. \n",
    "\n",
    "Back to the observation that reverse_complement can give cryptic error messages: let's add a `try` statement and place the `for` loop in it. We can then use `except` with our expected `KeyError` to print a different message and `raise` to both **fast fail** and return the full **traceback**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit the function below!\n",
    "def reverse_complement(dna_sequence):\n",
    "    \"\"\"Reverses the complement of a dna sequence\"\"\"\n",
    "    complements = {\"T\":\"A\", \"A\":\"T\", \"C\":\"G\", \"G\":\"C\"}\n",
    "    reverse = dna_sequence[::-1]\n",
    "    result = \"\"\n",
    "    # Add try - except - raise statements\n",
    "    for letter in reverse:\n",
    "        result = result + complements[letter]\n",
    "    return(result)\n",
    "\n",
    "help(reverse_complement)\n",
    "print(reverse_complement(\"CAAg\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer to the above exercise:\n",
    "def reverse_complement(dna_sequence):\n",
    "    \"\"\"Reverses the complement of a dna sequence\"\"\"\n",
    "    complements = {\"T\":\"A\", \"A\":\"T\", \"C\":\"G\", \"G\":\"C\"}\n",
    "    reverse = dna_sequence[::-1]\n",
    "    result = \"\"\n",
    "    # Add try - except - raise statements\n",
    "    try:\n",
    "        for letter in reverse:\n",
    "            result = result + complements[letter]\n",
    "    except KeyError:\n",
    "        #print(\"All letters must be capitalized (A, C, G, or T)\", file=sys.stderr)\n",
    "        raise ValueError(\"Input needs to be a capital A, C, G, or T\")\n",
    "    return(result)\n",
    "\n",
    "print(reverse_complement(\"CAAg\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `Error`s automatically output to `sys.stderr`.\n",
    "\n",
    "Alternatively, we can check the input for the dictionary and produce an error ourselves. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit the function below!\n",
    "def reverse_complement(dna_sequence):\n",
    "    \"\"\"Reverses the complement of a dna sequence\"\"\"\n",
    "    complements = {\"T\":\"A\", \"A\":\"T\", \"C\":\"G\", \"G\":\"C\"}\n",
    "    reverse = dna_sequence[::-1]\n",
    "    result = \"\"\n",
    "    for letter in reverse:\n",
    "        # Check that letter is valid, if not raise an Error.\n",
    "        result = result + complements[letter]\n",
    "    return(result)\n",
    "\n",
    "print(reverse_complement(\"CAAg\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer to the above exercise:\n",
    "def reverse_complement(dna_sequence):\n",
    "    \"\"\"Reverses the complement of a dna sequence\"\"\"\n",
    "    complements = {\"T\":\"A\", \"A\":\"T\", \"C\":\"G\", \"G\":\"C\"}\n",
    "    reverse = dna_sequence[::-1]\n",
    "    result = \"\"\n",
    "    for letter in reverse:\n",
    "        # Check that letter is valid, if not raise an Error.\n",
    "        if letter not in complements:\n",
    "            raise ValueError(f\"Input needs to be a capital A, C, G, or T - invalid input was \\\"{letter}\\\"\")\n",
    "        result = result + complements[letter]\n",
    "    return(result)\n",
    "\n",
    "print(reverse_complement(\"CAAg\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some notes on `if-else` versus `try-catch` for errors\n",
    "\n",
    "We just demonstrated two different ways to handle potential errors: `if-else` blocks and `try-catch` blocks. Both of them protected against invalid input, so what's the difference?\n",
    "* In the `if-else` block, we checked the value of `letter` *before* trying to use it in the `complements` dictionary. In the `try-catch` block, we immediately called `complements[letter]`, which threw an error.\n",
    "* If you already know what your inputs should or shouldn't contain, it's *probably* a good idea to code in that information beforehand througn an `if-else`.\n",
    "* `try-catch` immediately runs code and stops when an error happens. This means the things done in `try` aren't reversed when you hit an error.\n",
    "* If you don't know *how exactly* something might break, but know it *can* happen, use a `try-catch`. For example, throw a `KeyError` or `ValueError` for keys that don't exist in your dictionary, but you can't guess what keys the user might try to use.\n",
    "* For some more discussion on this, see this StackOverflow forum post: https://stackoverflow.com/questions/7604636/better-to-try-something-and-catch-the-exception-or-test-if-its-possible-first."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application: Sanitizing Input\n",
    "\n",
    "Another problem is when programs produce incorrect results instead of producing an error. Suppose we have a function that prints all k-mers of a given *k* from a sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Function to return all k-mers from a sequence\n",
    "def kmers_from_sequence(dna_sequence, k):\n",
    "    positions = len(dna_sequence) - k + 1\n",
    "    for i in range(positions):\n",
    "        kmer = dna_sequence[i:i + k]\n",
    "        print(kmer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmers_from_sequence(\"CACGTGACTAG\", 3)\n",
    "print(\"After the function\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmers_from_sequence(\"CACGTGACTAG\", -3)\n",
    "print(\"After the function\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can **sanitize** the inputs to solve this. The value, *k*, should be a number greater than 0 and not longer than the length of the sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Refactor the following function to check that the value of `k` is both:\n",
    "- A positive number\n",
    "- Not longer than the length of `dna_sequence`\n",
    "\n",
    "If there is a problem, `raise` a `ValueError` with an appropriate message. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit the function below\n",
    "def kmers_from_sequence(dna_sequence, k):\n",
    "    # Write code to check input here!\n",
    "    \n",
    "    positions = len(dna_sequence) - k + 1\n",
    "    for i in range(positions):\n",
    "        kmer = dna_sequence[i:i + k]\n",
    "        print(kmer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmers_from_sequence(\"CAATCGACGTA\", 12)  # Should return an error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer to the above exercise\n",
    "def kmers_from_sequence(dna_sequence, k):\n",
    "    # Write code to check input here!\n",
    "    if k < 1 or k > len(dna_sequence):\n",
    "        raise ValueError(\"k needs to be at least 1 and no longer than the length of dna_sequence\")\n",
    "    positions = len(dna_sequence) - k + 1\n",
    "    for i in range(positions):\n",
    "        kmer = dna_sequence[i:i + k]\n",
    "        print(kmer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Syntactical shortcut: Separate code with line breaks\n",
    "\n",
    "If you have lines that are long and hard to read, putting in line breaks can help. In Python, you can have line breaks inside parentheses. Let's demonstrate this on some code we wrote yesterday:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import data from yesterday\n",
    "import pandas as pd\n",
    "gapminder = pd.read_table(\"gapminderDataFiveYear_superDirty.txt\", sep = \"\\t\")\n",
    "gapminder['region'] = gapminder['region'].astype(str)\n",
    "\n",
    "# Method 1 for formatting the 'region' column:\n",
    "gapminder['region'] = gapminder['region'].str.lstrip() # Strip white space on left\n",
    "gapminder['region'] = gapminder['region'].str.rstrip() # Strip white space on right\n",
    "gapminder['region'] = gapminder['region'].str.lower() # Convert to lowercase\n",
    "\n",
    "# Method 2 for formatting the 'region' column:\n",
    "gapminder['region'] = gapminder['region'].str.lstrip().str.rstrip().str.lower() # Strip white space on left and right, and convert to lowercase\n",
    "\n",
    "print(gapminder['region'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three different transformations happening above: removing whitespace on the left, removing whitespace on the right, and converting the text to lowercase. We can make this one line more intuitive by breaking it up into three:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# New method of chaining functions\n",
    "gapminder_copy['region'] = (\n",
    "    gapminder['region']\n",
    "    .str.lstrip()\n",
    "    .str.rstrip()\n",
    "    .str.lower()\n",
    ")\n",
    "\n",
    "print(gapminder_copy['region'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same output as above! This code is functionally the same as methods 1 and 2. We benefit from explicitly delineating each step like in method 1, and we also get the nicer syntax of applying all cleaning steps at the same time with method 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outlining your code - Pseudocode\n",
    "\n",
    "Thinking about what you want your future code to do *before* coding anything reduces the time you spend physically coding. It forces you to think about the big pieces that go into solving your problem and how they'll fit together, revealing potential problems much earlier. Let's take an example from day 1 to illustrate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percent = 20\n",
    "if percent < 39:\n",
    "    print('Low')\n",
    "elif percent < 47:\n",
    "    print('Normal')\n",
    "else:\n",
    "    print('High')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make this code more relevant for us, let's introduce a real biological variable: *hematocrit*. Hematocrit is the percentage (by volume) of red blood cells in blood. The normal values for humans are:\n",
    "* Male: 41% - 50%\n",
    "* Female: 36% - 44%\n",
    "* Average: 39% - 47%  -  these percentages used in the code above\n",
    "\n",
    "Now let's go over the logic of our previous code:\n",
    "1. We first check if the percent is less than 39: if so, then label as \"Low\".\n",
    "2. We then check if the percent is less than 47: if so, then label as \"Normal\".\n",
    "3. Otherwise, label as \"High\".\n",
    "\n",
    "This code seems straightforward at first glance, but there's a conceptual oversight: the second case here isn't strictly correct. Values less than 47% are normal *only if* the value isn't already less than 39%. The logic presented in this code works by virtue of checking the 39% case before the 47% case. Suppose we unintentionally coded in case 2 before case 1: now we first check for values less than 47, and any values that fulfill that condition are \"Normal\", even if they're *also* less than 39. This is an easy and common accident that can lead to erroneous results.\n",
    "\n",
    "It's more accurate to explicitly state the values that fall under each category:\n",
    "* Values between 0% - 38% are \"Low\"\n",
    "* Values between 39% - 47% are \"Normal\"\n",
    "* Values between 48% - 100% are \"High\"\n",
    "\n",
    "How would we rewrite the previous logic to follow this biological meaning?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jot down some new logic here!\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "tags": []
   },
   "source": [
    "Answer to above prompt:\n",
    "1. Check if the percent is between 0 and 38: if so, then label as \"Low\".\n",
    "2. Then, check if the percent is between 39 and 47: if so, then label as \"Normal\".\n",
    "3. Then, check if the percent is between 48 and 100: if so, then label as \"High\".\n",
    "4. (new condition!) If it's any other integer, raise an error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is better! We've done a few defensive programming concepts here:\n",
    "* Written up **pseudo-code**: this is not actual code, but an outline of how you want your code to be structured.\n",
    "* Sanitized our input and guarded against a potential error.\n",
    "* More explicitly stated the biological meaning of our code, making it easier for others to follow.\n",
    "* Defended against the concept of a \"wrong order\" for our if/else statements - now it doesn't matter how the three conditions are ordered.\n",
    "\n",
    "Now, re-write the actual code to follow this new logic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# New code here!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer to above exercise:\n",
    "percent = 20\n",
    "if 0 <= percent <= 38:\n",
    "    print('Low')\n",
    "elif 39 <= percent <= 47:\n",
    "    print('Normal')\n",
    "elif 48 <= percent <= 100:\n",
    "    print('High')\n",
    "else:\n",
    "    raise ValueError(\"percent value must be between 0 and 100\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with files\n",
    "\n",
    "Real-world data will typically be in a file, **not** in your code. Python has functions (and libraries) to read files in various formats, including text, image, and binary formats. Sequence data is typically stored in plain text files, like FASTA. Let's walk through reading one of those."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open and close\n",
    "\n",
    "When you use a word processor or spreadsheet, you open files, work with them, and then close them when you're done. In Python, you do the same thing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Opening and closing a file\n",
    "f = open('ae.fa')\n",
    "for line in f:\n",
    "    print(line)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go through the steps we just did:\n",
    "1. We used the `open()` function on a string that represents a path to a file.\n",
    "    * The result of that function was saved to the variable `f`. This value is called a *file object*.\n",
    "2. We wrote a `for` loop. When you write a `for` loop for a file object, each time a loop happens, one line in the file is read and stored in the loop variable `line`.\n",
    "3. We printed the loop variable `line` for each loop.\n",
    "4. We used the `close()` function to close the file.\n",
    "\n",
    "Text files are generally read line-by-line in a loop like we just did with this FASTA file. Each line of text from the file is set into our loop variable, at which point we're free to do whatever we want with the loop variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading lines\n",
    "\n",
    "Notice the blank line in between each line from the file. Text files and programs indicate that there is a new line using a special newline character, `\\n`. In our previous example, each line in the file includes the `\\n` newline character at the end.\n",
    "\n",
    "Let's visualize this by adding each line to a single string, `all_lines`, and compare printing the string versus the raw data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Visualizing newline characters\n",
    "f = open('ae.fa')\n",
    "all_lines = ''  # Initializes a new variable to store all the strings\n",
    "for line in f:\n",
    "    print(line)\n",
    "    all_lines = all_lines + line  # Adds each successive line to the variable\n",
    "f.close()\n",
    "all_lines  # Prints the variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If desired, we can strip off the newline character with the `.strip()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Use strip() to remove newlines\n",
    "string_with_newline = \"There is a newline at the end\\n\"\n",
    "string_with_newline\n",
    "string_with_newline = string_with_newline.strip()\n",
    "string_with_newline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try this with our original code to print each line in a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "# Modifying our original code to remove newlines\n",
    "f = open('ae.fa')\n",
    "for line in f:\n",
    "    line = line.strip()\n",
    "    print(line)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, when reading text files, check that you're not accidentally including extra `\\n` characters (or excluding them when you need to separate lines)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Processing a fasta file\n",
    "\n",
    "Let's put some concepts together: write a function called `read_fasta` that takes a filename and returns a `list` of all the DNA sequences in the file.\n",
    "* Each element in the list should correspond to one FASTA entry.\n",
    "* Make sure to remove the header rows and only include the DNA sequences.\n",
    "* Test your code with the `ls_orchid.fasta` file.\n",
    "\n",
    "Some hints:\n",
    "* Look at `ls_orchid.fasta` before coding to know how your input data is formatted.\n",
    "* There are multiple FASTA sequences in this file - you'll need to know when to stop reading the previous sequence and start reading the next one.\n",
    "* Remember to strip the newlines off each FASTA sequence.\n",
    "* It's very helpful to write pseudocode first, to go through the logic of your function before writing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Answer to above exercise:\n",
    "def read_fasta(filename):\n",
    "    sequences = []\n",
    "    temp = ''\n",
    "    f = open(filename)\n",
    "    for line in f:\n",
    "        if line[0] in ['A','T','C','G','N']:\n",
    "            temp = temp + line.strip()\n",
    "        elif line[0] == '\\n':\n",
    "            sequences.append(temp)\n",
    "            temp = ''\n",
    "    f.close()\n",
    "    return sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(read_fasta('ls_orchid.fasta'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scripts\n",
    "\n",
    "We've done a good job of organizing our code into functions, but we've only been running them from this notebook. Suppose you want to run this program outside of Jupyter Notebook (e.g. on the command line) or give it to someone else to run: we want to nicely package this script in a way that makes it easy for the next person to run. We wouldn't just want to hand them this entire Jupyter Notebook and say \"the right function's in there somewhere, just keep looking\" or \"yeah, you just need to edit the arguments to use your own files\"!\n",
    "\n",
    "Let's take our `read_fasta` function and put it in a script that reads the `ls_orchid.fasta` file and prints its contents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo - copy/paste from above and edit\n",
    "%%writefile cool_functions.py\n",
    "def read_fasta(filename):\n",
    "    sequences = []\n",
    "    temp = ''\n",
    "    f = open(filename)\n",
    "    for line in f:\n",
    "        if line[0] in ['A','T','C','G','N']:\n",
    "            temp = temp + line.strip()\n",
    "        elif line[0] == '\\n':\n",
    "            sequences.append(temp)\n",
    "            temp = ''\n",
    "    f.close()\n",
    "    return sequences\n",
    "\n",
    "print(read_fasta('ls_orchid.fasta'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the first line contains a `%%` operator followed by the command `writefile` and a file name. This operator is specific to Jupyter notebooks, called a \"Cell Magic Command\", and copies the code written in a cell into a file.\n",
    "\n",
    "Our script reads the `ls_orchid.fa` file (and only the `ls_orchid.fa` file) every time it's run, an example of **hard-coding** an argument. Most programs let you specify the arguments you want to use, as opposed to being inflexible and pre-determined every time you run it. Let's change our script again to achieve this goal.\n",
    "\n",
    "In Python, our program can get these arguments, but we have to load the `sys` module. Remember from earlier when we imported `sys` to call `sys.stderr`? Well, `sys` can also grab arguments from the command line. Let's change our `read_fasta` function slightly to allow it to read in *any* file by using the `sys` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit function above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Live demo\n",
    "%%writefile cool_functions.py\n",
    "import sys\n",
    "\n",
    "def read_fasta(filename):\n",
    "    sequences = []\n",
    "    temp = ''\n",
    "    f = open(filename)\n",
    "    for line in f:\n",
    "        if line[0] in ['A','T','C','G','N']:\n",
    "            temp = temp + line.strip()\n",
    "        elif line[0] == '\\n':\n",
    "            sequences.append(temp)\n",
    "            temp = ''\n",
    "    f.close()\n",
    "    return sequences\n",
    "\n",
    "print(read_fasta(sys.argv[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what happens if we don't have an input file name? According to the `sys` documentation, `sys.argv` returns a list where the first item `sys.argv[0]` is the name of the script by default, and any additional items in the list are the command line arguments. If no arguments were passed, `sys.argv` should be a list with exactly one element, the script name. Let's add some more code to our script to check whether the user passed in a filename:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit function above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo\n",
    "%%writefile cool_functions.py\n",
    "import sys\n",
    "\n",
    "def read_fasta(filename):\n",
    "    sequences = []\n",
    "    temp = ''\n",
    "    f = open(filename)\n",
    "    for line in f:\n",
    "        if line[0] in ['A','T','C','G','N']:\n",
    "            temp = temp + line.strip()\n",
    "        elif line[0] == '\\n':\n",
    "            sequences.append(temp)\n",
    "            temp = ''\n",
    "    f.close()\n",
    "    return sequences\n",
    "\n",
    "if len(sys.argv) < 2:\n",
    "    print('Usage:', sys.argv[0], '<sequence.fa>')\n",
    "    sys.exit(1)\n",
    "\n",
    "print(read_fasta(sys.argv[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use yet another function of `sys`, `sys.exit`, which lets us stop a script or program before the whole script is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing your own scripts\n",
    "\n",
    "So far, we have used modules to help us work on our analyses such as:\n",
    "- Standard Library\n",
    "    - `sys`\n",
    "- Third Party\n",
    "    - `pandas`\n",
    "    - `numpy`\n",
    "    - `matplotlib`\n",
    "\n",
    "These are imported using the `import` keyword, which lets us use their functions. We also write functions for use in our own code. Having our own functions available to import into other scripts gives the benefit of:\n",
    "1. Letting us reuse code over multiple analyses (DRY).\n",
    "2. Letting others use our code in their own scripts without copy/pasting (DRY).\n",
    "\n",
    "While it may seem like going out of one's way to write a module for analysis, you can have the same Python file work in two different ways:\n",
    "1. A module that's imported into another Python notebook or script\n",
    "2. A standalone script that's run from the command line\n",
    "\n",
    "Let's demonstrate this by making a new Jupyter notebook file: `New -> Python3 Notebook`. In this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cool_functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is this giving an error when we import it into our Jupyter notebook? This is because we're using `sys.argv` to read arguments from the command line, but we're *not* on the command line right now, so there are no arguments/files to read in. We're attempting to do method #1 (import as a module within another Python file) with code that's meant for method #2 (run as a standalone script on the command line).\n",
    "\n",
    "So, how do we separate the `print` statement containing `sys.argv` such that it's only run when the script is run from the command line?\n",
    "\n",
    "Answer: Add this `if` statement before the `print` statement: `if __name__ == \"__main__\":`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit function above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Live demo\n",
    "%%writefile cool_functions.py\n",
    "import sys\n",
    "\n",
    "def read_fasta(filename):\n",
    "    sequences = []\n",
    "    temp = ''\n",
    "    f = open(filename)\n",
    "    for line in f:\n",
    "        if line[0] in ['A','T','C','G','N']:\n",
    "            temp = temp + line.strip()\n",
    "        elif line[0] == '\\n':\n",
    "            sequences.append(temp)\n",
    "            temp = ''\n",
    "    f.close()\n",
    "    return sequences\n",
    "\n",
    "if len(sys.argv) < 2:\n",
    "    print('Usage:', sys.argv[0], '<sequence.fa>')\n",
    "    sys.exit(1)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    print(read_fasta(sys.argv[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go back to your new Jupyter notebook file and re-run the `import` cell - it should no longer throw an error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Adding new functions to your module\n",
    "\n",
    "In your new notebook, you've now imported the `cool_functions.py` file, which contains the `read_fasta` function that returns a list of DNA sequences from a fasta file. Building off this module, how would you get a count of unique dinucleotides (AA, AT, AC, AG, TA, etc.) for each FASTA sequence? Write out some pseudocode first, then implement it.\n",
    "* Hint: you can loop over a string using `for letter in string`, giving you the actual letter.\n",
    "* Another hint: you can loop over the positions in a string using `for int in range(len(string))`, giving you integers for each letter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put your answer here!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# One potential answer - exercise is intentionally left somewhat ambiguous\n",
    "def count_dinucleotides(seq):\n",
    "    \"\"\"Given a DNA sequence, prints the counts of unique dinucleotides in the sequence.\n",
    "    \"\"\"\n",
    "    dnts = {}\n",
    "    for i in range(len(seq)-1):\n",
    "        dinucleotide = seq[i] + seq[i+1]\n",
    "        if dinucleotide not in dnts:\n",
    "            dnts[dinucleotide] = 1\n",
    "        else:\n",
    "            dnts[dinucleotide] += 1\n",
    "    print(dnts)\n",
    "    \n",
    "def count_all_dinucleotides(file):\n",
    "    \"\"\"Given a FASTA file, prints the count of unique dinucleotides for all sequences.\n",
    "    \"\"\"\n",
    "    for seq in read_fasta(file):\n",
    "        count_dinucleotides(seq)\n",
    "\n",
    "count_all_dinucleotides(\"ls_orchid.fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automation: Combining Python and Bash"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Automation saves both time and effort in large quantities. You save your own time by having the computer automatically loop over data files, performing multiple analyses in parallel, and generating output files and visualizations. By saving these directions in the form of code, you also save time by making your work documented and reproducible, saving you the effort of re-doing them manually should the need arise.\n",
    "\n",
    "Suppose you're given a few *hundred* fasta files you need to concatenate. You could manually type them all into a list in your Jupyter notebook and run it, or you could have `Bash` automate the `Python` script for you!\n",
    "\n",
    "First, a technical check: in your Git Bash or Terminal window, run `python --version`. Hopefully, you get some info about the version of Python you're running.\n",
    "\n",
    "Recall our bash lesson on the first day:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for filename in *.fa*\n",
    "do\n",
    "    echo $filename\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `python` command in the terminal, we can also run Python files without needing to open up Jupyter notebook! You use this `python` command just like any other bash command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "python cool_functions.py ls_orchid.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "And we can use it in a for loop in Bash, just like any other command:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for filename in *.fa*\n",
    "do\n",
    "    python cool_functions.py $filename\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get a lot of text outputted to your terminal. Better yet, let's redirect it to another file to save for later:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for filename in *.fa*\n",
    "do\n",
    "    python cool_functions.py $filename >> output_fastas.txt\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have seen the `>` bash symbol before, which means to redirect the output of a command into another file. Here, we use `>>`: it's is similar in that it also redirects output to a file, but this *concatenates* the result to whatever is in the file, instead of overwriting the whole file with the new results.\n",
    "\n",
    "Now we can open up and check `output_fastas.txt`. There's one line for every file processed.\n",
    "\n",
    "## Automating the Bash command\n",
    "\n",
    "Okay, time for one more level of automation: it's great that we can get Bash to loop over our Python script and our input files, but what if we don't want to type out the `for` loop in bash every single time we run an analysis? We can store that command in a Bash script, similar to how we stored our Python code in a Python script.\n",
    "\n",
    "In either `nano` or `vim`, copy and paste the above code into a new file called `script.sh`. Remember that we can create a new file by directly typing `nano script.sh` or `vim script.sh` on the command line. Once you've created your `script.sh` file, run it on the command line with:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "./script.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should work as intended (i.e. not output anything to the terminal), and you can open up your `output_fastas.txt` file to see the results.\n",
    "\n",
    "## Exercise: Optimizing your script\n",
    "\n",
    "How would you modify this `script.sh` so that it empties the contents of `output_fastas.txt` before running your program? Hint: there are multiple bash commands you could potentially use."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Edit your bash script!"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Solutions - any one of these will work:\n",
    "rm output_fastas.txt\n",
    "touch output_fastas.txt\n",
    "echo \"\" > output_fastas.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automating your time: Walking away from your computer\n",
    "\n",
    "Congratulations, you now have a `Bash` script that automates a `Python` script over however many fasta files you have! What if you were dealing with gigabytes or even terabytes of FASTA files? You'd be waiting forever for your script to finish, but you probably don't want to sit in front of your computer that long. What can you do?"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "nohup ./script.sh &"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* `nohup`: Stands for `no hangup` - even if you close your bash terminal, the script will continue to run in the background of your computer. Just make sure you don't shut down your computer before it finishes! On a computing cluster, this isn't really a problem since compute clusters generally stay online 24/7.\n",
    "  * The Duke Compute Cluster (DCC) is an example of a computing cluster: https://oit-rc.pages.oit.duke.edu/rcsupportdocs/dcc/. It provides free computational time and storage (specific amounts will vary) to Duke affiliates. Ask your PI for access if you're interested!\n",
    "* `&`: Run this program in the background of your terminal. This frees up your terminal so you can work on other things and run more commands. Note that this *does not* keep it running if the terminal is closed; you'll still need `nohup` for that.\n",
    "\n",
    "TL,DR: `&` *frees up* your terminal so you can do other things, `nohup` *lets you close* your terminal while still running your script."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}