orf_finder.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence = \"ATCCATGACTTGACTATGACGTTCAGCATCTGAAGACTTATGAAAATGATCTCATGGATCGACT\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`               ↑↑↑---↑↑↑  ↑↑↑------------↑↑↑      ↑↑↑----↑↑↑    ↑↑↑--------`\n",
    "\n",
    "If you double click on this in Jupyter, you can see the arrows line up with where we inserted ORFs for testing purposes.\n",
    "\n",
    "Please note, that the third \"ORF\" is not a **valid** ORF (hence the quotation marks), since there's 4 nucleotides between the start and the end codon, so it should NOT be picked up by our algorithm (stop codon isn't in the same reading frame!\n",
    "\n",
    "Also, the last start codon does **not** have a matching stop codon, so it shouldn't be found, either.\n",
    "\n",
    "Of course there are more things to test and consider -- for example an ORF inside another ORF. We encourage you to think about what might happen (step through the algorithm by hand to really try to understand what's happening) and finally verify your assumptions with an experiment.\n",
    "\n",
    "Finally, it's very well possible that there's still bugs in here we haven't caught. Please feel free to point these out to us :)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_atg(sequence, position):\n",
    "    \"\"\"Return the position of the first encountered start\n",
    "    codon in ``sequence``, starting from ``position``\"\"\"\n",
    "    for start_pos in range(position, len(sequence)):\n",
    "        if sequence[start_pos:(start_pos + 3)] == \"ATG\":\n",
    "            return start_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that there's definitely things we could do to still improve this approach! Here's some potentially interesting questions:\n",
    "\n",
    "- Do we need to traverse the *entire* string, or might it be sensible to stop a couple of letters before the end?\n",
    "- If we're taking the range up to `len(sequence)` but are then accessing `sequence[start_pos:(start_pos + 3)]`, why is it okay to even go beyond what should be the maximum index of our sequence?\n",
    "- What if we didn't *just* want to look for \"ATG\" as a start codon? Could you think of ways to make this more flexible?\n",
    "\n",
    "Also, if you've noticed the usage of the message in triple quotation marks (`\"\"\"Hello!\"\"\"`) right below the function definition, those are called \"docstrings\" and I recommend reading up on them :)! They're there for documentation purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_stop(sequence, position):\n",
    "    \"\"\"Return an ORF's sequence by finding the first stop codon\n",
    "    to appear after a start codon found at ``position``\"\"\"\n",
    "    for stop_pos in range(position, len(sequence), 3):\n",
    "        if sequence[stop_pos:(stop_pos + 3)] in (\"TAG\", \"TGA\", \"TAA\"):\n",
    "            return sequence[position:(stop_pos + 3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like with the previous function definition, think about how you might implement this differently, or even how to improve on this solution.\n",
    "\n",
    "Points to note:\n",
    "\n",
    "- The usage of the `in` keyword, which is a *very* handy tool to have in your belt. If you've not encountered it before, I'd strongly recommend reading up on it!\n",
    "- Also note that we're passing an unassuming `3` to the `range` function as well! What does this accomplish? (If you need a hint: think about something mentioned in the \"test case\" string at the top!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Begin with initialising your variables\n",
    "orfs = []\n",
    "position = 0\n",
    "\n",
    "# Loop through the sequence, i.e. step through every character\n",
    "while position < len(sequence):\n",
    "    \n",
    "    # Find first (if any) start codon\n",
    "    start = find_atg(sequence, position)\n",
    "    \n",
    "    # If there's no start codon in the entire sequence, just exit the loop\n",
    "    if not start:\n",
    "        break\n",
    "    \n",
    "    # Try to find a sequence by looking for stop codons after\n",
    "    # the start codon we just determined\n",
    "    orf = find_stop(sequence, start)\n",
    "    \n",
    "    # If we found one, add it to our list of codons, including where it started\n",
    "    if orf:\n",
    "        orfs.append((orf, start))\n",
    "        \n",
    "    # Finally, step to the next character in the sequence and repeat the process\n",
    "    position += start + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step through this by hand, if you're uncertain about what's going on. Whenever there's a function called, go to that definition and apply it, step by step, too.\n",
    "\n",
    "Things to note:\n",
    "\n",
    "- Variable initialisation\n",
    "- The use of a `while` loop – could/should we have used something different? What else needs to be taken care of? (Hint: why are we performing `position += start + 1` at the end?\n",
    "- `break` – another keyword to check up on (check out https://docs.python.org/3/library/keyword.html!)\n",
    "\n",
    "And as always: how could we do better/differently? Do you see any issues?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('ATGACTTGA', 4), ('ATGACGTTCAGCATCTGA', 15)]\n"
     ]
    }
   ],
   "source": [
    "# Let's check our results thus far:\n",
    "print(orfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_longest_orf(orfs):\n",
    "    \"\"\"Return the longest ORF in our list of ORFs\"\"\"\n",
    "    \n",
    "    # Initialise variable with the first ORF we found\n",
    "    longest = orfs[0]\n",
    "    \n",
    "    # Step through each remaining element in the list and see\n",
    "    # if any of these ORFs is longer than the first one\n",
    "    for orf in orfs[1:]:\n",
    "        if len(orf[0]) > len(longest[0]):\n",
    "            longest = orf\n",
    "    return longest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's an issue in this! Can you spot it? (Hint: what do you think happens, if the sequence we're checking doesn't have any ORFs in it? If you can't see it, try it out and figure out what's going wrong!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('ATGACGTTCAGCATCTGA', 15)\n"
     ]
    }
   ],
   "source": [
    "# Print our result:\n",
    "print(find_longest_orf(orfs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And done :)!"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}