{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practical 5: Understanding GQ and DP in sequence data\n",
    "\n",
    "## The cells of this practical can be entered (by cut and paste) into the IPython console.\n",
    "\n",
    "## Before entering the first cell, make sure you have changed to the directory hail-practical. Skip the first cell haven't closed IPython console since running the last practical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "from math import log, isnan\n",
    "from pprint import pprint\n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "from hail import *\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "def qqplot(pvals):\n",
    "    spvals = sorted([x for x in pvals if x and not(isnan(x))])\n",
    "    exp = [-log(float(i) / len(spvals), 10) for i in np.arange(1, len(spvals) + 1, 1)]\n",
    "    obs = [-log(p, 10) for p in spvals]\n",
    "    plt.scatter(exp, obs)\n",
    "    plt.plot(np.arange(0, max(max(exp), max(obs))), c=\"red\")\n",
    "    plt.xlabel(\"Expected p-value (-log10 scale)\")\n",
    "    plt.ylabel(\"Observed p-value (-log10 scale)\")\n",
    "    plt.xlim(xmin=0)\n",
    "    plt.ylim(ymin=0)\n",
    "\n",
    "hc = HailContext()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vds = hc.read('1kg.vds').cache()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Puzzle: Investigating the GQ distribution\n",
    "\n",
    "## GQ refresher\n",
    "\n",
    "GQ is \"genotype quality\", which is roughly the log-scaled probability that your genotype is called wrong.\n",
    "\n",
    "GQ 10 means 90% confidence.\n",
    "\n",
    "GQ 20 means 99% confidence.\n",
    "\n",
    "GQ 30 means 99.9% confidence.\n",
    "\n",
    "GQ is truncated at 99, which corresponds to a ~ 0.0000000001 chance that your call is wrong (if it's calibrated!)\n",
    "|\n",
    "If we plot a histogram of GQ values, what will it look like?\n",
    "\n",
    "### DP refresher\n",
    "\n",
    "DP is \"depth\", which is the total number of reads for a given sample at a given variant.\n",
    "\n",
    "### Produce a histogram of GQ and DP values for every genotype in our dataset\n",
    "\n",
    "### We'll use the `hist` functions again here. You _definitely_ can't collect all the GQ values for a full dataset: with 1K Genomes, you'd get 200 billion values back (and you probably don't have 1600 gigabytes of RAM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[gq_hist, dp_hist] = vds.query_genotypes(['gs.map(g => g.gq).hist(0, 100, 100)', \n",
    "                                          'gs.map(g => g.dp).hist(0, 30, 30)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.xlim(0, 100)\n",
    "plt.ylim(0, 2500000)\n",
    "plt.xlabel('GQ')\n",
    "plt.ylabel('Count')\n",
    "plt.title('GQ Histogram')\n",
    "plt.bar(gq_hist.binEdges[:-1], gq_hist.binFrequencies, width=1, label='GQ')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.xlim(0, 30)\n",
    "plt.ylim(0, 3500000)\n",
    "plt.xlabel('DP')\n",
    "plt.ylabel('Count')\n",
    "plt.title('DP Histogram')\n",
    "plt.bar(dp_hist.binEdges[:-1], dp_hist.binFrequencies, width=1, label='DP')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This GQ histogram is pretty strange. We're going to learn why.\n",
    "\n",
    "### By eye, we can identify at least 4 superimposed distributions. We need to pull these apart. Separating heterozygotes from homozygotes is a good place to start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[het_gq_hist, hom_gq_hist] = vds.query_genotypes(['gs.filter(g => g.isHet).map(g => g.gq).hist(0, 100, 100)', \n",
    "                                                  'gs.filter(g => !g.isHet).map(g => g.gq).hist(0, 100, 100)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.xlim(0, 100)\n",
    "plt.ylim(0, 2000000)\n",
    "plt.xlabel('GQ')\n",
    "plt.ylabel('Count')\n",
    "plt.title('GQ Histogram')\n",
    "plt.bar(het_gq_hist.binEdges[:-1], het_gq_hist.binFrequencies, width=1, color='red', label='heterozygotes')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.xlim(0, 100)\n",
    "plt.ylim(0, 2500000)\n",
    "plt.xlabel('GQ')\n",
    "plt.ylabel('Count')\n",
    "plt.title('GQ Histogram')\n",
    "plt.bar(hom_gq_hist.binEdges[:-1], hom_gq_hist.binFrequencies, width=1, label='homozygotes')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Separating the heterzogotes helped a bit, but doesn't address why we see three superimposed histograms with very different magnitudes.\n",
    "\n",
    "### Let's go back to depth to investigate -- remember that depth is a more 'primitive' piece of metadata, and is used to produce GQ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# argmax is the index of the largest value in a list\n",
    "from numpy import argmax\n",
    "gq_mode = int(hom_gq_hist.binEdges[argmax(hom_gq_hist.binFrequencies)])\n",
    "dp_mode = int(dp_hist.binEdges[argmax(dp_hist.binFrequencies)])\n",
    "\n",
    "print('GQ mode is %d' % gq_mode)\n",
    "print('DP mode is %d' % dp_mode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## There are 3 superimposed GQ distributions for homozygotes.\n",
    "\n",
    "## The ratio between the mode GQ and DP is 3.\n",
    "\n",
    "## We can visually assess correlation by looking at a histogram of the ratio between the two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gq_dp_hist = vds.query_genotypes('gs.filter(g => !g.isHet).map(g => g.gq / g.dp).hist(1, 7, 12)')\n",
    "\n",
    "plt.xlim(0, 8)\n",
    "plt.ylim(0, 20000000)\n",
    "plt.xlabel('GQ/DP')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('GQ / DP Histogram')\n",
    "plt.bar(gq_dp_hist.binEdges[:-1], gq_dp_hist.binFrequencies, width=0.5, label='homozygote GQ/DP')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This ratio is extremely consistent! Remember also that DP is inherently quantized due to the process of short-read sequencing.\n",
    "\n",
    "### To learn more about where this ratio comes from, continue to the code below!.\n",
    "\n",
    "### Since GQ is roughly the probability that the second-most-likely genotype is actually correct, for homozygotes it's the probabilitly that the genotype was actually heterozygous.\n",
    "\n",
    "### One simple model is a binomial: if we have _r_ reads, then the probability we call a homozygote wrong is the probability of flipping _r_ heads from a fair coin (this is how we model seeing a reference read at every read)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import binom\n",
    "from numpy import log10\n",
    "\n",
    "def binomial_p(num_reads):\n",
    "    return binom(num_reads, 0.5).pmf(0)\n",
    "\n",
    "# phred scaling is -10 * log10(x)\n",
    "def phred_scale(x):\n",
    "    return -10 * log10(x)\n",
    "                   \n",
    "\n",
    "xs = xrange(30)\n",
    "ys = [phred_scale(binomial_p(x)) for x in xs]\n",
    "observed_ratio = [3 * x for x in xs]\n",
    "\n",
    "plt.scatter(xs, ys, label='binomial probability of 0 reads')\n",
    "plt.plot(xs, map(lambda value: 3 * value, xs), color='k', label='slope=3')\n",
    "plt.xlim(0, 30)\n",
    "plt.ylim(0, 100)\n",
    "plt.xlabel('Reads')\n",
    "plt.ylabel('phred-scaled GQ')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why is the ratio 3?  One more read (+1 DP) forces us to flip one more coin, multiplying our probability by 0.5. When we phred scale this..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "- 10 * log10(0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus question (hard): Why do heterozygotes have such high GQ values?\n",
    "\n",
    "### GQ is truncated at 99, but you can see the true value by looking at the 2nd-smallest PL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "true_het_gq = vds.query_genotypes('gs.filter(g => g.isHet).map(g => min(g.pl[0], g.pl[2])).hist(0, 400, 100)')\n",
    "plt.xlim(0, 400)\n",
    "plt.ylim(0, 500000)\n",
    "plt.xlabel('True GQ')\n",
    "plt.ylabel('Count')\n",
    "plt.title('True GQ Histogram')\n",
    "plt.bar(true_het_gq.binEdges[:-1], true_het_gq.binFrequencies, width=4, color='red', label='heterozygotes')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is the binomial model appropriate?  What other sources of uncertainty are there?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: calculate the mean depth and mean GQ of the entire dataset.\n",
    "\n",
    "Fill in the `<?>` with code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[mean_gq, mean_dp] = vds.query_genotypes(['<?>', '<?>'])\n",
    "print(\"mean GQ is %f, mean DP is %f\" % (mean_gq, mean_dp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: what fraction of genotypes have GQ 20 or above?  What fraction of homozygotes?  What fraction of heterozygotes?\n",
    "\n",
    "Fill in the `<?>` with code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[gq20, hom_gq20, het_gq20] = vds.query_genotypes(['<?>', '<?>', '<?>'])\n",
    "print(\"Total fraction is %f, among homozygotes %f, and heterozygotes %f\" % (gq20, hom_gq20, het_gq20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}