GWAS Phase 1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hail's overview \n",
    "\n",
    "### https://hail.is/docs/devel/tutorials-landing.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PHASE 1 - 1k Genomes Project \n",
    "\n",
    "VCFs can be find in the public google bucket : \n",
    "https://console.cloud.google.com/storage/browser/genomics-public-data/1000-genomes/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hail environment and Python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/opt/conda/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "Running on Apache Spark version 2.2.1\n",
      "SparkUI available at http://10.142.0.5:4040\n",
      "Welcome to\n",
      "     __  __     <>__\n",
      "    / /_/ /__  __/ /\n",
      "   / __  / _ `/ / /\n",
      "  /_/ /_/\\_,_/_/_/   version devel-98f4d6a179b4\n",
      "NOTE: This is a beta version. Interfaces may change\n",
      "  during the beta period. We recommend pulling\n",
      "  the latest changes weekly.\n"
     ]
    }
   ],
   "source": [
    "import hail as hl\n",
    "import hail.expr.aggregators as agg\n",
    "hl.init()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <div class=\"bk-root\">\n",
       "        <a href=\"https://bokeh.pydata.org\" target=\"_blank\" class=\"bk-logo bk-logo-small bk-logo-notebook\"></a>\n",
       "        <span id=\"73f95e6b-ae73-42d7-82a0-150b14d74371\">Loading BokehJS ...</span>\n",
       "    </div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "(function(root) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  var force = true;\n",
       "\n",
       "  if (typeof (root._bokeh_onload_callbacks) === \"undefined\" || force === true) {\n",
       "    root._bokeh_onload_callbacks = [];\n",
       "    root._bokeh_is_loading = undefined;\n",
       "  }\n",
       "\n",
       "  var JS_MIME_TYPE = 'application/javascript';\n",
       "  var HTML_MIME_TYPE = 'text/html';\n",
       "  var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
       "  var CLASS_NAME = 'output_bokeh rendered_html';\n",
       "\n",
       "  /**\n",
       "   * Render data to the DOM node\n",
       "   */\n",
       "  function render(props, node) {\n",
       "    var script = document.createElement(\"script\");\n",
       "    node.appendChild(script);\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when an output is cleared or removed\n",
       "   */\n",
       "  function handleClearOutput(event, handle) {\n",
       "    var cell = handle.cell;\n",
       "\n",
       "    var id = cell.output_area._bokeh_element_id;\n",
       "    var server_id = cell.output_area._bokeh_server_id;\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null && id in Bokeh.index) {\n",
       "      Bokeh.index[id].model.document.clear();\n",
       "      delete Bokeh.index[id];\n",
       "    }\n",
       "\n",
       "    if (server_id !== undefined) {\n",
       "      // Clean up Bokeh references\n",
       "      var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
       "      cell.notebook.kernel.execute(cmd, {\n",
       "        iopub: {\n",
       "          output: function(msg) {\n",
       "            var id = msg.content.text.trim();\n",
       "            if (id in Bokeh.index) {\n",
       "              Bokeh.index[id].model.document.clear();\n",
       "              delete Bokeh.index[id];\n",
       "            }\n",
       "          }\n",
       "        }\n",
       "      });\n",
       "      // Destroy server and session\n",
       "      var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
       "      cell.notebook.kernel.execute(cmd);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when a new output is added\n",
       "   */\n",
       "  function handleAddOutput(event, handle) {\n",
       "    var output_area = handle.output_area;\n",
       "    var output = handle.output;\n",
       "\n",
       "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
       "    if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
       "      return\n",
       "    }\n",
       "\n",
       "    var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
       "\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
       "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
       "      // store reference to embed id on output_area\n",
       "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
       "    }\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
       "      var bk_div = document.createElement(\"div\");\n",
       "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
       "      var script_attrs = bk_div.children[0].attributes;\n",
       "      for (var i = 0; i < script_attrs.length; i++) {\n",
       "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
       "      }\n",
       "      // store reference to server id on output_area\n",
       "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function register_renderer(events, OutputArea) {\n",
       "\n",
       "    function append_mime(data, metadata, element) {\n",
       "      // create a DOM node to render to\n",
       "      var toinsert = this.create_output_subarea(\n",
       "        metadata,\n",
       "        CLASS_NAME,\n",
       "        EXEC_MIME_TYPE\n",
       "      );\n",
       "      this.keyboard_manager.register_events(toinsert);\n",
       "      // Render to node\n",
       "      var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
       "      render(props, toinsert[toinsert.length - 1]);\n",
       "      element.append(toinsert);\n",
       "      return toinsert\n",
       "    }\n",
       "\n",
       "    /* Handle when an output is cleared or removed */\n",
       "    events.on('clear_output.CodeCell', handleClearOutput);\n",
       "    events.on('delete.Cell', handleClearOutput);\n",
       "\n",
       "    /* Handle when a new output is added */\n",
       "    events.on('output_added.OutputArea', handleAddOutput);\n",
       "\n",
       "    /**\n",
       "     * Register the mime type and append_mime function with output_area\n",
       "     */\n",
       "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
       "      /* Is output safe? */\n",
       "      safe: true,\n",
       "      /* Index of renderer in `output_area.display_order` */\n",
       "      index: 0\n",
       "    });\n",
       "  }\n",
       "\n",
       "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
       "  if (root.Jupyter !== undefined) {\n",
       "    var events = require('base/js/events');\n",
       "    var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
       "\n",
       "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
       "      register_renderer(events, OutputArea);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  \n",
       "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
       "    root._bokeh_timeout = Date.now() + 5000;\n",
       "    root._bokeh_failed_load = false;\n",
       "  }\n",
       "\n",
       "  var NB_LOAD_WARNING = {'data': {'text/html':\n",
       "     \"<div style='background-color: #fdd'>\\n\"+\n",
       "     \"<p>\\n\"+\n",
       "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
       "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
       "     \"</p>\\n\"+\n",
       "     \"<ul>\\n\"+\n",
       "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
       "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
       "     \"</ul>\\n\"+\n",
       "     \"<code>\\n\"+\n",
       "     \"from bokeh.resources import INLINE\\n\"+\n",
       "     \"output_notebook(resources=INLINE)\\n\"+\n",
       "     \"</code>\\n\"+\n",
       "     \"</div>\"}};\n",
       "\n",
       "  function display_loaded() {\n",
       "    var el = document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\");\n",
       "    if (el != null) {\n",
       "      el.textContent = \"BokehJS is loading...\";\n",
       "    }\n",
       "    if (root.Bokeh !== undefined) {\n",
       "      if (el != null) {\n",
       "        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
       "      }\n",
       "    } else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(display_loaded, 100)\n",
       "    }\n",
       "  }\n",
       "\n",
       "\n",
       "  function run_callbacks() {\n",
       "    try {\n",
       "      root._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n",
       "    }\n",
       "    finally {\n",
       "      delete root._bokeh_onload_callbacks\n",
       "    }\n",
       "    console.info(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(js_urls, callback) {\n",
       "    root._bokeh_onload_callbacks.push(callback);\n",
       "    if (root._bokeh_is_loading > 0) {\n",
       "      console.log(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    root._bokeh_is_loading = js_urls.length;\n",
       "    for (var i = 0; i < js_urls.length; i++) {\n",
       "      var url = js_urls[i];\n",
       "      var s = document.createElement('script');\n",
       "      s.src = url;\n",
       "      s.async = false;\n",
       "      s.onreadystatechange = s.onload = function() {\n",
       "        root._bokeh_is_loading--;\n",
       "        if (root._bokeh_is_loading === 0) {\n",
       "          console.log(\"Bokeh: all BokehJS libraries loaded\");\n",
       "          run_callbacks()\n",
       "        }\n",
       "      };\n",
       "      s.onerror = function() {\n",
       "        console.warn(\"failed to load library \" + url);\n",
       "      };\n",
       "      console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "    }\n",
       "  };var element = document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\");\n",
       "  if (element == null) {\n",
       "    console.log(\"Bokeh: ERROR: autoload.js configured with elementid '73f95e6b-ae73-42d7-82a0-150b14d74371' but no matching script tag was found. \")\n",
       "    return false;\n",
       "  }\n",
       "\n",
       "  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-0.13.0.min.js\"];\n",
       "\n",
       "  var inline_js = [\n",
       "    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "    \n",
       "    function(Bokeh) {\n",
       "      \n",
       "    },\n",
       "    function(Bokeh) {\n",
       "      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.css\");\n",
       "      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.css\");\n",
       "      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.css\");\n",
       "      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.css\");\n",
       "      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.css\");\n",
       "      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.css\");\n",
       "    }\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    \n",
       "    if ((root.Bokeh !== undefined) || (force === true)) {\n",
       "      for (var i = 0; i < inline_js.length; i++) {\n",
       "        inline_js[i].call(root, root.Bokeh);\n",
       "      }if (force === true) {\n",
       "        display_loaded();\n",
       "      }} else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(run_inline_js, 100);\n",
       "    } else if (!root._bokeh_failed_load) {\n",
       "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
       "      root._bokeh_failed_load = true;\n",
       "    } else if (force !== true) {\n",
       "      var cell = $(document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\")).parents('.cell').data().cell;\n",
       "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  if (root._bokeh_is_loading === 0) {\n",
       "    console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(js_urls, function() {\n",
       "      console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(window));"
      ],
      "application/vnd.bokehjs_load.v0+json": "\n(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  var force = true;\n\n  if (typeof (root._bokeh_onload_callbacks) === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n  \n\n  \n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  var NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded() {\n    var el = document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\");\n    if (el != null) {\n      el.textContent = \"BokehJS is loading...\";\n    }\n    if (root.Bokeh !== undefined) {\n      if (el != null) {\n        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(display_loaded, 100)\n    }\n  }\n\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n    }\n    finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.info(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(js_urls, callback) {\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.log(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = js_urls.length;\n    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      var s = document.createElement('script');\n      s.src = url;\n      s.async = false;\n      s.onreadystatechange = s.onload = function() {\n        root._bokeh_is_loading--;\n        if (root._bokeh_is_loading === 0) {\n          console.log(\"Bokeh: all BokehJS libraries loaded\");\n          run_callbacks()\n        }\n      };\n      s.onerror = function() {\n        console.warn(\"failed to load library \" + url);\n      };\n      console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.getElementsByTagName(\"head\")[0].appendChild(s);\n    }\n  };var element = document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\");\n  if (element == null) {\n    console.log(\"Bokeh: ERROR: autoload.js configured with elementid '73f95e6b-ae73-42d7-82a0-150b14d74371' but no matching script tag was found. \")\n    return false;\n  }\n\n  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-0.13.0.min.js\"];\n\n  var inline_js = [\n    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\n    \n    function(Bokeh) {\n      \n    },\n    function(Bokeh) {\n      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.css\");\n      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.13.0.min.css\");\n      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.css\");\n      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.13.0.min.css\");\n      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.css\");\n      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-tables-0.13.0.min.css\");\n    }\n  ];\n\n  function run_inline_js() {\n    \n    if ((root.Bokeh !== undefined) || (force === true)) {\n      for (var i = 0; i < inline_js.length; i++) {\n        inline_js[i].call(root, root.Bokeh);\n      }if (force === true) {\n        display_loaded();\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      var cell = $(document.getElementById(\"73f95e6b-ae73-42d7-82a0-150b14d74371\")).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(js_urls, function() {\n      console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No module named 'dask'\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from collections import Counter\n",
    "from math import log, isnan\n",
    "from pprint import pprint\n",
    "from bokeh.io import show, output_notebook\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.layouts import gridplot\n",
    "from bokeh.models import Span\n",
    "output_notebook()\n",
    "import plotly\n",
    "import time\n",
    "import gcsfs\n",
    "import csv\n",
    "from bokeh.io import output_notebook, push_notebook, show\n",
    "import seaborn as sb\n",
    "import random\n",
    "import gcsfs\n",
    "fs = gcsfs.GCSFileSystem(project='avl-hail-ines')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the Matrix Table Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MT size:  (39706715, 1092)\n"
     ]
    }
   ],
   "source": [
    "# Phase 1 - 1k Genome Project\n",
    "mt = hl.read_matrix_table(\"gs://1k-genome/1000-genomes/VDS-of-all/ALL.chr.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.mt\")\n",
    "print('MT size: ', mt.count())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `cache` is used to optimize some of the downstream operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.variant_qc(mt).cache()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------\n",
      "Global fields:\n",
      "    None\n",
      "----------------------------------------\n",
      "Column fields:\n",
      "    's': str \n",
      "----------------------------------------\n",
      "Row fields:\n",
      "    'locus': locus<GRCh37> \n",
      "    'alleles': array<str> \n",
      "    'rsid': str \n",
      "    'qual': float64 \n",
      "    'filters': set<str> \n",
      "    'info': struct {\n",
      "        LDAF: float64, \n",
      "        AVGPOST: float64, \n",
      "        RSQ: float64, \n",
      "        ERATE: float64, \n",
      "        THETA: float64, \n",
      "        CIEND: array<int32>, \n",
      "        CIPOS: array<int32>, \n",
      "        END: int32, \n",
      "        HOMLEN: array<int32>, \n",
      "        HOMSEQ: array<str>, \n",
      "        SVLEN: int32, \n",
      "        SVTYPE: str, \n",
      "        AC: array<int32>, \n",
      "        AN: int32, \n",
      "        AA: str, \n",
      "        AF: array<float64>, \n",
      "        AMR_AF: float64, \n",
      "        ASN_AF: float64, \n",
      "        AFR_AF: float64, \n",
      "        EUR_AF: float64, \n",
      "        VT: str, \n",
      "        SNPSOURCE: array<str>\n",
      "    } \n",
      "    'variant_qc': struct {\n",
      "        AC: array<int32>, \n",
      "        AF: array<float64>, \n",
      "        AN: int32, \n",
      "        homozygote_count: array<int32>, \n",
      "        n_called: int64, \n",
      "        n_not_called: int64, \n",
      "        call_rate: float32, \n",
      "        n_het: int64, \n",
      "        n_non_ref: int64, \n",
      "        het_freq_hwe: float64, \n",
      "        p_value_hwe: float64\n",
      "    } \n",
      "----------------------------------------\n",
      "Entry fields:\n",
      "    'GT': call \n",
      "    'DS': float64 \n",
      "    'GL': array<float64> \n",
      "    'PL': array<int32> \n",
      "----------------------------------------\n",
      "Column key: ['s']\n",
      "Row key: ['locus', 'alleles']\n",
      "Partition key: ['locus']\n",
      "----------------------------------------\n"
     ]
    }
   ],
   "source": [
    "mt.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================\n",
      "Number of variants: 39706715\n",
      "==============================\n",
      "Alleles per variant\n",
      "-------------------\n",
      "  2 alleles: 39706715 variants\n",
      "==============================\n",
      "Variants per contig\n",
      "-------------------\n",
      "   1: 3007196 variants\n",
      "   2: 3307592 variants\n",
      "   3: 2763454 variants\n",
      "   4: 2736765 variants\n",
      "   5: 2530217 variants\n",
      "   6: 2424425 variants\n",
      "   7: 2215231 variants\n",
      "   8: 2183839 variants\n",
      "   9: 1652388 variants\n",
      "  10: 1882663 variants\n",
      "  11: 1894908 variants\n",
      "  12: 1828006 variants\n",
      "  13: 1373000 variants\n",
      "  14: 1258254 variants\n",
      "  15: 1130554 variants\n",
      "  16: 1210619 variants\n",
      "  17: 1046733 variants\n",
      "  18: 1088820 variants\n",
      "  19: 816115 variants\n",
      "  20: 855166 variants\n",
      "  21: 518965 variants\n",
      "  22: 494328 variants\n",
      "   X: 1487477 variants\n",
      "==============================\n",
      "Allele type distribution\n",
      "------------------------\n",
      "        SNP: 38248779 alternate alleles\n",
      "   Deletion: 872267 alternate alleles\n",
      "  Insertion: 577895 alternate alleles\n",
      "   Symbolic: 5488 alternate alleles\n",
      "    Complex: 2286 alternate alleles\n",
      "==============================\n"
     ]
    }
   ],
   "source": [
    "hl.summarize_variants(mt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Struct(GT=0|0, DS=0.2, GL=[-0.18, -0.47, -2.42], PL=None),\n",
       " Struct(GT=0|0, DS=0.15, GL=[-0.24, -0.44, -1.16], PL=None),\n",
       " Struct(GT=0|0, DS=0.15, GL=[-0.15, -0.54, -3.12], PL=None),\n",
       " Struct(GT=0|1, DS=0.6, GL=[-0.48, -0.48, -0.48], PL=None),\n",
       " Struct(GT=0|0, DS=0.55, GL=[-0.48, -0.48, -0.48], PL=None)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mt.entry.take(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+------------+---------+------+-------------+---------------------+\n",
      "| locus         | alleles    | s       | GT   |          DS | GL                  |\n",
      "+---------------+------------+---------+------+-------------+---------------------+\n",
      "| locus<GRCh37> | array<str> | str     | call |     float64 | array<float64>      |\n",
      "+---------------+------------+---------+------+-------------+---------------------+\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00096 | 0|0  | 2.00000e-01 | [-0.18,-0.47,-2.42] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00097 | 0|0  | 1.50000e-01 | [-0.24,-0.44,-1.16] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00099 | 0|0  | 1.50000e-01 | [-0.15,-0.54,-3.12] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00100 | 0|1  | 6.00000e-01 | [-0.48,-0.48,-0.48] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00101 | 0|0  | 5.50000e-01 | [-0.48,-0.48,-0.48] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00102 | 0|1  | 9.50000e-01 | [-1.92,-0.01,-2.5]  |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00103 | 0|0  | 5.00000e-02 | [-0.05,-0.93,-5.0]  |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00104 | 0|0  | 1.00000e-01 | [-0.11,-0.66,-4.22] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00106 | 0|1  | 5.50000e-01 | [-0.28,-0.43,-0.96] |\n",
      "| 1:10583       | [\"G\",\"A\"]  | HG00108 | 0|0  | 4.50000e-01 | [-0.48,-0.48,-0.48] |\n",
      "+---------------+------------+---------+------+-------------+---------------------+\n",
      "\n",
      "+--------------+\n",
      "| PL           |\n",
      "+--------------+\n",
      "| array<int32> |\n",
      "+--------------+\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "| null         |\n",
      "+--------------+\n",
      "showing top 10 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mt.entry.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+------------+\n",
      "| locus         | alleles    |\n",
      "+---------------+------------+\n",
      "| locus<GRCh37> | array<str> |\n",
      "+---------------+------------+\n",
      "| 1:10583       | [\"G\",\"A\"]  |\n",
      "| 1:10611       | [\"C\",\"G\"]  |\n",
      "| 1:13302       | [\"C\",\"T\"]  |\n",
      "| 1:13327       | [\"G\",\"C\"]  |\n",
      "| 1:13957       | [\"TC\",\"T\"] |\n",
      "+---------------+------------+\n",
      "showing top 5 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mt.rows().select().show(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Annotation file : encode the label Super_population "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------\n",
      "Global fields:\n",
      "    None\n",
      "----------------------------------------\n",
      "Row fields:\n",
      "    'Sample': str \n",
      "    'Family_ID': str \n",
      "    'Population': str \n",
      "    'Population_Description': str \n",
      "    'Gender': str \n",
      "    'Relationship': str \n",
      "    'Unexpected_Parent_Child': str \n",
      "    'Non_Paternity': str \n",
      "    'Siblings': str \n",
      "    'Grandparents': str \n",
      "    'Avuncular': str \n",
      "    'Half_Siblings': str \n",
      "    'Unknown_Second_Order': str \n",
      "    'Third_Order': str \n",
      "    'In_Low_Coverage_Pilot': str \n",
      "    'LC_Pilot_Platforms': str \n",
      "    'LC_Pilot_Centers': str \n",
      "    'In_High_Coverage_Pilot': str \n",
      "    'HC_Pilot_Platforms': str \n",
      "    'HC_Pilot_Centers': str \n",
      "    'In_Exon_Targetted_Pilot': str \n",
      "    'ET_Pilot_Platforms': str \n",
      "    'ET_Pilot_Centers': str \n",
      "    'Has_Sequence_in_Phase1': str \n",
      "    'Phase1_LC_Platform': str \n",
      "    'Phase1_LC_Centers': str \n",
      "    'Phase1_E_Platform': str \n",
      "    'Phase1_E_Centers': str \n",
      "    'In_Phase1_Integrated_Variant_Set': str \n",
      "    'Has_Phase1_chrY_SNPS': str \n",
      "    'Has_phase1_chrY_Deletions': str \n",
      "    'Has_phase1_chrMT_SNPs': str \n",
      "    'Main_project_LC_Centers': str \n",
      "    'Main_project_LC_platform': str \n",
      "    'Total_LC_Sequence': str \n",
      "    'LC_Non_Duplicated_Aligned_Coverage': str \n",
      "    'Main_Project_E_Centers': str \n",
      "    'Main_Project_E_Platform': str \n",
      "    'Total_Exome_Sequence': str \n",
      "    'X_Targets_Covered_to_20x_or_greater': str \n",
      "    'VerifyBam_E_Omni_Free': str \n",
      "    'VerifyBam_E_Affy_Free': str \n",
      "    'VerifyBam_E_Omni_Chip': str \n",
      "    'VerifyBam_E_Affy_Chip': str \n",
      "    'VerifyBam_LC_Omni_Free': str \n",
      "    'VerifyBam_LC_Affy_Free': str \n",
      "    'VerifyBam_LC_Omni_Chip': str \n",
      "    'VerifyBam_LC_Affy_Chip': str \n",
      "    'LC_Indel_Ratio': str \n",
      "    'E_Indel_Ratio': str \n",
      "    'LC_Passed_QC': str \n",
      "    'E_Passed_QC': str \n",
      "    'In_Final_Phase_Variant_Calling': str \n",
      "    'Has_Omni_Genotypes': str \n",
      "    'Has_Axiom_Genotypes': str \n",
      "    'Has_Affy_6_0_Genotypes': str \n",
      "    'Has_Exome_LOF_Genotypes': str \n",
      "    'EBV_Coverage': str \n",
      "    'DNA_Source_from_Coriell': str \n",
      "    'Has_Sequence_from_Blood_in_Index': str \n",
      "    'Super_Population': str \n",
      "    'Super_Population_Description': str \n",
      "----------------------------------------\n",
      "Key: ['Sample']\n",
      "----------------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 12:40:43 Hail: INFO: Reading table with no type imputation\n",
      "  Loading column 'Sample' as type 'str' (type not specified)\n",
      "  Loading column 'Family_ID' as type 'str' (type not specified)\n",
      "  Loading column 'Population' as type 'str' (type not specified)\n",
      "  Loading column 'Population_Description' as type 'str' (type not specified)\n",
      "  Loading column 'Gender' as type 'str' (type not specified)\n",
      "  Loading column 'Relationship' as type 'str' (type not specified)\n",
      "  Loading column 'Unexpected_Parent_Child' as type 'str' (type not specified)\n",
      "  Loading column 'Non_Paternity' as type 'str' (type not specified)\n",
      "  Loading column 'Siblings' as type 'str' (type not specified)\n",
      "  Loading column 'Grandparents' as type 'str' (type not specified)\n",
      "  Loading column 'Avuncular' as type 'str' (type not specified)\n",
      "  Loading column 'Half_Siblings' as type 'str' (type not specified)\n",
      "  Loading column 'Unknown_Second_Order' as type 'str' (type not specified)\n",
      "  Loading column 'Third_Order' as type 'str' (type not specified)\n",
      "  Loading column 'In_Low_Coverage_Pilot' as type 'str' (type not specified)\n",
      "  Loading column 'LC_Pilot_Platforms' as type 'str' (type not specified)\n",
      "  Loading column 'LC_Pilot_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'In_High_Coverage_Pilot' as type 'str' (type not specified)\n",
      "  Loading column 'HC_Pilot_Platforms' as type 'str' (type not specified)\n",
      "  Loading column 'HC_Pilot_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'In_Exon_Targetted_Pilot' as type 'str' (type not specified)\n",
      "  Loading column 'ET_Pilot_Platforms' as type 'str' (type not specified)\n",
      "  Loading column 'ET_Pilot_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Sequence_in_Phase1' as type 'str' (type not specified)\n",
      "  Loading column 'Phase1_LC_Platform' as type 'str' (type not specified)\n",
      "  Loading column 'Phase1_LC_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'Phase1_E_Platform' as type 'str' (type not specified)\n",
      "  Loading column 'Phase1_E_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'In_Phase1_Integrated_Variant_Set' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Phase1_chrY_SNPS' as type 'str' (type not specified)\n",
      "  Loading column 'Has_phase1_chrY_Deletions' as type 'str' (type not specified)\n",
      "  Loading column 'Has_phase1_chrMT_SNPs' as type 'str' (type not specified)\n",
      "  Loading column 'Main_project_LC_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'Main_project_LC_platform' as type 'str' (type not specified)\n",
      "  Loading column 'Total_LC_Sequence' as type 'str' (type not specified)\n",
      "  Loading column 'LC_Non_Duplicated_Aligned_Coverage' as type 'str' (type not specified)\n",
      "  Loading column 'Main_Project_E_Centers' as type 'str' (type not specified)\n",
      "  Loading column 'Main_Project_E_Platform' as type 'str' (type not specified)\n",
      "  Loading column 'Total_Exome_Sequence' as type 'str' (type not specified)\n",
      "  Loading column 'X_Targets_Covered_to_20x_or_greater' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_E_Omni_Free' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_E_Affy_Free' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_E_Omni_Chip' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_E_Affy_Chip' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_LC_Omni_Free' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_LC_Affy_Free' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_LC_Omni_Chip' as type 'str' (type not specified)\n",
      "  Loading column 'VerifyBam_LC_Affy_Chip' as type 'str' (type not specified)\n",
      "  Loading column 'LC_Indel_Ratio' as type 'str' (type not specified)\n",
      "  Loading column 'E_Indel_Ratio' as type 'str' (type not specified)\n",
      "  Loading column 'LC_Passed_QC' as type 'str' (type not specified)\n",
      "  Loading column 'E_Passed_QC' as type 'str' (type not specified)\n",
      "  Loading column 'In_Final_Phase_Variant_Calling' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Omni_Genotypes' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Axiom_Genotypes' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Affy_6_0_Genotypes' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Exome_LOF_Genotypes' as type 'str' (type not specified)\n",
      "  Loading column 'EBV_Coverage' as type 'str' (type not specified)\n",
      "  Loading column 'DNA_Source_from_Coriell' as type 'str' (type not specified)\n",
      "  Loading column 'Has_Sequence_from_Blood_in_Index' as type 'str' (type not specified)\n",
      "  Loading column 'Super_Population' as type 'str' (type not specified)\n",
      "  Loading column 'Super_Population_Description' as type 'str' (type not specified)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Import Annotation File as a hail.Table\n",
    "table = (hl.import_table('gs://1k-genome/1000-genomes/other/sample_info/sample_info.csv', delimiter=',', missing='', quote='\"').key_by('Sample'))\n",
    "table.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Annotation File Size :  (3500, 62)\n"
     ]
    }
   ],
   "source": [
    "# Import Annotation File as a pandas.DataFrame\n",
    "with fs.open('1k-genome/1000-genomes/other/sample_info/sample_info.csv',\"rt\") as f:\n",
    "    data = pd.read_csv(f, na_values = str, header = 0)\n",
    "print ('Annotation File Size : ', data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.astype(str)\n",
    "data[['In_Low_Coverage_Pilot', 'In_High_Coverage_Pilot', 'In_Exon_Targetted_Pilot', 'Has_Sequence_in_Phase1', \n",
    "      'In_Phase1_Integrated_Variant_Set',\"Has_Phase1_chrY_SNPS\",\n",
    "        'Has_phase1_chrY_Deletions','Has_phase1_chrMT_SNPs', 'Total_LC_Sequence','LC_Non_Duplicated_Aligned_Coverage', 'Total_Exome_Sequence',\n",
    "        'X_Targets_Covered_to_20x_or_greater','VerifyBam_E_Omni_Free','VerifyBam_E_Affy_Free','VerifyBam_E_Omni_Chip','VerifyBam_E_Affy_Chip','VerifyBam_LC_Omni_Free',                 \n",
    "        'VerifyBam_LC_Affy_Free', 'VerifyBam_LC_Omni_Chip','VerifyBam_LC_Affy_Chip','LC_Indel_Ratio','E_Indel_Ratio',                          \n",
    "        'LC_Passed_QC', 'E_Passed_QC', 'In_Final_Phase_Variant_Calling','Has_Omni_Genotypes','Has_Axiom_Genotypes','Has_Affy_6_0_Genotypes',                 \n",
    "        'Has_Exome_LOF_Genotypes','EBV_Coverage', \n",
    "        'Has_Sequence_from_Blood_in_Index' ]] = data[['In_Low_Coverage_Pilot', 'In_High_Coverage_Pilot', 'In_Exon_Targetted_Pilot', 'Has_Sequence_in_Phase1', 'In_Phase1_Integrated_Variant_Set',\"Has_Phase1_chrY_SNPS\",\n",
    "        'Has_phase1_chrY_Deletions','Has_phase1_chrMT_SNPs', 'Total_LC_Sequence','LC_Non_Duplicated_Aligned_Coverage', 'Total_Exome_Sequence',\n",
    "        'X_Targets_Covered_to_20x_or_greater','VerifyBam_E_Omni_Free','VerifyBam_E_Affy_Free','VerifyBam_E_Omni_Chip','VerifyBam_E_Affy_Chip','VerifyBam_LC_Omni_Free',                 \n",
    "        'VerifyBam_LC_Affy_Free',                 \n",
    "        'VerifyBam_LC_Omni_Chip',                 \n",
    "        'VerifyBam_LC_Affy_Chip',                 \n",
    "        'LC_Indel_Ratio',                         \n",
    "        'E_Indel_Ratio',                          \n",
    "        'LC_Passed_QC',                           \n",
    "        'E_Passed_QC',                            \n",
    "        'In_Final_Phase_Variant_Calling',         \n",
    "        'Has_Omni_Genotypes',                     \n",
    "        'Has_Axiom_Genotypes',                    \n",
    "        'Has_Affy_6_0_Genotypes',                 \n",
    "        'Has_Exome_LOF_Genotypes',                \n",
    "        'EBV_Coverage', 'Has_Sequence_from_Blood_in_Index' ]].astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the *LabelEncoder* method  from `scikit-learn` to encode the label Super Population which contains 4 groups for phase 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import LabelEncoder\n",
    "key = data['Super_Population'].values\n",
    "\n",
    "#create Label encoder for Super_Population\n",
    "le = LabelEncoder()\n",
    "le.fit(['SAS', 'EUR', 'AFR', 'AMR', 'EAS'])\n",
    "\n",
    "#append to the Dataframe\n",
    "data[\"Super_Population_Classification\"] = pd.Series(le.transform(key), index=data.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Annotate the MT file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 12:40:44 Hail: INFO: Coerced sorted dataset\n"
     ]
    }
   ],
   "source": [
    "#Add the label Super_Population_Classification \n",
    "table1 = hl.Table.to_pandas(table)\n",
    "table1[\"Super_Population_Classification\"] = pd.Series(le.transform(key).astype(float)) \n",
    "table1 = hl.Table.from_pandas(table1).key_by('Sample') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `annotate_cols` method to join the table with the MatrixTable containing our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 12:40:53 Hail: INFO: Coerced almost-sorted dataset\n",
      "2018-08-07 12:41:17 Hail: INFO: Ordering unsorted dataset with network shuffle\n"
     ]
    }
   ],
   "source": [
    "#Annotate the MF file\n",
    "mt = mt.annotate_cols(**table1[mt.s])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query functions and the Hail Expression Language"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `aggregate` method can be used to aggregate over rows of the table.\n",
    "`counter` is an aggregation function that counts the number of occurrences of each unique element. We can use this to pull out the population distribution by passing in a Hail Expression for the field that we want to count by."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 12:41:33 Hail: INFO: Coerced almost-sorted dataset\n",
      "2018-08-07 12:41:46 Hail: INFO: Ordering unsorted dataset with network shuffle\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'AFR': 1018, 'AMR': 535, 'EAS': 617, 'EUR': 669, 'SAS': 661}\n"
     ]
    }
   ],
   "source": [
    "pprint(table1.aggregate(agg.counter(table1.Super_Population)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "snp_counts = mt.aggregate_rows(agg.counter(hl.Struct(ref=mt.alleles[0], alt=mt.alleles[1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can list the counts in descending order using Python’s Counter class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(Struct(ref=G, alt=A), 7768570),\n",
       " (Struct(ref=C, alt=T), 7730953),\n",
       " (Struct(ref=A, alt=G), 5318273),\n",
       " (Struct(ref=T, alt=C), 5312221),\n",
       " (Struct(ref=C, alt=A), 1772424),\n",
       " (Struct(ref=G, alt=T), 1758628),\n",
       " (Struct(ref=C, alt=G), 1633912),\n",
       " (Struct(ref=G, alt=C), 1624685),\n",
       " (Struct(ref=A, alt=C), 1354055),\n",
       " (Struct(ref=T, alt=G), 1350151),\n",
       " (Struct(ref=T, alt=A), 1313500),\n",
       " (Struct(ref=A, alt=T), 1311407),\n",
       " (Struct(ref=A, alt=AT), 63401),\n",
       " (Struct(ref=T, alt=TA), 58720),\n",
       " (Struct(ref=AT, alt=A), 53873),\n",
       " (Struct(ref=TA, alt=T), 49001),\n",
       " (Struct(ref=C, alt=CA), 45240),\n",
       " (Struct(ref=C, alt=CT), 41699),\n",
       " (Struct(ref=G, alt=GA), 41038),\n",
       " (Struct(ref=CT, alt=C), 39935),\n",
       " (Struct(ref=TG, alt=T), 37299),\n",
       " (Struct(ref=G, alt=GT), 36619),\n",
       " (Struct(ref=AG, alt=A), 35852),\n",
       " (Struct(ref=CA, alt=C), 35025),\n",
       " (Struct(ref=GA, alt=G), 33063),\n",
       " (Struct(ref=AC, alt=A), 32616),\n",
       " (Struct(ref=TC, alt=T), 30579),\n",
       " (Struct(ref=T, alt=TG), 25565),\n",
       " (Struct(ref=T, alt=TC), 23599),\n",
       " (Struct(ref=GT, alt=G), 22966),\n",
       " (Struct(ref=A, alt=AG), 20299),\n",
       " (Struct(ref=A, alt=AC), 19724),\n",
       " (Struct(ref=GC, alt=G), 16872),\n",
       " (Struct(ref=CAT, alt=C), 16495),\n",
       " (Struct(ref=G, alt=GC), 11477),\n",
       " (Struct(ref=CAG, alt=C), 9354),\n",
       " (Struct(ref=CTT, alt=C), 8568),\n",
       " (Struct(ref=CTG, alt=C), 7393),\n",
       " (Struct(ref=ACT, alt=A), 7369),\n",
       " (Struct(ref=TTC, alt=T), 7325),\n",
       " (Struct(ref=AAG, alt=A), 7156),\n",
       " (Struct(ref=AAT, alt=A), 6782),\n",
       " (Struct(ref=CG, alt=C), 6662),\n",
       " (Struct(ref=AAC, alt=A), 6487),\n",
       " (Struct(ref=TTA, alt=T), 6267),\n",
       " (Struct(ref=C, alt=CAT), 6188),\n",
       " (Struct(ref=TAC, alt=T), 5390),\n",
       " (Struct(ref=CCT, alt=C), 5210),\n",
       " (Struct(ref=C, alt=CG), 5134),\n",
       " (Struct(ref=TAA, alt=T), 4979),\n",
       " (Struct(ref=ATG, alt=A), 4852),\n",
       " (Struct(ref=CAA, alt=C), 4671),\n",
       " (Struct(ref=TTG, alt=T), 4513),\n",
       " (Struct(ref=AAAG, alt=A), 4378),\n",
       " (Struct(ref=GTA, alt=G), 4250),\n",
       " (Struct(ref=T, alt=TTC), 4107),\n",
       " (Struct(ref=G, alt=GTA), 4014),\n",
       " (Struct(ref=T, alt=TTA), 3825),\n",
       " (Struct(ref=CTA, alt=C), 3717),\n",
       " (Struct(ref=TCA, alt=T), 3405),\n",
       " (Struct(ref=A, alt=AAT), 3342),\n",
       " (Struct(ref=T, alt=TAC), 3100),\n",
       " (Struct(ref=A, alt=ATG), 3018),\n",
       " (Struct(ref=C, alt=CTT), 3007),\n",
       " (Struct(ref=AAAAC, alt=A), 2911),\n",
       " (Struct(ref=CCTT, alt=C), 2852),\n",
       " (Struct(ref=TTTC, alt=T), 2820),\n",
       " (Struct(ref=TGA, alt=T), 2753),\n",
       " (Struct(ref=AAAC, alt=A), 2737),\n",
       " (Struct(ref=T, alt=TAA), 2729),\n",
       " (Struct(ref=ATC, alt=A), 2725),\n",
       " (Struct(ref=TCTC, alt=T), 2718),\n",
       " (Struct(ref=CCA, alt=C), 2607),\n",
       " (Struct(ref=T, alt=TTG), 2589),\n",
       " (Struct(ref=GAA, alt=G), 2562),\n",
       " (Struct(ref=GAC, alt=G), 2560),\n",
       " (Struct(ref=C, alt=CAG), 2556),\n",
       " (Struct(ref=TTTTG, alt=T), 2533),\n",
       " (Struct(ref=AAAAG, alt=A), 2513),\n",
       " (Struct(ref=AAAT, alt=A), 2479),\n",
       " (Struct(ref=A, alt=AAC), 2477),\n",
       " (Struct(ref=ATT, alt=A), 2468),\n",
       " (Struct(ref=A, alt=ACT), 2444),\n",
       " (Struct(ref=TTTG, alt=T), 2309),\n",
       " (Struct(ref=AGT, alt=A), 2295),\n",
       " (Struct(ref=GCA, alt=G), 2280),\n",
       " (Struct(ref=TAG, alt=T), 2245),\n",
       " (Struct(ref=C, alt=CTG), 2241),\n",
       " (Struct(ref=GTC, alt=G), 2163),\n",
       " (Struct(ref=C, alt=CAA), 2127),\n",
       " (Struct(ref=TATA, alt=T), 2099),\n",
       " (Struct(ref=C, alt=CTA), 2072),\n",
       " (Struct(ref=TTTTC, alt=T), 2063),\n",
       " (Struct(ref=GAT, alt=G), 2030),\n",
       " (Struct(ref=A, alt=ATT), 2028),\n",
       " (Struct(ref=GGA, alt=G), 1993),\n",
       " (Struct(ref=ATAT, alt=A), 1943),\n",
       " (Struct(ref=A, alt=AAG), 1884),\n",
       " (Struct(ref=CTCTT, alt=C), 1835),\n",
       " (Struct(ref=TTTA, alt=T), 1756),\n",
       " (Struct(ref=CATT, alt=C), 1741),\n",
       " (Struct(ref=GGT, alt=G), 1716),\n",
       " (Struct(ref=AAAAT, alt=A), 1716),\n",
       " (Struct(ref=G, alt=GAT), 1621),\n",
       " (Struct(ref=TAAG, alt=T), 1590),\n",
       " (Struct(ref=TAAC, alt=T), 1550),\n",
       " (Struct(ref=G, alt=GAA), 1528),\n",
       " (Struct(ref=A, alt=<DEL>), 1522),\n",
       " (Struct(ref=T, alt=<DEL>), 1504),\n",
       " (Struct(ref=TATC, alt=T), 1492),\n",
       " (Struct(ref=T, alt=TCA), 1465),\n",
       " (Struct(ref=GAGA, alt=G), 1444),\n",
       " (Struct(ref=T, alt=TTTC), 1394),\n",
       " (Struct(ref=GCT, alt=G), 1382),\n",
       " (Struct(ref=ACTT, alt=A), 1378),\n",
       " (Struct(ref=T, alt=TGA), 1370),\n",
       " (Struct(ref=C, alt=CCT), 1353),\n",
       " (Struct(ref=G, alt=GTT), 1335),\n",
       " (Struct(ref=TTTTA, alt=T), 1305),\n",
       " (Struct(ref=A, alt=ATC), 1291),\n",
       " (Struct(ref=TAAAC, alt=T), 1262),\n",
       " (Struct(ref=G, alt=GTC), 1251),\n",
       " (Struct(ref=GTT, alt=G), 1248),\n",
       " (Struct(ref=G, alt=<DEL>), 1241),\n",
       " (Struct(ref=TTAA, alt=T), 1227),\n",
       " (Struct(ref=C, alt=<DEL>), 1221),\n",
       " (Struct(ref=T, alt=TAG), 1220),\n",
       " (Struct(ref=TGC, alt=T), 1210),\n",
       " (Struct(ref=ATCT, alt=A), 1199),\n",
       " (Struct(ref=AATT, alt=A), 1189),\n",
       " (Struct(ref=A, alt=AGT), 1185),\n",
       " (Struct(ref=CTCT, alt=C), 1174),\n",
       " (Struct(ref=G, alt=GCA), 1173),\n",
       " (Struct(ref=TTCTC, alt=T), 1151),\n",
       " (Struct(ref=G, alt=GAC), 1133),\n",
       " (Struct(ref=CTGTT, alt=C), 1130),\n",
       " (Struct(ref=C, alt=CCA), 1094),\n",
       " (Struct(ref=AAGG, alt=A), 1089),\n",
       " (Struct(ref=AGAG, alt=A), 1080),\n",
       " (Struct(ref=CGT, alt=C), 1076),\n",
       " (Struct(ref=T, alt=TTTG), 1070),\n",
       " (Struct(ref=ATGT, alt=A), 1038),\n",
       " (Struct(ref=CAAT, alt=C), 1032),\n",
       " (Struct(ref=CAAAT, alt=C), 1018),\n",
       " (Struct(ref=ATTC, alt=A), 999),\n",
       " (Struct(ref=T, alt=TGC), 971),\n",
       " (Struct(ref=AGG, alt=A), 967),\n",
       " (Struct(ref=T, alt=TGG), 964),\n",
       " (Struct(ref=CAGA, alt=C), 945),\n",
       " (Struct(ref=T, alt=TCC), 943),\n",
       " (Struct(ref=G, alt=GGT), 926),\n",
       " (Struct(ref=ATTG, alt=A), 925),\n",
       " (Struct(ref=A, alt=AGG), 917),\n",
       " (Struct(ref=ATATT, alt=A), 917),\n",
       " (Struct(ref=G, alt=GGA), 903),\n",
       " (Struct(ref=TGG, alt=T), 900),\n",
       " (Struct(ref=TAAAG, alt=T), 900),\n",
       " (Struct(ref=TTAAA, alt=T), 889),\n",
       " (Struct(ref=GAGAA, alt=G), 881),\n",
       " (Struct(ref=CAAG, alt=C), 870),\n",
       " (Struct(ref=ATAAT, alt=A), 853),\n",
       " (Struct(ref=TGAA, alt=T), 845),\n",
       " (Struct(ref=AAGAG, alt=A), 845),\n",
       " (Struct(ref=A, alt=AAAC), 844),\n",
       " (Struct(ref=CCTTT, alt=C), 843),\n",
       " (Struct(ref=TCC, alt=T), 840),\n",
       " (Struct(ref=CAGG, alt=C), 837),\n",
       " (Struct(ref=T, alt=TTTTC), 835),\n",
       " (Struct(ref=TAGA, alt=T), 834),\n",
       " (Struct(ref=CTTG, alt=C), 822),\n",
       " (Struct(ref=TGAG, alt=T), 821),\n",
       " (Struct(ref=CTTCT, alt=C), 821),\n",
       " (Struct(ref=G, alt=GCT), 809),\n",
       " (Struct(ref=T, alt=TTTA), 806),\n",
       " (Struct(ref=A, alt=AAAAC), 801),\n",
       " (Struct(ref=ACG, alt=A), 798),\n",
       " (Struct(ref=T, alt=TTTTG), 791),\n",
       " (Struct(ref=A, alt=AAAT), 788),\n",
       " (Struct(ref=CTATT, alt=C), 784),\n",
       " (Struct(ref=GTTC, alt=G), 778),\n",
       " (Struct(ref=TAGAC, alt=T), 775),\n",
       " (Struct(ref=TATG, alt=T), 762),\n",
       " (Struct(ref=ACC, alt=A), 760),\n",
       " (Struct(ref=AAAAAG, alt=A), 758),\n",
       " (Struct(ref=CATA, alt=C), 756),\n",
       " (Struct(ref=AATG, alt=A), 755),\n",
       " (Struct(ref=A, alt=ACC), 754),\n",
       " (Struct(ref=CACAT, alt=C), 754),\n",
       " (Struct(ref=TAATA, alt=T), 736),\n",
       " (Struct(ref=A, alt=AAAG), 719),\n",
       " (Struct(ref=GATA, alt=G), 701),\n",
       " (Struct(ref=AGTT, alt=A), 698),\n",
       " (Struct(ref=TATTA, alt=T), 692),\n",
       " (Struct(ref=TAGAG, alt=T), 691),\n",
       " (Struct(ref=ATATG, alt=A), 689),\n",
       " (Struct(ref=C, alt=CGT), 688),\n",
       " (Struct(ref=TTTAA, alt=T), 685),\n",
       " (Struct(ref=CTTA, alt=C), 678),\n",
       " (Struct(ref=GGC, alt=G), 678),\n",
       " (Struct(ref=T, alt=TTTTA), 675),\n",
       " (Struct(ref=GTTA, alt=G), 673),\n",
       " (Struct(ref=GGAA, alt=G), 659),\n",
       " (Struct(ref=ATTTG, alt=A), 657),\n",
       " (Struct(ref=ATATAT, alt=A), 656),\n",
       " (Struct(ref=TCAAA, alt=T), 654),\n",
       " (Struct(ref=CTTTG, alt=C), 652),\n",
       " (Struct(ref=CACA, alt=C), 650),\n",
       " (Struct(ref=ATTAT, alt=A), 648),\n",
       " (Struct(ref=TATAC, alt=T), 644),\n",
       " (Struct(ref=CACTT, alt=C), 640),\n",
       " (Struct(ref=TGTG, alt=T), 631),\n",
       " (Struct(ref=C, alt=CCTT), 630),\n",
       " (Struct(ref=AAATT, alt=A), 626),\n",
       " (Struct(ref=A, alt=AAAAT), 618),\n",
       " (Struct(ref=CTGT, alt=C), 604),\n",
       " (Struct(ref=TAGAA, alt=T), 601),\n",
       " (Struct(ref=TACA, alt=T), 594),\n",
       " (Struct(ref=TAAA, alt=T), 594),\n",
       " (Struct(ref=A, alt=ACG), 593),\n",
       " (Struct(ref=GTGTA, alt=G), 590),\n",
       " (Struct(ref=ACTC, alt=A), 589),\n",
       " (Struct(ref=TACTC, alt=T), 586),\n",
       " (Struct(ref=CAATT, alt=C), 577),\n",
       " (Struct(ref=TTCC, alt=T), 573),\n",
       " (Struct(ref=TATATA, alt=T), 572),\n",
       " (Struct(ref=TAGG, alt=T), 572),\n",
       " (Struct(ref=TAGTC, alt=T), 572),\n",
       " (Struct(ref=CTCTG, alt=C), 570),\n",
       " (Struct(ref=GAAGA, alt=G), 567),\n",
       " (Struct(ref=AGC, alt=A), 562),\n",
       " (Struct(ref=TATAA, alt=T), 560),\n",
       " (Struct(ref=TTATC, alt=T), 558),\n",
       " (Struct(ref=CTTTA, alt=C), 556),\n",
       " (Struct(ref=CAGAG, alt=C), 552),\n",
       " (Struct(ref=CAAAG, alt=C), 551),\n",
       " (Struct(ref=CTAT, alt=C), 549),\n",
       " (Struct(ref=CTTT, alt=C), 547),\n",
       " (Struct(ref=CAGTT, alt=C), 544),\n",
       " (Struct(ref=A, alt=ACTT), 541),\n",
       " (Struct(ref=ATTGT, alt=A), 535),\n",
       " (Struct(ref=T, alt=TAAC), 522),\n",
       " (Struct(ref=GTCT, alt=G), 521),\n",
       " (Struct(ref=CACAA, alt=C), 519),\n",
       " (Struct(ref=C, alt=CATT), 517),\n",
       " (Struct(ref=AAAAAC, alt=A), 514),\n",
       " (Struct(ref=CATTT, alt=C), 513),\n",
       " (Struct(ref=AGAT, alt=A), 512),\n",
       " (Struct(ref=GAAC, alt=G), 511),\n",
       " (Struct(ref=TTGTC, alt=T), 510),\n",
       " (Struct(ref=GAAT, alt=G), 506),\n",
       " (Struct(ref=CAGAT, alt=C), 495),\n",
       " (Struct(ref=TCAA, alt=T), 494),\n",
       " (Struct(ref=GCTT, alt=G), 492),\n",
       " (Struct(ref=TTTTTC, alt=T), 492),\n",
       " (Struct(ref=G, alt=GGC), 485),\n",
       " (Struct(ref=TTAAG, alt=T), 482),\n",
       " (Struct(ref=A, alt=AATT), 478),\n",
       " (Struct(ref=CATG, alt=C), 478),\n",
       " (Struct(ref=TTTTTG, alt=T), 474),\n",
       " (Struct(ref=AAAGG, alt=A), 474),\n",
       " (Struct(ref=T, alt=TATC), 469),\n",
       " (Struct(ref=T, alt=TTAA), 468),\n",
       " (Struct(ref=AAAAAT, alt=A), 466),\n",
       " (Struct(ref=GTGTT, alt=G), 463),\n",
       " (Struct(ref=G, alt=GCC), 460),\n",
       " (Struct(ref=AATTT, alt=A), 459),\n",
       " (Struct(ref=ATAGT, alt=A), 457),\n",
       " (Struct(ref=ATTTC, alt=A), 456),\n",
       " (Struct(ref=CCAAA, alt=C), 454),\n",
       " (Struct(ref=GACA, alt=G), 454),\n",
       " (Struct(ref=GTTTA, alt=G), 451),\n",
       " (Struct(ref=CCCT, alt=C), 450),\n",
       " (Struct(ref=TACAC, alt=T), 449),\n",
       " (Struct(ref=AATAG, alt=A), 446),\n",
       " (Struct(ref=ACACT, alt=A), 442),\n",
       " (Struct(ref=ACAT, alt=A), 438),\n",
       " (Struct(ref=ATTT, alt=A), 437),\n",
       " (Struct(ref=TGTGA, alt=T), 435),\n",
       " (Struct(ref=GTGT, alt=G), 434),\n",
       " (Struct(ref=GGAAA, alt=G), 430),\n",
       " (Struct(ref=T, alt=TAAG), 429),\n",
       " (Struct(ref=GCC, alt=G), 426),\n",
       " (Struct(ref=CCTCT, alt=C), 422),\n",
       " (Struct(ref=AAGAC, alt=A), 419),\n",
       " (Struct(ref=TTGA, alt=T), 417),\n",
       " (Struct(ref=TTAAC, alt=T), 417),\n",
       " (Struct(ref=A, alt=AAAAG), 416),\n",
       " (Struct(ref=AAAGT, alt=A), 414),\n",
       " (Struct(ref=A, alt=AGC), 413),\n",
       " (Struct(ref=CCG, alt=C), 411),\n",
       " (Struct(ref=CTAA, alt=C), 410),\n",
       " (Struct(ref=AATC, alt=A), 408),\n",
       " (Struct(ref=AGAAG, alt=A), 407),\n",
       " (Struct(ref=GTAT, alt=G), 405),\n",
       " (Struct(ref=GTTGT, alt=G), 403),\n",
       " (Struct(ref=T, alt=TAAAC), 392),\n",
       " (Struct(ref=CTTAT, alt=C), 391),\n",
       " (Struct(ref=ATACT, alt=A), 391),\n",
       " (Struct(ref=TCTTA, alt=T), 384),\n",
       " (Struct(ref=GATAA, alt=G), 383),\n",
       " (Struct(ref=TAATC, alt=T), 383),\n",
       " (Struct(ref=TCTA, alt=T), 382),\n",
       " (Struct(ref=TTTCC, alt=T), 381),\n",
       " (Struct(ref=CCCA, alt=C), 381),\n",
       " (Struct(ref=CTAAT, alt=C), 380),\n",
       " (Struct(ref=GGTT, alt=G), 372),\n",
       " (Struct(ref=A, alt=ATTT), 372),\n",
       " (Struct(ref=GACAA, alt=G), 371),\n",
       " (Struct(ref=ACCT, alt=A), 370),\n",
       " (Struct(ref=T, alt=TAAA), 369),\n",
       " (Struct(ref=TCTCA, alt=T), 368),\n",
       " (Struct(ref=GATT, alt=G), 366),\n",
       " (Struct(ref=AAGT, alt=A), 365),\n",
       " (Struct(ref=ATAG, alt=A), 365),\n",
       " (Struct(ref=ATTCT, alt=A), 365),\n",
       " (Struct(ref=TTAAAA, alt=T), 365),\n",
       " (Struct(ref=TTAGA, alt=T), 362),\n",
       " (Struct(ref=C, alt=CAAAT), 360),\n",
       " (Struct(ref=TCTAA, alt=T), 359),\n",
       " (Struct(ref=ATATC, alt=A), 353),\n",
       " (Struct(ref=CTTAA, alt=C), 352),\n",
       " (Struct(ref=GTAA, alt=G), 352),\n",
       " (Struct(ref=TTTTTA, alt=T), 351),\n",
       " (Struct(ref=TATAG, alt=T), 348),\n",
       " (Struct(ref=C, alt=CTGTT), 342),\n",
       " (Struct(ref=TGATA, alt=T), 340),\n",
       " (Struct(ref=CACAG, alt=C), 339),\n",
       " (Struct(ref=TCTTC, alt=T), 338),\n",
       " (Struct(ref=CCAA, alt=C), 337),\n",
       " (Struct(ref=GGGA, alt=G), 336),\n",
       " (Struct(ref=ACTTT, alt=A), 335),\n",
       " (Struct(ref=ACAC, alt=A), 333),\n",
       " (Struct(ref=C, alt=CCG), 332),\n",
       " (Struct(ref=GAATA, alt=G), 332),\n",
       " (Struct(ref=A, alt=AGAT), 332),\n",
       " (Struct(ref=C, alt=CAAT), 329),\n",
       " (Struct(ref=ATGTT, alt=A), 328),\n",
       " (Struct(ref=CTTTT, alt=C), 327),\n",
       " (Struct(ref=CCAT, alt=C), 327),\n",
       " (Struct(ref=TTCA, alt=T), 327),\n",
       " (Struct(ref=A, alt=ATCT), 327),\n",
       " (Struct(ref=CCCTT, alt=C), 326),\n",
       " (Struct(ref=C, alt=CTATT), 325),\n",
       " (Struct(ref=TGTA, alt=T), 324),\n",
       " (Struct(ref=TGGA, alt=T), 322),\n",
       " (Struct(ref=CAGAA, alt=C), 321),\n",
       " (Struct(ref=TGTAA, alt=T), 320),\n",
       " (Struct(ref=AAGAT, alt=A), 318),\n",
       " (Struct(ref=CTTGT, alt=C), 317),\n",
       " (Struct(ref=GAAAC, alt=G), 315),\n",
       " (Struct(ref=GTTTC, alt=G), 314),\n",
       " (Struct(ref=A, alt=ATTG), 314),\n",
       " (Struct(ref=CTCTA, alt=C), 314),\n",
       " (Struct(ref=AGTTT, alt=A), 310),\n",
       " (Struct(ref=CAAA, alt=C), 309),\n",
       " (Struct(ref=ATAAG, alt=A), 309),\n",
       " (Struct(ref=CTTTCT, alt=C), 308),\n",
       " (Struct(ref=A, alt=ATGT), 308),\n",
       " (Struct(ref=CTCA, alt=C), 307),\n",
       " (Struct(ref=A, alt=ATTC), 306),\n",
       " (Struct(ref=TTAG, alt=T), 305),\n",
       " (Struct(ref=GTATT, alt=G), 304),\n",
       " (Struct(ref=C, alt=CAAA), 302),\n",
       " (Struct(ref=GAAAT, alt=G), 302),\n",
       " (Struct(ref=TAACA, alt=T), 302),\n",
       " (Struct(ref=CATCT, alt=C), 301),\n",
       " (Struct(ref=TTATA, alt=T), 300),\n",
       " (Struct(ref=ATCTG, alt=A), 300),\n",
       " (Struct(ref=C, alt=CAGG), 298),\n",
       " (Struct(ref=ATCTT, alt=A), 298),\n",
       " (Struct(ref=CTAAA, alt=C), 297),\n",
       " (Struct(ref=CGA, alt=C), 297),\n",
       " (Struct(ref=G, alt=GTTA), 297),\n",
       " (Struct(ref=GTCTC, alt=G), 295),\n",
       " (Struct(ref=C, alt=CTAT), 295),\n",
       " (Struct(ref=AACT, alt=A), 295),\n",
       " (Struct(ref=ACTGT, alt=A), 294),\n",
       " (Struct(ref=ATAAC, alt=A), 293),\n",
       " (Struct(ref=A, alt=AGTT), 293),\n",
       " (Struct(ref=GAATT, alt=G), 292),\n",
       " (Struct(ref=ATAC, alt=A), 291),\n",
       " (Struct(ref=AACAG, alt=A), 290),\n",
       " (Struct(ref=TCAC, alt=T), 289),\n",
       " (Struct(ref=AATAT, alt=A), 286),\n",
       " (Struct(ref=ACTCT, alt=A), 285),\n",
       " (Struct(ref=CATAA, alt=C), 285),\n",
       " (Struct(ref=AGAGT, alt=A), 285),\n",
       " (Struct(ref=GAGAC, alt=G), 284),\n",
       " (Struct(ref=TAGC, alt=T), 284),\n",
       " (Struct(ref=CCTAA, alt=C), 284),\n",
       " (Struct(ref=ACTAT, alt=A), 282),\n",
       " (Struct(ref=T, alt=TTTTTG), 282),\n",
       " (Struct(ref=TTAC, alt=T), 280),\n",
       " (Struct(ref=C, alt=CTTA), 280),\n",
       " (Struct(ref=CTGAG, alt=C), 279),\n",
       " (Struct(ref=GTCTT, alt=G), 278),\n",
       " (Struct(ref=TTCTG, alt=T), 277),\n",
       " (Struct(ref=C, alt=CTTT), 274),\n",
       " (Struct(ref=CAGT, alt=C), 273),\n",
       " (Struct(ref=T, alt=TGAA), 272),\n",
       " (Struct(ref=TCTG, alt=T), 272),\n",
       " (Struct(ref=CGTT, alt=C), 271),\n",
       " (Struct(ref=TAAAAC, alt=T), 271),\n",
       " (Struct(ref=ATCTC, alt=A), 271),\n",
       " (Struct(ref=C, alt=CTTTA), 270),\n",
       " (Struct(ref=TGTC, alt=T), 269),\n",
       " (Struct(ref=TAGTG, alt=T), 269),\n",
       " (Struct(ref=TTGG, alt=T), 268),\n",
       " (Struct(ref=TTTGA, alt=T), 268),\n",
       " (Struct(ref=ATGTG, alt=A), 268),\n",
       " (Struct(ref=TAAATA, alt=T), 265),\n",
       " (Struct(ref=CCTGT, alt=C), 264),\n",
       " (Struct(ref=TGTTA, alt=T), 263),\n",
       " (Struct(ref=G, alt=GTTT), 263),\n",
       " (Struct(ref=C, alt=CATA), 261),\n",
       " (Struct(ref=CTGA, alt=C), 260),\n",
       " (Struct(ref=A, alt=AATG), 260),\n",
       " (Struct(ref=ATAAAT, alt=A), 258),\n",
       " (Struct(ref=GTTAT, alt=G), 257),\n",
       " (Struct(ref=G, alt=GATA), 257),\n",
       " (Struct(ref=A, alt=AAGG), 257),\n",
       " (Struct(ref=GTCTA, alt=G), 256),\n",
       " (Struct(ref=TTCTA, alt=T), 256),\n",
       " (Struct(ref=GGTTT, alt=G), 255),\n",
       " (Struct(ref=TATATAC, alt=T), 254),\n",
       " (Struct(ref=C, alt=CAGA), 254),\n",
       " (Struct(ref=GCAA, alt=G), 254),\n",
       " (Struct(ref=TAAGA, alt=T), 252),\n",
       " (Struct(ref=C, alt=CTGT), 251),\n",
       " (Struct(ref=ATCC, alt=A), 250),\n",
       " (Struct(ref=CCTAT, alt=C), 250),\n",
       " (Struct(ref=AAACT, alt=A), 250),\n",
       " (Struct(ref=T, alt=TTAAA), 248),\n",
       " (Struct(ref=TTGAG, alt=T), 248),\n",
       " (Struct(ref=AAGC, alt=A), 247),\n",
       " (Struct(ref=T, alt=TTTTTC), 246),\n",
       " (Struct(ref=TCTGA, alt=T), 245),\n",
       " (Struct(ref=CTGAT, alt=C), 245),\n",
       " (Struct(ref=ATTAC, alt=A), 245),\n",
       " (Struct(ref=ATTAG, alt=A), 244),\n",
       " (Struct(ref=C, alt=CTTG), 244),\n",
       " (Struct(ref=TGAGA, alt=T), 243),\n",
       " (Struct(ref=AAAAAAG, alt=A), 243),\n",
       " (Struct(ref=GATTA, alt=G), 242),\n",
       " (Struct(ref=TACAA, alt=T), 242),\n",
       " (Struct(ref=GTTT, alt=G), 241),\n",
       " (Struct(ref=ATGG, alt=A), 240),\n",
       " (Struct(ref=C, alt=CATTT), 240),\n",
       " (Struct(ref=T, alt=TGAG), 240),\n",
       " (Struct(ref=GTTAA, alt=G), 239),\n",
       " (Struct(ref=CCTGA, alt=C), 239),\n",
       " (Struct(ref=TAAAA, alt=T), 238),\n",
       " (Struct(ref=T, alt=TTCTC), 238),\n",
       " (Struct(ref=GAAA, alt=G), 238),\n",
       " (Struct(ref=TGAAA, alt=T), 237),\n",
       " (Struct(ref=GTGA, alt=G), 237),\n",
       " (Struct(ref=TACC, alt=T), 237),\n",
       " (Struct(ref=TCAGA, alt=T), 237),\n",
       " (Struct(ref=CCTG, alt=C), 236),\n",
       " (Struct(ref=ACAGT, alt=A), 236),\n",
       " (Struct(ref=CCTA, alt=C), 236),\n",
       " (Struct(ref=C, alt=CTAA), 236),\n",
       " (Struct(ref=ACAG, alt=A), 235),\n",
       " (Struct(ref=GCAAA, alt=G), 235),\n",
       " (Struct(ref=AAGTT, alt=A), 232),\n",
       " (Struct(ref=T, alt=TCAAA), 232),\n",
       " (Struct(ref=AACTT, alt=A), 231),\n",
       " (Struct(ref=C, alt=CAAG), 231),\n",
       " (Struct(ref=TCG, alt=T), 229),\n",
       " (Struct(ref=GCCT, alt=G), 228),\n",
       " (Struct(ref=T, alt=TAGA), 228),\n",
       " (Struct(ref=TCCTC, alt=T), 226),\n",
       " (Struct(ref=CAAAA, alt=C), 226),\n",
       " (Struct(ref=C, alt=CACAT), 224),\n",
       " (Struct(ref=G, alt=GTAT), 224),\n",
       " (Struct(ref=GGAGA, alt=G), 223),\n",
       " (Struct(ref=CACT, alt=C), 222),\n",
       " (Struct(ref=A, alt=AAAAAC), 222),\n",
       " (Struct(ref=T, alt=TAAAG), 221),\n",
       " (Struct(ref=T, alt=TTTAA), 221),\n",
       " (Struct(ref=A, alt=ATTTG), 221),\n",
       " (Struct(ref=T, alt=TATAC), 219),\n",
       " (Struct(ref=GTATA, alt=G), 219),\n",
       " (Struct(ref=G, alt=GGAA), 217),\n",
       " (Struct(ref=GCTTT, alt=G), 217),\n",
       " (Struct(ref=GTGTC, alt=G), 216),\n",
       " (Struct(ref=A, alt=ACTC), 216),\n",
       " (Struct(ref=GATTT, alt=G), 216),\n",
       " (Struct(ref=ATTTAT, alt=A), 216),\n",
       " (Struct(ref=G, alt=GAAA), 214),\n",
       " (Struct(ref=T, alt=TACA), 214),\n",
       " (Struct(ref=C, alt=CCTTT), 213),\n",
       " (Struct(ref=TTGAC, alt=T), 212),\n",
       " (Struct(ref=T, alt=TTCC), 211),\n",
       " (Struct(ref=GGGT, alt=G), 211),\n",
       " (Struct(ref=CTCTTT, alt=C), 210),\n",
       " (Struct(ref=CATAT, alt=C), 210),\n",
       " (Struct(ref=TAAAAG, alt=T), 209),\n",
       " (Struct(ref=GCTCT, alt=G), 208),\n",
       " (Struct(ref=GTTCT, alt=G), 208),\n",
       " (Struct(ref=GAAGT, alt=G), 207),\n",
       " (Struct(ref=G, alt=GTGTA), 206),\n",
       " (Struct(ref=C, alt=CTTAT), 206),\n",
       " (Struct(ref=T, alt=TTGA), 206),\n",
       " (Struct(ref=G, alt=GTAA), 206),\n",
       " (Struct(ref=TAATAA, alt=T), 206),\n",
       " (Struct(ref=TTGTG, alt=T), 204),\n",
       " (Struct(ref=TTATAA, alt=T), 204),\n",
       " (Struct(ref=ACAAT, alt=A), 204),\n",
       " (Struct(ref=G, alt=GTTTA), 202),\n",
       " (Struct(ref=GAGAT, alt=G), 201),\n",
       " (Struct(ref=CTCCT, alt=C), 201),\n",
       " (Struct(ref=CTAAG, alt=C), 201),\n",
       " (Struct(ref=G, alt=GTTC), 201),\n",
       " (Struct(ref=T, alt=TAGAG), 200),\n",
       " (Struct(ref=TTTAG, alt=T), 200),\n",
       " (Struct(ref=GAACA, alt=G), 199),\n",
       " (Struct(ref=GTAAA, alt=G), 198),\n",
       " (Struct(ref=T, alt=TCAA), 197),\n",
       " (Struct(ref=GCTC, alt=G), 197),\n",
       " (Struct(ref=CAATA, alt=C), 195),\n",
       " (Struct(ref=GAGGA, alt=G), 194),\n",
       " (Struct(ref=A, alt=ATTTC), 193),\n",
       " (Struct(ref=GTCC, alt=G), 193),\n",
       " (Struct(ref=T, alt=TACAC), 193),\n",
       " (Struct(ref=GATC, alt=G), 193),\n",
       " (Struct(ref=GGGAA, alt=G), 192),\n",
       " (Struct(ref=GCTAA, alt=G), 192),\n",
       " (Struct(ref=CATTA, alt=C), 192),\n",
       " (Struct(ref=AAGGG, alt=A), 190),\n",
       " (Struct(ref=CTTAG, alt=C), 190),\n",
       " (Struct(ref=G, alt=GCTT), 190),\n",
       " (Struct(ref=G, alt=GAAT), 190),\n",
       " (Struct(ref=A, alt=AATC), 189),\n",
       " (Struct(ref=A, alt=ATATG), 189),\n",
       " (Struct(ref=CGG, alt=C), 189),\n",
       " (Struct(ref=C, alt=CGA), 189),\n",
       " (Struct(ref=TATATATA, alt=T), 189),\n",
       " (Struct(ref=ACTG, alt=A), 189),\n",
       " (Struct(ref=AGTG, alt=A), 187),\n",
       " (Struct(ref=GTCA, alt=G), 187),\n",
       " (Struct(ref=AACTC, alt=A), 186),\n",
       " (Struct(ref=AAATG, alt=A), 186),\n",
       " (Struct(ref=TTAAAG, alt=T), 186),\n",
       " (Struct(ref=A, alt=AAAGG), 184),\n",
       " (Struct(ref=ATATATG, alt=A), 184),\n",
       " (Struct(ref=CTGTG, alt=C), 183),\n",
       " (Struct(ref=AACC, alt=A), 183),\n",
       " (Struct(ref=C, alt=CATG), 183),\n",
       " (Struct(ref=CCAGA, alt=C), 183),\n",
       " (Struct(ref=C, alt=CTAAA), 183),\n",
       " (Struct(ref=C, alt=CGG), 182),\n",
       " (Struct(ref=TAGTA, alt=T), 182),\n",
       " (Struct(ref=C, alt=CTTTG), 181),\n",
       " (Struct(ref=TACTA, alt=T), 181),\n",
       " (Struct(ref=TTCCC, alt=T), 181),\n",
       " (Struct(ref=G, alt=GATT), 180),\n",
       " (Struct(ref=TTGC, alt=T), 180),\n",
       " (Struct(ref=CTATA, alt=C), 180),\n",
       " (Struct(ref=ATGAG, alt=A), 180),\n",
       " (Struct(ref=CATAG, alt=C), 180),\n",
       " (Struct(ref=GGATA, alt=G), 180),\n",
       " (Struct(ref=AATAC, alt=A), 179),\n",
       " (Struct(ref=CTTTAT, alt=C), 178),\n",
       " (Struct(ref=C, alt=CAATA), 177),\n",
       " (Struct(ref=GCCA, alt=G), 177),\n",
       " (Struct(ref=CAACA, alt=C), 177),\n",
       " (Struct(ref=G, alt=GTATT), 177),\n",
       " (Struct(ref=C, alt=CAGAT), 177),\n",
       " (Struct(ref=ATTATT, alt=A), 176),\n",
       " (Struct(ref=T, alt=TAGAC), 176),\n",
       " (Struct(ref=GGAT, alt=G), 175),\n",
       " (Struct(ref=T, alt=TTTCC), 175),\n",
       " (Struct(ref=G, alt=GAAC), 173),\n",
       " (Struct(ref=CACACAT, alt=C), 173),\n",
       " (Struct(ref=CAAGT, alt=C), 173),\n",
       " (Struct(ref=A, alt=AATTT), 172),\n",
       " (Struct(ref=ATATTT, alt=A), 171),\n",
       " (Struct(ref=C, alt=CCCT), 171),\n",
       " (Struct(ref=AAAAAAC, alt=A), 169),\n",
       " (Struct(ref=A, alt=AAATT), 169),\n",
       " (Struct(ref=T, alt=TAACA), 168),\n",
       " (Struct(ref=G, alt=GTCT), 168),\n",
       " (Struct(ref=C, alt=CCCA), 168),\n",
       " (Struct(ref=CTTTTA, alt=C), 168),\n",
       " (Struct(ref=T, alt=TAGG), 167),\n",
       " (Struct(ref=CAGTG, alt=C), 167),\n",
       " (Struct(ref=GTGTGTA, alt=G), 166),\n",
       " (Struct(ref=AATCT, alt=A), 166),\n",
       " (Struct(ref=T, alt=TAGAA), 166),\n",
       " (Struct(ref=AGGT, alt=A), 165),\n",
       " (Struct(ref=T, alt=TTAC), 165),\n",
       " (Struct(ref=T, alt=TCG), 165),\n",
       " (Struct(ref=G, alt=GGGA), 164),\n",
       " (Struct(ref=TTTTAA, alt=T), 163),\n",
       " (Struct(ref=CTTCTT, alt=C), 163),\n",
       " (Struct(ref=GCTA, alt=G), 163),\n",
       " (Struct(ref=C, alt=CAAAG), 162),\n",
       " (Struct(ref=CCAG, alt=C), 162),\n",
       " (Struct(ref=TCCC, alt=T), 162),\n",
       " (Struct(ref=CCTTTT, alt=C), 162),\n",
       " (Struct(ref=C, alt=CTTTT), 161),\n",
       " (Struct(ref=C, alt=CCAT), 161),\n",
       " (Struct(ref=AGAC, alt=A), 161),\n",
       " (Struct(ref=C, alt=CTGA), 161),\n",
       " (Struct(ref=C, alt=CAAAA), 161),\n",
       " (Struct(ref=ATTTT, alt=A), 161),\n",
       " (Struct(ref=TATTG, alt=T), 160),\n",
       " (Struct(ref=CAACT, alt=C), 160),\n",
       " (Struct(ref=GCTGT, alt=G), 160),\n",
       " (Struct(ref=TTTAC, alt=T), 159),\n",
       " (Struct(ref=A, alt=ATTTT), 158),\n",
       " (Struct(ref=TGTGC, alt=T), 158),\n",
       " (Struct(ref=GAGT, alt=G), 157),\n",
       " (Struct(ref=TCCA, alt=T), 157),\n",
       " (Struct(ref=C, alt=CAGAG), 157),\n",
       " (Struct(ref=C, alt=CAGTT), 157),\n",
       " (Struct(ref=TGAC, alt=T), 156),\n",
       " (Struct(ref=T, alt=TATA), 156),\n",
       " (Struct(ref=AAATC, alt=A), 156),\n",
       " (Struct(ref=TATTC, alt=T), 156),\n",
       " (Struct(ref=TATATATAC, alt=T), 155),\n",
       " (Struct(ref=T, alt=TCTC), 155),\n",
       " (Struct(ref=T, alt=TTAAG), 155),\n",
       " (Struct(ref=G, alt=GTTTT), 154),\n",
       " (Struct(ref=CCATT, alt=C), 154),\n",
       " (Struct(ref=C, alt=CAGT), 154),\n",
       " (Struct(ref=AGAGG, alt=A), 154),\n",
       " (Struct(ref=TTTCTC, alt=T), 154),\n",
       " (Struct(ref=A, alt=ACAT), 154),\n",
       " (Struct(ref=TCAG, alt=T), 153),\n",
       " (Struct(ref=G, alt=GTATA), 153),\n",
       " (Struct(ref=T, alt=TAAAA), 153),\n",
       " (Struct(ref=GACT, alt=G), 152),\n",
       " (Struct(ref=G, alt=GAAAT), 152),\n",
       " (Struct(ref=ATTTCT, alt=A), 151),\n",
       " (Struct(ref=TGGG, alt=T), 149),\n",
       " (Struct(ref=AGCT, alt=A), 149),\n",
       " (Struct(ref=TCTCC, alt=T), 149),\n",
       " (Struct(ref=A, alt=ATCTT), 149),\n",
       " (Struct(ref=AAATAT, alt=A), 149),\n",
       " (Struct(ref=A, alt=ATGTT), 148),\n",
       " (Struct(ref=A, alt=ATTGT), 148),\n",
       " (Struct(ref=C, alt=CCCTT), 148),\n",
       " (Struct(ref=ATAACT, alt=A), 147),\n",
       " (Struct(ref=TTATG, alt=T), 147),\n",
       " (Struct(ref=T, alt=TTCA), 147),\n",
       " (Struct(ref=A, alt=ATGTG), 146),\n",
       " (Struct(ref=CAAAAT, alt=C), 146),\n",
       " (Struct(ref=GTAGA, alt=G), 146),\n",
       " (Struct(ref=AATTG, alt=A), 146),\n",
       " (Struct(ref=GCTAT, alt=G), 145),\n",
       " (Struct(ref=C, alt=CATAT), 144),\n",
       " (Struct(ref=G, alt=GGAAA), 144),\n",
       " (Struct(ref=ATAAAG, alt=A), 143),\n",
       " (Struct(ref=AAAGAG, alt=A), 143),\n",
       " (Struct(ref=AGGAG, alt=A), 143),\n",
       " (Struct(ref=TAATG, alt=T), 141),\n",
       " (Struct(ref=TTTAAA, alt=T), 140),\n",
       " (Struct(ref=CCTCA, alt=C), 140),\n",
       " (Struct(ref=A, alt=AAGAG), 140),\n",
       " (Struct(ref=TATATC, alt=T), 140),\n",
       " (Struct(ref=A, alt=ACTTT), 140),\n",
       " (Struct(ref=GGGC, alt=G), 140),\n",
       " (Struct(ref=T, alt=TTAG), 139),\n",
       " (Struct(ref=A, alt=ATAT), 139),\n",
       " (Struct(ref=T, alt=TATG), 139),\n",
       " (Struct(ref=GAAAA, alt=G), 139),\n",
       " (Struct(ref=GCTGA, alt=G), 138),\n",
       " (Struct(ref=GAACT, alt=G), 138),\n",
       " (Struct(ref=A, alt=AATAT), 137),\n",
       " (Struct(ref=A, alt=ATATC), 136),\n",
       " (Struct(ref=AAAATT, alt=A), 136),\n",
       " (Struct(ref=C, alt=CAATT), 136),\n",
       " (Struct(ref=TATAAC, alt=T), 135),\n",
       " (Struct(ref=A, alt=ACCT), 135),\n",
       " (Struct(ref=GTATC, alt=G), 135),\n",
       " (Struct(ref=GAAAGA, alt=G), 135),\n",
       " (Struct(ref=T, alt=TTATA), 135),\n",
       " (Struct(ref=G, alt=GATAA), 135),\n",
       " (Struct(ref=GTTTT, alt=G), 134),\n",
       " (Struct(ref=AACTG, alt=A), 134),\n",
       " (Struct(ref=G, alt=GCCT), 134),\n",
       " (Struct(ref=TGCA, alt=T), 133),\n",
       " (Struct(ref=C, alt=CTCTG), 133),\n",
       " (Struct(ref=CCCG, alt=C), 133),\n",
       " (Struct(ref=TTTTTTG, alt=T), 133),\n",
       " (Struct(ref=ACTTTT, alt=A), 133),\n",
       " (Struct(ref=T, alt=TGGA), 132),\n",
       " (Struct(ref=C, alt=CTAAT), 131),\n",
       " (Struct(ref=C, alt=CCAA), 131),\n",
       " (Struct(ref=A, alt=ATAG), 130),\n",
       " (Struct(ref=AACAC, alt=A), 130),\n",
       " (Struct(ref=AAAAAAT, alt=A), 130),\n",
       " (Struct(ref=T, alt=TTTTTA), 130),\n",
       " (Struct(ref=A, alt=ATCTG), 129),\n",
       " (Struct(ref=CTGG, alt=C), 129),\n",
       " (Struct(ref=CATTTT, alt=C), 129),\n",
       " (Struct(ref=ATTAAT, alt=A), 129),\n",
       " (Struct(ref=GGCT, alt=G), 129),\n",
       " (Struct(ref=AATTTT, alt=A), 128),\n",
       " (Struct(ref=T, alt=TAGC), 128),\n",
       " (Struct(ref=CCTTA, alt=C), 128),\n",
       " (Struct(ref=G, alt=GACA), 128),\n",
       " (Struct(ref=CTTATT, alt=C), 127),\n",
       " (Struct(ref=AATAAT, alt=A), 127),\n",
       " (Struct(ref=G, alt=GGGT), 127),\n",
       " (Struct(ref=C, alt=CACTT), 127),\n",
       " (Struct(ref=C, alt=CTCTA), 127),\n",
       " (Struct(ref=G, alt=GCAA), 126),\n",
       " (Struct(ref=G, alt=GTGA), 126),\n",
       " (Struct(ref=A, alt=AACTT), 125),\n",
       " (Struct(ref=AAAAAAAG, alt=A), 125),\n",
       " (Struct(ref=TAACTC, alt=T), 125),\n",
       " (Struct(ref=C, alt=CTTAA), 125),\n",
       " (Struct(ref=ATAATT, alt=A), 124),\n",
       " (Struct(ref=T, alt=TGTA), 124),\n",
       " (Struct(ref=ATGC, alt=A), 124),\n",
       " (Struct(ref=T, alt=TCTA), 124),\n",
       " (Struct(ref=TATAAA, alt=T), 123),\n",
       " (Struct(ref=CATATT, alt=C), 123),\n",
       " (Struct(ref=T, alt=TTATC), 123),\n",
       " (Struct(ref=T, alt=TTGG), 123),\n",
       " (Struct(ref=A, alt=ATAC), 123),\n",
       " (Struct(ref=A, alt=ATTCT), 123),\n",
       " (Struct(ref=G, alt=GAATA), 122),\n",
       " (Struct(ref=T, alt=TGGG), 122),\n",
       " (Struct(ref=TACAG, alt=T), 121),\n",
       " (Struct(ref=G, alt=GAAAC), 121),\n",
       " (Struct(ref=TACTG, alt=T), 121),\n",
       " (Struct(ref=GTAAT, alt=G), 121),\n",
       " (Struct(ref=AATGT, alt=A), 120),\n",
       " (Struct(ref=AATTAT, alt=A), 120),\n",
       " (Struct(ref=G, alt=GTAAA), 120),\n",
       " (Struct(ref=TTTTTTC, alt=T), 120),\n",
       " (Struct(ref=A, alt=ATGG), 119),\n",
       " (Struct(ref=CCAAT, alt=C), 119),\n",
       " (Struct(ref=GAGTT, alt=G), 119),\n",
       " (Struct(ref=ACACG, alt=A), 119),\n",
       " (Struct(ref=ATACTT, alt=A), 119),\n",
       " (Struct(ref=GGTA, alt=G), 119),\n",
       " (Struct(ref=G, alt=GAGAT), 118),\n",
       " (Struct(ref=T, alt=TCAC), 118),\n",
       " (Struct(ref=G, alt=GTCA), 118),\n",
       " (Struct(ref=TTGAA, alt=T), 118),\n",
       " (Struct(ref=G, alt=GCTA), 117),\n",
       " (Struct(ref=GCAT, alt=G), 117),\n",
       " (Struct(ref=T, alt=TGCA), 117),\n",
       " (Struct(ref=ATTTTC, alt=A), 117),\n",
       " (Struct(ref=A, alt=AGTG), 117),\n",
       " (Struct(ref=C, alt=CCAAA), 117),\n",
       " (Struct(ref=TCTTTA, alt=T), 117),\n",
       " (Struct(ref=GAAGAA, alt=G), 116),\n",
       " (Struct(ref=TTTAAG, alt=T), 115),\n",
       " (Struct(ref=CTGAA, alt=C), 115),\n",
       " (Struct(ref=TATTTA, alt=T), 115),\n",
       " (Struct(ref=GGCA, alt=G), 115),\n",
       " (Struct(ref=C, alt=CTATA), 115),\n",
       " (Struct(ref=CAAATT, alt=C), 115),\n",
       " (Struct(ref=A, alt=AAGAT), 115),\n",
       " (Struct(ref=GACC, alt=G), 114),\n",
       " (Struct(ref=T, alt=TATCA), 114),\n",
       " (Struct(ref=A, alt=AAAAAT), 114),\n",
       " (Struct(ref=CACTA, alt=C), 114),\n",
       " (Struct(ref=G, alt=GTCTT), 113),\n",
       " (Struct(ref=T, alt=TAAAAC), 113),\n",
       " (Struct(ref=TCACA, alt=T), 113),\n",
       " (Struct(ref=TTGTA, alt=T), 113),\n",
       " (Struct(ref=T, alt=TCCC), 113),\n",
       " (Struct(ref=CTCAT, alt=C), 112),\n",
       " (Struct(ref=CTTTAA, alt=C), 112),\n",
       " (Struct(ref=GACTT, alt=G), 112),\n",
       " (Struct(ref=CATGT, alt=C), 112),\n",
       " (Struct(ref=CAGTA, alt=C), 111),\n",
       " (Struct(ref=GCAC, alt=G), 111),\n",
       " (Struct(ref=C, alt=CTGAT), 111),\n",
       " (Struct(ref=CTTCA, alt=C), 111),\n",
       " (Struct(ref=GACTA, alt=G), 111),\n",
       " (Struct(ref=TAAATC, alt=T), 110),\n",
       " (Struct(ref=A, alt=AGAG), 110),\n",
       " (Struct(ref=T, alt=TATAG), 110),\n",
       " (Struct(ref=C, alt=CATAA), 110),\n",
       " (Struct(ref=TATCA, alt=T), 110),\n",
       " (Struct(ref=T, alt=TGTAA), 110),\n",
       " (Struct(ref=TTATTA, alt=T), 109),\n",
       " (Struct(ref=ATGTC, alt=A), 109),\n",
       " (Struct(ref=ATGAT, alt=A), 109),\n",
       " (Struct(ref=T, alt=TACTC), 109),\n",
       " (Struct(ref=A, alt=ATAAG), 109),\n",
       " (Struct(ref=TGGC, alt=T), 108),\n",
       " (Struct(ref=A, alt=AACT), 108),\n",
       " (Struct(ref=ATATATT, alt=A), 108),\n",
       " (Struct(ref=A, alt=AGTTT), 108),\n",
       " (Struct(ref=G, alt=GTTAT), 108),\n",
       " (Struct(ref=GTTAC, alt=G), 107),\n",
       " (Struct(ref=GAGC, alt=G), 107),\n",
       " (Struct(ref=TTTCA, alt=T), 107),\n",
       " (Struct(ref=AAGTC, alt=A), 107),\n",
       " (Struct(ref=TGAAAA, alt=T), 107),\n",
       " (Struct(ref=AGTGT, alt=A), 107),\n",
       " (Struct(ref=A, alt=AGGG), 106),\n",
       " (Struct(ref=A, alt=ATCC), 106),\n",
       " (Struct(ref=G, alt=GGTT), 106),\n",
       " (Struct(ref=TTACTC, alt=T), 106),\n",
       " (Struct(ref=A, alt=ATAGT), 106),\n",
       " (Struct(ref=CTCAA, alt=C), 105),\n",
       " (Struct(ref=A, alt=ACTG), 105),\n",
       " (Struct(ref=CTCAG, alt=C), 105),\n",
       " (Struct(ref=GGTGA, alt=G), 105),\n",
       " (Struct(ref=A, alt=ACAG), 104),\n",
       " (Struct(ref=G, alt=GTTTC), 104),\n",
       " (Struct(ref=A, alt=AAAAAG), 104),\n",
       " (Struct(ref=A, alt=AAGT), 104),\n",
       " (Struct(ref=A, alt=AAGC), 104),\n",
       " (Struct(ref=GCAGA, alt=G), 104),\n",
       " (Struct(ref=G, alt=GTCTC), 104),\n",
       " (Struct(ref=CACTG, alt=C), 103),\n",
       " (Struct(ref=GTTTAT, alt=G), 103),\n",
       " (Struct(ref=CTTGA, alt=C), 103),\n",
       " (Struct(ref=ATATATATG, alt=A), 102),\n",
       " (Struct(ref=T, alt=TTGTC), 102),\n",
       " (Struct(ref=GGTAA, alt=G), 102),\n",
       " (Struct(ref=C, alt=CTCT), 102),\n",
       " (Struct(ref=AGGG, alt=A), 102),\n",
       " (Struct(ref=CAAGA, alt=C), 102),\n",
       " (Struct(ref=GTAAC, alt=G), 101),\n",
       " (Struct(ref=ATATATAT, alt=A), 101),\n",
       " (Struct(ref=A, alt=ATACT), 100),\n",
       " (Struct(ref=T, alt=TCCA), 100),\n",
       " (Struct(ref=GTAGT, alt=G), 100),\n",
       " (Struct(ref=C, alt=CTTGT), 100),\n",
       " (Struct(ref=G, alt=GGAT), 100),\n",
       " (Struct(ref=CTAGA, alt=C), 100),\n",
       " (Struct(ref=T, alt=TCAG), 100),\n",
       " (Struct(ref=AGATT, alt=A), 99),\n",
       " (Struct(ref=C, alt=CTGG), 99),\n",
       " (Struct(ref=T, alt=TTGC), 99),\n",
       " (Struct(ref=T, alt=TAGTC), 99),\n",
       " (Struct(ref=A, alt=ACCC), 99),\n",
       " (Struct(ref=C, alt=CTCA), 99),\n",
       " (Struct(ref=TGACA, alt=T), 98),\n",
       " (Struct(ref=G, alt=GATTT), 98),\n",
       " (Struct(ref=GTTTGT, alt=G), 98),\n",
       " (Struct(ref=CTTAAT, alt=C), 97),\n",
       " (Struct(ref=A, alt=ATCTC), 97),\n",
       " (Struct(ref=TTACA, alt=T), 97),\n",
       " (Struct(ref=T, alt=TTAAAA), 97),\n",
       " (Struct(ref=TTAATA, alt=T), 97),\n",
       " (Struct(ref=GATAAA, alt=G), 97),\n",
       " (Struct(ref=C, alt=CCTAT), 97),\n",
       " (Struct(ref=ATGAC, alt=A), 96),\n",
       " (Struct(ref=G, alt=GCTC), 96),\n",
       " (Struct(ref=C, alt=CTGAG), 96),\n",
       " (Struct(ref=G, alt=GGCT), 96),\n",
       " (Struct(ref=TAAAGA, alt=T), 96),\n",
       " (Struct(ref=T, alt=TGTC), 96),\n",
       " (Struct(ref=AGTC, alt=A), 95),\n",
       " (Struct(ref=TGTCA, alt=T), 95),\n",
       " (Struct(ref=T, alt=TCTAA), 95),\n",
       " (Struct(ref=CGAG, alt=C), 95),\n",
       " (Struct(ref=GATAT, alt=G), 95),\n",
       " (Struct(ref=C, alt=CACAG), 95),\n",
       " (Struct(ref=CTATTT, alt=C), 94),\n",
       " (Struct(ref=AGAAAG, alt=A), 94),\n",
       " (Struct(ref=T, alt=TTAGA), 93),\n",
       " (Struct(ref=GTAC, alt=G), 93),\n",
       " (Struct(ref=T, alt=TGTGA), 93),\n",
       " (Struct(ref=TATATAA, alt=T), 93),\n",
       " (Struct(ref=TTAAAC, alt=T), 93),\n",
       " (Struct(ref=TATAAG, alt=T), 93),\n",
       " (Struct(ref=TCTTTC, alt=T), 93),\n",
       " (Struct(ref=A, alt=AACC), 92),\n",
       " (Struct(ref=GTGAA, alt=G), 92),\n",
       " (Struct(ref=CAAACA, alt=C), 92),\n",
       " (Struct(ref=T, alt=TTCTA), 92),\n",
       " (Struct(ref=C, alt=CAGAA), 91),\n",
       " (Struct(ref=GCTTA, alt=G), 91),\n",
       " (Struct(ref=ACATT, alt=A), 91),\n",
       " (Struct(ref=A, alt=ACTCT), 91),\n",
       " (Struct(ref=T, alt=TACC), 91),\n",
       " (Struct(ref=G, alt=GGTTT), 90),\n",
       " (Struct(ref=C, alt=CCTG), 90),\n",
       " (Struct(ref=A, alt=AATAG), 90),\n",
       " (Struct(ref=A, alt=AAGGG), 90),\n",
       " (Struct(ref=T, alt=TAATC), 90),\n",
       " (Struct(ref=TTCAC, alt=T), 90),\n",
       " (Struct(ref=G, alt=GAATT), 90),\n",
       " (Struct(ref=G, alt=GGGC), 90),\n",
       " (Struct(ref=TAAGAA, alt=T), 89),\n",
       " (Struct(ref=G, alt=GCAT), 89),\n",
       " (Struct(ref=GCACA, alt=G), 89),\n",
       " (Struct(ref=C, alt=CCTGT), 89),\n",
       " (Struct(ref=TAAGG, alt=T), 89),\n",
       " (Struct(ref=CCCCA, alt=C), 89),\n",
       " (Struct(ref=T, alt=TGGC), 89),\n",
       " (Struct(ref=T, alt=TATTC), 89),\n",
       " (Struct(ref=ATTACT, alt=A), 88),\n",
       " (Struct(ref=CCAGT, alt=C), 88),\n",
       " (Struct(ref=G, alt=GCAAA), 88),\n",
       " (Struct(ref=T, alt=TGAAA), 88),\n",
       " (Struct(ref=C, alt=CACT), 88),\n",
       " (Struct(ref=A, alt=ATGC), 88),\n",
       " (Struct(ref=AAAAAAAC, alt=A), 87),\n",
       " (Struct(ref=C, alt=CTAAG), 87),\n",
       " (Struct(ref=G, alt=GTCTA), 87),\n",
       " (Struct(ref=G, alt=GACAA), 87),\n",
       " (Struct(ref=TTTTTTA, alt=T), 87),\n",
       " (Struct(ref=C, alt=CCATT), 87),\n",
       " (Struct(ref=T, alt=TCTG), 87),\n",
       " (Struct(ref=A, alt=AAATG), 86),\n",
       " (Struct(ref=C, alt=CTGAA), 86),\n",
       " (Struct(ref=G, alt=GTGTC), 86),\n",
       " (Struct(ref=TGCC, alt=T), 86),\n",
       " (Struct(ref=CCACT, alt=C), 86),\n",
       " (Struct(ref=A, alt=ACTGT), 86),\n",
       " (Struct(ref=ATAGG, alt=A), 86),\n",
       " (Struct(ref=G, alt=GATC), 86),\n",
       " (Struct(ref=A, alt=AGCT), 86),\n",
       " (Struct(ref=CTCTCT, alt=C), 86),\n",
       " (Struct(ref=T, alt=TGTGC), 85),\n",
       " (Struct(ref=CTAGT, alt=C), 85),\n",
       " (Struct(ref=TCAAAA, alt=T), 85),\n",
       " (Struct(ref=ATTAAG, alt=A), 85),\n",
       " (Struct(ref=A, alt=AAGAC), 85),\n",
       " (Struct(ref=TCATC, alt=T), 85),\n",
       " (Struct(ref=AAAGTT, alt=A), 85),\n",
       " (Struct(ref=TTCAA, alt=T), 85),\n",
       " (Struct(ref=G, alt=GCCA), 85),\n",
       " (Struct(ref=TATATATATAC, alt=T), 84),\n",
       " (Struct(ref=C, alt=CTCTT), 84),\n",
       " (Struct(ref=CGTG, alt=C), 84),\n",
       " (Struct(ref=CATATA, alt=C), 84),\n",
       " (Struct(ref=TAAGAG, alt=T), 84),\n",
       " (Struct(ref=TTTATA, alt=T), 84),\n",
       " (Struct(ref=TATATATATA, alt=T), 84),\n",
       " (Struct(ref=TAAGTA, alt=T), 84),\n",
       " (Struct(ref=CATAAA, alt=C), 84),\n",
       " (Struct(ref=GTGC, alt=G), 84),\n",
       " (Struct(ref=T, alt=TGAC), 83),\n",
       " (Struct(ref=C, alt=CCCG), 83),\n",
       " (Struct(ref=G, alt=GTTAA), 83),\n",
       " (Struct(ref=TATTAA, alt=T), 83),\n",
       " (Struct(ref=G, alt=GTAC), 82),\n",
       " (Struct(ref=AGAAT, alt=A), 82),\n",
       " (Struct(ref=T, alt=TAAGA), 81),\n",
       " (Struct(ref=C, alt=CCTAA), 81),\n",
       " (Struct(ref=GTTTTC, alt=G), 81),\n",
       " (Struct(ref=CTGTA, alt=C), 81),\n",
       " (Struct(ref=CTTAAA, alt=C), 81),\n",
       " (Struct(ref=ATATGT, alt=A), 81),\n",
       " (Struct(ref=AGATG, alt=A), 81),\n",
       " (Struct(ref=CATAAT, alt=C), 81),\n",
       " (Struct(ref=A, alt=AGGT), 81),\n",
       " (Struct(ref=AGGC, alt=A), 80),\n",
       " (Struct(ref=A, alt=ATAAC), 80),\n",
       " (Struct(ref=TAATTA, alt=T), 80),\n",
       " (Struct(ref=CAAAAG, alt=C), 80),\n",
       " (Struct(ref=ATATCT, alt=A), 80),\n",
       " (Struct(ref=GACAC, alt=G), 80),\n",
       " (Struct(ref=C, alt=CCAG), 80),\n",
       " (Struct(ref=AATATT, alt=A), 79),\n",
       " (Struct(ref=CGAT, alt=C), 79),\n",
       " (Struct(ref=GAGTA, alt=G), 79),\n",
       " (Struct(ref=G, alt=GATTA), 79),\n",
       " (Struct(ref=CTAG, alt=C), 79),\n",
       " (Struct(ref=CTACT, alt=C), 79),\n",
       " (Struct(ref=CTTTGT, alt=C), 79),\n",
       " (Struct(ref=C, alt=CTGTG), 79),\n",
       " (Struct(ref=GCCC, alt=G), 79),\n",
       " (Struct(ref=T, alt=TCTCA), 79),\n",
       " (Struct(ref=CTGTTT, alt=C), 78),\n",
       " (Struct(ref=TTAAGA, alt=T), 78),\n",
       " (Struct(ref=AAGAAAG, alt=A), 78),\n",
       " (Struct(ref=T, alt=TTGAA), 78),\n",
       " (Struct(ref=C, alt=CTTTAT), 78),\n",
       " (Struct(ref=GTGTGTGTA, alt=G), 78),\n",
       " (Struct(ref=ACCC, alt=A), 78),\n",
       " (Struct(ref=CAACTT, alt=C), 78),\n",
       " (Struct(ref=CGTGT, alt=C), 78),\n",
       " (Struct(ref=CTAATT, alt=C), 78),\n",
       " (Struct(ref=C, alt=CCTCT), 77),\n",
       " (Struct(ref=CCCCT, alt=C), 77),\n",
       " (Struct(ref=T, alt=TTTGA), 77),\n",
       " (Struct(ref=A, alt=AGGC), 77),\n",
       " (Struct(ref=TGAAG, alt=T), 77),\n",
       " (Struct(ref=CTATG, alt=C), 77),\n",
       " (Struct(ref=AAAAAAAAG, alt=A), 77),\n",
       " (Struct(ref=G, alt=GAGAC), 77),\n",
       " (Struct(ref=T, alt=TTAAC), 77),\n",
       " (Struct(ref=T, alt=TCTTA), 76),\n",
       " (Struct(ref=G, alt=GTCC), 76),\n",
       " (Struct(ref=AATAAG, alt=A), 76),\n",
       " (Struct(ref=G, alt=GGGAA), 75),\n",
       " (Struct(ref=AAAGC, alt=A), 75),\n",
       " (Struct(ref=GTTTTA, alt=G), 75),\n",
       " (Struct(ref=GAAATA, alt=G), 75),\n",
       " (Struct(ref=C, alt=CATCT), 75),\n",
       " (Struct(ref=GGTTA, alt=G), 75),\n",
       " (Struct(ref=AGACT, alt=A), 75),\n",
       " (Struct(ref=C, alt=CCTA), 75),\n",
       " ...]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "counts = Counter(snp_counts)\n",
    "counts.most_common()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Take a subset of the matrix"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "indices = list(range(mt.count_cols()))\n",
    "import random\n",
    "random.shuffle(indices)\n",
    "mt_subset = mt.head(18000)\n",
    "mt_subset.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QC is entirely based on the ability to understand the properties of a dataset. In Hail, the `sample_qc` method produces a set of useful metrics and stores them in a column field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.sample_qc(mt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "  <div class=\"bk-root\" id=\"0f5c6d6e-0045-4c2a-806d-737c17b2157d\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "    \n",
       "  var docs_json = {\"6f90481d-9ef1-4c0d-ba77-cb1cf9e67148\":{\"roots\":{\"references\":[{\"attributes\":{},\"id\":\"0a4c8494-a0fb-4f89-951d-a55b73a821f7\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{\"axis_label\":\"Frequency\",\"formatter\":{\"id\":\"3f7b6b33-4555-49f4-936a-afeb8a04fc82\",\"type\":\"BasicTickFormatter\"},\"plot\":{\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"},\"ticker\":{\"id\":\"af01b002-0cc2-4043-ac15-033a205960dd\",\"type\":\"BasicTicker\"}},\"id\":\"c846c93f-bb74-499e-bf8b-23be03d078b8\",\"type\":\"LinearAxis\"},{\"attributes\":{},\"id\":\"e7303869-d69a-422f-8128-de94e72aacd3\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"source\":{\"id\":\"9b35bd99-5f78-4cd4-a955-70a93c350f11\",\"type\":\"ColumnDataSource\"}},\"id\":\"d4ffd088-21c4-4134-850d-e05653bdf961\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"3f7b6b33-4555-49f4-936a-afeb8a04fc82\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"16d3b1da-4a76-452d-a869-076abc7cce62\",\"type\":\"ResetTool\"},{\"attributes\":{\"background_fill_color\":{\"value\":\"#EEEEEE\"},\"below\":[{\"id\":\"230daed0-530d-4f82-9a93-38ef53070e11\",\"type\":\"LinearAxis\"}],\"left\":[{\"id\":\"c846c93f-bb74-499e-bf8b-23be03d078b8\",\"type\":\"LinearAxis\"}],\"renderers\":[{\"id\":\"230daed0-530d-4f82-9a93-38ef53070e11\",\"type\":\"LinearAxis\"},{\"id\":\"155873b8-b3ed-43d0-aa92-527f5f2fc87e\",\"type\":\"Grid\"},{\"id\":\"c846c93f-bb74-499e-bf8b-23be03d078b8\",\"type\":\"LinearAxis\"},{\"id\":\"abc0a935-6d54-4aa1-a045-d8ff29f231d7\",\"type\":\"Grid\"},{\"id\":\"39765548-e403-4b15-8a82-140a11acec48\",\"type\":\"BoxAnnotation\"},{\"id\":\"50364971-b771-4a88-8a4d-bc799b773ae9\",\"type\":\"Legend\"},{\"id\":\"7e7e8f3c-3087-474b-bcfc-64fbfb98f6e3\",\"type\":\"GlyphRenderer\"}],\"title\":null,\"toolbar\":{\"id\":\"c0178539-8dc8-4d00-bf05-8d5e10b86d9d\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"691a968d-38f2-4326-a845-6672da7c84a6\",\"type\":\"DataRange1d\"},\"x_scale\":{\"id\":\"49a504ec-3f54-4d67-b3d7-940a042819b6\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"ce71a8e3-acf6-4a9e-8bce-709235989ef0\",\"type\":\"DataRange1d\"},\"y_scale\":{\"id\":\"4ce015a4-0f7f-46f4-8aa6-b3dcb31442b3\",\"type\":\"LinearScale\"}},\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{\"callback\":null},\"id\":\"ce71a8e3-acf6-4a9e-8bce-709235989ef0\",\"type\":\"DataRange1d\"},{\"attributes\":{\"overlay\":{\"id\":\"39765548-e403-4b15-8a82-140a11acec48\",\"type\":\"BoxAnnotation\"}},\"id\":\"8ea85a7d-7274-4c89-8f4e-6dc78eae0482\",\"type\":\"BoxZoomTool\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"plot\":null,\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"39765548-e403-4b15-8a82-140a11acec48\",\"type\":\"BoxAnnotation\"},{\"attributes\":{\"axis_label\":\"Call Rate\",\"formatter\":{\"id\":\"0a4c8494-a0fb-4f89-951d-a55b73a821f7\",\"type\":\"BasicTickFormatter\"},\"plot\":{\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"},\"ticker\":{\"id\":\"83a0698b-c844-49a7-853f-94d77229812d\",\"type\":\"BasicTicker\"}},\"id\":\"230daed0-530d-4f82-9a93-38ef53070e11\",\"type\":\"LinearAxis\"},{\"attributes\":{\"callback\":null,\"data\":{\"left\":[0.88,0.8824,0.8848,0.8872,0.8896000000000001,0.892,0.8944,0.8968,0.8992,0.9016,0.904,0.9064,0.9088,0.9112,0.9136,0.916,0.9184,0.9208,0.9232,0.9256,0.928,0.9304,0.9328,0.9352,0.9376,0.94,0.9424,0.9448,0.9472,0.9496,0.952,0.9544,0.9568,0.9592,0.9616,0.964,0.9664,0.9688,0.9712,0.9736,0.976,0.9783999999999999,0.9808,0.9832,0.9856,0.988,0.9904,0.9928,0.9952,0.9976],\"right\":[0.8824,0.8848,0.8872,0.8896000000000001,0.892,0.8944,0.8968,0.8992,0.9016,0.904,0.9064,0.9088,0.9112,0.9136,0.916,0.9184,0.9208,0.9232,0.9256,0.928,0.9304,0.9328,0.9352,0.9376,0.94,0.9424,0.9448,0.9472,0.9496,0.952,0.9544,0.9568,0.9592,0.9616,0.964,0.9664,0.9688,0.9712,0.9736,0.976,0.9783999999999999,0.9808,0.9832,0.9856,0.988,0.9904,0.9928,0.9952,0.9976,1.0],\"top\":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43359732780]},\"selected\":{\"id\":\"3d6269c6-f267-45e3-9016-d61f15ccc184\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"e7303869-d69a-422f-8128-de94e72aacd3\",\"type\":\"UnionRenderers\"}},\"id\":\"9b35bd99-5f78-4cd4-a955-70a93c350f11\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"638b93bf-c19b-43ac-ae1b-43bd256611d6\",\"type\":\"HelpTool\"},{\"attributes\":{\"bottom\":{\"value\":0},\"fill_color\":{\"value\":\"#1f77b4\"},\"left\":{\"field\":\"left\"},\"right\":{\"field\":\"right\"},\"top\":{\"field\":\"top\"}},\"id\":\"7689762c-3eb4-4f37-a364-4c0653bdb9ea\",\"type\":\"Quad\"},{\"attributes\":{},\"id\":\"3d6269c6-f267-45e3-9016-d61f15ccc184\",\"type\":\"Selection\"},{\"attributes\":{\"dimension\":1,\"plot\":{\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"},\"ticker\":{\"id\":\"af01b002-0cc2-4043-ac15-033a205960dd\",\"type\":\"BasicTicker\"}},\"id\":\"abc0a935-6d54-4aa1-a045-d8ff29f231d7\",\"type\":\"Grid\"},{\"attributes\":{},\"id\":\"3f8c95f7-5dab-437d-ab0f-73677c521373\",\"type\":\"WheelZoomTool\"},{\"attributes\":{},\"id\":\"af01b002-0cc2-4043-ac15-033a205960dd\",\"type\":\"BasicTicker\"},{\"attributes\":{\"items\":[{\"id\":\"81c125b2-9b32-47ed-aeef-9f0878666295\",\"type\":\"LegendItem\"}],\"plot\":{\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"}},\"id\":\"50364971-b771-4a88-8a4d-bc799b773ae9\",\"type\":\"Legend\"},{\"attributes\":{\"plot\":{\"id\":\"067941eb-9019-4ca4-960d-bd1aa37d3872\",\"subtype\":\"Figure\",\"type\":\"Plot\"},\"ticker\":{\"id\":\"83a0698b-c844-49a7-853f-94d77229812d\",\"type\":\"BasicTicker\"}},\"id\":\"155873b8-b3ed-43d0-aa92-527f5f2fc87e\",\"type\":\"Grid\"},{\"attributes\":{},\"id\":\"4ce015a4-0f7f-46f4-8aa6-b3dcb31442b3\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"3e657cd6-83f4-4b6c-9cfb-778b42ace5e6\",\"type\":\"PanTool\"},{\"attributes\":{},\"id\":\"83a0698b-c844-49a7-853f-94d77229812d\",\"type\":\"BasicTicker\"},{\"attributes\":{\"data_source\":{\"id\":\"9b35bd99-5f78-4cd4-a955-70a93c350f11\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"7689762c-3eb4-4f37-a364-4c0653bdb9ea\",\"type\":\"Quad\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"972782eb-8c02-4ac0-921d-96ebf84c4bc5\",\"type\":\"Quad\"},\"selection_glyph\":null,\"view\":{\"id\":\"d4ffd088-21c4-4134-850d-e05653bdf961\",\"type\":\"CDSView\"}},\"id\":\"7e7e8f3c-3087-474b-bcfc-64fbfb98f6e3\",\"type\":\"GlyphRenderer\"},{\"attributes\":{},\"id\":\"76d0a161-5e08-482f-aa43-e9ffc7662884\",\"type\":\"SaveTool\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"3e657cd6-83f4-4b6c-9cfb-778b42ace5e6\",\"type\":\"PanTool\"},{\"id\":\"3f8c95f7-5dab-437d-ab0f-73677c521373\",\"type\":\"WheelZoomTool\"},{\"id\":\"8ea85a7d-7274-4c89-8f4e-6dc78eae0482\",\"type\":\"BoxZoomTool\"},{\"id\":\"76d0a161-5e08-482f-aa43-e9ffc7662884\",\"type\":\"SaveTool\"},{\"id\":\"16d3b1da-4a76-452d-a869-076abc7cce62\",\"type\":\"ResetTool\"},{\"id\":\"638b93bf-c19b-43ac-ae1b-43bd256611d6\",\"type\":\"HelpTool\"}]},\"id\":\"c0178539-8dc8-4d00-bf05-8d5e10b86d9d\",\"type\":\"Toolbar\"},{\"attributes\":{\"bottom\":{\"value\":0},\"fill_alpha\":{\"value\":0.1},\"fill_color\":{\"value\":\"#1f77b4\"},\"left\":{\"field\":\"left\"},\"line_alpha\":{\"value\":0.1},\"line_color\":{\"value\":\"#1f77b4\"},\"right\":{\"field\":\"right\"},\"top\":{\"field\":\"top\"}},\"id\":\"972782eb-8c02-4ac0-921d-96ebf84c4bc5\",\"type\":\"Quad\"},{\"attributes\":{},\"id\":\"49a504ec-3f54-4d67-b3d7-940a042819b6\",\"type\":\"LinearScale\"},{\"attributes\":{\"callback\":null},\"id\":\"691a968d-38f2-4326-a845-6672da7c84a6\",\"type\":\"DataRange1d\"},{\"attributes\":{\"label\":{\"value\":\"Call Rate\"},\"renderers\":[{\"id\":\"7e7e8f3c-3087-474b-bcfc-64fbfb98f6e3\",\"type\":\"GlyphRenderer\"}]},\"id\":\"81c125b2-9b32-47ed-aeef-9f0878666295\",\"type\":\"LegendItem\"}],\"root_ids\":[\"067941eb-9019-4ca4-960d-bd1aa37d3872\"]},\"title\":\"Bokeh Application\",\"version\":\"0.13.0\"}};\n",
       "  var render_items = [{\"docid\":\"6f90481d-9ef1-4c0d-ba77-cb1cf9e67148\",\"roots\":{\"067941eb-9019-4ca4-960d-bd1aa37d3872\":\"0f5c6d6e-0045-4c2a-806d-737c17b2157d\"}}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    var attempts = 0;\n",
       "    var timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        embed_document(root);\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "      attempts++;\n",
       "      if (attempts > 100) {\n",
       "        console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\")\n",
       "        clearInterval(timer);\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "067941eb-9019-4ca4-960d-bd1aa37d3872"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = hl.plot.histogram(mt.sample_qc.call_rate, range=(.88,1), legend='Call Rate')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GWAS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After filtering for only 2 SupPopulations, 427/1092 samples remaining\n",
      "The SuperPopulations are:  ['AFR' 'AMR']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/sklearn/preprocessing/label.py:151: DeprecationWarning:\n",
      "\n",
      "The truth value of an empty array is ambiguous. Returning False, but in future this will result in an error. Use `array.size > 0` to check that an array is not empty.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Filter for only 2 sub populations then do logistic regression\n",
    "mt_filt1 = mt.filter_cols(((mt.col.Super_Population_Classification == 0) | (mt.col.Super_Population_Classification == 1)))\n",
    "\n",
    "print(\"After filtering for only 2 SupPopulations, %d/1092 samples remaining\" % mt_filt1.count_cols())\n",
    "print (\"The SuperPopulations are: \", le.inverse_transform([0, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example above considers a model of the form\n",
    "\t``` Prob(population) = sigmoid(β0+β1gt+ε),ε∼N(0,σ2) ``` \n",
    "where sigmoid is the sigmoid function, the genotype gt is coded as 0 for HomRef, 1 for Het, and 2 for HomVar, and the covariate is.female is coded as for 1 for True (female) and 0 for False (male). The null model sets β1=0.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 13:01:18 Hail: INFO: logistic_regression: running wald on 427 samples for response variable y,\n",
      "    with input variable x, intercept, and 0 additional covariates...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------------\n",
      "Type:\n",
      "    struct {\n",
      "        locus: locus<GRCh37>, \n",
      "        alleles: array<str>, \n",
      "        rsid: str, \n",
      "        qual: float64, \n",
      "        filters: set<str>, \n",
      "        info: struct {\n",
      "            LDAF: float64, \n",
      "            AVGPOST: float64, \n",
      "            RSQ: float64, \n",
      "            ERATE: float64, \n",
      "            THETA: float64, \n",
      "            CIEND: array<int32>, \n",
      "            CIPOS: array<int32>, \n",
      "            END: int32, \n",
      "            HOMLEN: array<int32>, \n",
      "            HOMSEQ: array<str>, \n",
      "            SVLEN: int32, \n",
      "            SVTYPE: str, \n",
      "            AC: array<int32>, \n",
      "            AN: int32, \n",
      "            AA: str, \n",
      "            AF: array<float64>, \n",
      "            AMR_AF: float64, \n",
      "            ASN_AF: float64, \n",
      "            AFR_AF: float64, \n",
      "            EUR_AF: float64, \n",
      "            VT: str, \n",
      "            SNPSOURCE: array<str>\n",
      "        }, \n",
      "        variant_qc: struct {\n",
      "            AC: array<int32>, \n",
      "            AF: array<float64>, \n",
      "            AN: int32, \n",
      "            homozygote_count: array<int32>, \n",
      "            n_called: int64, \n",
      "            n_not_called: int64, \n",
      "            call_rate: float32, \n",
      "            n_het: int64, \n",
      "            n_non_ref: int64, \n",
      "            het_freq_hwe: float64, \n",
      "            p_value_hwe: float64\n",
      "        }, \n",
      "        logreg: struct {\n",
      "            beta: float64, \n",
      "            standard_error: float64, \n",
      "            z_stat: float64, \n",
      "            p_value: float64, \n",
      "            fit: struct {\n",
      "                n_iterations: int32, \n",
      "                converged: bool, \n",
      "                exploded: bool\n",
      "            }\n",
      "        }\n",
      "    }\n",
      "--------------------------------------------------------\n",
      "Source:\n",
      "    <hail.matrixtable.MatrixTable object at 0x7fd957610c18>\n",
      "Index:\n",
      "    ['row']\n",
      "--------------------------------------------------------\n",
      "----------------------------------------\n",
      "Global fields:\n",
      "    None\n",
      "----------------------------------------\n",
      "Column fields:\n",
      "    's': str \n",
      "    'Family_ID': str \n",
      "    'Population': str \n",
      "    'Population_Description': str \n",
      "    'Gender': str \n",
      "    'Relationship': str \n",
      "    'Unexpected_Parent_Child': str \n",
      "    'Non_Paternity': str \n",
      "    'Siblings': str \n",
      "    'Grandparents': str \n",
      "    'Avuncular': str \n",
      "    'Half_Siblings': str \n",
      "    'Unknown_Second_Order': str \n",
      "    'Third_Order': str \n",
      "    'In_Low_Coverage_Pilot': str \n",
      "    'LC_Pilot_Platforms': str \n",
      "    'LC_Pilot_Centers': str \n",
      "    'In_High_Coverage_Pilot': str \n",
      "    'HC_Pilot_Platforms': str \n",
      "    'HC_Pilot_Centers': str \n",
      "    'In_Exon_Targetted_Pilot': str \n",
      "    'ET_Pilot_Platforms': str \n",
      "    'ET_Pilot_Centers': str \n",
      "    'Has_Sequence_in_Phase1': str \n",
      "    'Phase1_LC_Platform': str \n",
      "    'Phase1_LC_Centers': str \n",
      "    'Phase1_E_Platform': str \n",
      "    'Phase1_E_Centers': str \n",
      "    'In_Phase1_Integrated_Variant_Set': str \n",
      "    'Has_Phase1_chrY_SNPS': str \n",
      "    'Has_phase1_chrY_Deletions': str \n",
      "    'Has_phase1_chrMT_SNPs': str \n",
      "    'Main_project_LC_Centers': str \n",
      "    'Main_project_LC_platform': str \n",
      "    'Total_LC_Sequence': str \n",
      "    'LC_Non_Duplicated_Aligned_Coverage': str \n",
      "    'Main_Project_E_Centers': str \n",
      "    'Main_Project_E_Platform': str \n",
      "    'Total_Exome_Sequence': str \n",
      "    'X_Targets_Covered_to_20x_or_greater': str \n",
      "    'VerifyBam_E_Omni_Free': str \n",
      "    'VerifyBam_E_Affy_Free': str \n",
      "    'VerifyBam_E_Omni_Chip': str \n",
      "    'VerifyBam_E_Affy_Chip': str \n",
      "    'VerifyBam_LC_Omni_Free': str \n",
      "    'VerifyBam_LC_Affy_Free': str \n",
      "    'VerifyBam_LC_Omni_Chip': str \n",
      "    'VerifyBam_LC_Affy_Chip': str \n",
      "    'LC_Indel_Ratio': str \n",
      "    'E_Indel_Ratio': str \n",
      "    'LC_Passed_QC': str \n",
      "    'E_Passed_QC': str \n",
      "    'In_Final_Phase_Variant_Calling': str \n",
      "    'Has_Omni_Genotypes': str \n",
      "    'Has_Axiom_Genotypes': str \n",
      "    'Has_Affy_6_0_Genotypes': str \n",
      "    'Has_Exome_LOF_Genotypes': str \n",
      "    'EBV_Coverage': str \n",
      "    'DNA_Source_from_Coriell': str \n",
      "    'Has_Sequence_from_Blood_in_Index': str \n",
      "    'Super_Population': str \n",
      "    'Super_Population_Description': str \n",
      "    'Super_Population_Classification': float64 \n",
      "    'sample_qc': struct {\n",
      "        call_rate: float64, \n",
      "        n_called: int64, \n",
      "        n_not_called: int64, \n",
      "        n_hom_ref: int64, \n",
      "        n_het: int64, \n",
      "        n_hom_var: int64, \n",
      "        n_non_ref: int64, \n",
      "        n_singleton: int64, \n",
      "        n_snp: int64, \n",
      "        n_insertion: int64, \n",
      "        n_deletion: int64, \n",
      "        n_transition: int64, \n",
      "        n_transversion: int64, \n",
      "        n_star: int64, \n",
      "        r_ti_tv: float64, \n",
      "        r_het_hom_var: float64, \n",
      "        r_insertion_deletion: float64\n",
      "    } \n",
      "----------------------------------------\n",
      "Row fields:\n",
      "    'locus': locus<GRCh37> \n",
      "    'alleles': array<str> \n",
      "    'rsid': str \n",
      "    'qual': float64 \n",
      "    'filters': set<str> \n",
      "    'info': struct {\n",
      "        LDAF: float64, \n",
      "        AVGPOST: float64, \n",
      "        RSQ: float64, \n",
      "        ERATE: float64, \n",
      "        THETA: float64, \n",
      "        CIEND: array<int32>, \n",
      "        CIPOS: array<int32>, \n",
      "        END: int32, \n",
      "        HOMLEN: array<int32>, \n",
      "        HOMSEQ: array<str>, \n",
      "        SVLEN: int32, \n",
      "        SVTYPE: str, \n",
      "        AC: array<int32>, \n",
      "        AN: int32, \n",
      "        AA: str, \n",
      "        AF: array<float64>, \n",
      "        AMR_AF: float64, \n",
      "        ASN_AF: float64, \n",
      "        AFR_AF: float64, \n",
      "        EUR_AF: float64, \n",
      "        VT: str, \n",
      "        SNPSOURCE: array<str>\n",
      "    } \n",
      "    'variant_qc': struct {\n",
      "        AC: array<int32>, \n",
      "        AF: array<float64>, \n",
      "        AN: int32, \n",
      "        homozygote_count: array<int32>, \n",
      "        n_called: int64, \n",
      "        n_not_called: int64, \n",
      "        call_rate: float32, \n",
      "        n_het: int64, \n",
      "        n_non_ref: int64, \n",
      "        het_freq_hwe: float64, \n",
      "        p_value_hwe: float64\n",
      "    } \n",
      "    'logreg': struct {\n",
      "        beta: float64, \n",
      "        standard_error: float64, \n",
      "        z_stat: float64, \n",
      "        p_value: float64, \n",
      "        fit: struct {\n",
      "            n_iterations: int32, \n",
      "            converged: bool, \n",
      "            exploded: bool\n",
      "        }\n",
      "    } \n",
      "----------------------------------------\n",
      "Entry fields:\n",
      "    'GT': call \n",
      "    'DS': float64 \n",
      "    'GL': array<float64> \n",
      "    'PL': array<int32> \n",
      "----------------------------------------\n",
      "Column key: ['s']\n",
      "Row key: ['locus', 'alleles']\n",
      "Partition key: ['locus']\n",
      "----------------------------------------\n"
     ]
    }
   ],
   "source": [
    "#Logistic Regression\n",
    "gwas_logistic = hl.logistic_regression(\n",
    "           test='wald',\n",
    "           y=mt_filt1.col.Super_Population_Classification,\n",
    "           x=mt_filt1.GT.n_alt_alleles()\n",
    "       )\n",
    "gwas_logistic.row.describe()\n",
    "gwas_logistic.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qqplot(pvals):\n",
    "    spvals = sorted(filter(lambda x: x and not (isnan(x)), pvals))\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",
    "    p = figure(title=\"Q Q Plot\",\n",
    "              x_axis_label = 'Expected P-Value (-log10 scale)',\n",
    "              y_axis_label = 'Observe P-Value (-log10 scale)')\n",
    "    p.scatter(x=exp, y=obs, color = 'black')\n",
    "    bound = max(max(exp), max(obs)) * 1.1\n",
    "    p.line([0, bound], [0, bound], color = 'red')\n",
    "    show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvals = gwas_logistic.logreg.p_value.collect()\n",
    "qqplot(pvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-08-07 13:31:15 Hail: INFO: linear_regression: running on 427 samples for 1 response variable y,\n",
      "    with input variable x, intercept, and 0 additional covariates...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------------\n",
      "Type:\n",
      "    struct {\n",
      "        locus: locus<GRCh37>, \n",
      "        alleles: array<str>, \n",
      "        rsid: str, \n",
      "        qual: float64, \n",
      "        filters: set<str>, \n",
      "        info: struct {\n",
      "            LDAF: float64, \n",
      "            AVGPOST: float64, \n",
      "            RSQ: float64, \n",
      "            ERATE: float64, \n",
      "            THETA: float64, \n",
      "            CIEND: array<int32>, \n",
      "            CIPOS: array<int32>, \n",
      "            END: int32, \n",
      "            HOMLEN: array<int32>, \n",
      "            HOMSEQ: array<str>, \n",
      "            SVLEN: int32, \n",
      "            SVTYPE: str, \n",
      "            AC: array<int32>, \n",
      "            AN: int32, \n",
      "            AA: str, \n",
      "            AF: array<float64>, \n",
      "            AMR_AF: float64, \n",
      "            ASN_AF: float64, \n",
      "            AFR_AF: float64, \n",
      "            EUR_AF: float64, \n",
      "            VT: str, \n",
      "            SNPSOURCE: array<str>\n",
      "        }, \n",
      "        variant_qc: struct {\n",
      "            AC: array<int32>, \n",
      "            AF: array<float64>, \n",
      "            AN: int32, \n",
      "            homozygote_count: array<int32>, \n",
      "            n_called: int64, \n",
      "            n_not_called: int64, \n",
      "            call_rate: float32, \n",
      "            n_het: int64, \n",
      "            n_non_ref: int64, \n",
      "            het_freq_hwe: float64, \n",
      "            p_value_hwe: float64\n",
      "        }, \n",
      "        linreg: struct {\n",
      "            n: int32, \n",
      "            sum_x: float64, \n",
      "            y_transpose_x: float64, \n",
      "            beta: float64, \n",
      "            standard_error: float64, \n",
      "            t_stat: float64, \n",
      "            p_value: float64\n",
      "        }\n",
      "    }\n",
      "--------------------------------------------------------\n",
      "Source:\n",
      "    <hail.matrixtable.MatrixTable object at 0x7fd95c14c160>\n",
      "Index:\n",
      "    ['row']\n",
      "--------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "#Linear Regression\n",
    "gwas_linear = hl.linear_regression(\n",
    "           test='wald',\n",
    "           y=mt_filt1.Super_Population_Classification,\n",
    "           x=mt_filt1.GT.n_alt_alleles(),\n",
    "       )\n",
    "gwas_linear.row.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.qq(gwas_linear.linreg.p_value)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear_regression method can also take column fields to use as covariates. We already annotated our samples with reported ancestry, but it is good to be skeptical of these labels due to human error. Genomes don’t have that problem! Instead of using reported ancestry, we will use genetic ancestry by including computed principal components in our model.\n",
    "\n",
    "The pca method produces eigenvalues as a list and sample PCs as a Table, and can also produce variant loadings when asked. The hwe_normalized_pca method does the same, using HWE-normalized genotypes for the PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_eigenvalues, pca_scores, _ = hl.hwe_normalized_pca(mt_filt1.GT, k=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = hl.plot.scatter(pca_scores.scores[0], pca_scores.scores[1],\n",
    "                    label=common_mt.cols()[pca_scores.s].SuperPopulation,\n",
    "                    title='PCA', xlabel='PC1', ylabel='PC2')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can rerun our linear regression, controlling by the first few principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#New table for pca data \n",
    "pcatable = pca_score.select(Super_Population = mt_filt1.cols()[pca_score.s].population,\n",
    "                             PC1 = pca_score.scores[0],\n",
    "                             PC2 = pca_score.scores[1]\n",
    "                           )\n",
    "\n",
    "\n",
    "\n",
    "pca_df = hl.Table.to_pandas(pcatable)\n",
    "mt_filt1 = mt_filt1.annotate_cols(pca = pca_score[mt_filt1.s])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas_logistic2 = hl.logistic_regression(\n",
    "           test='lrt',\n",
    "           y=mt_filt1.col.Super_Population_Classification,\n",
    "           x=mt_filt1.GT.n_alt_alleles(),\n",
    "           covariates=[mt_filt1.pca.scores[0], mt_filt1.pca.scores[1], mt_filt1.pca.scores[2],mt_filt1.col.Super_Population_Classification])\n",
    "       \n",
    "pvals = gwas_logistic2.logreg.p_value.collect()\n",
    "qqplot(pvals)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas_linear2 = hl.linear_regression(\n",
    "           test='wald',\n",
    "           y=mt_filt1.Super_Population_Classification,\n",
    "           x=mt_filt1.GT.n_alt_alleles(),\n",
    "           covariates=[mt_filt1.pca.scores[0], mt_filt1.pca.scores[1], mt_filt1.pca.scores[2],mt_filt1.col.Super_Population_Classification])\n",
    "\n",
    "\n",
    "p = hl.plot.qq(gwas_linear2.linreg.p_value)\n",
    "show(p)\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## MANHATTAN PLOT\n",
    "from bokeh.plotting import figure, output_file, show\n",
    "\n",
    "filtered = gwas_linear2.filter_rows(gwas_linear2.row.logreg.p_value < .005, keep = True)\n",
    "\n",
    "#filtered.count()\n",
    "x = hl.plot.manhattan(filtered.row.logreg.p_value, locus=filtered.row.locus, title=\"Practice\", size=4, hover_fields=None)\n",
    "show(x)"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}