SnpArrays/SnpArraysTutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SnpArrays.jl\n",
    "\n",
    "Data from [*genome-wide association studies (GWAS)*](https://en.wikipedia.org/wiki/Genome-wide_association_study) are often saved as a [**PLINK binary biallelic genotype table**](https://www.cog-genomics.org/plink2/formats#bed) or `.bed` file. To be useful, such files should be accompanied by a `.fam` file, containing metadata on the rows of the table, and a `.bim` file,\n",
    "containing metadata on the columns. The `.fam` and `.bim` files are in tab-separated format.\n",
    "\n",
    "The table contains the observed allelic type at `n` [*single nucleotide polymorphism*](https://en.wikipedia.org/wiki/Single-nucleotide_polymorphism) (SNP) positions for `m` individuals. A SNP corresponds to a nucleotide position on the genome where some degree of variation has been observed in a population, with each individual have one of two possible *alleles* at that position on each of a pair of chromosomes. Three possible genotypes and corresponding coding are\n",
    "\n",
    "| Genotype | Plink/SnpArray |  \n",
    "|:---:|:---:|  \n",
    "| A1,A1 | 0x00 |  \n",
    "| missing | 0x01 |\n",
    "| A1,A2 | 0x02 |  \n",
    "| A2,A2 | 0x03 |  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "This package requires Julia v0.7 or later, which can be obtained from\n",
    "<https://julialang.org/downloads/> or by building Julia from the sources in the\n",
    "<https://github.com/JuliaLang/julia> repository.\n",
    "\n",
    "The package has not yet been registered and must be installed using the repository location.\n",
    "Start julia and use the `]` key to switch to the package manager REPL\n",
    "```julia\n",
    "(v1.0) pkg> add https://github.com/OpenMendel/SnpArrays.jl.git\n",
    "```\n",
    "Use the backspace key to return to the Julia REPL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 1.1.0\n",
      "Commit 80516ca202 (2019-01-21 21:24 UTC)\n",
      "Platform Info:\n",
      "  OS: macOS (x86_64-apple-darwin14.5.0)\n",
      "  CPU: Intel(R) Core(TM) i7-6567U CPU @ 3.30GHz\n",
      "  WORD_SIZE: 64\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n",
      "Environment:\n",
      "  JULIA_NUM_THREADS = 8\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SnpArrays, BenchmarkTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example data\n",
    "\n",
    "An example data set, from a study published in 2006, is about 5 Mb in size. Data from recent studies, which have samples from tens of thousands of individuals at over a million SNP positions, would be in the tens or even hundreds of Gb range.\n",
    "\n",
    "Example data location can be retrieved by using the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filename = SnpArrays.datadir(\"mouse\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r--  1 janets  staff  4922753 Feb  8 00:25 /Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.bed\n",
      "-rw-r--r--  1 janets  staff   306000 Feb  8 00:25 /Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.bim\n",
      "-rw-r--r--  1 janets  staff    92060 Feb  8 00:25 /Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.fam\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Process(`\u001b[4mls\u001b[24m \u001b[4m-l\u001b[24m \u001b[4m/Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.bed\u001b[24m \u001b[4m/Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.bim\u001b[24m \u001b[4m/Users/janets/.julia/packages/SnpArrays/c2tz4/src/../data/mouse.fam\u001b[24m`, ProcessExited(0))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run(`ls -l $filename.bed $filename.bim $filename.fam`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructor\n",
    "\n",
    "There are various ways to initialize a SnpArray.\n",
    "\n",
    "### Intitialize from Plink files\n",
    "\n",
    "All three files (`.bed`, `.fam`, `.bim`) need to be present."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×10150 SnpArray:\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x02  0x02  0x02  0x02  0x02\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x02  0x02  0x02  0x02  0x02\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03  …  0x00  0x00  0x00  0x00  0x00  0x00\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x00  0x00  0x00  0x00  0x00  0x00\n",
       "    ⋮                             ⋮  ⋱           ⋮                        \n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x03  0x02  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02     0x01  0x01  0x01  0x01  0x01  0x01\n",
       " 0x00  0x00  0x00  0x00  0x03  0x00     0x03  0x03  0x03  0x03  0x03  0x03"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const mouse = SnpArray(SnpArrays.datadir(\"mouse.bed\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The virtual size of the GWAS data is 1940 observations at each of 10150 SNP positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1940, 10150)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the file is memory-mapped opening the file and accessing the data is fast, even for very large .bed files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  114.169 μs (79 allocations: 390.73 KiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(SnpArray(SnpArrays.datadir(\"mouse.bed\")));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the memory-mapped file is read only, changing entries is not allowed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ReadOnlyMemoryError",
     "evalue": "ReadOnlyMemoryError()",
     "output_type": "error",
     "traceback": [
      "ReadOnlyMemoryError()",
      "",
      "Stacktrace:",
      " [1] | at ./int.jl:320 [inlined]",
      " [2] setindex!(::SnpArray, ::UInt8, ::Int64, ::Int64) at /Users/janets/.julia/packages/SnpArrays/c2tz4/src/snparray.jl:109",
      " [3] top-level scope at In[9]:1"
     ]
    }
   ],
   "source": [
    "mouse[1, 1] = 0x00"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To possibly change genoytpes in a bed file, open with write permission\n",
    "```julia\n",
    "mouse = SnpArray(SnpArrays.datadir(\"mouse.bed\"), \"w\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize and create bed file\n",
    "\n",
    "Initialize 5 rows and 3 columns with all (A1, A1) genotype (0x00) and memory-map to a bed file `tmp.bed`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpbf = SnpArray(\"tmp.bed\", 5, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change entries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x03  0x03  0x00\n",
       " 0x03  0x03  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpbf[1:2, 1:2] .= 0x03\n",
    "tmpbf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x02  0x02  0x02\n",
       " 0x02  0x02  0x02\n",
       " 0x02  0x02  0x02\n",
       " 0x02  0x02  0x02\n",
       " 0x02  0x02  0x02"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fill!(tmpbf, 0x02)\n",
    "tmpbf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(\"tmp.bed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize 5 rows and 3 columns with undefined genotypes without memory-mapping to any file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x01  0x00  0x01\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpbf = SnpArray(undef, 5, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a bed file corresponding to an existing SnpArray and memory-map it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x01  0x00  0x01\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpbf = SnpArray(\"tmp.bed\", tmpbf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 SnpArray:\n",
       " 0x02  0x00  0x01\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00\n",
       " 0x00  0x00  0x00"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmpbf[1, 1] = 0x02\n",
    "tmpbf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(\"tmp.bed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `convert` and `copyto!`\n",
    "\n",
    "Most common usage of SnpArray is to convert genotypes to numeric values for statistical analysis. The conversion rule depends on genetic models (additive, dominant, or recessive), centering, scaling, or imputation.\n",
    "\n",
    "### `convert`\n",
    "\n",
    "`convert` function has 4 keyword arguments: `model`, `center`, `scale`, and `impute`.\n",
    "\n",
    "`model` keyword specifies the SNP model for conversion. By default `convert` function translates genotypes according to the *additive* SNP model, which essentially counts the number of **A2** allele (0, 1 or 2) per genotype. Other SNP models are *dominant* and *recessive*, both in terms of the **A2** allele.\n",
    "\n",
    "| Genotype | `SnpArray` | `model=ADDITIVE_MODEL` | `model=DOMINANT_MODEL` | `model=RECESSIVE_MODEL` |    \n",
    "|:---:|:---:|:---:|:---:|:---:|  \n",
    "| A1,A1 | 0x00 | 0 | 0 | 0 |  \n",
    "| missing | 0x01 | NaN | NaN | NaN |\n",
    "| A1,A2 | 0x02 | 1 | 1 | 0 |  \n",
    "| A2,A2 | 0x03 | 2 | 1 | 1 |  \n",
    "\n",
    "`center=true` tells `convert` to center each column by its mean. Default is `false`.\n",
    "\n",
    "`scale=true` tells `convert` to scale each column by its standard deviation. Default is `false`.\n",
    "\n",
    "`impute=true` tells `convert` to impute missing genotypes (0x01) by column mean. Default is `false`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the whole SnpArray to a Float64 matrix using defaults (`model=ADDITIVE_MODEL`, `center=false`, `scale=false`, `impute=false`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×10150 Array{Float64,2}:\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0  …    2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  2.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       1.0    1.0    1.0    1.0    1.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0  …    2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  2.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  2.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       1.0    1.0    1.0    1.0    1.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  …    0.0    0.0    0.0    0.0    0.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       0.0    0.0    0.0    0.0    0.0\n",
       " ⋮                        ⋮              ⋱    ⋮                              \n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  …    2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0       2.0    2.0    2.0    2.0    2.0\n",
       " 2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  2.0  1.0  2.0  1.0  …    2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  2.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  2.0  1.0  1.0  1.0  1.0  2.0       2.0    2.0    2.0    2.0    2.0\n",
       " 1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0     NaN    NaN    NaN    NaN    NaN  \n",
       " 0.0  0.0  0.0  0.0  2.0  0.0  2.0  0.0       2.0    2.0    2.0    2.0    2.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Matrix{Float64}, mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "!!! note  \n",
    "\n",
    "    When `convert` or `copyto!` is used a slice or subarray of SnpArray, using `view`, `@view` or `views` is necessary for both correctness and efficiency. Without view, it's simply converting the UInt8 coding in original bed file.\n",
    "    \n",
    "!!!\n",
    "\n",
    "Convert a column to Float64 vector using defaults (`model=ADDITIVE_MODEL`, `center=false`, `scale=false`, `impute=false`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " ⋮  \n",
       " 2.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 0.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# convert(Vector{Float64}, view(mouse, :, 1)) # alternative syntax\n",
    "# @views convert(Vector{Float64}, mouse[:, 1]) # alternative syntax\n",
    "convert(Vector{Float64}, @view(mouse[:, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert a subarray of SnpArray to a Float64 matrix using defaults (`model=ADDITIVE_MODEL`, `center=false`, `scale=false`, `impute=false`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 1.0  1.0  2.0  2.0  1.0\n",
       " 2.0  2.0  2.0  2.0  2.0\n",
       " 2.0  2.0  2.0  2.0  2.0\n",
       " 1.0  1.0  2.0  2.0  1.0\n",
       " 1.0  2.0  1.0  1.0  1.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Matrix{Float64}, @view(mouse[1:2:10, 1:2:10]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a matrix of different SNP models (`ADDITIVE_MODEL` vs `DOMINANT_MODEL` vs `RECESSIVE_MODEL`) for the same snp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×3 Array{Float64,2}:\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " ⋮            \n",
       " 2.0  1.0  1.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 2.0  1.0  1.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 1.0  1.0  0.0\n",
       " 0.0  0.0  0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@views [convert(Vector{Float64}, mouse[:, 1], model=ADDITIVE_MODEL) convert(Vector{Float64}, mouse[:, 1], model=DOMINANT_MODEL) convert(Vector{Float64}, mouse[:, 1], model=RECESSIVE_MODEL)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center and scale (last column) while uaing `convert`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "  -1.8819155626127624\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "  -1.8819155626127624\n",
       "  -4.2359771983402785\n",
       "   0.4721460731147541\n",
       "  -4.2359771983402785\n",
       "   ⋮                 \n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       "   0.4721460731147541\n",
       " NaN                 \n",
       "   0.4721460731147541"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Vector{Float64}, @view(mouse[:, end]), center=true, scale=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center, scale, and impute (last column) while using `convert`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       " -1.8819155626127624\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       " -1.8819155626127624\n",
       " -4.2359771983402785\n",
       "  0.4721460731147541\n",
       " -4.2359771983402785\n",
       "  ⋮                 \n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.4721460731147541\n",
       "  0.0               \n",
       "  0.4721460731147541"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Vector{Float64}, @view(mouse[:, end]), center=true, scale=true, impute=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `copyto!`\n",
    "\n",
    "`copyto!` is the in-place version of `convert`. It takes the same keyword arguments (`model`, `center`, `scale`, `impute`) as `convert`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy a column to a Float64 vector using defaults (`model=:additive`, `center=false`, `scale=false`, `impute=false`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " ⋮  \n",
       " 2.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 2.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0\n",
       " 0.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = zeros(size(mouse, 1))\n",
    "copyto!(v, @view(mouse[:, 1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.297 μs (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime(copyto!($v, $@view(mouse[:, 1])));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy columns using defaults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×2 Array{Float64,2}:\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " ⋮       \n",
       " 2.0  2.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 2.0  2.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 1.0  1.0\n",
       " 0.0  0.0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v2 = zeros(size(mouse, 1), 2)\n",
    "copyto!(v2, @view(mouse[:, 1:2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  11.083 μs (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "# roughly double the cost of copying 1 column\n",
    "@btime(copyto!($v2, $@view(mouse[:, 1:2])));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center and scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       "  ⋮                  \n",
       "  1.2624897581484626 \n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       "  1.2624897581484626 \n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       " -0.16084075452851265\n",
       " -1.584171267205488  "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "copyto!(v, @view(mouse[:, 1]), center=true, scale=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.926 μs (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "# more cost becoz of extra pass for center, scale, and/or impute\n",
    "@btime(copyto!($v, $(@view(mouse[:, 1])), center=true, scale=true));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looping over all columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  42.784 ms (10150 allocations: 475.78 KiB)\n"
     ]
    }
   ],
   "source": [
    "v = Vector{Float64}(undef, size(mouse, 1))\n",
    "function loop_test(v, s)\n",
    "    for j in 1:size(s, 2)\n",
    "        copyto!(v, @view(s[:, j]))\n",
    "    end\n",
    "end\n",
    "@btime(loop_test($v, $mouse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy the whole SnpArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  43.859 ms (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "M = similar(mouse, Float64)\n",
    "@btime(copyto!($M, $mouse));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summaries\n",
    "\n",
    "### Counts\n",
    "\n",
    "Counts of each the four possible values for each column are returned by `counts`.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×10150 Array{Int64,2}:\n",
       "  358   359  252   358    33   359  …    56    56    56    56    56    56\n",
       "    2     0    4     3     4     1      173   173   162   173   174   175\n",
       " 1003  1004  888  1004   442  1004      242   242   242   242   242   242\n",
       "  577   577  796   575  1461   576     1469  1469  1480  1469  1468  1467"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts(mouse, dims=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Column 2 has no missing values (code `0x01`, the second row in the column-counts table).\n",
    "In that SNP position for this sample, 359 individuals are homozygous allele 1 (`G` according to the `.bim` file), 1004 are heterozygous, and 577 are homozygous allele 2 (`A`).\n",
    "\n",
    "The counts by column and by row are cached in the `SnpArray` object. Accesses after the first are extremely fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.845 ns (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime(counts($mouse, dims=1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minor allele frequencies\n",
    "\n",
    "Minor allele frequencies (MAF) for each SNP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10150-element Array{Float64,1}:\n",
       " 0.4434984520123839  \n",
       " 0.4438144329896907  \n",
       " 0.359504132231405   \n",
       " 0.4439855446566856  \n",
       " 0.13119834710743805 \n",
       " 0.44404332129963897 \n",
       " 0.1412706611570248  \n",
       " 0.30299123259412064 \n",
       " 0.4445018069179143  \n",
       " 0.44424367578729995 \n",
       " 0.43427835051546393 \n",
       " 0.14075413223140498 \n",
       " 0.304639175257732   \n",
       " ⋮                   \n",
       " 0.0527624309392265  \n",
       " 0.052980132450331174\n",
       " 0.08079096045197742 \n",
       " 0.08253250423968339 \n",
       " 0.08253250423968339 \n",
       " 0.10022650056625138 \n",
       " 0.10016977928692694 \n",
       " 0.10016977928692694 \n",
       " 0.09955005624296964 \n",
       " 0.10016977928692694 \n",
       " 0.10022650056625138 \n",
       " 0.10028328611898019 "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maf(mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minor allele (`false` means A1 is the minor allele; `true` means A2 is the minor allele) for each SNP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10150-element BitArray{1}:\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       "     ⋮\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false\n",
       " false"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minorallele(mouse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `mean` and `var`\n",
    "\n",
    "The package provides methods for the generics `mean` and `var` from the `Statistics` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10150 Array{Float64,2}:\n",
       " 1.113  1.11237  1.28099  1.11203  …  1.8009  1.79966  1.79955  1.79943"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean(mouse, dims=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10150 Array{Float64,2}:\n",
       " 0.815273  0.814948  0.869835  0.815178  …  0.968308  0.96829  0.968272"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean(mouse, dims=1, model=DOMINANT_MODEL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10150 Array{Float64,2}:\n",
       " 0.469929  0.470089  0.462605  0.469365  …  0.223714  0.223818  0.223923"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var(mouse, dims=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These methods make use of the cached column or row counts and thus are very fast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  14.688 μs (2 allocations: 79.39 KiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(mean($mouse, dims=1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The column-wise or row-wise standard deviations are returned by `std`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×1 Array{Float64,2}:\n",
       " 0.6504997290784408\n",
       " 0.6379008244533891\n",
       " 0.6558172726141286\n",
       " 0.6532675479248437\n",
       " 0.6744432174014563\n",
       " 0.6519092298111158\n",
       " 0.6779881845456428\n",
       " 0.6955814098050999\n",
       " 0.6437566832989493\n",
       " 0.6505283141088536\n",
       " 0.665444994623426 \n",
       " 0.659392039592328 \n",
       " 0.6641674726999468\n",
       " ⋮                 \n",
       " 0.6599158250006595\n",
       " 0.688387450736178 \n",
       " 0.6664063015924304\n",
       " 0.6613451651895259\n",
       " 0.6659810347614777\n",
       " 0.6274577846909379\n",
       " 0.6823658517777204\n",
       " 0.6695299551061924\n",
       " 0.710756592739754 \n",
       " 0.6387913736114869\n",
       " 0.6736492722732016\n",
       " 0.688855476425891 "
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std(mouse, dims=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing rate\n",
    "\n",
    "Proportion of missing genotypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10150-element Array{Float64,1}:\n",
       " 0.0010309278350515464\n",
       " 0.0                  \n",
       " 0.002061855670103093 \n",
       " 0.0015463917525773195\n",
       " 0.002061855670103093 \n",
       " 0.0005154639175257732\n",
       " 0.002061855670103093 \n",
       " 0.0005154639175257732\n",
       " 0.0015463917525773195\n",
       " 0.0015463917525773195\n",
       " 0.0                  \n",
       " 0.002061855670103093 \n",
       " 0.0                  \n",
       " ⋮                    \n",
       " 0.06701030927835051  \n",
       " 0.06597938144329897  \n",
       " 0.08762886597938144  \n",
       " 0.08814432989690722  \n",
       " 0.08814432989690722  \n",
       " 0.08969072164948454  \n",
       " 0.08917525773195877  \n",
       " 0.08917525773195877  \n",
       " 0.08350515463917525  \n",
       " 0.08917525773195877  \n",
       " 0.08969072164948454  \n",
       " 0.09020618556701031  "
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "missingrate(mouse, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940-element Array{Float64,1}:\n",
       " 0.00019704433497536947\n",
       " 0.0                   \n",
       " 0.018423645320197045  \n",
       " 0.0007881773399014779 \n",
       " 0.0                   \n",
       " 0.004236453201970443  \n",
       " 0.0051231527093596055 \n",
       " 0.00039408866995073894\n",
       " 0.005517241379310344  \n",
       " 0.0016748768472906405 \n",
       " 0.0                   \n",
       " 9.852216748768474e-5  \n",
       " 0.0004926108374384236 \n",
       " ⋮                     \n",
       " 0.000689655172413793  \n",
       " 0.004729064039408867  \n",
       " 0.0004926108374384236 \n",
       " 0.001083743842364532  \n",
       " 0.00019704433497536947\n",
       " 0.0025615763546798028 \n",
       " 0.0038423645320197044 \n",
       " 0.001379310344827586  \n",
       " 0.0064039408866995075 \n",
       " 0.002857142857142857  \n",
       " 0.0011822660098522167 \n",
       " 0.00029556650246305416"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "missingrate(mouse, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Location of the missing values\n",
    "\n",
    "The positions of the missing data are evaluated by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×10150 SparseArrays.SparseMatrixCSC{Bool,Int32} with 33922 stored entries:\n",
       "  [702  ,     1]  =  true\n",
       "  [949  ,     1]  =  true\n",
       "  [914  ,     3]  =  true\n",
       "  [949  ,     3]  =  true\n",
       "  [1604 ,     3]  =  true\n",
       "  [1891 ,     3]  =  true\n",
       "  [81   ,     4]  =  true\n",
       "  [990  ,     4]  =  true\n",
       "  [1882 ,     4]  =  true\n",
       "  [81   ,     5]  =  true\n",
       "  [676  ,     5]  =  true\n",
       "  [990  ,     5]  =  true\n",
       "  ⋮\n",
       "  [1791 , 10150]  =  true\n",
       "  [1795 , 10150]  =  true\n",
       "  [1846 , 10150]  =  true\n",
       "  [1848 , 10150]  =  true\n",
       "  [1851 , 10150]  =  true\n",
       "  [1853 , 10150]  =  true\n",
       "  [1860 , 10150]  =  true\n",
       "  [1873 , 10150]  =  true\n",
       "  [1886 , 10150]  =  true\n",
       "  [1894 , 10150]  =  true\n",
       "  [1897 , 10150]  =  true\n",
       "  [1939 , 10150]  =  true"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp = missingpos(mouse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  36.420 ms (19273 allocations: 1.80 MiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(missingpos($mouse));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, for example, the number of missing data values in each column can be evaluated as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10150 Array{Int64,2}:\n",
       " 2  0  4  3  4  1  4  1  3  3  0  4  0  …  174  173  173  162  173  174  175"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(mp, dims=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a faster, but somewhat more obscure, way to calculate the same sums"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×10150 view(::Array{Int64,2}, 2:2, :) with eltype Int64:\n",
       " 2  0  4  3  4  1  4  1  3  3  0  4  0  …  174  173  173  162  173  174  175"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "view(counts(mouse, dims=1), 2:2, :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Genetic relationship matrix\n",
    "\n",
    "`grm` function computes the empirical kinship matrix using either the classical genetic relationship matrix, `grm(A, model=:GRM)`, or the method of moment method, `grm(A, model=:MoM)`, or the robust method, `grm(A, model=:Robust)`. \n",
    "\n",
    "Classical genetic relation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×1940 Array{Float64,2}:\n",
       "  0.478301    -0.0331304    0.0135612    …  -0.0347737   -0.0129443 \n",
       " -0.0331304    0.422771    -0.0389227        0.0457987    0.00556832\n",
       "  0.0135612   -0.0389227    0.509248        -0.0356689   -0.0608705 \n",
       "  0.0198205    0.00728645  -0.00935362      -0.0302404   -0.0102152 \n",
       "  0.056747    -0.0163418   -0.00495283      -0.0413347   -0.0415659 \n",
       " -0.0165628   -0.0191127   -0.0112181    …   0.0177118   -0.0193087 \n",
       "  0.123771    -0.0404167    0.00442739       0.00880649  -0.0437565 \n",
       " -0.0628362    0.172552    -0.0728312        0.0640027   -0.0281429 \n",
       "  0.0605018   -0.0260505    0.00398852      -0.00277754  -0.0607773 \n",
       "  0.108886    -0.0204594   -0.00767711      -0.0210501    0.00343526\n",
       " -0.0142307    0.00270989  -0.0235504    …  -0.0223563   -0.028408  \n",
       " -0.0306022    0.197743    -0.00244269       0.0213998   -0.0478472 \n",
       " -0.0131463   -0.0226707    0.0223522       -0.037288     0.0493662 \n",
       "  ⋮                                      ⋱                          \n",
       "  0.0176725   -0.0165609    0.0378308        0.0238751   -0.0420143 \n",
       "  0.0024949   -0.0411137    0.0154847       -0.0380656   -0.0650806 \n",
       "  0.0952286    0.00894298  -0.0163446    …  -0.0202633   -0.0219594 \n",
       " -0.0309488   -0.0228342   -0.0478253       -0.014896     0.261623  \n",
       " -0.004804    -0.0375168   -0.0211418       -0.0172572    0.0359166 \n",
       "  0.0076296    0.0481887   -0.0328968        0.0920425   -0.0292548 \n",
       "  0.070045    -0.0302138    0.000647283      0.00892069  -0.00632566\n",
       "  0.0378132   -6.59565e-5   0.00888932   …   0.00230815  -0.0291622 \n",
       " -0.00132837   0.00223654   0.0495928       -0.00936248   0.0299075 \n",
       "  0.0640864   -0.0241218    0.00602283       0.00403413   0.00689551\n",
       " -0.0347737    0.0457987   -0.0356689        0.509228    -0.035215  \n",
       " -0.0129443    0.00556832  -0.0608705       -0.035215     0.552712  "
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# grm(mouse, method=:MoM)\n",
    "# grm(mouse, method=:Robust)\n",
    "grm(mouse, method=:GRM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  664.864 ms (30 allocations: 28.95 MiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(grm($mouse, method=:GRM));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Float32 (single precision) potentially saves memory usage and computation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×1940 Array{Float32,2}:\n",
       "  0.478301    -0.0331304    0.0135612    …  -0.0347737   -0.0129443 \n",
       " -0.0331304    0.422771    -0.0389227        0.0457987    0.00556833\n",
       "  0.0135612   -0.0389227    0.509248        -0.0356689   -0.0608705 \n",
       "  0.0198205    0.00728645  -0.00935361      -0.0302404   -0.0102152 \n",
       "  0.056747    -0.0163418   -0.00495284      -0.0413347   -0.0415659 \n",
       " -0.0165628   -0.0191127   -0.0112181    …   0.0177117   -0.0193087 \n",
       "  0.123771    -0.0404167    0.0044274        0.0088065   -0.0437565 \n",
       " -0.0628363    0.172552    -0.0728312        0.0640027   -0.0281429 \n",
       "  0.0605018   -0.0260505    0.00398853      -0.00277754  -0.0607773 \n",
       "  0.108886    -0.0204594   -0.00767711      -0.0210501    0.00343524\n",
       " -0.0142307    0.00270989  -0.0235504    …  -0.0223563   -0.028408  \n",
       " -0.0306022    0.197743    -0.00244268       0.0213998   -0.0478472 \n",
       " -0.0131464   -0.0226707    0.0223522       -0.037288     0.0493662 \n",
       "  ⋮                                      ⋱                          \n",
       "  0.0176725   -0.016561     0.0378308        0.0238751   -0.0420143 \n",
       "  0.00249491  -0.0411137    0.0154847       -0.0380656   -0.0650806 \n",
       "  0.0952286    0.00894298  -0.0163446    …  -0.0202633   -0.0219594 \n",
       " -0.0309488   -0.0228342   -0.0478253       -0.014896     0.261623  \n",
       " -0.00480401  -0.0375167   -0.0211418       -0.0172572    0.0359166 \n",
       "  0.00762961   0.0481887   -0.0328968        0.0920425   -0.0292547 \n",
       "  0.070045    -0.0302138    0.000647269      0.00892068  -0.00632566\n",
       "  0.0378132   -6.59475e-5   0.00888932   …   0.00230815  -0.0291622 \n",
       " -0.00132838   0.00223653   0.0495928       -0.00936246   0.0299075 \n",
       "  0.0640864   -0.0241219    0.00602283       0.00403413   0.00689551\n",
       " -0.0347737    0.0457987   -0.0356689        0.509228    -0.035215  \n",
       " -0.0129443    0.00556833  -0.0608705       -0.035215     0.552712  "
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grm(mouse, method=:GRM, t=Float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  369.065 ms (31 allocations: 14.60 MiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(grm($mouse, method=:GRM, t=Float32));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, `grm` excludes SNPs with minor allele frequency below 0.01. This default can be changed by the keyword argument `minmaf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×1940 Array{Float64,2}:\n",
       "  0.478556    -0.0331783    0.013541     …  -0.0348225   -0.0129761 \n",
       " -0.0331783    0.422993    -0.0389741        0.0457975    0.00554753\n",
       "  0.013541    -0.0389741    0.50952         -0.0357183   -0.0609305 \n",
       "  0.0203209    0.00777944  -0.00887047      -0.0297696   -0.00972836\n",
       "  0.0567523   -0.0163798   -0.00498406      -0.0413874   -0.0416146 \n",
       " -0.0166009   -0.0191523   -0.0112531    …   0.0176939   -0.0193442 \n",
       "  0.123816    -0.0404689    0.00440171       0.0087834   -0.0438065 \n",
       " -0.0629017    0.172626    -0.0729026        0.0640123   -0.0281836 \n",
       "  0.0605093   -0.0260942    0.00396257      -0.00280748  -0.0608373 \n",
       "  0.108922    -0.0204998   -0.00770996      -0.0210909    0.00341321\n",
       " -0.0142674    0.00268319  -0.0235927    …  -0.0223978   -0.0284489 \n",
       " -0.0306486    0.197832    -0.00247243       0.0213842   -0.0478996 \n",
       " -0.0131824   -0.0227124    0.0223371       -0.0373384    0.0493713 \n",
       "  ⋮                                      ⋱                          \n",
       "  0.0176546   -0.016599     0.0378249        0.0238609   -0.0420633 \n",
       "  0.00246808  -0.0411663    0.0154656       -0.0381165   -0.0651432 \n",
       "  0.0952566    0.00891997  -0.0163826    …  -0.0203036   -0.0219965 \n",
       " -0.0309912   -0.0228718   -0.0478777       -0.0149289    0.261754  \n",
       " -0.00483514  -0.0375673   -0.0211827       -0.0172957    0.0359138 \n",
       "  0.00770862   0.0482917   -0.0328417        0.0921714   -0.0292961 \n",
       "  0.0700582   -0.03026      0.000619365      0.00889767  -0.00635348\n",
       "  0.0378313   -7.02155e-5   0.00889036   …   0.0023053   -0.0291795 \n",
       " -0.00133338   0.00223364   0.0496179       -0.00937223   0.0299252 \n",
       "  0.0641201   -0.0241403    0.00602217       0.0040323    0.00689958\n",
       " -0.0348225    0.0457975   -0.0357183        0.509501    -0.0352599 \n",
       " -0.0129761    0.00554753  -0.0609305       -0.0352599    0.553015  "
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute GRM excluding SNPs with MAF≤0.05 \n",
    "grm(mouse, minmaf=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To specify specific SNPs for calculating the empirical kinship, use the `cinds` keyword (default is `nothing`). When `cinds` is specified, `minmaf` is ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1940×1940 Array{Float64,2}:\n",
       "  0.477       -0.0307774     0.0118026   …  -0.0320301    -0.0125113 \n",
       " -0.0307774    0.425085     -0.0367459       0.0480442     0.00519065\n",
       "  0.0118026   -0.0367459     0.505038       -0.0385129    -0.0631557 \n",
       "  0.0166017    0.00614789   -0.00919695     -0.0399744    -0.0104884 \n",
       "  0.05724     -0.0122148    -0.00543377     -0.0395663    -0.0372998 \n",
       " -0.0193129   -0.0224378    -0.009277    …   0.0153785    -0.0220184 \n",
       "  0.12194     -0.0410682     0.00274307      0.00796748   -0.0441578 \n",
       " -0.0624031    0.173985     -0.0724784       0.0663191    -0.0294243 \n",
       "  0.0627626   -0.0288615     0.00265615     -0.00449877   -0.0579702 \n",
       "  0.110878    -0.0232715    -0.00881604     -0.021272      0.00169016\n",
       " -0.00800735  -0.00149824   -0.019791    …  -0.024124     -0.0289397 \n",
       " -0.0272944    0.19894      -0.00534771      0.0209384    -0.0511051 \n",
       " -0.011388    -0.0281003     0.0273853      -0.0360047     0.0459359 \n",
       "  ⋮                                      ⋱                           \n",
       "  0.0169431   -0.0136989     0.0340794       0.0272811    -0.041189  \n",
       "  0.00201325  -0.0426611     0.0124353      -0.0387982    -0.0656181 \n",
       "  0.097587     0.0058123    -0.0160698   …  -0.021457     -0.023226  \n",
       " -0.0342014   -0.0211246    -0.0490112      -0.0129575     0.256552  \n",
       " -0.00324255  -0.0423482    -0.0192699      -0.0149015     0.0339388 \n",
       "  0.00575353   0.0464237    -0.0294694       0.0924759    -0.0275451 \n",
       "  0.0748725   -0.0258461    -0.00141068      0.0115232    -0.00486589\n",
       "  0.0386555    0.000612169   0.00959997  …  -0.000357284  -0.0334687 \n",
       " -0.00343056   0.0120673     0.0455375      -0.0103798     0.0336959 \n",
       "  0.0656909   -0.0193469     0.00600815      0.00188545    0.00726181\n",
       " -0.0320301    0.0480442    -0.0385129       0.513285     -0.0317963 \n",
       " -0.0125113    0.00519065   -0.0631557      -0.0317963     0.54471   "
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# GRM using every other SNP\n",
    "grm(mouse, cinds=1:2:size(mouse, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitering by missing rate\n",
    "\n",
    "Count the number of rows and columns that have proportion of missingness < 0.001."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rowmr = count(missingrate(mouse, 2) .< 0.001) = 1157\n",
      "colmr = count(missingrate(mouse, 1) .< 0.001) = 5726\n"
     ]
    }
   ],
   "source": [
    "@show rowmr = count(missingrate(mouse, 2) .< 0.001)\n",
    "@show colmr = count(missingrate(mouse, 1) .< 0.001);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter according to minimum success rates (1 - proportion of missing genotypes) per row and column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Bool[true, true, false, true, true, false, false, true, false, false  …  true, false, true, false, false, false, false, false, false, true], Bool[false, true, false, false, false, true, false, true, false, false  …  false, false, false, false, false, false, false, false, false, false])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rowmask, colmask =  SnpArrays.filter(mouse, 0.999, 0.999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1157, 5726)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count(rowmask), count(colmask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  90.812 ms (8 allocations: 96.45 KiB)\n"
     ]
    }
   ],
   "source": [
    "@btime(SnpArrays.filter($mouse, 0.999, 0.999));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One may use the `rowmask` and `colmask` to filter and save the filtering results as Plink files.\n",
    "```julia\n",
    "SnpArrays.filter(SnpArrays.datadir(\"mouse\"), rowmask, colmask)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter Plink files\n",
    "\n",
    "Filter a set of Plink files according to row indices and column indices. By result, filtered Plink files are saved as `srcname.filtered.bed`, `srcname.filtered.fam`, and `srcname.filtered.bim`, where `srcname` is the source Plink file name. You can also specify the destination file name using keyword `des`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 SnpArray:\n",
       " 0x02  0x02  0x02  0x02  0x03\n",
       " 0x02  0x02  0x03  0x02  0x02\n",
       " 0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x02\n",
       " 0x03  0x03  0x03  0x03  0x03"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SnpArrays.filter(SnpArrays.datadir(\"mouse\"), 1:5, 1:5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.bed\")\n",
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.fam\")\n",
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.bim\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter a set of Plink files according to logical vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1157×5726 SnpArray:\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02  …  0x03  0x02  0x00  0x02  0x03  0x02\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x00  0x03  0x03  0x03  0x03  0x00\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x03  0x03  0x00  0x03  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x03  0x03  0x02  0x03  0x02\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x00  0x03  0x03  0x00  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03  …  0x03  0x02  0x00  0x03  0x03  0x03\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02     0x00  0x03  0x03  0x03  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x00  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x03  0x03  0x03  0x03  0x00\n",
       " 0x00  0x00  0x02  0x02  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03  …  0x03  0x02  0x02  0x03  0x03  0x00\n",
       " 0x02  0x02  0x02  0x02  0x02  0x02     0x00  0x03  0x03  0x02  0x03  0x00\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x02  0x03  0x02  0x03  0x03  0x02\n",
       "    ⋮                             ⋮  ⋱     ⋮                             ⋮\n",
       " 0x00  0x00  0x02  0x00  0x02  0x02  …  0x02  0x03  0x02  0x03  0x02  0x02\n",
       " 0x00  0x00  0x02  0x02  0x02  0x02     0x02  0x03  0x03  0x02  0x03  0x03\n",
       " 0x00  0x00  0x00  0x00  0x00  0x00     0x02  0x03  0x03  0x03  0x00  0x00\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x00  0x00  0x02  0x02  0x02  0x02     0x00  0x03  0x03  0x02  0x02  0x03\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03  …  0x02  0x03  0x03  0x03  0x03  0x02\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x00  0x03  0x03  0x02  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x03  0x02  0x03  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x02  0x02  0x03  0x03  0x03\n",
       " 0x02  0x02  0x03  0x02  0x03  0x03     0x00  0x03  0x03  0x03  0x03  0x00\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03  …  0x00  0x03  0x03  0x03  0x03  0x03\n",
       " 0x00  0x00  0x00  0x00  0x00  0x00     0x03  0x03  0x03  0x03  0x03  0x03"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SnpArrays.filter(SnpArrays.datadir(\"mouse\"), rowmask, colmask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ls: ../data/mouse.filtered.bed: No such file or directory\n",
      "ls: ../data/mouse.filtered.bim: No such file or directory\n",
      "ls: ../data/mouse.filtered.fam: No such file or directory\n"
     ]
    }
   ],
   "source": [
    ";ls -l ../data/mouse.filtered.bed ../data/mouse.filtered.fam ../data/mouse.filtered.bim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.bed\")\n",
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.fam\")\n",
    "rm(SnpArrays.datadir(\"mouse\") * \".filtered.bim\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear algebra with SnpArray\n",
    "\n",
    "In some applications we want to perform linear algebra using SnpArray directly without expanding it to numeric matrix. This is achieved by the `SnpBitMatrix` type. The implementation assumes:\n",
    "\n",
    "1. the SnpArray does not have missing genotypes, and\n",
    "2. the matrix corresponding to SnpArray is the matrix of A2 allele counts.\n",
    "\n",
    "First let's load a data set without missing genotypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "379×54051 SnpArray:\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x02  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x02  0x02  0x03  0x03  0x02\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x02  0x02  0x02  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03  …  0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x02  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03  …  0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x03  0x02  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       "    ⋮                             ⋮  ⋱     ⋮                             ⋮\n",
       " 0x03  0x03  0x03  0x00  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x02  0x02  0x02  0x03  0x02  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x03  0x03  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x00  0x03     0x02  0x02  0x02  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03  …  0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x00  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const EUR = SnpArray(SnpArrays.datadir(\"EUR_subset.bed\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To instantiate a SnpBitMatrix based on SnpArray,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "const EURbm = SnpBitMatrix{Float64}(EUR, model=ADDITIVE_MODEL, center=true, scale=true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constructor shares the same keyword arguments as the `convert` or `copyto!` functions. The type parameter, `Float64` in this example, indicates the SnpBitMatrix acts like a Float64 matrix.\n",
    "\n",
    "The memory usage of the SnpBitMatrix should be similar to the SnpArray, or equivalently bed file size, if `model=ADDITIVE_MODEL` or half of that of SnpArray if `model=DOMINANT_MODEL` or `model=RECESSIVE_MODEL`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6876757, 6421960)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Base.summarysize(EUR), Base.summarysize(EURbm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A SnpBitMatrix acts similar to a regular matrix and responds to `size`, `eltype`, and SnpBitMatrix-vector multiplications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size(EURbm) = (379, 54051)\n",
      "eltype(EURbm) = Float64\n",
      "typeof(EURbm) <: AbstractMatrix = true\n"
     ]
    }
   ],
   "source": [
    "@show size(EURbm)\n",
    "@show eltype(EURbm)\n",
    "@show typeof(EURbm) <: AbstractMatrix;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SnpBitMatrix-vector multiplication is mathematically equivalent to the corresponding Float matrix contained from `convert` or `copyto!` a SnpArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.892068045077235e-11"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "v1 = randn(size(EUR, 1))\n",
    "v2 = randn(size(EUR, 2))\n",
    "A = convert(Matrix{Float64}, EUR, model=ADDITIVE_MODEL, center=true, scale=true)\n",
    "norm(EURbm * v2 -  A * v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.505289475808036e-12"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(EURbm' * v1 - A' * v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the Float64 matrix fits into memory so Matrix{Float64}-vector multiplication (highly optimized BLAS) is much faster than the SnpBitMatrix-vector multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  85.600 ms (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "out1 = Vector{Float64}(undef, size(EUR, 1))\n",
    "out2 = Vector{Float64}(undef, size(EUR, 2))\n",
    "@btime(mul!($out1, $EURbm, $v2));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  7.319 ms (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime(mul!($out1, $A, $v2));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  87.429 ms (1 allocation: 16 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime(mul!($out2, $transpose($EURbm), $v1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.646 ms (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime(mul!($out2, $transpose($A), $v1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In another test example with a ~1GB bed file, SnpBitMatrix-vector multiplication is about 3-5 fold faster than the corresponding Matrix{Float64}-vector multiplication, because the Matrix{Float64} matrix cannot fit into the memory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`SnpBitMatrix` can be created from a subarray of SnpArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "EURsub = @view EUR[1:2:100, 1:2:100]\n",
    "EURsubbm = SnpBitMatrix{Float64}(EURsub, model=ADDITIVE_MODEL, center=true, scale=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2600"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Base.summarysize(EURsubbm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size(EURsubbm) = (50, 50)\n",
      "eltype(EURsubbm) = Float64\n",
      "typeof(EURsubbm) <: AbstractMatrix = true\n"
     ]
    }
   ],
   "source": [
    "@show size(EURsubbm)\n",
    "@show eltype(EURsubbm)\n",
    "@show typeof(EURsubbm) <: AbstractMatrix;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.18686710289806e-14"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "v1 = randn(size(EURsub, 1))\n",
    "v2 = randn(size(EURsub, 2))\n",
    "A = convert(Matrix{Float64}, EURsub, model=ADDITIVE_MODEL, center=true, scale=true)\n",
    "norm(EURsubbm * v2 -  A * v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.516016712755663e-14"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(EURsubbm' * v1 - A' * v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SnpData"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a `SnpData` dataframe, which is a `SnpArray` with information on SNP and subject appended."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SnpData(379, 54051, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x03; … ; 0x02 0x03 … 0x03 0x03; 0x03 0x03 … 0x03 0x03], 54051×6 DataFrames.DataFrame. Omitted printing of 2 columns\n",
       "│ Row   │ chromosome │ snpid       │ genetic_distance │ position │\n",
       "│       │ \u001b[90mString\u001b[39m     │ \u001b[90mString\u001b[39m      │ \u001b[90mFloat64\u001b[39m          │ \u001b[90mInt64\u001b[39m    │\n",
       "├───────┼────────────┼─────────────┼──────────────────┼──────────┤\n",
       "│ 1     │ 17         │ rs34151105  │ 0.0              │ 1665     │\n",
       "│ 2     │ 17         │ rs143500173 │ 0.0              │ 2748     │\n",
       "│ 3     │ 17         │ rs113560219 │ 0.0              │ 4702     │\n",
       "│ 4     │ 17         │ rs1882989   │ 5.6e-5           │ 15222    │\n",
       "│ 5     │ 17         │ rs8069133   │ 0.000499         │ 32311    │\n",
       "│ 6     │ 17         │ rs112221137 │ 0.000605         │ 36405    │\n",
       "│ 7     │ 17         │ rs34889101  │ 0.00062          │ 36975    │\n",
       "│ 8     │ 17         │ rs35840960  │ 0.000668         │ 38827    │\n",
       "│ 9     │ 17         │ rs144918387 │ 0.000775         │ 42965    │\n",
       "│ 10    │ 17         │ rs62057022  │ 0.000948         │ 49640    │\n",
       "⋮\n",
       "│ 54041 │ 22         │ rs6010070   │ 0.750874         │ 51180959 │\n",
       "│ 54042 │ 22         │ rs6010072   │ 0.750878         │ 51181519 │\n",
       "│ 54043 │ 22         │ rs6009960   │ 0.750879         │ 51181568 │\n",
       "│ 54044 │ 22         │ rs56807126  │ 0.750879         │ 51181624 │\n",
       "│ 54045 │ 22         │ rs6010073   │ 0.750882         │ 51181986 │\n",
       "│ 54046 │ 22         │ rs184517959 │ 0.750893         │ 51183421 │\n",
       "│ 54047 │ 22         │ rs113391741 │ 0.750923         │ 51187440 │\n",
       "│ 54048 │ 22         │ rs151247655 │ 0.750938         │ 51189403 │\n",
       "│ 54049 │ 22         │ rs187225588 │ 0.751001         │ 51197602 │\n",
       "│ 54050 │ 22         │ rs9616967   │ 0.75109          │ 51209343 │\n",
       "│ 54051 │ 22         │ rs148755559 │ 0.751156         │ 51218133 │, 379×6 DataFrames.DataFrame\n",
       "│ Row │ fid       │ iid       │ father    │ mother    │ sex       │ phenotype │\n",
       "│     │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │\n",
       "├─────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┤\n",
       "│ 1   │ 1         │ HG00096   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 2   │ 2         │ HG00097   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 3   │ 3         │ HG00099   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 4   │ 4         │ HG00100   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 5   │ 5         │ HG00101   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 6   │ 6         │ HG00102   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 7   │ 7         │ HG00103   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 8   │ 8         │ HG00104   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 9   │ 9         │ HG00106   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 10  │ 10        │ HG00108   │ 0         │ 0         │ 1         │ 1         │\n",
       "⋮\n",
       "│ 369 │ 369       │ NA20810   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 370 │ 370       │ NA20811   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 371 │ 371       │ NA20812   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 372 │ 372       │ NA20813   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 373 │ 373       │ NA20814   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 374 │ 374       │ NA20815   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 375 │ 375       │ NA20816   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 376 │ 376       │ NA20818   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 377 │ 377       │ NA20819   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 378 │ 378       │ NA20826   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 379 │ 379       │ NA20828   │ 0         │ 0         │ 2         │ 1         │)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EUR_data = SnpData(SnpArrays.datadir(\"EUR_subset\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can split `SnpData` by SNP's choromosomes using `split_plink`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{AbstractString,SnpData} with 6 entries:\n",
       "  \"21\" => SnpData(379, 5813, UInt8[0x02 0x02 … 0x03 0x03; 0x03 0x03 … 0x03 0x03…\n",
       "  \"17\" => SnpData(379, 11041, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x0…\n",
       "  \"19\" => SnpData(379, 9690, UInt8[0x02 0x02 … 0x03 0x02; 0x03 0x03 … 0x00 0x03…\n",
       "  \"20\" => SnpData(379, 9327, UInt8[0x03 0x03 … 0x02 0x03; 0x03 0x03 … 0x02 0x03…\n",
       "  \"22\" => SnpData(379, 5938, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x03…\n",
       "  \"18\" => SnpData(379, 12242, UInt8[0x03 0x02 … 0x02 0x03; 0x03 0x03 … 0x03 0x0…"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "splitted = SnpArrays.split_plink(SnpArrays.datadir(\"EUR_subset\"); prefix=\"tmp.chr.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a SnpArray for chromosome 17."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SnpData(379, 11041, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x02; … ; 0x02 0x03 … 0x03 0x02; 0x03 0x03 … 0x00 0x03], 11041×6 DataFrames.DataFrame. Omitted printing of 2 columns\n",
       "│ Row   │ chromosome │ snpid       │ genetic_distance │ position │\n",
       "│       │ \u001b[90mString\u001b[39m     │ \u001b[90mString\u001b[39m      │ \u001b[90mFloat64\u001b[39m          │ \u001b[90mInt64\u001b[39m    │\n",
       "├───────┼────────────┼─────────────┼──────────────────┼──────────┤\n",
       "│ 1     │ 17         │ rs34151105  │ 0.0              │ 1665     │\n",
       "│ 2     │ 17         │ rs143500173 │ 0.0              │ 2748     │\n",
       "│ 3     │ 17         │ rs113560219 │ 0.0              │ 4702     │\n",
       "│ 4     │ 17         │ rs1882989   │ 5.6e-5           │ 15222    │\n",
       "│ 5     │ 17         │ rs8069133   │ 0.000499         │ 32311    │\n",
       "│ 6     │ 17         │ rs112221137 │ 0.000605         │ 36405    │\n",
       "│ 7     │ 17         │ rs34889101  │ 0.00062          │ 36975    │\n",
       "│ 8     │ 17         │ rs35840960  │ 0.000668         │ 38827    │\n",
       "│ 9     │ 17         │ rs144918387 │ 0.000775         │ 42965    │\n",
       "│ 10    │ 17         │ rs62057022  │ 0.000948         │ 49640    │\n",
       "⋮\n",
       "│ 11031 │ 17         │ rs7501742   │ 1.28534          │ 81050982 │\n",
       "│ 11032 │ 17         │ rs4056      │ 1.28534          │ 81063675 │\n",
       "│ 11033 │ 17         │ rs77182059  │ 1.28534          │ 81072675 │\n",
       "│ 11034 │ 17         │ rs71193771  │ 1.28534          │ 81082286 │\n",
       "│ 11035 │ 17         │ rs189547061 │ 1.28534          │ 81082531 │\n",
       "│ 11036 │ 17         │ rs187145448 │ 1.28534          │ 81087107 │\n",
       "│ 11037 │ 17         │ rs141554018 │ 1.28534          │ 81089800 │\n",
       "│ 11038 │ 17         │ rs144646847 │ 1.28534          │ 81103302 │\n",
       "│ 11039 │ 17         │ rs11489019  │ 1.28534          │ 81106522 │\n",
       "│ 11040 │ 17         │ rs139247601 │ 1.28534          │ 81130897 │\n",
       "│ 11041 │ 17         │ rs181245287 │ 1.28534          │ 81170925 │, 379×6 DataFrames.DataFrame\n",
       "│ Row │ fid       │ iid       │ father    │ mother    │ sex       │ phenotype │\n",
       "│     │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │\n",
       "├─────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┤\n",
       "│ 1   │ 1         │ HG00096   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 2   │ 2         │ HG00097   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 3   │ 3         │ HG00099   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 4   │ 4         │ HG00100   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 5   │ 5         │ HG00101   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 6   │ 6         │ HG00102   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 7   │ 7         │ HG00103   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 8   │ 8         │ HG00104   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 9   │ 9         │ HG00106   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 10  │ 10        │ HG00108   │ 0         │ 0         │ 1         │ 1         │\n",
       "⋮\n",
       "│ 369 │ 369       │ NA20810   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 370 │ 370       │ NA20811   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 371 │ 371       │ NA20812   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 372 │ 372       │ NA20813   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 373 │ 373       │ NA20814   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 374 │ 374       │ NA20815   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 375 │ 375       │ NA20816   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 376 │ 376       │ NA20818   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 377 │ 377       │ NA20819   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 378 │ 378       │ NA20826   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 379 │ 379       │ NA20828   │ 0         │ 0         │ 2         │ 1         │)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "piece = splitted[\"17\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can merge the split dictionary back into one SnpData using `merge_plink`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SnpData(379, 54051, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x03; … ; 0x03 0x03 … 0x03 0x02; 0x03 0x03 … 0x03 0x03], 54051×6 DataFrames.DataFrame. Omitted printing of 2 columns\n",
       "│ Row   │ chromosome │ snpid       │ genetic_distance │ position │\n",
       "│       │ \u001b[90mString\u001b[39m     │ \u001b[90mString\u001b[39m      │ \u001b[90mFloat64\u001b[39m          │ \u001b[90mInt64\u001b[39m    │\n",
       "├───────┼────────────┼─────────────┼──────────────────┼──────────┤\n",
       "│ 1     │ 17         │ rs34151105  │ 0.0              │ 1665     │\n",
       "│ 2     │ 17         │ rs143500173 │ 0.0              │ 2748     │\n",
       "│ 3     │ 17         │ rs113560219 │ 0.0              │ 4702     │\n",
       "│ 4     │ 17         │ rs1882989   │ 5.6e-5           │ 15222    │\n",
       "│ 5     │ 17         │ rs8069133   │ 0.000499         │ 32311    │\n",
       "│ 6     │ 17         │ rs112221137 │ 0.000605         │ 36405    │\n",
       "│ 7     │ 17         │ rs34889101  │ 0.00062          │ 36975    │\n",
       "│ 8     │ 17         │ rs35840960  │ 0.000668         │ 38827    │\n",
       "│ 9     │ 17         │ rs144918387 │ 0.000775         │ 42965    │\n",
       "│ 10    │ 17         │ rs62057022  │ 0.000948         │ 49640    │\n",
       "⋮\n",
       "│ 54041 │ 22         │ rs6010070   │ 0.750874         │ 51180959 │\n",
       "│ 54042 │ 22         │ rs6010072   │ 0.750878         │ 51181519 │\n",
       "│ 54043 │ 22         │ rs6009960   │ 0.750879         │ 51181568 │\n",
       "│ 54044 │ 22         │ rs56807126  │ 0.750879         │ 51181624 │\n",
       "│ 54045 │ 22         │ rs6010073   │ 0.750882         │ 51181986 │\n",
       "│ 54046 │ 22         │ rs184517959 │ 0.750893         │ 51183421 │\n",
       "│ 54047 │ 22         │ rs113391741 │ 0.750923         │ 51187440 │\n",
       "│ 54048 │ 22         │ rs151247655 │ 0.750938         │ 51189403 │\n",
       "│ 54049 │ 22         │ rs187225588 │ 0.751001         │ 51197602 │\n",
       "│ 54050 │ 22         │ rs9616967   │ 0.75109          │ 51209343 │\n",
       "│ 54051 │ 22         │ rs148755559 │ 0.751156         │ 51218133 │, 379×6 DataFrames.DataFrame\n",
       "│ Row │ fid       │ iid       │ father    │ mother    │ sex       │ phenotype │\n",
       "│     │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │\n",
       "├─────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┤\n",
       "│ 1   │ 1         │ HG00096   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 2   │ 2         │ HG00097   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 3   │ 3         │ HG00099   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 4   │ 4         │ HG00100   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 5   │ 5         │ HG00101   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 6   │ 6         │ HG00102   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 7   │ 7         │ HG00103   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 8   │ 8         │ HG00104   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 9   │ 9         │ HG00106   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 10  │ 10        │ HG00108   │ 0         │ 0         │ 1         │ 1         │\n",
       "⋮\n",
       "│ 369 │ 369       │ NA20810   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 370 │ 370       │ NA20811   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 371 │ 371       │ NA20812   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 372 │ 372       │ NA20813   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 373 │ 373       │ NA20814   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 374 │ 374       │ NA20815   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 375 │ 375       │ NA20816   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 376 │ 376       │ NA20818   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 377 │ 377       │ NA20819   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 378 │ 378       │ NA20826   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 379 │ 379       │ NA20828   │ 0         │ 0         │ 2         │ 1         │)"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "merged = SnpArrays.merge_plink(\"tmp.merged\", splitted) # write_plink is included here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also merge the plink formatted files based on their common prefix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SnpData(379, 54051, UInt8[0x03 0x03 … 0x03 0x03; 0x03 0x02 … 0x03 0x03; … ; 0x03 0x03 … 0x03 0x02; 0x03 0x03 … 0x03 0x03], 54051×6 DataFrames.DataFrame. Omitted printing of 2 columns\n",
       "│ Row   │ chromosome │ snpid       │ genetic_distance │ position │\n",
       "│       │ \u001b[90mString\u001b[39m     │ \u001b[90mString\u001b[39m      │ \u001b[90mFloat64\u001b[39m          │ \u001b[90mInt64\u001b[39m    │\n",
       "├───────┼────────────┼─────────────┼──────────────────┼──────────┤\n",
       "│ 1     │ 17         │ rs34151105  │ 0.0              │ 1665     │\n",
       "│ 2     │ 17         │ rs143500173 │ 0.0              │ 2748     │\n",
       "│ 3     │ 17         │ rs113560219 │ 0.0              │ 4702     │\n",
       "│ 4     │ 17         │ rs1882989   │ 5.6e-5           │ 15222    │\n",
       "│ 5     │ 17         │ rs8069133   │ 0.000499         │ 32311    │\n",
       "│ 6     │ 17         │ rs112221137 │ 0.000605         │ 36405    │\n",
       "│ 7     │ 17         │ rs34889101  │ 0.00062          │ 36975    │\n",
       "│ 8     │ 17         │ rs35840960  │ 0.000668         │ 38827    │\n",
       "│ 9     │ 17         │ rs144918387 │ 0.000775         │ 42965    │\n",
       "│ 10    │ 17         │ rs62057022  │ 0.000948         │ 49640    │\n",
       "⋮\n",
       "│ 54041 │ 22         │ rs6010070   │ 0.750874         │ 51180959 │\n",
       "│ 54042 │ 22         │ rs6010072   │ 0.750878         │ 51181519 │\n",
       "│ 54043 │ 22         │ rs6009960   │ 0.750879         │ 51181568 │\n",
       "│ 54044 │ 22         │ rs56807126  │ 0.750879         │ 51181624 │\n",
       "│ 54045 │ 22         │ rs6010073   │ 0.750882         │ 51181986 │\n",
       "│ 54046 │ 22         │ rs184517959 │ 0.750893         │ 51183421 │\n",
       "│ 54047 │ 22         │ rs113391741 │ 0.750923         │ 51187440 │\n",
       "│ 54048 │ 22         │ rs151247655 │ 0.750938         │ 51189403 │\n",
       "│ 54049 │ 22         │ rs187225588 │ 0.751001         │ 51197602 │\n",
       "│ 54050 │ 22         │ rs9616967   │ 0.75109          │ 51209343 │\n",
       "│ 54051 │ 22         │ rs148755559 │ 0.751156         │ 51218133 │, 379×6 DataFrames.DataFrame\n",
       "│ Row │ fid       │ iid       │ father    │ mother    │ sex       │ phenotype │\n",
       "│     │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │ \u001b[90mAbstract…\u001b[39m │\n",
       "├─────┼───────────┼───────────┼───────────┼───────────┼───────────┼───────────┤\n",
       "│ 1   │ 1         │ HG00096   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 2   │ 2         │ HG00097   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 3   │ 3         │ HG00099   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 4   │ 4         │ HG00100   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 5   │ 5         │ HG00101   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 6   │ 6         │ HG00102   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 7   │ 7         │ HG00103   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 8   │ 8         │ HG00104   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 9   │ 9         │ HG00106   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 10  │ 10        │ HG00108   │ 0         │ 0         │ 1         │ 1         │\n",
       "⋮\n",
       "│ 369 │ 369       │ NA20810   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 370 │ 370       │ NA20811   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 371 │ 371       │ NA20812   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 372 │ 372       │ NA20813   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 373 │ 373       │ NA20814   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 374 │ 374       │ NA20815   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 375 │ 375       │ NA20816   │ 0         │ 0         │ 1         │ 1         │\n",
       "│ 376 │ 376       │ NA20818   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 377 │ 377       │ NA20819   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 378 │ 378       │ NA20826   │ 0         │ 0         │ 2         │ 1         │\n",
       "│ 379 │ 379       │ NA20828   │ 0         │ 0         │ 2         │ 1         │)"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "merged_from_splitted_files = merge_plink(\"tmp.chr\"; des = \"tmp.merged.2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "isfile(\"tmp.merged.bim\") && rm(\"tmp.merged.bim\") # cleanup\n",
    "isfile(\"tmp.merged.fam\") && rm(\"tmp.merged.fam\")\n",
    "isfile(\"tmp.merged.bed\") && rm(\"tmp.merged.bed\")\n",
    "isfile(\"tmp.merged.2.bim\") && rm(\"tmp.merged.2.bim\")\n",
    "isfile(\"tmp.merged.2.fam\") && rm(\"tmp.merged.2.fam\")\n",
    "isfile(\"tmp.merged.2.bed\") && rm(\"tmp.merged.2.bed\")\n",
    "for k in keys(splitted)\n",
    "    isfile(\"tmp.chr.$(k).bim\") && rm(\"tmp.chr.$(k).bim\")\n",
    "    isfile(\"tmp.chr.$(k).fam\") && rm(\"tmp.chr.$(k).fam\")\n",
    "    isfile(\"tmp.chr.$(k).bed\") && rm(\"tmp.chr.$(k).bed\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "SnpArrays provides users with the ability to rapidly and accurately manipulate very large genetic data sets. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "84px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_position": {
    "height": "755px",
    "left": "0px",
    "right": "1099px",
    "top": "66px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}