ORF-Open_reading_frames.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open reading frames, from [Rosalind.info](https://www.rosalind.info)\n",
    "\n",
    "(Text copied from http://rosalind.info/problems/orf/)\n",
    "\n",
    "<div class=\"problem-statement problem-statement-bordered\" problem=\"122\">\n",
    "<h2 id=\"problem\">Problem</h2>\n",
    "<p>Either strand of a DNA double helix can serve as the <a class=\"term\" href=\"/glossary/coding-strand/\" id=\"term-322\" rel=\"tooltip\" title=\"\n",
    "The strand of a double-stranded DNA molecule that is copied or transcribed into RNA.\">coding strand</a> for RNA transcription.\n",
    "Hence, a given DNA string implies six total <a class=\"term new\" href=\"/glossary/reading-frame/\" id=\"term-267\" rel=\"tooltip\" title=\"New term: \n",
    "One of three possible ways to read a given strand of DNA, depending upon the starting position.\">reading frames</a>, or ways in which the same region of DNA can be translated into amino acids:\n",
    "three reading frames result from reading the string itself, whereas three more result\n",
    "from reading its <a class=\"term\" href=\"/glossary/reverse-complement/\" id=\"term-252\" rel=\"tooltip\" title=\"\n",
    "The DNA string formed by reversing and complementing each symbol.\">reverse complement</a>.</p>\n",
    "<p>An <a class=\"term new\" href=\"/glossary/open-reading-frame/\" id=\"term-268\" rel=\"tooltip\" title=\"New term: \n",
    "A sequence in DNA or RNA potentially able to encode the protein.\">open reading frame</a> (ORF) is one which starts from the <a class=\"term\" data-math=\"true\" href=\"/glossary/start-codon/\" id=\"term-260\" rel=\"tooltip\" title=\"\n",
    "The RNA codon $\\textrm{AUG}$, which codes for the amino acid methionine and indicates the beginning\n",
    "of translation into protein.\">start codon</a> and ends by <a class=\"term\" href=\"/glossary/stop-codon/\" id=\"term-261\" rel=\"tooltip\" title=\"\n",
    "One of three possible RNA codons that indicate the termination of protein translation.\">stop codon</a>, without\n",
    "any other <a class=\"term\" href=\"/glossary/stop-codon/\" id=\"term-261\" rel=\"tooltip\" title=\"\n",
    "One of three possible RNA codons that indicate the termination of protein translation.\">stop codons</a> in between. Thus, a candidate protein string is derived by translating an open reading\n",
    "frame into amino acids until a stop codon is reached.</p>\n",
    "<p><span class=\"given-return\">Given:</span> A <a class=\"term\" href=\"/glossary/dna-string/\" id=\"term-349\" rel=\"tooltip\" title=\"\n",
    "A string constructed from the alphabet {A, C, G, T}.\">DNA string</a> <mathjax>$s$</mathjax> of length at most 1 <a class=\"term\" href=\"/glossary/kbp/\" id=\"term-394\" rel=\"tooltip\" title=\"\n",
    "1 kbp = 1000 base pairs\">kbp</a> in <a class=\"term\" href=\"/glossary/fasta-format/\" id=\"term-759\" rel=\"tooltip\" title=\"\n",
    "A text format used for naming genetic strings in databases.\">FASTA format</a>.</p>\n",
    "<p><span class=\"given-return\">Return:</span> Every distinct candidate protein string that can be translated from ORFs of <mathjax>$s$</mathjax>.\n",
    "Strings can be returned in any order.</p>\n",
    "<h2 id=\"sample-dataset\">Sample Dataset</h2>\n",
    "<div class=\"codehilite\"><pre>&gt;Rosalind_99\n",
    "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG\n",
    "</pre></div>\n",
    "\n",
    "\n",
    "<h2 id=\"sample-output\">Sample Output</h2>\n",
    "<div class=\"codehilite\"><pre>MLLGSFRLIPKETLIQVAGSSPCNLS\n",
    "M\n",
    "MGMTPRLGLESLLE\n",
    "MTPRLGLESLLE\n",
    "</pre></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## My interpretation/reasoning\n",
    "\n",
    "1. This problem seems straightforward: return all possible six frame translations for a given DNA sequence\n",
    "\n",
    "2. Several tools exist that predict and translate open reading frames, such as [this useful webtool by ExPASy](https://web.expasy.org/translate/), or [Prodigal](https://github.com/hyattpd/Prodigal) also seems to use something like six frame translations. \n",
    "\n",
    "3. However, for this particular exercise, these tools do not seem so handy: the  output needs to be formatted as above: One translation per line, each possible translation is requested, and nothing else.\n",
    "\n",
    "So...\n",
    "\n",
    "To put this into a script, I want to have:\n",
    " - A function to read the sequence from a fasta file\n",
    " - Translations for the DNA sequences in all 6 possible open reading frames, _not_ stopping at any one particular stop codon. That is, stop the current translation, but keep looking in the rest of the sequence.\n",
    " - This bit of Biopython's documentation can help me do this: https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc295\n",
    " \n",
    "In particular, this piece of code may be of use:\n",
    "\n",
    "```python\n",
    ">>> from Bio import SeqIO\n",
    ">>> record = SeqIO.read(\"NC_005816.fna\", \"fasta\")\n",
    ">>> table = 11\n",
    ">>> min_pro_len = 100\n",
    "\n",
    "Here is a neat trick using the Seq object’s split method to get a list of all the possible ORF translations in the six reading frames:\n",
    "\n",
    ">>> for strand, nuc in [(+1, record.seq), (-1, record.seq.reverse_complement())]:\n",
    "...     for frame in range(3):\n",
    "...         length = 3 * ((len(record)-frame) // 3) #Multiple of three\n",
    "...         for pro in nuc[frame:frame+length].translate(table).split(\"*\"):\n",
    "...             if len(pro) >= min_pro_len:\n",
    "...                 print(\"%s...%s - length %i, strand %i, frame %i\" \\\n",
    "...                       % (pro[:30], pro[-3:], len(pro), strand, frame))\n",
    "GCLMKKSSIVATIITILSGSANAASSQLIP...YRF - length 315, strand 1, frame 0\n",
    "KSGELRQTPPASSTLHLRLILQRSGVMMEL...NPE - length 285, strand 1, frame 1\n",
    "GLNCSFFSICNWKFIDYINRLFQIIYLCKN...YYH - length 176, strand 1, frame 1\n",
    "VKKILYIKALFLCTVIKLRRFIFSVNNMKF...DLP - length 165, strand 1, frame 1\n",
    "NQIQGVICSPDSGEFMVTFETVMEIKILHK...GVA - length 355, strand 1, frame 2\n",
    "RRKEHVSKKRRPQKRPRRRRFFHRLRPPDE...PTR - length 128, strand 1, frame 2\n",
    "TGKQNSCQMSAIWQLRQNTATKTRQNRARI...AIK - length 100, strand 1, frame 2\n",
    "QGSGYAFPHASILSGIAMSHFYFLVLHAVK...CSD - length 114, strand -1, frame 0\n",
    "IYSTSEHTGEQVMRTLDEVIASRSPESQTR...FHV - length 111, strand -1, frame 0\n",
    "WGKLQVIGLSMWMVLFSQRFDDWLNEQEDA...ESK - length 125, strand -1, frame 1\n",
    "RGIFMSDTMVVNGSGGVPAFLFSGSTLSSY...LLK - length 361, strand -1, frame 1\n",
    "WDVKTVTGVLHHPFHLTFSLCPEGATQSGR...VKR - length 111, strand -1, frame 1\n",
    "LSHTVTDFTDQMAQVGLCQCVNVFLDEVTG...KAA - length 107, strand -1, frame 2\n",
    "RALTGLSAPGIRSQTSCDRLRELRYVPVSL...PLQ - length 119, strand -1, frame 2\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_single_sequence_file(input_file):\n",
    "    \"\"\"\n",
    "Read a single sequence from a file.\n",
    "    \"\"\"\n",
    "    return(SeqIO.read(input_file, \"fasta\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def translate_six_frames(seq_record):\n",
    "    \"\"\"\n",
    "Given a Biopython sequence record with a DNA (or RNA?) sequence,\n",
    "translate into amino acid (protein) sequences in six frames.\n",
    "Returns translations as list of strings.\n",
    "    \"\"\"\n",
    "    translation_list = []\n",
    "    \n",
    "    for strand, nuc in [(+1, seq_record.seq), (-1,  seq_record.seq.reverse_complement())]:\n",
    "        for frame in range(3):\n",
    "            length = 3 * ((len(seq_record)-frame) // 3)\n",
    "            for pro in nuc[frame:frame+length].translate().split(\"*\"):\n",
    "                translation_list.append(pro)\n",
    "                \n",
    "    return(translation_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_record = read_single_sequence_file(\"data/Example_open_reading_frames.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_list = translate_six_frames(test_record)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SHVANSGYMGMTPRLGLESLLE\n",
      "A\n",
      "MIRVASQ\n",
      "AM\n",
      "LTQVTWG\n",
      "PRDLD\n",
      "SLFWNKPE\n",
      "SE\n",
      "HL\n",
      "PCS\n",
      "LRLHGDDPATWIRVSFGISLNDPSSIS\n",
      "LRCYSDHSGLFQKRL\n",
      "SKSRGHPHVT\n",
      "VSYMA\n",
      "\n",
      "DATRIIQAYSKRDSNPSRGVIPM\n",
      "PELATW\n",
      "EMLLGSFRLIPKETLIQVAGSSPCNLS\n",
      "LHG\n"
     ]
    }
   ],
   "source": [
    "for translation in test_list: print(translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not quite what I expected... The list is way too long and does not always start with start codons (or Ms in the translation).\n",
    "Let's debug the function a bit by using print statements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_translate_six_frames(seq_record):\n",
    "    \"\"\"\n",
    "Given a Biopython sequence record with a DNA (or RNA?) sequence,\n",
    "translate into amino acid (protein) sequences in six frames.\n",
    "Returns translations as list of strings.\n",
    "    \"\"\"\n",
    "    translation_list = []\n",
    "    \n",
    "    for strand, nuc in [(+1, seq_record.seq), (-1,  seq_record.seq.reverse_complement())]:\n",
    "        print(\"Strand: %s\\nNuc: %s\" % (strand, nuc))\n",
    "        for frame in range(3):\n",
    "            print(\"Frame: %s\" % frame)\n",
    "            length = 3 * ((len(seq_record)-frame) // 3)\n",
    "            print(\"Length: %s\" % length)\n",
    "            print(\"Possible translations: %s\" % nuc[frame:frame+length].translate())\n",
    "            for pro in nuc[frame:frame+length].translate().split(\"*\"):\n",
    "                translation_list.append(pro)\n",
    "                \n",
    "    return(translation_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Strand: 1\n",
      "Nuc: AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG\n",
      "Frame: 0\n",
      "Length: 96\n",
      "Possible translations: SHVANSGYMGMTPRLGLESLLE*A*MIRVASQ\n",
      "Frame: 1\n",
      "Length: 93\n",
      "Possible translations: AM*LTQVTWG*PRDLD*SLFWNKPE*SE*HL\n",
      "Frame: 2\n",
      "Length: 93\n",
      "Possible translations: PCS*LRLHGDDPATWIRVSFGISLNDPSSIS\n",
      "Strand: -1\n",
      "Nuc: CTGAGATGCTACTCGGATCATTCAGGCTTATTCCAAAAGAGACTCTAATCCAAGTCGCGGGGTCATCCCCATGTAACCTGAGTTAGCTACATGGCT\n",
      "Frame: 0\n",
      "Length: 96\n",
      "Possible translations: LRCYSDHSGLFQKRL*SKSRGHPHVT*VSYMA\n",
      "Frame: 1\n",
      "Length: 93\n",
      "Possible translations: *DATRIIQAYSKRDSNPSRGVIPM*PELATW\n",
      "Frame: 2\n",
      "Length: 93\n",
      "Possible translations: EMLLGSFRLIPKETLIQVAGSSPCNLS*LHG\n",
      "SHVANSGYMGMTPRLGLESLLE\n",
      "A\n",
      "MIRVASQ\n",
      "AM\n",
      "LTQVTWG\n",
      "PRDLD\n",
      "SLFWNKPE\n",
      "SE\n",
      "HL\n",
      "PCS\n",
      "LRLHGDDPATWIRVSFGISLNDPSSIS\n",
      "LRCYSDHSGLFQKRL\n",
      "SKSRGHPHVT\n",
      "VSYMA\n",
      "\n",
      "DATRIIQAYSKRDSNPSRGVIPM\n",
      "PELATW\n",
      "EMLLGSFRLIPKETLIQVAGSSPCNLS\n",
      "LHG\n"
     ]
    }
   ],
   "source": [
    "for translation in test_translate_six_frames(test_record): print(translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think I will need to do some searching, for example with Regular Expressions (`regex`), to filter out the open reading frames, which start with start codons and end with stop codons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import regex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('', 'TAG'), ('GAA', 'TAA'), ('GAA', 'TAA')]\n",
      "[('AGT', 'TAG'), ('', 'TAA')]\n"
     ]
    }
   ],
   "source": [
    "#Let's try a little test with the RegEx: ATG([ACGT]{3})*(TAA|TAG|TGA)\n",
    "# This should start with \"ATG\" (start), followed by any number of triplets consisting of A, C, G and/or T, and end with one of TAA, TAG or TGA (stop codons)\n",
    "\n",
    "test_forward = \"AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG\"\n",
    "\n",
    "test_reverse = \"CTGAGATGCTACTCGGATCATTCAGGCTTATTCCAAAAGAGACTCTAATCCAAGTCGCGGGGTCATCCCCATGTAACCTGAGTTAGCTACATGGCT\"\n",
    "\n",
    "pattern = \"ATG([ACGT]{3})*?(TAA|TAG|TGA)\"\n",
    "\n",
    "#Explanation:\n",
    "# 1. Start with \"ATG\" literally (start codon)\n",
    "# 2.1. Then match any number of triplets of A, C, G and/or G\n",
    "# 2.2. However, since '*' is a 'greedy' match, it will look for the longest possible match\n",
    "#     and I want to stop at the first stop codon (= point 3). Therefore, I added the '?'\n",
    "#     to make the matching method 'lazy'. (Also see https://regex101.com/)\n",
    "# 3. Stop the match at either \"TAA\", \"TAG\" or \"TGA\" (stop codons)\n",
    "\n",
    "print(regex.findall(r\"%s\" % pattern, test_forward, overlapped=True))\n",
    "\n",
    "print(regex.findall(r\"%s\" % pattern, test_reverse, overlapped=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m\n",
       "\u001b[0mregex\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfindall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mflags\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[0mpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mendpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0moverlapped\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mconcurrent\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mignore_unused\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\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[0mfindall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mendpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moverlapped\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m  \u001b[0mconcurrent\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mignore_unused\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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\"\"\"Return a list of all matches in the string. The matches may be overlapped\u001b[0m\n",
       "\u001b[0;34m    if overlapped is True. If one or more groups are present in the pattern,\u001b[0m\n",
       "\u001b[0;34m    return a list of groups; this will be a list of tuples if the pattern has\u001b[0m\n",
       "\u001b[0;34m    more than one group. Empty matches are included in the result.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0m_compile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpattern\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mignore_unused\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfindall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpos\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m      \u001b[0mendpos\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moverlapped\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconcurrent\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/miniconda3/envs/jupyter/lib/python3.8/site-packages/regex/regex.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "??regex.findall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<regex.Match object; span=(4, 88), match='ATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAG'>\n",
      "<regex.Match object; span=(24, 75), match='ATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGA'>\n",
      "<regex.Match object; span=(30, 75), match='ATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGA'>\n"
     ]
    }
   ],
   "source": [
    "for match in regex.finditer(r\"%s\" % pattern, test_forward, overlapped=True): print(match)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<regex.Match object; span=(5, 86), match='ATGCTACTCGGATCATTCAGGCTTATTCCAAAAGAGACTCTAATCCAAGTCGCGGGGTCATCCCCATGTAACCTGAGTTAG'>\n",
      "<regex.Match object; span=(70, 76), match='ATGTAA'>\n"
     ]
    }
   ],
   "source": [
    "for match in regex.finditer(r\"%s\" % pattern, test_reverse, overlapped=True): print(match)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks a bit better. This would give 5 ORFs, which is almost consistent with the example answer, which has 4 sequences.\n",
    "\n",
    "Let's see what happens if I translate these:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Seq import Seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG\n",
      "\n",
      "ATGTAG\n",
      "M\n",
      "ATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAA\n",
      "MGMTPRLGLESLLE\n",
      "ATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAA\n",
      "MTPRLGLESLLE\n"
     ]
    }
   ],
   "source": [
    "print(test_forward)\n",
    "print()\n",
    "\n",
    "for match in regex.finditer(r\"%s\" % pattern, test_forward, overlapped=True):\n",
    "    test_sequence = Seq(match.captures()[0])\n",
    "    print(test_sequence)\n",
    "    print(test_sequence.translate().strip(\"*\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, that looks pretty nice. Now on the other strand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATGCTACTCGGATCATTCAGGCTTATTCCAAAAGAGACTCTAATCCAAGTCGCGGGGTCATCCCCATGTAACCTGAGTTAG\n",
      "MLLGSFRLIPKETLIQVAGSSPCNLS\n",
      "ATGTAA\n",
      "M\n"
     ]
    }
   ],
   "source": [
    "for match in regex.finditer(r\"%s\" % pattern, test_reverse, overlapped=True):\n",
    "    test_sequence = Seq(match.captures()[0])\n",
    "    print(test_sequence)\n",
    "    print(test_sequence.translate().strip(\"*\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That also looks pretty good. Now I probably want to throw all these sequences into a set, to remove duplicates. (The example shows only one time 'M'.) So the final function should look something like...:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pattern = \"ATG([ACGT]{3})*?(TAA|TAG|TGA)\"\n",
    "\n",
    "#Explanation:\n",
    "# 1. Start with \"ATG\" literally (start codon)\n",
    "# 2.1. Then match any number of triplets of A, C, G and/or G\n",
    "# 2.2. However, since '*' is a 'greedy' match, it will look for the longest possible match\n",
    "#     and I want to stop at the first stop codon (= point 3). Therefore, I added the '?'\n",
    "#     to make the matching method 'lazy'. (Also see https://regex101.com/)\n",
    "# 3. Stop the match at either \"TAA\", \"TAG\" or \"TGA\" (stop codons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def translate_six_frames2(seq_record):\n",
    "    \"\"\"\n",
    "Given a Biopython sequence record with a DNA (or RNA?) sequence,\n",
    "translate into amino acid (protein) sequences in six frames.\n",
    "Returns translations as (deduplicated) list of strings.\n",
    "    \"\"\"\n",
    "    translation_list = []\n",
    "    \n",
    "    for strand, nuc in [(+1, seq_record.seq), (-1,  seq_record.seq.reverse_complement())]:\n",
    "        for match in regex.finditer(r\"%s\" % pattern, str(nuc), overlapped=True):\n",
    "            sequence = Seq(match.captures()[0])\n",
    "            translation = sequence.translate().strip(\"*\")\n",
    "            \n",
    "            translation_list.append(translation)\n",
    "            \n",
    "    return(list(set(translation_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_list2 = translate_six_frames2(test_record)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MTPRLGLESLLE\n",
      "MLLGSFRLIPKETLIQVAGSSPCNLS\n",
      "M\n",
      "MGMTPRLGLESLLE\n"
     ]
    }
   ],
   "source": [
    "for translation in test_list2: print(translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, let's see if this is enough to pass the test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MLDPGLRTGCCIPTQRFWALASRHRFRSGGLSCQQLR\n",
      "MQQPVLSPGSNMVC\n",
      "MPWLIYSYFHRAMEMTTGQGLAR\n",
      "MNIRHS\n",
      "MFIFQRVYLGYLSSCAEVDSLENAYSARRLTPTYHRSNYCVCSLYRQDDVAVAPRSMTTGQGLAR\n",
      "MLTVITGKDGAVSLGST\n",
      "MSSVGDL\n",
      "MSLSTPGKHVGR\n",
      "MSVNEHPSRLSWPIPPTVIAEQVFWLSCEGATMFGLFRLIRFAKSAPCSPSSYYAGSSEHTSNASSSWPSIPARS\n",
      "MVC\n",
      "MEMTTGQGLAR\n",
      "MLRGATATSS\n",
      "MSQHNNY\n",
      "MFGLFRLIRFAKSAPCSPSSYYAGSSEHTSNASSSWPSIPARS\n",
      "MCRG\n",
      "MTTGQGLAR\n",
      "MLAHRNIPVTRVVLGRRSPPGRNSSGVMSSVGDL\n",
      "MLICLGEIQRVLWLQSDGDRDIALQTMSVNEHPSRLSWPIPPTVIAEQVFWLSCEGATMFGLFRLIRFAKSAPCSPSSYYAGSSEHTSNASSSWPSIPARS\n",
      "MTLQLPLAA\n",
      "MVQQMPWLIYSYFHRAMEMTTGQGLAR\n",
      "MPKIVAWGCNNQFLVLDLTWYANGYHG\n",
      "MFR\n",
      "MSTRAVYPGQYLQQL\n",
      "MAPK\n",
      "MASTQVLPRLTAPSLPVITVSIPC\n",
      "MSVEM\n",
      "MSTGQCPKSLRGDATTSS\n",
      "MV\n",
      "MSDVHLPTCLPGVLKLMCRG\n",
      "MMASTQVLPRLTAPSLPVITVSIPC\n",
      "M\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sam/miniconda3/envs/jupyter/lib/python3.8/site-packages/Bio/Seq.py:175: BiopythonWarning: Biopython Seq objects now use string comparison. Older versions of Biopython used object comparison. During this transition, please use hash(id(my_seq)) or my_dict[id(my_seq)] if you want the old behaviour, or use hash(str(my_seq)) or my_dict[str(my_seq)] for the new string hashing behaviour.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "seq_record = read_single_sequence_file(\"data/rosalind_orf.txt\")\n",
    "\n",
    "translation_list = translate_six_frames2(seq_record)\n",
    "\n",
    "for translation in translation_list:\n",
    "    print(translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Success!!\n",
    "\n",
    "It worked. I passed the test."
   ]
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}