ATGUworkshop_CVanalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATGU Workshop -- Common Variant Analysis\n",
    "\n",
    "In this section, we are learning how to \n",
    "\n",
    "1) Use simple python code  and Jupyter notebook\n",
    "\n",
    "2) Setting Hail up in your environment\n",
    "\n",
    "3) Simple descriptive analysis with Hail\n",
    "\n",
    "4) Common variant analysis with Hail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hail on Jupyter\n",
    "\n",
    "Part of what we think is so exciting about Hail is that it has coincided with a larger shift in the data science community.\n",
    "\n",
    "Three years ago, most computational biologists at Broad analyzed genetic data using command-line tools, and took advantage of research compute clusters by explicitly using scheduling frameworks like LSF or Sun Grid Engine.\n",
    "\n",
    "Now, they have the option to use Hail in interactive Python notebooks backed by thousands of cores on public compute clouds like [Google Cloud](https://cloud.google.com/), [Amazon Web Services](https://aws.amazon.com/), or [Microsoft Azure](https://azure.microsoft.com/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Jupyter\n",
    "### Running cells\n",
    "Evaluate cells using SHIFT + ENTER. Select the next cell and run it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Hello, world')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modes\n",
    "\n",
    "Jupyter has two modes, a **navigation mode** and an **editor mode**.\n",
    "\n",
    "#### Navigation mode:\n",
    "\n",
    " - <font color=\"blue\"><strong>BLUE</strong></font> cell borders\n",
    " - `UP` / `DOWN` move between cells\n",
    " - `ENTER` while a cell is selected will move to **editing mode**.\n",
    " - Many letters are keyboard shortcuts! This is a common trap.\n",
    " \n",
    "#### Editor mode:\n",
    "\n",
    " - <font color=\"green\"><strong>GREEN</strong></font> cell borders\n",
    " - `UP` / `DOWN`/ move within cells before moving between cells.\n",
    " - `ESC` will return to **navigation mode**.\n",
    " - `SHIFT + ENTER` will evaluate a cell and return to **navigation mode**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell types\n",
    "\n",
    "There are several types of cells in Jupyter notebooks. The two you will see here are **Markdown** (text) and **Code**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a code cell\n",
    "my_variable = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This is a markdown cell**, so even if something looks like code (as below), it won't get executed!\n",
    "\n",
    "my_variable += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tips and tricks\n",
    "\n",
    "Keyboard shortcuts:\n",
    "\n",
    " - `SHIFT + ENTER` to evaluate a cell\n",
    " - `ESC` to return to navigation mode\n",
    " - `y` to turn a markdown cell into code\n",
    " - `m` to turn a code cell into markdown\n",
    " - `a` to add a new cell **above** the currently selected cell\n",
    " - `b` to add a new cell **below** the currently selected cell\n",
    " - `d, d` (repeated) to delete the currently selected cell\n",
    " - `TAB` to activate code completion\n",
    " \n",
    "To try this out, create a new cell below this one using `b`, and print `my_variable` by starting with `print(my` and pressing `TAB`!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up our Python environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to Hail, we import a few methods from the Hail plotting library. We'll see examples soon!\n",
    "For further installation instructions, please visit our Getting Started page (https://hail.is/docs/0.2/getting_started.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hail as hl\n",
    "from hail.plot import output_notebook, show"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we initialize Hail and set up plotting to display inline in the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hl.init()\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The demonstration materials are designed to work on a small (~20MB) downsampled chunk of the public 1000 Genomes dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to call command-line utilities from Jupyter by prefixing a line with a `!`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls -1 resources/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explore genetic data with Hail\n",
    "\n",
    "#### Learning Objectives:\n",
    "\n",
    "- To be comfortable exploring Hail data structures.\n",
    "- To understand categories of functionality for performing QC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data from VCF\n",
    "\n",
    "The [Variant Call Format (VCF)](https://en.wikipedia.org/wiki/Variant_Call_Format) is a common file format for representing genetic data collected on multiple individuals (samples).\n",
    "\n",
    "Hail's [import_vcf](https://hail.is/docs/0.2/methods/impex.html#hail.methods.import_vcf) function can read this format.\n",
    "\n",
    "However, VCF is a text format that is easy for humans to read, but very inefficient to process on a computer. \n",
    "\n",
    "The first thing we do is import (`import_vcf`) and convert the `VCF` file into a Hail native file format. This is done by using the `write` method below. The resulting file is **much** faster to process because it is scalable and easily parallelizable.\n",
    "\n",
    "Let's read in a chunk of data, specifically from 1000 Genomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hl.import_vcf('resources/1kg.vcf.bgz', min_partitions=4).write('resources/1kg.mt', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read 1KG into Hail\n",
    "\n",
    "We represent genetic data as a Hail [`MatrixTable`](https://hail.is/docs/0.2/overview/matrix_table.html), and name our variable `mt` to indicate this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.read_matrix_table('resources/1kg.mt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is a `MatrixTable`?\n",
    "\n",
    "Let's explore it!\n",
    "\n",
    "You can see:\n",
    " - **numeric** types:\n",
    "     - integers (`int32`, `int64`), e.g. `5`\n",
    "     - floating point numbers (`float32`, `float64`), e.g. `5.5` or `3e-8`\n",
    " - **strings** (`str`), e.g. `\"Foo\"`\n",
    " - **boolean** values  (`bool`) e.g. `True`\n",
    " - **collections**:\n",
    "     - arrays (`array`), e.g. `[1,1,2,3]`\n",
    "     - sets (`set`), e.g. `{1,3}`\n",
    "     - dictionaries (`dict`), e.g. `{'Foo': 5, 'Bar': 10}`\n",
    " - **genetic data types**:\n",
    "     - loci (`locus`), e.g. `[GRCh37] 1:10000` or `[GRCh38] chr1:10024`\n",
    "     - genotype calls (`call`), e.g. `0/2` or `1|0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `show`\n",
    "\n",
    "You can show individual fields like the sample ID (`s`), "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.s.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the locus (`locus`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.locus.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or the called genotype (`GT`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.GT.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `summarize`\n",
    "`summarize` Prints (potentially) useful information about any field or object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.DP.summarize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.AD.summarize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `count`\n",
    "\n",
    "`MatrixTable.count` returns a tuple with the number of rows (variants) and number of columns (samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hail has functions built for genetics\n",
    "\n",
    "For example, `hl.summarize_variants` prints useful statistics about the genetic variants in the dataset. These are not part of the generic `summarize()` function, which must support all kinds of data, not just variant data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hl.summarize_variants(mt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Annotation and quality control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrate sample information\n",
    "\n",
    "This is a text file containing phenotype information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! head resources/1kg_annotations.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can import it as a [Hail Table](https://hail.is/docs/0.2/overview/table.html) with [hl.import_table](https://hail.is/docs/0.2/methods/impex.html?highlight=import_table#hail.methods.import_table).\n",
    "\n",
    "We call it \"sa\" for \"sample annotations\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sa = hl.import_table('resources/1kg_annotations.txt', \n",
    "                      impute=True, \n",
    "                      key='s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we can see the names and types of fields in the logging messages, we can also `show` this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sa.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can `summarize` each field in `sa`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sa.summarize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add sample metadata into our 1KG `MatrixTable`\n",
    "\n",
    "It just takes one line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.annotate_cols(pheno = sa[mt.s])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What's going on here?\n",
    "\n",
    "Understanding what's going on here is a bit more difficult. To understand, we need to understand a few pieces:\n",
    "\n",
    "#### 1. `annotate` methods\n",
    "\n",
    "In Hail, `annotate` methods refer to **adding new fields**. \n",
    "\n",
    " - `MatrixTable`'s `annotate_cols` adds new column (**sample**) fields.\n",
    " - `MatrixTable`'s `annotate_rows` adds new row (**variant**) fields.\n",
    " - `MatrixTable`'s `annotate_entries` adds new entry (**genotype**) fields.\n",
    " - `Table`'s `annotate` adds new row fields.\n",
    "\n",
    "In the above cell, we are adding a new column (**sample**) field called \"pheno\". This field should be the values in our table `sa` associated with the sample ID `s` in our `MatrixTable` - that is, this is performing a **join**.\n",
    "\n",
    "Python uses square brackets to look up values in dictionaries:\n",
    "\n",
    "    d = {'foo': 5, 'bar': 10}\n",
    "    d['foo']\n",
    "\n",
    "You should think of this in much the same way - for each column of `mt`, we are looking up the fields in `sa` using the sample ID `s`.\n",
    "\n",
    "Let's look at where does this go into the `MatrixTable`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Query the phenotype fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What’s the fraction of samples with `purple_hair`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.aggregate_cols(hl.agg.fraction(mt.pheno.purple_hair))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many people are in each self-reported major ancestry group?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.aggregate_cols(hl.agg.counter(mt.pheno.super_population))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample QC\n",
    "\n",
    "We'll start with examples of sample QC.\n",
    "\n",
    "Hail has the function [hl.sample_qc](https://hail.is/docs/0.2/methods/genetics.html#hail.methods.sample_qc) to compute a list of useful statistics about samples from sequencing data. This function adds a new column field, `sample_qc`, with the computed statistics.\n",
    "\n",
    "**Click the link** above to see the documentation, which lists the fields and their descriptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.sample_qc(mt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hail includes a plotting library built on [bokeh](https://bokeh.pydata.org/en/latest/index.html) that makes it easy to visualize fields of Hail tables and matrix tables.\n",
    "\n",
    "Let's visualize the distribution of `Mean DP` (`DP` = Read Depth) to `Call Rate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.scatter(x=mt.sample_qc.dp_stats.mean,\n",
    "                    y=mt.sample_qc.call_rate,\n",
    "                    xlabel='Mean DP',\n",
    "                    ylabel='Call Rate',\n",
    "                    hover_fields={'ID': mt.s},\n",
    "                    size=8)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter columns using generated QC statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before filtering samples, we should compute a raw sample count:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.count_cols()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`filter_cols` removes entire columns from the matrix table. Here, we keep columns (samples) where the `call_rate` is over 95%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.filter_cols(mt.sample_qc.call_rate >= 0.95)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute a final sample count:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.count_cols()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_hist = mt.aggregate_entries(hl.expr.aggregators.hist(mt.DP, 0, 30, 30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del dp_hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.histogram(mt.DP, legend='DP', title='DP Histogram')\n",
    "show(p)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variant QC\n",
    "\n",
    "Hail has the function [hl.variant_qc](https://hail.is/docs/0.2/methods/genetics.html#hail.methods.variant_qc) to compute a list of useful statistics about **variants** from sequencing data.\n",
    "\n",
    "Once again, **Click the link** above to see the documentation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.variant_qc(mt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can `show()` the computed information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.variant_qc.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metrics like `call_rate` are important for QC. Let's plot the cumulative density function of call rate per variant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.filter_rows(hl.agg.call_stats(mt.GT, mt.alleles).AF[1] > 0.01).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "show(hl.plot.cdf(mt.variant_qc.call_rate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before filtering variants, we should compute a raw variant count:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pre-qc variant count\n",
    "mt.count_rows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`filter_rows` removes entire rows of the matrix table. Here, we keep rows where the `call_rate` is over 95%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.filter_rows(mt.variant_qc.call_rate > 0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After filtering, we can see more resolution of the top end of the call rate distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "show(hl.plot.cdf(mt.variant_qc.call_rate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then compute the final sample and variant count:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Association Testing and PCA\n",
    "\n",
    "We will first filter to common variants (those with an alternate allele frequency over 1%). GWAS cannot detect signal from extremely rare variants, like those only observed in one individual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.filter_rows(hl.agg.call_stats(mt.GT, mt.alleles).AF[1] > 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that in a GWAS we're running independent association tests on each variant.\n",
    "\n",
    "In Hail, the method we use is [hl.linear_regression_rows](https://hail.is/docs/0.2/methods/stats.html#hail.methods.linear_regression_rows). Why isn't this called `hl.gwas`? Modularity!\n",
    "\n",
    "We use the phenotype `caffeine_consumption` as our dependent variable, the number of alternate alleles as our independent variable, and no covariates besides an intercept term (that's the `1.0`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas = hl.linear_regression_rows(y=mt.pheno.caffeine_consumption, \n",
    "                                 x=mt.GT.n_alt_alleles(), \n",
    "                                 covariates=[1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas.describe(widget=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "Let’s visualize our association test results from the linear regression. We can do so by creating 2 common plots: a [Manhattan plot](https://en.wikipedia.org/wiki/Manhattan_plot) and a [Q-Q plot](https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot).\n",
    "\n",
    "We'll start with the Manhattan plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.manhattan(gwas.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other common plot is the Q-Q (quantile-quantile) plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.qq(gwas.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confounded!\n",
    "\n",
    "The Q-Q plot indicates **extreme** inflation of p-values.\n",
    "\n",
    "What do you think is a good range for genomic inflation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've done a GWAS before, you've probably included a few other covariates that might be confounders -- age, sex, and principal components.\n",
    "\n",
    "Principal components are a measure of genetic ancestry, and can be used to control for [population stratification](https://en.wikipedia.org/wiki/Population_stratification).\n",
    "\n",
    "Principal component analysis (PCA) is a very general statistical method for reducing high dimensional data to a small number of dimensions which capture most of the variation in the data. Hail has the function [pca](https://hail.is/docs/0.2/methods/stats.html#hail.methods.pca) for performing generic PCA.\n",
    "\n",
    "PCA typically works best on normalized data (e.g. mean centered). Hail provides the specialized function [hwe_normalized_pca](https://hail.is/docs/0.2/methods/genetics.html#hail.methods.hwe_normalized_pca) which first normalizes the genotypes according to the Hardy-Weinberg Equilibium model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_eigenvalues, pca_scores, pca_loadings = hl.hwe_normalized_pca(mt.GT, compute_loadings=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_scores.write('resources/pca_scores.ht', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pca function returns three things:\n",
    "* **eigenvalues**: an array of doubles\n",
    "* **scores**: a table keyed by sample\n",
    "* **loadings**: a table keyed by variant\n",
    "\n",
    "The **loadings** are the *principal directions*, each of which is a weighted sum of variants (a \"virtual variant\"). By default, `hwe_normalized_pca` gives us 10 principal directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.annotate_cols(pca = pca_scores[mt.s])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Control confounders and run another GWAS\n",
    "\n",
    "Now that we have computed principal components and saved it into our `mt`, let’s attempt to adjust the inflation that we saw in our initial Q-Q plot e.g.  `mt.pheno.is_female`. We will now add PCs as `covariates` in `linear_regression_rows`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas = hl.linear_regression_rows(\n",
    "    y=mt.pheno.caffeine_consumption, \n",
    "    x=mt.GT.n_alt_alleles(),\n",
    "    covariates=[1.0, mt.pca.scores[0], mt.pca.scores[1], mt.pca.scores[2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.qq(gwas.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above Q-Q plot indicates much better (but not good!) genomic control. Let's try the Manhattan plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.manhattan(gwas.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What other covariates can you think off that could possibly clean up this analysis?\n",
    "\n",
    "#### Zoom Breakout rooms Activity\n",
    "\n",
    "We will assign you into TWO breakout rooms. \n",
    "\n",
    "**Team/Room _Purple Hair_**\n",
    "\n",
    "Create a model with **purple hair** as the outcome\n",
    "\n",
    "\n",
    "**Team/Room _Polydactylism_**\n",
    "\n",
    "Create a model with **six toes** as the outcome\n",
    "\n",
    "Your assignment would be to :\n",
    "\n",
    "1) What is the distribution of people who have the phenotype? A simple list with do from `count()` or `show()`! \n",
    "\n",
    "2) Create a logistic model with the given phenotype outcome using [Hail documentation](https://hail.is/docs/0.2/methods/stats.html#hail.methods.logistic_regression_rows). Use the search function at the top of the documentation page if you need more information!  \n",
    "\n",
    "3) Print QQ plot\n",
    "\n",
    "4) Print Manhattan plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write QC'ed final dataset to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.write('output/post_qc.mt', overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "242px",
    "width": "283px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}