_1_Basic_poretools_command.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This is jupyter\n",
    "\n",
    "* In Jupyter, if you see the symbol In \\[\\*\\], that means that the process is running.\n",
    "* Or it means that the process is waiting to be run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import os\n",
    "import sys\n",
    "import subprocess\n",
    "import h5py\n",
    "from IPython.display import Image\n",
    "from IPython.core.display import HTML\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We will be using poretools. Let's see which poretools we are using.\n",
    "poretools_path = subprocess.check_output(\"which poretools\", shell=True).rstrip()\n",
    "poretools_version = subprocess.check_output(\"poretools --version\", shell=True, stderr=subprocess.STDOUT).rstrip()\n",
    "print \"The poretools path is %s.\\nThe version of poretools is %s\" % (poretools_path, poretools_version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Now we need to set our directories:\n",
    "HOME_DIR = os.environ.get('PORECAMPAU_ANALYSIS_PATH') + \"/\"\n",
    "PORECAMPAU_DATA_DIR = os.environ.get('PORECAMPAU_DATA_PATH') + \"/\"\n",
    "\n",
    "SAMPLE_NAME = \"e_coli_R9\"\n",
    "FAST5_DIR = PORECAMPAU_DATA_DIR + \"fast5/%s/\" % SAMPLE_NAME\n",
    "QUALITY_METRICS_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"quality_metrics/\"\n",
    "PLOTS_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"plots/\"\n",
    "METADATA_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"meta/\"\n",
    "STATS_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"stats/\"\n",
    "SQUIGGLE_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"squiggle/\"\n",
    "FASTQ_DIR = os.path.abspath(HOME_DIR) + \"/\" + \"fastq/\"\n",
    "    \n",
    "sub_dirs = (QUALITY_METRICS_DIR, PLOTS_DIR, METADATA_DIR, STATS_DIR, SQUIGGLE_DIR, FASTQ_DIR)\n",
    "for directory in sub_dirs:\n",
    "    if not os.path.isdir(directory):\n",
    "        os.mkdir(directory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this first part of the tutorial, we will be using poretools to generate some plots of our data to go in the plot folder.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hist_file = PLOTS_DIR + SAMPLE_NAME + \".hist.png\"\n",
    "hist_command_options = [\"hist\"]\n",
    "hist_command_options.append(\"--saveas %s\" % hist_file)\n",
    "hist_command_options.append(FAST5_DIR)\n",
    "hist_command = \"poretools %s\" % ' '.join(hist_command_options)\n",
    "hist_stderr = subprocess.check_output(hist_command, shell=True, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's have a look at the plot we have generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Image(filename=hist_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm.... this is because we've only got a subset of files in the fast5 folder.\n",
    "Let's use a different file from the list below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# List the set of images in the images directory\n",
    "hist_dir = PORECAMPAU_DATA_DIR + \"hist/\"\n",
    "for dirpath, dirnames, filenames in os.walk(hist_dir):\n",
    "    if len(filenames) == 0:  # empty folder\n",
    "        continue\n",
    "    for filename in filenames:\n",
    "        if not filename.endswith(\".png\"):  # Not a histogram image\n",
    "            continue\n",
    "        print dirpath + \"/\" + filename\n",
    "        \n",
    "# Copy this filename into the image loader in the previous code block."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see a peak at around 7kb which is what we would expect for a nanopore run.  \n",
    "You should expect the histogram plot to look similar to that of your tapestation prior to your library preparation.\n",
    "Excessive smaller fragments suggest that unintentional shearing occurred during your library preparation.\n",
    "Remember that long DNA is fragile :)\n",
    "We will now look at generating a yield plot, so that we can see the data generated over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "yield_file = PLOTS_DIR + SAMPLE_NAME + \".yield.png\"\n",
    "yield_command_options = [\"yield_plot\"]\n",
    "yield_command_options.append(\"--saveas %s\" % yield_file)\n",
    "yield_command_options.append(\"--plot-type basepairs\")\n",
    "yield_command_options.append(FAST5_DIR)\n",
    "yield_command = \"poretools %s\" % ' '.join(yield_command_options)\n",
    "\n",
    "yield_stderr = subprocess.check_output(yield_command, shell=True, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Image(yield_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's collect some quantitative stats for our data. We will again use poretools to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Do the same thing for the yield file:\n",
    "yield_dir = PORECAMPAU_DATA_DIR + \"yield/\"\n",
    "for dirpath, dirnames, filenames in os.walk(yield_dir):\n",
    "    if len(filenames) == 0:  # empty folder\n",
    "        continue\n",
    "    for filename in filenames:\n",
    "        if not filename.endswith(\".png\"):  # Not a histogram image\n",
    "            continue\n",
    "        print dirpath + \"/\" + filename\n",
    "        \n",
    "# Copy this filename into the previous image loader."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stats_file = STATS_DIR + SAMPLE_NAME + \".stats.txt\"\n",
    "stats_command_options = [\"stats\"]\n",
    "stats_command_options.append(FAST5_DIR)\n",
    "stats_command = \"poretools %s > %s\" % (' '.join(stats_command_options), stats_file)\n",
    "stats_stderr = subprocess.check_output(stats_command, shell=True, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use the cat command to print the output of the file to the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stats_output = subprocess.check_output(\"cat %s\" % stats_file, shell=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(stats_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the N25, N50 and N75 mean in this context?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality Metrics with Poretools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "qualpos_file = QUALITY_METRICS_DIR + SAMPLE_NAME + \".qualpos.jpg\"\n",
    "qualpos_command_options = [\"qualpos\"]\n",
    "qualpos_command_options.append(\"--saveas %s\" % qualpos_file)\n",
    "qualpos_command_options.append(FAST5_DIR)\n",
    "\n",
    "qualpos_command = \"poretools %s\" % ' '.join(qualpos_command_options)\n",
    "print(qualpos_command)\n",
    "qualpos_stderr = subprocess.check_output(qualpos_command, shell=True, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now look at the image generated by qualpos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Image(qualpos_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting Fastq data from fast5 files\n",
    "There are quite a few software packages to do this.\n",
    "It should be noted that this process can be exceptionally slow.\n",
    "We incorporate a parallel process here so that we can accelerate the process. If you are unfamiliar with parallel processing, please run through the Jupyter notebook on parallel processing first. Yes, you can run this without doing that, but the following commands won't make much sense and you won't have learned very much.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What the Fast5!\n",
    "\n",
    "Now that we know the tools required to extract and analyse the data. So far it's all been a bit of magic. Let's have a look inside the fast5 files that hold all the information that we've been plotting.\n",
    "If you are unfamiliar with exploring hdf5 files in python I suggest you look at the previous tutorial for guidance on some definitions of attributes and datasets etc.\n",
    "\n",
    "Installed along with anaconda h5py package is many useful command line tools for exploring hdf5 data types. Freeware for all operating systems exists so that the hdf5 file can also be visualised without a commandline interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# In a given MinION reads folder, over 100'000 files can exist. \n",
    "# This can cause many delays and headaches even for those with experience on the commandline.\n",
    "# Caution must be taken when entering directories with many folders. Wild characters, such as * can often breakdown\n",
    "# with so many files in the directory.\n",
    "# We can play it safe with the pipe symbol and the head command.\n",
    "print(FAST5_DIR)\n",
    "p = subprocess.Popen(\"ls %s | head -n 10\" % FAST5_DIR, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "p.wait()\n",
    "stdout, stderr = p.communicate()\n",
    "print(stdout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the pipe allows what has been output by the ls command into the head command before the ls has even finished.\n",
    "Very useful given we only want ten files, not ten thousand! In fact the ls would have been so long that the head command\n",
    "Finished before the ls command throwing an error which we see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(stderr) # Nothing to worry about."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# An easy way to handle python files is to store them in an array or dictionary.\n",
    "# The following command adds the directory onto the start of each filename and stores the file paths in an array\n",
    "\n",
    "fast5files = [FAST5_DIR + fast5file for fast5file in stdout.splitlines()]\n",
    "\n",
    "for fast5file in fast5files[0:3]:  # print the first three commands\n",
    "    print(fast5file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the first fast5file to delve into the hdf5 file structure.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fast5file = fast5files[0]\n",
    "# Firstly we open the fast5 file using a handler. Subsequent analysis will use the handler, not the filename.\n",
    "f = h5py.File(fast5file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "hdf5 files are a bit like a tree, where our filename is the root.\n",
    "To see the first set of branches, we use the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for group in f:\n",
    "    print group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These names aren't too intuitive. We can use the following command to then list the second set of branches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for group in f:\n",
    "    for subgroup in f[group]:\n",
    "        print \"Group: %s \\t Subgroup: %s \" % (group, subgroup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how to iterate through a group we need to use the square brackets to indicate which group we are highlighting.\n",
    "To continue further, we don't add another set of square brackets, but append the subgroup to the group within the current set of brackets. See the example below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for subsubgroup in f[\"Analyses/Basecall_1D_000\"]:\n",
    "    print(subsubgroup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the id parameter to see if each node is a dataset or another group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for subsubsubgroup in f[\"Analyses/Basecall_1D_000\"]:\n",
    "    subsubsubtype = f[\"Analyses/Basecall_1D_000/%s\" % subsubsubgroup].id\n",
    "    print(\"%s \\t %s\" % (subsubsubgroup, subsubsubtype))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what else is in the fast5 file?\n",
    "We'd rather not keep naming variables with infinite subs. Given we don't know how deep this goes, it may be best to generate a recursive list to see an overview of the file structure. Although there exists ways to do this already, by running h5ls on the command line.. why not come up with our own, which is better!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def recursive_list_hdf5(f):\n",
    "    for name in f:\n",
    "        path = str(f).split('\"')[1]  \n",
    "        #print(path)\n",
    "        depth = path.count(\"/\")\n",
    "        #print(str(f[name]))\n",
    "        info = str(f[name]).split(\" \") # Example: <HDF5 dataset \"Events\": shape (7063,), type \"|V78\">\n",
    "        object_type = info[1]          # Or <HDF5 group \"/Analyses/Basecall_1D_000\" (5 members)>\n",
    "        print_order = []\n",
    "        \n",
    "        # We need to indent to show directory structure.\n",
    "        tab_depth = []\n",
    "        for tab in range(0, depth):\n",
    "            tab_depth.append(\"\\t\\t\")\n",
    "        tab_depth = ''.join(tab_depth)\n",
    "        \n",
    "        if object_type == \"group\":\n",
    "            object_name = info[2].split(\"/\")[-1] # We only want the last attribute.\n",
    "            object_member_count = re.sub(\"\\(\",\"\", info[-2])  # Always just before the last space\n",
    "            print_order.append(\"Name: %s\" % re.sub(\"\\\"\", \"\", object_name))\n",
    "            print_order.append(\"Type: Group\")\n",
    "            print_order.append(\"Members: %s \" % object_member_count)\n",
    "        elif object_type == \"dataset\":\n",
    "            object_name  = re.sub(\"\\\"|\\\"|:$\", \"\", info[2])  # Remove any colons from the end of dataset names.\n",
    "            object_length = re.sub(\",$\",\"\", info[4])\n",
    "            print_order.append(\"Name: %s\" % object_name)\n",
    "            print_order.append(\"Type: Dataset\")\n",
    "            print_order.append(\"Shape: %s\" % ''.join(object_length))\n",
    "        \n",
    "        print_order = ', '.join(print_order)\n",
    "        print(tab_depth + print_order)\n",
    "        \n",
    "        # A group will have child nodes, so we iterate through these.\n",
    "        if object_type == \"group\":\n",
    "            recursive_list_hdf5(f[name])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "recursive_list_hdf5(f)\n",
    "#for name in f[\"/Analyses/Basecall_1D_000/BaseCalled_complement\"]:\n",
    "#   print name.parent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aha! We can find the fastq files under: \n",
    "* Analyses/Basecall_2D_000/BaseCalled_2D/Fastq or \n",
    "* Analyses/Basecall_1D_000/BaseCalled_complement/Fastq or\n",
    "* Analyses/Basecall_1D_000/BaseCalled_template/Fastq\n",
    "\n",
    "Note that for extracting datasets we do something a little different to groups.\n",
    "Here we go down to the last group the dataset is in then add the dataset name in its own square brakets.\n",
    "The '[()]' ending allows us to extract the actual dataset, rather than a reference to the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fastq = f[\"Analyses/Basecall_1D_000/BaseCalled_complement/\"]['Fastq'][()]\n",
    "fastq = fastq.replace(\"\\\\n\", \"\\n\")\n",
    "print(fastq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Phoa! You can comment out the print statement above and then re-run to suppress the your output. Extracting fastq data seems quite complex, is fortunate that we have software packages that can extract this data for us. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Running poretools fastq\n",
    "p = subprocess.Popen(\"poretools fastq %s > %s%s.fastq\" % (FAST5_DIR, FASTQ_DIR, SAMPLE_NAME), \n",
    "                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "p.wait()\n",
    "stdout, stderr = p.communicate()\n",
    "print(stdout, stderr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was pretty quick! In future, you can run parallel -j < threads \\> | poretools fastq < dir \\> \n",
    "when running on hundreds of thousands of files!  \n",
    "(It seems to run faster than any thing else I've seen - would be delighted to be disproven!) "
   ]
  }
 ],
 "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
}