ngs-5x-digital-normalization.ipynb

{
 "metadata": {
  "name": "ngs-5x-digital-normalization"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "source": [
      "You'll need to install khmer and screed; see ngs-02 and ngs-03 notebooks."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "cd /mnt"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## First, grab some code from the digital normalization paper for making simulated genomes & reads."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!git clone git://github.com/ged-lab/2012-paper-diginorm.git"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "mkdir kmer"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "cd kmer"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Now, make a random genome of length 8kb (4 x 2000)"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import random",
      "random.seed(1)",
      "",
      "x = [\"A\"] + [\"G\"] + [\"C\"] + [\"T\"]",
      "x = x*2000",
      "random.shuffle(x)",
      "x = \"\".join(x)",
      "x"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "fp = open('genome.fa', 'w')",
      "fp.write('>genome\\n')",
      "fp.write(x)",
      "fp.close()"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Make reads with a 1% error rate to a coverage of about 200"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /mnt/2012-paper-diginorm/pipeline/make-reads.py genome.fa > reads.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!head reads.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Now, use khmer to load all the 20-mers in the reads into a counting data structure.",
      "",
      "The counting data structure is saved as counts.kh.  Other parameters explained in:",
      "",
      "http://readthedocs.org/docs/khmer/en/latest/scripts.html"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/scripts/load-into-counting.py -x 1e8 -N 4 -k 20 counts.kh reads.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Get the abundance distribution of the k-mers in the reads"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/scripts/abundance-dist.py -s counts.kh reads.fa reads.dist"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## The output is in 'reads.dist'",
      "",
      "Column 0 is the k-mer count; column 1 is the number of k-mers with that count; column 2 is the running total of column 1; and column 3 is the fraction of totla in column 2."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!head reads.dist"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Now, plot the abundance distribution"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x = []",
      "y = []",
      "for line in open('reads.dist'):",
      "    count, num, total, frac = line.split()",
      "    x.append(int(count))",
      "    y.append(int(num))"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Reads with errors have many k-mers with abundance=1."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "plot(x, y)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## There will also be a bunch of high-count k-mers correlated with coverage"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "plot(x, y)",
      "axis(ymax=1000)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Apply digital normalization to downsample to a coverage of 20 @ k=20"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/scripts/normalize-by-median.py -x 2e8 -N 4 -k 20 -C 20 reads.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Count the k-mers in the digitally normalized reads"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/scripts/load-into-counting.py -x 1e8 -N 4 -k 20 counts-dn.kh reads.fa.keep"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/scripts/abundance-dist.py -s counts-dn.kh reads.fa.keep reads.dn.dist"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!head reads.dn.dist"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Compare the raw and diginormed reads"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x2 = []",
      "y2 = []",
      "for line in open('reads.dn.dist'):",
      "    count, num, total, frac = line.split()",
      "    x2.append(int(count))",
      "    y2.append(int(num))"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "plot(x2, y2)",
      "axis(ymax=2000)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Now let's look at mapping-based measures of coverage",
      "",
      "### First, install bowtie"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!curl -O -L http://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12.7/bowtie-0.12.7-linux-x86_64.zip",
      "!unzip bowtie-0.12.7-linux-x86_64.zip",
      "!cp bowtie-0.12.7/bowtie* /usr/local/bin"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "### Index the reference genome"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!bowtie-build genome.fa genome"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Map the raw and abundance filtered reads to the reference genome"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!bowtie -f -a --strata --best genome reads.fa > reads.map",
      "!bowtie -f -a --strata --best genome reads.fa.keep > reads.map.keep"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!head -3 reads.map"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Build a function to read in the mapping output file"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "# define a function that we can use to read in the mapping output file",
      "def read_mapping_file(filename):",
      "    for line in open(filename):",
      "        _, orient, refname, pos, read, _, _, mismatches = line.split('\\t')[:8]",
      "        pos = int(pos)",
      "        mismatches = mismatches.split(',')",
      "        x = []",
      "        for posn in mismatches:",
      "            posn = posn.strip()",
      "            if posn:",
      "                posn = posn.split(':')[0]",
      "                posn = int(posn)",
      "                x.append(posn)",
      "        yield orient, refname, pos, read, x"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Plot *coverage* distribution (using 'hist')"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "genome_pos = {}",
      "",
      "for (orient, refname, pos, read, mismatches) in read_mapping_file('reads.map'):",
      "    for readi in range(len(read)):",
      "        genome_pos[pos + readi] = genome_pos.get(pos + readi, 0) + 1",
      "",
      "genome_pos_keep = {}",
      "",
      "for (orient, refname, pos, read, mismatches) in read_mapping_file('reads.map.keep'):",
      "    for readi in range(len(read)):",
      "        genome_pos_keep[pos + readi] = genome_pos_keep.get(pos + readi, 0) + 1",
      "        "
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Coverage distribution of original reads"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "hist(genome_pos.values(), bins=100)",
      "xlabel('coverage')",
      "ylabel('number of bases with that coverage')"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Coverage distribution of normalized reads"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "hist(genome_pos_keep.values(), bins=100)",
      "xlabel('coverage')"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      ""
     ],
     "language": "python",
     "outputs": []
    }
   ]
  }
 ]
}