python-package/Local_Ancestry_Example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial illustrates a simple usage of Loter for local ancestry inference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "* [I) Installation](#installation)\n",
    "* [II) Quick overview](#loter_usage)\n",
    "* [III) Details about local ancestry inference functions](#details)\n",
    "* [IV) Dealing with more than 2 ancestral populations](#3pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I) Installation <a class=\"anchor\" id=\"installation\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install Loter follow the procedure explained in the [README.md](../README.md) file.\n",
    "If you have already installed Loter, jump to the [next section](#loter_usage)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# II) Quick overview <a class=\"anchor\" id='loter_usage'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We show here how to perform local ancestry inference using Loter. In this section, we consider haplotypes from diploid admixed individuals. Loter can be used on non-diploid organisms, see [next section](#details)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a) Haplotype data <a  class=\"anchor\" id='hap_data'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use human haplotype data generated from the Hapmap dataset (see <https://www.ncbi.nlm.nih.gov/probe/docs/projhapmap/>). In particular, we use haplotypes of the 50,000 first SNPs from both copies of chromosome 1 for individuals sampled from different populations: 44 European individuals (*CEU*), 50 African inviduals (*YRI*) and 23 Mexican individuals (*MEX*).\n",
    "\n",
    "We detail in the next section how to load data from VCF files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "\n",
    "H_ceu = np.load(os.path.join(os.pardir, \"data\", \"H_ceu.npy\"), allow_pickle=True)\n",
    "H_yri = np.load(os.path.join(os.pardir, \"data\", \"H_yri.npy\"), allow_pickle=True)\n",
    "H_mex = np.load(os.path.join(os.pardir, \"data\", \"H_mex.npy\"), allow_pickle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`H_ceu`, `H_yri` and `H_mex` each represent a matrix of haplotypes sampled from a population (resp. Europe, Africa and Mexico). The matrix $H$ of dimension $(2\\times n, p)$ represents the haplotypes of $n$ individuals in rows and $p$ SNPs in columns. In diploid species, there are two haplotypes per individual (hence $2\\times n$ haplotypes for $n$ individuals). Haplotype matrices must be encoded with $\\{0,1\\}$ values. Each entry $H_{ij}$ represents the allele of SNP $j$ in haplotype $i$: 0 being the ancestral allele and 1 the derived (mutated) allele."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(88, 50000)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_ceu.shape # CEU haplotypes: 2*44 hapolotypes of the 1st 50,000 SNPs of chromosome 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 50000)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_yri.shape # YRI haplotypes: 2*50 hapolotypes of the 1st 50,000 SNPs of chromosome 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(46, 50000)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_mex.shape # MEX haplotypes: 2*23 hapolotypes of the 1st 50,000 SNPs of chromosome 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Haplotype matrices are generally obtained as a result of a **phasing algorithm** applied to genotype matrices.\n",
    "\n",
    "A genotype matrix $G$ of dimension $(n,p)$ represents genotype information regarding $n$ individuals in rows and $p$ SNPs in columns. Genotype matrices are encoded with ${0,1,2}$ values. Each entry $G_{ij}$ represent the alleles of a SNP $j$ on both copies of the corresponding chromosome in individual $i$: 0 = homozygous ancestral, 2 = homozygous mutated, 1 = heterozygous (ancestral on one chromosome and derived on the other).\n",
    "\n",
    "Data in this example (`H_ceu`, `H_yri` and `H_mex`) are the results of the phasing software [Beagle](https://faculty.washington.edu/browning/beagle/beagle.html) applied to Hapmap data and stored in Numpy array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b) Loading data from VCF Files <a  class=\"anchor\" id='load_data_vcf'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Haplotype data are generally stored in VCF format. We show here how lo load data from such files to be used in Loter.\n",
    "\n",
    "To load a VCF, we use [scikit allel](https://scikit-allel.readthedocs.io/en/latest/#). To install it, you can run `pip install scikit-allel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import allel\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`vcf2npy` convert a VCF File to a numpy matrix with values 0, 1 and 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vcf2npy(vcfpath):\n",
    "    callset = allel.read_vcf(vcfpath)\n",
    "    haplotypes_1 = callset['calldata/GT'][:,:,0]\n",
    "    haplotypes_2 = callset['calldata/GT'][:,:,1]\n",
    "    \n",
    "    m, n = haplotypes_1.shape\n",
    "    mat_haplo = np.empty((2*n, m))\n",
    "    mat_haplo[::2] = haplotypes_1.T\n",
    "    mat_haplo[1::2] = haplotypes_2.T\n",
    "    \n",
    "    return mat_haplo.astype(np.uint8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These 3 files are the results of the phasing software [Beagle](https://faculty.washington.edu/browning/beagle/beagle.html). Hence, we only get $\\{0,1\\}$ valued matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "H_ceu = vcf2npy(os.path.join(os.pardir, 'data', 'ceu.vcf'))\n",
    "H_yri = vcf2npy(os.path.join(os.pardir, 'data', 'yri.vcf'))\n",
    "H_mex = vcf2npy(os.path.join(os.pardir, 'data', 'mex.vcf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c) Run Loter <a class=\"anchor\" id='run_loter'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run Loter with default parameters.\n",
    "\n",
    "We assume, in this example, that Mexicans (*MEX*) are the result of an admixture between Europeans and Africans. We use present-day samples from European (*CEU*) and African (*YRI*) populations as a proxy for haplotypes of ancestral populations.\n",
    "\n",
    "The function `loter.locanc.local_ancestry.loter_smooth` takes at least two arguments (when using default values for the method hyper-parameters):\n",
    "* `l_H`: a list of \"ancestral\" or reference haplotypes matrices. Its length is equal to the number of ancestral populations.\n",
    "* `h_adm`: a matrix of admixed haplotypes\n",
    "\n",
    "and returns a matrix indicating the origin (regarding the ancetry) of each SNP in each haplotype: each entry $(i,j)$ corresponds to the index of the ancestral population in the list of ancestral haplotypes matrices from which the SNPs $j$ in haplotype $i$ originates (0 for the first population in `l_H` and $1$ for the second population in `l_H`).\n",
    "\n",
    "**Note:** The function `loter.locanc.local_ancestry.loter_smooth` uses a module that is able to correct phasing error (c.f. [next section](#details)). This module is specifically designed for diploid organism (i.e. that `h_adm` contains an even number of haplotypes).\n",
    "\n",
    "This computation may take several minutes to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import loter.locanc.local_ancestry as lc\n",
    "\n",
    "res_loter = lc.loter_smooth(l_H=[H_ceu, H_yri], h_adm=H_mex, num_threads=8) ## set the number of threads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(46, 50000)\n"
     ]
    }
   ],
   "source": [
    "print(res_loter.shape) # 2*23 haplotype and 50000 SNPs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## d) Visualize ancestry chunks <a class=\"anchor\" id='plot_result'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For data visualization, we use the `matplotlib` package. To install it, you can run `pip install matplotlib`.\n",
    "\n",
    "We display admixture chunks for each admixed haplotypes (in rows) and SNPs (in columns). Chunks of 0 correspond to CEU (first element of the list) and chunks of 1 correspond to YRI (second element of the list)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f0cc8f4fac8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD8CAYAAABJsn7AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAHeZJREFUeJzt3X/MJdV93/H3h+XHBhsH8NrWmiVlra4rry0H3BVgUTUE7LIQC/wHqdi2KUlpVmpMa8tpHagr0tBKtR21dq2i2KsYmUSJMXGSsqLrbigBJWlrYAk/zELWPGyQ2S7ylhgcV5YNu8+3f9x5YPZy750zc2fumTt8XtboufPrzDlz7x7G55zvHEUEZma2eCfkzoCZ2euVK2Azs0xcAZuZZeIK2MwsE1fAZmaZuAI2M8vEFbCZWQVJt0o6IunxKfsl6fOSViQ9Jul9Kem6AjYzq/ZlYPuM/ZcDW4plJ/AbKYm6AjYzqxARfwJ8d8YhVwG/FSPfAE6XtLEq3RPnyZSk7cB/BtYBvxkRn5p1/Mk6JdbzhuO2vfO9PwDgW4+dmnTNd773B5XHphzT5vmzjh/fVy7vPNeZlO54WrOOr1OWuunMe/+baPuas763PlrL7yK/nya/rzppfp8Xno+ItzRN67KffkP81XePJR370GM/2g/8sLRpV0TsqnG5s4BnS+uHim3PzTpJTUORJa0DvgV8sLjYg8COiHhi2jlv0plxgS49btvew48AcNnbz0267t7Dj1Qem3JMm+fPOn58X7m881xnUrrjac06vk5Z6qYz7/1vou1rzvre+mgtv4v8fpr8vuqk+T/iaw9FxLamaW37yfXxwN6fSDp23canKq8l6Rzgroh4z4R9/w34DxHxZ8X6PcAnIuKhWWnO0wRxPrASEQcj4iXgdkaP4WZm2QWwmvi/FhwCzi6tbwIOV500TwU87ZH7OJJ2Stonad/L/GiOy5mZpQuCl+NY0tKC3cA/LkZDXAh8LyJmNj/AfG3AmrDtNe0ZRTvKLhg1QcxxPTOzWlp6ukXSV4CLgQ2SDgG/CpwEEBFfAPYAVwArwA+AX0hJd54KuNEj97i+tKnN04ZV57y2yjuezqR0m1yrL99H3/i+zJbaB9I0nSaC4FhLr9uNiB0V+wP4SN1052mCeBDYImmzpJOBaxg9hpuZ9cIqkbTk0vgJOCKOSroe2MtoGNqtEbG/tZyZmc0hgGMZK9cUc40Djog9jNo+zMx6J+fTbYpeRMLtPfzIK+Ms6xxfXiYdU0e57WneMcSp+7tuV6x7X9fOafP6fUyraXpt5yElzdRrrv2WJo1dbqrq/C5+v23e4wBejkhacpnrCdjMrK+CGHYThJlZbwUc63f96wrYzIZpFAnXb66AzWygxLGJ8WL90YsKuG5jvgfF15d6z/p6b/uQry7y0IdyTdPnvKUYdcK5AjYzW7jROGBXwGZmWaz6CdjMbPH8BGxmlkkgjvUj1myqXuRuVkRbU2uzAzRNs+l5VR0XdWfAmMdlbz/3uBkLml6zaSdp2zNUtJnGvN9TF1Fxda6Vcv1l7ERrO8+roaQlFz8Bm9kgBeKlWJc7GzO5AjazQRoFYvTi/+RP5QrYzAbLnXBmZhlEiGPhJ+BKXXYWLGKaoWXQJBJuka/OrNLG9XOXoYlpeV7GsuSw6idgM7PFG3XC9buK63fuzMwaciecmVlGx3oeityL/zyMB0zMGmSe0vbVRSDANIuYpmbSvRn/20XemgRyNFUVNDNr6qlZ26dNXVU3uGF837yBGqnBR9P2p+R1vLx9mNZpkWmvRcKlLLn4CdjMBmvVoyDMzBZv9DIeV8BmZgsXiJcdimxmtngROBDj9a5JR+B4J2I5jVn7urKoazTpfG0aqDBt/7TOtrbvwbxTRNXpjK5zvTrq5mHx5EAMM7McAj8Bm5ll4044M7MMgrwvW0/hCtjMBmk0LX2/q7jePJ933WHQlroRUW1qcq2UiLmUNOqe3+R689zLJhFjk64/vjTVpymNcmr6u2nn35U4lrjk0u//PJiZNRQ4Es7MLJu+z4jR7/88mJk1FCFW44SkJYWk7ZIOSFqRdMOE/T8h6V5JD0t6TNIVVWn6CdjMBmnUCddOKLKkdcAtwAeBQ8CDknZHxBOlw/4NcEdE/IakrcAe4JxZ6boCHrhFTeWz9/AjvexI7Us+xvU1X7m1e19anRPufGAlIg4CSLoduAooV8ABvKn4/OPA4apEXQGb2SCNOuGS24A3SNpXWt8VEbtK62cBz5bWDwEXjKXxb4E/kvTPgTcAH6i6qCtgMxusGpFwz0fEthn7J9XkMba+A/hyRPxHSe8HflvSeyJidVqiroDNbJBajoQ7BJxdWt/Ea5sYrgO2A0TE/5a0HtgAHJmWaOV/HiTdKumIpMdL286UdLekp4q/Z9QoCDB92p1ZU8nM+jwrrS5VTXVTtT4pvSb7JmnSJlt3qp7x9OeZZmeRqsq5qDy2NSVR+d/K+DRSTaZ7mvRvre5vGfK3da9yQtKS4EFgi6TNkk4GrgF2jx3zbeBSAEnvAtYD/3dWoilX/jJFrV5yA3BPRGwB7inWzcx6IwJeXj0haalOK44C1wN7gScZjXbYL+lmSVcWh/0y8IuSHgW+Avx8RIw3UxynsgkiIv5E0jljm68CLi4+3wbcB/xKZSnMzBZk1ATRXqhDROxhNLSsvO2m0ucngIvqpNm0DfhtEfFccdHnJL112oGSdgI7AdZzasPLmZnV1/dIuM474YqhHLsA3qQzZz6Om5m1peYwtCyaVsDfkbSxePrdyIxevmmmdRDNarQvT1tTdVwT48EEfdN13troeOr/NDVp+p7Hcv6q/k00mZ6pj0E19bXbBNGFprnbDVxbfL4WuLOd7JiZtWe1mBeuasml8glY0lcYdbhtkHQI+FXgU8Adkq5jNPTiZ7vMpJlZXaNREEs+LX1E7Jiy69KW82Jm1hpPSWRmllHfp6XvTQt1atROVYdBOa2608pMS6fqnPHIo3HlMjXpEKmTn2nHzxu1VCetNqZBmpRmm+mm/I7Gr9+m1CjN1E7NOp1yVflKlfqWvCbXb+N+r42CSFly8ROwmQ1W30dBuAI2s0GKEEddAZuZ5eFOODOzDIYcCdeqcoN7n6Ju6ualaSTSvOkvQpPOzDbzmzsyq+1rdlWGPrzqsw2j+7MydzqugM3MMvA4YDOzjPo+DtgVsJkNUgQcTXjZek69yV35TWep7VjTpkoZ3zdpaqN55Gpnm+e6XbXHziPHFECpZk3R0+V16u4fN+93syxt3akciGFmloHbgM3MMgpXwGZmebgTzswsgwiPA04yPnA/teF+1tvLupiyqMn5fez8mqZvHWGT5Jo2Kndn0rLKO82XONbzURC9qIDNzLrgNmAzswz8Lggzs1xi1A7cZ66AzWyw+j4Kot8t1BNUTVc0aWqWtq43KZKuy46rSdFYbaU1Sddvc5uVbtOO1670pdOtST7KU2RN+t67iGRsS5v3PYpOuJQlFz8Bm9lguQnCzCwTj4IwM8sgwhWwmVk2Hoa25HJPhdOWZc67NdP1d74Mvym3AZuZZRCIVYcim5nl0fMH4OUbB2xmlqTohEtZUkjaLumApBVJN0w55u9LekLSfkm/W5Vmbyrg1EHdKe1ObU2LPm2qo7r5GT+/ramRZuULaBSUMq08TYJPZt2/vr95bVH5q/r95JgGqc10mvz7aFUkLhUkrQNuAS4HtgI7JG0dO2YLcCNwUUS8G/hYVbq9qYDNzNrW4hPw+cBKRByMiJeA24Grxo75ReCWiHhhdO04UpWoK2AzG6QAVleVtAAbJO0rLTvHkjsLeLa0fqjYVvZO4J2S/qekb0jaXpVHd8KZ2TAFkD4O+PmI2DZj/6SExhsvTgS2ABcDm4A/lfSeiHhxWqJ+AjazwYpIWxIcAs4urW8CDk845s6IeDki/hI4wKhCnqo3FXAfB3X3LU9dTrM0S9OOkVz5NXtFS51wwIPAFkmbJZ0MXAPsHjvmvwI/DSBpA6MmiYOzEnUThJkNVPoQsyoRcVTS9cBeYB1wa0Tsl3QzsC8idhf7/p6kJ4BjwL+KiL+ala4rYDMbrhYjMSJiD7BnbNtNpc8BfLxYkrgCNrNhCojVfr+Mp7INWNLZku6V9GQR3fHRYvuZku6W9FTx94zus2tmVocSlzxSOuGOAr8cEe8CLgQ+UkSA3ADcExFbgHuK9VbM6vQpR1KNR5VNW580tU/TiLTxzqPx9KumERqfNqnuteqoU87x+9bUpO9g1jGT8lrn+vPmdzytNakRanXzmnL+pN9U1e8q9befksdpeZv0vVZdp+71W++cba8TrhOVFXBEPBcRf158/j7wJKMByFcBtxWH3QZ8uKtMmpk10vMKuFYbsKRzgPOA+4G3RcRzMKqkJb11yjk7gZ0A6zl1nryamaWrF4iRRXIFLOmNwO8DH4uIv5bSChYRu4BdAG/SmX1/O5yZDUjfX8ieFIgh6SRGle/vRMQfFJu/I2ljsX8jUPniCTOzhVpV2pJJ5ROwRo+6XwKejIj/VNq1G7gW+FTx9862MpUSQTV+zKxXL846tup6VS57+7kzOxr6EA22lsdJ15v1Gsuqjp9p1xq/7rT942k2uR/T8tq1uq88Tb0P0zppq64zfh9Sv+tZac36d7IskYzq+RNwShPERcDPAd+UtPYN/mtGFe8dkq4Dvg38bDdZNDNrIHMHW4rKCjgi/ozpA+UubTc7ZmZt0XA64czMlk7Pn4B78za0snnb88ptWE3bCJsGZ7QRPNHEtPOq2irX7tG0dsdJ97JOkMKkduDywP3xtvtpQRxVgQl1VZ0/7/5xVQE4ddt+q77PtTyOt7GntiWnHLMUUxKtJi6Z+AnYzIZpSOOAzcyWzRBGQZiZLaeeV8C9bAM2M3s96OUTcBuDvOcNtsg10Lzpdfs0MH68s66q0yal86ntIIA2OqP6bhnL0Hae3QRhZpZDkDXMOIUrYDMbLj8Bm5nl4SYIM7Ncel4BL90oiJRImSbRNLOiteqc2zRKq4uoonK6bUUYVb39rerc8t+yeabRqZL63ZbL1uRNcF2p82azNqcHakPO+wYMa0YMM7NloXAThJlZPh4FYWaWh5+Al0STWQhS0mo6Bf0k02a1SJVyftV9aLMtOWVbF9caL8+k76vNWVPq5K1v6c17rezBIK6AzcwycBuwmVlGroDNzPJQxpetp1i6ccBmZkPRiwq46fQ/8xwzTx6mXS+1k6yNjqwmU8Q3zdukqW2aqjvFTtdydxLNG9RSTqcqeKTtsnYVFDXPea/hQAwzswzcCWdmlpErYDOzTFwBm5ktnvAoiCS5O0LaykObUW+5rtWH78LqS+1kfV2JV1/IU7WkkLRd0gFJK5JumHHc1ZJC0raqNHtRAZuZdaKlURCS1gG3AJcDW4EdkrZOOO404F8A96dkzxWwmQ1Xe8PQzgdWIuJgRLwE3A5cNeG4fwd8BvhhSqKugM1ssGo0QWyQtK+07BxL6izg2dL6oWLbq9eSzgPOjoi7UvPnTjgzG670URDPR8SsNttJLxZ+JXVJJwCfBX4++Yr04Am47SlzxtNuK9JonjzMWq9zbt08NSn/tOPb7ODJNU1Nm9ddtg6vLu75pGi8qmPWjlvIbyBGoyBSlgSHgLNL65uAw6X104D3APdJega4ENhd1RHnJ2AzG672xgE/CGyRtBn4P8A1wD945TIR3wM2rK1Lug/4lxGxb1ai2Z+Azcy60tYwtIg4ClwP7AWeBO6IiP2SbpZ0ZdP8+QnYzIarxUi4iNgD7BnbdtOUYy9OSTN7BTxrqvK20l5W8+a/rfIv87Q5fbhuHyyi7PMGB7XeLpz5TWcpslfAZmZdEH4bmplZNq6Azcxy6XkFXDkKQtJ6SQ9IelTSfkm/VmzfLOl+SU9J+qqkk7vPrplZDT2fESNlGNqPgEsi4ieBc4Htki4EPg18NiK2AC8A13WVyXKwRrmhfm19UrDD+DLp+CaN/lXnLCrIoM50MPMOfM8VOLFoOcpZ1XE1LU+pU1JN+7ewKHWnGysvc2v5bWhdqKyAY+T/FasnFUsAlwBfK7bfBny4kxyamTU1gCdgJK2T9AhwBLgbeBp4sRicDBNeTGFmlluLocidSOqEi4hjwLmSTgf+EHjXpMMmnVu8VWgnwHpObZhNM7P6+j4KolYockS8CNzH6EUTp0taq8DHX0xRPmdXRGyLiG0ncco8eTUzS5fa/NDnaeklvQV4OSJelPRjwAcYdcDdC1zN6MXE1wJ3dpXJWdFys7ZNe5vY2vYckVF9jU6zdu09/MhCp56adt3cv5Pc1+/7MLSUJoiNwG3FlBwnMHoJxV2SngBul/TvgYeBL3WYTzOzWgYRCRcRjwHnTdh+kNE0HWZmvaTVftfAjoQzs2Hyy3jMzPLpexNEL17I3uW0QZOiaubpGKg6t2lkU8qx41GAqWlMugd1ph5q0hHUlbaiuZpENNaN6mpD3Smi1qLfZp27yGi4tqbEamzZR0GYmS2rvj8BuwI2s+FyBWxmlkHkDTNO4Qr4dSr7APmGljXflqbN73cQ44DNzJZW9LsGdgVsZoPlJ2AzsxwciGFmlk/fO+F6EYjR1aD1OukuanB6nTzNCiDpqjNqWgDC2pIa0JEy7VPKdDl1p9OpSiMlj9O2l9+m18X0PpMCbZp8z7PeHvh6M4gXspuZLZ3AnXBmZrm4E87MLBdXwGZmi+dAjAVru9OhrWllcup6WpxpnXLTOoKavt1r2v6UjteqbeXP0zq/ytu7+k0s+2+tdyL8QnYzs2z6Xf+6Ajaz4XIThJlZDgG4CcLMLJN+17/9iIQbjwBqK8KoTlrzdIDUmSqozjQxdaYkanrdlP11vp/UaaLW/qZ0alV1ArbZebWW3qTOvTYjNqsiBcej7uqm3xe5p7NSpC1JaUnbJR2QtCLphgn7Py7pCUmPSbpH0t+oSrMXFbCZWRe0GklLZTrSOuAW4HJgK7BD0taxwx4GtkXEe4GvAZ+pStcVsJkNU+qEnGlPwOcDKxFxMCJeAm4HrjruchH3RsQPitVvAJuqEnUbsJkN0igQI7kReIOkfaX1XRGxq7R+FvBsaf0QcMGM9K4Dvl510V5UwF295WvetFLPrxOo0Fb5HAzQjnKQx6LKXuf3/nr7Plovb/qbzp6PiG0z9mvCtom1u6R/BGwDfqrqor2ogM3MulDjCbjKIeDs0vom4PBrrid9APgk8FMR8aOqRN0GbGbD1G4b8IPAFkmbJZ0MXAPsLh8g6Tzgi8CVEXEkJVE/AZvZQLX3LoiIOCrpemAvsA64NSL2S7oZ2BcRu4FfB94I/J4kgG9HxJWz0nUFbGbD1eIL2SNiD7BnbNtNpc8fqJvmIJsgFjl4fd4AiDWpb/Fqcn7KVELlY5tcK9WkDq8upvepMs+UP13kY/zzPNqe4isleKhPwR+vCE9JZGaWj6ckMjPLpN/1rytgMxsurfZ7XnpXwGY2TEGdQIwsBtkJB/k7V/pu3g6+XHrZ2dOSNu5zm/enjQjBnN+XCBRpSy5+Ajaz4XInnJlZJq6AzcwyGFIbsKR1kh6WdFexvlnS/ZKekvTVIj7azKw3tLqatORSpxPuo8CTpfVPA5+NiC3AC4zef5nVpOimtqckqkqvzQ6reTswyvejKqqpKp0uo55ydAjm7Fgs38+qfNS917nK1aeO2lfFqAkiZckkqQKWtAn4GeA3i3UBlzCadgPgNuDDXWTQzKyRoPcVcGob8OeATwCnFetvBl6MiKPF+iFGb4x/DUk7gZ0A6zm1eU7NzOpa9jZgSR8CjkTEQ+XNEw6d+J+RiNgVEdsiYttJnNIwm2Zm9Q1hHPBFwJWSrgDWA29i9ER8uqQTi6fgiW+HX7R53yjW1jW7SruNt5A1Oa7tt2vl+J5ymvQ2ulRDvi8L0fNhaJVPwBFxY0RsiohzGL0F/o8j4h8C9wJXF4ddC9zZWS7NzOqKgGOraUsm84Qi/wrwcUkrjNqEv9ROlszMWjKQTjgAIuI+4L7i80Hg/PazZGbWkp43QTgSzsyGKYCW5oTrymDfhjbJrKlv1vaVl/K+8t+qc2dda9L2JtcdP27SMZOm/ZnWqTMr39OmL6pr1nldvTWrTrqzjp0nGGVSZ2pbgT9tndNEShnq/lbmCRh6rYBYTVsy8ROwmQ1TkLWDLYUrYDMbLrcBm5ll4grYzCyHvEPMUryuKuBZHQIp++qcv7Y+reOmqtNs1vYmnWJ1IuFyTCPTVcRXWx2GbUYlzpOPrs4ZpAA8KaeZWSZ+AjYzyyE8CsLMLIuAyDjGN4UrYDMbLkfCVauKvJo3vdRzUvJRp/Ns2r61SKjxZfw60zrfUiOpqiL/pl2r7WlyUvI1HkU47ZyUaMZJ58zKS2oa0yIe51F1fp8j4VLy3jQvrXUkDullPGZmSyPCoyDMzLLxKAgzsxyCOHYsdyZm6kUFnGuQ+6LNO/16m/elTptvF5q0mzdNsyqQZdJxqdPGz3vv2r73OYJopskV1POKJXgdZS8qYDOzTvR8GFovRkGYmbUtgFiNpCWFpO2SDkhakXTDhP2nSPpqsf9+SedUpekK2MyGKdp7IbukdcAtwOXAVmCHpK1jh10HvBARfxP4LPDpqnRdAZvZYMWxY0lLgvOBlYg4GBEvAbcDV40dcxVwW/H5a8ClkjQrUcUCh2lI+j5wYGEXXIwNwPO5M9EBl2t5DLFMAH8rIk5rerKk/87o3qRYD/ywtL4rInaV0roa2B4R/7RY/znggoi4vnTM48Uxh4r1p4tjpn43i+6EOxAR2xZ8zU5J2je0MoHLtUyGWCYYlWue8yNie1t5ASY9yY4/vaYccxw3QZiZVTsEnF1a3wQcnnaMpBOBHwe+OytRV8BmZtUeBLZI2izpZOAaYPfYMbuBa4vPVwN/HBVtvItugthVfcjSGWKZwOVaJkMsE/SoXBFxVNL1wF5gHXBrROyXdDOwLyJ2A18CflvSCqMn32uq0l1oJ5yZmb3KTRBmZpm4AjYzy2QhFXBVCF8fSLpV0pFiLN/atjMl3S3pqeLvGcV2Sfp8UZ7HJL2vdM61xfFPSbq2tP1vS/pmcc7nqwZot1SmsyXdK+lJSfslfXQg5Vov6QFJjxbl+rVi++YiBPSpIiT05GL71BBRSTcW2w9Iuqy0PctvVtI6SQ9LumtAZXqm+I08sja0bNl/g62JiE4XRg3WTwPvAE4GHgW2dn3dBvn8u8D7gMdL2z4D3FB8vgH4dPH5CuDrjMb9XQjcX2w/EzhY/D2j+HxGse8B4P3FOV8HLl9AmTYC7ys+nwZ8i1EY5bKXS8Abi88nAfcX+b0DuKbY/gXgnxWffwn4QvH5GuCrxeetxe/xFGBz8Ttdl/M3C3wc+F3grmJ9CGV6Btgwtm2pf4Ot3ZsF3Pz3A3tL6zcCN+Yu+JS8nsPxFfABYGPxeSOjQBKALwI7xo8DdgBfLG3/YrFtI/AXpe3HHbfA8t0JfHBI5QJOBf4cuIBRNNiJ4787Rj3X7y8+n1gcp/Hf4tpxuX6zjMaW3gNcAtxV5HGpy1Rc6xleWwEP5jc4z7KIJoizgGdL64eKbcvgbRHxHEDx963F9mllmrX90ITtC1P8X9TzGD0tLn25iv+r/ghwBLib0dPdixFxdEJeXsl/sf97wJupX96ufQ74BLD2dpg3s/xlglE02B9JekjSzmLb0v8G27CIccC1w/OWwLQy1d2+EJLeCPw+8LGI+OsZTWRLU66IOAacK+l04A+Bd83IS938T3ow6bRckj4EHImIhyRdvLZ5Rj56X6aSiyLisKS3AndL+osZxy7Nb7ANi3gCTgnh66vvSNoIUPw9UmyfVqZZ2zdN2N45SScxqnx/JyL+oNi89OVaExEvAvcxai88XaMQ0PG8TAsRrVveLl0EXCnpGUZv2rqE0RPxMpcJgIg4XPw9wug/luczoN/gXBbQ/nMiowbzzbza+P/u3G0vU/J6Dse3Af86x3cUfKb4/DMc31HwQLH9TOAvGXUSnFF8PrPY92Bx7FpHwRULKI+A3wI+N7Z92cv1FuD04vOPAX8KfAj4PY7vsPql4vNHOL7D6o7i87s5vsPqIKPOqqy/WeBiXu2EW+oyAW8ATit9/l/A9mX/DbZ2fxb0g7qCUQ/808Ancxd6Sh6/AjwHvMzov6rXMWpTuwd4qvi79oWL0cuZnwa+CWwrpfNPgJVi+YXS9m3A48U5/4UiCrHjMv0dRv937DHgkWK5YgDlei/wcFGux4Gbiu3vYNQjvlJUXKcU29cX6yvF/neU0vpkkfcDlHrPc/5mOb4CXuoyFfl/tFj2r1132X+DbS0ORTYzy8SRcGZmmbgCNjPLxBWwmVkmroDNzDJxBWxmlokrYDOzTFwBm5ll8v8BncAQIu44b8kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.imshow(res_loter, interpolation='nearest', aspect='auto')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can save the result. Each line correspond to an haplotype. The SNPs are space delimited."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"result_mex.txt\", res_loter, fmt=\"%i\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# III) Details about local ancestry inference functions <a class=\"anchor\" id='details'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section concerns more advanced users that would like to change default parameter values and understand the algorithm behing the local ancestry inference approach in Loter. We recommend to use default function with default parameter values, designed to give the best results for a broad range of species.\n",
    "\n",
    "All details regarding local ancestry inference can be found in this reference: Dias-Alves, T., Mairal, J., Blum, M.G.B., 2018. Loter: A Software Package to Infer Local Ancestry for a Wide Range of Species. Mol Biol Evol 35, 2318–2326. https://doi.org/10.1093/molbev/msy126\n",
    "\n",
    "We use the same example as previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "\n",
    "H_ceu = np.load(os.path.join(os.pardir, \"data\", \"H_ceu.npy\"))\n",
    "H_yri = np.load(os.path.join(os.pardir, \"data\", \"H_yri.npy\"))\n",
    "H_mex = np.load(os.path.join(os.pardir, \"data\", \"H_mex.npy\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a) Loter with phase correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm for local ancestry inference in Loter is based on a penalty parameter $\\lambda$ (`lambda`). To avoid choosing it, we use a bagging approach: the model is trained for different values of $\\lambda$ over mutliple bootstrap resamplings, and the result is determined by a majority vote among the different trained models.\n",
    "\n",
    "In addition, the standard approach for local ancestry inference in Loter integrates a post-processing step based on a phase correction module, which accounts for phasing error (i.e. switches of homologous loci between the two copies of a chromosome in an individual). The phase correction step corresponds to a smoothing of the bagging results and is specifically designed for data sampled from diploid organisms.\n",
    "\n",
    "The function `loter.locanc.local_ancestry.loter_smooth` implements the method with both bagging and phase correction. Besides `l_H` and `h_adm` (c.f. [previously](#run_loter)), it takes several additional hyper-parameters:\n",
    "* `range_lambda`: list or 1-d array of candidate values (>0) for $\\lambda$\n",
    "* `threshold`: smoothing parameter in [0,1] for the phase correction module (=1 corresponds to no smoothing)\n",
    "* `rate_vote`: minimal bagging vote parameter in [0,1] (SNP with a majority vote lower than this threshold are considered indecisive and reassigned based on the reference population for nearest SNP with decisive vote)\n",
    "* `nb_bagging`: number of resampling in the bagging (positive integer)\n",
    "* `num_threads`: number of threads for parallel computations (requires OpenMP, c.f. [README.md](../README.md))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import loter.locanc.local_ancestry as lc\n",
    "\n",
    "res_loter = lc.loter_smooth(\n",
    "    l_H=[H_ceu, H_yri], h_adm=H_mex, range_lambda=np.arange(1.5, 5.5, 0.5),\n",
    "    threshold=0.90, rate_vote=0.5, nb_bagging=20, num_threads=8\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b) Loter without phase correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to run the local ancestry inference approach in Loter without the phase correction module. To do so, you should use the function `loter.locanc.local_ancestry.loter_local_ancestry`.\n",
    "\n",
    "Besides `l_H` and `h_adm` (c.f. [previously](#run_loter)), it takes several additional hyper-parameters:\n",
    "* `range_lambda`: list or 1-d array of candidate values (>0) for $\\lambda$\n",
    "* `rate_vote`: irrelevant in this mode\n",
    "* `nb_bagging`: number of resampling in the bagging (positive interger)\n",
    "* `num_threads`: number of threads for parallel computations (requires OpenMP, c.f. [README.md](../README.md))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import loter.locanc.local_ancestry as lc\n",
    "\n",
    "res_loter = lc.loter_local_ancestry(\n",
    "    l_H=[H_ceu, H_yri], h_adm=H_mex, range_lambda=np.arange(1.5, 5.5, 0.5),\n",
    "    nb_bagging=20, num_threads=8\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result containing haplotypic ancestry can be found in `res_loter`.\n",
    "\n",
    "In particular, `res_loter[0]` contains the ancestry and `res_loter[1]` contains the number of time that ancestry was picked in the bagging procedure.\n",
    "\n",
    "Again, we display admixture chunks. Chunks of 0 correspond to CEU (first element of the list) and chunks of 1 correspond to YRI (second element of the list)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.imshow(res_no_impute[0], interpolation='nearest', aspect='auto')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** The default version of the function `lc.loter_local_ancestry` can be used on non-diploid organisms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to get additional information about genotypic ancestry, only in the case of diploid organisms, you can change the `default` parameter value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_impute, res_no_impute = lc.loter_local_ancestry(\n",
    "    l_H=[H_ceu, H_yri], h_adm=H_mex, range_lambda=np.arange(1.5, 5.5, 0.5),\n",
    "    rate_vote=0.5, nb_bagging=20, num_threads=8, default=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the result containing haplotypic ancestry can be found in the second variable `res_no_impute`. In particular, `res_no_impute[0]` contains the ancestry and `res_no_impute[1]` contains the number of time that ancestry was picked in the bagging procedure.\n",
    "\n",
    "In addition, `res_impute` contains the genotypic ancestry. The genotypic ancestry corresponds to the paired haplotypic ancestries without order. For instance, if there are 3 ancestral populations, there are 6 possible ancestry values for genotypic ancestry.\n",
    "\n",
    "Denoting by $h1$ and $h2$ the haplotypic ancestry at a locus, the mapping function is the following:\n",
    "\n",
    "$$g: (h1, h2) -> (max(h1, h2) * (max(h1, h2) + 1))/2 + min(h1, h2)$$\n",
    "\n",
    "For example, let's consider 3 populations:\n",
    "* g((0, 0)) = (0 * 1)/2 + 0 = 0\n",
    "\n",
    "* g((0, 1)) = g((1, 0)) = (0 * 1)/2 + 1 = 1\n",
    "\n",
    "* g((1, 1)) = (1 * 2)/2 + 1 = 2\n",
    "\n",
    "* g((2, 0)) = g((0, 2)) = (2 * 3)/2 + 0 = 3\n",
    "\n",
    "* g((2, 1)) = g((1, 2)) = (2 * 3)/2 + 1 = 4\n",
    "\n",
    "* g((2, 2)) = (2 * 3)/2 + 2 = 5\n",
    "\n",
    "So genotypic ancestry compresses the pair information into one single number.\n",
    "It is useful for the computation of the genotypic ancestry error rate.\n",
    "\n",
    "In this mode, the `rate_vote` input parameter is the minimal bagging vote parameter in [0,1]. The SNP with a majority vote lower than this threshold are considered indecisive and reassigned based on the reference population for nearest SNP with decisive vote.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IV) Dealing with more than 2 ancestral populations <a class=\"anchor\" id='3pop'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is not currently possible to account for phasing errors in Loter local ancestry inference procedure when there are more than two ancestral populations. Haplotypic length of ancestry tracts will be affected and should not be used when there are more than 2 ancestral populations. Other statistics at the scale of genotypes (mean ancestry in the population, mean ancestry per individual) are accurately inferred.\n",
    "\n",
    "Thus, you can used the bagging-only procedure `loter.locanc.local_ancestry.loter_local_ancestry` described [previously](#details) (in default mode or not).\n",
    "\n",
    "Here, we show a toy example where Mexican haplotypes are both the admixed haplotypes and in the reference list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import loter.locanc.local_ancestry as lc\n",
    "\n",
    "res_loter = lc.loter_local_ancestry([H_ceu, H_yri, H_mex], H_mex, num_threads=8) ## set the number of threads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [previous section](#details) for details about the structure of `res_loter`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}