_7_mux_stats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the filenames from the MinION to generate\n",
    "some interesting statistics about our run!\n",
    "\n",
    "The first thing we'll do is test to see if the QC held true.\n",
    "If the initial QC says 1041 pores, then how many of those actually\n",
    "produce at least one viable sequence?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set paths\n",
    "# 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",
    "PASS_DIR = PORECAMPAU_DATA_DIR + \"fast5/e_coli_R9/\"\n",
    "MY_MUX_DIR = HOME_DIR + \"my_mux_data/\"\n",
    "\n",
    "if not os.path.isdir(MY_MUX_DIR):\n",
    "    os.mkdir(MY_MUX_DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Collect the fast5 files\n",
    "fast5_files = [PASS_DIR + fast5_file for fast5_file in os.listdir(PASS_DIR)\n",
    "               if fast5_file.endswith(\".fast5\")]\n",
    "\n",
    "# We will now write to file for each fastfile, the channel id, and the mux id.\n",
    "mux_file = MY_MUX_DIR + \"e_coli_R7.3.mux.txt\"\n",
    "mux_handle = open(mux_file, 'w')\n",
    "\n",
    "for fast5file in fast5_files:\n",
    "    channel = fast5file.split('_')[-3]  # 'channel_512'\n",
    "    channel_num = channel.replace('ch', \"\")  # '512'\n",
    "    read = fast5file.split('_')[-2]  # 'read_23042'\n",
    "    read_num = read.replace(\"read\", \"\")  # '23042'\n",
    "    # Open fast5 file\n",
    "    f = h5py.File(fast5file, 'r')\n",
    "    mux_path = \"/Raw/Reads/Read_\" + str(read_num)\n",
    "    mux = f[mux_path].attrs.values()[2]\n",
    "    mux_handle.write(channel_num + \"\\t\" + str(mux) + \"\\n\")\n",
    "\n",
    "mux_handle.close()  # Don't forget to close the gate or all your bytes will escape!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The list that we import has as many rows as we have fast5 files, and two columns.\n",
    "The first column is the channel id, the second is the mux id.\n",
    "Open up the mux file to get a better look. It will only be 100 lines long as we only have 100 reads\n",
    "in this folder. For your benefit I have a similar text file using all the pass reads from a run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SHARED_MUX_DIRECTORY = PORECAMPAU_DATA_DIR + \"mux_counts/\"\n",
    "# List files in the mux directory that we can play with\n",
    "mux_files = [mux_file for mux_file in os.listdir(SHARED_MUX_DIRECTORY)]\n",
    "print '\\n'.join(mux_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R -i SHARED_MUX_DIRECTORY,MY_MUX_DIR -o mux_table,mux_stats\n",
    "\n",
    "# copy and paste one of the mux files printed above into to quotation marks below.\n",
    "mux_file = paste(SHARED_MUX_DIRECTORY, \"e_coli_R9.4.mux.txt\", sep=\"/\") \n",
    "\n",
    "mux_data <- read.table(mux_file, header=FALSE, sep=\"\\t\")\n",
    "names(mux_data) <- c(\"channel\", \"mux\")\n",
    "mux_data$channel = factor(mux_data$channel, levels=c(1:512))\n",
    "\n",
    "THRESHOLD = 1;  # How many pass reads need to pass through the pore\n",
    "                # to consider this a pass read.\n",
    "    \n",
    "total_viable_pores = 0;  # Initialise value then increment.\n",
    "\n",
    "mux_per_channel = NULL  # array which we will append to then summarise\n",
    "                      # for each channel. Will be given the number\n",
    "                      # 0:4 for a given channel.\n",
    "        \n",
    "reads_per_pore = NULL  # for each pore a value of the number of reads\n",
    "                          # produced by this channel.\n",
    "\n",
    "mux_table <- as.matrix(table(mux_data$channel, mux_data$mux))\n",
    "\n",
    "for(channel in 1:512){\n",
    "    mux_in_this_channel = sum(mux_table[channel, ] > THRESHOLD)  # number between 0 and 4.\n",
    "    mux_per_channel = c(mux_per_channel, mux_in_this_channel)\n",
    "}\n",
    "\n",
    "\n",
    "# Create a histogram using the mux_per_channel variable.\n",
    "breaks = seq(-0.5, 4.5, by = 1)\n",
    "xlab = \"Number of pores used in channel.\"\n",
    "ylab = \"Number of channels.\"\n",
    "main = \"Histogram of the number of pores used.\"\n",
    "hist(mux_per_channel, breaks=breaks,\n",
    "    xlab=xlab, ylab=ylab, main=main)\n",
    "\n",
    "# Save plot using cairo\n",
    "plot_name = paste(MY_MUX_DIR, \"e_coli_R7.3.png\", sep=\"/\")\n",
    "png(plot_name, type=\"cairo\")\n",
    "hist(mux_per_channel, breaks=breaks,\n",
    "    xlab=xlab, ylab=ylab, main=main)\n",
    "dev.off()\n",
    "\n",
    "# Create a log histogram of the number of pass reads\n",
    "# each channel produced.\n",
    "\n",
    "mux_stats = table(mux_per_channel)  # Return this value so we can print\n",
    "# mux_stats will produce a \n",
    "# a final QC using python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use the output from R to mimic the initial QC stage.\n",
    "# How does this plot compare to the initial QC you performed.\n",
    "for i in range(0,5):\n",
    "    try:\n",
    "        print \"The number of channels with %d viable pores: %d\" % (i, mux_stats[i])\n",
    "    except IndexError:\n",
    "        print(\"No channels present with %d viable pores\" % i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examining half lives:\n",
    "So now we know how many pores in a flowcells actually produced reads.\n",
    "But of those that produced reads, how many reads did they produce?\n",
    "A really cool plot to create is to build a half life plot of the reads.\n",
    "This can also be useful for run diagnostics, good DNA generally allows a pore to last longer ==>\n",
    "more reads can go through the pore. Poor DNA can clog up a pore and we should see a much steeper slope in the number of reads each pore produces.\n",
    "\n",
    "In this next plot we should see a nice logarithmic curve showing the half-life of pores.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R -i mux_table,MY_MUX_DIR\n",
    "\n",
    "# Given we're wanting to see the longevity of a pore, we'll take out those pores\n",
    "# that did not produce at least one read.\n",
    "mux_table <- mux_table[which(mux_table != 0)]\n",
    "\n",
    "run_name = \"e_coli_R7.3\"\n",
    "plot_filename = paste(MY_MUX_DIR, run_name, \".pore_durability.png\", sep=\"\")\n",
    "\n",
    "# Create the plot below and also print it to file\n",
    "hist(mux_table,  main = paste(\"Pore Quality Distribution of\", run_name),\n",
    "        xlab = \"Number of pass reads produced by pore\", ylab = \"Number of pores\")\n",
    "\n",
    "png(plot_filename, type = \"cairo\")\n",
    "hist(mux_table,  main = paste(\"Pore Quality Distribution of\", run_name),\n",
    "        xlab = \"Number of pass reads produced by pore\", ylab = \"Number of pores\")\n",
    "dev.off()"
   ]
  }
 ],
 "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
}