tutorial_local.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",
    "assuming that readers do not have a cloud computing environment set up.\n",
    "\n",
    "As a first step, let's download the EMS file in tsv format to local (see the [finucane lab page](http://finucanelab.org/data) for instructions).\n",
    "\n",
    "Once we have the tsv (block-zipped) in local, we can read that in pandas:\n",
    "(Note: Here we are directly reading the `tsv.gz` in the cloud using hail and hadoop, but this can be done simpler by `df = pd.read_csv(your_local_data_dir)`), as shown in the comment below."
   ]
  },
  {
   "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:4042\n",
      "Welcome to\n",
      "     __  __     <>__\n",
      "    / /_/ /__  __/ /\n",
      "   / __  / _ `/ / /\n",
      "  /_/ /_/\\_,_/_/_/   version 0.2.26-2dcc3d963867\n",
      "LOGGING: writing to /home/hail/hail-20200914-2207-0.2.26-2dcc3d963867.log\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>v</th>\n",
       "      <th>g</th>\n",
       "      <th>rf_score_raw</th>\n",
       "      <th>ems</th>\n",
       "      <th>ems_normalized</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr10_100006504_T_C_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.98046</td>\n",
       "      <td>0.009245</td>\n",
       "      <td>761.83</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr10_100006605_T_C_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.87271</td>\n",
       "      <td>0.009629</td>\n",
       "      <td>793.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr10_100006744_C_G_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.94024</td>\n",
       "      <td>0.001634</td>\n",
       "      <td>134.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr10_100009635_T_G_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.91919</td>\n",
       "      <td>0.001803</td>\n",
       "      <td>148.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr10_100167711_G_A_b38</td>\n",
       "      <td>ENSG00000107566.13</td>\n",
       "      <td>0.92195</td>\n",
       "      <td>0.002227</td>\n",
       "      <td>183.52</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         v                   g  rf_score_raw       ems  \\\n",
       "0  chr10_100006504_T_C_b38  ENSG00000107554.16       0.98046  0.009245   \n",
       "1  chr10_100006605_T_C_b38  ENSG00000107554.16       0.87271  0.009629   \n",
       "2  chr10_100006744_C_G_b38  ENSG00000107554.16       0.94024  0.001634   \n",
       "3  chr10_100009635_T_G_b38  ENSG00000107554.16       0.91919  0.001803   \n",
       "4  chr10_100167711_G_A_b38  ENSG00000107566.13       0.92195  0.002227   \n",
       "\n",
       "   ems_normalized  \n",
       "0          761.83  \n",
       "1          793.46  \n",
       "2          134.62  \n",
       "3          148.62  \n",
       "4          183.52  "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#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\n",
    "\n",
    "tissue_name = \"Whole_Blood\"\n",
    "\n",
    "###\n",
    "#import and initialize hail \n",
    "import hail as hl\n",
    "import hail.expr.aggregators as agg\n",
    "hl.init(default_reference='GRCh38')\n",
    "#and read the file with hail hadoop\n",
    "fn = \"gs://expression-modifier-score/public/tsv/ems_top_{0}.tsv.bgz\".format(tissue_name)\n",
    "with hl.hadoop_open(fn, 'r') as f:\n",
    "    df = pd.read_csv(f, sep=\"\\t\")\n",
    "###\n",
    "\n",
    "#In local, this ### can be replaced with the following one line:\n",
    "# df = pd.read_csv(your_local_file_path, sep=\"\\t\")\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the EMS table as a pandas dataframe, filtering, aggregation and other major analyses are straightforward:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering\n",
    "Example of filtering to variant-genes with `ems_normalized`>500:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "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>v</th>\n",
       "      <th>g</th>\n",
       "      <th>rf_score_raw</th>\n",
       "      <th>ems</th>\n",
       "      <th>ems_normalized</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr10_100006504_T_C_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.98046</td>\n",
       "      <td>0.009245</td>\n",
       "      <td>761.83</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr10_100006605_T_C_b38</td>\n",
       "      <td>ENSG00000107554.16</td>\n",
       "      <td>0.87271</td>\n",
       "      <td>0.009629</td>\n",
       "      <td>793.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>chr10_100168530_TGA_T_b38</td>\n",
       "      <td>ENSG00000213341.10</td>\n",
       "      <td>0.87273</td>\n",
       "      <td>0.009629</td>\n",
       "      <td>793.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>chr10_100184325_G_C_b38</td>\n",
       "      <td>ENSG00000107566.13</td>\n",
       "      <td>0.89993</td>\n",
       "      <td>0.008127</td>\n",
       "      <td>669.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>chr10_100184325_G_C_b38</td>\n",
       "      <td>ENSG00000236308.1</td>\n",
       "      <td>0.89710</td>\n",
       "      <td>0.008082</td>\n",
       "      <td>665.98</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            v                   g  rf_score_raw       ems  \\\n",
       "0     chr10_100006504_T_C_b38  ENSG00000107554.16       0.98046  0.009245   \n",
       "1     chr10_100006605_T_C_b38  ENSG00000107554.16       0.87271  0.009629   \n",
       "9   chr10_100168530_TGA_T_b38  ENSG00000213341.10       0.87273  0.009629   \n",
       "11    chr10_100184325_G_C_b38  ENSG00000107566.13       0.89993  0.008127   \n",
       "12    chr10_100184325_G_C_b38   ENSG00000236308.1       0.89710  0.008082   \n",
       "\n",
       "    ems_normalized  \n",
       "0           761.83  \n",
       "1           793.46  \n",
       "9           793.46  \n",
       "11          669.68  \n",
       "12          665.98  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfhigh = df[df.ems_normalized>500]\n",
    "dfhigh.head()"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [
    {
     "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>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>chr10_100006504_T_C_b38</th>\n",
       "      <td>0.009245</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr10_100006605_T_C_b38</th>\n",
       "      <td>0.009629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr10_100006744_C_G_b38</th>\n",
       "      <td>0.001634</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr10_100009635_T_G_b38</th>\n",
       "      <td>0.001803</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr10_100167711_G_A_b38</th>\n",
       "      <td>0.002227</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              max\n",
       "v                                \n",
       "chr10_100006504_T_C_b38  0.009245\n",
       "chr10_100006605_T_C_b38  0.009629\n",
       "chr10_100006744_C_G_b38  0.001634\n",
       "chr10_100009635_T_G_b38  0.001803\n",
       "chr10_100167711_G_A_b38  0.002227"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfv = df.groupby(\"v\").ems.agg([\"max\"])\n",
    "dfv.head()"
   ]
  },
  {
   "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",
    "Let's assume we have a tsv file with 10000 rows of variant and genes:\n",
    "(We will make that by a random sampling of the `df`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "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>v</th>\n",
       "      <th>g</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr21_41362487_T_G_b38</td>\n",
       "      <td>ENSG00000183486.12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr17_73233212_C_T_b38</td>\n",
       "      <td>ENSG00000141219.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr2_74120467_G_C_b38</td>\n",
       "      <td>ENSG00000217702.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr22_20015967_C_CAG_b38</td>\n",
       "      <td>ENSG00000183597.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr17_78157921_T_A_b38</td>\n",
       "      <td>ENSG00000108639.7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          v                   g\n",
       "0    chr21_41362487_T_G_b38  ENSG00000183486.12\n",
       "1    chr17_73233212_C_T_b38  ENSG00000141219.15\n",
       "2     chr2_74120467_G_C_b38   ENSG00000217702.2\n",
       "3  chr22_20015967_C_CAG_b38  ENSG00000183597.15\n",
       "4    chr17_78157921_T_A_b38   ENSG00000108639.7"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vgs = df.sample(n=10000, random_state=1).reset_index()[[\"v\",\"g\"]]\n",
    "vgs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can annotate the EMS by join:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "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>v</th>\n",
       "      <th>g</th>\n",
       "      <th>ems</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>chr21_41362487_T_G_b38_ENSG00000183486.12</th>\n",
       "      <td>chr21_41362487_T_G_b38</td>\n",
       "      <td>ENSG00000183486.12</td>\n",
       "      <td>0.001962</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr17_73233212_C_T_b38_ENSG00000141219.15</th>\n",
       "      <td>chr17_73233212_C_T_b38</td>\n",
       "      <td>ENSG00000141219.15</td>\n",
       "      <td>0.012074</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr2_74120467_G_C_b38_ENSG00000217702.2</th>\n",
       "      <td>chr2_74120467_G_C_b38</td>\n",
       "      <td>ENSG00000217702.2</td>\n",
       "      <td>0.001641</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr22_20015967_C_CAG_b38_ENSG00000183597.15</th>\n",
       "      <td>chr22_20015967_C_CAG_b38</td>\n",
       "      <td>ENSG00000183597.15</td>\n",
       "      <td>0.004186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chr17_78157921_T_A_b38_ENSG00000108639.7</th>\n",
       "      <td>chr17_78157921_T_A_b38</td>\n",
       "      <td>ENSG00000108639.7</td>\n",
       "      <td>0.004465</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                                    v  \\\n",
       "chr21_41362487_T_G_b38_ENSG00000183486.12      chr21_41362487_T_G_b38   \n",
       "chr17_73233212_C_T_b38_ENSG00000141219.15      chr17_73233212_C_T_b38   \n",
       "chr2_74120467_G_C_b38_ENSG00000217702.2         chr2_74120467_G_C_b38   \n",
       "chr22_20015967_C_CAG_b38_ENSG00000183597.15  chr22_20015967_C_CAG_b38   \n",
       "chr17_78157921_T_A_b38_ENSG00000108639.7       chr17_78157921_T_A_b38   \n",
       "\n",
       "                                                              g       ems  \n",
       "chr21_41362487_T_G_b38_ENSG00000183486.12    ENSG00000183486.12  0.001962  \n",
       "chr17_73233212_C_T_b38_ENSG00000141219.15    ENSG00000141219.15  0.012074  \n",
       "chr2_74120467_G_C_b38_ENSG00000217702.2       ENSG00000217702.2  0.001641  \n",
       "chr22_20015967_C_CAG_b38_ENSG00000183597.15  ENSG00000183597.15  0.004186  \n",
       "chr17_78157921_T_A_b38_ENSG00000108639.7      ENSG00000108639.7  0.004465  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#index the both df with v and g\n",
    "df.index = df.v + \"_\" + df.g\n",
    "vgs.index = vgs.v + \"_\" + vgs.g\n",
    "#join\n",
    "vgs = vgs.join(df[\"ems\"], how=\"left\")\n",
    "vgs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# End notes\n",
    "\n",
    "This tutorial explained basic properties of EMS and how to annotate EMS to a set of variant-gene pairs in local python environment. Such approach would be simple and fast when the size of the variant-gene pairs are small and/or people are interested in a small subset of variant-genes with high EMS.\n",
    "\n",
    "\n",
    "For users who are familiar with cloud computing using hail and seeking of large scale computation using EMS, we recommend using the cloud; the corresponding tutorial is available at `tutorial.ipynb`."
   ]
  }
 ],
 "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
}