07-Impute-Chu/MendelImpute_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genotype imputation and phasing\n",
    "\n",
    "This notebook showcases our software [MendelImpute.jl](https://github.com/biona001/MendelImpute.jl), the fastest and most memory-efficient phasing/imputation software as of 2020. Our paper will be available on arxiv soon.\n",
    "\n",
    "## Package installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Uncomment to add necessary packages\n",
    "# using Pkg\n",
    "# Pkg.add(PackageSpec(url=\"https://github.com/OpenMendel/SnpArrays.jl.git\"))\n",
    "# Pkg.add(PackageSpec(url=\"https://github.com/OpenMendel/VCFTools.jl.git\"))\n",
    "# Pkg.add(PackageSpec(url=\"https://github.com/OpenMendel/MendelImpute.git\"))\n",
    "# Pkg.add([\"GeneticVariation\", \"Random\", \"SparseArrays\", \"Plots\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load necessary packages\n",
    "using MendelImpute\n",
    "using VCFTools\n",
    "using GeneticVariation\n",
    "using Random\n",
    "using SparseArrays\n",
    "using UnicodePlots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is imputation and phasing?\n",
    "\n",
    "+ Making standard SNP array (GWAS) data more like sequence data by increasing the number of variants that are available for study\n",
    "+ Inputs: \n",
    "    + Unphased (entries 0, 1, 2) genotypes at ~1 million markers \n",
    "    + A reference haplotype panel (entries 0, 1) at $\\ge 40$ million markers\n",
    "+ Output:\n",
    "    + Phased genotypes at all markers (if genotype is 1, you know which chromosome it is on)\n",
    "    \n",
    "# Method overview\n",
    "\n",
    "We take a data driven approach for genotype imputation and phasing. For each target genotype, we try to select 2 haplotype segments from the reference panel to match it up as best as possible within a small genomic window. Numerous interacting tactics are summarized in the cartoon illustration below. \n",
    "\n",
    "\n",
    "\n",
    "![fdsa](method.png)\n",
    "\n",
    "\n",
    "\n",
    "Overview of *MendelImpute*'s algorithm. (a) After alignment, imputation and phasing are carried out on short, non-overlapping windows of the typed SNPs. (b) Based on a least squares criterion, we find two unique haplotypes whose vector sum approximates the genotype vector on the current window. Once this is done, all reference haplotypes corresponding to these two unique haplotypes are assembled into two sets of candidate haplotypes. (c) We intersect candidate haplotype sets window by window, carrying along the surviving set and switching orientations if the result generates more surviving haplotypes. (d) After three windows the top extended chromosome possess no surviving haplotypes, but a switch to the second orientation in the current window allows $h_5$ to survive on the top chromosome. Eventually we must search for a break point separating $h_1$ from $h_2$ or $h_6$ between windows 3 and 4 (bottom panel)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example data\n",
    "\n",
    "Using [msprime](https://msprime.readthedocs.io/en/stable/) and our [VCFTools.jl](https://github.com/OpenMendel/VCFTools.jl) package, the `simulate_data.ipynb` notebook simulated data for this demonstration:\n",
    "\n",
    "+ **`ref.excludeTarget.vcf.gz`**: haplotype reference panel with 2400 samples on 36063 SNPs.\n",
    "+ **`target.typedOnly.masked.vcf.gz`**: unphased target genotypes with 100 samples, 5000 SNPs, and 1% genotypes missing at random.\n",
    "+ **`target.full.vcf.gz`**: complete genotype information (for checking imputation accuracy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imputation and phasing\n",
    "\n",
    "Now we illustrate how to use `MendelImpute.jl` for phasing and imputation. \n",
    "\n",
    "## Step 1: Generate `.jlso` compressed reference panel\n",
    "\n",
    "MendelImpute requires one to pre-process the reference panel for faster reading. This is achieved via the [compress_haplotypes](https://openmendel.github.io/MendelImpute.jl/dev/man/api/#MendelImpute.compress_haplotypes) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32mimporting reference data...100%|████████████████████████| Time: 0:00:13\u001b[39m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 27.782844 seconds (212.40 M allocations: 14.952 GiB, 9.69% gc time)\n"
     ]
    }
   ],
   "source": [
    "reffile = \"./data/ref.excludeTarget.vcf.gz\"\n",
    "tgtfile = \"./data/target.typedOnly.masked.vcf.gz\"\n",
    "outfile = \"./data/ref.excludeTarget.jlso\"\n",
    "@time compress_haplotypes(reffile, tgtfile, outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Run MendelImpute!\n",
    "\n",
    "Note: First run will be slower. Run twice for more accurate timing result. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of threads = 1\n",
      "Importing reference haplotype data...\n",
      "Total windows = 16, averaging ~ 708 unique haplotypes per window.\n",
      "\n",
      "Timings: \n",
      "    Data import                     = 0.585041 seconds\n",
      "        import target data             = 0.384378 seconds\n",
      "        import compressed haplotypes   = 0.200663 seconds\n",
      "    Computing haplotype pair        = 0.298163 seconds\n",
      "        BLAS3 mul! to get M and N      = 0.0282044 seconds per thread\n",
      "        haplopair search               = 0.250691 seconds per thread\n",
      "        initializing missing           = 0.00429557 seconds per thread\n",
      "        allocating and viewing         = 0.0147918 seconds per thread\n",
      "        index conversion               = 5.539e-5 seconds per thread\n",
      "    Phasing by win-win intersection = 0.0617189 seconds\n",
      "        Window-by-window intersection  = 0.00103244 seconds per thread\n",
      "        Breakpoint search              = 0.0587399 seconds per thread\n",
      "        Recording result               = 0.0017482 seconds per thread\n",
      "    Imputation                     = 0.318495 seconds\n",
      "        Imputing missing               = 0.017883 seconds\n",
      "        Writing to file                = 0.300612 seconds\n",
      "\n",
      "    Total time                      = 1.26407 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "tgtfile = \"./data/target.typedOnly.masked.vcf.gz\"\n",
    "reffile = \"./data/ref.excludeTarget.jlso\"\n",
    "outfile = \"./data/imputed.vcf.gz\"\n",
    "phase(tgtfile, reffile, outfile);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-analysis: Import, manipulate, and visualize data\n",
    "\n",
    "**Import:** OpenMendel's `VCFTools.jl` allow us to import VCF files directly as numeric matrices. Thus it supports all standard matrix operations and can be directly plugged into any Julia package that accepts matrices as inputs. For simplicity, we will calculate the error rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error rate = 0.0007129190583146162\n"
     ]
    }
   ],
   "source": [
    "# import genotypes as double-precision matrices\n",
    "Ximputed = convert_gt(Float64, \"./data/imputed.vcf.gz\")\n",
    "Xtrue = convert_gt(Float64, \"./data/target.full.vcf.gz\")\n",
    "err = sum(Xtrue .!= Ximputed) / size(Xtrue, 1) / size(Xtrue, 2)\n",
    "println(\"error rate = $err\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Manipulate:** Genotype matrices permit all standard matrix operations (BLAS, LAPACK...etc), so users gain access to statistical analysis *immediately* after imputation. \n",
    "\n",
    "Here, we multiply the first 10 rows and 10 columns of the imputed matrix by a random vector, getting another random vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10-element Array{Union{Missing, Float64},1}:\n",
       " 1.0201571201911694\n",
       " 0.9522113232558804\n",
       " 1.6069507335118702\n",
       " 1.6069507335118702\n",
       " 1.4507189839267838\n",
       " 1.2795810283838753\n",
       " 1.13084507188393\n",
       " 2.037512597247485\n",
       " 1.6069507335118702\n",
       " 1.6069507335118702"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ximputed[1:10, 1:10] * rand(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Visualize**: Julia is shipped with a rich plotting ecosystem. One application would be to see if there are specific regions of the chromosomes that seem to have more errors than other regions. Below we plot the error distribution for the first 200 SNPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[1m                           Visualized imputation error\u001b[22m\n",
       "\u001b[90m       ┌───────────────────────────────────────────────────────────────────┐\u001b[39m    \n",
       "     \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[31m⠤\u001b[39m\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠤\u001b[39m\u001b[31m⠠\u001b[39m\u001b[31m⠠\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠁\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠂\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠠\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠈\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠠\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[31m⠒\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠐\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠐\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠈\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠂\u001b[39m\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠠\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠄\u001b[39m\u001b[0m⠀\u001b[31m⠤\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[31m⠠\u001b[39m\u001b[31m⠄\u001b[39m\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠐\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[31m⠂\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "      \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[31m⠁\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "   \u001b[90m100\u001b[39m\u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m    \n",
       "\u001b[90m       └───────────────────────────────────────────────────────────────────┘\u001b[39m    \n",
       "\u001b[90m       1\u001b[39m\u001b[90m                                 \u001b[39m\u001b[90m                                200\u001b[39m\n",
       "\u001b[0m                                     nz = 41"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disagreeing_entries = sparse(Xtrue .!= Ximputed)\n",
    "spy(disagreeing_entries[:, 1:200], title=\"Visualized imputation error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure, each row is a sample genotype, each column is a SNP, and each dot represent an imputation error. Thus, in the top left panel, the line of dots (in sample 28) suggest the given sample is not imputed well in that region. \n",
    "\n",
    "**Post-imputation quality control**: Naturally, we seek a post-imputation quality score for each SNP (or each sample). This is not demonstrated in our talk, but [our documentation](https://openmendel.github.io/MendelImpute.jl/dev/man/Phasing_and_Imputation/#Post-imputation:-per-SNP-Imputation-Quality-Score) details this process. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion \n",
    "\n",
    "+ OpenMendel offers a fast and intuitive pipeline for genotype imputation and phasing\n",
    "+ After imputation, your data is immediately ready for statistical analysis and visualization\n",
    "\n",
    "# Other package features\n",
    "\n",
    "+ Built-in support for imputing genotypes stored in VCF files (.vcf, .vcf.gz) or PLINK files.\n",
    "+ Out-of-the-box multithreaded (shared memory) parallelism.\n",
    "+ Admixture estimation, with code examples to make pretty plots!\n",
    "+ Ultra-compressed file for phased genotypes.\n",
    "+ Imputation on dosage data\n",
    "\n",
    "**All code is on GitHub:** https://github.com/OpenMendel/MendelImpute.jl\n",
    "\n",
    "**Latest documentation:** https://openmendel.github.io/MendelImpute.jl/dev/\n",
    "\n",
    "**Our paper is going live on bioarxiv within a week.** Once it is up we will update our github page, please stay tuned!"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.0",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}