explore_imputed_genotypes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explore your results from *Imputer* on Open Humans.  \n",
    "\n",
    "If you'd like to learn how to explore you're imputed genotypes, please follow along in the cells below.  \n",
    "Plese see [openimpute.com/about](http://openimpute.com/about) or the [GitHub page](https://github.com/arvkevi/openhumansimputer) for questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import io\n",
    "import bz2\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import requests\n",
    "from ohapi import api\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "user = api.exchange_oauth2_member(os.environ.get('OH_ACCESS_TOKEN'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imputation results were returned in [vcf](https://samtools.github.io/hts-specs/VCFv4.3.pdf) format, the next cell will download your imputed genotype file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found your imputed genotype data in Open Humans!\n"
     ]
    }
   ],
   "source": [
    "for record in user['data']:\n",
    "    if record['basename'] == 'member.imputed.vcf.bz2':\n",
    "        print('Found your imputed genotype data in Open Humans!')\n",
    "        datafile = requests.get(record['download_url']).content\n",
    "        datafile = bz2.decompress(datafile)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select a Chromosome to analyze  \n",
    "Not all chromosomes will fit into memory on the server, so in this example, choose one.  \n",
    "The next three cells will subset your `vcf` file to a single chromosome and will convert the data into a [pandas](https://pandas.pydata.org/) `DataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# {1-22, 'X'}\n",
    "chrom = '5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished loading chromosome 5, stopping early\n"
     ]
    }
   ],
   "source": [
    "vcf_array = []\n",
    "for line in datafile.splitlines():\n",
    "\n",
    "    # this will stop searching the file if the chromsome has been read.\n",
    "    if len(vcf_array) > 0 and not line.startswith(chrom.encode('utf-8')):\n",
    "        print('Finished loading chromosome {}, stopping early'.format(chrom))\n",
    "        break\n",
    "    \n",
    "    # add the genetic variants from the user-defined chrom\n",
    "    if line.startswith(chrom.encode('utf-8')):\n",
    "        vcf_array.append(line.decode('utf-8'))\n",
    "del datafile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(io.StringIO('\\n'.join(vcf_array)),\n",
    "            sep='\\t',\n",
    "            comment='#',\n",
    "            header=None,\n",
    "            names=[\n",
    "                'CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE'\n",
    "            ])\n",
    "df.set_index(['CHROM', 'POS', 'REF', 'ALT'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 4 columns which serve as a unique index for genomic position:  \n",
    "`CHROM`, `POS`, `REF`, and `ALT`  \n",
    "\n",
    "To clean up the `DataFrame`, drop some columns that aren't used:  \n",
    "`ID` is a unique string, essentially a concatenated form of the unique index.  \n",
    "`QUAL` is `.` for every position, we did not store this data type.  \n",
    "`FILTER` is also `.` for every position, we did not store this data type.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>QUAL</th>\n",
       "      <th>FILTER</th>\n",
       "      <th>INFO</th>\n",
       "      <th>FORMAT</th>\n",
       "      <th>SAMPLE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">5</th>\n",
       "      <th>10056</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>5_10056_A_C</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>PR;INFO=0.001</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:0.998,0.002,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10062</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>5_10062_A_C</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>PR;INFO=0.0</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10063</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>5_10063_A_C</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>PR;INFO=0.052</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:0.885,0.115,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10067</th>\n",
       "      <th>T</th>\n",
       "      <th>A</th>\n",
       "      <td>5_10067_T_A</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>PR;INFO=0.0</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10099</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>5_10099_A_C</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>PR;INFO=0.97</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/1:0.015,0.985,0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              ID QUAL FILTER           INFO FORMAT  \\\n",
       "CHROM POS   REF ALT                                                  \n",
       "5     10056 A   C    5_10056_A_C    .      .  PR;INFO=0.001  GT:GP   \n",
       "      10062 A   C    5_10062_A_C    .      .    PR;INFO=0.0  GT:GP   \n",
       "      10063 A   C    5_10063_A_C    .      .  PR;INFO=0.052  GT:GP   \n",
       "      10067 T   A    5_10067_T_A    .      .    PR;INFO=0.0  GT:GP   \n",
       "      10099 A   C    5_10099_A_C    .      .   PR;INFO=0.97  GT:GP   \n",
       "\n",
       "                                  SAMPLE  \n",
       "CHROM POS   REF ALT                       \n",
       "5     10056 A   C    0/0:0.998,0.002,0.0  \n",
       "      10062 A   C        0/0:1.0,0.0,0.0  \n",
       "      10063 A   C    0/0:0.885,0.115,0.0  \n",
       "      10067 T   A        0/0:1.0,0.0,0.0  \n",
       "      10099 A   C    0/1:0.015,0.985,0.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop(['ID', 'QUAL', 'FILTER'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>INFO</th>\n",
       "      <th>FORMAT</th>\n",
       "      <th>SAMPLE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">5</th>\n",
       "      <th>10056</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>PR;INFO=0.001</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:0.998,0.002,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10062</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>PR;INFO=0.0</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10063</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>PR;INFO=0.052</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:0.885,0.115,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10067</th>\n",
       "      <th>T</th>\n",
       "      <th>A</th>\n",
       "      <td>PR;INFO=0.0</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/0:1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10099</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>PR;INFO=0.97</td>\n",
       "      <td>GT:GP</td>\n",
       "      <td>0/1:0.015,0.985,0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              INFO FORMAT               SAMPLE\n",
       "CHROM POS   REF ALT                                           \n",
       "5     10056 A   C    PR;INFO=0.001  GT:GP  0/0:0.998,0.002,0.0\n",
       "      10062 A   C      PR;INFO=0.0  GT:GP      0/0:1.0,0.0,0.0\n",
       "      10063 A   C    PR;INFO=0.052  GT:GP  0/0:0.885,0.115,0.0\n",
       "      10067 T   A      PR;INFO=0.0  GT:GP      0/0:1.0,0.0,0.0\n",
       "      10099 A   C     PR;INFO=0.97  GT:GP  0/1:0.015,0.985,0.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell will transform the remaining columns into a more usable format.  \n",
    "\n",
    "Specifically:  \n",
    "`INFO` will be parsed and float value after `INFO=` will be extracted.  \n",
    "`FORMAT` will be dropped, it is used to inform the reader what is contained in the `SAMPLE` field.  \n",
    "`SAMPLE` will be split and parsed into two columns as indicated by `FORMAT` -- `GT` and `GP`.  \n",
    "`GT` is your genotype at that genomic position.  \n",
    "`GP` is the probability output by the *Imputer* pipeline that the `GT` is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = df['INFO'].str.split(';', expand=True)\n",
    "df['INFO'] = info[1].str.split('=', expand=True)[1].astype(float)\n",
    "del info\n",
    "\n",
    "df.drop(['FORMAT'], axis=1, inplace=True)\n",
    "\n",
    "sample = df['SAMPLE'].str.split(':', expand=True)\n",
    "df['GT'] = sample[0]\n",
    "df['GP'] = sample[1]\n",
    "df.drop(['SAMPLE'], inplace=True, axis=1)\n",
    "del sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>INFO</th>\n",
       "      <th>GT</th>\n",
       "      <th>GP</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">5</th>\n",
       "      <th>10056</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0/0</td>\n",
       "      <td>0.998,0.002,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10062</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0/0</td>\n",
       "      <td>1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10063</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.052</td>\n",
       "      <td>0/0</td>\n",
       "      <td>0.885,0.115,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10067</th>\n",
       "      <th>T</th>\n",
       "      <th>A</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0/0</td>\n",
       "      <td>1.0,0.0,0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10099</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.970</td>\n",
       "      <td>0/1</td>\n",
       "      <td>0.015,0.985,0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      INFO   GT               GP\n",
       "CHROM POS   REF ALT                             \n",
       "5     10056 A   C    0.001  0/0  0.998,0.002,0.0\n",
       "      10062 A   C    0.000  0/0      1.0,0.0,0.0\n",
       "      10063 A   C    0.052  0/0  0.885,0.115,0.0\n",
       "      10067 T   A    0.000  0/0      1.0,0.0,0.0\n",
       "      10099 A   C    0.970  0/1  0.015,0.985,0.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## post-impuation QC  \n",
    "The INFO metric serves as a \"confidence\" that the genotype is called accurately. From the [IMPUTE2 homepage](http://mathgen.stats.ox.ac.uk/impute/impute_v2.html):  \n",
    "\n",
    ">Our metric typically takes values between 0 and 1, where values near 1 indicate that a SNP has been imputed with high certainty. The metric can occasionally take negative values when the imputation is very uncertain, and we automatically assign a value of -1 when the metric is undefined (e.g., because it wasn't calculated).\n",
    "\n",
    ">Investigators often use the info metric to remove poorly imputed SNPs from their association testing results. There is no universal cutoff value for post-imputation SNP filtering; various groups have used cutoffs of 0.3 and 0.5, for example, but the right threshold for your analysis may differ. One way to assess different info thresholds is to see whether they produce sensible Q-Q plots, although we emphasize that Q-Q plots can look bad for many reasons besides your post-imputation filtering scheme.\n",
    "\n",
    "Let's take an inclusive approach and set a cutoff of **0.3**  \n",
    "How many sites were considered to be imputed with high confidence?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All snps: 867428\n",
      "High confidence snps: 163374\n",
      "% hi conf snps: 18.83%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF6NJREFUeJzt3X+QXeV93/H318jYCjGWALPDSDQiE8U1gTGGHVDGM+nGpELgDuIP04EhkWDUqkOxx6mZNqL9Qy3UHdIOdQPjKFGLipQhxiqtK40trGpk7rjtIIyICTIQRmusoK1UFCOhsGZsV+63f9xn6WV9d++j1eoere77NXPnnvM9zznP81wJffb82EtkJpIk1Xhf0wOQJM0dhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGrzmh7AbLvoootyyZIlM9r3Rz/6Eeedd97sDugM55wHg3M++53qfJ9//vkfZuZHerU760JjyZIl7N27d0b7tlotRkZGZndAZzjnPBic89nvVOcbEX9Z087LU5KkaoaGJKmaoSFJqtYzNCLioxHxQsfrryPidyPigojYFRH7y/vC0j4i4uGIGI2IFyPi6o5jrS7t90fE6o76NRGxr+zzcEREqXftQ5LUjJ6hkZmvZuZVmXkVcA3wDvA1YB2wOzOXArvLOsCNwNLyWgtsgHYAAOuB64BrgfUdIbChtJ3Yb0WpT9WHJKkBJ3t56nrg+5n5l8BKYHOpbwZuKcsrgS3ZtgdYEBGXADcAuzLzaGYeA3YBK8q28zPzmWz/H6G2TDpWtz4kSQ042UdubwO+UpaHMvMwQGYejoiLS30RcLBjn7FSm64+1qU+XR/vERFraZ+pMDQ0RKvVOslptY2Pj89437nKOQ8G53z269d8q0MjIs4Fbgbu69W0Sy1nUK+WmRuBjQDDw8M502eVB+25bnDOg8I5n/36Nd+TuTx1I/BnmflGWX+jXFqivB8p9THg0o79FgOHetQXd6lP14ckqQEnc3nqdv7/pSmA7cBq4MHyvq2j/tmIeIL2Te/j5dLSTuBfddz8Xg7cl5lHI+LtiFgGPAusAh7p0cdpse9/HefOdd84nV10deDBT/e9T0maiarQiIhfAP428A86yg8CWyNiDfA6cGup7wBuAkZpP2l1F0AJhweA50q7+zPzaFm+G3gMmA88VV7T9SFJakBVaGTmO8CFk2pv0n6aanLbBO6Z4jibgE1d6nuBK7rUu/YhSWqGvxEuSapmaEiSqhkakqRqhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSapmaEiSqhkakqRqhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqlYVGhGxICKejIi/iIhXIuLXI+KCiNgVEfvL+8LSNiLi4YgYjYgXI+LqjuOsLu33R8Tqjvo1EbGv7PNwRESpd+1DktSM2jONPwC+mZl/E/g48AqwDtidmUuB3WUd4EZgaXmtBTZAOwCA9cB1wLXA+o4Q2FDaTuy3otSn6kOS1ICeoRER5wO/ATwKkJk/zcy3gJXA5tJsM3BLWV4JbMm2PcCCiLgEuAHYlZlHM/MYsAtYUbadn5nPZGYCWyYdq1sfkqQGzKto88vAXwH/MSI+DjwPfB4YyszDAJl5OCIuLu0XAQc79h8rtenqY13qTNPHe0TEWtpnKgwNDdFqtSqm9fOG5sO9V56Y0b6nYqbjnQ3j4+ON9t8E5zwYBm3O/ZpvTWjMA64GPpeZz0bEHzD9ZaLoUssZ1Ktl5kZgI8Dw8HCOjIyczO7veuTxbTy0r+YjmV0H7hjpe58TWq0WM/285irnPBgGbc79mm/NPY0xYCwzny3rT9IOkTfKpSXK+5GO9pd27L8YONSjvrhLnWn6kCQ1oGdoZOb/Bg5GxEdL6XrgZWA7MPEE1GpgW1neDqwqT1EtA46XS0w7geURsbDcAF8O7Czb3o6IZeWpqVWTjtWtD0lSA2qvxXwOeDwizgVeA+6iHThbI2IN8Dpwa2m7A7gJGAXeKW3JzKMR8QDwXGl3f2YeLct3A48B84GnygvgwSn6kCQ1oCo0MvMFYLjLpuu7tE3gnimOswnY1KW+F7iiS/3Nbn1Ikprhb4RLkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSapmaEiSqhkakqRqhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSapmaEiSqlWFRkQciIh9EfFCROwttQsiYldE7C/vC0s9IuLhiBiNiBcj4uqO46wu7fdHxOqO+jXl+KNl35iuD0lSM07mTOM3M/OqzBwu6+uA3Zm5FNhd1gFuBJaW11pgA7QDAFgPXAdcC6zvCIENpe3Efit69CFJasCpXJ5aCWwuy5uBWzrqW7JtD7AgIi4BbgB2ZebRzDwG7AJWlG3nZ+YzmZnAlknH6taHJKkB8yrbJfDfIiKBP87MjcBQZh4GyMzDEXFxabsIONix71ipTVcf61Jnmj7eIyLW0j5TYWhoiFarVTmt9xqaD/deeWJG+56KmY53NoyPjzfafxOc82AYtDn3a761ofHJzDxU/tHeFRF/MU3b6FLLGdSrlRDbCDA8PJwjIyMns/u7Hnl8Gw/tq/1IZs+BO0b63ueEVqvFTD+vuco5D4ZBm3O/5lt1eSozD5X3I8DXaN+TeKNcWqK8HynNx4BLO3ZfDBzqUV/cpc40fUiSGtAzNCLivIj40MQysBz4HrAdmHgCajWwrSxvB1aVp6iWAcfLJaadwPKIWFhugC8HdpZtb0fEsvLU1KpJx+rWhySpATXXYoaAr5WnYOcBf5qZ34yI54CtEbEGeB24tbTfAdwEjALvAHcBZObRiHgAeK60uz8zj5blu4HHgPnAU+UF8OAUfUiSGtAzNDLzNeDjXepvAtd3qSdwzxTH2gRs6lLfC1xR24ckqRn+RrgkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSapmaEiSqhkakqRqhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpWnVoRMQ5EfHdiPh6Wb8sIp6NiP0R8dWIOLfUP1DWR8v2JR3HuK/UX42IGzrqK0ptNCLWddS79iFJasbJnGl8HnilY/33gS9l5lLgGLCm1NcAxzLzV4AvlXZExOXAbcCvASuAPyxBdA7wZeBG4HLg9tJ2uj4kSQ2oCo2IWAx8GvgPZT2ATwFPliabgVvK8sqyTtl+fWm/EngiM3+SmT8ARoFry2s0M1/LzJ8CTwAre/QhSWrAvMp2/w74J8CHyvqFwFuZeaKsjwGLyvIi4CBAZp6IiOOl/SJgT8cxO/c5OKl+XY8+3iMi1gJrAYaGhmi1WpXTeq+h+XDvlSd6N5xlMx3vbBgfH2+0/yY458EwaHPu13x7hkZE/B3gSGY+HxEjE+UuTbPHtqnq3c52pmv/88XMjcBGgOHh4RwZGenWrKdHHt/GQ/tqc3T2HLhjpO99Tmi1Wsz085qrnPNgGLQ592u+Nf9CfhK4OSJuAj4InE/7zGNBRMwrZwKLgUOl/RhwKTAWEfOADwNHO+oTOvfpVv/hNH1IkhrQ855GZt6XmYszcwntG9nfysw7gKeBz5Rmq4FtZXl7Wads/1ZmZqnfVp6uugxYCnwHeA5YWp6UOrf0sb3sM1UfkqQGnMrvafwe8IWIGKV9/+HRUn8UuLDUvwCsA8jMl4CtwMvAN4F7MvNn5Szis8BO2k9nbS1tp+tDktSAk7qAn5ktoFWWX6P95NPkNj8Gbp1i/y8CX+xS3wHs6FLv2ockqRn+RrgkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSapmaEiSqhkakqRqhoYkqZqhIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGo9QyMiPhgR34mIP4+IlyLiX5T6ZRHxbETsj4ivRsS5pf6Bsj5ati/pONZ9pf5qRNzQUV9RaqMRsa6j3rUPSVIzas40fgJ8KjM/DlwFrIiIZcDvA1/KzKXAMWBNab8GOJaZvwJ8qbQjIi4HbgN+DVgB/GFEnBMR5wBfBm4ELgduL22Zpg9JUgN6hka2jZfV95dXAp8Cniz1zcAtZXllWadsvz4iotSfyMyfZOYPgFHg2vIazczXMvOnwBPAyrLPVH1IkhpQdU+jnBG8ABwBdgHfB97KzBOlyRiwqCwvAg4ClO3HgQs765P2map+4TR9SJIaMK+mUWb+DLgqIhYAXwM+1q1ZeY8ptk1V7xZc07X/ORGxFlgLMDQ0RKvV6tasp6H5cO+VJ3o3nGUzHe9sGB8fb7T/JjjnwTBoc+7XfKtCY0JmvhURLWAZsCAi5pUzgcXAodJsDLgUGIuIecCHgaMd9Qmd+3Sr/3CaPiaPayOwEWB4eDhHRkZOZlrveuTxbTy076Q+kllx4I6Rvvc5odVqMdPPa65yzoNh0Obcr/nWPD31kXKGQUTMB34LeAV4GvhMabYa2FaWt5d1yvZvZWaW+m3l6arLgKXAd4DngKXlSalzad8s3172maoPSVIDan6svgTYXJ5yeh+wNTO/HhEvA09ExL8Evgs8Wto/CvxJRIzSPsO4DSAzX4qIrcDLwAngnnLZi4j4LLATOAfYlJkvlWP93hR9SJIa0DM0MvNF4BNd6q/RfvJpcv3HwK1THOuLwBe71HcAO2r7kCQ1w98IlyRVMzQkSdUMDUlSNUNDklTN0JAkVTM0JEnVDA1JUjVDQ5JUzdCQJFUzNCRJ1QwNSVI1Q0OSVM3QkCRVMzQkSdUMDUlSNUNDklTN0JAkVTM0JEnVDA1JUjVDQ5JUzdCQJFXrGRoRcWlEPB0Rr0TESxHx+VK/ICJ2RcT+8r6w1CMiHo6I0Yh4MSKu7jjW6tJ+f0Ss7qhfExH7yj4PR0RM14ckqRk1ZxongHsz82PAMuCeiLgcWAfszsylwO6yDnAjsLS81gIboB0AwHrgOuBaYH1HCGwobSf2W1HqU/UhSWpAz9DIzMOZ+Wdl+W3gFWARsBLYXJptBm4pyyuBLdm2B1gQEZcANwC7MvNoZh4DdgEryrbzM/OZzExgy6RjdetDktSAk7qnERFLgE8AzwJDmXkY2sECXFyaLQIOduw2VmrT1ce61JmmD0lSA+bVNoyIXwT+M/C7mfnX5bZD16ZdajmDerWIWEv78hZDQ0O0Wq2T2f1dQ/Ph3itPzGjfUzHT8c6G8fHxRvtvgnMeDIM2537Ntyo0IuL9tAPj8cz8L6X8RkRckpmHyyWmI6U+Blzasfti4FCpj0yqt0p9cZf20/XxHpm5EdgIMDw8nCMjI92a9fTI49t4aF91js6aA3eM9L3PCa1Wi5l+XnOVcx4Mgzbnfs235umpAB4FXsnMf9uxaTsw8QTUamBbR31VeYpqGXC8XFraCSyPiIXlBvhyYGfZ9nZELCt9rZp0rG59SJIaUPNj9SeB3wH2RcQLpfZPgQeBrRGxBngduLVs2wHcBIwC7wB3AWTm0Yh4AHiutLs/M4+W5buBx4D5wFPlxTR9SJIa0DM0MvN/0P2+A8D1XdoncM8Ux9oEbOpS3wtc0aX+Zrc+JEnN8DfCJUnVDA1JUjVDQ5JUzdCQJFUzNCRJ1QwNSVI1Q0OSVM3QkCRVMzQkSdUMDUlSNUNDklSt/98DLklnsSXrvtFIv4+tOK8v/XimIUmqZmhIkqoZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSarWMzQiYlNEHImI73XULoiIXRGxv7wvLPWIiIcjYjQiXoyIqzv2WV3a74+I1R31ayJiX9nn4YiI6fqQJDWn5kzjMWDFpNo6YHdmLgV2l3WAG4Gl5bUW2ADtAADWA9cB1wLrO0JgQ2k7sd+KHn1IkhrSMzQy89vA0UnllcDmsrwZuKWjviXb9gALIuIS4AZgV2YezcxjwC5gRdl2fmY+k5kJbJl0rG59SJIaMtN7GkOZeRigvF9c6ouAgx3txkptuvpYl/p0fUiSGjLb33IbXWo5g/rJdRqxlvYlLoaGhmi1Wid7CACG5sO9V56Y0b6nYqbjnQ3j4+ON9t8E5zwYmppzE/+GQP/mO9PQeCMiLsnMw+US05FSHwMu7Wi3GDhU6iOT6q1SX9yl/XR9/JzM3AhsBBgeHs6RkZGpmk7rkce38dC+/n9b/IE7Rvre54RWq8VMP6+5yjkPhqbmfGeDX43ej/nO9PLUdmDiCajVwLaO+qryFNUy4Hi5tLQTWB4RC8sN8OXAzrLt7YhYVp6aWjXpWN36kCQ1pOeP1RHxFdpnCRdFxBjtp6AeBLZGxBrgdeDW0nwHcBMwCrwD3AWQmUcj4gHgudLu/sycuLl+N+0ntOYDT5UX0/QhSWpIz9DIzNun2HR9l7YJ3DPFcTYBm7rU9wJXdKm/2a0PSVJz/I1wSVI1Q0OSVM3QkCRVMzQkSdUMDUlSNUNDklTN0JAkVTM0JEnVDA1JUjVDQ5JUzdCQJFUzNCRJ1QwNSVI1Q0OSVM3QkCRVMzQkSdUMDUlSNUNDklTN0JAkVTM0JEnVDA1JUrUzPjQiYkVEvBoRoxGxrunxSNIgO6NDIyLOAb4M3AhcDtweEZc3OypJGlzzmh5AD9cCo5n5GkBEPAGsBF5udFSzbMm6bzTW971XnuDOBvo/8OCn+96npFN3pofGIuBgx/oYcF1DY9EsGsSgbJJz1myJzGx6DFOKiFuBGzLz75X13wGuzczPTWq3FlhbVj8KvDrDLi8CfjjDfecq5zwYnPPZ71Tn+0uZ+ZFejc70M40x4NKO9cXAocmNMnMjsPFUO4uIvZk5fKrHmUuc82Bwzme/fs33jL4RDjwHLI2IyyLiXOA2YHvDY5KkgXVGn2lk5omI+CywEzgH2JSZLzU8LEkaWGd0aABk5g5gR5+6O+VLXHOQcx4Mzvns15f5ntE3wiVJZ5Yz/Z6GJOkMMpCh0eurSSLiAxHx1bL92YhY0v9Rzq6KOX8hIl6OiBcjYndE/FIT45xNtV9BExGfiYiMiDn9pE3NfCPi75Y/55ci4k/7PcbZVvH3+m9ExNMR8d3yd/umJsY5myJiU0QciYjvTbE9IuLh8pm8GBFXz+oAMnOgXrRvqH8f+GXgXODPgcsntfmHwB+V5duArzY97j7M+TeBXyjLdw/CnEu7DwHfBvYAw02P+zT/GS8FvgssLOsXNz3uPsx5I3B3Wb4cOND0uGdh3r8BXA18b4rtNwFPAQEsA56dzf4H8Uzj3a8mycyfAhNfTdJpJbC5LD8JXB8R0ccxzraec87MpzPznbK6h/bvxMxlNX/OAA8A/xr4cT8HdxrUzPfvA1/OzGMAmXmkz2OcbTVzTuD8svxhuvye11yTmd8Gjk7TZCWwJdv2AAsi4pLZ6n8QQ6PbV5MsmqpNZp4AjgMX9mV0p0fNnDutof2TylzWc84R8Qng0sz8ej8HdprU/Bn/KvCrEfE/I2JPRKzo2+hOj5o5/3PgtyNijPZTmJ/j7Hey/72flDP+kdvToNsZw+RHyGrazCXV84mI3waGgb91Wkd0+k0754h4H/Al4M5+Deg0q/kznkf7EtUI7TPJ/x4RV2TmW6d5bKdLzZxvBx7LzIci4teBPylz/r+nf3iNOa3/fg3imUbNV5O82yYi5tE+rZ3udPBMV/V1LBHxW8A/A27OzJ/0aWynS685fwi4AmhFxAHa1363z+Gb4bV/r7dl5v/JzB/Q/o62pX0a3+lQM+c1wFaAzHwG+CDt72g6m1X99z5TgxgaNV9Nsh1YXZY/A3wryx2mOarnnMulmj+mHRhz/Vo39JhzZh7PzIsyc0lmLqF9H+fmzNzbzHBPWc3f6/9K+4EHIuIi2perXuvrKGdXzZxfB64HiIiP0Q6Nv+rrKPtvO7CqPEW1DDiemYdn6+ADd3kqp/hqkoi4H9ibmduBR2mfxo7SPsO4rbkRn7rKOf8b4BeB/1Tu+b+emTc3NuhTVDnns0blfHcCyyPiZeBnwD/OzDebG/WpqZzzvcC/j4h/RPsSzZ1z/AdAIuIrtC8xXlTu1awH3g+QmX9E+97NTcAo8A5w16z2P8c/P0lSHw3i5SlJ0gwZGpKkaoaGJKmaoSFJqmZoSJKqGRqSpGqGhiSpmqEhSar2/wCNfxRQ8niG6gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_snps = df.shape[0]\n",
    "n_hiconf_snps = df.loc[df['INFO']>=0.3].shape[0]\n",
    "print('All snps: {}'.format(n_snps))\n",
    "print('High confidence snps: {}'.format(n_hiconf_snps))\n",
    "print('% hi conf snps: {0:.2f}%'.format((n_hiconf_snps / n_snps)*100))\n",
    "\n",
    "df['INFO'].hist();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the possible values for the `GT` field:  \n",
    "`0/0` indicates you are homozygous reference at this position.  \n",
    "`0/1` indicates you are heterozygous at this position.  \n",
    "`1/1` indicates you are homozygous alternate at this position.  \n",
    "`./.` indicates that the *GP* threshold of `0.6` was not met at this position, so no call was made.  \n",
    "*The `./.` threshold was arbitrarily chosen.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0/0    567302\n",
       "0/1    155154\n",
       "1/1    138541\n",
       "./.      6431\n",
       "Name: GT, dtype: int64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['GT'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell will split the genotype probabilities into three columns:  \n",
    "1. `GP_homref` for homozygous reference calls\n",
    "2. `GP_het` for heterozygous calls\n",
    "3. `GP_homalt` for homozygous alterante calls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp = df['GP'].str.split(',', expand=True)\n",
    "df['GP_homref'] = gp[0].astype(float)\n",
    "df['GP_het'] = gp[1].astype(float)\n",
    "df['GP_homalt'] = gp[2].astype(float)\n",
    "del gp\n",
    "df.drop(['GP'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>INFO</th>\n",
       "      <th>GT</th>\n",
       "      <th>GP_homref</th>\n",
       "      <th>GP_het</th>\n",
       "      <th>GP_homalt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">5</th>\n",
       "      <th>10056</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0/0</td>\n",
       "      <td>0.998</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10062</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0/0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10063</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.052</td>\n",
       "      <td>0/0</td>\n",
       "      <td>0.885</td>\n",
       "      <td>0.115</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10067</th>\n",
       "      <th>T</th>\n",
       "      <th>A</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0/0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10099</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <td>0.970</td>\n",
       "      <td>0/1</td>\n",
       "      <td>0.015</td>\n",
       "      <td>0.985</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      INFO   GT  GP_homref  GP_het  GP_homalt\n",
       "CHROM POS   REF ALT                                          \n",
       "5     10056 A   C    0.001  0/0      0.998   0.002        0.0\n",
       "      10062 A   C    0.000  0/0      1.000   0.000        0.0\n",
       "      10063 A   C    0.052  0/0      0.885   0.115        0.0\n",
       "      10067 T   A    0.000  0/0      1.000   0.000        0.0\n",
       "      10099 A   C    0.970  0/1      0.015   0.985        0.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a closer look at the genotype probabilities, `GP`, for the uncalled sites, they are very uncertain, with a maximum genotype probability of much lower than `1.0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CHROM  POS     REF  ALT\n",
       "5      16120   A    T      0.555\n",
       "       83380   A    C      0.555\n",
       "       92468   G    GTA    0.555\n",
       "       262306  C    G      0.555\n",
       "       272367  C    T      0.555\n",
       "dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.loc[df['GT'] == './.'][['GP_homref', 'GP_het', 'GP_homalt']].astype(float).max(axis=1).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully this notebook will help you explore your imputed genotype data, if you have any questions or suggestions please open an issue on the [GitHub page](https://github.com/arvkevi/openhumansimputer)."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}