23_python_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data science in Python\n",
    "\n",
    "- Course GitHub repo: https://github.com/pycam/python-data-science\n",
    "- Python website: https://www.python.org/ \n",
    "\n",
    "## Session 2.3: Biological data with BioPython\n",
    "\n",
    "\n",
    "- [BioPython](#BioPython)\n",
    "- [Working with sequences](#Working-with-sequences)\n",
    "- [Connecting with biological databases](#Connecting-with-biological-databases)\n",
    "- [Exercise 2.3.1](#Exercise-2.3.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BioPython\n",
    "\n",
    "The goal of Biopython is to make it as easy as possible to use Python for bioinformatics by creating high-quality, reusable modules and classes. Biopython features include parsers for various Bioinformatics file formats (BLAST, Clustalw, FASTA, Genbank,...), access to online services (NCBI, Expasy,...), interfaces to common and not-so-common programs (Clustalw, DSSP, MSMS...), a standard sequence class, various clustering modules, a KD tree data structure etc. and documentation as well as a tutorial: http://biopython.org/DIST/docs/tutorial/Tutorial.html.\n",
    "\n",
    "The Biopython Application Programming Interface (API) documentation: http://biopython.org/DIST/docs/api/."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a sequence by defining a `Seq` object with strings. `Bio.Seq()` takes as input a string and converts in into a Seq object. We can print the sequences, individual residues, lengths and use other functions to get summary statistics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGTACATGACTGGTTTAG\n"
     ]
    }
   ],
   "source": [
    "from Bio.Seq import Seq\n",
    "\n",
    "my_dna = Seq(\"AGTACATGACTGGTTTAG\")\n",
    "print(my_dna)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Seq class *inherets* from string, and so all string methods work with it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGTACATGACTGGTTTAG\n",
      "T\n",
      "GTAC\n",
      "18\n",
      "5\n",
      "AGTACATGACTGGTTTAGAGTACATGACTGGTTTAG\n"
     ]
    }
   ],
   "source": [
    "print(my_dna)\n",
    "print(my_dna[10])\n",
    "print(my_dna[1:5])\n",
    "print(len(my_dna))\n",
    "print(my_dna.count(\"A\"))\n",
    "print(my_dna+my_dna)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TCATGTACTGACCAAATC\n"
     ]
    }
   ],
   "source": [
    "my_dna.complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dna.reverse_complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('ST*LV*', HasStopCodon(IUPACProtein(), '*'))"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_dna.translate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use functions from `Bio.SeqUtils` to get idea about a sequence "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "38.888888888888886\n",
      "5633.600099999999\n"
     ]
    }
   ],
   "source": [
    "# Calculate the molecular weight\n",
    "from Bio.SeqUtils import GC, molecular_weight\n",
    "print(GC(my_dna))\n",
    "print(molecular_weight(my_dna))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One letter code protein sequences can be converted into three letter codes using `seq3` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AlaGlyThrAlaCysAlaCysThrGlyGlyThr\n"
     ]
    }
   ],
   "source": [
    "from Bio.SeqUtils import seq3\n",
    "print(seq3(my_dna))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parsing sequence file format: FASTA files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequence files can be parsed and read the same way we read other files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">swissprot|P02724|GLPA_HUMAN Glycophorin-A;\n",
      "MYGKIIFVLLLSEIVSISASSTTGVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAH\n",
      "EVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKK\n",
      "SPSDVKPLPSPDTDVPLSSVEIENPETSDQ\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with open( \"data/glpa.fa\" ) as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Biopython provides specific functions to allow parsing/reading sequence files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "swissprot|P02724|GLPA_HUMAN\n",
      "MYGKIIFVLLLSEIVSISASSTTGVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ\n"
     ]
    }
   ],
   "source": [
    "# Reading FASTA files\n",
    "from Bio import SeqIO\n",
    "\n",
    "with open(\"data/glpa.fa\") as f:\n",
    "    for protein in SeqIO.parse(f, 'fasta'):\n",
    "        print(protein.id)\n",
    "        print(protein.seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequence objects can be written into files using file handles with the function `SeqIO.write()`. We need to provide the name of the output sequence file and the sequence file format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">THEID a description\n",
      "MYGKIIFVLLLSEIVSISASSTTGVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAH\n",
      "EVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFG\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Writing FASTA files\n",
    "from Bio import SeqIO\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Seq import Seq\n",
    "\n",
    "sequence = 'MYGKIIFVLLLSEIVSISASSTTGVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFG'\n",
    " \n",
    "seq = Seq(sequence)\n",
    "protein = [SeqRecord(seq, id=\"THEID\", description='a description'),]\n",
    "\n",
    "with open( \"biopython.fa\", \"w\") as f:\n",
    "    SeqIO.write(protein, f, 'fasta')\n",
    "\n",
    "with open( \"biopython.fa\" ) as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connecting with biological databases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequences can be searched and downloaded from public databases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DQ091202.1 with 6 features\n",
      "TTCTGGGCCTCAGTTTCCTCATTTGTATAATAACAGAATTGGAGAGTAAATTCTTAAGAGGCTTACCAGGCTGTAATTCTAAAAAAAATGCATAAATAAACTTGCCAAGGCAGATGTTTTTAGCAGCAATTCCTGAAAGAAACGGGACCAGGAGATAAGTAGAGAAAGAGTGAAGGTCTGAAATCAAACTAATAAGACAGTCCCAGACTGTCAAGGAGAGGTATGGCTGTCATCATTCAGGCCTCACCCTGCAGAACCACACCCTGGCCTTGGCCAATCTGCTCACAAGAGCAAAAAGGGCAGGACCAGGGTTGGGCATATAAGGAAGAGTAGTGCCAGCTGCTGTTTACACTCACTTCTGACACAACTGTGTTGACTAGCAACTACCCAATCAGACACCATGGTGAATCTGACTGCTGCTGAGAAGACACAAGTCACCAACCTGTGGGGCAAGGTGAATGTGAAAGAGCTTGGTGGTGAGGCCCTGAGCAGGTTTGTATCTAGGTTGCAAGGTAGACTTAAGGAGGGTTGAGTGGGGCTGGGCATGTGGAGACAGAACAGTCTCCCAGTTTCTGACAGGCACTGACTTCCTCTGCACCSTGTGGTGCTTTCACCTTCAGGCTGCTGGTGGTCTACCCATGGACCCGGAGGTTCTTTGAACACTTTGGGGACCTGTCCACTGCTGACGCTGTCCTGCACAACGCTAAAGTGCTGGCCCATGGCGAGAAAGTGTTGACCTCCTTTGGTGAGGGCCTGAAGCACCTGGACAACCTCAAGGGCACCTTTGCCGATCTGAGCGAGCTGCACTGTGACAAGCTGCACGTGGATCCTGAGAATTTCAGGGTGAGTCTAGGAGACACTCTATTTTTTCTTTTCACTTTGTAGTCTTTCACTGTGATTATTTTGCTTATTTGAATTTCCTCTGTATCTCTTTTTACTCGACTATGTTTCATCATTTAGTGTTTTTTCAACTTATACCATTTTGTATTACTTTTCTTTCAATATTCTTCCTTTTTTCCTGACTCACATTCTTGCTTTATATCATGCTCTTTATTTAATTTCCTACGTTTTTGCTCTTGCTCTCCCTTTCTCCTAGTTTCCTTCCCTCTGAACAGTACCCAAATTGTGCATACCACCTCTCGTCCACTATTTCTGCACTGGGGCAAATCCCCACCCCTCCTCCATATGAGGGTTGGAAAGGACTGAATCAAAGAGGAGAGGATCATGGTGCTGTTCTAGAGTATGTGATTCATTTCAGACTTGAAGGATAACTTGAATAATATAAAATCAGGAGTAAATGGAGAGGAAAGTCAGTATCTGAGAATGAAAGATCAGAAGGTCATAGACGAGATGGGGAGCAGAAGTTACTAAGAAACTGACCATTGTGGCTATAATTAATCACTTAATTAGTTAATTAATATGTTTGTTATTTATTCACGTTTTTCATTTTGGTGGGAGTAAATTTGGGCTAGTGTGTGGGCAACATAAATGGGTTTCACCCCATTGTCTCAGAGGCCAAGCTGGATTGCTTTGTTAACCATGTCTGTGTATGTATCTACCTCTTCCCCATAGCTCCTGGGCAATGTGCTGGTGATTGTCCTGGCCCGCCACTTTGGCAAGGAATTCACCCCAGATGTTCAGGCTGCCTATGAGAAGGTTGTGGCAGGTGTGGCGAATGCCCTGGCTCACAAATACCACTGAGATCCTGGCCTGTTCCTGGTATCCATCGGAAGCCCCATTTCCCGAGATGCTATCTCTGAATTTGGGAAAATAATGCCAACTCTCAAGGGCATCTCTTCTGCCTAATAAAGTACTTTCAGCTCAACTTTCTGATTCATTTATTTTTTTTCTCAGTCACTCTTGTGGTGGGGGAAGTTCCCAAGGCTCTATGGACAGAGAGCTCTTGTGCCTTATAGGAAAAGTTCAAGGGAAATTGGAAAATAAAGGGAACCATACACAGATATTAATGGGAACAATTCTACTTCAAAGGCATAAAGATTGGGAAGGTTTGGCAAATAGGATACTGGTACTACAGGGATTCCATGGGCCTCAGGCCTAAGACATAGCCCCAGGGCTAACTTTCAGATTCAATTCCAGAAATTACTCACAAAATAATGGA\n",
      ">DQ091202.1 Elephas maximus HBB/D gene, complete cds\n",
      "TTCTGGGCCTCAGTTTCCTCATTTGTATAATAACAGAATTGGAGAGTAAATTCTTAAGAG\n",
      "GCTTACCAGGCTGTAATTCTAAAAAAAATGCATAAATAAACTTGCCAAGGCAGATGTTTT\n",
      "TAGCAGCAATTCCTGAAAGAAACGGGACCAGGAGATAAGTAGAGAAAGAGTGAAGGTCTG\n",
      "AAATCAAACTAATAAGACAGTCCCAGACTGTCAAGGAGAGGTATGGCTGTCATCATTCAG\n",
      "GCCTCACCCTGCAGAACCACACCCTGGCCTTGGCCAATCTGCTCACAAGAGCAAAAAGGG\n",
      "CAGGACCAGGGTTGGGCATATAAGGAAGAGTAGTGCCAGCTGCTGTTTACACTCACTTCT\n",
      "GACACAACTGTGTTGACTAGCAACTACCCAATCAGACACCATGGTGAATCTGACTGCTGC\n",
      "TGAGAAGACACAAGTCACCAACCTGTGGGGCAAGGTGAATGTGAAAGAGCTTGGTGGTGA\n",
      "GGCCCTGAGCAGGTTTGTATCTAGGTTGCAAGGTAGACTTAAGGAGGGTTGAGTGGGGCT\n",
      "GGGCATGTGGAGACAGAACAGTCTCCCAGTTTCTGACAGGCACTGACTTCCTCTGCACCS\n",
      "TGTGGTGCTTTCACCTTCAGGCTGCTGGTGGTCTACCCATGGACCCGGAGGTTCTTTGAA\n",
      "CACTTTGGGGACCTGTCCACTGCTGACGCTGTCCTGCACAACGCTAAAGTGCTGGCCCAT\n",
      "GGCGAGAAAGTGTTGACCTCCTTTGGTGAGGGCCTGAAGCACCTGGACAACCTCAAGGGC\n",
      "ACCTTTGCCGATCTGAGCGAGCTGCACTGTGACAAGCTGCACGTGGATCCTGAGAATTTC\n",
      "AGGGTGAGTCTAGGAGACACTCTATTTTTTCTTTTCACTTTGTAGTCTTTCACTGTGATT\n",
      "ATTTTGCTTATTTGAATTTCCTCTGTATCTCTTTTTACTCGACTATGTTTCATCATTTAG\n",
      "TGTTTTTTCAACTTATACCATTTTGTATTACTTTTCTTTCAATATTCTTCCTTTTTTCCT\n",
      "GACTCACATTCTTGCTTTATATCATGCTCTTTATTTAATTTCCTACGTTTTTGCTCTTGC\n",
      "TCTCCCTTTCTCCTAGTTTCCTTCCCTCTGAACAGTACCCAAATTGTGCATACCACCTCT\n",
      "CGTCCACTATTTCTGCACTGGGGCAAATCCCCACCCCTCCTCCATATGAGGGTTGGAAAG\n",
      "GACTGAATCAAAGAGGAGAGGATCATGGTGCTGTTCTAGAGTATGTGATTCATTTCAGAC\n",
      "TTGAAGGATAACTTGAATAATATAAAATCAGGAGTAAATGGAGAGGAAAGTCAGTATCTG\n",
      "AGAATGAAAGATCAGAAGGTCATAGACGAGATGGGGAGCAGAAGTTACTAAGAAACTGAC\n",
      "CATTGTGGCTATAATTAATCACTTAATTAGTTAATTAATATGTTTGTTATTTATTCACGT\n",
      "TTTTCATTTTGGTGGGAGTAAATTTGGGCTAGTGTGTGGGCAACATAAATGGGTTTCACC\n",
      "CCATTGTCTCAGAGGCCAAGCTGGATTGCTTTGTTAACCATGTCTGTGTATGTATCTACC\n",
      "TCTTCCCCATAGCTCCTGGGCAATGTGCTGGTGATTGTCCTGGCCCGCCACTTTGGCAAG\n",
      "GAATTCACCCCAGATGTTCAGGCTGCCTATGAGAAGGTTGTGGCAGGTGTGGCGAATGCC\n",
      "CTGGCTCACAAATACCACTGAGATCCTGGCCTGTTCCTGGTATCCATCGGAAGCCCCATT\n",
      "TCCCGAGATGCTATCTCTGAATTTGGGAAAATAATGCCAACTCTCAAGGGCATCTCTTCT\n",
      "GCCTAATAAAGTACTTTCAGCTCAACTTTCTGATTCATTTATTTTTTTTCTCAGTCACTC\n",
      "TTGTGGTGGGGGAAGTTCCCAAGGCTCTATGGACAGAGAGCTCTTGTGCCTTATAGGAAA\n",
      "AGTTCAAGGGAAATTGGAAAATAAAGGGAACCATACACAGATATTAATGGGAACAATTCT\n",
      "ACTTCAAAGGCATAAAGATTGGGAAGGTTTGGCAAATAGGATACTGGTACTACAGGGATT\n",
      "CCATGGGCCTCAGGCCTAAGACATAGCCCCAGGGCTAACTTTCAGATTCAATTCCAGAAA\n",
      "TTACTCACAAAATAATGGA\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Read FASTA file from NCBI GenBank\n",
    "from Bio import Entrez\n",
    "\n",
    "Entrez.email = 'A.N.Other@example.com' # Always tell NCBI who you are\n",
    "handle = Entrez.efetch(db=\"nucleotide\", id=\"71066805\", rettype=\"gb\")\n",
    "seq_record = SeqIO.read(handle, \"gb\")\n",
    "handle.close()\n",
    "\n",
    "print(seq_record.id, 'with', len(seq_record.features), 'features')\n",
    "print(seq_record.seq)\n",
    "print(seq_record.format(\"fasta\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-1558745ff4e7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mhandle\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExPASy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_sprot_raw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'HBB_HUMAN'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mprot_record\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSeqIO\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"swiss\"\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      6\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SeqIO/__init__.py\u001b[0m in \u001b[0;36mread\u001b[0;34m(handle, format, alphabet)\u001b[0m\n\u001b[1;32m    724\u001b[0m     \u001b[0miterator\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mformat\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[1;32m    725\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 726\u001b[0;31m         \u001b[0mfirst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0miterator\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    727\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mStopIteration\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    728\u001b[0m         \u001b[0mfirst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SeqIO/__init__.py\u001b[0m in \u001b[0;36mparse\u001b[0;34m(handle, format, alphabet)\u001b[0m\n\u001b[1;32m    659\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Unknown format '%s'\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    660\u001b[0m         \u001b[0;31m# This imposes some overhead... wait until we drop Python 2.4 to fix it\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 661\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mr\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mi\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    662\u001b[0m             \u001b[0;32myield\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    663\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SeqIO/SwissIO.py\u001b[0m in \u001b[0;36mSwissIterator\u001b[0;34m(handle)\u001b[0m\n\u001b[1;32m     83\u001b[0m     \"\"\"\n\u001b[1;32m     84\u001b[0m     \u001b[0mswiss_records\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSwissProt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhandle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 85\u001b[0;31m     \u001b[0;32mfor\u001b[0m \u001b[0mswiss_record\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mswiss_records\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     86\u001b[0m         \u001b[0;31m# Convert the SwissProt record to a SeqRecord\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m         \u001b[0mseq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSeq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSeq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mswiss_record\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msequence\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mAlphabet\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgeneric_protein\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SwissProt/__init__.py\u001b[0m in \u001b[0;36mparse\u001b[0;34m(handle)\u001b[0m\n\u001b[1;32m    138\u001b[0m     \"\"\"\n\u001b[1;32m    139\u001b[0m     \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 140\u001b[0;31m         \u001b[0mrecord\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhandle\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    141\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mrecord\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    142\u001b[0m             \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SwissProt/__init__.py\u001b[0m in \u001b[0;36m_read\u001b[0;34m(handle)\u001b[0m\n\u001b[1;32m    253\u001b[0m             \u001b[0m_read_kw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrecord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    254\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'FT'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 255\u001b[0;31m             \u001b[0m_read_ft\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrecord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mline\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    256\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'SQ'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    257\u001b[0m             \u001b[0mcols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\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;32m~/miniconda3/envs/general/lib/python3.7/site-packages/Bio/SwissProt/__init__.py\u001b[0m in \u001b[0;36m_read_ft\u001b[0;34m(record, line)\u001b[0m\n\u001b[1;32m    582\u001b[0m         \u001b[0mdescription\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m29\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m70\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    583\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# is continuation of last one\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 584\u001b[0;31m         \u001b[0;32massert\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mfrom_res\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mto_res\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    585\u001b[0m         \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfrom_res\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mto_res\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mold_description\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mold_ft_id\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrecord\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfeatures\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    586\u001b[0m         \u001b[0;32mdel\u001b[0m \u001b[0mrecord\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfeatures\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Read SWISSPROT record\n",
    "from Bio import ExPASy\n",
    "\n",
    "handle = ExPASy.get_sprot_raw('HBB_HUMAN')\n",
    "prot_record = SeqIO.read(handle, \"swiss\")\n",
    "handle.close()\n",
    "\n",
    "print(prot_record.description)\n",
    "print(prot_record.seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2.3.1\n",
    "\n",
    "-  a FASTA file named `data/sample.fa` using BioPython and answer the following questions:\n",
    "  - How many sequences are in the file?\n",
    "  - What are the IDs and the lengths of the longest and the shortest sequences?\n",
    "  - Select sequences longer than 500bp. What is the average length of these sequences?\n",
    "  - Calculate and print the percentage of GC in each of the sequences.\n",
    "  - Write the newly created sequences into a FASTA file named `long_sequences.fa` "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieving more biological data (live coding session)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next session\n",
    "\n",
    "Go to our next notebook: [Session 2.4: Data project report in Jupyter](24_python_data.ipynb)"
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}