Plot_insertsizes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting insert size distributions\n",
    "====================\n",
    "\n",
    "We will use the following tools:\n",
    "\n",
    "* ipython notebook\n",
    "* the python *pysam* module for reading (and writing) SAM/BAM files\n",
    "* matplotlib for plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check where we are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what files we have in the folder. You should have at least one BAM file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pysam\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set some parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "infile = 'map_pe.sorted.bam' #SAM or BAM file\n",
    "max_inserts = 10000 # collect no more than this insert sizes (for speed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the file, go through all alignments, collect the insert size only when it is larger than 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with pysam.Samfile(infile, 'rb' ) as samfile:\n",
    "    # NB 'rb' for BAM file, 'r' for SAM file\n",
    "    # collect all alignments in the SAM/BAM file\n",
    "    alignments = samfile.fetch()\n",
    "    insert_sizes = []\n",
    "    # go through all alignments\n",
    "    for alignment in alignments:\n",
    "        insert_size = alignment.isize\n",
    "        if insert_size > 0:\n",
    "            # collect insert size\n",
    "            insert_sizes.append(insert_size)\n",
    "        if len(insert_sizes) >= max_inserts:\n",
    "            # stop\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a simple histogram of alignments  \n",
    "Use 100 bins and a suitable range for the insert size, e.g. 0-500 for paired end, 0-10000 for mate pair (experiment with this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(insert_sizes, bins=100, range = (0,500))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beautify by adding\n",
    "\n",
    "* limits for the x-axis\n",
    "* a label for the data\n",
    "* a legend\n",
    "\n",
    "Experiment with limits, labels etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(insert_sizes, bins=100, range = (0,500), label = 'Paired end')\n",
    "plt.xlim(200,400)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add labels and save as file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(insert_sizes, bins=100, range = (0,500), label = 'Paired end')\n",
    "plt.xlim(200,400)\n",
    "plt.legend()\n",
    "plt.xlabel('Insert size')\n",
    "plt.ylabel('Count')\n",
    "plt.show()\n",
    "plt.savefig('Inserts.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}