ATGUworkshop_RVanalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rare Variant Analysis\n",
    "\n",
    "Objective:\n",
    "\n",
    "- Understand basic principles behind simple variant aggregation and burden tests.\n",
    "\n",
    "GWAS is a great tool for finding associations between **common variants** and disease, but is underpowered to detect rare-variant associations, because rare variants by definition have small sample sizes.\n",
    "\n",
    "It is possible to find associations between rare variants and disease by **grouping variants of similar effect**, and testing each group.\n",
    "\n",
    "One possible solution is to sum variant counts according to some genomic interval (for instance, gene), and then association with these intervals. A version of this kind of test is called a burden test. \n",
    "\n",
    "We'll do a burden test that associates rare variant burden with our `caffeine_consumption` phenotype. We shouldn't hope to find anything here -- especially because we've only got a few thousand rare variants!"
   ]
  },
  {
   "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",
    "# make plots display inline, rather than creating files\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=\"#1a0dab\">Step 1:</font> Import variant data\n",
    "\n",
    "First, we'll need to start again from the QC'ed matrix table on disk -- `mt` has been filtered to include only common variants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.read_matrix_table('resources/post_qc.mt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will keep variants with an allele frequency of under 1%. Including common variants will only reduce the power of a burden test."
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "# <font color=\"#1a0dab\">Step 2:</font> Group by gene\n",
    "\n",
    "\n",
    "To assign variants to genes, we'll use a tab-separated file that contains genomic intervals and corresponding genes.\n",
    "\n",
    "Additionally, you can also use our vep annotation tool which works like magic with the correct Google Cloud Platform (GCP) settings. More information [here](https://hail.is/docs/0.2/annotation_database_ui.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_ht = hl.import_table('resources/ensembl_gene_annotations.txt', impute=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_ht.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many intervals (genes) are there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_ht.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotate variants with genes\n",
    "\n",
    "In order join our two tables, we need to create a field of type `interval` so that Hail knows how to execute a join.\n",
    "\n",
    "We'll use the [transmute](https://hail.is/docs/0.2/hail.Table.html?highlight=transmute#hail.Table.transmute) function, which is like `annotate`, but drops any fields referenced in the computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('before transmute')\n",
    "gene_ht.describe()\n",
    "\n",
    "gene_ht = gene_ht.transmute(\n",
    "    interval = hl.locus_interval(gene_ht.chromosome,\n",
    "                                 gene_ht.start,\n",
    "                                 gene_ht.end))\n",
    "\n",
    "print('')\n",
    "print('after transmute')\n",
    "gene_ht.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This field needs to be the key of the table, so we will use [key_by](https://hail.is/docs/0.2/hail.Table.html?highlight=key_by#hail.Table.key_by) to assign this computed field as the table key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keyed_gene_table = gene_ht.key_by('interval')\n",
    "\n",
    "keyed_gene_table.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall how we annotated sample phenotypes earlier in the common variant tutorial -- this join looks very similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = mt.annotate_rows(gene = keyed_gene_table[mt.locus].gene_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's `show` the resulting annotations on the matrix table. How do they differ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.gene.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=\"#1a0dab\">Step 3:</font> Aggregate by gene\n",
    "\n",
    "Hail's modularity makes it easy to perform non-kernel-based burden tests.\n",
    "\n",
    "We'll compose two general tools:\n",
    " - [group_rows_by](https://hail.is/docs/0.2/hail.MatrixTable.html#hail.MatrixTable.group_rows_by) / [aggregate](https://hail.is/docs/0.2/hail.GroupedMatrixTable.html#hail.GroupedMatrixTable.aggregate)\n",
    " - [hl.linear_regression_rows](https://hail.is/docs/0.2/methods/stats.html#hail.methods.linear_regression_rows).\n",
    " \n",
    "This means that you can flexibly specify the way genotypes are summarized per gene. Using other tools, you may have a few ways to aggregate, but if you want to do something different you are out of luck!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burden_mt = (\n",
    "    mt\n",
    "    .group_rows_by('gene')\n",
    "    .aggregate(n_variants = hl.agg.count_where(mt.GT.n_alt_alleles() > 0))\n",
    ")\n",
    "\n",
    "# filter to genes with at least one rare variant!\n",
    "burden_mt = burden_mt.filter_rows(hl.agg.sum(burden_mt.n_variants) > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burden_mt.describe(widget=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burden_mt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=\"#1a0dab\">Step 4:</font> Run linear regression per gene\n",
    "\n",
    "This should look familiar! We can reuse the same modular components (like `linear_regression_rows`) for many different purposes."
   ]
  },
  {
   "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": [
    "burden_mt = burden_mt.annotate_cols(pca = pca_scores[burden_mt.s])\n",
    "\n",
    "burden_results = hl.linear_regression_rows(\n",
    "    y=burden_mt.pheno.caffeine_consumption, \n",
    "    x=burden_mt.n_variants,\n",
    "    covariates=[1.0, \n",
    "                burden_mt.pheno.is_female, \n",
    "                burden_mt.pca.scores[0], \n",
    "                burden_mt.pca.scores[1], \n",
    "                burden_mt.pca.scores[2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sorry, no `hl.plot.manhattan` for genes!\n",
    "\n",
    "Manhattan plots are really only useful for standard GWAS. Instead, we can simply sort by p-value using [order_by](https://hail.is/docs/0.2/hail.Table.html#hail.Table.order_by), and print:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burden_results.order_by(burden_results.p_value).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we use a QQ plot to help us with what we are expecting from our data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.qq(burden_results.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With relatively few points, it'll be a little unstable.\n",
    "\n",
    "RVAS QQ plots tend to be a bit lower for the same sample size.\n",
    "\n",
    "Deflation would imply an underpowered study and and this RVAS is definitely underpowered."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Any questions, team?\n",
    "\n",
    "\n",
    "### What other covariates can you think off that could possibly clean up this analysis? It's the same dataset that we played with a few weeks ago\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What do you have to do?\n",
    "\n",
    "1) Introduce yourselves! \n",
    "\n",
    "2) Identify a note-taker (and a back up, just in case). This person will also share their screen with the group for code reviewing.\n",
    "\n",
    "3) Identify a reporter who will share your group’s responses with the larger group.\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) Which genes are ranked highest? What do you think of the results? \n",
    "\n",
    "&emsp; Kumar and Arcturus will pop in and out of your rooms to check in; please use the “Ask for Help” button to bring Kumar or Arcturus into your group as and when needed "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# If you have questions, ask them! We may have answers :)"
   ]
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}