python_fm_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with Python: functions and modules\n",
    "\n",
    "## Session 4: Using third party libraries\n",
    "\n",
    "- [Matplotlib](#Matplotlib)\n",
    "- [Exercise 4.1](#Exercise-4.1)\n",
    "- [BioPython](#BioPython)\n",
    "- [Working with sequences](#Working-with-sequences)\n",
    "- [Connecting with biological databases](#Connecting-with-biological-databases)\n",
    "- [Exercise 4.2](#Exercise-4.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matplotlib\n",
    "\n",
    "[matplotlib](http://matplotlib.org/) is probably the single most used Python package for graphics. It provides both a very quick way to visualize data from Python and publication-quality figures in many formats.\n",
    "\n",
    "matplotlib.pyplot is a collection of command style functions that make matplotlib work like MATLAB. Each pyplot function makes some change to a figure: e.g., creates a figure, creates a plotting area in a figure, plots some lines in a plotting area, decorates the plot with labels, etc. \n",
    "\n",
    "Let's start with a very simple plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as mpyplot\n",
    "mpyplot.plot([1,2,3,4])\n",
    "mpyplot.ylabel('some numbers')\n",
    "mpyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plot()` is a versatile command, and will take an arbitrary number of arguments. For example, to plot x versus y, you can issue the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpyplot.plot([1,2,3,4], [1,4,9,16])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For every x, y pair of arguments, there is an **optional third argument** which is the format string that indicates the color and line type of the plot. The letters and symbols of the format string are from MATLAB, and you concatenate a color string with a line style string. The default format string is `b-`, which is a solid blue line. For example, to plot the above with red circles, you would chose `ro`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as mpyplot\n",
    "mpyplot.plot([1,2,3,4], [1,4,9,16], 'ro')\n",
    "mpyplot.axis([0, 6, 0, 20])\n",
    "mpyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`matplotlib` has a few methods in the **`pyplot` module** that make creating common types of plots faster and more convenient because they automatically create a Figure and an Axes object. The most widely used are:\n",
    "\n",
    "- [mpyplot.bar](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.bar) – creates a bar chart.\n",
    "- [mpyplot.boxplot](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.boxplot) – makes a box and whisker plot.\n",
    "- [mpyplot.hist](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist) – makes a histogram.\n",
    "- [mpyplot.plot](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot) – creates a line plot.\n",
    "- [mpyplot.scatter](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.scatter) – makes a scatter plot.\n",
    "\n",
    "Calling any of these methods will automatically setup `Figure` and `Axes` objects, and draw the plot. Each of these methods has different parameters that can be passed in to modify the resulting plot.\n",
    "\n",
    "The [Pyplot tutorial](http://matplotlib.org/users/pyplot_tutorial.html) is where these simple examples above are coming from. More could be learn from it if you wish during your own time.\n",
    "\n",
    "Let's now try to plot the GC content along the chain we have calculated during the previous session, while solving the Exercises 3.3 and 3.4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq = 'ATGGTGCATCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTG'\n",
    "gc = [40.0, 60.0, 80.0, 60.0, 40.0, 60.0, 40.0, 40.0, 40.0, 60.0, \n",
    "      40.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, \n",
    "      60.0, 40.0, 40.0, 40.0, 40.0, 40.0, 60.0, 60.0, 80.0, 80.0, \n",
    "      80.0, 60.0, 40.0, 40.0, 20.0, 40.0, 60.0, 80.0, 80.0, 80.0, \n",
    "      80.0, 60.0, 60.0, 60.0, 80.0, 80.0, 100.0, 80.0, 60.0, 60.0, \n",
    "      60.0, 40.0, 60.0]\n",
    "window_ids = range(len(gc))\n",
    "\n",
    "import matplotlib.pyplot as mpyplot\n",
    "mpyplot.plot(window_ids, gc, '--' )\n",
    "mpyplot.xlabel('5 bases window id along the sequence')\n",
    "mpyplot.ylabel('%GC')\n",
    "mpyplot.title('GC plot for sequence\\n' + seq)\n",
    "mpyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.1\n",
    "\n",
    "Re-use the GapMinder dataset to plot, in Jupyter using Matplotlib, from the world data the life expectancy against GDP per capita for 1957 and 2007 using a scatter plot, add title to your graph as well as a legend."
   ]
  },
  {
   "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."
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating sequence\n",
    "from Bio.Seq import Seq\n",
    "my_seq = Seq(\"AGTACACTGGT\")\n",
    "print(my_seq)\n",
    "print(my_seq[10])\n",
    "print(my_seq[1:5])\n",
    "print(len(my_seq))\n",
    "print(my_seq.count(\"A\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use functions from `Bio.SeqUtils` to get idea about a sequence "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the molecular weight\n",
    "from Bio.SeqUtils import GC, molecular_weight\n",
    "print(GC(my_seq))\n",
    "print(molecular_weight(my_seq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One letter code protein sequences can be converted into three letter codes using `seq3` utility "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.SeqUtils import seq3\n",
    "print(seq3(my_seq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alphabets defines how the strings are going to be treated as sequence object. `Bio.Alphabet` module defines the available alphabets for Biopython. `Bio.Alphabet.IUPAC` provides basic definition for DNA, RNA and proteins. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Alphabet import IUPAC\n",
    "my_dna = Seq(\"AGTACATGACTGGTTTAG\", IUPAC.unambiguous_dna)\n",
    "print(my_dna)\n",
    "print(my_dna.alphabet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dna.complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dna.reverse_complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dna.translate()"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Writing FASTA files\n",
    "from Bio import SeqIO\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Seq import Seq\n",
    "from Bio.Alphabet import IUPAC\n",
    "\n",
    "sequence = 'MYGKIIFVLLLSEIVSISASSTTGVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFG'\n",
    " \n",
    "seq = Seq(sequence, IUPAC.protein)\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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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 4.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "- Retrieve 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": [
    "## Congratulation! You reached the end! "
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}