tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro \n",
    "\n",
    "In this notebook, we will learn how to use hail functions to annotate the Expression Modifier Score (EMS) for a given list of variant-gene pairs (vgs).\n",
    "\n",
    "Most of the work in this notebook is performed in google cloud using hail. Specifically, we ran\n",
    "\n",
    "`hailctl dataproc start qbwcluster-mpl -w 8 --pkgs matplotlib,seaborn,sklearn`\n",
    "\n",
    "Followed by\n",
    "\n",
    "`hailctl dataproc connect qbwcluster-mpl notebook`\n",
    "\n",
    "To start the jupyter notebook in the google cloud cluster.\n",
    "\n",
    "(For more info on setting up the cluster and notebook, please visit the relevant pages: [hail](https://hail.is/index.html) and the [hail dataproc](https://github.com/Nealelab/cloudtools/blob/master/README.md))\n",
    "\n",
    "After building the cluster and starting this notebook, the first step is to import hail and relevant tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Running on Apache Spark version 2.4.5\n",
      "SparkUI available at http://qbwcluster-mpl-m.c.encode-uk-biobank.internal:4041\n",
      "Welcome to\n",
      "     __  __     <>__\n",
      "    / /_/ /__  __/ /\n",
      "   / __  / _ `/ / /\n",
      "  /_/ /_/\\_,_/_/_/   version 0.2.26-2dcc3d963867\n",
      "LOGGING: writing to /home/hail/hail-20200914-2112-0.2.26-2dcc3d963867.log\n"
     ]
    }
   ],
   "source": [
    "#import and initialize hail \n",
    "import hail as hl\n",
    "import hail.expr.aggregators as agg\n",
    "hl.init(default_reference='GRCh38')\n",
    "\n",
    "#import other useful tools\n",
    "import time as tm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now assuming that we have the environment built, first, let's start taking a look at the EMS in Whole Blood by reading it as a hail table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>rf_score_raw</td><td>ems</td><td>ems_normalized</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>float64</td><td>float64</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000014919.12&quot;</td><td>3.90e-02</td><td>1.78e-06</td><td>1.47e-01</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000023839.10&quot;</td><td>4.48e-02</td><td>2.42e-06</td><td>1.99e-01</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>1.77e-02</td><td>5.31e-07</td><td>4.37e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000075290.7&quot;</td><td>2.07e-02</td><td>5.08e-07</td><td>4.19e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000075826.16&quot;</td><td>1.84e-02</td><td>5.31e-07</td><td>4.37e-02</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+----------+\n",
       "| v                         | g                    | rf_score_raw |      ems |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| str                       | str                  |      float64 |  float64 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000014919.12\" |     3.90e-02 | 1.78e-06 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000023839.10\" |     4.48e-02 | 2.42e-06 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000055950.16\" |     1.77e-02 | 5.31e-07 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000075290.7\"  |     2.07e-02 | 5.08e-07 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000075826.16\" |     1.84e-02 | 5.31e-07 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "\n",
       "+----------------+\n",
       "| ems_normalized |\n",
       "+----------------+\n",
       "|        float64 |\n",
       "+----------------+\n",
       "|       1.47e-01 |\n",
       "|       1.99e-01 |\n",
       "|       4.37e-02 |\n",
       "|       4.19e-02 |\n",
       "|       4.37e-02 |\n",
       "+----------------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tissue_name = \"Whole_Blood\"\n",
    "ems = hl.read_table(\"gs://expression-modifier-score/public/ht/ems_{0}.ht\".format(tissue_name))\n",
    "ems.head(5).show() #show the first 5 entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `v` column shows the variant (in hg38), \n",
    "\n",
    "the `g` column shows the gene, \n",
    "\n",
    "the `rf_score_raw` is the raw output score from the random forest model used in building the EMS (ranging from 0 to 1), \n",
    "\n",
    "the `ems` is the EMS, corresponding to the estimated probability of the variant-gene pairs being a putative causal eQTL, and\n",
    "\n",
    "the `ems_normalized` is the EMS divided by the probability that a random variant-gene pair being a putative causal eQTL."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since EMS is defined for all the variant-gene pairs in GTEx, the table is quite large:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "149429947"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ems.count() #counts the number of rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this number is quite large, in this tutorial we will cut down and focus on the first 100,000 entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ems = ems.head(100000) #cut down to the first 100,000 entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EMS in different tissues\n",
    "\n",
    "We can simply change the tissue name to get the EMS for a different tissue:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>rf_score_raw</td><td>ems</td><td>ems_normalized</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>float64</td><td>float64</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000014919.12&quot;</td><td>2.57e-02</td><td>9.74e-07</td><td>8.97e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000023839.10&quot;</td><td>5.29e-02</td><td>2.22e-06</td><td>2.05e-01</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>7.14e-03</td><td>3.53e-07</td><td>3.25e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000075290.7&quot;</td><td>1.00e-02</td><td>4.27e-07</td><td>3.93e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000075826.16&quot;</td><td>8.57e-03</td><td>4.12e-07</td><td>3.80e-02</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+----------+\n",
       "| v                         | g                    | rf_score_raw |      ems |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| str                       | str                  |      float64 |  float64 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000014919.12\" |     2.57e-02 | 9.74e-07 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000023839.10\" |     5.29e-02 | 2.22e-06 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000055950.16\" |     7.14e-03 | 3.53e-07 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000075290.7\"  |     1.00e-02 | 4.27e-07 |\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000075826.16\" |     8.57e-03 | 4.12e-07 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "\n",
       "+----------------+\n",
       "| ems_normalized |\n",
       "+----------------+\n",
       "|        float64 |\n",
       "+----------------+\n",
       "|       8.97e-02 |\n",
       "|       2.05e-01 |\n",
       "|       3.25e-02 |\n",
       "|       3.93e-02 |\n",
       "|       3.80e-02 |\n",
       "+----------------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tissue_name = \"Muscle_Skeletal\"\n",
    "ems_muscle = hl.read_table(\"gs://expression-modifier-score/public/ht/ems_{0}.ht\".format(tissue_name))\n",
    "ems_muscle.head(5).show() #show the first 5 entries for the EMS in skeletal muscle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering\n",
    "Here we will show an example of filtering the EMS table with specific condition of interest. \n",
    "\n",
    "One might be interested in only the variants that are more likely to be putative causal than by chance. \n",
    "Then we can filter the table with the value of `ems_normalized`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>rf_score_raw</td><td>ems</td><td>ems_normalized</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>float64</td><td>float64</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>5.91e-01</td><td>1.21e-04</td><td>9.93e+00</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>4.25e-01</td><td>4.71e-05</td><td>3.88e+00</td></tr>\n",
       "<tr><td>&quot;chr10_100000979_T_C_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>6.52e-01</td><td>2.27e-04</td><td>1.87e+01</td></tr>\n",
       "<tr><td>&quot;chr10_100002012_T_C_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>6.38e-01</td><td>1.80e-04</td><td>1.49e+01</td></tr>\n",
       "<tr><td>&quot;chr10_100002038_G_A_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>4.85e-01</td><td>8.12e-05</td><td>6.69e+00</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+----------+\n",
       "| v                         | g                    | rf_score_raw |      ems |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| str                       | str                  |      float64 |  float64 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000107554.16\" |     5.91e-01 | 1.21e-04 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000107554.16\" |     4.25e-01 | 4.71e-05 |\n",
       "| \"chr10_100000979_T_C_b38\" | \"ENSG00000107554.16\" |     6.52e-01 | 2.27e-04 |\n",
       "| \"chr10_100002012_T_C_b38\" | \"ENSG00000107554.16\" |     6.38e-01 | 1.80e-04 |\n",
       "| \"chr10_100002038_G_A_b38\" | \"ENSG00000107554.16\" |     4.85e-01 | 8.12e-05 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "\n",
       "+----------------+\n",
       "| ems_normalized |\n",
       "+----------------+\n",
       "|        float64 |\n",
       "+----------------+\n",
       "|       9.93e+00 |\n",
       "|       3.88e+00 |\n",
       "|       1.87e+01 |\n",
       "|       1.49e+01 |\n",
       "|       6.69e+00 |\n",
       "+----------------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ems_sub = ems.filter(ems.ems_normalized>1)\n",
    "ems_sub.head(5).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000\n",
      "8548\n"
     ]
    }
   ],
   "source": [
    "print (ems.count())\n",
    "print (ems_sub.count())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only ~=8% of the vgs have `ems_normalized`>1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aggregation\n",
    "One might be interested in only the variants. We can aggregate the table to get the maximum EMS for a variant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:12:54 Hail: INFO: Ordering unsorted dataset with network shuffle\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>max_ems</td></tr>\n",
       "<tr><td>str</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>1.21e-04</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>4.71e-05</td></tr>\n",
       "<tr><td>&quot;chr10_100000979_T_C_b38&quot;</td><td>2.27e-04</td></tr>\n",
       "<tr><td>&quot;chr10_100002012_T_C_b38&quot;</td><td>1.80e-04</td></tr>\n",
       "<tr><td>&quot;chr10_100002038_G_A_b38&quot;</td><td>8.12e-05</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------+\n",
       "| v                         |  max_ems |\n",
       "+---------------------------+----------+\n",
       "| str                       |  float64 |\n",
       "+---------------------------+----------+\n",
       "| \"chr10_100000235_C_T_b38\" | 1.21e-04 |\n",
       "| \"chr10_100000943_G_A_b38\" | 4.71e-05 |\n",
       "| \"chr10_100000979_T_C_b38\" | 2.27e-04 |\n",
       "| \"chr10_100002012_T_C_b38\" | 1.80e-04 |\n",
       "| \"chr10_100002038_G_A_b38\" | 8.12e-05 |\n",
       "+---------------------------+----------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ems_max = ems.group_by(\"v\").aggregate(max_ems = hl.agg.max(ems.ems))\n",
    "ems_max.head(5).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Annotating EMS for a set of vgs\n",
    "Here we will show an example of annotating EMS for a set of variant and genes of interest.\n",
    "\n",
    "\n",
    "Let's assume we have a `tsv` file with 10000 rows of variant and genes, with the TSS distance annotated. \n",
    "\n",
    "We can read that as a hail table by `import_table`:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:12:55 Hail: INFO: Reading table to impute column types\n",
      "2020-09-14 21:12:56 Hail: INFO: Finished type imputation\n",
      "  Loading column 'v' as type 'str' (imputed)\n",
      "  Loading column 'g' as type 'str' (imputed)\n",
      "  Loading column 'tss_distance' as type 'int32' (imputed)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>tss_distance</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>int32</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>-9684</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986423</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000107816.17&quot;</td><td>-995675</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000119906.11&quot;</td><td>-911797</td></tr>\n",
       "<tr><td>&quot;chr10_100000979_T_C_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986387</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+\n",
       "| v                         | g                    | tss_distance |\n",
       "+---------------------------+----------------------+--------------+\n",
       "| str                       | str                  |        int32 |\n",
       "+---------------------------+----------------------+--------------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000107554.16\" |        -9684 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000055950.16\" |      -986423 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000107816.17\" |      -995675 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000119906.11\" |      -911797 |\n",
       "| \"chr10_100000979_T_C_b38\" | \"ENSG00000055950.16\" |      -986387 |\n",
       "+---------------------------+----------------------+--------------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = hl.import_table(\"gs://qingbowang/vg_head10000_example.tsv\", impute=True)\n",
    "df.head(5).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can annotate the EMS using hail language. Below are two examples of doing so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<StructExpression of type struct{v: str, g: str}>\n",
      "<StructExpression of type struct{v: str, g: str}>\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:12:57 Hail: INFO: Coerced sorted dataset\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>tss_distance</td><td>ems</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>int32</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>-9684</td><td>1.21e-04</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986423</td><td>2.74e-07</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000107816.17&quot;</td><td>-995675</td><td>2.74e-07</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000119906.11&quot;</td><td>-911797</td><td>2.74e-07</td></tr>\n",
       "<tr><td>&quot;chr10_100000979_T_C_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986387</td><td>5.16e-07</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+----------+\n",
       "| v                         | g                    | tss_distance |      ems |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| str                       | str                  |        int32 |  float64 |\n",
       "+---------------------------+----------------------+--------------+----------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000107554.16\" |        -9684 | 1.21e-04 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000055950.16\" |      -986423 | 2.74e-07 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000107816.17\" |      -995675 | 2.74e-07 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000119906.11\" |      -911797 | 2.74e-07 |\n",
       "| \"chr10_100000979_T_C_b38\" | \"ENSG00000055950.16\" |      -986387 | 5.16e-07 |\n",
       "+---------------------------+----------------------+--------------+----------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:12:59 Hail: INFO: Coerced sorted dataset\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>v</td><td>g</td><td>tss_distance</td><td>rf_score_raw</td><td>ems</td><td>ems_normalized</td></tr>\n",
       "<tr><td>str</td><td>str</td><td>int32</td><td>float64</td><td>float64</td><td>float64</td></tr>\n",
       "</thead><tbody><tr><td>&quot;chr10_100000235_C_T_b38&quot;</td><td>&quot;ENSG00000107554.16&quot;</td><td>-9684</td><td>5.91e-01</td><td>1.21e-04</td><td>9.93e+00</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986423</td><td>8.61e-03</td><td>2.74e-07</td><td>2.26e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000107816.17&quot;</td><td>-995675</td><td>8.61e-03</td><td>2.74e-07</td><td>2.26e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000943_G_A_b38&quot;</td><td>&quot;ENSG00000119906.11&quot;</td><td>-911797</td><td>9.38e-03</td><td>2.74e-07</td><td>2.26e-02</td></tr>\n",
       "<tr><td>&quot;chr10_100000979_T_C_b38&quot;</td><td>&quot;ENSG00000055950.16&quot;</td><td>-986387</td><td>1.87e-02</td><td>5.16e-07</td><td>4.25e-02</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+---------------------------+----------------------+--------------+\n",
       "| v                         | g                    | tss_distance |\n",
       "+---------------------------+----------------------+--------------+\n",
       "| str                       | str                  |        int32 |\n",
       "+---------------------------+----------------------+--------------+\n",
       "| \"chr10_100000235_C_T_b38\" | \"ENSG00000107554.16\" |        -9684 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000055950.16\" |      -986423 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000107816.17\" |      -995675 |\n",
       "| \"chr10_100000943_G_A_b38\" | \"ENSG00000119906.11\" |      -911797 |\n",
       "| \"chr10_100000979_T_C_b38\" | \"ENSG00000055950.16\" |      -986387 |\n",
       "+---------------------------+----------------------+--------------+\n",
       "\n",
       "+--------------+----------+----------------+\n",
       "| rf_score_raw |      ems | ems_normalized |\n",
       "+--------------+----------+----------------+\n",
       "|      float64 |  float64 |        float64 |\n",
       "+--------------+----------+----------------+\n",
       "|     5.91e-01 | 1.21e-04 |       9.93e+00 |\n",
       "|     8.61e-03 | 2.74e-07 |       2.26e-02 |\n",
       "|     8.61e-03 | 2.74e-07 |       2.26e-02 |\n",
       "|     9.38e-03 | 2.74e-07 |       2.26e-02 |\n",
       "|     1.87e-02 | 5.16e-07 |       4.25e-02 |\n",
       "+--------------+----------+----------------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#make sure that two tables of interest have the same keys\n",
    "df = df.key_by(\"v\", \"g\")\n",
    "print (ems.key)\n",
    "print (df.key)\n",
    "\n",
    "#1. simply `annotate` the column of interest\n",
    "df1 = df.annotate(ems = ems[df.key].ems)\n",
    "df1.head(5).show()\n",
    "\n",
    "#2. use `join` to annotate all the columns together\n",
    "df2 = df.join(ems, how=\"left\")\n",
    "df2.head(5).show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aggregation and Enrichment analysis\n",
    "\n",
    "In our manuscript, we quantified the enrichment of different genomic annotations in different EMS bins.\n",
    "\n",
    "As an example to replicate such analysis, here, we will investigate the relationship between distance to TSS and EMS.\n",
    "\n",
    "First, let's bin EMS and TSS in different windows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:13:00 Hail: INFO: Coerced sorted dataset\n",
      "2020-09-14 21:13:02 Hail: INFO: Coerced sorted dataset\n",
      "2020-09-14 21:13:02 Hail: INFO: Coerced dataset with out-of-order partitions.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><thead style=\"font-weight: bold;\"><tr><td>tss_distance_bin</td><td>ems_normalized_bin</td><td>n</td></tr>\n",
       "<tr><td>int32</td><td>int32</td><td>int64</td></tr>\n",
       "</thead><tbody><tr><td>1</td><td>1</td><td>146</td></tr>\n",
       "<tr><td>1</td><td>2</td><td>5389</td></tr>\n",
       "<tr><td>2</td><td>1</td><td>78</td></tr>\n",
       "<tr><td>2</td><td>2</td><td>4387</td></tr>\n",
       "</tbody></table>"
      ],
      "text/plain": [
       "+------------------+--------------------+-------+\n",
       "| tss_distance_bin | ems_normalized_bin |     n |\n",
       "+------------------+--------------------+-------+\n",
       "|            int32 |              int32 | int64 |\n",
       "+------------------+--------------------+-------+\n",
       "|                1 |                  1 |   146 |\n",
       "|                1 |                  2 |  5389 |\n",
       "|                2 |                  1 |    78 |\n",
       "|                2 |                  2 |  4387 |\n",
       "+------------------+--------------------+-------+"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df2 = df2.filter(hl.is_defined(df2.ems))#remove small number of NAs\n",
    "df2 = df2.annotate(tss_distance_bin = hl.case()\n",
    "                                      .when(df2.tss_distance <10**4, 1) \n",
    "                                      .default(2), \n",
    "                   ems_normalized_bin = hl.case()\n",
    "                                       .when(df2.ems_normalized>10, 1)\n",
    "                                       .default(2) )\n",
    "\n",
    "summary = df2.group_by(\"tss_distance_bin\", \"ems_normalized_bin\").aggregate(n = hl.agg.count())\n",
    "summary.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make small pandas dataframe to compare the distribution of distance to TSS for vgs with high vs low EMS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-09-14 21:13:02 Hail: INFO: Coerced sorted dataset\n",
      "2020-09-14 21:13:03 Hail: INFO: Coerced sorted dataset\n",
      "2020-09-14 21:13:03 Hail: INFO: Coerced dataset with out-of-order partitions.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>&gt;10</th>\n",
       "      <th>&lt;10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>&lt;10000</th>\n",
       "      <td>146</td>\n",
       "      <td>5389</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&gt;10000</th>\n",
       "      <td>78</td>\n",
       "      <td>4387</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        >10   <10\n",
       "<10000  146  5389\n",
       ">10000   78  4387"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = summary.to_pandas()\n",
    "summary = pd.pivot_table(summary, values='n', index=\"tss_distance_bin\",columns=\"ems_normalized_bin\")\n",
    "summary.index.names = [\"tss_distance\"]\n",
    "summary.columns.names = [\"ems_normalized\"]\n",
    "summary.index = [\"<10000\", \">10000\"]\n",
    "summary.columns = [\">10\", \"<10\"]\n",
    "summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And not surprisingly, variants with distance to TSS <10kb are enriched in high (>10x the random) EMS variants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEOCAYAAAApCQiDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlAUlEQVR4nO3deZhcVZ3/8feHsCQYMAIhLBqD6LAjkRZli4EYcZlxEEdRQQioiKKgCPwAAQkKkdVlZAsjBESQAdmEKHtQJAQ6BoFhU4EEE0KCEBHIRvj+/ji3oLqorj5d3VVd6f68nuc+XXXuPeeezPjUl3PvOd+jiMDMzKwVrNLXHTAzMytxUDIzs5bhoGRmZi3DQcnMzFqGg5KZmbWMVfu6Ayu79dZbL0aNGtXX3TAzW2nMnDnzuYgYXu2cg1IPjRo1ivb29r7uhpnZSkPS7M7O+fGdmZm1DAclMzNrGQ5KZmbWMhyUzMysZTgomZlZy3BQMjOzluGgZGZmLcNByczMWoaDUh/Z+/zp7H3+9L7uhplZS3FQMjOzltG0oCTpfZIukdReHL+Q9L5m3d/MzFpfU4KSpH2A+4ANganFMQK4V9K+zeiDmZm1vmYlZD0ZOD4iTikvlHQM8APg0ib1w8zMWlizHt8NB/63SvmVwPpN6oOZmbW4ZgWlO4CxVcrHAnc2qQ9mZtbiGvb4TtJeZV9/C0yS1AbcU5R9ENgLOLFRfTAzs5VLI98pXVWl7KDiKPffwDkN7IeZma0kGhaUIsJroMzMrFuaNSV80xrnxjWjD2Zm1vqaNZq5WdKIykJJHwaubVIfzMysxTUrKP0OuEXSW0sFRUC6BjiiSX0wM7MW16yg9A3g/4DfShoiaTxwHfCdiDi/SX0wM7MW15SgFBEBfBH4J2nN0jXAtyJicjPub2ZmK4dGrlOqlmx1Eiml0CXAzNI1EfGnRvXDzMxWHo1cp9QOBKCystL3g4GvFp8DGNTAfpiZ2UqikUFpkwa2bWZm/VAjF8/OblTbZmbWPzV86wpJbwG+AOwEbFAUzwf+CFweES83ug9mZrZyaOjsO0lbAo8DZ5K2r5hXHMOLsseKa8zMzBo+UjobuAvYPyKWlJ+QNBiYUlyzW4P7YWZmK4FGB6UPAG2VAQkgIpZI+gFwb4P7YGZmK4m6Ht8VWRk+LOmdXVz6AvBvNc6/p7jGzMwsLyhJmiLp68Xn1Umjm5tJ74Q+VqPqBcDFko6WtL2ktxfH9pKOBi4EnGbIzMyA/JHSHryxY+wngbVIM+lOpMbOsRFxInAKcChwHzC7OO4ryiZFxEm5nZU0RtL1kuZKCkkTurh+rKTrJD0j6RVJD0g6sMp1H5I0U9ISSU9IOji3T2Zm1ntyg9LbgAXF548Cv46IBcCvgJqz5yLi1IjYCNgU2KU4No2IjSLitG72dyjwEHAYsDjj+p2AB4H/ArYGzgUmS/pC6QJJmwBTgbuB0aRUSP8t6dPd7JuZmfVQ7kSH+cDWkp4hjZpKW5oPBZbnNBARTwJPdruHHduYSgogSJqScf0pFUXnStoN+DRwWVF2MDAvIr5ZfH9E0gdIW2r8uif9NTOz7skdKV0IXEEapawAbivKPwA8Wu/NJb1D0oX11q/T2nScXLEj6f1YuZuANkmrVWtA0kGS2iW1L1y4sEHdNDMbeLKCUvHe50BgMrBLRCwrTr0KnNqD+68D7N+D+t0i6d+BcaR/R8kGwLMVlz5LGkWuV62diJgcEW0R0TZ8+PCG9NXMbCDKenwnaQxwXUS8WnHql6T3Np3V26+Lpkfm3L83SNqZ9Mju0Ijw2igzsxaU+07pDmBD3pjsUPLW4lxnW09MAV4hbU9RTVM2GZS0C+ld1AkRcW7F6fnAiIqyEaRR4HNN6J6ZmRVyg0Jp36NK6wK1EqrOA/aLiLWqHcDO3exvtxWjvN8CJ0bEj6tcMh0YX1E2HmiPiKxJHGZm1jtqjpQkXV98DOBSSUvLTg8iTbO+u0YTM4H3AVd3cr5yE8CaJA0F3l18XQUYKWk74PmImCNpErBDRIwrrh8L3AicA1wmqZSlfEVElGYonAd8Q9KPSQt5dwYmAJ/P7ZeZmfWOrh7f/aP4K9KMtfK1QctIyVYvqFH/DNK08c78le4lY20jPS4smVgcF5MCyYak9VAlE4A1SdO7jygrnw2MgjRVXdLHgR8BXyON7g6NCE8HNzNrsppBKSIOAJD0FHBGd/c+iog/dHH+ZeDObrQ3jRojq4iYUOX7hGrXVlx3J2lEZ2ZmfShrokNETGx0R8zMzHITsq4j6VxJj0taJOnF8qOLuhtKOr5I5Fpevoqk70pq2rRwMzNrbblTwn9Oygs3mfTOpbMp3tUsAA4B/o+OEx7GA98BTu9GW2Zm1o/lBqVxwPiImNHdG0TECkmXAfvRMSjtC1xZlh3CzMwGuNygtAB4qQf3uQS4R9I6EfG8pLcAe5KSu5qZmQH5i2e/C5xUrBPqtoi4n5S4de+i6FPA/IiotcbJzMwGmNyR0nGkdT0LJM2mYruKiNg2o41LSI/wzgW+WHw3MzN7XW5QuqoX7nUpMEnSh4CxvLEn04Bz7ay5zJqziGUrXmPnH97OkXtsxp6jN+7rbpmZ9bmmrVOKiAWSbgF+AUyPiNk9bXNldO2suRxz9YMsW/EaAHMXLeaYqx8EcGAyswGvKVm6y1wCvJ2UFmhAOv2mx1i8fEWHssXLV3D6TY/1UY/MzFpH7uLZ1SVNLBbPLpG0ovzoxv2uIeW6+1Vdve0H5i1a3K1yM7OBJHek9H3SDrFnAq8BRwJnkxK2fj33ZhGxPCLujIgB+wu80bAh3So3MxtIcoPSZ4GDI+J8YAVpF9pDge/x5r2IrIYj99iMIat13BNxyGqDOHKPzfqoR2ZmrSM3KI0AHi4+vwQMKz7/DvhIL/epX9tz9MZM2msbVh+U/k+/8bAhTNprG09yMDMjf0r4HGCj4u9fSZkYZgI70nGPJcuw5+iNufzeOQBc8dUd+7g3ZmatI3ekdA0p/x3AT4CJkp4EpgD/04B+mZnZAJS7TumYss9XSfo7sBPweETc0KjOmZnZwFLXOqWIuCcizuoqIEl6r6TdKsr2kfSEpAWSzqvcZ8nMzAauTkdKkvYCfhMRy4vPnYqIqzs59QPgXuCOos0tgYuK748CBwJzSVPOzcxsgKv1+O4qYAPSthW1ct8FMKiTc+8jBaaSzwEPR8QeAJIeAL6Ng5KZmVHj8V1ErBIRC8o+d3Z0FpAA1iXtVFsyBvhN2fdpgLdDNzMzIOOdkqTVJF0hadM62l8IbFy0MwjYHijfvXZ1UoYIMzOzroNSRCwnLZCNOtqfBnxP0ruA7xRld5Sd3xJ4qo52zcysH8qdfXc1UHOyQyeOB95DWnB7MnBURLxcdv6LwG11tGtmZv1QdzI6HCdpV6AdKA8sRMRZ1SpFxFOSNge2AhZGxLyKS74H/L17XTYzs/4qNyhNAF4Ati2OcgFUDUoAEfEq8OfyMkmrAoMj4s/Va5mZ2UCUm9Fhk3oalzQOWDci/res7GjgRGBVSbcCn4uIRfW0b2Zm/Uujd549mrTTLACSdgBOIW2JfhTwXuC7De6DmZmtJHIf3yHp34D/Iq0r6pAaKCIO7KTaNqTAVPIZ4O6I+ErR5tOkxbVHZvZhDHAEaWr5RsABETGlxvWDgfNIi3i3AP4YEWMrrhlLxxmBJVtExKM5/TIzs96RFZQkfQL4NTCLFBDuAzYF1gD+UKPqMFJGiJKdgall3++jWMeUaSjwEHBJcXRlELAE+Bnwcd7YB6qarYDny74v7Ea/zMysF+Q+vjsJmBgROwJLSVO5RwG3ktYideYZUvBC0hrAaGB62fm1ivayRMTUiDg2Iq4iY9FtRLwcEQdHxGS6nuW3ICLmlx0rcvtlZma9IzcobQZcUXxeDqwZEUtIwepbNer9FjhN0u7AqaSp5OUjq21Ja5haQbukZyTdVpnZ3MzMmiM3KP0LGFx8fgZ4d/F5VeBtNeqdQHp8dispI/hXImJZ2fkDgVuye9sYzwBfAz5NWiD8GHBbsSarKkkHSWqX1L5woZ/ymZn1ltyJDjOAXYCHgRuBMyW9F/gUHR/HdRARzwFjJL0VeKnKI7HPAC91u9e9KCIeIwWikumSRpEmX1R9X1Y8DpwM0NbWVk/6JTMzqyI3KB1OmmQAaY3RWqSRxePFuaokrQA2LGUbrxQRz1crbwEzSNtsmJlZE+Uunn2i7PMrpMddOVRPp1rAdqTHemZm1kS5U8KvJS14/U3FO6GmkjSUN95nrQKMlLQd8HxEzJE0CdghIsaV1dmStK5qPWBocT0RcX9x/lukTOX/V1y3L7AnaSRoZmZNlPv47hXgYmC5pF8Dv4iIOzPrflbSi7UuiIicNUcAbXRc6DqxOC4m5efbkGIKepmpwDvLvs8q/pZGcasDp5MyTywmBadPRET5eiozM2uC3Md3X5D0FtLEhi8At0h6BrgcuDQiHqpR/YfU3ospyFsIS0RMo8YjwYiYUKVsVBdtngaclnN/MzNrrOw0Q8U+SJcCl0oaDuwNHExK+1OrnXd1NtHBzMysXLcTshb55HYH9gD+DXi6xuWeLm1mZtmygpKSj0i6GHgWOBeYB4zrYluLlXX2nZmZ9YHcx3fPAGuT0gZNAG7MnIV3MWnygJmZWZdyg9LxwJXd3YwvIg7odo/MzGzAyp19d0GjO2JmZtbonWfNzMyyOSiZmVnLaHpQkjS0WIhrZmbWQadBSdLtkoYVn/crdo6tm6RDJM0B/gm8KGm2pK/3pE0zM+tfao2UdgbWLD5fBLy13ptIOpaUbujnwEeK4yLgh5KOrrddMzPrX2rNvnsUOEXSHaRFsJ0mVs1IqHowcFBEXF5WdpukvwCnkAKWmZkNcLWC0teAnwD/SUoX1Fli1ZyEqusD91UpvxcY0XU3zcxsIOj08V1E3B0R74+It5FGSu+KiLWqHGtn3OdxUnbxSl+g41bkZmY2gOVmdNgEWNiD+5wI/K+kMcAfi7KdgQ8Bn+lBu2Zm1o/kZnSYLWmEpEOALUmP7B4GzomIZzPqXy3pA8C3gX8vih8h7RI7q/OaZmY2kORuh74z8DtShvDpRfE+wLcl7RER0zutXIiImaStxs3MzKrKfXx3BmmX2YMj4jUASasA5wFnAjvVqixpBbBh5WZ/ktYFFkTEoO52fGV3xVd37OsumJm1nNygtB0woRSQACLiNUlnATmP3zrbV2kNIGcLDDMzGwByg9I/SZMdKmfKbQIs6qySpMOLjwEcLOmlstODgF1J66HMzMyyg9KvgJ9LOgq4uyjbGTiV9FivM98s/gr4MrCi7Nwy4CnSwlozM7PsoHQUKbBcWFZnOWlb9E7TBJW2Si+yQuwVES/U31UzM+vvcqeELwMOk3QMsGlR/LeIeCWz/m519s/MzAaQ3JESAEUQerBBfTEzswHOm/yZmVnLcFAyM7OW4aBkZmYtIysoSRop6U0LYJWMzGxjhKQjJJ0rab2ibGdJm3Svy2Zm1l/ljpSeBIZXKV+nOFeTpO1JC2/3Ab4ElLa7GA+cnNkHMzPr53KDkqi+wd9QYElG/TOAn0TEaGBpWflNpEW4eZ2Qxki6XtJcSSFpQhfXD5Y0RdIDkpZLmtbJdR+SNFPSEklPSPKCXjOzPlBzSriknxYfA5gkqXxd0iBgB+D+jPtsTxohVXqG7u08OxR4iLTTbVe73ULq4xLgZ8DHgWGVFxSPD6eSFgbvC+wCnCNpYUT8uht9MzOzHupqndI2xV8BW9Axeeoy4E+kUVBXFgNvq1K+ObCgSnlVETGVFECQNCXj+pcp0hhJ2pYqQak4Py8iSimRHin2fjoCcFAyM2uimkGplIlB0kXAYRHxYp33uQ74nqTSLrMhaRQpd15f//DvCNxcUXYTsL+k1SJieR/0ycxsQMp6pxQRB/QgIEEadaxD2lJ9TeAu4K+kDOPH9aDd3rABafPCcs+SAvZ61SpIOkhSu6T2hQt7sku8mZmVy915djBwGDAOWJ+KYBYR29aqXwS0XSTtDryvqP+niLi1nk73tYiYDEwGaGtrqzYBxMzM6pCb++4c4FPAlaStK+r6IY6I24Hb66nbQPN582SLEcCrwHPN746Z2cCVG5T2BD5T78imeCf1UEScWVF+OLBlRHy5nnZ7yXRSwC03Hmj3+yQzs+bKXaf0CvB0D+7zMaqPkG4nTdXOImmopO0kbUfq+8ji+8ji/CRJt1XU2bK4fj2gvH7JecDGkn4saQtJXwYmkDer0MzMelFuUDoNOLxaqqFMw4CXqpS/TJoAkasNmFUcQ4CJxeeTivMb8sZ+TyVTi2v2Jq2XKtUHICKeJAXGMaQ1V98FDvUaJTOz5st9fDce2BX4qKSHSbvOvi4iPtlF/cdJP/w/qSj/BGkWXpaImEZaM9XZ+QlVykZltHsnaQKGmZn1odyg9BxwTQ/ucyZwnqT1eeMx3jjgW8AhPWjXzMz6kdzt0A/oyU0i4uJiWvlxwDFF8Vzg8Ii4qCdtm5lZ/9Gt7dAltZHe2dwQES9LeguwNCJe7apuRJwPnC9pePHdq07NzKyD3MWzI0ipgnYgrVF6D/AEcBYp4elhuTd0MDIzs87kzr77ESn1zrqk6eElVwIf6aqypHWKzf0el7RI0ovlR/e7bWZm/VHu47txwLiIeKFiVvjfgJydZ38OjCal5plHnRkhzMysf8sNSkPouG1FyXDyNvkbB4yPiBm5HTMzs4En9/Hd70lZDkpC0iDg/wG3Va3R0QKqL541MzN7Xe5I6SjgTknvB9YgrTvaCngreduZfxc4SdL+EeHgZGZmVeWuU3pY0jbA14ClwGDSJIezI+KZjCaOA0YBCyTN5s0ZIWpufWFmZgND9jqliJgPfK/O+1xVZz0zMxtAsoOSpNWBram+yd/UWnUjYmJdvTMzswEld/HseOAXpIBUKYBBvdkpMzMbmHJn350N3ABsAqxJmiJeOtbsqrKk1SVNLBbPLpG0ovyot/NmZta/5AalDYFTImJ2RCyJiKXlR0b97wP7k2btvQYcSQp0/wC+Xk/Hzcys/8kNSjcAO/XgPp8FDi6Ssq4ArouIQ0kTJ8b3oF0zM+tHcic6HAz8UtL2wEO8eUr3JV3UHwE8XHx+ibQTLcDvgFMz+2BmZv1cblDag5Qq6OOkhKzluesC6CoozQE2Kv7+tWhvJrAjsLgb/TUzs34s9/HdGcDPgLUiYmhErFV2rJ1R/xpSUIO0JfpESU8CU4D/6W6nzcysf8odKQ0DzouIl+u5SUQcU/b5KklPk9ITPR4RN9TTppmZ9T+5I6VfAx+u9yaSxkh6PQBGxIyIOAv4naQx9bZrZmb9S+5I6Qng5CKAPMCbJzqc1UX9O0jTyhdUlL+1OOfFt2Zmlh2UDgT+RZoWXjk1PEjbotciqm/sty5Q1yNBMzPrf3KzhG9ST+OSri81AVwqqXyh7SBSLr2762nbzMz6n+yErHX6R/FXwAt0nP69DLgLuKDBfTAzs5VEd7KE702a1l0tS/gnq9WJiAOKuk8BZ9Q7e8/MzAaG3CzhpwPfIk1KmEf190O1fL+ivQ2Afwcejgg/vjMzMyB/pLQf8PmIqHezvhtJKYV+Imko0A68BRgq6UsZaYrMzGwAyF2ntApwfw/u0wbcXnzeC3iR9BjwK8ARPWjXzMz6kdygNBnYtwf3GQosKj5/BLgmIpaTAtWmuY0Ui3CvlzRXUkiakFFnG0l3Slpc1DtBksrOTyjaqjwGd+tfaGZmPdbp4ztJPy37ugqwT7EDbbXFs4d2cZ85wM6SfkNKxvqZonwdUoLXXENJWcovoesksEhaG7gF+D3wfmBz4CLS2qgzyy59hYrgGBFLutEvMzPrBbXeKW1T8f3+4u/mFeU5kx7OIm2n/hIwmxQkAMYAD2bUTzeKmApMBZA0JaPKPqSdcfePiMXAQ5I2Bw6XdFZElPoeETE/tx9mZtYYnQaliNitt24SEedLmgm8A7glIl4rTv0NOL637lPFjsAfioBUchNpNuAo4MmibIik2aQFvfcDx0fErAb2y8zMqsh6pyRpA0lvr1L+dkkjctqIiPaIuCYiXioruzEi/pjf3W7bAHi2ouzZsnMAj5HSKP0n8HlgCfBHSe/prFFJB0lql9S+cOHCXu6ymdnAlTsl/FLgCt6cfWEPYG/S5IUOJB0OnBMRS4rPncpI6NowETEdmF76Lulu0mjpm0DVd2URMZk0+YO2trburtkyM7NO5AalNuCQKuV/AE7vpM43gYtJI49v1mg7J6FrveaTtmIvN6Ls3Js7E7FCUjvQ6UjJzMwaIzcorQqsUaV8cCflHZK41pvQtRdMB06VNLhsNt14UlaKp6pVKKaLbwv8uSk9NDOz1+WuU5oBfK1K+SHAfb3XndokDZW0naTtSH0fWXwfWZyfJOm2siqXkaZ7T5G0taS9gKOB12feSfqepD0kvato9+ekoHRes/5dZmaW5I6UvgvcLmlb3sjMsDswmk52pJV0Qm4nIuKkzEvbSPn3SiYWx8XABNJGgq+vN4qIfxZrq84mpTZ6gbQ+qfxx4TDS+6ENgH8Cs4AxEXFvbv/NzKx36I2lOl1cKL0XOJIUiCD9eJ8eEVUfc0mqXH/0TtKaoXnF941Io5inImLbbva7ZbS1tUV7e3tfd8PMbKUhaWZEtFU7l711RRF8slMNRcTri28lHUBK6rp/RMwpykaSsiv8MrdNMzPr32qlGVonIp4vfa7VSOm6Gk4A9iwFpKLOHEnfAa4DLszvspmZ9Ve1RkoLJW0YEQuA56ieTkhF+aAu7jMCGFKlfDCwXk5Hzcys/6sVlHYHSiOgnqYcugW4QNJXSLP1AtgBOL84Z2ZmVjP33Z0AklYFtgKujYh5nV3fhS+TZsjdDawoylYh5aH7Sp1tmplZP9PlRIeIeLXYDv3Gem8SEQuBjxf55LYoih+NiMfrbdPMzPqf3Nl39wDbk7adqFtE/AX4S0/aMDOz/is3KF0AnFFM455J2iTvdRHxp97umJmZDTy5Qemy4m+1xKk5s+/MzMy6lBuU+iqhqpmZDSBZQSkievQuyczMLEd2mqFi59kxwPpUZBfP3aRP0kad1Pc7KTMzywtKkvYhpQJ6FVhIx+wOXW7SJ2k0affazUlZIMr5nZSZmQH5I6WTSFs+HB8RK7q6uIrJwNOkhbLzqJ6yyMzMBrjcoDQC+J86AxLAlsBoL5Y1M7NacneenQp8oAf3eZC0iZ6ZmVmnckdKtwCnStqKFGCWl5+MiKsrK1Rsd3EscJqk4zqp39XWF2ZmNgDkBqXzi7/HVjnX2USFyu0uBNxcpcwTHczMDMhfp5T7mK9cT7e7MDOzASZ7nVJ3lba+gNe3Pn86IjrMupMk4B2N6oOZma1cao6AJN0taVjZ90nl74okrSdpTtXKHT0JDK9Svk5xzszMrMvZdx8EVi/7fggwrOz7IGDjjPuU3h1VGgosyahvZmYDQHcf31VmY6h9sfTT4mMAkyS9UnZ6EGlL9Pu72QczM+unGvZOqbBN8VekHWeXlZ1bBvwJOKPBfTAzs5VEV0EpePNjt+wUQRGxG4Cki4DDIuLF7nXPzMwGkq6CkoBLJS0tvg8GLih7DLdGzk0i4oA6+2dmZgNIV0Hp4orvl1a55pJqFSVdD+wbES8WnzsVEZ/soh9mZjYA1AxKPRzh/IM3HvX9owftmJnZANHIxbMHVPtsZmYrt73Pnw7AFV/dsdfbrid9ULdJ2klSjwOgpDGSrpc0V1JImpBRZxtJd0paXNQ7ocgkUX7NpyU9LGlp8fdTPe2rmZl1X1OCEnA78IKkmyUd24MgNRR4CDgMWNzVxZLWJmU4fxZ4f1HvSODwsmt2BK4AfglsV/y9UlJPtuowM7M6NCsovQ34FDAD+Bgdg9QxuY1ExNSIODYirgJey6iyD7AmsH9EPFTUOxU4vGy09C3gjog4OSIeiYiTgWlFuZmZNVFTglJELI6IWyPi+IjYFdgKuBIYC/yggbfeEfhDRJSPqm4CNgJGlV1zc0W9m4CdOmtU0kGS2iW1L1y4sBe7a2Y2sDXrndL6kj4r6VxJjwAPAJsAJwO7N/DWG5Ae3ZV7tuxcrWs63Sk3IiZHRFtEtA0fXi3PrJlZ/3TtrLnMmrOIGU8+z84/vJ1rZ83t1fYbnWaoZD6wkLRZ4FeBGRGxtHYVMzNrJdfOmssxVz/IshXp7cncRYs55uoHAdhzdE5u7q41653SZcBS0kSDo4BvSNq+chZcA8wHRlSUjSg7V+ua+ZiZ2etOv+kxFi9f0aFs8fIVnH7TY712j2a9U9o3IkYC7wOuIc1yuxp4XtJ1Dbz1dGBXSYPLysYD84Cnyq4ZX1FvPHB3A/tlZrbSmbeo+qTnzsrr0ayRUsmTpCndDwOPAW8BPppbWdJQSdtJ2o7U95HF95HF+UmSbiurchnwCjBF0taS9gKOBs4q2wX3J8Duko6WtHkxG3A34Mc9+YeamfU3Gw0b0q3yejRrosNRkqYCi4DfA/8BzCz+rlOjaqU2YFZxDAEmFp9PKs5vCGxaujgi/kka9WwEtANnA2cCZ5VdczfwOWACaQLGfsDeETGje/9KM7P+7cg9NmPIaoM6lA1ZbRBH7rFZr91DbwwYGkfSdNLan2nAXRHxcsNv2iRtbW3R3t7e190wM2uKa2fN5airHmDZitfYeNgQjtxjs25PcpA0MyLaqp1ryuy7iOj9BElmZtZ0e47emMvvnQOsxLnvzMzMcjgomZlZy2jW4lkzM+snGvHYrsQjJTMzaxkOSmZm1jIclMzMrGU4KJmZWctwUDIzs5bhoGRmZi3DQcnMzFqGg5KZmbUMByUzM2sZTckS3p9JWgjMrrP6esBzvdgdM7Nm6cnv1zsjYni1Ew5KfUhSe2fp283MWlmjfr/8+M7MzFqGg5KZmbUMB6W+NbmvO2BmVqeG/H75nZKZmbUMj5TMzKxlOCiZmVnLcFBqAkmjJIWk7OmTkiZIeqmR/TIzazUOSj0kaYqkG6qUtxWBaBTwNLAhcH+Tu2dmlqWz37JmW7WvOzAQRMQKYH5f98PMrNV5pNQE1R7fSfqEpMckLZH0e0mfKxtZldcdJ+khSS9LukPSJk3/B5jZgCZpjKQZxe/Vs5J+JGn14txHJf1L0qrF93cXv2XnldX/gaRbc+7loNQHJI0ErgZuBN4L/BQ4rcqlawDHAAcCOwLDgPOqXGdm1hCSNgZ+C8wCRgNfAj4PTCouuQsYDJT+o3ssKSfe2LJmxgLTcu7noNQ7PirppfID+H2N678GPBERh0fEYxFxFdWDzarAIRFxb0Q8AJwBjJWk3v8nmJlV9XVgHvD1iHgkIm4Ajga+IWnNiHgJmAnsVlw/FvgZ8E5JG0paE3g/DkpN9Xtgu4rjCzWu3xy4r6JsRpXrlkbEY2Xf5wGrA2+rr5tmZt22BXBPRLxWVnYX6bfo3cX3abwxMvoQaWQ1oyjbCXgVuDfnZp7o0DteiYi/lhdIGtYL7b5a8b2UfsP/MWFmraD0mzSNNHLaAlibNHKaRho9LQCmR8SynAb949Y3HuWN568lO/RFR8zMuvAI8EFJ5fFiF2AZ8Lfi+12kd+BHAXcVM46nkYLSWDIf3YGDUl85D9hU0hmSNpO0F/DV4pyTEZpZX1lb0nblBzAV2Ag4R9IWkj4B/BD4WUS8AlD2Xmlf4I6irXuAtwMfxEGptUXEbODTwCeBPwPfBiYWp5f0Vb/MbMDblTTLrvz4NvAx0sy7+4ELgcuBYyvqTiO9EpoGEBFLSO+VlpL5PgmcJbxlSDoMOAkYFv5/ipkNUJ7o0EckHUKagbeQNLw9HpjigGRmA5mDUt95N2n4uy7wd9J7ppP6tEdmZn3Mj+/MzKxleKKDmZm1DAclMzNrGQ5KZmbWMhyUzAYgSSdKeqin15j1NgclsyYqdveMKsc9Zdc8VZR9sUr9GcW5I8rKNpF0qaS/S1oqaZ6kGyWN7mF3zyAl1zRrGk8JN2u+W4HKgFOZrPJp0j5avygVSNoa2Br4R1nZasAtpBxknwXmklLCfARYpyedLFLHvNSTNsy6yyMls+ZbGhHzK47nK665DNhR0rvKyr4EXEXHQLEVsClp3627I2J2REyPiIkRcVtXHZH0ZUlzJC2WdK2k9crOdXh8V4zybpB0mKS5kl6QdFGxX45Zr3BQMmtNzwG/AQ4AKLae3hf4ecV1C4HXgE+XtqPuhlFFm/8JfBh4DymvWS27kkZrHwb2Bj4FHNbN+5p1ykHJrPnetFOxpFOrXHchsH+xZcAngUUR0WFH44iYCxwKnAAsknSnpO9L2iqjH0OA/SJiVkT8kZSp/j8kvadGnReBg4sdSG8GrgTGZdzLLIuDklnzVdup+PQq190ECBhPenRXdRQTEWcDG5B2O76LNPK5v9pEiQpzI2JO2fcZpFHXFjXqPFzslVMyD1i/i/uYZfNEB7Pme9NOxdVExGuSLiblSPwgKTB1du2/gOuB6yUdRwpo36dsokQvWV55a/wft9aL/D8ms9Z2Iek9zi0RMS+nQpFp/lFgaBeXbizpHWXfdyD9JjxST0fNeoNHSmbNt4akDSrKVkTEwsoLI+KJYkbc4moNFTuDTiSNiB4mTS3/EGk6+eVd9GMxcLGkw0nvl84DboyIv3Tj32LWqxyUzJrvw8AzFWVzSVtHv0mV6eLl/g48QZroMIo00plDWvj6wy768RTwK9Isv/WAm4Evd1HHrKG8dYWZmbUMv1MyM7OW4aBkZmYtw0HJzMxahoOSmZm1DAclMzNrGQ5KZmbWMhyUzMysZTgomZlZy/j/GUoLEjbf2I4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rd = summary.iloc[0,:].sum() / summary.sum().sum() #prob(close to TSS) for a random variant\n",
    "enr = summary.iloc[0,:] / summary.sum(axis=0) #prob(close to TSS | high or low EMS)\n",
    "err = np.sqrt(enr*(1-enr) / summary.sum(axis=0)) #sem\n",
    "enr = enr/rd #divide by the random to get an enrichment\n",
    "err = err/rd #divide by the random, for an estimate of the error of the enrichment\n",
    "\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.rcParams['font.sans-serif'] = \"Arial\"\n",
    "plt.rcParams['font.family'] = \"sans-serif\"\n",
    "plt.errorbar([0,1], enr, err, fmt=\"o\")\n",
    "plt.xlabel(\"EMS bin\")\n",
    "plt.ylabel(\"Enrichment of variants \\n with distance to TSS <10kb\")\n",
    "plt.xticks([0,1],[\"High\",\"Low\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# End notes\n",
    "\n",
    "This tutorial explained basic properties of EMS, how to annotate EMS to a set of variant-gene pairs, and an example of enrichment analysis using EMS, in a cloud environment. \n",
    "\n",
    "For users who prefer file formats other than hail table, we have also prepared `tsv` files for a subset of variant-gene pairs (those with `ems_normalized>100`), and the corresponding tutorial is available at `tutorial_local.ipynb`.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Hail",
   "language": "python",
   "name": "hail"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}