4 - Host genotypes - Data preparation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Warning: usually don't do this kind of inclusion : using *, you 'merge' the namespaces\n",
    "#          but in this case we actually want to mimic the imports as if they were imported here\n",
    "from src.setup import *\n",
    "\n",
    "from src.host import run_plink\n",
    "from src import host\n",
    "from src import common\n",
    "\n",
    "# This prevents to restart the kernel whenever an external file is edited\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# Comutation management\n",
    "DATA_THINNING = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "**14'738'460 variants, 622 people.**\n",
    "\n",
    "## Data processing plan\n",
    "\n",
    "### Preprocessing \n",
    "\n",
    "1. Data thinning (if `DATA_THINNING == True`)\n",
    "1. Create SNP map (SNP id and SNP position)\n",
    "\n",
    "### Filtering\n",
    "1. Remove missing values (variants)\n",
    "1. Remove missing values (individuals)\n",
    "1. Filter out rare variants and variants at Hardy-Weinberg disequilibrium\n",
    "\n",
    "## Computation and path management\n",
    "\n",
    "### Idea\n",
    "\n",
    "Sometimes we don't want to re-do all the computations if the resulting data is already calculated. The features provided below aim to easily manage the paths and the computations:\n",
    "1. A way of easily running shell commands\n",
    "1. A dedicated function to easily run plink commands on the shell\n",
    "1. A way of checking if the results already exist, and loading data if this is the case\n",
    "    * Unless we specifically specify to re-do the calculations.\n",
    "\n",
    "###### DOCUMENTATION \n",
    "\n",
    "Running shell commands and managing the outputs : https://janakiev.com/til/python-shell-commands/, https://www.pythonforbeginners.com/os/subprocess-for-system-administrators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General principle\n",
    "\n",
    "* `run_shell_command`: manage to run shell commands with a specific package (`subprocess`) that helps to manage the outputs. Returns a 2-tuple (stdout, stderr).\n",
    "* `run_plink`: manage plink version, binary files location, output location, display errors, check if results already exist.\n",
    "\n",
    "*NOTE: possible improvement: automatically find the file extension. Also, if the command results in several output files, we should check that all of them are present. So **be careful, we assume the user doesn't  only a fraction of the output files**.*\n",
    "\n",
    "###### Workflow\n",
    "1. Call `run_plink` command with the specific command (e.g. `--freq`), the input, the output and the file extension to look for (to avoid re-doing the computations). \n",
    "1. The function checks if the result already exist (needs the output and the extension resulting from the command)\n",
    "    * already exists: display a message, abort the process (returns 2-tuple of empty strings). This occurs if `force==False` (which is the default value).\n",
    "1. The command is runned with plink2 by default\n",
    "1. The error messages are displayed (by default, `verbose_err=True`)\n",
    "1. The outputs are returned as a tuple (stdout, stderr)\n",
    "\n",
    "###### Modify behaviour\n",
    "\n",
    "Change the boolean value defined in the first cell of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Path management"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_raw = PATH_HOST_RAW_PLINK_DATA\n",
    "path_out = PATH_PLINK_DATA\n",
    "file_prefix = 'host_geno'\n",
    "\n",
    "if not os.path.exists(PATH_PLINK_LOG):\n",
    "    os.system('mkdir '+PATH_PLINK_LOG)\n",
    "    print(\"Created the directory\", PATH_PLINK_LOG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview on raw data\n",
    "\n",
    "### Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 6.6G\r\n",
      "-rwxr-xr-x 1 mazeller SV-BA4-unit  245  9 aoû 16:11 dictionary_wes_vcf.txt\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit 2.2G 13 aoû 15:09 hbv_gilead_bcftools.bed\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit 431M 13 aoû 15:09 hbv_gilead_bcftools.bim\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit  18K 13 aoû 15:09 hbv_gilead_bcftools.fam\r\n",
      "-rwxr-xr-x 1 mazeller SV-BA4-unit 2.2G  9 aoû 16:11 hbv_gilead.bed\r\n",
      "-rwxr-xr-x 1 mazeller SV-BA4-unit 431M  9 aoû 16:11 hbv_gilead.bim\r\n",
      "-rwxr-xr-x 1 mazeller SV-BA4-unit  18K  9 aoû 16:11 hbv_gilead.fam\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit 1.4G 11 aoû 01:59 hbv_gilead.frq\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit  930 11 aoû 01:59 hbv_gilead.log\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit  13K 11 aoû 01:59 hbv_gilead.nosex\r\n",
      "-rw-r--r-- 1 mazeller SV-BA4-unit  677 11 aoû 17:56 hbv_gilead_validate.log\r\n"
     ]
    }
   ],
   "source": [
    "path = PATH_HOST_RAW_DATA\n",
    "!ls -lh $path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those are **plink binary files**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dictionnary of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info about WES_VCF by Ondrej.Podlaha@gilead.com on Dec 12 2018.\r\n",
      "\r\n",
      "1. Human reference genome: GRCh37.87\r\n",
      "2. Variant calling was performed using GATK (version nightly-2016-08-27-g667f78b; HaplotypeCaller)\r\n",
      "3. dbSNP v138 was used by HaplotypeCaller\r\n",
      " \r\n"
     ]
    }
   ],
   "source": [
    "!cat $path/dictionary_wes_vcf.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data integrity\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "run_plink: command '--validate', the file 'data/plink/host_geno_validate.log' already exists (force is set to False).\n",
      "PLINK v2.00a2LM 64-bit Intel (31 Jul 2019)     www.cog-genomics.org/plink/2.0/\r\n",
      "(C) 2005-2019 Shaun Purcell, Christopher Chang   GNU General Public License v3\r\n",
      "Logging to data/plink/host_geno_validate.log.\r\n",
      "Options in effect:\r\n",
      "  --bfile data/raw/wes_plink/hbv_gilead_bcftools\r\n",
      "  --out data/plink/host_geno_validate\r\n",
      "  --validate\r\n",
      "\r\n",
      "Start time: Tue Aug 13 17:31:45 2019\r\n",
      "64232 MiB RAM detected; reserving 32116 MiB for main workspace.\r\n",
      "Using up to 16 threads (change this with --threads).\r\n",
      "622 samples (0 females, 0 males, 622 ambiguous; 622 founders) loaded from\r\n",
      "data/raw/wes_plink/hbv_gilead_bcftools.fam.\r\n",
      "14738460 variants loaded from data/raw/wes_plink/hbv_gilead_bcftools.bim.\r\n",
      "Validating data/raw/wes_plink/hbv_gilead_bcftools.bed... done.\r\n",
      "End time: Tue Aug 13 17:31:51 2019\r\n"
     ]
    }
   ],
   "source": [
    "run_plink('--validate', file=path_raw, out=path_out+file_prefix+'_validate', \n",
    "               extension='log', log_name='validate')\n",
    "!cat $PATH_PLINK_LOG'validate.log'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## .bim file\n",
    "\n",
    "This is an extended .map file, \n",
    "1. Chromosode code\n",
    "1. Variant id (SNP)\n",
    "1. Position in centimorgan\n",
    "1. Baise-pair coordinate\n",
    "1. ALT allele code\n",
    "1. REF allele code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\trs779258992\t0\t10146\tA\tAC\r\n",
      "1\t1:10400:C:T\t0\t10400\tT\tC\r\n",
      "1\trs55998931\t0\t10492\tT\tC\r\n",
      "1\trs58108140\t0\t10583\tA\tG\r\n",
      "1\trs376342519\t0\t10616\tC\tCCGCCGTTGCAAAGGCGCGCCG\r\n",
      "1\t1:13012:A:G\t0\t13012\tA\tG\r\n",
      "1\t1:13054:C:T\t0\t13054\tT\tC\r\n",
      "1\t1:13250:G:T\t0\t13250\tT\tG\r\n",
      "1\trs531730856\t0\t13273\tC\tG\r\n",
      "1\trs538791886\t0\t13289\tC\tCCT\r\n"
     ]
    }
   ],
   "source": [
    "!head -10 $path_raw'.bim'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arbitrary thinning of the data\n",
    "\n",
    "We use the plink2 `--thin <p>` option to **reduce the amount of data** (keep entries with probability p). We'll make sure that the code runs first. When everything works, we can simply set `DATA_THINNING = False`.\n",
    "\n",
    "Management of the **data being taken into account**: we simply put `path_raw` to target the reduced data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# We change the default raw data for subsequent analyses\n",
    "if DATA_THINNING == True :\n",
    "    file_prefix = 'host_reduced'\n",
    "    run_plink(\"--thin 0.001 --make-bed\", file=path_raw, out=path_out+file_prefix, \n",
    "              extension='bed', log_name='thinning')\n",
    "    path_raw = path_out+file_prefix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note on plink\n",
    "\n",
    "| Feature | Summary statistic | Inclusion criteria  |\n",
    "|----------------------------|-----------|-----------|\n",
    "| Missingness per individual | --missing | --mind N  |\n",
    "| Missingness per marker     | --missing | --geno N  |\n",
    "| Allele frequency           | --freq    | --maf N   |\n",
    "| Hardy-Weinberg equilibrium | --hardy   | --hwe     |\n",
    "| Mendel error rates         | --mendel  |  --me N M |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SNP map generation\n",
    "\n",
    "We create a Python dictionnary object and store it to easily convert the SNP references to their genome location. Those information are extracted from the `.bim` file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "manage_pickle: the file 'data/SNP_dict' has been loaded\n"
     ]
    }
   ],
   "source": [
    "# Create dictionnary\n",
    "# Use custom manage_pickle function from common.py\n",
    "def create_dict() :\n",
    "    # Create DataFrame from .bim file\n",
    "    df_bim = pd.read_csv(path_raw+'.bim', sep='\\s+', header=None,\n",
    "                         usecols=[0, 1, 3])\n",
    "    df_bim.columns = ['#CHROM', 'SNP', 'COORD']\n",
    "    print(df_bim.head(1))\n",
    "    return dict(zip(df_bim['SNP'], df_bim['COORD']))\n",
    "\n",
    "SNP_dic = common.manage_pickle(path=PATH_SNP_MAP, fun=create_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Missing values \n",
    "\n",
    "`.smiss` for sample-based missing data report:\n",
    "* MISSING_CT: Number of missing hardcalls, not counting het haploids\n",
    "* OBS_CT = Denominator (#males on chrY, otherwise #samples)\n",
    "* F_MISS = Missing hardcall rate, not counting het haploids\n",
    "\n",
    "###### Pandas - memory efficiency\n",
    "\n",
    "When reading big files, one must provide to `pandas.read_csv` a map of the data types that are read to improve the memory usage:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtype = {'F_MISS':float, 'ID':str, '#CHROM':str}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Run plink missing report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink2 --bfile data/raw/wes_plink/hbv_gilead_bcftools --out data/plink/host_geno_clean_missing --missing '\n",
      "Log written to 'data/plink/log/missing.log'\n"
     ]
    }
   ],
   "source": [
    "path_missing = path_out+file_prefix+'_missing'\n",
    "o,e = run_plink(\"--missing\", file=path_raw, out=path_missing, \n",
    "                extension='vmiss', log_name='missing', force=True)\n",
    "df_miss_ind = pd.read_csv(path_missing+'.smiss', sep='\\s+', \n",
    "                          usecols=['F_MISS'], dtype=dtype)\n",
    "df_miss_var = pd.read_csv(path_missing+'.vmiss', sep='\\s+', \n",
    "                          usecols=['#CHROM', 'ID', 'F_MISS'], dtype=dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Plot function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_missing(df_ind, df_var):\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(15, 3))\n",
    "    sns.distplot(df_ind['F_MISS'], norm_hist=False, kde=False, ax=ax[0])\n",
    "    ax[0].set_ylabel('number of individuals'); ax[0].set_xlabel('rate of missing variants');\n",
    "    sns.distplot(df_var.F_MISS.dropna(), kde=False, norm_hist=False, ax=ax[1])\n",
    "    ax[1].set_xlabel('rate of missing individuals'); ax[1].set_ylabel('number of variants');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the rate of missing variants for each individual (this tells us about how to filter individuals out). We also plot the rate of missing individuals for each variants (this tells us about how to filter individuals out)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "75.2% variants with missing rate of 100%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAADbCAYAAADK8LgjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debgsZXnv/e8PkEFBBtkxyOBGRRNiUHELGjWAQ8QhYCIqoAYQw4mzcYh6zKuo8T2occAhGlRkuGSS6HEnIg7IFo/KsBFk8iCEyY1EEBUFIuN9/qha2HvtNdQaunv1Wt/PddXVNfdd1b3W00/VU/eTqkKSJEmSNJrWG3YAkiRJkqTZs1InSZIkSSPMSp0kSZIkjTArdZIkSZI0wqzUSZIkSdIIs1InSZIkSSPMSp0kSSMsydFJbkxySYd1P5Lkwnb4SZJfDyJGSVJ/xX7qJEkaXUn+HLgVOK6qHj2D7V4LPK6qXt634CRJA+GdOkmSRlhVnQX8sndekocnOT3J+Um+m+SPJtj0AODEgQQpSeqrDYYdgCRJmndHAX9XVVck2R34F+BpYwuTPBTYEfj2kOKTJM0jK3WSJC0iSTYF/gz4YpKx2RuNW21/4NSqumeQsUmS+sNKnSRJi8t6wK+r6rFTrLM/8OoBxSNJ6jOfqZMkaRGpqt8AVyd5IUAajxlb3j5ftyXwgyGFKEmaZ1bqJEkaYUlOpKmgPSrJmiSHAi8BDk3yI+BSYN+eTfYHTirTX0vSomGXBpIkSZI0wrxTJ0mSJEkjzEqdJEmSJI2wkch+ufXWW9fy5cuHHYYkqc/OP//8X1TVsmHHMSosHyVp6ZiqjByJSt3y5ctZvXr1sMOQJPVZkmuHHcMosXyUpKVjqjLS5peSJPVZkqOT3JjkkkmWvyTJRUkuTvL93i4IJEmajpU6SZL67xhg7ymWXw3sUVV/CrwXOGoQQUmSFoeRaH4pSdIoq6qzkiyfYvn3eybPBrbrd0ySpMXDO3WSJC0shwJfm2xhksOSrE6y+qabbhpgWJKkhcpKnSRJC0SSvWgqdW+dbJ2qOqqqVlTVimXLTBQqSbL5pTRjJ5xzXaf1Dtx9hz5HImkxSbIL8Fng2VV187DjkSTNzFS/Efv9u9A7dZIkDVmSHYAvAS+rqp8MOx5J0mjxTp0kSX2W5ERgT2DrJGuAdwH3A6iqTwPvBB4E/EsSgLurasVwopUkjRordZIk9VlVHTDN8lcArxhQOJKkRcbml5IkSZI0wqzUSZIkSdIIs1InSZIkSSPMSp0kSZIkjTArdZIkSZI0wqzUSZIkSdIIs1InSZIkSSPMSp0kSZIkjTArdZIkSZI0wqzUSZIkSdIIs1InSZIkSSPMSp0kSZIkjTArdZIkSZI0wqzUSZIkSdIIs1InSZIkSSPMSp0kSZIkjbC+VuqS/H2SS5NckuTEJBsn2THJOUmuTHJykg37GYMkSZIkLWbTVuqSPDzJRu34nklel2SLDtttC7wOWFFVjwbWB/YH3g98pKoeAfwKOHQuByBJ0ihIcnSSG5NcMsnyJPlYe9HzoiS7DjpGSdJo6nKn7t+Ae5I8AjgK2B44oeP+NwA2SbIBcH/gBuBpwKnt8mOB588oYkmSRtMxwN5TLH82sFM7HAZ8agAxSZIWgS6Vunur6m7gr4CPV9VbgG2m26iqrgf+GbiOpjJ3C3A+8Ot2fwBrgG1nE7gkSaOkqs4CfjnFKvsCx1XjbGCLJNOWt5IkdanU3ZXkAOAg4D/aefebbqMkW9IUUDsCDwEewNRXKMdvf1iS1UlW33TTTV03kyRpVG0L/LRn2gufkqROulTqDgGeBLyvqq5OsiNwfIftngFcXVU3VdVdwJeAJ9NcedygXWc74PqJNq6qo6pqRVWtWLZsWYe3kyRp8fOipyRpvGkrdVV1WVW9rqpObKevrqr3d9j3dcATk9w/SYCnA5cBZwL7tescBHxldqFLkrSoXE/z3PqYCS98etFTkjTeBpMtSHIxUJMtr6pdptpxVZ2T5FTgh8DdwAU0iVa+CpyU5J/aeZ+bRdySJC02K4HXJDkJ2B24papuGHJMkqQRMGmlDnjeXHdeVe8C3jVu9lXAbnPdtyRJg5bk9cDngd8CnwUeB7ytqr7RYdsTgT2BrZOsoSkf7wdQVZ8GTgOeA1wJ3E7z+IMkSdOatFJXVdcOMhBJkkbAy6vqyCTPArYEXkbznPm0lbqqOmCa5QW8el6ilCQtKV06H39ikvOS3JrkziT3JPnNIIKTJGmBSfv6HOD4qrq0Z54kSUPRJfvlJ4ADgCuATYBXAJ/sZ1CSJC1Q5yf5Bk2l7utJNgPuHXJMkqQlrkuljqq6Eli/qu6pqs8zg/7mJElaRA4F3gY8oapuBzbEZ98kSUM2VaKUMbcn2RC4MMkHgBvoWBmUJGmR+WZVPX1soqpuTnIKTbc9kiQNRZfK2cuA9YHXALfR9KHzgn4GJUnSQpJk4yRb0WSu3DLJVu2wHNh2uNFJkpa6ae/U9WTB/G/g3f0NR5KkBel/AG8AHgKcz++To/yG5tlzSZKGZtpKXZKrmaAT8qp6WF8ikiRpgamqI4Ejk7y2qj4+7HgkSerV5Zm6FT3jGwMvBLbqTziSJC1cVfXxJH8GLKenDK2q44YWlCRpyevS/PLmcbM+muR84J39CUmSpIUpyfHAw4ELgXva2QVYqZMkDU2X5pe79kyuR3PnrssdPkmSFpsVwM5Vtc5jCZIkDUuXytmHesbvBq4BXtSXaCRJWtguAf6QpnsfSZIWhC7NL/caRCCSJI2ArYHLkpwL3DE2s6r2GV5IkqSlbtJKXZI3TrVhVX14/sORJGlBO3zYAUiSNN5Ud+o2a18fBTwBWNlO/yVwbj+DkgbthHOuG3YIkkZAVX1n2DFIkjTepJW6qno3QJKzgF2r6rft9OHAVwcSnSRJC0iSJwIfB/4Y2BBYH7itqh441MAkSUvaeh3WeTBwZ8/0ne08SZKWmk8ABwBXAJsArwA+OdSIJElLXpfsl8cB5yb5cjv9fOCYvkUkSdICVlVXJlm/qu4BPp/kAuDtw45LkrR0TXunrqreBxwC/KodDqmq/9XvwCRJWoBuT7IhcGGSDyT5e7q1eiHJ3kkuT3JlkrdNsHyHJGcmuSDJRUmeM9/BS5IWp6myXz6wqn6TZCuavumu6Vm2VVX9sv/hSXNjAhRJ8+xlNM/RvQb4e2B74AXTbZRkfZpmms8E1gDnJVlZVZf1rPaPwClV9akkOwOnAcvnN3xJ0mI0VfPLE4DnAecD1TM/7fTD+hiXJEkLTlVd247+N/DuGWy6G3BlVV0FkOQkYF+gt1JXwFjClc2Bn80tWknSUjFV9svnta87Di4cSZIWniSnVNWLklzM2hc6AaiqXabZxbbAT3um1wC7j1vncOAbSV4LPAB4xuwjliQtJdMmSkmyEjgR+EpV3d7/kCRJWnBe374+r4/vcQBwTFV9KMmTgOOTPLqq7u1dKclhwGEAO+ywQx/DkSSNii4Pd38IeCrw4ySnJtkvycZ9jkuSpAWjqm5on4s7pqquHT902MX1NM/fjdmundfrUOCU9v1+AGwMbD1BLEdV1YqqWrFs2bJZHY8kaXHpkv3yO1X1Kppn6P4VeBFwY5edJ9mirQj+3yQ/TvKkJFsl+WaSK9rXLed2CJIk9V/bhcG9STafxebnATsl2bHNnrk/sHLcOtcBTwdI8sc0lbqb5hCyJGmJ6NJPHUk2Af4SeDGwK3Bsx/0fCZxeVfu1hdj9gf8JnFFVR7Qpnd8GvHXGkUuSNHi3Ahcn+SZw29jMqnrdVBtV1d1JXgN8nSZ75tFVdWmS9wCrq2ol8CbgM203CQUcXFXrPL8nSdJ4XZ6pO4Uma9fpwCeA74xv3z/JdpsDfw4cDFBVdwJ3JtkX2LNd7VhgFVbqJEmj4UvtMGNVdRpNNwW9897ZM34Z8OQ5RSdJWpK63Kn7HHBA2+xkJnakaTby+SSPoeka4fXAg6vqhnad/wIePMP9SpI0FFXVtaWKJEkDM1Xn40+rqm/TpFXeN8lay6tquiuVG9A01XxtVZ2T5Eiappa9+6gkEzYtMbuXRt1MOj4/cHe/49IoSLIT8L+AnWmeeQOgquy7VZI0NFPdqdsD+DbNs3TjFdM3P1kDrKmqc9rpU2kqdT9Psk2bSWwbJkm6UlVHAUcBrFixwmcKJEkLweeBdwEfAfYCDqFbJmlJkvpmqs7H39W+HjKbHVfVfyX5aZJHVdXlNBm9LmuHg4Aj2tevzGb/kiQNwSZVdUaStF0ZHJ7kfOCd020oSVK/TNX88o1TbVhVH+6w/9cCX2gzX17F769onpLkUOBami4SJEkaBXckWQ+4os1meT2w6ZBjkiQtcVM1v9ysfX0U8AR+35/OXwLndtl5VV0IrJhg0dO7BihJ0gLyeprueV4HvJemCeZBQ41IkrTkTdX88t0ASc4Cdq2q37bThwNfHUh0kiQtLPdU1a00/dXN6vEESZLmW5eHux8M3NkzfSd2QyBJWpo+lOTHSd6b5NHDDkaSJOjWT91xwLlJvtxOP5+m03BJkpaUqtoryR/SPA/+r0keCJxcVf805NAkSUvYtJW6qnpfktOBp7SzDqmqC/oblrS0dO3Tzv7spOGrqv8CPpbkTOAfaDJfWqmTJA1Nlzt1ABcCN4ytn2SHqures7IkSYtAkj8GXgy8ALgZOBl401CDkiQtedNW6pK8lqaj1Z8D9wCh6Xx8l/6GJknSgnM0cBLwrKr62bCDkSQJut2pez3wqKq6ud/BSJK0kFXVk4YdgyRJ43XJfvlT4JZ+ByJJkiRJmrkud+quAlYl+Spwx9jMqvpw36KSJEmSJHXS5U7ddcA3gQ2BzXoGSZKWhCTHt6+vH3YskiSN16VLg3cPIhBJkhawxyd5CPDyJMfRJA27T1X9cjhhSZI0RaUuyUer6g1J/p0m2+VaqmqfvkYmSdLC8WngDOBhwPmsXamrdr4kSUMx1Z2649vXfx5EIJIkLVRV9TGaDsc/VVWvnM0+kuwNHAmsD3y2qo6YYJ0XAYfTVBR/VFUHzj5qSdJSMWmlrqrOb1+/M7hwJElauKrqlUkeAzy1nXVWVV003XZJ1gc+CTwTWAOcl2RlVV3Ws85OwNuBJ1fVr5L8wfwfgSRpMeqSKEWSJAFJXgd8AfiDdvhCktd22HQ34Mqquqqq7qTpwHzfcev8LfDJqvoVQFXdOH+RS5IWsy5dGkiSpMYrgN2r6jaAJO8HfgB8fJrttqXp93XMGmD3ces8st3n92iaaB5eVaeP31GSw4DDAHbYYYdZHIIkabGZ9E6d6ZslSVpHgHt6pu9hXCbMOdgA2AnYEzgA+EySLcavVFVHVdWKqlqxbNmyeXprSdIom+pOnembJUla2+eBc5J8uZ1+PvC5DttdD2zfM71dO6/XGuCcqroLuDrJT2gqeefNLWRJ0mI3VaXO9M1akE4457phhyBpiaqqDydZBTylnXVIVV3QYdPzgJ2S7EhTmdsfGJ/Z8n/T3KH7fJKtaZpjXjUvgUuSFrWpsl/OOX2zJEmLTVX9EPjhDLe5O8lrgK/TPC93dFVdmuQ9wOqqWtku+4skl9E063xLVd08z+FLkhahaROlzDZ9syRJ+r2qOg04bdy8d/aMF/DGdpAkqbNpuzSYQ/pmSZIkSVKfdenSYLbpmyVJWjTaDsS/VVV7DTsWSZJ6danU9TN9s3QfE6BIWsiq6p4k9ybZvKpuGXY8kiSN6VKpm236ZuC+K5urgeur6nlt5q+TgAfRZNV8WVXdObOwJUkailuBi5N8E7htbGZVvW54IUmSlrouiVJmm755zOuBHwMPbKffD3ykqk5K8mngUOBTM9ifJEnD8qV2kCRpwehyp25W6ZsBkmwHPBd4H/DGJAGexu/75jkWOBwrdZKkEVBVxybZBNihqi4fdjySJEGH7Jdz9FHgH4B72+kHAb+uqrvb6TXAtn2OQZKkeZHkL4ELgdPb6ccmWTncqCRJS13fKnVJngfcWFXnz3L7w5KsTrL6pptumufoJEmalcOB3YBfA1TVhcDDhhmQJElTVuqSrJ/kzFnu+8nAPkmuoUmM8jTgSGCLJGPNPrcDrp9o46o6qqpWVNWKZcuWzTIESZLm1V0TZL68d8I1JUkakCkrdVV1D3Bvks1nuuOqentVbVdVy4H9gW9X1UuAM4H92tUOAr4y031LkjQklyY5EFg/yU5JPg58f9hBSZKWti6JUuY7ffNbgZOS/BNwATPoHkGSpCF7LfAO4A7gRODrwHuHGpEkacnrUqmbc/rmqloFrGrHr6J5HkGSpJFSVbcD70jy/mayfjvsmCRJ6tJPnembJUkCkjwBOBrYrJ2+BXj5bJOCSZI0H6bNfmn6ZkmS7vM54FVVtbx9ZvzVwOeHG5Ikaanr0qXB4Zi+WZIkgHuq6rtjE1X1f4C7p1hfkqS+6/JM3V1VdUuS3nmmb5aG4IRzruu87oG779DHSKSlJcmu7eh3kvwrTZKUAl5M+8y4JEnD0qVSt1b6ZuB1mL5ZkrS0fGjc9Lt6xmuQgUiSNF6XSp3pmyVJS1pV7TXXfSTZGzgSWB/4bFUdMcl6LwBOBZ5QVavn+r6SpMWvS/ZL0zdLkgQk2QL4G2A5PWXodH23Jlkf+CTwTGANcF6SlVV12bj1NgNeD5wzv5FLkhazLtkvn5DkYuAimk7If5Tk8f0PTZKkBec0mgrdxcD5PcN0dgOurKqrqupO4CRg3wnWey/wfuB38xKtJGlJ6NL8cix983cBkjyFJn3zLv0MTJKkBWjjqnrjLLbbFvhpz/QaYPfeFdpkLNtX1VeTvGWyHSU5DDgMYIcdTIgkSerWpYHpmyVJahyf5G+TbJNkq7FhrjtNsh7wYeBN061bVUdV1YqqWrFs2bK5vrUkaRGY9E6d6ZslSVrHncAHaRKIjWW9LKbvv/V6YPue6e3aeWM2Ax4NrGq7EPpDYGWSfUyWIkmazlTNL03fLEnS2t4EPKKqfjHD7c4DdkqyI01lbn/gwLGFVXULsPXYdJJVwJut0EmSupi0Ujcf6ZslSVpkrgRun+lGVXV3ktfQdAu0PnB0VV2a5D3A6qpaOc9xSpKWkGkTpcw2fbMkSYvQbcCFSc6k6b8V6FYmVtVpNNkze+e9c5J195xbmJKkpaRL9svTgLNp0jff299wJEla0P53O0iStGB0qdTNNn2zJEmLSlUdO+wYJEkar0ul7vgkfwv8B2s3Nfll36KSJGkBSnI1EyQLq6rpsl9KktQ3XSp1s03fLEnSYrOiZ3xj4IXAnPupkyRpLrpU6mabvlnihHOuG3YIkjRvqurmcbM+muR8YMKEJ5IkDUKXSt2s0jdLkrTYJNm1Z3I9mjt3XcpSSZL6pktBNOv0zZIkLTIf6hm/G7gGeNFwQpEkqdGlUmf6ZmkE9aPp64G77zDv+5RGSVXtNewYJEkab9pKnembJUlqJNkIeAGwnJ4ytKreM6yYJEmatlJn+mZJku7zFeAW4Hx6HkmQJGmYujS/nFX65iTbA8cBD6apFB5VVUcm2Qo4meYq5zXAi6rqVzMLW5KkodiuqvYedhCSJPVab7oVqurmnuH6qvoo8NwO+74beFNV7Qw8EXh1kp2BtwFnVNVOwBnttCRJo+D7Sf502EFIktSrS/PLWaVvrqobgBva8d8m+TGwLbAvsGe72rHAKuCtMwlakqQheQpwcPtowh1AgKqqXYYbliRpKevS/HLO6ZuTLAceB5wDPLit8AH8F03zTEmSRsGzhx2AJEnjdbnjNqf0zUk2Bf4NeENV/SZJ774ryTpJWNrtDgMOA9hhB9OoS5KGr6quHXYMkiSN16X55azTNye5H02F7gtV9aV29s+TbFNVNyTZBrhxom2r6ijgKIAVK1ZMWPGTJEmSpKVu2kQpNOmb96VpenlbzzClNLfkPgf8uKo+3LNoJXBQO35Qu39JkiRJ0ix0eaZutumbnwy8DLg4yYXtvP8JHAGckuRQ4Fpm+HyeJEmjKMnewJHA+sBnq+qIccvfCLyC5iLqTcDLbe4pSeqiS6Xu+0n+tKounsmOq+r/0GQFm8jTZ7IvSZJGWZL1gU8CzwTWAOclWVlVl/WsdgGwoqpuT/JK4APAiwcfrSRp1HRpfvkU4Pwklye5KMnFSS7qd2CSJC0iuwFXVtVVVXUncBLNow33qaozq+r2dvJsYLsBxyhJGlFd7tSZvlmSpLnZFvhpz/QaYPcp1j8U+FpfI5IkLRpdujSwPb8kSQOS5KXACmCPSZbb5Y8kaS1d7tRJEgAnnHNdp/UO3N0fmtI41wPb90xv185bS5JnAO8A9qiqOybakV3+SJLG6/JMnSRJmpvzgJ2S7JhkQ2B/mi5+7pPkccC/AvtU1YR9uEqSNBErdZIk9VlV3Q28Bvg68GPglKq6NMl7kuzTrvZBYFPgi0kuTLJykt1JkrQWm19KkjQAVXUacNq4ee/sGX/GwIOSJC0K3qmTJEmSpBFmpU6SJEmSRpiVOkmSJEkaYVbqJEmSJGmEmShFs9K1vzJpOvZ9J0mSNDfeqZMkSZKkEeadOknzrh93cmeyT+/qSZKkpcQ7dZIkSZI0wqzUSZIkSdIIs/mlpEXH5CvSzEz1N+PfiSQtfN6pkyRJkqQRZqVOkiRJkkaYlTpJkiRJGmFW6iRJkiRphFmpkyRJkqQRZqVOkiRJkkaYXRpIUgddu0noh36klJ/J8ZjSXpKkhW0olbokewNHAusDn62qI/r9nv34QTYqP3Tss0uShm+6si/JRsBxwOOBm4EXV9U1g45TkjR6Bl6pS7I+8EngmcAa4LwkK6vqskHHIknSIHQs+w4FflVVj0iyP/B+4MWDj3Ztk10Y9EKgJC0cw7hTtxtwZVVdBZDkJGBfwErdDM333cdhNi+ThsHvfDfe7Z8XXcq+fYHD2/FTgU8kSVXVIAPtaqrvhd8FSYvVQv3tMIxK3bbAT3um1wC7DyEOSZIGpUvZd986VXV3kluABwG/GEiE82ih/uiRpMVqwSZKSXIYcFg7eWuSy+dp11szTwXkS+ZjJ7Mzb8cwJMY/fKN+DKMeP8zgGIb4v2Yq68Q/T3E+dH52s3j1qXxcDH9T881zsi7Pybo8J+vynEzgJfNzXiYtI4dRqbse2L5nert23lqq6ijgqPl+8ySrq2rFfO93kEb9GIx/+Eb9GEY9fhj9Yxj1+IegS9k3ts6aJBsAm9MkTFlLP8pHP891eU7W5TlZl+dkXZ6TifX7vAyjn7rzgJ2S7JhkQ2B/YOUQ4pAkaVC6lH0rgYPa8f2Aby/U5+kkSQvLwO/Utc8JvAb4Ok1a56Or6tJBxyFJ0qBMVvYleQ+wuqpWAp8Djk9yJfBLmoqfJEnTGsozdVV1GnDaMN6bPjTpHIJRPwbjH75RP4ZRjx9G/xhGPf6Bm6jsq6p39oz/DnjhoONq+Xmuy3OyLs/Jujwn6/KcTKyv5yW27JAkSZKk0TWMZ+okSZIkSfNkUVXqkuyd5PIkVyZ52wTLd0hyZpILklyU5Dk9y97ebnd5kmcNNvL7YphV/EmemeT8JBe3r08bfPT3xTjrz6Bn+a1J3jy4qNd6/7l8h3ZJ8oMkl7afxcaDjX5O36H7JTm2jfvHSd4+6Nh7YpzuGB6a5Iw2/lVJtutZdlCSK9rhoPHbDsJs40/y2J7vz0VJXjz46O+LcdafQbv8gUnWJPnE4KLWdDp8rhslObldfk6S5YOPcvA6nJc3Jrms/b6fkWTRd7sx3TnpWe8FSSrJos902OWcJHlR+125NMkJg45x0Ob6m28xSnJ0khuTXDLJ8iT5WHvOLkqy67y9eVUtioHmwfP/BB4GbAj8CNh53DpHAa9sx3cGrukZ/xGwEbBju5/1Ryj+xwEPaccfDVw/ap9Bz/JTgS8Cbx6l+GmeT70IeEw7/aAR+w4dCJzUjt8fuAZYvkA/gy8CB7XjTwOOb8e3Aq5qX7dsx7ccofgfCezUjj8EuAHYYpQ+g57lRwInAJ8YdPwOc/pcXwV8uh3fHzh52HEvkPOyF3D/dvyVi/28dDkn7XqbAWcBZwMrhh33sM8JsBNwwVi5A/zBsONeAOdkyt98i3EA/hzYFbhkkuXPAb4GBHgicM58vfdiulO3G3BlVV1VVXcCJwH7jlungAe245sDP2vH96X5QXtHVV0NXNnub5BmHX9VXVBVY8dyKbBJko0GEPN4c/kMSPJ84GqaYxiGucT/F8BFVfUjgKq6uaruGUDMveYSfwEPSNM31ibAncBv+h/yOrocw87At9vxM3uWPwv4ZlX9sqp+BXwT2HsAMfeadfxV9ZOquqId/xlwI7BsIFGvbS6fAUkeDzwY+MYAYlV3XT7XfYFj2/FTgacnyQBjHIZpz0tVnVlVt7eTZ9P0MbiYdfmuALwXeD/wu0EGNyRdzsnfAp9syx+q6sYBxzhoc/rNt1hV1Vk02Ysnsy9wXDXOBrZIss18vPdiqtRtC/y0Z3pNO6/X4cBLk6yhyUD22hls229zib/XC4AfVtUd/QhyGrM+hiSbAm8F3t3/MCc1l8/gkUAl+XqSHyb5h34HO4G5xH8qcBvN3aHrgH+uqqn+KfVLl2P4EfDX7fhfAZsleVDHbfttLvHfJ8luNFc+/7NPcU5l1seQZD3gQ8BQmk9rSl0+1/vWqaq7gVtoWh0sZjP9v3EozVX2xWzac9I2Gdu+qr46yMCGqMv35JHAI5N8L8nZSQZ9UXHQ5ut361LTt98qi6lS18UBwDFVtR3N7c/j2x8ho2LK+JP8Cc1Vs/8xpPi6mOwYDgc+UlW3DjO4DiaLfwPgKcBL2te/SvL04YU5qcni3w24h6bZ347Am5I8bHhhTunNwB5JLgD2AK6niX1UTBl/e8XueOCQqrp3OCFOa7JjeBVwWlWtGWZwUj8keSmwAvjgsGMZprbM+DDwpmHHssBsQNMEc0+asvYzSbYYakTDN+q/u0fKUPqp65Prge17prdr5/U6lLY5VlX9IE0ii607bttvc4n/xjZRwZeBv6mqYVzdh7kdw+7Afkk+AGwB3Jvkd1U1yEQLc4l/DXBWVf0CIMlpNG2qz+h30D3mEv+BwOlVdRfN9+l7ND9erl1VPFcAAAjcSURBVOp71Gub9hjapol/Dffd4X1BVf06yfU0hWnvtqv6GewEZh1/O/1A4KvAO9pmGcMwl8/gScBTk7wK2BTYMMmtVTVpogUNTJf/D2PrrGmbYm8O3DyY8IamU/mf5BnAO4A9htQSZpCmOyeb0Ty/v6ptnfuHwMok+1TV6oFFOVhdvidraJ6Pugu4OslPaCp55w0mxIGb0+/WgUS4MPWvztGPhwSHMdBUUK+iucsw9sDmn4xb52vAwe34H9O07Q3wJ6ydKOUqBp/kYi7xb9Gu/9ej+hmMW+dwhpMoZS6fwZbAD2mSjGwAfAt47gjF/1bg8+38BwCXAbss0M9ga2C9dvx9wHva8a1onsncsh2uBrYaofg3pLkI8IZBn/f5OoZx6xyMiVIWzNDxc301aydKOWXYcS+Q8/I4mqbQOw073oVyTsatv4rFnyily/dkb+DYdnxrmiZ2Dxp27EM+J9P+5luMA7CcyROlPJe1E6WcO2/vO+wDn+eT+BzgJ+0/33e0894D7NOO7wx8r/3iXQj8Rc+272i3uxx49ijFD/wjzfNQF/YMQ8m6NJfPoGcfhzOESt08fIdeSpPk5RLgA6MUP81dlS+28V8GvGUY8Xc8hv2AK9p1Pgts1LPty2kSHV1J03xxZOJvvz93jfs7fuwoHcO4fRyMlboFNXT4XDdu/w9cCZwLPGzYMS+Q8/It4Oc9f5crhx3zsM/JuHVXscgrdR2/J6FplnoZcDGw/7BjXgDnZNrffIttAE6kyU9wF83d20OBvwP+rud78sn2nF08n387ad9AkiRJkjSCfFhRkiRJkkaYlTpJkiRJGmFW6iRJkiRphFmpkyRJkqQRZqVOkiRJkkaYlTotWknekOT+87SvP0pyYZILkjx8FtuvSPKxWWx3WpItZrpdPyT5/hy2PTjJQ+YzHknS1CwH19rPQ5KcOsNtDk7yiXb875L8zTTrfzbJzlPtZzaS3DrbbbV02KWBRlaS0HyH751k+TU0/X/8Yh7e623ABlX1T3Pd16hJskFV3T3Hfayi6Xtw9fxEJUmyHOyvJAfTnL/XDHM/SW6tqk3nEoMWP+/UaaQkWZ7k8iTH0XTyvX2STyVZneTSJO9u13sd8BDgzCRntvP+IskPkvwwyReTrPMPMsljk5yd5KIkX06yZZLnAG8AXjm2r3Hb3Jrkg+37fyvJbklWJbkqyT7tOnsm+Y92fI/2aufYFc/NkmyT5Kx23iVJntque02Srdvj/nGSz7Tv840km7TrPKGN98I2jksmiPGkJM/tmT4myX7tfr/bnpMfJvmznni/m2QlTUeq910pTLJpkjPa9S9Osm/PZ7NOjEn2A1YAX2hj3CTJEUkua+P+51l9GSRpCbIcnHU5uHxsfpo7Z19KcnqSK5J8oGe9Q5L8JMm5wJN75h+e5M1p7lieO26/F7fjq5KsmGY/x7Tl4n3nrn2dsGwddwwTniMJYOg9rzs4zGQAlgP3Ak/smbdV+7o+sArYpZ2+Bti6Hd8aOAt4QDv9VuCdE+z/ImCPdvw9wEfb8cNp7jRNFFMBz27Hvwx8A7gf8Bjgwnb+nsB/tOP/Djy5Hd8U2AB4E/COnuPYrPcY2uO+G3hsO/8U4KXt+CXAk9rxI4BLJojxr4Bj2/ENgZ8CmwD3BzZu5+8ErO6J9zZgx5593Nq+bgA8sOe8XglkmhhX0VylBHgQcDm/bymwxbC/Vw4ODg6jMlgOzrocXD42HzgYuArYHNgYuBbYHtgGuA5Y1paV3wM+Mf74gQvHysf2PP5jO76K5iLmVPs5BtivJ64py9Zx60x4jhwcqso7dRpJ11bV2T3TL0ryQ+AC4E+AddqzA09s538vyYXAQcBDe1dIsjlNBeM77axjgT/vEM+dwOnt+MXAd6rqrnZ8+QTrfw/4cHsVdYtqmjaeBxyS5HDgT6vqtxNsd3VVXdiOnw8sT/OcwWZV9YN2/gmTxPg1YK8kGwHPBs6qqv+mKXQ/015l/CJrn7tzq+rqCfYV4P9PchHwLWBb4MGTxTjB9rcAvwM+l+SvgdsniVmSNDHLwZmXg+OdUVW3VNXvaFqkPBTYHVhVVTdV1Z3AyZNsewrw4nb8xROs13U/vaYqW8d0OUdaoqzUaRTdNjaSZEfgzcDTq2oX4Ks0V93GC/DNqnpsO+xcVYfOUzx3VdXYw6n3AncAVPOMwwbjV66qI4BX0Nwp+16SP6qqs2gKzuuBYzLxw9h39IzfM9G+J9MWWquAZ7F2AfT3wM9prqauoLmiOOY2JvYSmquPj6+qx7bbj53zaWNsC+/dgFOB5/H7HwKSpG4sB2dYDs7zvk6mqUg/EqiqumIG295N+/s7yXr8vtydqmyF5o26nCMtUVbqNOoeSFO43ZLkwTR3ocb8FtisHT8beHKSRwAkeUD7z/g+VXUL8KueNuovA77DPEvy8Kq6uKreT3PV7Y+SPBT4eVV9BvgssGuXfVXVr4HfJtm9nbX/FKufDBwCPJXfV6Q2B25oC96X0TTnmM7mwI1VdVeSvRh3pXcS930W7TMcm1fVaTSVysd02F6SNDHLwe7l4HTOAfZI8qAk9wNeOMl7/idNRfD/Y+K7cFPt5xrg8e34PjQtZqBD2Trbc6SlYS5XOKShq6ofJbkA+L80z4l9r2fxUcDpSX5WVXulyT51YtsEEeAfgZ+M2+VBwKfTpIC+iqYSNN/e0P7Dvhe4lKZp5P7AW5LcBdwKzOTq26E0TSjvpSl8b5lkvW8AxwNfaZuDAPwL8G/t1b7TmfzuXK8vAP/eNtlcTXPup3MMzXn9b5ofHF9JsjHNleM3dthekjQBy0Ggezk4paq6oW3a+APg1zTPzk3mZOCDwI4z3M9naMrAH7F2udulbN2T2Z8jLXJ2aSCNuCSbVtVY9qy3AdtU1euHHJYkSQNhOSh5p05aDJ6b5O00f8/X0mT1kiRpqbAc1JLnnTpJkiRJGmEmSpEkSZKkEWalTpIkSZJGmJU6SZIkSRphVuokSZIkaYRZqZMkSZKkEWalTpIkSZJG2P8DyQQFK1ynSbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_missing(df_miss_ind, df_miss_var)\n",
    "only_missing = df_miss_var.F_MISS[df_miss_var.F_MISS == 1.0].index\n",
    "print(\"{:.3}% variants with missing rate of 100%\"\n",
    "      .format(len(only_missing)/df_miss_var.shape[0]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### First report of missing values\n",
    "\n",
    "Looking at the plot on the right, we see the many variants have 100% of missing individuals. This is expected since the raw data has changed from the beginning. The provided data is now more adapted to exome sequencing data, where variants of poor sequencing quality (ex some intronic regions that are close to exonic regions) have been set to only missing values for each individual. \n",
    "\n",
    "This allows us to remove variants with 100% missingness before looking at the NaN's again, since missingness is biased here by those bad quality variants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intermediate data\n",
    "\n",
    "We first make a new look at the missing report *after* removing variants with 100% missingness (we can't simply drop the variants in the DataFrame since the plot of showing the rate of missing variants per individuals would be biased, we have to run the plink missing report again).\n",
    "\n",
    "*Note: we actually would like to run `--geno 1.0`, but plink tests if the missingness is greater than the provided value (and not greater or equal). Thus we provide a value that is very close to 1.0.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink2 --bfile data/raw/wes_plink/hbv_gilead_bcftools --out data/plink/host_geno_clean_intermed --geno 0.999999999 --make-bed '\n"
     ]
    }
   ],
   "source": [
    "path_intermed=path_out+file_prefix+'_intermed'\n",
    "o, e= run_plink('--geno 0.999999999 --make-bed', file=path_raw, out=path_intermed,\n",
    "         extension='bed', force=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log reports 3'661'936 remaining variants (11'076'524 removed).\n",
    "\n",
    "## New plink missing report\n",
    "\n",
    "At this point we have a dataset in which there is no variant with 100% missingness. We want to have an unbiased look at this new dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink2 --bfile data/plink/host_geno_intermed --out data/plink/host_geno_clean_missing --missing '\n",
      "Log written to 'data/plink/log/missing.log'\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAADQCAYAAACgN0jgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5hsVXnn8e9PEEHlKoyjXHJQiYYkavAoZMxF4jyKl4AzXoIkCgR14t1oMuIkUdQ4o0m8oUaDglweEQjR8URRJAqSGAEPQkRgDCeAXIJKQFFAub7zx14txaG7end3Vdep7u/nefZTe6/al3fvU6dXrdprvTtVhSRJkiRpOt1v0gFIkiRJkhbPRp0kSZIkTTEbdZIkSZI0xWzUSZIkSdIUs1EnSZIkSVPMRp0kSZIkTbHNJx1AHzvuuGOtWbNm0mFIksbs/PPP/4+q2mnScUwL60dJWj2G1ZFT0ahbs2YN69evn3QYkqQxS/KdSccwTawfJWn1GFZH2v1SkiRJkqaYjTpJkiRJmmI26iRJkiRpitmokyRJkqQpZqNOkiRJkqbYVGS/lDYlJ557Va/1Dtp7tzFHIkmjMezvmn/LJGnT5506SZIkSZpiNuokSZIkaYrZqJMkSZKkKWajTpIkSZKmmI06SZIkSZpiNuokSZIkaYrZqJMkSZKkKWajTpIkSZKmmI06SZIkSZpiNuokSZIkaYrZqJMkSZKkKWajTpIkSZKmmI06SZIkSZpiNuokSRqQ5A+TXJzkW0k+mWTLJLsnOTfJhiQnJ9mirfuAtryhvb9mYD9vauXfTvL0gfL9WtmGJIcPlM96DEmS5mOjTpKkJsnOwGuAtVX1S8BmwIHAu4D3VtWjgB8Ah7VNDgN+0Mrf29YjyZ5tu18E9gP+OslmSTYDPgQ8A9gTeGFblyHHkCRpKBt1kiTd2+bAVkk2Bx4IXAf8FnBqe/844Dlt/oC2THv/qUnSyk+qqtuq6gpgA/CkNm2oqsur6nbgJOCAts1cx5AkaSgbdZIkNVV1LfBXwFV0jbmbgPOBH1bVnW21a4Cd2/zOwNVt2zvb+g8ZLN9om7nKHzLkGJIkDWWjTpKkJsn2dHfZdgceDjyIrvvkJiPJy5KsT7L++uuvn3Q4kqRNgI06SZLu8V+BK6rq+qq6A/gU8GRgu9YdE2AX4No2fy2wK0B7f1vghsHyjbaZq/yGIce4l6o6qqrWVtXanXbaaSnnKklaIcbaqFtIBjFJkjYBVwH7JHlgG+f2VOAS4EzgeW2dg4HPtPl1bZn2/perqlr5gS075u7AHsB5wNeBPVpduAVdMpV1bZu5jiFJ0lCbz7dCkkcC11TVbUmeAjwWOL6qfjjPdjMZxPasqp8kOYWu8nomXXavk5J8hC6714eXeB6SJC1ZVZ2b5FTgG8CdwAXAUcDngJOS/HkrO7ptcjRwQpINwI109RxVdXGr9y5p+3llVd0FkORVwOl0mTWPqaqL277eOMcxJElT4MRzr5rzvYP23m2sx563UQf8HbA2yaPoKrbPACfSNc767H+rJHdw7wxiB7X3jwOOwEadJGkTUVVvAd6yUfHldJkrN173p8Dz59jPO4B3zFJ+GnDaLOWzHkOSpPn06X55d8vG9d+AD1TVHwMPm2+jRWQQkyRJkiQtUJ9G3R1JXkjXv/+zrez+82201AxiZveSJEmSpPn1adQdCvwq8I6quqIN+D6hx3YLzSB2L2b3kiRJkqT5zduoq6pLquo1VfXJtnxFVb2rx74XmkFMkiRJkrRAcyZKSXIRUHO9X1WPHbbjRWQQkyRJkiQt0LDsl89e6s4XkkFMkqRRSvJa4OPAj4GPAb8CHF5VX5xoYJIkjdicjbqq+s5yBiJJ0oj9flW9P8nTge2BF9GNCbdRJ0laUeYdU5dknyRfT3JzktuT3JXkR8sRnCRJS5D2+kzghPaQ7wxZX5KkqdQn++UHgRcClwFbAS8BPjTOoCRJGoHzk3yRrlF3epKtgbsnHJMkSSPXp1FHVW0ANququ6rq4yzgeXOSJE3IYcDhwBOr6lZgC7rH9EiStKIMS5Qy49YkWwAXJvkL4Dp6NgYlSZqgM6rqqTMLVXVDklPoHrEjSdKK0adx9iJgM+BVwC3ArsBzxxmUJEmLlWTLJDsAOybZPskObVoD7DzZ6CRJGr1579QNZMH8CfDW8YYjSdKS/Q/gdcDDgfO5JznKj+jGiUuStKLM26hLcgWzPIS8qh4xlogkSVqCqno/8P4kr66qD0w6HkmSxq3PmLq1A/NbAs8HdhhPOJIkjUZVfSDJfwHWMFDfVdXxEwtKkqQx6NP98oaNit6X5HzgzeMJSZKkpUtyAvBI4ELgrlZcgI06SdKK0qf75V4Di/eju3PX5w6fJEmTtBbYs6ruM4RAkqSVpE/j7N0D83cCVwIvGEs0kiSNzreA/0z3KB5JklasPt0v912OQCRJGrEdgUuSnAfcNlNYVftPLiRJkkZvzkZdktcP27Cq3jP6cCRJGpkjJh2AJEnLYdiduq3b66OBJwLr2vJvA+eNMyhpuZ147lWTDkHSiFXVVyYdgyRJy2HORl1VvRUgydnAXlX147Z8BPC5ZYlOkqRFSrIP8AHgF4AtgM2AW6pqm4kGJknSiN2vxzoPBW4fWL69lUmStCn7IPBC4DJgK+AlwIcmGpEkSWPQJ/vl8cB5ST7dlp8DHDu2iCRJGpGq2pBks6q6C/h4kguAN006LkmSRqlP9st3JPk88Out6NCqumC8YUmStGS3JtkCuDDJX9A92qBPDxVJkqbKsOyX21TVj5LsQPdsuisH3tuhqm4cf3jS0pgARVrVXkQ3ju5VwB8CuwLPnWhEkiSNwbBfLE9sr+cD6wemmWVJkjZZVfWdqvpJVf2oqt5aVa+vqg3zbZdkuySnJvl/SS5N8qtJdkhyRpLL2uv2bd0kOTLJhiTfTLLXwH4ObutfluTggfInJLmobXNkkrTyWY8hSdJ85mzUVdWz2+vuVfWIgWn3qnrE8oUoSVJ/SU5prxe1hta9ph67eD/whap6DPA44FLgcOBLVbUH8KW2DPAMYI82vQz4cDv2DsBbgL2BJwFvGWikfRh46cB2+7XyuY4hSdJQ846pS7IO+CTwmaq6dfwhSZK0JK9tr89e6IZJtgV+AzgEoKpuB25PcgDwlLbaccBZwBuBA4Djq6qAc9pdvoe1dc+YGaqQ5AxgvyRnAdtU1Tmt/Hi6BGSfb/ua7RiSJA3VZ8D4u+mSpFzauqM8L8mWY45LkqRFqarrkmwGHNu6YN5rmmfz3YHraZkyk3wsyYOAh1bVdW2d73LPo312Bq4e2P6aVjas/JpZyhlyjHtJ8rIk65Osv/766+c5HUnSajBvo66qvlJVrwAeAfwN8ALg+312vpBxCZIkjUp7hMHd7c7bQmwO7AV8uKp+BbiFjbpBtrtyNZJA5zDsGFV1VFWtraq1O+200zjDkCRNiV6pnZNsRZcx7A+AJ9J1C+ljIeMSJEkapZuBi5Ic3RKSHJnkyHm2uQa4pqrObcun0jXyvte6VdJeZ37cvJYuq+aMXVrZsPJdZilnyDEkSRpq3kZdG3B+KfBbwAeBR1bVq3tsNzMu4WjoxiVU1Q/pxgzMNAqPoxtLIEnSqH0K+DPgbLrMzTPTnKrqu8DVSR7dip4KXAKsA2YyWB4MfKbNrwNe3LJg7gPc1LpQng48Lcn2rUfK04DT23s/SrJPy3r54o32NdsxJEkaat5EKXSNshe2riwLMTgu4XF0Felr6TlmQJKkpaiqvr1KNvZq4BPtweWXA4fS/Qh6SpLDgO/QDUUAOA14JrABuLWtS1XdmOTtwNfbem8beL7rK4Bjga3oEqR8vpW/c45jSJI01LCHj/9WVX0ZeBBwQHuMzs9U1ad67Hsv4NVVdW6S9zPLuIQks44ZSPIyuvTQ7LbbbvOdh7TJWciDzw/a28+4NGpJ9gD+D7An8LMEX/M9lqeqLgTWzvLWU2dZt4BXzrGfY4BjZilfD/zSLOU3zHYMSZLmM+xO3W8CXwZ+e5b3iq5byzCzjUs4nDZmoGUnm3PMQFUdBRwFsHbt2rEOSJckrUgfp3tW3HuBfbnnjpskSSvKnI26qnpLez10MTuuqu8muTrJo6vq29wzLuESurEC78QxA5Kk8dmqqr6UJO1RBkckOR9486QDkyRplIZ1v3z9sA2r6j099r+QcQmSJI3SbUnuB1yW5FV0WSYfPOGYJEkauWHdL7dur4+me4zBurb828B5fXa+kHEJkiSN2GuBBwKvAd5O1wXz4KFbSJI0hYZ1v3wrQJKzgb2q6sdt+Qjgc8sSnSRJi3dXVd1M97y6RQ0lkCRpGvQZMP5Q4PaB5dvxMQSSpE3fu5NcmuTtSe6TbVKSpJWiz3PqjgfOS/Lptvwc7nl4uCRJm6Sq2jfJf6Ybu/03SbYBTq6qP59waJIkjdS8jbqqekeSLwC/1ooOraoLxhuWtLr0faadz7OTFqaqvgscmeRM4H/SZb60USdJWlH63KkDuBC4bmb9JLtVVf8nK0uStMyS/ALwO8BzgRuAk4E3TDQoSZLGYN5GXZJX0z289XvAXUDoHj7+2PGGJknSkhwDnAQ8var+fdLBSJI0Ln3u1L0WeHRV3TDuYCRJGpWq+tVJxyBJ0nLok/3yauCmcQciSZIkSVq4PnfqLgfOSvI54LaZwqp6z9iikiRJkiT10udO3VXAGcAWwNYDkyRJm5wkJ7TX1046FkmSlkOfRxq8dTkCkSRpRJ6Q5OHA7yc5ni7B189U1Y2TCUuSpPGYs1GX5H1V9bokf0+X7fJeqmr/sUYmSdLifAT4EvAI4Hzu3airVi5J0oox7E7dCe31r5YjEEmSRqGqjqR74PiHq+rlk45HkqRxm7NRV1Xnt9evLF84kiSNRlW9PMnjgF9vRWdX1TcnGZMkSePQJ1GKJElTJ8lrgE8A/6lNn0jy6slGJUnS6PV5pIEkSdPoJcDeVXULQJJ3AV8DPjDRqCRJGrE579SZElqSNOUC3DWwfBcbZcKUJGklGHanzpTQkqRp9nHg3CSfbsvPAY6eYDySJI3FsEadKaG1STrx3KsmHYKkKVBV70lyFvBrrejQqrpggiFJkjQWw7JfmhJakjTVquobwDcmHYckSeM0b6IUU0JLkiRJ0qZr3kcamBJakrTaJNksyQVJPtuWd09ybpINSU5OskUrf0Bb3tDeXzOwjze18m8nefpA+X6tbEOSwwfKZz2GJEnz6fOcupmU0G+uqjcD+wAvHW9YkiQtXmuUnbmEXbwWuHRg+V3Ae6vqUcAPgMNa+WHAD1r5e9t6JNkTOBD4RWA/4K9bTJsBHwKeAewJvLCtO+wYkiQN1ec5daaE1rIwAYqkUamqu5LcnWTbqrppIdsm2QV4FvAO4PVJAvwWcFBb5TjgCODDwAFtHuBU4INt/QOAk6rqNuCKJBuAJ7X1NlTV5e1YJwEHJLl0yDEkSRqqT6NuSSmh26+S64Frq+rZSXYHTgIeQpdV80VVdfvCwpYkaV43AxclOQO4Zaawql4zz3bvA/4nsHVbfgjww6q6sy1fA+zc5ncGrm77vTPJTW39nYFzBvY5uM3VG5XvPc8xJEkaat7ul1X1HuBQ4MY2HVpV71vAMfp2YZEkaZQ+BfwZcDbdj4gz05ySPBv4flUNXW+Skrwsyfok66+//vpJhyNJ2gT0uVO36JTQC+zCIknSyFTVcUm2Anarqm/33OzJwP5JnglsCWwDvB/YLsnm7U7aLsC1bf1rgV2Ba5JsDmwL3DBQPmNwm9nKbxhyjI3P6yjgKIC1a9dWz/OSJK1gfRKlLMVMF5a727LdSyRJyyLJbwMXAl9oy49Psm7YNlX1pqraparW0CU6+XJV/S5wJvC8ttrBwGfa/Lq2THv/y1VVrfzAlh1zd2AP4Dzg68AeLdPlFu0Y69o2cx1DkqShxtaoW2oXFruXSJKW6Ai65CQ/BKiqC4FHLHJfb6TrcbKB7gfKmbHlRwMPaeWvBw5vx7oYOAW4hK5R+cqquqv9qPkq4HS6oQmntHWHHUOSpKGGdr9sSU7+oar2XcS+F9qF5V7sXiJJWqI7quqmruf/z9w918obq6qzgLPa/OXck71ycJ2fAs+fY/t30A0/2Lj8NOC0WcpnPYYkSfMZeqeuqu4C7k6y7UJ3vIguLJIkjdLFSQ4CNkuyR5IPAP886aAkSRq1PolSFpsSei5vBE5K8ufABdi9RJI0Hq8G/gS4DfgkXZfHt080IkmSxqBPo+5TbVq0Pl1YJEkapaq6FfiTJO/qFuvHk45JkqRxmLdRt8iU0JIkTVSSJwLH0B4i3h4M/vub8jPoJElajHmzXy4mJbQkSZuAo4FXVNWaNr77lcDHJxuSJEmj1+eRBkcwupTQkiQtl7uq6h9nFqrqn4A7h6wvSdJU6jOmbkkpoSWNzonnXtV73YP23m2MkUibriR7tdmvJPkbuiQpBfwObXy3JEkrSZ9G3b1SQgOvwZTQkqRN17s3Wn7LwLzPPZUkrTh9GnWmhJYkTY2q2nfSMUiStJz6ZL80JbQkaeok2Q54MbCGgfpuCc9ZlSRpkzRvo86U0JKkKXUacA5wEY4FlyStYH26X86khP5HgCS/RpcS+rHjDEySpCXasqpeP+kgJEkatz6PNDAltCRpGp2Q5KVJHpZkh5lp0kFJkjRqc96pMyW0JGnK3Q78JV2yr5msl4XPWpUkrTDDul+aElqSNM3eADyqqv5j0oFIkjROczbqTAktSZpyG4BbJx2EJEnj1if7pSmhJUnT6BbgwiRn0j1rFbD+kiStPH2yX5oSWpI0jf5vmyRJWtH6NOpMCS1JmjpVddykY5AkaTn0adSdkOSlwGe5d/eVG8cWlSRJS5TkCmZJ7FVVZr+UJK0ofRp1poSWJE2jtQPzWwLPB3xOnSRpxenTqDMltBbtxHOvmnQIklapqrpho6L3JTkfePMk4pEkaVz6NOpMCS1JmjpJ9hpYvB/dnbs+9Z4kSVOlT+VmSmhJ0jR698D8ncCVwAsmE4okSePTp1FnSmhpCo2j6+tBe+828n1K41JV+046BkmSlsO8jTpTQkuSplGSBwDPBdYwUN9V1dsmFZMkSeNwv/lWSHJFkss3npYjOEmSluAzwAF0XS9vGZjmlGTXJGcmuSTJxUle28p3SHJGksva6/atPEmOTLIhyTcHx/ElObitf1mSgwfKn5DkorbNkUky7BiSJM2nT/fLRaWETrIrcDzwULpHIBxVVe9PsgNwMt0vp1cCL6iqHywsbEmS5rVLVe23wG3uBN5QVd9IsjVwfpIzgEOAL1XVO5McDhwOvBF4BrBHm/YGPgzs3eq6t9DVodX2s67Vdx8GXgqcC5wG7Ad8vu1ztmNIkjTUvHfqquqGgenaqnof8Kwe+56pGPcE9gFemWRP7qm09gC+1JYlSRq1f07yywvZoKquq6pvtPkfA5cCO9Pd8ZsZjnAc8Jw2fwBwfHXOAbZL8jDg6cAZVXVja8idAezX3tumqs6pqqL78XNwX7MdQ5Kkoea9U7fYlNBVdR1wXZv/cZLBivEpbbXjgLPwl0hJ0uj9GnBIkivosjcHqKp6bJ+Nk6wBfoXujtpDW70G8F26XijQ1WtXD2x2TSsbVn7NLOUMOcbGcb0MeBnAbruNP3nRXEmXTJwkSZuOPt0vl5wSumfFKEnSKD1jsRsmeTDwd8DrqupHbdgb0LUKk9QI4pvTsGNU1VHAUQBr164daxySpOnQ547bklJCL7ZiXO5fIiVJK0tVfWcx2yW5P1299Ymq+lQr/l6Sh1XVda0L5fdb+bXArgOb79LKruWeXikz5We18l1mWX/YMSRJGqpP9ssHJDkoyf9K8uaZqc/Oh1WM7f05K62qOqqq1lbV2p122qnf2UiStAQtE+XRwKVV9Z6Bt9YBMxksD6bLrDlT/uKWBXMf4KbWG+V04GlJtm9ZLJ8GnN7e+1GSfdqxXrzRvmY7hiRJQ/XpfvkZ4CbgfLoxCb30qBjfiZWWJGnT8mTgRcBFSS5sZf+Lrs46JclhwHe4ZxjCacAzgQ3ArcChAFV1Y5K3A19v672tqm5s868AjgW2ost6+flWPtcxJEkaqk+jbjEpoWHhFaMkSRNVVf9El1BlNk+dZf0CXjnHvo4BjpmlfD3wS7OU3zDbMSRJmk+fRt0/J/nlqrpoITteaMUoSZIkSVq4Po26JaWEliRJkiSNT59G3aJTQkuSJEmSxqvPIw0WlRJakiRJkjR+fe7USRIAJ557Va/1DtrbZ0tKkiQtl3mfUydJkiRJ2nTZqJMkSZKkKWb3S0mStGDDumPbBVuSlpd36iRJkiRpitmokyRJkqQpZqNOkiRJkqaYY+okSZIkqYe+j3dabjbqtCib6gda08dn30mSJC2N3S8lSZIkaYp5p07SyI3jTu5C9uldPWmyfNyBJC0v79RJkiRJ0hTzTp0kSVo23sWTtKmbxtwRNuokrTgmX5EkSauJ3S8lSZIkaYp5p06SJG0S7JopSYtjo06SJG3ybPBJGqVpHDc3jI06SZI01WzwSZrNSmu4DWOjTpIkrViL/VJnY1CaDqup4TaMjTpJkqSNjOOLog1FrQaLvXNu42xpbNRJUg+TrGzG8UVwIefjF1FpNKb9S+tiv5CPY7th5trncv8tHXX88+1zGj5f0xDjtJpIoy7JfsD7gc2Aj1XVO8d9zNX8i5vP7JKk6TGJOlLqY7HfpZZ7u+Xa3ySOZ6NIc1n259Ql2Qz4EPAMYE/ghUn2XO44JEna1FhHSpIWYxJ36p4EbKiqywGSnAQcAFwygVim2rT/giVNmp/5frzbv6ysIyVJC7bsd+qAnYGrB5avaWWSJK121pGSpAXbZBOlJHkZ8LK2eHOSb4/5kDsC/7GQDX53TIEskwWf7xRbTecKnu+Ks9Hfmk32fEf0N/HnRrOblWtM9eMm+7maIK/JfXlNZud1uS+vyUZ+dzTXZM46chKNumuBXQeWd2ll91JVRwFHLVdQSdZX1drlOt6krabzXU3nCp7vSrfazncVmreOHEf96Ofqvrwm9+U1mZ3X5b68Jvc17msyie6XXwf2SLJ7ki2AA4F1E4hDkqRNjXWkJGnBlv1OXVXdmeRVwOl06ZqPqaqLlzsOSZI2NdaRkqTFmMiYuqo6DThtEsceYtm6em4iVtP5rqZzBc93pVtt57vqTKiO9HN1X16T+/KazM7rcl9ek/sa6zVJVY1z/5IkSZKkMZrEmDpJkiRJ0oisqkZdkv2SfDvJhiSHz/L+65NckuSbSb6UZKpTa/c43z9IclGSC5P8U5I9JxHnqMx3vgPrPTdJJZnqrEw9/n0PSXJ9+/e9MMlLJhHnqPT5903ygvZ/+OIkJy53jKPU49/3vQP/tv+a5IeTiFPTo8dn6gFJTm7vn5tkzfJHufxW23eDPlZbfdrHaquD+urx/2e3JGcmuaD9H3rmJOJcLkmOSfL9JN+a4/0kObJdr28m2WtkB6+qVTHRDTj/N+ARwBbAvwB7brTOvsAD2/zLgZMnHfeYz3ebgfn9gS9MOu5xnm9bb2vgbOAcYO2k4x7zv+8hwAcnHesynu8ewAXA9m35P0067nGe70brv5ouocbEY3faNKee/4deAXykzR84zXXgiK/LivluMKpr0tZbEfXpCD8nK6YOGvF1OQp4eZvfE7hy0nGP+Zr8BrAX8K053n8m8HkgwD7AuaM69mq6U/ckYENVXV5VtwMnAQcMrlBVZ1bVrW3xHLrnA02rPuf7o4HFBwHTPMBy3vNt3g68C/jpcgY3Bn3Pd6Xoc74vBT5UVT8AqKrvL3OMo7TQf98XAp9clsg0rfp8pg4AjmvzpwJPTZJljHESVtt3gz5WW33ax2qrg/rqc10K2KbNbwv8+zLGt+yq6mzgxiGrHAAcX51zgO2SPGwUx15NjbqdgasHlq9pZXM5jK4lPa16nW+SVyb5N+AvgNcsU2zjMO/5tlvcu1bV55YzsDHp+3l+bru9f2qSXWd5f1r0Od+fB34+yVeTnJNkv2WLbvR6/71qXcF2B768DHFpevX5TP1snaq6E7gJeMiyRDc5q+27QR+rrT7tY7XVQX31uS5HAL+X5Bq6rL6vXp7QNlkL/ZvT22pq1PWW5PeAtcBfTjqWcauqD1XVI4E3An866XjGJcn9gPcAb5h0LMvo74E1VfVY4Azu+QV+pdqcrvvLU+juXH00yXYTjWh5HAicWlV3TToQaSVbTd8Nhlml9Wkfq7UOms8LgWOrahe6rocntM+QRmw1XdRrgcE7Fbu0sntJ8l+BPwH2r6rblim2ceh1vgNOAp4z1ojGa77z3Rr4JeCsJFfS9WNeN8WDu+f9962qGwY+wx8DnrBMsY1Dn8/zNcC6qrqjqq4A/pWugp1GC/n/eyB2vdT8+nymfrZOks3pukrdsCzRTc5q+27Qx2qrT/tYbXVQX32uy2HAKQBV9TVgS2DHZYlu07TQ7+e9raZG3deBPZLsnmQLui9C6wZXSPIrwN/Q/dGe9r7Qfc538I/Ns4DLljG+URt6vlV1U1XtWFVrqmoN3biI/atq/WTCXbI+/76DfbT3By5dxvhGbd7zBf4v3S+kJNmRrivM5csZ5Aj1OV+SPAbYHvjaMsen6dPnM7UOOLjNPw/4crWR/SvYavtu0Mdqq0/7WG11UF99rstVwFMBkvwCXaPu+mWNctOyDnhxy4K5D3BTVV03ih1vPoqdTIOqujPJq4DT6bL1HFNVFyd5G7C+qtbRdal4MPC3bWz4VVW1/8SCXoKe5/uq9uvjHcAPuKcynzo9z3fF6Hm+r0myP3An3aDdQyYW8BL1PN/TgacluQS4C/jjqprKuwwL+DwfCJy0Cr54a4l6fqaOpusatYHub8aBk4t4eay27wZ9rLb6tI/VVgf11fO6vIGuK+of0iVNOWQl11lJPknXuN+xjSN8C3B/gKr6CN24wmcCG4BbgUNHduwVfF0lSZIkacVbTd0vJUmSJGnFsVEnSZIkSVPMRp0kSZIkTTEbdZIkSZI0xWzUSZIkSdIUs1GnFSvJ65I8cET7ekySC5NckOSRi9h+bZIjF7HdaUm2W+h245Dkn5ew7SFJHg7nFl0AAAYrSURBVD7KeCRJw1kP3ms/D09y6gK3OSTJB9v8HyR58TzrfyzJnsP2sxhJbl7stlo9fKSBpla6Bwalqu6e4/0rgbVV9R8jONbhwOZV9edL3de0SbJ5Vd25xH2cBfzRCn84rSQtK+vB8UpyCN31e9Uk95Pk5qp68FJi0MrnnTpNlSRrknw7yfHAt4Bdk3w4yfokFyd5a1vvNcDDgTOTnNnKnpbka0m+keRvk9znD2SSxyc5J8k3k3w6yfZJngm8Dnj5zL422ubmJH/Zjv8PSZ6U5Kwkl7eHf5PkKUk+2+Z/s/3aOfOL59ZJHpbk7Fb2rSS/3ta9MsmO7bwvTfLRdpwvJtmqrfPEFu+FLY5vzRLjSUmeNbB8bJLntf3+Y7sm30jyXwbi/cck64BLZs6zvT44yZfa+hclOWDg3+Y+MSZ5HrAW+ESLcask70xySYv7rxb1YZCkVch6cNH14JqZ8nR3zj6V5AtJLkvyFwPrHZrkX5OcBzx5oPyIJH+U7o7leRvt96I2f1aStfPs59hWL/7s2rXXWevWjc5h1mskAVBVTk5TMwFrgLuBfQbKdmivmwFnAY9ty1cCO7b5HYGzgQe15TcCb55l/98EfrPNvw14X5s/gu5O02wxFfCMNv9p4IvA/YHHARe28qcAn23zfw88uc0/GNgceAPwJwPnsfXgObTzvhN4fCs/Bfi9Nv8t4Ffb/DuBb80S438DjmvzWwBXA1sBDwS2bOV7AOsH4r0F2H1gHze3182BbQau6wYg88R4Ft2vlAAPAb7NPT0Ftpv058rJyclpWibrwUXXg2tmyoFDgMuBbYEtge8AuwIPA64Cdmp15VeBD258/sCFM/Vju45/2ubPovsRc9h+jgWeNxDX0Lp1o3VmvUZOTlXlnTpNpe9U1TkDyy9I8g3gAuAXgfv0Zwf2aeVfTXIhcDDwc4MrJNmWroHxlVZ0HPAbPeK5HfhCm78I+EpV3dHm18yy/leB97RfUberrmvj14FDkxwB/HJV/XiW7a6oqgvb/PnAmnTjDLauqq+18hPniPHzwL5JHgA8Azi7qn5CV+l+tP3K+Lfc+9qdV1VXzLKvAP87yTeBfwB2Bh46V4yzbH8T8FPg6CT/Hbh1jpglSbOzHlx4PbixL1XVTVX1U7oeKT8H7A2cVVXXV9XtwMlzbHsK8Dtt/ndmWa/vfgYNq1tn9LlGWqVs1Gka3TIzk2R34I+Ap1bVY4HP0f3qtrEAZ1TV49u0Z1UdNqJ47qiqmcGpdwO3AVQ3xmHzjVeuqncCL6G7U/bVJI+pqrPpKs5rgWMz+2Ds2wbm75pt33NpldZZwNO5dwX0h8D36H5NXUv3i+KMW5jd79L9+viEqnp8237mms8bY6u8nwScCjybe74ISJL6sR5cYD044n2dTNeQ/nmgquqyBWx7J+37d5L7cU+9O6xuhe5Afa6RVikbdZp229BVbjcleSjdXagZPwa2bvPnAE9O8iiAJA9qf4x/pqpuAn4w0Ef9RcBXGLEkj6yqi6rqXXS/uj0myc8B36uqjwIfA/bqs6+q+iHw4yR7t6IDh6x+MnAo8Ovc05DaFriuVbwvouvOMZ9tge9X1R1J9mWjX3rn8LN/izaGY9uqOo2uUfm4HttLkmZnPdi/HpzPucBvJnlIkvsDz5/jmP9G1xD8M2a/CzdsP1cCT2jz+9P1mIEedetir5FWh6X8wiFNXFX9S5ILgP9HN07sqwNvHwV8Icm/V9W+6bJPfbJ1QQT4U+BfN9rlwcBH0qWAvpyuETRqr2t/sO8GLqbrGnkg8MdJ7gBuBhby69thdF0o76arfG+aY70vAicAn2ndQQD+Gvi79mvfF5j77tygTwB/37psrqe79vM5lu66/oTuC8dnkmxJ98vx63tsL0mahfUg0L8eHKqqrmtdG78G/JBu7NxcTgb+Eth9gfv5KF0d+C/cu97tU7c+hcVfI61wPtJAmnJJHlxVM9mzDgceVlWvnXBYkiQtC+tByTt10krwrCRvovv//B26rF6SJK0W1oNa9bxTJ0mSJElTzEQpkiRJkjTFbNRJkiRJ0hSzUSdJkiRJU8xGnSRJkiRNMRt1kiRJkjTFbNRJkiRJ0hT7//ll6NBKxca9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "o,e = run_plink(\"--missing\", file=path_intermed, out=path_missing, \n",
    "                extension='vmiss', log_name='missing', force=True)\n",
    "df_miss_ind = pd.read_csv(path_missing+'.smiss', sep='\\s+', \n",
    "                          usecols=['F_MISS'], dtype=dtype)\n",
    "df_miss_var = pd.read_csv(path_missing+'.vmiss', sep='\\s+', \n",
    "                          usecols=['#CHROM', 'ID', 'F_MISS'], dtype=dtype)\n",
    "plot_missing(df_miss_ind, df_miss_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Report\n",
    "\n",
    "If we carefully look at the plot on the left, we reassuringly observe that the histogram was **heavily shifted to the left** - i.e. to lower missing rates. This is expected since the raw data had many fully missing variants, which increased the missing rate of variants per individual. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Median missingness per chromosome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZQU5bnH8e8T0KjgzqCERVDRJC5BQdLeqJi4sGhEhRhIMCoiV+/lxuiNW7yHm+WYuMUsJ14NIq6JBEENyMCgiaJZWhh0QDYVQQWCLAYlYKKiz/2jaoam1+qZqekZ6vc5h0NVve9b9XR1Tz391vK2uTsiIpJcn6p0ACIiUllKBCIiCadEICKScEoEIiIJp0QgIpJw7SsdQGN06tTJe/bsWekwRETalAULFmxy96rs5W0yEfTs2ZPa2tpKhyEi0qaY2Zv5luvUkIhIwikRiIgknBKBiEjCKRGIiCScEoGISMIpEYiIJJwSgYhIwikRtJC6u75a6RAKuv+BMysdQrMZ8sT3Kh1CIox4bFWlQ2hWL0/YUHabdbeuiyGSylAiaAH1SaA1JoP6JLArJIP6JKBkEK/6JLCrJIP6JFBOMqhPArtKMog1EZjZJDPbYGaLC5Sbmf3SzFaY2SIzOz7OeCoh++DfmpJB9sG/LSeD7IO/kkF0KzZvjVw3++Df1pNB9sE/SjLIPvjvCskg7h7B/cCgIuWDgd7hv7HAXTHHIyIZqles57o/LmXW6+srHYpUUKyJwN2fA/5epMpQ4EEPpIH9zKxLnDG1pELf/ltDr6DQt/+22Cso9O1fvYLiPtz+CfcvCoaeuW/hm3y4/ZOi9Qt9+2+rvYJC3/6L9QoKfftv672CSg861xVYnTG/JlyWs1fNbCxBr4EePXrkrOjd3z/VML3f0DMibXzz76c0TO8/9IJIbTY+dgsAVedfF6k+wF4HfRaA99cvj9zmpSevB+C4s2+O3Ob5WdcCcPLgWyPVP/igYxum316/KFKbqU9f0zA9/PTbIrWZ+OyONmNOjdbmquevbpj+2cl3RGpzTKdeALy8KfqB6brngs/ALadEe/8Brn/uSQBuPuXs6G3mPh20GXB6pPo3zH2+YfonA06O1ObGuTsGYrxpQL+S9e+Yt4KPw58s/9jhZ/NWcN2/HVG0Tae9cg8Zm97fXrTNxJqNOcvGDMwZADPH3LnvATBgwL4l66Yfz/99M3XeAQXb7NYx93vwR1uLJ8N2+7TLWfbxlo+Lttl4zys5y6ouO7Jom02T5ucs6zT6hKJtmqLSiSAyd58ATADo16+fZ5a98/DjZM8fOOq8ouvb9PDEnPlOo8YUbbPugWtyprtcVPjAtlvHTgB8tG3TzvNbNxVs89fJl+SdP3HEfQXb1Ey5KGd+4AUPFKwP0LHDQWzdun7n+W3FTw9M+v2FOfOjhz5UtM1tM7+VM3/NWQ8WbXPhnG/lzD90ZvE2nffcn/Xb3m2YBtjwz80F618ye2Le+fsGFf4MXDLrkbzz9w0eWbDN6OoncuYnDTm3YH2AS6trcubvHTKwaJvLZv0pZ/6ewScVrP/Gu9uYv+7dnZbNW/cub773Pofsu1fRbcmup9KJYC3QPWO+W7hsl9bnihmVDqGgiy+aU+kQylZ97o+5uCY3IVef++MKRNM2TF6a/89s8pI1BXsFk8/vxbjZq/Mub26zZm3OmR48eP9m304+x4ztXHabLtc2/Yz2+l/+KXL5Qd8unOQbo9KJYDowzswmA18E3nP3tn2yTaQNGHlUt5weAcCIo7pVIJqWVerOoPryxiSEtirWRGBmjwCnAp3MbA3wv8BuAO5+N1ANDAFWAO8Dl+Rfk7SEttgbkMY5ZN+9OKHLfjslg/5d9tNpoQJKXQyuL2+OnkElxJoI3L3wydOg3IH/jDMGEcnv6v6HM2p6LR87tDO4qv/hzbbuOx5/O1L51ecd3KTtPPtw7oXofOWnjip9cbo5vX1H3ken8pYffPXRcYdTUqVPDYlIheze/lNcfOwh3LvwTS75wiHs3j7/3eSlbg+tL4/jWoG0DCUCkTbq3KlPRS5/Ynj+W6qHHH4QRxzYgcP379issUnbokQgiXLWY8WfSagvn3n+1UXr7UqUBESJQJrdjY8WG1Vk5/KbvjYbgEHThxRtU18++5zqhmWlnhyuL2/qbaRnTbsnUvnMYZc1aTtfnfp4pLIZw4s/IyNSLiUCKeoXvy3+IFNm+ZXfqClSU0RaKw1DLSKScOoRtBLPTDwrUvmXx8xs0naKDSqXWaZnCkSSQ4kgBqVGF60vb+pQE49PGhyp7LzRs5q0HSnf2VN/E6nsyeHfbIlwRIpKVCLYcHfxO0bqyztfnpw7RkREEpUIGmPNr0ZHKu82blJLhCMi0ux0sVhEJOGUCEREEq7NnhraeNfDkcurrhgVdzgi0kTTphX+wabM8mHDOrVEOImiHoGISMIpEYiIJJwSgYhIwikRiIgknBKBiEjCKRGIiCScEoFIwi3ZuKXSIUiFKRGIiCScEoGISMIpEYiIJJwSgYhIwikRiIgknBKBiEjCKRGIiCScEoFIwm14/8NKhyAVpkQgIpJwSgQiIgkXeyIws0Fm9oqZrTCz6/OU9zCzZ8zsJTNbZGZD4o5JRER2iDURmFk74E5gMPB5YKSZfT6r2v8AU9z9OGAE8H9xxiQiIjuLu0fQH1jh7ivd/UNgMjA0q44D+4TT+wJ/izkmERHJEHci6AqszphfEy7L9H1glJmtAaqB/8q3IjMba2a1Zla7cePGOGIVEUmk1nCxeCRwv7t3A4YAD5lZTlzuPsHd+7l7v6qqqhYPUkRkVxV3IlgLdM+Y7xYuy3QpMAXA3f8K7AF0ijkuEREJxZ0I5gO9zayXme1OcDF4eladt4DTAMzscwSJQOd+RERaSKyJwN23A+OAGmAZwd1BS8zsh2Z2Tljtv4HLzGwh8Ahwsbt7nHGJiMgO7ePegLtXE1wEzlw2PmN6KfCluOMQEWmt0uk0K1+bF7n+oen2pFKpZtt+a7hYLCIVtGHbB5UOQSos9h6BiIgUl0ql6DVve+T6BzVjbwDUIxARSTwlAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQSbgN739Q6RCkwpQIREQSTolARCThlAhERBJOiUBEJOGUCEREEi72RGBmg8zsFTNbYWbXF6hzgZktNbMlZvbbuGMSEZEd2se5cjNrB9wJnAGsAeab2XR3X5pRpzdwA/Ald99sZp3jjElERHYWd4+gP7DC3Ve6+4fAZGBoVp3LgDvdfTOAu2+IOSYREckQdyLoCqzOmF8TLst0BHCEmf3ZzNJmNijfisxsrJnVmlntxo0bYwpXRCR5WsPF4vZAb+BUYCRwj5ntl13J3Se4ez9371dVVdXCIYqI7LrKTgRmtpuZHRfxXP5aoHvGfLdwWaY1wHR3/8jdVwGvEiQGEYnZ+LnL2LjtA8bPXVbpUKSCSiYCM7vbzI4Kp/cFFgIPAi+Z2cgSzecDvc2sl5ntDowApmfVeYKgN4CZdSI4VbSynBchIiKNF6VHcLK7LwmnLwFedfdjgL7AtcUauvt2YBxQAywDprj7EjP7oZmdE1arAd4xs6XAM8A17v5OI16LiIg0QpTbRz/MmD4DeBTA3d82s5KN3b0aqM5aNj5j2oGrw38iItLCovQI3jWzs83sOOBLwGwAM2sP7BlncCIiEr8oPYJ/B34JHAx8x93fDpefBsyMKzAREWkZJROBu78K5Nzb7+41BOf3RUSkDYty19Bl4TAQWOA+M9tiZovC00UiItKGRblGcCXwRjg9EjgW6EVwcfeX8YQlInG7fFYdG97/oOHf5bPqKh2SVEiURLDd3T8Kp88GHnT3d9z9aaBDfKGJiEhLiJIIPjGzLma2B8EF4qczynTXkIhIGxflrqHxQC3QjmAoiCUAZjYAPQEsItLmRblr6EkzOwTYu36o6FAt8PXYIhMRkRZRMhGY2fkZ0/mqPNacAYmISMuKcmroq0XKHCUCEZE2LcqpoUtaIhAREamMKKeGig4G5+53NF84IiLS0qKcGrodqANmAR8ApYccFRGRNiNKIjiO4Inis4AFwCPAH8Lho0VEpI0r+UCZuy909+vdvQ9wLzAUWJrxwzIiItKGRekRAGBmVQS9g2MIfmd4Q1xBiUg8hk2bF7l82rD+cYcjrUSUi8WjgQuAPYCpwAXuriQgIrKLiNIjmAgsBt4EBgJnZj5Y5u46RSQi0oZFSQRfjj0KERGpmCgPlM2NsiIzm+buw5oekoiItKQow1BHdWgzrktERFpIcyYCPVcgItIGNWciEBGRNqg5E4GGnhARaYOiPEfQw93firCu65ohHhGRNi2dTrNyxcLI9Q9NbyWVSsUYUWlRegRP1E+Y2bRCldx9TrNEJCIiLSrKcwSZp3x0Z5CIJEo6neb1la9HqntY+jBSqRQ9/9Ix8voPTh3d2NCaTZQegReYFhGRXUCUHsEXzGwLQc9gz3CacN7dfZ9ijc1sEPALoB0w0d1vLlBvGMFYRie4e23UFyAiElU6naZu+apIdbele5FKpUilUhzy3CGR2nRJdWlKeBUT5cnido1duZm1A+4EziAYsXS+mU1396VZ9fYGrgReaOy2RKTtSqfTLFr0RqS6Xbv2JJVKkU6nmb8kWps90j0rfkG2NYs8DHUj9QdWuPtKADObTPh7Bln1fgTcAlwTczwiErN0Os2KhSuj1e1yaIseoFOpFB0WRbvUeUyqc8zRtB5xJ4KuwOqM+TXAFzMrmNnxQHd3n2lmBROBmY0FxgL06NEjhlBFpFJSqRRr1x4esW6nhjb/WnFYxDZVjY4tCeJOBEWZ2aeAO4CLS9V19wnABIB+/frporVIK5VKpfjLup4R6x4cbzASSdxDTKwFumfMdwuX1dsbOBp41szeAFLAdDPrF3NcIiISirtHMB/obWa9CBLACOAb9YXu/h7QqX7ezJ4Fvqu7hkRah3Q6zeZFy6PV/cx6XZBto2LtEbj7dmAcUAMsA6a4+xIz+6GZ6ZfNRERagdivEbh7NVCdtWx8gbqnxh2PiESXSqXY/28HRazbK+ZoJC4ahlpEJOEqeteQiATS6TT/WBRtxMp0t+Dp1S2LFkRff7e9AXhv0bzobbp+onP+CaFEINLMgoP6kmh1u6V1sJWKUyIQaQVSqRR7r1kXuS7APmv+Udb6AfZdG/1scCrVP3JdaduUCESKCL7dR7x9smu6YZCyvddEG7ZYvQFpDXSxWEQk4dQjkMQIvt2viFb3Mxnf7te+HKmNvt1LW6VEIBWXTqfZWrc1Wt3O6YYhiLcuXF26AZA+WBdkRYpRIpDESKVS7P23v0SuK5IUSgRScalUio4bov3Ga/0BOpVK0fHt6WW1EZH8lAikoHQ6zWuLt5VVH2D1y9FO8wCku+u0jUil6a4hEZGEU49ACkqlUrywskNZ9QFmrI52miezjYhUjhJBgqTTaZYujXaqJ53WKRuRpNCpIRGRhFOPoI1Kp9MsXBbtomyX9I6Ho5a/Eu1Uj3oDIsmhHoGISMKpR9AKpNNpapdHG0lyz4xv9+uWlnfvvYhIPkoERaTTaV57dWOkur3TO4Y+ePGVaAf1f+mCrIi0AolJBOl0mlWvRBtwrFcLH6BTqRT/XLx35LoiIs0pMYmgMVKpFN1qqyLV7ZYx9MEeL0U7qPfRQV1EWoHEJIJUKsWhddEGHOusA7SIJIjuGhIRSbg22SPYtm0bf3gl2o+DA/TURVkRkYLUIxARSbg22SPo0KEDpx15VOT6VeoNiIgUpB6BiEjCKRGIiCScEoGISMIpEYiIJFzsicDMBpnZK2a2wsyuz1N+tZktNbNFZvYHMzsk7phERGSHWBOBmbUD7gQGA58HRprZ57OqvQT0c/djganArXHGJCIiO4u7R9AfWOHuK939Q2AyMDSzgrs/4+7vh7NpoFvMMYmISIa4E0FXYHXG/JpwWSGXArPyFZjZWDOrNbPajRujDQ0tIiKltZqLxWY2CugH3Jav3N0nuHs/d+9XVRVtRFARESkt7ieL1wLdM+a7hct2YmanAzcCA9z9g5hjEhGRDHH3COYDvc2sl5ntDowApmdWMLPjgF8D57j7hpjjERGRLLEmAnffDowDaoBlwBR3X2JmPzSzc8JqtwEdgUfNrM7MphdYnYiIxCD2QefcvRqozlo2PmP69LhjEBGRwlrNxWIREakMJQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSbjYf49ARFqPacP6N0xfPqsup/zuwX1aMhxpJdQjEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOD1HICLSChz07ZMapjdNmp9T3mn0CbFtWz0CEZGEU49ARFqFYcM6NUzPmrU5p3zw4P1bMpxEUY9ARCTh1CMQEWlmB199dMP0xnteySmvuuzIlgynJCUCkRg8OfybAIyufiKnbNKQc/O2mTH8PAAura7JKbt3yMBmjC5+V593cMP0xJqNOeVjBlY1y3ZOHRWsJ/343/OWp847IGfZMWM7A7D8t5tyyj77jU45ywC6XNsFgA13b8gp63x552jBtmKJSgSdL7+6YXrTwxNzyjuNGtMs2+lzxQwAljx0Sd7yoy68L2fZl8fMBOCvk/O3OXFEbpvzRs8CoGbKRTllAy94IO96Lr5oDgBTp16YUzZ8+EM5y678xo6D0qTf57YZPTS3zU1fm90wfdvMb+WUX3PWgznLZp9T3TB94ZzcNg+dmdum+twfN0xfXHNbTvn9A6/JWTbz/OAzcMns3Pcf4L5BuZ+BmcMuC9rMeiR/m8Ej8y6P2xPDz2iYvmzWn3LK7xl8Us6yxph8fq+G6XGzV+eU/2pQ92bZjlRO7InAzAYBvwDaARPd/eas8k8DDwJ9gXeAr7v7G6XWW3XFqIbpdx5+PKf8wFHnNSXsBt3GTWqYXvdA7oGly0W5ByARkbYk1kRgZu2AO4EzgDXAfDOb7u5LM6pdCmx298PNbARwC/D1OOOSxhk99CGmPr0jGQ4/vfUkwfsHXsN1f7oHgFtOuqzC0exQfxro+rlPA3DzgNNLtrl3yEBumPt8w/xPBpxcss09g0/ixrm1DfM3DehXsk1jfnvgV4O684Pn1jXM/+8pXUq2GTOwiinPB6duLjg591RNPoMH78/cue8BMGDAviXr158CqnvqXQD6nLFfyTb1p4FWzgjuUDr0q6XvSup8eWfeeeSdhvkDRx5Ysk3VZUfy90dWNswfMPLQkm06jT6Bv09ZuKPNBV8o2aZJ3D22f8CJQE3G/A3ADVl1aoATw+n2wCbAiq23b9++nm3zE3Ma/pXjvT/MKqv+hmk3+4ZpN5fVZlXNzb6qprw2L864zl+ccV1ZbWqfu72s+u7uzz9/W9ltav5SfptH55Xf5qb5N5Xd5qcLHi2vfu0s/2lteZ+BO+Y/63fMf7a8NvP+UlZ9d/efzastu83P5y8uu01j3Dl/Q9ltZtW+G0MkO1v25y2+7M9bymqz+pn3yt7O5pmby28z462y27w7c1nZbYoBaj3PMdWCsniY2XBgkLuPCecvBL7o7uMy6iwO66wJ518P62zKWtdYYCxAjx49+r755puxxS0isisyswXuntNdbDPPEbj7BHfv5+79qqqa544DERGJPxGsBTJvKegWLstbx8zaA/sSXDQWEZEWEHcimA/0NrNeZrY7MAKYnlVnOlB//+Nw4I8e5/kqERHZSax3Dbn7djMbR3BBuB0wyd2XmNkPCS5aTAfuBR4ysxXA3wmShYiItJDYnyNw92qgOmvZ+IzpfwFfizsOERHJr81cLBYRkXgoEYiIJJwSgYhIwikRiIgkXKxPFsfFzDYC+R4t7kQwREU51KZl2rTWuNSm9calNs2/jUPcPfeJ3HzjTrTVfxQYR0NtKt+mtcalNq03LrVpubh0akhEJOGUCEREEm5XSwQT1KbVtmmtcalN641LbVoorjZ5sVhERJrPrtYjEBGRMikRiIgkXJtNBGY2ycw2hL9wVr/sd2ZWF/57w8zqIrT5mpktMbNPzCznl3sKtPmCmf3VzF42sxlmtk+ROLub2TNmtjTczpWNeW0l6u9hZvPMbGG4jR9EXa+Z3WZmy81skZk9bmb7lWqTUfbfZuZm1qlEfG+E+6rOzGoL1MkX2/fNbG3Gezokq03efWtmB5jZU2b2Wvh/wR+jNbP9zGxquA+WmdmJUWILl/9X2G6Jmd0aIa6i+zqj/ZEZr7nOzLaY2XcKvYaMdlea2eJwmzn1w/fq4Yz59ma20cyeNLNLMrb3Ycb7dXOB9ybyPs5o087MXjKzJ0vU625mq8zsgHB+/3C+Z566+WL7UbiP68xsjpl9JmvdZb83BbbTx8zS9Z9rM+tf5DVdFW5vsZk9YmZ7FHjd+WIreowK65iZ/cnMBmcs+5qZzS4U007Kvd+0tfwDTgGOBxYXKP8pML5UG+BzwJHAs0C/KNsh+J2FAeH0aOBHReLsAhwfTu8NvAp8vimvLU99AzqG07sBLwCpiK/lTKB9OH0LcEuUWAh+TKiG4MG+TiXieyNCnXyxfR/4brn7FrgVuD5cfn32a8paxwPAmHB6d2C/iLF9GXga+HQ43zlCXEX3dYH42gFvEzwIVKze0cBiYC+CUYWfBg7PqrMVqAP2DOcHh/NPFnu/Crz+yPs4o83VwG+zt1eg7rXAhHD612T91nmJ2PbJmP42cHdT35sC25kDDA6nhwDPFoixK7AqY79PAS4u4/Nc9BiV9RlYBuwBdAReAw4rta/d2/BzBO7+HMHvF+QwMwMuAB4p1cbdl7n7K2Vu5wjguXD6KWBYkfbr3P3FcPofBG9U10L1i2yzWH13963h7G7hv5y7AAq8/jnuvj2cTRP8ilyUWH5G8MfaLHcblPuawzaF9u1QggM84f/n5mtvZvsS/IHfG67jQ3d/N2JsVwA3u/sHYZ0NpeIqta8LOA143d1L/Uj354AX3P39cBtzgfPz1KsGzgqnR5L1N5JPgdcfaR/XM7Nu4XYnltpe6GdAKuzZnATcHjU2d9+SMduBjM9oY9+bAvvAgfqzAfsCfyvyetoDe1rwK4x75atbJLaix6iM9ouBGcB1wHjgQXd/vVQ7aMOnhko4GVjv7q/FtP4lBH8IEPyWQvcidRuEXdvjCL6xN6uw210HbACecvfGbGM0MCvCtoYCa919YcT1OjDHzBaY2dgyYxoXdtcnlTjF05Md+/Ygd18XFr0NHFSgWS9gI3BfeMpiopl1iBjXEcDJZvaCmc01sxMixJUp0r4m+KGmkgdrgt7AyWZ2oJntRfANNd/ncjIwIjw1cWyeuKKKuo/r/Zzgi8MnUVbu7h8B1xAkhO+E85GZ2U1mthr4JsFBMV+dnjTtvfkOcFu4nduBG/JVcve1YflbwDrgPXefUyL+QrGV8gPgGwS9vVtL1G2wqyaCSN90mmA08B9mtoCgC/dhqQZm1hGYRvCh3lKqfrnc/WN370PwTaa/mR1dTnszuxHYDvymRL29gO9R4I+rgJPc/XiCD+d/mtkpEdvdBRwG9CH4A/ppgZgK7lsP+syFei3tCbr7d7n7ccA2gtMcUbQHDgBSBAesKWFPtGRcZezr3YFzgEdLBePuywhOacwBZhOc8vk4T71FQE+Cv5Hq7PLGKLGPMbOzgQ3uvqDMVQ8meN/L+iyHMd3o7t0J9vG4PDE16b0JXQFcFW7nKsKeZZ5t7U/wxbEX8Bmgg5mNKrTSphwr3H0b8DvgofreahS7XCIIu17nE+yMWLj7cnc/0937EiScot0vM9uN4I39jbs/FldcYWzvAs8Ag6K2MbOLgbOBb4Z/1MUcRvCBXmhmbxAknhfN7OAiMa0N/98APA4UvKiW1W59mOA+Ae7J167Avl1vZl3C8i4EvaR81gBrMnpPUwkSQxRrgMfC03LzCL7pNlw0L/Sel7mvBwMvuvv6KAG5+73u3tfdTwE2E5xjzmc6wTfUpnxZirqPAb4EnBN+XiYDX7GMi9b5mFkf4AyCRHtV/bYa4TdknbptpvcGgt9ar2//KIU/16cDq9x9Y9izeQz4t3wVm+lY8QkRe171drlEQLDTl7v7mrg2YGadw/8/BfwPcHeRukbwTWGZu98RUzxV9Xc5mNmeBH9AyyO2HUTQZT/H3d8vVd/dX3b3zu7e0917EhwQj3f3twusv4OZ7V0/TXBRLurdUJl//Odltyuyb6cT/JES/v/7Aq/lbWC1mR0ZLjoNWBolNuAJggvGmNkRBBeaNxWLq9x9TZk924zPZQ+CL0O/LVB1EvADd3856rrziLSPAdz9BnfvFn5eRgB/dPdi34iNoDf4HXd/C7iNAtcICrTvnTE7lIy/hWZ8byA4zz8gnP4KwcXZfN4iuN6xV7j90wjO/2fHHfuxoiCPcEW5Nf4j+ANZB3xEcDC6NFx+P3B51DYEB5g1wAfAeqAmQpsrCb5tvQrcTPiEdoFtnkTQbV5E0F2vA4Y05rUVqX8s8FK4jcVk3S1V4rWsAFZnxHZ3ObFQ4o4g4FBgYfhvCXBjGbE9BLwcvq7pQJco+xY4EPgDwR/m08ABReLrA9SG63gC2D9ibLsDD4f7+0XgKxHiKrqvs7bZAXgH2LeMv4nnCRLZQuC0POVb8yw7ldJ3DeV7/ZH3cant5akzFvhdxny7cB8PiPjeTAvfl0UEF0+7NvW9KbCdk4AF4f5+Aehb5DX9gCAhLQ4/15+OeqygxDEqz3q+T5G77fL90xATIiIJtyueGhIRkTIoEYiIJJwSgYhIwikRiIgknBKBiEjCKRFI4pnZT8zsy2Z2rpndkLH8uxaMSFlnZvPN7Fvh8mczR4E0s54WjkppZqea2Xthm+VmdnvWts4Nh8xYZsEIn+dmlN1vZu/XP3cRLvu5RRjhVaQplAhE4IsEA40NIGlb8O4AAAIESURBVBxM0MwuJ3gwr78HQ3ecRjDKaxTPh22OA842sy+F6/wCwYNRQ939cwTDR9xuZsdmtF1BOI5V+MDiV4C1TXt5IsUpEUhiWTAG/SLgBOCvwBjgLjMbTzCe0hUejvXi7lvc/YHCa8vl7v8keCiofrTZ7wI/dvdVYfkq4CcEYxXVmwx8PZw+Ffgzwdg3IrFRIpDEcvdrCJ4QvZ8gGSxy92MJRsrc291XFmn+m/D0Tx0FBm8LBxvrzY4hy48ieBI1U224vN6rQFXYdiRBYhCJlRKBJN3xBEMEfJY8478U8U137xOeAhqSVXaymS0kOKVT4wXGYSriMYIxeb5IMGyESKzaVzoAkUoIR7e8n2D01E0EPxZi4Tf8E4GtZnZoiV5BIc+7+9lm1gtIm9kUd68jGAeoL0HiqdeXYAymTL8j6Dk84O6fZI1uLdLs1COQRHL3uvDbfP3PAf4RGBh+y/8nwbn7Oy38PWoz61h/11AZ21hFMCjhdeGi24Ebwh8dqf/xke+R9TsLHvwa2Y3A/zXmtYmUSz0CSSwzqwI2h9+6P+vumUNQ30Xwu6/zzewjglEn8/4wTgl3A981s57uXmdm1wEzwnHnPwKuDXsLO3H3XzdiWyKNotFHRUQSTqeGREQSTolARCThlAhERBJOiUBEJOGUCEREEk6JQEQk4ZQIREQS7v8BYoe53rEE3e8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_miss_chr = df_miss_var[['#CHROM', 'F_MISS']].groupby('#CHROM').median()\n",
    "sns.boxenplot(x='#CHROM', y='F_MISS', data=df_miss_var,\n",
    "             order=df_miss_chr.sort_values(by='F_MISS').index);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigate missingness patterns w.r.t. some clinical variables\n",
    "\n",
    "We want to make sure there are no patterns among variables such as sex, race, ... \n",
    "\n",
    "<br />\n",
    "\n",
    "**TO BE DONE**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final missingness processing\n",
    "\n",
    "Note: we need to separately filter missingness by individual and by variant, since plink first removes individuals (and we don't want that). Thus we first create a dataset with `--geno`, then with `--mind`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Managing paths and files\n",
    "file_prefix = file_prefix + '_clean'\n",
    "path_clean = path_out+file_prefix\n",
    "path_clean_variants = path_clean+'_var'\n",
    "path_clean_ind = path_clean+'_ind'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "run_plink: command '--geno 0.01 --make-bed', the file 'data/plink/host_geno_clean_var.bed' already exists (force is set to False).\n",
      "run_plink: command '--mind 0.1 --make-bed', the file 'data/plink/host_geno_clean_ind.bed' already exists (force is set to False).\n",
      "run_plink: command '--missing', the file 'data/plink/host_geno_clean_ind.vmiss' already exists (force is set to False).\n"
     ]
    }
   ],
   "source": [
    "run_plink('--geno ' + str(THRESHOLD_MISSING_VARIANTS) + ' --make-bed', file=path_intermed, \n",
    "          out=path_clean_variants, extension='bed', log_name='filter_variants');\n",
    "run_plink('--mind ' + str(THRESHOLD_MISSING_INDIVIDUALS) + ' --make-bed', \n",
    "          file=path_clean_variants, \n",
    "          out=path_clean_ind, extension='bed', log_name='filter_individuals');\n",
    "run_plink('--missing', file=path_clean_ind, out=path_clean_ind,\n",
    "         extension='vmiss', log_name='missing_final')\n",
    "df_miss_ind = pd.read_csv(path_clean_ind+'.smiss', sep='\\s+', \n",
    "                          usecols=['F_MISS'], dtype=dtype)\n",
    "df_miss_var = pd.read_csv(path_clean_ind+'.vmiss', sep='\\s+', \n",
    "                          usecols=['#CHROM', 'ID', 'F_MISS'], dtype=dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA38AAADQCAYAAAC3FdgpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3defxtdV3v8ddbZpVRiIsMHUxuhqaIR6C0G8MNAU30qgiWApFUgkNaCdVNnLpaKYoahYEMVwRCuZChSAhSJsMBjowSRwbhdBRkkiFB4HP/WN8f7PPjN+zfsH/j6/l4rMde67umz15nn993f/f6rs83VYUkSZIkaWF7xmwHIEmSJEkaPBt/kiRJkrQI2PiTJEmSpEXAxp8kSZIkLQI2/iRJkiRpEbDxJ0mSJEmLwJqzHcBUbLrpprVkyZLZDkOSNAOuuOKKH1fVZrMdx3xhHSlJi8NE6sd53fhbsmQJy5Ytm+0wJEkzIMltsx3DfGIdKUmLw0TqR7t9SpIkSdIiYONPkiRJkhYBG3+SJEmStAjY+JMkSZKkRcDGnyRJkiQtAvM62+d0OPXSH4y67i07bzODkUiSNLdYR0rSwuKdP0mSJElaBGz8SZIkSdIiYONPkiRJkhYBG3+SJEmStAjY+JMkSZKkRcDGnyRJfUiybpLLknw3yXVJPtjKT0xyS5LlbdqhlSfJMUlWJLk6yY49xzowyU1tOrCn/GVJrmn7HJMkrXyTJOe37c9PsvFMv39J0vxn40+SpP48AuxeVS8BdgD2SrJLW/fHVbVDm5a3sr2B7dp0KHAsdA054APAzsBOwAd6GnPHAm/v2W+vVn4EcEFVbQdc0JYlSZoQG3+SJPWhOg+2xbXaVGPssi9wctvvEmCjJFsArwLOr6p7qupe4Hy6huQWwAZVdUlVFXAy8LqeY53U5k/qKZckqW82/iRJ6lOSNZIsB+6ka8Bd2lZ9tHXtPDrJOq1sS+D2nt3vaGVjld8xQjnA5lW1qs3/ENh8lPgOTbIsybK77rprcm9SkrRg2fiTJKlPVfV4Ve0AbAXslORFwJHAC4CXA5sA7x9wDMUodxyr6riqWlpVSzfbbLNBhiFJmods/EmSNEFVdR9wIbBXVa1qXTsfAb5A9xwfwEpg657dtmplY5VvNUI5wI9at1Da653T+44kSYuBjT9JkvqQZLMkG7X59YDfAL7X0ygL3bN417ZdzgHe1rJ+7gLc37pungfsmWTjluhlT+C8tu4nSXZpx3obcHbPsYaygh7YUy5JUt/WnO0AJEmaJ7YATkqyBt2Pp2dU1VeTfDPJZkCA5cDvt+3PBfYBVgAPAwcDVNU9ST4MXN62+1BV3dPm3wGcCKwHfK1NAB8DzkhyCHAbsN/A3qUkacEaaOMvya3AA8DjwGNVtbSluD4dWALcCuxXVfe2Xzk/TVdRPgwcVFVXDjI+SZL6VVVXAy8doXz3UbYv4LBR1p0AnDBC+TLgRSOU3w3sMcGQJUlazUx0+9ytjXu0tC2PNlbRiOMhSZIkSZKmbjae+RttrKLRxkOSJEmSJE3RoBt/BXwjyRVJDm1lo41VNNq4R6txDCNJkiRJmrhBJ3x5ZVWtTPJzwPlJvte7sqoqyYhjFY2mqo4DjgNYunTphPaVJEmSpMVqoHf+qmple70TOItu7KPRxioabdwjSZIkSdIUDazxl+RZSdYfmqcbx+haRh+raLTxkCRJkiRJUzRut88kvwDcUVWPJNkVeDFdYpb7xtl1c+CsbgQH1gROraqvJ7mckccqGnE8JEmSJEnS1PXzzN+XgaVJnk/3rN3ZwKl0DbVRVdXNwEtGKB9xrKKxxkOSJEmSJE1NP90+n6iqx4DXA5+pqj8GHIJBkiRJkuaRfhp/P0tyAN3zeV9tZWsNLiRJkiRJ0nTrp/F3MPArwEer6pYk2wKnDDYsSZIkSdJ0GrfxV1XXV9W7qupLbfmWqvr44EOTJGnuSLJuksuSfDfJdUk+2Mq3TXJpkhVJTk+yditfpy2vaOuX9BzryFZ+Y5JX9ZTv1cpWJDmip3zEc0iSNBGjNv6SXJPk6tGmmQxSkqQ54BFg96p6CbADsFcbmujjwNFV9XzgXuCQtv0hwL2t/Oi2HUm2B/YHXgjsBfxtkjWSrAF8Dtgb2B44oG3LGOeQJKlvY935ew3wm2NMkiTNO0nenWSDNq7s8UmuTLLnePtV58G2uFabCtgdOLOVnwS8rs3v25Zp6/dIN/7RvsBpVfVIVd1CN8TRTm1aUVU3V9WjwGnAvm2f0c4hSVLfRh3qoapum8lAJEmaIb9TVZ9u3S03Bt5K9yz7N8bbsd2duwJ4Pt1duu8D97Ws2AB3AFu2+S2B2wGq6rEk9wPPaeWX9By2d5/bh5Xv3PYZ7RySJPVt3Gf+kuyS5PIkDyZ5NMnjSX4yE8FJkjQAaa/7AKdU1XU9ZWOqqseragdgK7o7dS8YTIiTk+TQJMuSLLvrrrtmOxxJ0hzTT7bPzwIHADcB6wG/S/drpyRJ89EVSb5B1/g7L8n6wBMTOUBV3QdcSJcNe6MkQz1ptgJWtvmVwNYAbf2GwN295cP2Ga387jHOMTyu46pqaVUt3WyzzSbyliRJi0A/jT+qagWwRvvF8wt0D6hLkjQfHQIcAby8qh4G1qYb1mhMSTZLslGbXw/4DeAGukbgG9tmBwJnt/lz2jJt/Terqlr5/i0b6LbAdsBlwOXAdi2z59p0SWHOafuMdg5Jkvo26jN/PR5uldDyJH8FrKLPRqMkSXPQ+VW1x9BCVd2d5AxgjzH2AdgCOKk99/cM4Iyq+mqS64HTknwEuAo4vm1/PHBKkhXAPXSNOarquna+64HHgMOq6nGAJIcD5wFrACe0LqkA7x/lHJIk9a2fxt9b6Sqhw4E/pOuS8oZBBiVJ0nRLsi7wTGDTJBvz1HN+G9BHApWquhp46QjlN9M9/ze8/KfAm0Y51keBj45Qfi5wbr/nkCRpIsZt/PVk/fwv4IODDUeSpIH5PeA9wHPpMnYONf5+Qvd8uyRJC9q4jb8kt9CNY7SaqnreQCKSJGkAqurTwKeTvLOqPjPb8UiSNNP66fa5tGd+XbouLJsMJhxJkgarqj6T5FeBJfTUg1V18qwFJUnSDOin2+fdw4o+leQK4C8GE5IkSYOT5BTgF4DlwOOtuAAbf5KkBa2fbp879iw+g+5OYD93DCVJmouWAtu3IRQkSVo0+mnEfaJn/jHgVmC/gUQjSdLgXQv8N7qhiyRJWjT66fa521RO0MZDWgasrKrXtAFtTwOeQ5dt7a1V9WiSdei63LwMuBt4c1XdOpVzS5I0gk2B65NcBjwyVFhVr529kCRJGrxRG39J3jvWjlX1yT7P8W7gBrpxlAA+DhxdVacl+TvgEODY9npvVT0/yf5tuzf3eQ5Jkvp11GwHIEnSbBjrzt/67fUXgZcD57Tl3wQu6+fgSbYCXk03kO17kwTYHXhL2+Qkukr4WGBfnqqQzwQ+myQ+kyFJmk5V9a3ZjkGSpNkwauOvqj4IkORiYMeqeqAtHwX8c5/H/xTwJzzVkHwOcF9VPdaW7wC2bPNbAre3cz+W5P62/Y/7fTOSJI0nyS7AZ4BfAtYG1gAeqqoNxtxRkqR57hl9bLM58GjP8qOtbExJXgPcWVVXTDK20Y57aJJlSZbddddd03loSdLi8FngAOAmYD3gd4HPzWpEkiTNgH6yfZ4MXJbkrLb8OuDEPvZ7BfDaJPvQDQ6/AfBpYKMka7a7f1sBK9v2K4GtgTuSrAlsSJf4ZTVVdRxwHMDSpUvtEipJmrCqWpFkjap6HPhCkquAI2c7LkmSBmncO39V9VHgYODeNh1cVf+nj/2OrKqtqmoJsD/wzar6LeBC4I1tswOBs9v8OW2Ztv6bPu8nSRqAh5OsDSxP8ldJ/pA+6sMkWye5MMn1Sa5L8u5WflSSlUmWt2mfnn2OTLIiyY1JXtVTvlcrW5HkiJ7ybZNc2spPb3GSZJ22vKKtXzJ9l0OStFiMWtkl2aC9bkI3tt8pbbqtlU3W++mSv6yge6bv+FZ+PPCcVv5e4IhR9pckaSreSvec3+HAQ3S9Tt7Qx36PAe+rqu2BXYDDkmzf1h1dVTu06VyAtm5/4IXAXsDfJlmjDYH0OWBvYHvggJ7jDGXEfj7dD66HtPInM2IDR7ftJEmakLG6fZ4KvIZuLL7eO3Bpy8/r9yRVdRFwUZu/GdhphG1+Cryp32NKkjQZVXVbm/0v4IMT2G8VbWD4qnogyQ08lbRsJPsCp1XVI8At7cfNofpvRasPSXIasG87nhmxJUkDM1a2z9e0121nLhxJkgYjyRlVtV+Sa1j9R00AqurFEzjWEuClwKV0z7gfnuRtwDK6u4P30jUML+nZrTfD9e3DynfGjNiSpAEbN+FLknOALwFnV9XDgw9JkqSBeHd7fc1UDpLk2cCXgfdU1U+SHAt8mK5B+WHgE8DvTOUcU4jtUOBQgG222WY2QpAkzWH9DPXwCeDXgBuSnJnkjUnWHXBckiRNq6pa1Z63O7Gqbhs+9XOMJGvRNfy+WFVfacf9UVU9XlVPAJ/nqa6dQ1mshwxluB6t/G5aRuxh5asda7yM2FW1tKqWbrbZZv28JUnSItJPts9vVdU76J7x+3tgP+DOQQcmSdJ0a0M7PJFkw4numyR0ycluqKpP9pRv0bPZ64Fr2/w5wP4tU+e2wHbAZcDlwHYts+fadElhzmnP75kRW5I0MP2M80eS9YDfBN4M7Ej3ELokSfPRg8A1Sc6ny/YJQFW9a5z9XkGXKfSaJMtb2Z/SZevcga7b563A77XjXZfkDOB6ukyhh7XGJ0kOB86jyzp6QlVd1473fuC0JB8BrmL1jNintKQx99A1GCVJmpB+nvk7g64Ly9eBzwLfal1bJEmaj77Spgmpqn+jy3g93Llj7PNR4KMjlJ870n5mxJYkDVI/d/6OBw4Y+rVSkqT5rKrsvSJJWpRGbfwl2b2qvgk8i278odXWDz3oLknSfJJkO+D/0A2w/mQCs6rqe/xaSZLmo7Hu/P068E26Z/2GKybRZUaSpDngC8AHgKOB3YCD6S/7tSRJ89pYg7x/oL0ePHPhSJI0cOtV1QVJ0oZ4OCrJFcBfzHZgkiQN0ljdPt871o69aa4lSZpHHknyDOCmlnVzJfDsWY5JkqSBG6vb5/rt9ReBl9ONMQRdN9DLBhmUJEkD9G7gmcC7gA/Tdf08cMw9JElaAMbq9vlBgCQXAztW1QNt+Sjgn2ckOkmSpt/jVfUg3Xh/PtogSVo0+nnAfXPg0Z7lR1uZJEnz0SeS3JDkw0leNNvBSJI0U/oZ5+9k4LIkZ7Xl1wGOkSRJmpeqarck/w3YD/j7JBsAp1fVR2Y5NEmSBmrcO39V9VHgd4B723RwVf3loAOTJGlQquqHVXUM8PvAcsz0KUlaBPq58wddxbhqaPsk21TVDwYWlSRJA5Lkl4A3A28A7gZOB943q0FJkjQDxm38JXkn3WC4PwIeB0I3yPuLBxuaJEkDcQJwGvCqqvrP2Q5GkqSZ0k/Cl3cDv1hVL6yqF1fVL1eVDT9J0rxUVb9SVZ+eaMMvydZJLkxyfZLrkry7lW+S5PwkN7XXjVt5khyTZEWSq5Ps2HOsA9v2NyU5sKf8ZUmuafsckyRjnUOSpInop/F3O3D/RA+cZN0klyX5bqskh4aO2DbJpa1iOz3J2q18nba8oq1fMtFzSpI0QI8B76uq7YFdgMOSbA8cAVxQVdsBF7RlgL2B7dp0KHAsdA05uh41OwM7AR/oacwdC7y9Z7+9Wvlo55AkqW/9NP5uBi5KcmSS9w5Nfez3CLB7Vb0E2AHYK8kuwMeBo6vq+XQJZA5p2x8C3NvKj27bSZI0J1TVqqq6ss0/ANwAbAnsy1NZsE+iy4pNKz+5OpcAGyXZAngVcH5V3VNV9wLn09WRWwAbVNUlVVV02bZ7jzXSOSRJ6ls/jb8f0FVMawPr90xjapXdg21xrTYVsDtwZisfXkkOVWxnAnsMdXeRJGmqkpzSXt89DcdaArwUuBTYvKpWtVU/5KmxcLek6z0z5I5WNlb5HSOUM8Y5JEnq27gJX6rqg5M9eJI1gCuA5wOfA74P3FdVj7VNeiu2JyvDqnosyf3Ac4AfT/b8kiT1eFmS5wK/k+RkugRmT6qqe/o5SJJnA18G3lNVP+n9nbKqKklNY8xPM9Y5khxK18WUbbbZZpBhSJLmoVEbf0k+VVXvSfJPdHfsVlNVrx3v4FX1OLBDko2As4AXTCXYFpcVmyRpMv6O7nm559H9MNnb+KtWPqYka9E1/L5YVV9pxT9KskVVrWpdN+9s5SuBrXt236qVrQR2HVZ+USvfaoTtxzrHaqrqOOA4gKVLlw60ESpJmn/G6vZ5Snv9G+ATI0x9q6r7gAuBX6F75mGo0dlbsT1ZSbb1G9KNvzT8WMdV1dKqWrrZZptNJAxJ0iJWVcdU1S8BJ1TV86pq256pn4ZfgOOBG6rqkz2rzgGGMnYeCJzdU/62lvVzF+D+1nXzPGDPJBu3RC97Aue1dT9Jsks719uGHWukc0iS1LdR7/xV1RXt9VuTOXCSzYCfVdV9SdYDfoMuicuFwBvpxlgaXkkeCHynrf9me+BdkqRpU1V/kOQlwK+1oour6uo+dn0F8FbgmiTLW9mfAh8DzkhyCHAbsF9bdy6wD7ACeBg4uJ3/niQfBi5v232op8vpO4ATgfWAr7WJMc4hSVLfxn3mbwq2AE5qz/09Azijqr6a5HrgtCQfAa6i+xWV9npKkhXAPcD+A4xNkrRIJXkX3eMDQ902v5jkuKr6zFj7VdW/Mew5wR57jLB9AYeNcqwT6AabH16+DHjRCOV3j3QOSZImYmCNv/Yr6ktHKL+Zblyj4eU/Bd40qHgkSWp+F9i5qh4CSPJxul4nYzb+JEma70Z95m86U2JLkjSHBHi8Z/lxRr+jJ0nSgjHWnb9pSYktSdIc8wXg0iRnteXX8dQjCJIkLVhjNf6mnBJbkqS5pqo+meQi4JWt6OCqumoWQ5IkaUaMle3zGOCYJMdW1R/MYEySJA1UVV0JXDnbcUiSNJPGTfgyhZTYkiRJkqQ5YqxB3oEnU2J/Efi5Nn0xyTsHHZgkSZIkafr0M9SDKbElSQtCG3v2X6pqt9mORZKkmTbunT9MiS1JWiCq6nHgiSQbznYskiTNtH7u/JkSW5K0kDwIXJPkfOChocKqetfshSRJ0uD1k/DFlNiSpIXkK22SJGlR6efOnymxJUkLRlWdlGQ9YJuqunG245Ekaab088yfJEkLRpLfBJYDX2/LOyQ5p4/9TkhyZ5Jre8qOSrIyyfI27dOz7sgkK5LcmORVPeV7tbIVSY7oKd82yaWt/PQka7fyddryirZ+yXRcB0nS4mPjT5K02BwF7ATcB1BVy4Hn9bHficBeI5QfXVU7tOlcgCTbA/sDL2z7/G2SNVq20c8BewPbAwe0bQE+3o71fOBe4JBWfghwbys/um0nSdKEjdn4axXVhTMVjCRJM+BnVXX/sLInxtupqi4G7unzHPsCp1XVI1V1C7CCrsG5E7Ciqm6uqkeB04B9kwTYHTiz7X8SXYK1oWOd1ObPBPZo20uSNCFjNv5MiS1JWoCuS/IWYI0k2yX5DPDvUzje4Umubt1CN25lWwK392xzRysbrfw5wH1V9diw8tWO1dbf37aXJGlC+un2OZQS+/gkxwxNgw5MkqQBeSddd8xHgC8BPwHeM8ljHQv8ArADsAr4xHQEOFlJDk2yLMmyu+66azZDkSTNQf1k+zQltiRpwaiqh4E/S/LxbrEemMKxfjQ0n+TzwFfb4kpg655Nt2pljFJ+N7BRkjXb3b3e7YeOdUeSNYEN2/YjxXMccBzA0qVLa7LvS5K0MPUzzp8psSVJC0aSlwMnAOu35fuB36mqKyZxrC2qalVbfD0wlAn0HODUJJ8EngtsB1wGBNguybZ0jbr9gbdUVbVn7N9I9xzggcDZPcc6EPhOW//NqrJhJ0masHEbfy0l9t8AawPbJtkB+FBVvXbQwUmSNADHA++oqn8FSPJK4AvAi8faKcmXgF2BTZPcAXwA2LXViwXcCvweQFVdl+QM4HrgMeCw9hw9SQ4HzgPWAE6oquvaKd4PnJbkI8BVLc6heE9JsoIu4cz+U70AglMv/cGo696y8zYzGIkkzZx+un0eRZed7CLoUmInGTcldpKtgZOBzekqxeOq6tNJNgFOB5bQVZT7VdW9LXPZp4F9gIeBg9rg8pIkTafHhxp+AFX1b0keG2uHtt0BIxQfP0LZ0PYfBT46Qvm5wLkjlN9MV98OL/8p8Kbx4pMkaTz9JHyZVEpsul8631dV2wO7AIe1sYyOAC6oqu2AC9oydGMebdemQ+keopckaVok2THJjsC3kvx9kl2T/HqSv6X9wClJ0kLWz52/1VJiA++ij5TY7RmIVW3+gSQ30KWr3peu2wx04xZdRNfVZV/g5PYcwyVJNhr2LIUkSVMxPBPnB3rmfYZOkrTg9dP4eyfwZzyVEvs84MMTOUmSJcBLgUuBzXsadD+k6xYKo499tFrjL8mhdHcG2WYb++RLkvpTVbvNdgySJM2mfrJ9TikldpJnA18G3lNVP+ke7Xvy2JVkQr+2msZakjQVSTYC3kb37PmT9WBVvWu2YpIkaSb0k+1z0imxk6xF1/D7YlUNjRX4o6HunEm2AO5s5WONiSRJ0nQ5F7gEuIb+nmGXJGlB6Kfb52RTYqfte0NVfbJn1dB4RR/j6eMYHZ7kNGBn4H6f95MkDcC6VfXe2Q5CkqSZ1k/jb1IpsYFXAG8FrkmyvJX9KV2j74wkhwC3Afu1defSDfOwgm6oh4P7ewuSJE3IKUneDnyV7nl2AKrqntkLSZKkwRu18dfSYUNLiU2X7KWAN9NHSuyq+jcgo6zeY4TtCzhsvONKkjRFjwJ/TZfMbOjZ8QLGHcNWkqT5bKw7f6bEliQtRO8Dnl9VP57tQCRJmkmjNv5MiS1JWqCGHi+QJGlR6SfbpymxJUkLyUPA8iQXsvozf9ZrkqQFrZ+EL6bEliQtJP+vTZIkLSr9NP5MiS1JWjCq6qTJ7JfkBOA1wJ1V9aJWtglwOl3vmFuB/arq3jbc0afpslg/DBxUVVe2fQ4E/rwd9iND8SR5GXAisB7dD6/vrqoa7RyTeQ+SpMXtGX1sc0qStyfZIskmQ9PAI5MkaQCS3JLk5uFTH7ueCOw1rOwI4IKq2g64oC0D7A1s16ZDgWPbuTehS6C2M7AT8IEkG7d9jgXe3rPfXuOcQ5KkCennzp8psSVJC8nSnvl1gTcB4/6oWVUXJ1kyrHhfYNc2fxLdUEjvb+Unt2GMLkmyUZIt2rbnD40pmOR8YK8kFwEbVNUlrfxk4HXA18Y4hyRJE9JP48+U2JKkBaOq7h5W9KkkVwB/MYnDbV5Vq9r8D4HN2/yWwO09293RysYqv2OE8rHOIUnShPTT+DMltiRpwUiyY8/iM+juBPZTH46pPZ830HFwxztHkkPpupmyzTbbDDIUSdI81E9lZ0psSdJC8ome+cdoSVQmeawfJdmiqla1bp13tvKVwNY9223VylbyVBfOofKLWvlWI2w/1jmepqqOA44DWLp06UAbopKk+aefxp8psSVJC0ZV7TaNhzsHOBD4WHs9u6f88CSn0SV3ub813s4D/rInycuewJFVdU+SnyTZBbiUbnzdz4xzDkmSJmTcxt9kU2JLkjQXJVkHeAPd0AlP1oNV9aFx9vsS3V27TZPcQZe182PAGUkOAW7jqTuI59IN8zD06MTB7Rz3JPkwcHnb7kNDyV+Ad/DUUA9faxNjnEOSpAkZt/GX5BaeyvL5pKoy26ckaT46G7gfuIKexxnGU1UHjLJqjxG2LeCwUY5zAnDCCOXLgBeNUH73SOeQJGmi+un2OamU2JIkzVFbVdXw8fokSVrwxh3kvaru7plWVtWngFfPQGySJA3Cvyf55dkOQpKkmdZPt8+BpMSWJGmWvBI4qD3W8AgQup6aL57dsCRJGqx+GnHTmRJbkqTZtvdsByBJ0mzoJ9vndKbEliRpVlXVbbMdgyRJs6Gfbp+TSoktSZIkSZo7+un2OamU2ElOAF4D3FlVL2plmwCn0zUkbwX2q6p7kwT4NN2YSA8DB1XVlf2/DUmSJA059dIfjLruLTtvM4ORSJpL+mn8TTYl9onAZ4GTe8qOAC6oqo8lOaItv5/u+Yvt2rQzcGx7lSRJkiRNg3GHemCSKbGr6mLgnmHF+wIntfmTgNf1lJ9cnUuAjZJsMdFzSpIkSZJG1s+dv+lMib15Va1q8z8ENm/zWwK392x3RytbxTBJDgUOBdhmG7stSJIkSVI/+mn8DSQldlVVkprEfscBxwEsXbp0wvtLkiRJ0mLUz1AP05kS+0dJtqiqVa1b552tfCWwdc92W7UySZIkSdI06OeZv+l0DnBgmz+QLpPoUPnb0tkFuL+ne6gkSXNakluTXJNkeZJlrWyTJOcnuam9btzKk+SYJCuSXJ1kx57jHNi2vynJgT3lL2vHX9H2zcy/S0nSfDewxl+SLwHfAX4xyR1JDgE+BvxGkpuA/9mWAc4FbgZWAJ8H3jGouCRJGpDdqmqHqlralocyXG8HXNCWYfUM14fSZbgeGg7pA3TZrncCPjDUYGzbvL1nv8lk4ZYkLXL9PPM3KVV1wCir9hhh2wIOG1QskiTNgn2BXdv8ScBFdMMbPZnhGrgkyVCG612B86vqHoAk5wN7JbkI2KBlwybJyXTZsr82Y+9EkrQgzHS3T0mSFqICvpHkipaVGiae4Xqs8jtGKJckaUIGdudPkqRF5JVVtTLJzwHnJ/le78rJZrieKIdDkiSNxTt/kiRNUVWtbK93AmfRPbP3o9adkz4zXI9VvtUI5SPFcVxVLa2qpZttttlU35YkaYGx8SdJ0hQkeVaS9YfmgT2Ba5l4huvzgD2TbNwSvewJnNfW/STJLi3L59t6jiVJUt/s9ilJ0tRsDpzVRl9YEzi1qr6e5HLgjJbt+jZgv7b9ucA+dBmuHwYOBqiqe5J8GLi8bfehoeQvdFmwTwTWo0v0YrIXaQSnXvqDUde9ZWe7Qks2/kM9WLYAAA/HSURBVCRJmoKquhl4yQjldzPBDNdVdQJwwgjly4AXTTlYSdKiZrdPSZIkSVoEbPxJkiRJ0iJg40+SJEmSFgEbf5IkSZK0CNj4kyRJkqRFwMafJEmSJC0CNv4kSZIkaRGw8SdJkiRJi4CNP0mSJElaBNac7QAkSZIkje3US38w6rq37LzNDEai+cw7f5IkSZK0CHjnbwz+wiJJkiRpoZhTd/6S7JXkxiQrkhwx2/FIkjRXWEdKkqZqzjT+kqwBfA7YG9geOCDJ9rMblSRJs886UpI0HeZSt8+dgBVVdTNAktOAfYHrZzWqUYzVJXQsdheVJE3CvKojJWmiRvtu7Xfn6TWXGn9bArf3LN8B7DxLsQzMZBuNY/E/heYbn6eVJmxR1JGSpMGaS42/viQ5FDi0LT6Y5MYpHnJT4MdTPMZMelq8vzVLgUzAvL/Gc9yCincOfp4X1PWdo/qN+ecHHch8N5N15Bz8vzpt+nhv8/H/2ZMG/G83Z6/NHPjMDuzazIH3NhXz7XvBTJr2+nEuNf5WAlv3LG/VylZTVccBx03XSZMsq6ql03W8QZtv8cL8i9l4B8t4B2u+xQvzM+ZZYB05h3hdRue1GZ3XZmRel9EN4trMmYQvwOXAdkm2TbI2sD9wzizHJEnSXGAdKUmasjlz56+qHktyOHAesAZwQlVdN8thSZI066wjJUnTYc40/gCq6lzg3Bk+7bR1j5kh8y1emH8xG+9gGe9gzbd4YX7GPOOsI+cUr8vovDaj89qMzOsyumm/Nqmq6T6mJEmSJGmOmUvP/EmSJEmSBmRBN/6S7JXkxiQrkhwxwvp1kpze1l+aZEnPuiNb+Y1JXjWX403ynCQXJnkwyWdnItYpxvsbSa5Ick173X2Ox7tTkuVt+m6S189EvFOJuWf9Nu1z8UdzOd4kS5L8V891/ru5HG9b9+Ik30lyXfssrztX403yWz3XdnmSJ5LsMIfjXSvJSe263pDkyEHHuhgMok4c75jzxXRfmyRbp6uXr29/I949c+9meg3ic9PWrZHkqiRfHfy7mH4D+v+0UZIzk3yv/e37lZl5N9NrQNfmD9v/pWuTfGkm6tzpNoU6cdTv+Ule1urKFUmOSZJxA6mqBTnRPRD/feB5wNrAd4Hth23zDuDv2vz+wOltfvu2/TrAtu04a8zheJ8FvBL4feCz8+D6vhR4bpt/EbByjsf7TGDNNr8FcOfQ8lyNuWf9mcA/An80l+MFlgDXzsRnd5riXRO4GnhJW37OXP4bMWybXwa+P8ev71uA09r8M4FbgSUz+flYaNMU/z1GrBP7OeZ8mAZ0bbYAdmzbrA/8h9dm9e9SwHuBU4Gvzvb7nCvXBTgJ+N02vzaw0Wy/17lwbYAtgVuA9dp2ZwAHzfZ7ncHrMur3fOAyYBcgwNeAvceLZSHf+dsJWFFVN1fVo8BpwL7DttmX7j8adF+S92gt5n3pvng8UlW3ACva8eZkvFX1UFX9G/DTAcfYayrxXlVV/9nKrwPWS7LOHI734ap6rJWvC8zUg7JT+QyT5HV0fyxnKiPglOKdBVOJd0/g6qr6LkBV3V1Vj8/heHsd0PYdtKnEW8CzkqwJrAc8CvxkBmJeyAZRJ/ZzzPlg2q9NVa2qqisBquoB4Aa6L7DzzUC+SyXZCng18A8z8B4GYdqvS5INgf8BHA9QVY9W1X0z8F6m26C+f69J931xTbofBf+T+WXav+cn2QLYoKouqa4leDLwuvECWciNvy2B23uW7+Dpf3if3KZ9ub+f7hf8fvadblOJdzZMV7xvAK6sqkcGFOfTYmkmFG+SnZNcB1wD/H5PY3CQJh1zkmcD7wc+OANxPi2WZqKfiW1bF6BvJfm1QQfL1OL970AlOS/JlUn+ZI7H2+vNwJcGFOOIsTQTifdM4CFgFfAD4G+q6p5BB7zADaJOnI26chAG+n2hdd16KXDpNMY8UwZ1bT4F/AnwxPSHPCMGcV22Be4CvtDqwn9I8qzBhD9Q035tqmol8Dd09cEq4P6q+sZAoh+cQXzP37IdZ6xjPs1CbvxpjkvyQuDjwO/NdizjqapLq+qFwMuBI+dBX/OjgKOr6sHZDqRPq4BtquqltK5ASTaY5ZjGsiZdF4zfaq+vT7LH7IY0viQ7Aw9X1bWzHcs4dgIeB55L94XofUmeN7shSRPXfoj7MvCeqvLuNZDkNcCdVXXFbMcyx6wJ7Agc2+rCh4B5+xztdEqyMd1dsW3p6oVnJfnt2Y1q/lrIjb+VwNY9y1u1shG3abeRNwTu7nPf6TaVeGfDlOJtXT7OAt5WVd8feLTTdH2r6gbgQbpnFQdtKjHvDPxVkluB9wB/mm6A6DkZb+vicTdA+0Lwfbq7a3MyXrpf1y6uqh9X1cN0Y6/tOIfjHbI/M3PXb7VYmonE+xbg61X1s6q6E/g2sHTgES9sg6gTZ6OuHISBfF9IshZdw++LVfWVgUQ+eIO4Nq8AXtvqp9OA3ZP830EEP0CDuC53AHdU1dAd4jMZfL0yCIO4Nv8TuKWq7qqqnwFfAX51INEPziC+569sxxnrmE833kOB83Wi+wXlZrpfCYYerHzhsG0OY/UHK89o8y9k9QdOb2bwyRwmHW/P+oOYuYQvU7m+G7Xt/9c8+Txsy1MJX36erp/5pnM55mHbHMXMJHyZyjXejKceeH8e3R+vTeZwvBsDV9KSAQH/Arx6rsbblp/RruvzBv1ZmIbr+37gC23+WcD1wItnIu6FOk3x32PEOrGfY86HaUDXJnTP33xqtt/fXLs2w/bdlfmZ8GUg1wX4V+AX2/xRwF/P9nudC9eG7gft6+jq3NA9F/fO2X6vM3VdetYfxPgJX/YZN5bZvhgDvtD70GXY+j7wZ63sQ8Br2/y6dJkQV7SL97yeff+s7XcjfWTOmQPx3grcQ3dX6g5mIKvYZOMF/pyuO8Pynunn5nC8b21/dJbTfeF/3Xz4DPcc4yhmoPE3xWv8hmHX+Dfncrxt3W+3mK8F/moexLsrcMlMfXan+Hl4diu/jq7h98czGfdCnab4+RmxThzpmPNxmu5rQ9cdvOiyAg/Vc+N+KZuL0yA+Nz3rd2UeNv4GdV2AHYBl7XPz/4CNZ/t9zqFr80Hge3R17inAOrP9Pmf4utzKCN/z6XrFXNuO+Vkg48WRtqMkSZIkaQFbyM/8SZIkSZIaG3+SJEmStAjY+JMkSZKkRcDGnyRJkiQtAjb+JEmSJGkRsPGnRSfJe5I8c5qO9YIky5NcleQXJrH/0iTHTGK/c5NsNNH9BiHJv09h34OSPHc645Ekjc16cLXjPDfJmRPc56Akn23zv5/kbeNs/w9Jth/rOJOR5MHJ7qvFy6EetOAkCd1n+4lR1t8KLK2qH0/DuY6gGwD+I1M91nyTZM2qemyKx7iIbgzCZdMTlSTJenCwkhxEd/0On83jJHmwqp49lRi0+HjnTwtCkiVJbkxyMt1gl1snOTbJsiTXJflg2+5dwHOBC5Nc2Mr2TPKdJFcm+cckT/tDmmSHJJckuTrJWUk2TrIP8B7gD4aONWyfB5P8dTv/vyTZKclFSW5O8tq2za5Jvtrmf739ejr0C+r6SbZIcnEruzbJr7Vtb02yaXvfNyT5fDvPN5Ks17Z5eYt3eYvj2hFiPC3Jq3uWT0zyxnbcf23X5Mokv9oT778mOYdu8O0nf3lM8uwkF7Ttr0myb8+/zdNiTPJGusFJv9hiXC/Jx5Jc3+L+m0l9GCRpEbIenHQ9uGSoPN2duK8k+XqSm5L8Vc92Byf5jySXAa/oKT8qyR+luwN62bDjXtPmL0qydJzjnNjqxSevXXsdsW4d9h5GvEbSiGZ7tHsnp+mYgCXAE8AuPWWbtNc1gIuAF7flW4FN2/ymwMXAs9ry+4G/GOH4VwO/3uY/BHyqzR9Fd+dqpJgK2LvNnwV8A1gLeAmwvJXvCny1zf8T8Io2/2xgTeB9wJ/1vI/1e99De9+PATu08jOA327z1wK/0uY/Blw7QoyvB05q82sDtwPrAc8E1m3l2wHLeuJ9CNi25xgPttc1gQ16rusKIOPEeBHdr54AzwFu5KkeCRvN9ufKycnJab5M1oOTrgeXDJUDBwE3AxsC6wK3AVsDWwA/ADZrdeW3gc8Of//A8qH6sV3HP2/zF9H92DnWcU4E3tgT15h167BtRrxGTk4jTd7500JyW1Vd0rO8X5IrgauAFwJP628P7NLKv51kOXAg8PO9GyTZkK4h8q1WdBLwP/qI51Hg623+GuBbVfWzNr9khO2/DXyy/Sq7UXVdKi8HDk5yFPDLVfXACPvdUlXL2/wVwJJ0z0GsX1XfaeWnjhLj14DdkqwD7A1cXFX/RVc5f779avmPrH7tLquqW0Y4VoC/THI18C/AlsDmo8U4wv73Az8Fjk/yv4CHR4lZkjQy68GJ14PDXVBV91fVT+l6uPw8sDNwUVXdVVWPAqePsu8ZwJvb/JtH2K7f4/Qaq24d0s81kgC7fWpheWhoJsm2wB8Be1TVi4F/pvsVb7gA51fVDm3avqoOmaZ4flZVQw/VPgE8AlDdMxhrDt+4qj4G/C7dnbdvJ3lBVV1MV8GuBE7MyA+VP9Iz//hIxx5Nq9wuAl7F6hXVHwI/ovt1dindL5RDHmJkv0X3a+bLqmqHtv/QNR83xlbJ7wScCbyGp74wSJL6Yz04wXpwmo91Ol2D+78DVVU3TWDfx2jfy5M8g6fq3bHqVuhO1M81kgAbf1q4NqCrBO9PsjndXa0hDwDrt/lLgFckeT5Akme1P9pPqqr7gXt7+tC/FfgW0yzJL1TVNVX1cbpf8V6Q5OeBH1XV54F/AHbs51hVdR/wQJKdW9H+Y2x+OnAw8Gs81eDaEFjVKui30nUjGc+GwJ1V9bMkuzHsl+NRPPlv0Z4x2bCqzqVrfL6kj/0lSSOzHuy/HhzPpcCvJ3lOkrWAN41yzu/TNRj/NyPf1RvrOLcCL2vzr6XrgQN91K2TvUZanKbyy4g0Z1XVd5NcBXyP7jm2b/esPg74epL/rKrd0mXb+lLr+gjw58B/DDvkgcDfpUuNfTNdY2m6vaf9YX8CuI6uS+b+wB8n+RnwIDCRX/MOoeu6+QRdJX3/KNt9AzgFOLt1QwH4W+DL7dfDrzP63b5eXwT+qXUVXUZ37cdzIt11/S+6LyZnJ1mX7pfo9/axvyRpBNaDQP/14JiqalXrUvkd4D66Z/tGczrw18C2EzzO5+nqwO+yer3bT926K5O/RlpkHOpBWqCSPLuqhrKFHQFsUVXvnuWwJEmaEdaD0tN5509auF6d5Ei6/+e30WUxkyRpsbAelIbxzp8kSZIkLQImfJEkSZKkRcDGnyRJkiQtAjb+JEmSJGkRsPEnSZIkSYuAjT9JkiRJWgRs/EmSJEnSIvD/AVhJBYSV4UV6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_missing(df_miss_ind, df_miss_var)\n",
    "# Free memory\n",
    "df_miss_ind = None\n",
    "df_miss_var = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last log reports 676'575 remaining variants. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minor allele frequencies\n",
    "\n",
    "## plink1.9 report\n",
    "\n",
    "* A1: minor allele (usually)\n",
    "* A2: major allele (usually)\n",
    "* MAF: frequency of A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "run_plink: command '--freq', the file 'data/plink/host_geno_clean_ind.frq' already exists (force is set to False).\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASf0lEQVR4nO3df6zdd33f8eeLuIG0IyQkbhTZpg6tpSpElJLb2FWnaSUjcdIpjkqEErrFRSlWS9DoQBqhVMoKkxa2qhnRaLaIWDhrg0npWDwGeF4IqvjDJjdtSEhSlkv4EVtAXDskBTZY6Ht/nI/Tw+V87j227z33Xt/nQzq63+/7+/l+P5+PTnJf9/vjHKeqkCRplBct9QAkScuXISFJ6jIkJEldhoQkqcuQkCR1rVnqASy0c889tzZu3LjUw5CkFeXBBx/8m6paO7t+yoXExo0bmZ6eXuphSNKKkuRro+pebpIkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpy5CQJHWdcp+4Phl3H/j6yPqbNr9iwiORpOXBMwlJUpchIUnqMiQkSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSusYKiSRfTfJIkoeSTLfay5PsS/JE+3l2qyfJbUlmkjyc5LVDx9ne2j+RZPtQ/eJ2/Jm2b+bqQ5I0GcdzJvGrVfWaqppq6zcB91XVJuC+tg5wBbCpvXYAt8PgFz5wM7AZuAS4eeiX/u3AW4b22zpPH5KkCTiZy03bgF1teRdw9VD9rhrYD5yV5HzgcmBfVR2tqmeAfcDWtu3MqtpfVQXcNetYo/qQJE3AuCFRwP9M8mCSHa12XlV9oy1/EzivLa8Dnhra92CrzVU/OKI+Vx8/IsmOJNNJpg8fPjzmlCRJ8xn3C/7+YVUdSvLTwL4kfz28saoqSS388Mbro6ruAO4AmJqaWtRxSNJqMtaZRFUdaj+fBj7O4J7Ct9qlItrPp1vzQ8CGod3Xt9pc9fUj6szRhyRpAuYNiSQ/leSlx5aBy4AvAnuAY08obQfubct7gOvbU05bgGfbJaO9wGVJzm43rC8D9rZtzyXZ0p5qun7WsUb1IUmagHEuN50HfLw9lboGuLuqPp3kAeCeJDcAXwPe2Np/ErgSmAG+B7wZoKqOJnkf8EBr996qOtqW3wp8GDgD+FR7AdzS6UOSNAHzhkRVPQn8woj6EeDSEfUCbuwcayewc0R9Grho3D4kSZPhJ64lSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVKXISFJ6jIkJEldhoQkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnLkJAkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSusYOiSSnJfmrJJ9o6xckOZBkJslHk5ze6i9u6zNt+8ahY7y71b+U5PKh+tZWm0ly01B9ZB+SpMk4njOJtwOPD62/H7i1qn4OeAa4odVvAJ5p9VtbO5JcCFwLvArYCvxxC57TgA8CVwAXAte1tnP1IUmagLFCIsl64NeAD7X1AK8DPtaa7AKubsvb2jpt+6Wt/TZgd1V9v6q+AswAl7TXTFU9WVU/AHYD2+bpQ5I0AeOeSfwH4F8Bf9fWzwG+XVXPt/WDwLq2vA54CqBtf7a1f6E+a59efa4+fkSSHUmmk0wfPnx4zClJkuYzb0gk+afA01X14ATGc0Kq6o6qmqqqqbVr1y71cCTplLFmjDa/AlyV5ErgJcCZwAeAs5KsaX/prwcOtfaHgA3AwSRrgJcBR4bqxwzvM6p+ZI4+JEkTMO+ZRFW9u6rWV9VGBjeeP1NVvwHcD1zTmm0H7m3Le9o6bftnqqpa/dr29NMFwCbg88ADwKb2JNPprY89bZ9eH5KkCTiZz0m8C3hHkhkG9w/ubPU7gXNa/R3ATQBV9ShwD/AY8Gngxqr6YTtLeBuwl8HTU/e0tnP1IUmagAz+YD91TE1N1fT09Ante/eBr4+sv2nzK05mSJK07CV5sKqmZtf9xLUkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnLkJAkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVKXISFJ6jIkJEldhoQkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlS17whkeQlST6f5AtJHk3yB61+QZIDSWaSfDTJ6a3+4rY+07ZvHDrWu1v9S0kuH6pvbbWZJDcN1Uf2IUmajHHOJL4PvK6qfgF4DbA1yRbg/cCtVfVzwDPADa39DcAzrX5ra0eSC4FrgVcBW4E/TnJaktOADwJXABcC17W2zNGHJGkC5g2JGvhOW/2J9irgdcDHWn0XcHVb3tbWadsvTZJW311V36+qrwAzwCXtNVNVT1bVD4DdwLa2T68PSdIEjHVPov3F/xDwNLAP+DLw7ap6vjU5CKxry+uApwDa9meBc4brs/bp1c+Zo4/Z49uRZDrJ9OHDh8eZkiRpDGOFRFX9sKpeA6xn8Jf/zy/qqI5TVd1RVVNVNbV27dqlHo4knTKO6+mmqvo2cD/wy8BZSda0TeuBQ235ELABoG1/GXBkuD5rn179yBx9SJImYJynm9YmOastnwG8HnicQVhc05ptB+5ty3vaOm37Z6qqWv3a9vTTBcAm4PPAA8Cm9iTT6Qxubu9p+/T6kCRNwJr5m3A+sKs9hfQi4J6q+kSSx4DdSf4N8FfAna39ncB/STIDHGXwS5+qejTJPcBjwPPAjVX1Q4AkbwP2AqcBO6vq0Xasd3X6kCRNQAZ/sJ86pqamanp6+oT2vfvA10fW37T5FSczJEla9pI8WFVTs+t+4lqS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVKXISFJ6jIkJEldhoQkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnLkJAkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpa96QSLIhyf1JHkvyaJK3t/rLk+xL8kT7eXarJ8ltSWaSPJzktUPH2t7aP5Fk+1D94iSPtH1uS5K5+pAkTcY4ZxLPA++sqguBLcCNSS4EbgLuq6pNwH1tHeAKYFN77QBuh8EvfOBmYDNwCXDz0C/924G3DO23tdV7fUiSJmDekKiqb1TVX7blvwUeB9YB24Bdrdku4Oq2vA24qwb2A2clOR+4HNhXVUer6hlgH7C1bTuzqvZXVQF3zTrWqD4kSRNwXPckkmwEfhE4AJxXVd9om74JnNeW1wFPDe12sNXmqh8cUWeOPmaPa0eS6STThw8fPp4pSZLmMHZIJPkHwJ8Dv1tVzw1va2cAtcBj+xFz9VFVd1TVVFVNrV27djGHIUmrylghkeQnGATEn1bVf23lb7VLRbSfT7f6IWDD0O7rW22u+voR9bn6kCRNwDhPNwW4E3i8qv5oaNMe4NgTStuBe4fq17ennLYAz7ZLRnuBy5Kc3W5YXwbsbdueS7Kl9XX9rGON6kOSNAFrxmjzK8A/Bx5J8lCr/R5wC3BPkhuArwFvbNs+CVwJzADfA94MUFVHk7wPeKC1e29VHW3LbwU+DJwBfKq9mKMPSdIEzBsSVfU5IJ3Nl45oX8CNnWPtBHaOqE8DF42oHxnVhyRpMvzEtSSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVKXISFJ6jIkJEldhoQkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnLkJAkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVLXvCGRZGeSp5N8caj28iT7kjzRfp7d6klyW5KZJA8nee3QPttb+yeSbB+qX5zkkbbPbUkyVx+SpMkZ50ziw8DWWbWbgPuqahNwX1sHuALY1F47gNth8AsfuBnYDFwC3Dz0S/924C1D+22dpw9J0oTMGxJV9RfA0VnlbcCutrwLuHqoflcN7AfOSnI+cDmwr6qOVtUzwD5ga9t2ZlXtr6oC7pp1rFF9SJIm5ETvSZxXVd9oy98EzmvL64CnhtodbLW56gdH1Ofq48ck2ZFkOsn04cOHT2A6kqRRTvrGdTsDqAUYywn3UVV3VNVUVU2tXbt2MYciSavKiYbEt9qlItrPp1v9ELBhqN36Vpurvn5Efa4+JEkTcqIhsQc49oTSduDeofr17SmnLcCz7ZLRXuCyJGe3G9aXAXvbtueSbGlPNV0/61ij+pAkTcia+Rok+Qjwj4Fzkxxk8JTSLcA9SW4Avga8sTX/JHAlMAN8D3gzQFUdTfI+4IHW7r1Vdexm+FsZPEF1BvCp9mKOPiRJEzJvSFTVdZ1Nl45oW8CNnePsBHaOqE8DF42oHxnVhyRpcvzEtSSpy5CQJHUZEpKkLkNCktRlSEiSugwJSVKXISFJ6jIkJEldhoQkqcuQkCR1GRKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnLkJAkdRkSkqQuQ0KS1GVISJK6DAlJUpchIUnqMiQkSV2GhCSpa81SD2AluPvA10fW37T5FRMeiSRNlmcSkqQuQ0KS1GVISJK6vCdxErxXIelUt+xDIslW4APAacCHquqWJR7SvAwPSaeKZR0SSU4DPgi8HjgIPJBkT1U9trQjOzG98DgRBo6kSVjWIQFcAsxU1ZMASXYD24AVGRILaSEDR9LKt1h/OC73kFgHPDW0fhDYPLtRkh3Ajrb6nSRfOsH+zgX+5gT3Xamc8+rgnE9xv3Hy8/2ZUcXlHhJjqao7gDtO9jhJpqtqagGGtGI459XBOZ/6Fmu+y/0R2EPAhqH19a0mSZqA5R4SDwCbklyQ5HTgWmDPEo9JklaNZX25qaqeT/I2YC+DR2B3VtWji9jlSV+yWoGc8+rgnE99izLfVNViHFeSdApY7pebJElLyJCQJHWtypBIsjXJl5LMJLlpxPYXJ/lo234gycbJj3JhjTHnf5TkL5M8n+SapRjjQhpjvu9I8liSh5Pcl2TkM+IryRhz/u0kjyR5KMnnkly4FONcSPPNeajdG5JUkhX/SOwY7/NvJjnc3ueHkvzWSXVYVavqxeAG+JeBVwKnA18ALpzV5q3Af2rL1wIfXepxT2DOG4FXA3cB1yz1mCcw318FfrIt/84qeY/PHFq+Cvj0Uo97sefc2r0U+AtgPzC11OOewPv8m8B/XKg+V+OZxAtf9VFVPwCOfdXHsG3Arrb8MeDSJJngGBfavHOuqq9W1cPA3y3FABfYOPO9v6q+11b3M/gMzko2zpyfG1r9KWClP7Uyzv/LAO8D3g/830kObpGMO+cFsxpDYtRXfazrtamq54FngXMmMrrFMc6cTyXHO98bgE8t6ogW31hzTnJjki8D/w74FxMa22KZd85JXgtsqKr/McmBLaJx/9t+Q7uU+rEkG0ZsH9tqDAnpBUn+GTAF/PulHsskVNUHq+pngXcBv7/U41lMSV4E/BHwzqUey4T9d2BjVb0a2MffXxU5IasxJMb5qo8X2iRZA7wMODKR0S2O1fb1JmPNN8k/Ad4DXFVV35/Q2BbL8b7Hu4GrF3VEi2++Ob8UuAj4bJKvAluAPSv85vW873NVHRn67/lDwMUn0+FqDIlxvupjD7C9LV8DfKbaHaEVarV9vcm8803yi8B/ZhAQTy/BGBfaOHPeNLT6a8ATExzfYphzzlX1bFWdW1Ubq2ojg3tPV1XV9NIMd0GM8z6fP7R6FfD4SfW41Hfrl+gJgSuB/83gKYH3tNp7GfwHBPAS4M+AGeDzwCuXeswTmPMvMbi++V0GZ02PLvWYF3m+/wv4FvBQe+1Z6jFPYM4fAB5t870feNVSj3mx5zyr7WdZ4U83jfk+/9v2Pn+hvc8/fzL9+bUckqSu1Xi5SZI0JkNCktRlSEiSugwJSVKXISFJ6jIkpAXQvmH0T4bW17Rv4vzErHb/Lcn+WbV/neTQ0Ld23jKpcUvzWdb/fKm0gnwXuCjJGVX1f4DXM+uTsEnOYvDp1+8keWVVPTm0+daq+sPJDVcaj2cS0sL5JINPMgNcB3xk1vZfZ/C9OrsZfFJWWvYMCWnh7AauTfISBv82x4FZ248Fx0fa8rB/OXS56fLFH6o0Hi83SQukqh5u/4rhdQzOKl6Q5DxgE/C5qqok/y/JRVX1xdbEy01aljyTkBbWHuAP+fFLTW8Ezga+0r6RdCM/fjYhLTuGhLSwdgJ/UFWPzKpfB2ytv/9G0ovxvoRWAENCWkBVdbCqbhuutUtQP8Pgq6qPtfsK8GySzRMdoHSc/BZYSVKXZxKSpC5DQpLUZUhIkroMCUlSlyEhSeoyJCRJXYaEJKnr/wPVSJ1UVQFazgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "path_freq = path_clean_ind\n",
    "run_plink('--freq', extension='frq', file=path_clean_ind, \n",
    "          out=path_clean_ind, plink2=False, log_name='freq_plink1');\n",
    "df_freq = pd.read_csv(path_freq+\".frq\", sep=\"\\s+\")\n",
    "df_freq.head(1)\n",
    "sns.distplot(df_freq.MAF.dropna(), kde=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage of variants to drop for different frequency thresholds:\n",
      "0.1\t86.9%\n",
      "0.05\t84.5%\n",
      "0.01\t78.6%\n",
      "0.005\t74.7%\n"
     ]
    }
   ],
   "source": [
    "frequencies = [0.1, 0.05, 0.01, 0.005]\n",
    "print(\"Percentage of variants to drop for different frequency thresholds:\")\n",
    "N = len(df_freq.MAF)\n",
    "for f in frequencies :\n",
    "    variants = df_freq.MAF[df_freq.MAF < f]\n",
    "    print(\"{}\\t{:.3}%\".format(f, len(variants)/N*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plink2 report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtype={'#CHROM':str, 'ID':str, 'REF':str, 'ALT':str, 'ALT_FREQS':float, 'OBS_CT':int}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "run_plink: command '--freq', the file 'data/plink/host_geno_clean_ind.afreq' already exists (force is set to False).\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEHCAYAAAC3Ph1GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARsElEQVR4nO3de5BkZX3G8e8TbppAgsiE2gJ0UfBCabnoiFheoiiGkJRgQgx4CSkpVwymvJXGS1XUJFZJJUJMitKsQliNKARFtoySEFyKYMnqICssSxQENIsrO15QyAVl+eWPPkvGYWa7d6a7Z172+6nqmnPec7rP752eeebM6XPOm6pCktSeX1rqAiRJC2OAS1KjDHBJapQBLkmNMsAlqVF7jnNjBx54YK1cuXKcm5Sk5l133XU/qKqJ2e1jDfCVK1cyNTU1zk1KUvOSfGeudg+hSFKjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo8Z6JeZiXLjhu3O2v+JZjxlzJZK0PLgHLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUqIEDPMkeSa5P8vlu/rAkG5LcmuSiJHuPrkxJ0my7sgf+RuDmGfNnAedU1eHAj4HTh1mYJGnnBgrwJIcAvw18rJsPcCxwSbfKWuCkURQoSZrboHvgfwO8HXigm380cHdV3d/NbwEOHnJtkqSd6BvgSX4H2FZV1y1kA0lWJ5lKMjU9Pb2Ql5AkzWGQPfDnAC9NcgfwaXqHTj4E7J9kx6X4hwB3zvXkqlpTVZNVNTkx8ZBBlSVJC9Q3wKvqnVV1SFWtBE4BvlRVrwTWAyd3q50GXDayKiVJD7GY88D/FHhLklvpHRM/bzglSZIGsUt3I6yqq4CruunbgKOHX5IkaRBeiSlJjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJatQggxo/IslXk3wjyU1J3te1X5Dk9iQbu8eq0ZcrSdphkBF57gOOrap7k+wFXJPki92yt1XVJaMrT5I0n74BXlUF3NvN7tU9apRFSZL6G+gYeJI9kmwEtgFXVNWGbtH7k9yQ5Jwk+8zz3NVJppJMTU9PD6lsSdJAAV5V26tqFXAIcHSSpwDvBJ4EPBM4gN4o9XM9d01VTVbV5MTExJDKliTt0lkoVXU3sB44vqq2Vs99wD/gCPWSNFaDnIUykWT/bvqRwHHAfyRZ0bUFOAnYNMpCJUm/aJCzUFYAa5PsQS/wL66qzyf5UpIJIMBG4IwR1ilJmmWQs1BuAI6ao/3YkVQkSRqIV2JKUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckho1yJBqj0jy1STfSHJTkvd17Ycl2ZDk1iQXJdl79OVKknYYZA/8PuDYqnoasAo4PskxwFnAOVV1OPBj4PTRlSlJmq1vgHcjz9/bze7VPQo4Frika19Lb2BjSdKYDHQMPMkeSTYC24ArgG8Dd1fV/d0qW4CD53nu6iRTSaamp6eHUbMkiQEDvKq2V9Uq4BDgaOBJg26gqtZU1WRVTU5MTCywTEnSbLt0FkpV3Q2sB54N7J9kx6j2hwB3Drk2SdJODHIWykSS/bvpRwLHATfTC/KTu9VOAy4bVZGSpIfas/8qrADWJtmDXuBfXFWfT7IZ+HSSvwSuB84bYZ2SpFn6BnhV3QAcNUf7bfSOh0uSloBXYkpSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoQUbkOTTJ+iSbk9yU5I1d+3uT3JlkY/c4YfTlSpJ2GGREnvuBt1bV15PsB1yX5Ipu2TlV9dejK0+SNJ9BRuTZCmztpu9JcjNw8KgLkyTt3C4dA0+ykt7wahu6pjckuSHJ+UkeNeTaJEk7MXCAJ9kX+Azwpqr6KfBh4PHAKnp76B+c53mrk0wlmZqenh5CyZIkGDDAk+xFL7w/WVWfBaiqu6pqe1U9AHyUeQY4rqo1VTVZVZMTExPDqluSdnuDnIUS4Dzg5qo6e0b7ihmrvQzYNPzyJEnzGeQslOcArwZuTLKxa3sXcGqSVUABdwCvG0mFkqQ5DXIWyjVA5lj0heGXI0kalFdiSlKjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaNciQaocmWZ9kc5Kbkryxaz8gyRVJbum+Oiq9JI3RIHvg9wNvraojgWOAM5McCbwDuLKqjgCu7OYlSWPSN8CramtVfb2bvge4GTgYOBFY2622FjhpVEVKkh5ql46BJ1kJHAVsAA6qqq3dou8DB83znNVJppJMTU9PL6JUSdJMAwd4kn2BzwBvqqqfzlxWVUVvdPqHqKo1VTVZVZMTExOLKlaS9P8GCvAke9EL709W1We75ruSrOiWrwC2jaZESdJcBjkLJcB5wM1VdfaMReuA07rp04DLhl+eJGk+ew6wznOAVwM3JtnYtb0L+ABwcZLTge8ALx9NiZKkufQN8Kq6Bsg8i1803HIkSYPySkxJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMGGVLt/CTbkmya0fbeJHcm2dg9ThhtmZKk2QbZA78AOH6O9nOqalX3+MJwy5Ik9dM3wKvqauBHY6hFkrQLFnMM/A1JbugOsTxqvpWSrE4ylWRqenp6EZuTJM200AD/MPB4YBWwFfjgfCtW1ZqqmqyqyYmJiQVuTpI024ICvKruqqrtVfUA8FHg6OGWJUnqZ0EBnmTFjNmXAZvmW1eSNBp79lshyaeAFwAHJtkCvAd4QZJVQAF3AK8bYY2SpDn0DfCqOnWO5vNGUIskaRd4JaYkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmN6hvg3aDF25JsmtF2QJIrktzSfZ13UGNJ0mgMsgd+AXD8rLZ3AFdW1RHAld28JGmM+gZ4VV0N/GhW84nA2m56LXDSkOuSJPWx0GPgB1XV1m76+8BBQ6pHkjSgRX+IWVVFb3DjOSVZnWQqydT09PRiNydJ6iw0wO9KsgKg+7ptvhWrak1VTVbV5MTExAI3J0mabaEBvg44rZs+DbhsOOVIkgY1yGmEnwK+AjwxyZYkpwMfAI5Lcgvw4m5ekjRGe/ZboapOnWfRi4ZciyRpF3glpiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUX0HdNiZJHcA9wDbgfuranIYRUmS+ltUgHdeWFU/GMLrSJJ2gYdQJKlRiw3wAv41yXVJVs+1QpLVSaaSTE1PTy9yc5KkHRYb4M+tqqcDvwWcmeT5s1eoqjVVNVlVkxMTE4vcnCRph0UFeFXd2X3dBlwKHD2MoiRJ/S04wJP8SpL9dkwDLwE2DaswSdLOLeYslIOAS5PseJ0Lq+ryoVQlSeprwQFeVbcBTxtiLZKkXeBphJLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0axqDGS+rCDd+ds/0Vz3rMmCuRpPFyD1ySGmWAS1KjFnUIJcnxwIeAPYCPVdUHhlLVEHhoRdLD3YIDPMkewLnAccAW4GtJ1lXV5mEVNwrzBftC+MdA0kw7y5dR5MVi9sCPBm7thlYjyaeBE4FlHeDDNMw/Btp9zfeL7X+R6mcxAX4w8J8z5rcAz5q9UpLVwOpu9t4k31zg9g4EfrDA57bKPu8GXrmLfX7lCGsZI9/nXfPYuRpHfhphVa0B1iz2dZJMVdXkEEpqhn3ePdjn3cMo+ryYs1DuBA6dMX9I1yZJGoPFBPjXgCOSHJZkb+AUYN1wypIk9bPgQyhVdX+SNwD/Qu80wvOr6qahVfZQiz4M0yD7vHuwz7uHofc5VTXs15QkjYFXYkpSowxwSWrUsgvwJMcn+WaSW5O8Y47l+yS5qFu+IcnK8Vc5XAP0+S1JNie5IcmVSeY8J7Ql/fo8Y73fS1JJmj7lbJD+Jnl59z7flOTCcdc4bAP8XD8myfok13c/2ycsRZ3DlOT8JNuSbJpneZL8bfc9uSHJ0xe1wapaNg96H4Z+G3gcsDfwDeDIWev8MfCRbvoU4KKlrnsMfX4h8Mvd9Ot3hz536+0HXA1cC0wudd0jfo+PAK4HHtXN//pS1z2GPq8BXt9NHwncsdR1D6HfzweeDmyaZ/kJwBeBAMcAGxazveW2B/7g5flV9TNgx+X5M50IrO2mLwFelCRjrHHY+va5qtZX1X93s9fSO+e+ZYO8zwB/AZwF/O84ixuBQfr7WuDcqvoxQFVtG3ONwzZInwv41W7614DvjbG+kaiqq4Ef7WSVE4GPV8+1wP5JVix0e8stwOe6PP/g+dapqvuBnwCPHkt1ozFIn2c6nd5f8Jb17XP3r+WhVfXP4yxsRAZ5j58APCHJl5Nc293ps2WD9Pm9wKuSbAG+APzJeEpbUrv6+75TzY/IsztJ8ipgEviNpa5llJL8EnA28EdLXMo47UnvMMoL6P2HdXWSp1bV3Uta1WidClxQVR9M8mzgE0meUlUPLHVhrVhue+CDXJ7/4DpJ9qT3r9cPx1LdaAx0S4IkLwbeDby0qu4bU22j0q/P+wFPAa5Kcge9Y4XrGv4gc5D3eAuwrqp+XlW3A9+iF+itGqTPpwMXA1TVV4BH0Lvh08PZUG9BstwCfJDL89cBp3XTJwNfqu7TgUb17XOSo4C/pxferR8bhT59rqqfVNWBVbWyqlbSO+7/0qqaWppyF22Qn+vP0dv7JsmB9A6p3DbOIodskD5/F3gRQJIn0wvw6bFWOX7rgD/szkY5BvhJVW1d8Kst9ae283xK+y16n2C/u2v7c3q/wNB7k/8JuBX4KvC4pa55DH3+N+AuYGP3WLfUNY+6z7PWvYqGz0IZ8D0OvcNGm4EbgVOWuuYx9PlI4Mv0zlDZCLxkqWseQp8/BWwFfk7vv6rTgTOAM2a8z+d235MbF/tz7aX0ktSo5XYIRZI0IANckhplgEtSowxwSWqUAS5JjTLAJalRBriWtSQndbeTfVI3v3L2rTqTfCTJxu5WrP/TTW9M8rJ5XvMfk9w+Y70zu/YtSW7sbvO5PsmDV/wm2T5j/Y1J3tYt2yfJ3yX5dneL0HU7ntct/7Pu9rA3dLdNfeaovlfa/XgvFC13pwLXdF/fM9cKVXUGQJLDgUuqatUAr/vmqvrcHO3Pq6q7k7wfeBe92/cC3DPP654F7AM8oaq2J3ktcGkX1M8FXgIcVVU/SzKBv3MaIvfAtWwl2ZdeCJ5O71LscfoKfe4Sl2Q/4FXAW6pqO0BVfZTeVXgvBFYA09W7nSpVNV2LuWxamsUA13J2InB5VX0L+GGSZwzxtc+ZcTjkyDmW/ya9+5PssN+sQygn07vZ1O1Vde+s507Ru0z8cuDx3ag05yZ53hDrlwxwLWun0hsIgO7rqUN87TdX1arusXlG+78n+R69myxdNKP9nhnrr6qqS/ptoKp+Sm90ljPo3THzkiSvHmIftJvzeJyWpSQHAMcCT01S9IboKno3Ahql5wH/Re+mRO8B3r6TdW8BDkuy76y98GcAn4QHBx1ZD6xPshn4A+AToyhcux/3wLVcnQx8oqoeW73byh4K3M4v3kt5JKrq58CbgNck2X8n690DXAj8VTcIBUleA9xdVdcmeXL3weoOq4DvjLB07WYMcC1XpwKXzmr7DPBO4IndKX87Hr8/7I1X1RZ6ty3ecRbK7GPg7+/a3w5sB27pDr2cCfxut2xfeqPM3JTkRuBwerdTlYbC28lKQ9INTns58KGqOn+p69HDnwEuSY3yQ0w9bCX5CL3xNGc6u6o+vhT1SMPmHrgkNcoPMSWpUQa4JDXKAJekRhngktSo/wNIxQognVnTawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "path_freq = path_clean_ind\n",
    "run_plink('--freq', extension='afreq', file=path_clean_ind, out=path_freq, log_name='freq_plink2');\n",
    "df_freq = pd.read_csv(path_freq+\".afreq\", sep=\"\\s+\", dtype=dtype)\n",
    "sns.distplot(df_freq.ALT_FREQS.dropna(), kde=False, norm_hist=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_freq = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hardy-Weinberg equilibrium\n",
    "\n",
    "* A1: Tested allele\n",
    "* AX: Non-A1 alleles (comma-separated)\n",
    "* HOM_A1_CT: Homozygous-A1 genotype count\n",
    "* HET_A1_CT: Heterozygous-A1 genotype count\n",
    "* TWO_AX_CT: number of NONmissing calls with no A1 copies\n",
    "* O(HET_A1): Observed heterozygous-major frequency\n",
    "* E(HET_A1): Expected heterozygous-major frequency\n",
    "* P: Hardy-Weinberg equilibrium exact test p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink2 --bfile data/plink/host_geno_clean_ind --out data/plink/host_geno_clean_ind_hardy --hardy '\n",
      "Log written to 'data/plink/log/hardy.log'\n"
     ]
    }
   ],
   "source": [
    "path_hardy = path_clean_ind+'_hardy'\n",
    "run_plink('--hardy', file=path_clean_ind, out=path_hardy, \n",
    "          extension='hardy', log_name='hardy');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "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>#CHROM</th>\n",
       "      <th>ID</th>\n",
       "      <th>A1</th>\n",
       "      <th>AX</th>\n",
       "      <th>HOM_A1_CT</th>\n",
       "      <th>HET_A1_CT</th>\n",
       "      <th>TWO_AX_CT</th>\n",
       "      <th>O(HET_A1)</th>\n",
       "      <th>E(HET_A1)</th>\n",
       "      <th>P</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1:14773:C:T</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>607</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0.004918</td>\n",
       "      <td>0.004906</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>rs571121669</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>605</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0.006568</td>\n",
       "      <td>0.006547</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>rs201013861</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>609</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.001639</td>\n",
       "      <td>0.001638</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   #CHROM           ID A1 AX  HOM_A1_CT  HET_A1_CT  TWO_AX_CT  O(HET_A1)  \\\n",
       "0       1  1:14773:C:T  C  T        607          3          0   0.004918   \n",
       "1       1  rs571121669  C  T        605          4          0   0.006568   \n",
       "2       1  rs201013861  G  A        609          1          0   0.001639   \n",
       "\n",
       "   E(HET_A1)    P  \n",
       "0   0.004906  1.0  \n",
       "1   0.006547  1.0  \n",
       "2   0.001638  1.0  "
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_hardy = pd.read_csv(path_hardy+'.hardy', sep=\"\\s+\")\n",
    "df_hardy.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtgAAADQCAYAAAA9HxZBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZRl1Vnn8e8vdEjIC4FAyyBNplvTxukwmpAW0KgTQ4QOxjQzxgyYkU7sFdQQJ050GVDXEPPiJDqKwUQUAwayDC+iMb2UBHsIDMuYJjQv4VWkQ17oHhJaIETEJBKf+ePsgktRL7e6Tlfdqvv9rHVXnbPPPmfvXbdq11P7nr1PqgpJkiRJ/XjKYldAkiRJWk4MsCVJkqQeGWBLkiRJPTLAliRJknpkgC1JkiT1yABbkiRJ6tGKxa5A3w4++OBavXr1YldDkubs+uuv/8eqWrnY9VhI9tmSlqqZ+uxlF2CvXr2a7du3L3Y1JGnOknxxseuw0OyzJS1VM/XZ3iIiSZIk9cgAW5IkSeqRAbYkSZLUIwNsSZIkqUcG2JIkSVKPlt0qIpK0kD5y7ZeelPZTRz9vEWoyXvy+SxpljmBLkiRJPTLAliRJknpkgC1JYyLJ+UnuS3LrFMd+KUklObjtJ8nZSXYkuTnJkQN5NyW5q702DaS/JMkt7Zyzk2RhWiZJo8UAW5LGx4eADZMTkxwOHAcM3tj8SmBte50KnNPyPhc4EzgaOAo4M8mB7ZxzgDcOnPeksiRpHBhgS9KYqKprgAemOHQW8CtADaRtBC6szjbggCSHAscDW6vqgap6ENgKbGjH9q+qbVVVwIXAiXuzPZI0qgywJWmMJdkI7Kqqz046dBhwz8D+zpY2U/rOKdIlaey4TJ8kjakkzwB+le72kIUs91S620543vNcWk/S8jPUCHaSL7SJKzcl2d7Snptka5vksnXiHrw+J8ZMV4YkqRffCawBPpvkC8Aq4IYk/w7YBRw+kHdVS5spfdUU6U9SVedW1fqqWr9y5cqemiJJo2Mut4j8SFW9qKrWt/3TgSurai1wZduHfifGTFeGJGmequqWqvq2qlpdVavpbus4sqq+DGwBTmmDJscAD1XVvcAVwHFJDmx9+HHAFe3Y15Ic0wZJTgE+tigNk6RFNp97sDcCF7TtC3h8MkufE2OmK0OSNEdJLgI+Dbwgyc4km2fIfjlwN7AD+GPgTQBV9QDwTuC69npHS6Pl+WA753PAx/dGOyRp1A17D3YBf5OkgD+qqnOBQ9qIBcCXgUPadp8TY6Yr4wm8n0+SZldVJ89yfPXAdgGnTZPvfOD8KdK3A0fMr5aStPQNG2D/YFXtSvJtwNYkfz94sKqqBd97zUxltID/XID169fv1XpIkiRJMxnqFpGq2tW+3gd8lO4e6q+02ztoX+9r2fucGDNdGZIkSdJImjXATvLMJM+e2Kab0HIr3QSYiZVANvH4ZJY+J8ZMV4YkSZI0koa5ReQQ4KNt5bwVwEeq6hNJrgMubZNkvgi8tuW/HDiBbpLLI8AboJsYk2RiYgw8eWLMh4D96CbFTEyMec80ZUiSJEkjadYAu6ruBr53ivT7gWOnSO9tYsx0ZUiSJEmjykelS5IkST0ywJYkSZJ6ZIAtSZIk9cgAW5IkSeqRAbYkSZLUIwNsSZIkqUcG2JIkSVKPDLAlaUwkOT/JfUluHUj77SR/n+TmJB9NcsDAsTOS7EhyZ5LjB9I3tLQdSU4fSF+T5NqWfkmSfReudZI0OgywJWl8fAjYMCltK3BEVX0P8A/AGQBJ1gEnAS9s5/xBkn2S7AN8AHglsA44ueUFeC9wVlU9H3gQ2Lx3myNJo8kAW5LGRFVdAzwwKe1vqurRtrsNWNW2NwIXV9U3qurzwA7gqPbaUVV3V9U3gYuBjUkCvBy4rJ1/AXDiXm2QJI0oA2xJ0oSfAT7etg8D7hk4trOlTZd+EPDVgWB9Il2Sxo4BtiSJJL8GPAr86QKUdWqS7Um27969e28XJ0kLzgBbksZcktcDrwJeV1XVkncBhw9kW9XSpku/HzggyYpJ6U9SVedW1fqqWr9y5cre2iFJo8IAW5LGWJINwK8Ar66qRwYObQFOSvK0JGuAtcBngOuAtW3FkH3pJkJuaYH5VcBr2vmbgI8tVDskaZQYYEvSmEhyEfBp4AVJdibZDLwfeDawNclNSf4QoKpuAy4Fbgc+AZxWVd9q91i/GbgCuAO4tOUFeBvw1iQ76O7JPm8BmydJI2PF7FkkSctBVZ08RfK0QXBVvRt49xTplwOXT5F+N90qI5I01hzBliRJknpkgC1JkiT1yABbkiRJ6pEBtiRJktSjoQPsJPskuTHJX7X9NUmuTbIjySVtuSbakk6XtPRrk6weuMYZLf3OJMcPpG9oaTuSnD6QPmUZkiRJ0qiaywj2W+iWZJrwXuCsqno+8CCwuaVvBh5s6We1fCRZR7de6guBDcAftKB9H+ADwCuBdcDJLe9MZUiSJEkjaagAO8kq4MeAD7b9AC8HLmtZLgBObNsb2z7t+LEt/0bg4qr6RlV9HthBt5zTUcCOqrq7qr4JXAxsnKUMSZIkaSQNO4L9e3RP+vq3tn8Q8NX2wAGAncBhbfsw4B6Advyhlv+x9EnnTJc+UxlPkOTUJNuTbN+9e/eQTZIkSZL6N2uAneRVwH1Vdf0C1GePVNW5VbW+qtavXLlysasjSZKkMTbMkxxfCrw6yQnA04H9gfcBByRZ0UaYVwG7Wv5dwOHAziQrgOcA9w+kTxg8Z6r0+2coQ5IkSRpJs45gV9UZVbWqqlbTTVL8ZFW9DrgKeE3Ltgn4WNve0vZpxz9ZVdXST2qrjKwB1gKfAa4D1rYVQ/ZtZWxp50xXhiRJkjSS5rMO9tuAtybZQXe/9Hkt/TzgoJb+VuB0gKq6DbgUuB34BHBaVX2rjU6/GbiCbpWSS1vemcqQJEmSRtIwt4g8pqquBq5u23fTrQAyOc/XgZ+c5vx3A++eIv1y4PIp0qcsQ5IkSRpVPslRksZEkvOT3Jfk1oG05ybZmuSu9vXAlp4kZ7cHfd2c5MiBcza1/Hcl2TSQ/pIkt7Rzzm7LrUrS2DHAlqTx8SG6B30NOh24sqrWAle2fege/rW2vU4FzoEuIAfOBI6m+4TxzImgvOV548B5k8uSpLFggC1JY6KqrgEemJQ8+HCwyQ8Nu7A62+hWdToUOB7YWlUPVNWDwFZgQzu2f1Vta5PUL8SHg0kaUwbYkjTeDqmqe9v2l4FD2vZcHw52WNuenP4kPhxM0nJngC1JAqCNPNcClOPDwSQtawbYkjTevtJu76B9va+lT/dwsJnSV02RLkljxwBbksbb4MPBJj807JS2msgxwEPtVpIrgOOSHNgmNx4HXNGOfS3JMW31kFPw4WCSxtSc1sGWJC1dSS4CXgYcnGQn3Wog7wEuTbIZ+CLw2pb9cuAEYAfwCPAGgKp6IMk76Z7CC/COqpqYOPkmupVK9gM+3l6SNHYMsCVpTFTVydMcOnaKvAWcNs11zgfOnyJ9O3DEfOooScuBt4hIkiRJPTLAliRJknpkgC1JkiT1yABbkiRJ6pEBtiRJktQjA2xJkiSpRwbYkiRJUo8MsCVJkqQeGWBLkiRJPTLAliRJknpkgC1JkiT1aNYAO8nTk3wmyWeT3JbkN1r6miTXJtmR5JIk+7b0p7X9He346oFrndHS70xy/ED6hpa2I8npA+lTliFJkiSNqmFGsL8BvLyqvhd4EbAhyTHAe4Gzqur5wIPA5pZ/M/BgSz+r5SPJOuAk4IXABuAPkuyTZB/gA8ArgXXAyS0vM5QhSepRkv/RBlFuTXJRG1zpbSBFksbJrAF2dR5uu09trwJeDlzW0i8ATmzbG9s+7fixSdLSL66qb1TV54EdwFHttaOq7q6qbwIXAxvbOdOVIUnqSZLDgP8OrK+qI4B96AZEehlIWci2SNIoGOoe7DbSfBNwH7AV+Bzw1ap6tGXZCRzWtg8D7gFoxx8CDhpMn3TOdOkHzVDG5PqdmmR7ku27d+8epkmSpCdaAeyXZAXwDOBe+htIkaSxMlSAXVXfqqoXAavoOsvv3qu1mqOqOreq1lfV+pUrVy52dSRpSamqXcD/Br5EF1g/BFxPfwMpkjRW5rSKSFV9FbgK+H7ggDbSAV3gvatt7wIOB2jHnwPcP5g+6Zzp0u+foQxJUk+SHEg3+rwG+HbgmXS3eOyt8vzUUdKyNswqIiuTHNC29wN+FLiDLtB+Tcu2CfhY297S9mnHP1lV1dJPapNj1gBrgc8A1wFr22Safenu39vSzpmuDElSf14BfL6qdlfVvwJ/AbyU/gZSnsBPHSUtd8OMYB8KXJXkZrpgeGtV/RXwNuCtSXbQfTR4Xst/HnBQS38rcDpAVd0GXArcDnwCOK3devIo8GbgCrrA/dKWlxnKkCT150vAMUme0e6lPpaur+5rIEWSxsqK2TJU1c3Ai6dIv5spJq9U1deBn5zmWu8G3j1F+uXA5cOWIUnqT1Vdm+Qy4AbgUeBG4Fzgr4GLk7yrpQ0OpHy4DX48QPfJI1V1W5KJgZRHaQMpC9oYSRoBswbYkqTlr6rOBM6clNzbQIokjRMflS5JkiT1yABbkiRJ6pEBtiRJktQjA2xJkiSpRwbYkiRJUo8MsCVJkqQeGWBLkiRJPTLAliRJknpkgC1JkiT1yABbkiRJ6pEBtiRJktQjA2xJkiSpRwbYkiSSHJDksiR/n+SOJN+f5LlJtia5q309sOVNkrOT7Ehyc5IjB66zqeW/K8mmxWuRJC0eA2xJEsD7gE9U1XcD3wvcAZwOXFlVa4Er2z7AK4G17XUqcA5AkucCZwJHA0cBZ04E5ZI0TgywJWnMJXkO8MPAeQBV9c2q+iqwEbigZbsAOLFtbwQurM424IAkhwLHA1ur6oGqehDYCmxYwKZI0kgwwJYkrQF2A3+S5MYkH0zyTOCQqrq35fkycEjbPgy4Z+D8nS1tunRJGisG2JKkFcCRwDlV9WLgn3n8dhAAqqqA6qOwJKcm2Z5k++7du/u4pCSNFANsSdJOYGdVXdv2L6MLuL/Sbv2gfb2vHd8FHD5w/qqWNl36E1TVuVW1vqrWr1y5steGSNIomDXATnJ4kquS3J7ktiRvaem9zS5P8pIkt7Rzzk6SmcqQJPWnqr4M3JPkBS3pWOB2YAsw0VdvAj7WtrcAp7T+/hjgoXYryRXAcUkObP31cS1NksbKMCPYjwK/VFXrgGOA05Kso9/Z5ecAbxw4b2JSzHRlSJL69QvAnya5GXgR8JvAe4AfTXIX8Iq2D3A5cDewA/hj4E0AVfUA8E7guvZ6R0uTpLGyYrYMbVTi3rb9T0nuoJu0shF4Wct2AXA18DYGZpcD29raqoe2vFsnOtskW4ENSa4G9m8z0UlyId1M9Y/PUIYkqUdVdROwfopDx06Rt4DTprnO+cD5/dZOkpaWOd2DnWQ18GLgWvqbXX5Y256czgxlSJIkSSNp6AA7ybOAPwd+saq+Nnisz9nl05mpDGekS5IkaVQMFWAneSpdcP2nVfUXLbmv2eW72vbk9JnKeAJnpEuSJGlUDLOKSOie7nVHVf3uwKFeZpe3Y19Lckwr65RJ15qqDEmSJGkkzTrJEXgp8NPALUluamm/Sjeb/NIkm4EvAq9txy4HTqCbXf4I8AboZpcnmZhdDk+cXf4m4EPAfnSTGz/e0qcrQ5IkSRpJw6wi8rdApjncy+zyqtoOHDFF+v1TlSFJkiSNKp/kKEmSJPXIAFuSJEnqkQG2JEmS1CMDbEmSJKlHBtiSJElSjwywJUmSpB4ZYEuSJEk9MsCWJJFknyQ3Jvmrtr8mybVJdiS5JMm+Lf1pbX9HO7564BpntPQ7kxy/OC2RpMVngC1JAngLcMfA/nuBs6rq+cCDwOaWvhl4sKWf1fKRZB1wEvBCYAPwB0n2WaC6S9JIMcCWpDGXZBXwY8AH236AlwOXtSwXACe27Y1tn3b82JZ/I3BxVX2jqj4P7ACOWpgWSNJoMcCWJP0e8CvAv7X9g4CvVtWjbX8ncFjbPgy4B6Adf6jlfyx9inOeIMmpSbYn2b579+4+2yFJI8EAW5LGWJJXAfdV1fULVWZVnVtV66tq/cqVKxeqWElaMCsWuwKSpEX1UuDVSU4Ang7sD7wPOCDJijZKvQrY1fLvAg4HdiZZATwHuH8gfcLgOZI0VhzBlqQxVlVnVNWqqlpNN0nxk1X1OuAq4DUt2ybgY217S9unHf9kVVVLP6mtMrIGWAt8ZoGaIUkjxRFsSdJU3gZcnORdwI3AeS39PODDSXYAD9AF5VTVbUkuBW4HHgVOq6pvLXy1JWnxGWBLkgCoqquBq9v23UyxCkhVfR34yWnOfzfw7r1XQ0laGrxFRJIkSeqRAbYkSZLUIwNsSZIkqUcG2JIkSVKPZg2wk5yf5L4ktw6kPTfJ1iR3ta8HtvQkOTvJjiQ3Jzly4JxNLf9dSTYNpL8kyS3tnLPbI3enLUOSJEkaZcOMYH8I2DAp7XTgyqpaC1zZ9gFeSbf26VrgVOAc6IJl4EzgaLpZ6WcOBMznAG8cOG/DLGVIkiRJI2vWALuqrqFb63TQRuCCtn0BcOJA+oXV2Ub3JLBDgeOBrVX1QFU9CGwFNrRj+1fVtvagggsnXWuqMiRJkqSRtaf3YB9SVfe27S8Dh7Ttw4B7BvLtbGkzpe+cIn2mMiRJkqSRNe9Jjm3kuXqoyx6XkeTUJNuTbN+9e/ferIokSZI0oz0NsL/Sbu+gfb2vpe8CDh/It6qlzZS+aor0mcp4kqo6t6rWV9X6lStX7mGTJEmSpPnb00elbwE2Ae9pXz82kP7mJBfTTWh8qKruTXIF8JsDExuPA86oqgeSfC3JMcC1wCnA789ShiRJQ/vItV96UtpPHf28RaiJpHExa4Cd5CLgZcDBSXbSrQbyHuDSJJuBLwKvbdkvB04AdgCPAG8AaIH0O4HrWr53VNXExMk30a1Ush/w8fZihjIkSZKkkTVrgF1VJ09z6Ngp8hZw2jTXOR84f4r07cARU6TfP1UZkqR+JTmcbhWnQ+jmu5xbVe9rS6xeAqwGvgC8tqoebM8reB/dgMojwOur6oZ2rU3Ar7dLv6uqLkCSxoxPcpQkPQr8UlWtA44BTkuyjn6feSBJY2NP78GWJC0TbUnUe9v2PyW5g27J1I10twhC9zyCq4G3MfDMA2BbkolnHryM9swDgCRb6R4edtGCNWaeprpfG7xnW9LcOIItSXpMktXAi+kmnvf1zIPJZbi0qqRlzQBbkgRAkmcBfw78YlV9bfBYn888cGlVScudAbYkiSRPpQuu/7Sq/qIl9/XMA0kaKwbYkjTm2qog5wF3VNXvDhyaeB4BPPmZB6ekcwztmQfAFcBxSQ5skxuPa2mSNFac5ChJeinw08AtSW5qab9Kv888kKSxYYAtSWOuqv4WyDSHe3nmwUKYbgWQUeDqJNJ48RYRSZIkqUcG2JIkSVKPDLAlSZKkHnkPtiRJS5T3dkujyRFsSZIkqUeOYEuSpKEs9Ij5VOU5Oq+lwABbkjR2RnlJP0lLnwF243/JkiRJ6oMBtiRJ0hw4uVSzMcCWJGkW3lKipWKpBf9Lrb7DMsCWJGmReHuitDwZYEuSNEJGYbS8jzos15FJaRgjH2An2QC8D9gH+GBVvWeRqyRJmoZ99miYS4C8nAPe5R7kL+Q/Y8v9e9m3kQ6wk+wDfAD4UWAncF2SLVV1++LWTJI0mX320rTQI+ajcFvMXOowl+9PH//YjML3Z1QtpSB/pANs4ChgR1XdDZDkYmAjsCCd9VJ6IyVpBCxqn63FM98gfW8FsaNsodu8t75vo3DdPurQd2w36gH2YcA9A/s7gaMXqS6PWS6/3JreXEcy5vKLubdGTqa6xt76J3FU/vmc70jPqLRjGRnJPlsahn/b1adRD7CHkuRU4NS2+3CSO/fgMgcD/9hfrUaKbZuj1+3l/HM4f+j2zaUO861vT9cdifdub1zjdXvetn8/z6KXBPvsWS3ntsHybp9tW6L2sN+ets8e9QB7F3D4wP6qlvYEVXUucO58CkqyvarWz+cao8q2LV3LuX22bVmyz+7Bcm4bLO/22balq+/2PaWvC+0l1wFrk6xJsi9wErBlkeskSZqafbYkMeIj2FX1aJI3A1fQLfl0flXdtsjVkiRNwT5bkjojHWADVNXlwOULUNS8Pq4ccbZt6VrO7bNty5B9di+Wc9tgebfPti1dvbYvVdXn9SRJkqSxNur3YEuSJElLylgF2Ek2JLkzyY4kp09x/GlJLmnHr02yeuFrueeGaN9bk9ye5OYkVyZZMkuCzda2gXw/kaSSLJmZzsO0Lclr23t3W5KPLHQd52OIn8vnJbkqyY3tZ/OExajnnkhyfpL7ktw6zfEkObu1/eYkRy50HZei+fTVSc5o6XcmOX4h6z2MPW1bkoPa78nDSd6/0PUexjza9qNJrk9yS/v68oWu+zDm0b6jktzUXp9N8p8Xuu6zmW981Prxh5P88kLVeVjzeN9WJ/mXgffuD+dUcFWNxYtuws3ngO8A9gU+C6yblOdNwB+27ZOASxa73j2370eAZ7Ttn18q7RumbS3fs4FrgG3A+sWud4/v21rgRuDAtv9ti13vntt3LvDzbXsd8IXFrvcc2vfDwJHArdMcPwH4OBDgGODaxa7zqL/m01e3n5/PAk8D1rTr7LPYbeqpbc8EfhD4OeD9i92Wntv2YuDb2/YRwK7Fbk/P7XsGsKJtHwrcN7E/Cq/5tG3g+GXAnwG/vNjt6fF9Wz1d3z7Ma5xGsB97hG9VfROYeITvoI3ABW37MuDYJFnAOs7HrO2rqquq6pG2u41ujdqlYJj3DuCdwHuBry9k5eZpmLa9EfhAVT0IUFX3LXAd52OY9hWwf9t+DvD/FrB+81JV1wAPzJBlI3BhdbYBByQ5dGFqt2TNp6/eCFxcVd+oqs8DO9r1RsUet62q/rmq/pbR7d/m07Ybq2ri9/42YL8kT1uQWg9vPu17pKoebelPp+vzRsm84qMkJwKfp3vvRs2ixX7jFGBP9Qjfw6bL034ZHgIOWpDazd8w7Ru0mW5kbSmYtW3to/fDq+qvF7JiPRjmffsu4LuSfCrJtiQbFqx28zdM+94O/LckO+lWn/iFhanagpjr76Xm11eP+vd7Of8d6qttPwHcUFXf2Ev13FPzal+So5PcBtwC/NxAwD0K9rhtSZ4FvA34jQWo556Y78/lmnb74v9N8kNzKXjkl+lT/5L8N2A98J8Wuy59SPIU4HeB1y9yVfaWFXS3ibyM7lOHa5L8x6r66qLWqj8nAx+qqt9J8v3Ah5McUVX/ttgVk7RwkryQ7lPI4xa7Ln2rqmuBFyb5D8AFST5eVaP6acRcvB04q6oeXjof+A/tXuB5VXV/kpcAf5nkhVX1tWFOHqcR7GEe4ftYniQr6D6uvn9Bajd/Qz2iOMkrgF8DXj2CIwTTma1tz6a7b+/qJF+gu9d1S5bGRMdh3redwJaq+tf2sfc/0AXcS8Ew7dsMXApQVZ+m+wj14AWp3d431O+lnmA+ffWof7+X89+hebUtySrgo8ApVfW5vV7buevlvauqO4CH6f5mjYr5tO1o4Lfa395fBH413cOmRsUet63danY/QFVdT3cv93cNW/A4BdjDPMJ3C7Cpbb8G+GS1O92XgFnbl+TFwB/RBddL6T7eGdtWVQ9V1cFVtbqqVtPdX/7qqtq+ONWdk2F+Lv+SbvSaJAfT/YLfvZCVnIdh2vcl4FiANrrzdGD3gtZy79kCnJLOMcBDVXXvYldqxM2nr94CnNRWBVhD94/oZxao3sNYzn+H9rhtSQ4A/ho4vao+tWA1npv5tG9NC9xIt3rXdwNfWJhqD2WP21ZVPzTwt/f3gN+sqlFa5WY+79vKJPsAJPkOuv5k+L+9ezo7cim+6Gb0/wPdfyG/1tLeQReMQfeH/c/oJsZ8BviOxa5zz+37P8BXgJvaa8ti17mvtk3KezVLZBWRId+30N0Cczvd/XsnLXade27fOuBTdLO7bwKOW+w6z6FtF9F9jPivdJ80bKZb5eHnBt67D7S237KUfi5H/Gdm2r6a7hO6zwF3Aq9c7Lb03LYv0E2qfbj9vD1pNaWl2Dbg14F/HvjbdBMjuFrSPNr303QTAG8CbgBOXOy29PlzOXCNtzNiq4jM8337iUnv24/PpVyf5ChJkiT1aJxuEZEkSZL2OgNsSZIkqUcG2JIkSVKPDLAlSZKkHhlgS5IkST0ywNZekeThSfuvTzKvtTGTXD3sw2PausP/mOTAtn9okkrygwN5dieZ9hHESV6d5PRZynlZkr8atg2jZqnXX1I/7LOXhqVe/3FigK2RNLEo/56qbv3JbcD3t6QfAG5sX0nyAronNU37hLSq2lJV75lPPWYzsYi9JC1l9tnSExlga8El+fEk1ya5Mcn/SXJIS397kg8n+RTw4ST7Jbk4yR1JPgrs1/L9TJLfG7jeG5OcNUVRf0frnNvXs3hi5/2pdv7KJH+e5Lr2emlLf2wEJ8mHkpyd5O+S3J3kNQPl7J/kr5PcmeQPkzylnXNckk8nuSHJnyV5Vkv/QpL3JrkB+Mkk35fk5iQ3JfntJLdO8T17WZJrpipnUr5tSV44sH91kvVJjmp1ubG14QWzvlGShH22fbb2hAG29pb9WudzU5Kb6J6aNOFvgWOq6sXAxcCvDBxbB7yiqk4Gfh54pKr+A3Am8JKW51Lgx5M8te2/ATh/ijp8isc766OAjwKHt/0foOvMAd4HnFVV30f35KYPTtOmQ4EfBF4FDI6SHAX8Qqv7dwL/Jd0jzX+9teVIYDvw1oFz7q+qI6vqYuBPgJ+tqhcB35qm7CnLmSLPJcBrofuIFTi0ukfG/z3wQ+17/j+B35yhHEnjxz7bPls9mtdHOtIM/qV1PkA3sgBM3Iu3CrikdSb7Ap8fOG9LVf1L2/5h4GyAqro5yc1t++EknwReleQO4KlVdcsUdbgOeHGSZ7Y8D7eRjOfTdda/0/K9AliXZOK8/SdGLib5y6r6N+D2iRGc5jNVdXdr50V0HfrXaY8Ab9fdF/j0wDmXtPwHAO6MVS0AAAJYSURBVM+uqoljH6H7YzCVqcq5bFKeS4G/ofvj9tqB488BLkiyFijgqUjS4+yz7bPVI0ewtRh+H3h/Vf1H4GeBpw8c++chr/FB4PV0IyF/ApDktIERmG+vqkeAu4CfAW5o520DTgC+DbizpT2FbnTmRe11WFU9YcJP842B7Qxs16R81Y5vHbjmuqravAftnHzdJ+wn+c8DbV5fVbuA+5N8D/BfaX8UgHcCV1XVEcCP88TvuSTNxD7bPltzZICtxfAcYFfb3jRDvmuAnwJIcgTwPRMHqupauo8Ofwq4qKV9YKBz/H8t698Bv8jjIxGfBt4CbGuTaqAbPfiFiWsneWwUZ0hHJVnT7q/7r3Qfp24DXtpGXkjyzCTfNfnEqvoq8E9Jjm5JJ82lnKr66ECbt7d8l9B9hPucqrq5pQ1+z18/x/ZJGm/22Y+3wz5bQzHA1mJ4O/BnSa4H/nGGfOcAz2ofKb4DuH7S8UuBT1XVgzNc41PAd/B4Z30D3cedfzeQ578D69uklduBnxu2Ic11wPuBO+g+Ov1oVe2m6xQvah+Tfhr47mnO3wz8cbvv8ZnAQ8OWM02+y+g6/UsH0n4L+F9JbsRbwyTNzduxzx5kn61Z5fF/CKWlJd1aoGdV1ZWLXZf5SPKsiY83063hemhVvWVSnpcBv1xV093rJ0kjzT5b48QRbC05SQ5I8g90k3KWdEfd/Fi7H+9W4IeAdy12hSSpL/bZGkeOYEuSJEk9cgRbkiRJ6pEBtiRJktQjA2xJkiSpRwbYkiRJUo8MsCVJkqQeGWBLkiRJPfr//27iGo3nIEYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(12, 3))\n",
    "sns.distplot(df_hardy.P, kde=False, ax=ax[0])\n",
    "ax[0].set_xlabel('Hardy-Weinberg p-val');\n",
    "\n",
    "sns.distplot(df_hardy.P[df_hardy.P <= 0.05], kde=False, ax=ax[1])\n",
    "ax[1].set_xlabel('Hardy-Weinberg p-val');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#df_hardy['COORD'] = df_hardy.ID.map(SNP_dic)\n",
    "#sns.scatterplot(x=df_hardy.COORD, y=df_hardy.P);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mendelian error rates\n",
    "\n",
    "[Quote](http://zzz.bwh.harvard.edu/plink/summary.shtml#mendel):\n",
    "> The `*.mendel` file contains all Mendel errors (i.e. one line per error); the `*.imendel` file contains a summary of per-individual error rates; the `*.fmendel` file contains a summary of per-family error rates; the `*.lmendel` file contains a summary of per-SNP error rates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink --bfile data/plink/host_geno_clean_ind --out data/plink/host_geno_clean_ind --mendel '\n",
      "Warning: Nonmissing nonmale Y chromosome genotype(s) present; many commands\n",
      "treat these as missing.\n",
      "Warning: Skipping --me/--mendel since there are no trios.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "run_plink(\"--mendel\", extension='mendel', file=path_clean_ind,\n",
    "          out=path_clean_ind, plink2=False, log_name='mendel');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there are **no trio**, there isn't any mendelian error rate to compute."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter MAF and variants in HW disequilibrium\n",
    "\n",
    "After this last step of filtering, we want to target to `path_clean`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 'plink2 --bfile data/plink/host_geno_clean_ind --out data/plink/host_geno_clean --maf 0.01 --hwe 1e-6 --make-bed '\n",
      "Log written to 'data/plink/log/filter_maf_hardy.log'\n"
     ]
    }
   ],
   "source": [
    "run_plink(\"--maf 0.01 --hwe 1e-6 --make-bed\", file=path_clean_ind, out=path_clean,\n",
    "          extension='bim', force=True, log_name='filter_maf_hardy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PCA pre-processing\n",
    "\n",
    "## Find intersecting individuals\n",
    "\n",
    "We first want to include only individuals that are in all the three datasets. We run the custom function `common.find_individuals`. We'll later use the command `--keep-fam <file>` of plink to filter the individuals. We could actually use `--keep <file>`, but in our case those two ids are the same (and the file for the former option has to have two columns).\n",
    "\n",
    "[Quote](https://www.cog-genomics.org/plink/1.9/filter):\n",
    "> --keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column, and removes all unlisted samples from the current analysis. --remove does the same for all listed samples.\n",
    ">\n",
    ">Similarly, --keep-fam and --remove-fam accept text files with family IDs in the first column, and keep or remove entire families."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "725 in clinical\n",
      "789 in viral\n",
      "611 in host genotypes\n",
      "Viral and clinical ids:\t611 in intersection, 725 in union\n",
      "Geno and clinical ids :\t724 in intersection, 790 in union\n",
      "Global intersection : 610\n",
      "Written to data/individuals/intersecting_ids\n"
     ]
    }
   ],
   "source": [
    "common.find_individuals(fam=PATH_HOST_CLEAN_DATA,\n",
    "                        clinical=PATH_CLINICAL_DATA,\n",
    "                        viral=PATH_VIRAL_DATA, \n",
    "                        output_path=PATH_INTERSECTING_INDIVIDUALS)\n",
    "# Write them into a file that can be read by the --keep command of plink\n",
    "#with open(PATH_INTERSECTING_INDIVIDUALS, 'w') as file:\n",
    "#    file.write(df_inds[ID_IGM_CLINICAL_DF].to_csv(header=False, index=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***\n",
    "<br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br /><br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Asian individuals\n",
    "\n",
    "We aim here to generate a dataset with only asian people with excluded X, Y and chromosome 6. We first compute the list of asian individuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "§"
   ]
  }
 ],
 "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
}