_4_IGV_In_Jupyter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### IGV, Interactive Genome Viewer.\n",
    " Have you used IGV before? Probably. It's an awesome way to see your reads stacked up against a reference!\n",
    " \n",
    " Using the E coli reference genome we will compare the output of the different chemistries, the basecalling algorithms used for the same chemistry and the reads of a nanopolished file.\n",
    " \n",
    " IGV is a great way to see what type of errors are seen in the nanopore. Homopolymers can cause issues as you should be able to observe. Most of the errors in nanopore reads can be seen to be insertion/deletion errors, an issue when determining if a change in the current is due to a nucleotide or just noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os, subprocess, re\n",
    "from igv import IGV, Reference, Track"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Show paths to the bam files\n",
    "reference_file = \"references/e_coli_k12_mg1655/NC_000913.fna\"\n",
    "bam_dir = \"/mnt/shared/PoreCampAU/data/alignment/e_coli_R9/\"\n",
    "my_alignment_dir = \"/home/researcher/alignment/\"\n",
    "# List the set of files in the bam directory\n",
    "for dirpath, dirnames, filenames in os.walk(bam_dir):\n",
    "    if len(filenames) == 0:  # empty folder\n",
    "        continue\n",
    "    for filename in filenames:\n",
    "        if not filename.endswith(\".bam\"):  # Not a bam file, maybe an index file.\n",
    "            continue\n",
    "        print dirpath + \"/\" + filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "IGV(locus=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's easy to visually see the differences in quality between each chemistry and alignment algorithm.\n",
    "However, quantitative metrics are often easier to explain to someone.\n",
    "To do this, we'll use the stats module of samtools to generate a stats report from the bam file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bam_file = \"/mnt/shared/PoreCampAU/data/alignment/e_coli_R9/nanonet2d/2016-11-15_E_COLI_R9_bwa-mem.sorted.bam\"\n",
    "stats_file = my_alignment_dir + \"e_coli_R9_metrichor_stats.txt\"  # rename this for each bam file to stop overwriting.\n",
    "samtools_stats_command = \"samtools stats %s > %s\" % (bam_file, stats_file)\n",
    "stderr = subprocess.check_call(samtools_stats_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr==\"\":\n",
    "    print \"Stderr = %s\" % stderr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, now this file is particularly big for a summary sheet.\n",
    "Fortunately it's sorted into components that we can extract using the grep command.\n",
    "...we could also use python, because python is beautiful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stats_file_handler = open(stats_file, 'r')\n",
    "for line in stats_file_handler:\n",
    "    if line.startswith(\"SN\\t\"):\n",
    "        print(line.rstrip()) # rstrip gets rid of the \\n at the end of the line.\n",
    "stats_file_handler.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}