chp3-seq-objects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 2 and 3\n",
    "\n",
    "## Reading fasta files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gi|2765658|emb|Z78533.1|CIZ78533\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet())\n",
      "740\n",
      "gi|2765657|emb|Z78532.1|CCZ78532\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAG...GGC', SingleLetterAlphabet())\n",
      "753\n",
      "gi|2765656|emb|Z78531.1|CFZ78531\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet())\n",
      "748\n",
      "gi|2765655|emb|Z78530.1|CMZ78530\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAAACAACAT...CAT', SingleLetterAlphabet())\n",
      "744\n",
      "gi|2765654|emb|Z78529.1|CLZ78529\n",
      "Seq('ACGGCGAGCTGCCGAAGGACATTGTTGAGACAGCAGAATATACGATTGAGTGAA...AAA', SingleLetterAlphabet())\n",
      "733\n"
     ]
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "fastas = list(SeqIO.parse(\"files/ls_orchid.fasta\",\"fasta\"))  # generator to list\n",
    "for seq_record in fastas[:5]:\n",
    "    print(seq_record.id)\n",
    "    print(repr(seq_record.seq))\n",
    "    print(len(seq_record))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequence objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('AGTACACTGGT')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio.Seq import Seq\n",
    "my_seq = Seq(\"AGTACACTGGT\")\n",
    "my_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('AGTACACTGGT', IUPACProtein())"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio.Alphabet import IUPAC\n",
    "my_prot = Seq(\"AGTACACTGGT\", IUPAC.protein)\n",
    "my_prot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequences as strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 A\n",
      "1 G\n",
      "2 T\n",
      "3 A\n",
      "4 C\n",
      "5 A\n",
      "6 C\n",
      "7 T\n",
      "8 G\n",
      "9 G\n",
      "10 T\n"
     ]
    }
   ],
   "source": [
    "for index, letter in enumerate(my_seq):\n",
    "    print(\"%i %s\" % (index, letter))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11\n"
     ]
    }
   ],
   "source": [
    "print(len(my_seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq.count(\"GT\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "46.875"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio.Alphabet import IUPAC\n",
    "from Bio.SeqUtils import GC\n",
    "my_seq = Seq(\"GATCGATGGGCCTATATAGGATCGAAAATCGC\", IUPAC.unambiguous_dna)\n",
    "GC(my_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slicing sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('GATGGGCC', IUPACUnambiguousDNA())"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq[4:12]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below prints every other 3rd nucleotide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('GCTGTAGTAAG', IUPACUnambiguousDNA())"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq[0::3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python trick, `[::-1]` can be used to get reverse of a string (sequence in our case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('CGCTAAAAGCTAGGATATATCCGGGTAGCTAG', IUPACUnambiguousDNA())"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concatenation\n",
    "\n",
    "protein and dna sequences can not be concatenated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Incompatible alphabets IUPACProtein() and IUPACUnambiguousDNA()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-28-0b91ff3c1125>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmy_prot\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mmy_seq\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/srv/conda/lib/python3.6/site-packages/Bio/Seq.py\u001b[0m in \u001b[0;36m__add__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    309\u001b[0m                 raise TypeError(\n\u001b[1;32m    310\u001b[0m                     \"Incompatible alphabets {0!r} and {1!r}\".format(\n\u001b[0;32m--> 311\u001b[0;31m                         self.alphabet, other.alphabet))\n\u001b[0m\u001b[1;32m    312\u001b[0m             \u001b[0;31m# They should be the same sequence type (or one of them is generic)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m             \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAlphabet\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_consensus_alphabet\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0malphabet\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0malphabet\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Incompatible alphabets IUPACProtein() and IUPACUnambiguousDNA()"
     ]
    }
   ],
   "source": [
    "my_prot + my_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sum` is acting funny, thus `+` is better way to concatenate strings (*in the book, the `sum` was used with list elements their order was okay*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('TTTTAGTC')"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seq1 = Seq(\"AGTC\")\n",
    "seq2 = Seq(\"TTTT\")\n",
    "sum(seq1,seq2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case and (reverse) complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('gatcgatgggcctatataggatcgaaaatcgc', DNAAlphabet())"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq.lower()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq.reverse_complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('ctagctacccggatatatcctagcttttagcg', DNAAlphabet())"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq.lower().complement()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transcription and translation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG\n",
      "CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT\n",
      "AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG\n"
     ]
    }
   ],
   "source": [
    "coding_dna = Seq(\"ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG\", IUPAC.unambiguous_dna)\n",
    "template_dna = coding_dna.reverse_complement()\n",
    "messenger_rna = coding_dna.transcribe()\n",
    "print(\"%s\\n%s\\n%s\" % (coding_dna, template_dna, messenger_rna))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "template_dna.reverse_complement().transcribe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coding_dna.transcribe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('MAIVMGR*KGAR*', HasStopCodon(IUPACProtein(), '*'))"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coding_dna.translate()    # no need for coding_dna.transcribe().translate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('MAIVMGR', IUPACProtein())"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coding_dna.translate(to_stop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MutableSeq objects\n",
    "\n",
    "Seq objects are not mutable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'Seq' object does not support item assignment",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-53-50e848a55f5c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmy_seq\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"G\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: 'Seq' object does not support item assignment"
     ]
    }
   ],
   "source": [
    "my_seq[5] = \"G\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GATCGATGGGCCTATATAGGATCGAAAATCGC\n",
      "GATCGCTGGGCCTATATAGGATCGAAAATCGC\n"
     ]
    }
   ],
   "source": [
    "mutable_seq = my_seq.tomutable()\n",
    "mutable_seq[5] = \"C\"\n",
    "print(\"%s\\n%s\" % (my_seq,mutable_seq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mutating actions, such as `remove()` and `reverse()` can be used with MutableSeq objects.\n",
    "\n",
    "MutableSeq objects can be transferred to Seq (immutable) objects by `toseq()` method. (e.g. `new_seq = mutable_seq.toseq()`)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}