MRNA-Inferring_mRNA_from_protein.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferring mRNA from protein, from [Rosalind.info](https://www.rosalind.info)\n",
    "\n",
    "(Text copied from http://rosalind.info/problems/mrna/)\n",
    "\n",
    "<div class=\"problem-statement problem-statement-bordered\" problem=\"130\">\n",
    "    <blockquote>\n",
    "<h2 id=\"pitfalls-of-reversing-translation\">Pitfalls of Reversing Translation</h2>\n",
    "<p>When researchers discover a new <a class=\"term\" href=\"/glossary/protein/\" id=\"term-206\" rel=\"tooltip\" title=\"The functional unit of the cell.\">protein</a>, they would like to infer\n",
    "the strand of <a class=\"term\" href=\"/glossary/messenger-rna/\" id=\"term-266\" rel=\"tooltip\" title=\"\n",
    "An RNA molecule that serves as the blueprint for translation into protein.\">mRNA</a> from which this protein could have been <a class=\"term\" href=\"/glossary/translation/\" id=\"term-247\" rel=\"tooltip\" title=\"\n",
    "The process by which mRNA is converted into a peptide chain for the creation of a protein.\">translated</a>,\n",
    "thus allowing them to locate <a class=\"term\" href=\"/glossary/gene/\" id=\"term-363\" rel=\"tooltip\" title=\"\n",
    "An interval of DNA whose nucleotides are translated into a polypeptide for protein creation.\">genes</a> associated with this protein on the <a class=\"term\" href=\"/glossary/genome/\" id=\"term-360\" rel=\"tooltip\" title=\"\n",
    "The collection of all of an organism's DNA taken from all its chromosomes.\">genome</a>.</p>\n",
    "<p>Unfortunately, although any <a class=\"term\" href=\"/glossary/rna-string/\" id=\"term-350\" rel=\"tooltip\" title=\"\n",
    "A string constructed from the alphabet {A, C, G, U}.\">RNA string</a> can be translated into a unique <a class=\"term\" href=\"/glossary/protein-string/\" id=\"term-371\" rel=\"tooltip\" title=\"\n",
    "A string composed of symbols taken from the English alphabet less B, J, O, U, X, and Z;\n",
    "representing a peptide chain formed from amino acids.\">protein string</a>,\n",
    "reversing the process yields a huge number of possible RNA strings from a single\n",
    "protein string because most amino acids correspond to multiple RNA <a class=\"term\" href=\"/glossary/codon/\" id=\"term-259\" rel=\"tooltip\" title=\"\n",
    "A triplet of contiguous nucleotides.\">codons</a>\n",
    "(see the <a class=\"term\" href=\"/glossary/rna-codon-table/\" id=\"term-217\" rel=\"tooltip\" title=\"A table indicating the translation of individual RNA codons into amino acids for the\n",
    "purpose of protein creation.\">RNA Codon Table</a>).</p>\n",
    "<p>Because of memory considerations, most data formats that are built into languages have\n",
    "upper bounds on how large an integer can be: in some versions of Python,\n",
    "an \"int\" variable may be required to be no larger than <mathjax>$2^{31} -1$</mathjax>, or 2,147,483,647.\n",
    "As a result, to deal with very large numbers in Rosalind, we need to devise a system\n",
    "that allows us to manipulate large numbers without actually having to store large numbers.</p>\n",
    "</blockquote>\n",
    "<h2 id=\"problem\">Problem</h2>\n",
    "<p>For positive integers <mathjax>$a$</mathjax> and <mathjax>$n$</mathjax>, <mathjax>$a$</mathjax> <a class=\"term new\" href=\"/glossary/modular-arithmetic/\" id=\"term-577\" rel=\"tooltip\" title=\"New term: \n",
    "The study of arithmetic on integer remainders.\">modulo</a> <mathjax>$n$</mathjax> (written <mathjax>$a\\mod n$</mathjax> in shorthand) is the remainder\n",
    "when <mathjax>$a$</mathjax> is divided by <mathjax>$n$</mathjax>.  For example, <mathjax>$29 \\mod 11 = 7$</mathjax> because <mathjax>$29 = 11 \\times 2 + 7$</mathjax>.</p>\n",
    "<p><a class=\"term new\" href=\"/glossary/modular-arithmetic/\" id=\"term-577\" rel=\"tooltip\" title=\"New term: \n",
    "The study of arithmetic on integer remainders.\">Modular arithmetic</a> is the study of addition, subtraction, multiplication, and division\n",
    "with respect to the modulo operation.  We say that <mathjax>$a$</mathjax> and <mathjax>$b$</mathjax> are <a class=\"term new\" href=\"/glossary/modular-arithmetic/\" id=\"term-577\" rel=\"tooltip\" title=\"New term: \n",
    "The study of arithmetic on integer remainders.\">congruent</a> modulo <mathjax>$n$</mathjax>\n",
    "if <mathjax>$a \\mod n = b \\mod n$</mathjax>; in this case, we use the notation <mathjax>$a \\equiv b \\mod n$</mathjax>.</p>\n",
    "<p>Two useful facts in modular arithmetic are that if <mathjax>$a \\equiv b \\mod n$</mathjax> and <mathjax>$c \\equiv d \\mod n$</mathjax>,\n",
    "then <mathjax>$a+c \\equiv b+d \\mod n$</mathjax> and <mathjax>$a \\times c \\equiv b \\times d \\mod n$</mathjax>.  To check your understanding of these rules,\n",
    "you may wish to verify these relationships for <mathjax>$a = 29$</mathjax>, <mathjax>$b = 73$</mathjax>, <mathjax>$c = 10$</mathjax>, <mathjax>$d = 32$</mathjax>, and <mathjax>$n = 11$</mathjax>.</p>\n",
    "<p>As you will see in this exercise, some Rosalind problems will ask for a (very large)\n",
    "integer solution modulo a smaller number to avoid the computational pitfalls that arise with\n",
    "storing such large numbers.</p>\n",
    "<p><span class=\"given-return\">Given:</span> A <a class=\"term\" href=\"/glossary/protein-string/\" id=\"term-371\" rel=\"tooltip\" title=\"\n",
    "A string composed of symbols taken from the English alphabet less B, J, O, U, X, and Z;\n",
    "representing a peptide chain formed from amino acids.\">protein string</a> of length at most 1000 <a class=\"term\" href=\"/glossary/amino-acid/\" id=\"term-198\" rel=\"tooltip\" title=\"\n",
    "The monomer unit for proteins; the same 20 amino acids commonly occur in most species.\">aa</a>.</p>\n",
    "<p><span class=\"given-return\">Return:</span> The total number of different RNA strings from which the protein could have been translated, modulo 1,000,000.\n",
    "(Don't neglect the importance of the <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> in protein translation.)</p>\n",
    "<h2 id=\"sample-dataset\">Sample Dataset</h2>\n",
    "<div class=\"codehilite\"><pre>MA\n",
    "</pre></div>\n",
    "\n",
    "\n",
    "<h2 id=\"sample-output\">Sample Output</h2>\n",
    "<div class=\"codehilite\"><pre>12\n",
    "</pre></div>\n",
    "\n",
    "\n",
    "<blockquote>\n",
    "<h2 id=\"hint\">Hint</h2>\n",
    "<p>What does it mean intuitively to take a number modulo 1,000,000?</p>\n",
    "</blockquote>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## My interpretation/reasoning\n",
    "\n",
    "1. First the hint: Intuitively, to me this means that the final number is smaller than the modulo number. So a (big) number modulo 1,000,000 will return a number smaller than 1,000,000, which is safely below the integer limit of $2^{31}-1$.\n",
    "\n",
    "2. This problem is about 'backtranslating' amino acid sequences (proteins) into nucleotide sequences (mRNA).\n",
    "\n",
    "3. However, the exercise does not ask to make actual mRNA sequences.\n",
    "\n",
    "4. Instead, I get to calculate how many options exist for the backtranslation.\n",
    "\n",
    "5. For some amino acid there is only one codon, for some there are multiple. (Note that there are also multiple stop codons! This is not written in the amino acid sequence, but does add 3 possibilities.)\n",
    "\n",
    "6. The final answer should be the product of the possibilities per amino acid, e.g. \"MA\" reads: 1 codon for M * 4 codons for A * 3 stop codons = 1 * 4 * 3 = 12 possible mRNA sequences.\n",
    "\n",
    "_Note: Rosalind has included this RNA codon table as reminder:_\n",
    "\n",
    "```\n",
    "UUU F      CUU L      AUU I      GUU V\n",
    "UUC F      CUC L      AUC I      GUC V\n",
    "UUA L      CUA L      AUA I      GUA V\n",
    "UUG L      CUG L      AUG M      GUG V\n",
    "UCU S      CCU P      ACU T      GCU A\n",
    "UCC S      CCC P      ACC T      GCC A\n",
    "UCA S      CCA P      ACA T      GCA A\n",
    "UCG S      CCG P      ACG T      GCG A\n",
    "UAU Y      CAU H      AAU N      GAU D\n",
    "UAC Y      CAC H      AAC N      GAC D\n",
    "UAA Stop   CAA Q      AAA K      GAA E\n",
    "UAG Stop   CAG Q      AAG K      GAG E\n",
    "UGU C      CGU R      AGU S      GGU G\n",
    "UGC C      CGC R      AGC S      GGC G\n",
    "UGA Stop   CGA R      AGA R      GGA G\n",
    "UGG W      CGG R      AGG R      GGG G\n",
    "```\n",
    "\n",
    "To put this into a script, I want to have:\n",
    " - A list of codons per amino acid\n",
    " - A function to read the amino acid sequence (from a text file)\n",
    " - Then for each amino acid in the sequence, I want to multiply the number of possibilities by the corresponding number of possibilities (starting with 1)\n",
    " - The final product should be _first_ multiplied by 3 to add in the stop codon, _then_ take the modulo 1,000,000 to produce the final answer.\n",
    " \n",
    "Sounds not too difficult, right? Let's hope it work out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_protein_sequence(input_file):\n",
    "    \"\"\"\n",
    "Read a protein sequence from an input file.\n",
    "Return the sequence as string (text).\n",
    "\n",
    "Note: I am assuming exactly 1 line of text.\n",
    "    \"\"\"\n",
    "    with open(input_file, 'r') as read_file:\n",
    "        for line in read_file:\n",
    "            sequence = line.strip()\n",
    "            \n",
    "    return(sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MA\n"
     ]
    }
   ],
   "source": [
    "test_sequence = read_protein_sequence(\"data/Example_Inferring_mRNA_from_protein.txt\")\n",
    "\n",
    "print(test_sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64\n"
     ]
    }
   ],
   "source": [
    "number_of_codons = {\n",
    "    \"A\": 4,\n",
    "    \"C\": 2,\n",
    "    \"D\": 2,\n",
    "    \"E\": 2,\n",
    "    \"F\": 2,\n",
    "    \"G\": 4,\n",
    "    \"H\": 2,\n",
    "    \"I\": 3,\n",
    "    \"K\": 2,\n",
    "    \"L\": 6,\n",
    "    \"M\": 1,\n",
    "    \"N\": 2,\n",
    "    \"P\": 4,\n",
    "    \"Q\": 2,\n",
    "    \"R\": 6,\n",
    "    \"S\": 6,\n",
    "    \"T\": 4,\n",
    "    \"V\": 4,\n",
    "    \"W\": 1,\n",
    "    \"Y\": 2,\n",
    "    \"Stop\": 3\n",
    "}\n",
    "\n",
    "#As a little check: the sum of these numbers should be 4 * 16 = 64:\n",
    "\n",
    "print(sum(number_of_codons.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "from operator import mul\n",
    "#Use these functions to multiple the contents of a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_mrna_possibilities(aa_sequence):\n",
    "    \"\"\"\n",
    "Given a protein sequence (amino acids),\n",
    "calculate the number of mRNA (nucleotide) sequences\n",
    "that can translate into that protein sequence.\n",
    "    \"\"\"\n",
    "    possibilities = 3 * reduce(mul, [\n",
    "        number_of_codons[amino_acid] for amino_acid in aa_sequence\n",
    "        ], 1)\n",
    "    \n",
    "    return(possibilities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12\n"
     ]
    }
   ],
   "source": [
    "test_possibilities = calculate_mrna_possibilities(test_sequence)\n",
    "\n",
    "print(test_possibilities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far so good. Ready to test on the real sequence!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "OverflowError",
     "evalue": "int too large to convert to float",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOverflowError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-0b57bab204ac>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mpossibilities\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_mrna_possibilities\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msequence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0manswer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpossibilities\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;36m1e6\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0manswer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mOverflowError\u001b[0m: int too large to convert to float"
     ]
    }
   ],
   "source": [
    "sequence = read_protein_sequence(\"data/rosalind_mrna.txt\")\n",
    "\n",
    "possibilities = calculate_mrna_possibilities(sequence)\n",
    "\n",
    "answer = possibilities % 1e6\n",
    "\n",
    "print(answer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right, an error. Isn't that just what we were warned for?\n",
    "Maybe we need the modulo earlier!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_mrna_possibilities_modulo(aa_sequence):\n",
    "    \"\"\"\n",
    "Given a protein sequence (amino acids),\n",
    "calculate the number of mRNA (nucleotide) sequences\n",
    "that can translate into that protein sequence.\n",
    "    \"\"\"\n",
    "    possibilities = 3 * reduce(mul, [\n",
    "        number_of_codons[amino_acid] for amino_acid in aa_sequence\n",
    "        ], 1) % 1e6\n",
    "    \n",
    "    return(possibilities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "OverflowError",
     "evalue": "int too large to convert to float",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOverflowError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-4e2e193a05be>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0msequence\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mread_protein_sequence\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"data/rosalind_mrna.txt\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mpossibilities\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_mrna_possibilities_modulo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msequence\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpossibilities\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-13-5df04834278b>\u001b[0m in \u001b[0;36mcalculate_mrna_possibilities_modulo\u001b[0;34m(aa_sequence)\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mthat\u001b[0m \u001b[0mcan\u001b[0m \u001b[0mtranslate\u001b[0m \u001b[0minto\u001b[0m \u001b[0mthat\u001b[0m \u001b[0mprotein\u001b[0m \u001b[0msequence\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \"\"\"\n\u001b[0;32m----> 7\u001b[0;31m     possibilities = 3 * reduce(mul, [\n\u001b[0m\u001b[1;32m      8\u001b[0m         \u001b[0mnumber_of_codons\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mamino_acid\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mamino_acid\u001b[0m \u001b[0;32min\u001b[0m \u001b[0maa_sequence\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m         ], 1) % 1e6\n",
      "\u001b[0;31mOverflowError\u001b[0m: int too large to convert to float"
     ]
    }
   ],
   "source": [
    "sequence = read_protein_sequence(\"data/rosalind_mrna.txt\")\n",
    "\n",
    "possibilities = calculate_mrna_possibilities_modulo(sequence)\n",
    "\n",
    "print(possibilities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this is not yet a solution... Let me think how to handle this then...\n",
    "\n",
    "Right! I think I need to spend some more attention to **modular arithmetic** and **congruent modulo**.  \n",
    "I should be able to use those to solve this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or else try this first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_mrna_possibilities_modulo2(aa_sequence):\n",
    "    \"\"\"\n",
    "Given a protein sequence (amino acids),\n",
    "calculate the number of mRNA (nucleotide) sequences\n",
    "that can translate into that protein sequence.\n",
    "    \"\"\"\n",
    "    possibilities = 3 * reduce(mul, [\n",
    "        number_of_codons[amino_acid] for amino_acid in aa_sequence\n",
    "        ], 1) % int(1e6)\n",
    "    \n",
    "    return(possibilities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "623808\n"
     ]
    }
   ],
   "source": [
    "print(calculate_mrna_possibilities_modulo2(sequence))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So when the big `int` does not have to be converted to a `float`, which was necessary in this case because `1e6` is a `float` and to modulo those numbers, they need to be the same type. So that big product had to be converted. But `1e6` can be an `int` just as well, so if that is converted to an `int`, the calculation can be done and I get an answer.\n",
    "\n",
    "This answer seems a lot easier than figuring out how to incorporate modular arithmetic. Perhaps it feels a bit like cheating, but the main goal was to make a working script, right? Not necessarily to do it by the one path that may have been intended by whoever made this exercise.\n",
    "\n",
    "Let's go again!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35456\n"
     ]
    }
   ],
   "source": [
    "new_sequence = read_protein_sequence(\"data/rosalind_mrna2.txt\")\n",
    "\n",
    "new_possibilities = calculate_mrna_possibilities_modulo2(new_sequence)\n",
    "\n",
    "print(new_possibilities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Success!!\n",
    "\n",
    "That worked."
   ]
  }
 ],
 "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
}