kmer-abundance.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Remember to go to 'File' and 'Make a copy'\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd /mnt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Change the name below to something unique to you"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "NetID='titus'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mkdir -p $NetID"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd $NetID"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You should now be in '/mnt/USERNAME', where USERNAME is *not* CHANGEME"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Make two random genomes with a 10:1 abundance"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# configure some stuff...\n",
      "import sys\n",
      "sys.path.insert(0, '/usr/local/share/khmer/python')\n",
      "import os\n",
      "os.environ['PYTHONPATH'] = '/usr/local/share/khmer/python'\n",
      "import khmer"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random\n",
      "random.seed(1)\n",
      "\n",
      "x = [\"A\"] + [\"G\"] + [\"C\"] + [\"T\"]\n",
      "x = x*1000\n",
      "random.shuffle(x)\n",
      "x = \"\".join(x)\n",
      "\n",
      "y = [\"A\"] + [\"G\"] + [\"C\"] + [\"T\"]\n",
      "y = y*1000\n",
      "random.shuffle(y)\n",
      "y = \"\".join(y)\n",
      "\n",
      "print 'x is', x[:100]\n",
      "print 'y is', y[:100]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "outfp = open('metagenome.fa', 'w')\n",
      "print >>outfp, \">x 1\"\n",
      "print >>outfp, x\n",
      "print >>outfp, \">y 2\"\n",
      "print >>outfp, y\n",
      "outfp.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!python /usr/local/share/2012-paper-diginorm/pipeline/make-biased-reads.py metagenome.fa | head -100000 > reads.fa"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(Yes, you should see an error.)\n",
      "\n",
      "## Now, calculate a k-mer abundance histogram for k=20"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!/usr/local/share/khmer/scripts/abundance-dist-single.py -k 20 -x 1e8 reads.fa reads.hist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "histdata = numpy.loadtxt('reads.hist')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(histdata[:,0], histdata[:,1])\n",
      "xlabel(\"k-mer abundance\")\n",
      "ylabel(\"N of k-mers with that abundance\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(histdata[:,0], histdata[:,1])\n",
      "axis(ymax=500)\n",
      "xlabel(\"k-mer abundance\")\n",
      "ylabel(\"N of k-mers with that abundance\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(histdata[:,0], histdata[:,1])\n",
      "axis(xmax=10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Questions to answer\n",
      "\n",
      "How do the 3 peaks (1, 100, 800) shift relative to each other with (a) changing coverage and (b) different k values?\n",
      "\n",
      "The '-k' parameter is set up above in the abundance-dist-single command; the coverage is set by how many sequences you keep ('head -10000') in the make-biased-reads.py command."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}