{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Donor match\n",
    "\n",
    "<span style=\"color:red\">Please updated to the latest version on GitHub by the following command</span>\n",
    "\n",
    "```bash\n",
    "pip install -U git+https://github.com/single-cell-genetics/vireo\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we show an example on how to align donors by their genotypes: \n",
    "1) between vireo and other 'omics data (e.g., SNP array, bulk RNA-seq, Exome-seq)\n",
    "2) between multiple batches all estimated by vireo.\n",
    "\n",
    "The idea is the same for the two cases: align the donors with giving least genotype difference, either using categorical genotype value or genotype probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vireoSNP version: 0.5.8\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import vireoSNP\n",
    "print(\"vireoSNP version: %s\" %vireoSNP.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Option 1: using wrap function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simple option is to use the wrap function `vireoSNP.vcf.match_VCF_samples`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape for Geno Prob in VCF1: (3784, 4, 3)\n",
      "Shape for Geno Prob in VCF2: (3784, 4, 3)\n",
      "n_variants in VCF1, VCF2 and matched: 3784, 3784, 3783\n",
      "aligned donors:\n",
      "['MantonCB1' 'MantonCB2' 'MantonCB3' 'MantonCB4']\n",
      "['donor2' 'donor1' 'donor3' 'donor0']\n"
     ]
    }
   ],
   "source": [
    "res = vireoSNP.vcf.match_VCF_samples('../data/donors.cellSNP.vcf.gz', \n",
    "                                     '../data/outs/cellSNP_noGT/GT_donors.vireo.vcf.gz',\n",
    "                                     GT_tag1 = 'PL', GT_tag2='PL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3gUlEQVR4nO3deXgV5dnH8e+dDUIWkpCDJhESgojKYhBQEFAKKiJoBURUUMCt1Ver4k6t4r5htW6taK0gKIKCArIoKgoFl4CAqFhZiywlmD1CCDn3+8cM4SQECGByJsn9ua5c5jzzzJx7nkzOb56ZIYqqYowxxnhNSLALMMYYYypjAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsqYaiIir4vIwzX4fmNEZGJNvZ8x1c0CyniCiFwqIl+KSJGIbHe/v0FEpIbrGCMiJSJSKCK5IrJYRLrWwPumiYi671soIv8TkVkics5Rbi/sN6zxUxHJEpF8EVkhIr8PWDY6oPZCEdkpIn4RSXSXJ4jI2yLyi4jsEJFJIhJblW1XUkeciLwmIttEpEBE/iMidwcsVxH5VkRCAtoeFpHXK4zN3lo3BK5vvMMCygSdiNwG/A14CjgWOAb4I9ANiAhCSW+rajTgAxYB0yoLShEJrYb3jnPf+xTgI2C6iIyohvc5EjcDSaoaC1wHTBSRJABVfVRVo/d+AU8AC1R1h7vuw0A80AJoifMzHlOVbVfiGSAaOAloDFwIrKnQJxm49BD7s3esLwPuE5HzDtHf1DALKBNUItIYeBC4QVXfUdUCdXyjqkNVtdjt10BExorIf93ZxT9EJNJd1lNEfhaR29zZ11YRGRn4HiIywT1D3ygi9waeXR+IqpYA43FCs4l7ye7vIjJbRIqA34nISSKywJ1tfSciF1bYTKKIfOSe6X8mIqlVGRdV3aaqf8P5EH9ib70ikiwi77r7sl5E/nSATXzu/jfXnSV0FZGWIvJJhVlMXFXqcWtaqap79r4EwoFmFfu5YX4lztjt1QJ4T1XzVTUPmA60OdxtuzoDb6pqjqr6VXW1qr5Toc+TwANVmUGq6hLgO6CtiCS6M9dcEckWkYVVOVZM9bCBN8HWFWgAvH+Ifo8DJwAZwPFACnBfwPJjcc6mU4CrgRdFJN5d9ry7LB04C+fDcySHICINgBHApoCZwOXAI0AM8CUwE/gQaArcBEwSkdYBmxkKPAQkAsuBSYd63wqmudtu7X5QzgRWuPvZG7hFRPpUst6Z7n/j3FnNEkCAx3BmFyfhBMCYgP19SUReOlgx7of3Lpx9XwBkVtKth1vzuwFtLwL9RSTe/bkMAuYcwbYBvgAeEZGRItLqAH2mAfk4P7+D7Y+ISDecsPwGuA34GWf2fAwwGicwTTCoqn3ZV9C+gGHAtgpti4FcYCfOB60ARUDLgD5dgfXu9z3dvmEBy7cDXYBQYDdwcsCyP+BcfqqsnjFu/1x3G58AHd1lrwMTAvr2ALYBIQFtbwFjAvpPDlgWDZQCzSp53zScD8KwCu0N3fZuwOnAfyssvwf4V0DtEw+2vQrrXgR8cwQ/s3CgLzDqAMv/CbxeoS0ZmA/43a+PgIjD3bbbJxInOJYCJTiX9/oGLFeck5jzgY04l4kf3ltTwNjkAjnAD8Cf3GUP4pwsHR/s3w37UptBmaD7BecyWNmlGFU9Q1Xj3GUhOGezjYCl7qWXXGCu2162Hd13iQjgV5xASMT50NsYsGwjzgzkQKaoapyqNlXVXqq6NGDZpoDvk3FmV/6DbLusv6oWAtnuelW1d1vZQCqQvHcM3HEYjXOmf0gicoyITBaRzSKSD0zEGZ/DoqolqjoHOLfiJU0RaQQMpvzlPYApwH9wZp6xwFr3/au87YA+O9W559URaOJue6qIJFToNxtnNvSHA+xKoqrGq+pJqvqc2/YUTuB9KCLr7OGJ4LKAMsG2BCgGDvjUFrADZ4bUxg2OOFVtrM4N7kPZgXOWHXjvpzmw+QjrDbzcswVoVuEeRcVtl91HEZFoIMFdr6oG4MzkfsQJu/UBYxCnqjGqev4h6tzrUbe9nToPIwzDmZ0eqTCcBx4q1puNc4kuUAbwsqoWuUH9D5wZzuFsez+qmo+zX1E497kq+jNOiDc61Lbc7RWo6m2qmo7z8MUoEeldlXXNb88CygSVquYCDwAvicjFIhIjIiEikoHzoYM7Q3kFeEZEmgKISMoB7r1U3H4pzhn2I+62U4FRVHL2fgS+xJmp3Ski4SLSE7gAmBzQ53wR6S4iETj3or5Q1U37bakCd7ZzI3A/cI87Bl8BBSJyl4hEikioiLQVkc6VbCIL51JaekBbDFAI5IlICnBHVXdURE4Ukb7u+4aLyDCcy6+fVeg6HOcyaMWA/Bq4xl0/EudJvZWHue29tfxFRDqLSISINMR5AjAXJ8TLUdUFwCq3rqrsZ38ROd590CMP55Ks/xCrmWpiAWWCTlWfxAmNO4H/uV8vA3fh3I/C/X4N8IV7eWo+0Hr/rVXqJpx7WOtwHht/E3jtN6h7N04g9cWZqb0EXKmqqwO6vYkTMtlAR5xZy8Hkuk8Ifoszwxisqq+571cK9MeZjax33/NVnAdAKtb2K87DHP92Lwd2wTkROBXng/cDnAcJyojzZOQ/DlCX4Nzj2o4TfjcDQ1R1WcD6KUAvYEIl61+Fc+/nZ5wZZjr7QuOQ2664e8C/3P3fApwD9HNnZpW5F2fmWhWtcI6tQpzZ/Uuq+mkV1zW/Mdn/RMcYY4wJPptBGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ40m/2l47rqobRjTSmyX4PSZmjFVqjf6S8fvHbg0/V4diEY4NdQp216pvlO1TVV7HdAuoQYpo0ZuDoq4JdRt0TG4w/Ul5PFJUEu4I66e4htwW7hDorPTZhY2XtdonPGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeFJYsAswVbPpu7UsnvIR6ldO7HYKGeedUWm/dctWM3/cNAbcMxJfahLb129h4aTZAKhCx/49aNGhdU2W7mmbVvzE4glzUb+fE393KhkX9qi037qvvmf+s1MY8PC1+NJT2FXwKx/9bQpZazdzwpkZdB/Zr4Yr975Nq9ayePKHzjHbI4OMvgc4ZpeuZv4/3mXAn0fiS0tm+/rNLJzgHrNAxwt60OLUE2uwcm/77KP5PHjXaPylpVwy/AquH3VLueWvvvAiU8a/QWhYGAmJiTz54vOkNG/G9yu/5S+33kZhQQEhoaH83+2j6D9oYHB2oooOGVAiosAkVR3mvg4DtgJfqmr/I3lTERmtqo8eyboB27gduAbYBZQAz6vqBBFZACQBO4EGwDOqOq7CujOAdFVtezQ11BS/38+it+bR7+bLiIqPZfpj/yK1fSvik33l+u3eVcyqT76maYvksraEFB8D7rmKkNAQfs0r5J2HXyW1fStCQm3y7Pf7WfSv2fS75wqimsQy/d5XSD21NfHHNS3Xb/fOYlbN/YKmx6eUtYWGh9H54t+R/fN2sjdtr+nSPc/v97Pozbn0u/Vy55h95DVSTznAMfvxV+WP2eSmDLj3aueYzS3gnQdfJfWUE+yYBUpLS7n/tjuZ8P40jk1J5qKevTn7/PNodeK+AG/Tvj3vf/YJkY0aMfHV13j8vvt5/vXXaBgZydiX/06L41vyv61bufDMXpzZuzexcY2DuEcHV5WfeBHQVkQi3dfnAJuP8n1HH83KIvJHt47TVDUD6A1IQJehbns34AkRiQhYdyBQeDTvX9OyNmyhcdN4Yn3xhIaF0rLzyWxY+dN+/TJnfE5Gn66Ehu077wiLCC/7xd5TsqfcINV3WWs20/iYBGKPSSA0LIyWXduyYemP+/XLnPoJGRd0JzR837iGN4zg2BNTy7WZfbLWb6GxL6H8Mbv8P/v1y3zvMzLO61puHMMaBB6zpYgdtWVWZC4lNb0FzVukERERQf9BA/nogznl+nQ9sweRjRoB0KFzJ7Zt3gJAeqvjaXF8SwCOSUqiiS+RX3bsqNkdOExVPSWZDey9hnEZ8NbeBSJymogsEZFvRGSxiLR220eIyDQRmSsiP4nIk27740CkiCwXkUlu2ygRWeV+3eK2pYnIDyLyioh8JyIfBoTkaOB6Vc0HUNV8VR1fSd3ROAFb6m4zGhgFPFzF/faEopwCouJjy15HxcVQlFNQrs+O/26jMCef5u2O32/97es3M/WBcbzz0Ct0v7yvnYm6inLyiWoSMK4JsRRl55frs2P9Fgp/yad5hxNqurxarSi3gKiEmLLXUfGxFOVWOGY3bnWO2fat9lt/+7rNTL3vZd55YBzdh51nx6xr29atJB23byaflJzM/7ZsPWD/KRMmctY5Z+/XviJzKSW7d5Oa3qJa6vytVPWnPhm4VEQaAu2BLwOWrQZ6qGoH4D4g8NJdBjAEaAcMEZFmqno3sFNVM1R1qIh0BEYCpwNdgGtFpIO7fivgRVVtA+QCg0QkFohR1XUHqXeSiKwEfgQeUtVSt/0h4Gng14PtrIhcJyKZIpK5q/CgXT1B/cqSqfPpOqh3pcubtkhh8P3XMeDukSyfu5g9JXtquMLaSf1+lkycR9dh5wa7lDpH/cqSKfPpOnj/D0+ApukpDH7wDwz481Usn2PH7JF4b/IUvv3mG669+aZy7du3bWPUddfz5EsvEBLi7eCv0vUJVV0pImk4s6fZFRY3BsaLSCuce5rhAcs+VtU8ABH5HkgFNlVYvzswXVWL3H7TgB7ADGC9qi53+y0F0qq0V84lvkwR8QGLRWQuEA+0VNVb3X052P6OA8YB+FKTtIrvWW2i4mMoytl3Zl+UW0BU/L6z05LiYrK3ZDHzr5MA2JlfyLyXptLnhsH4UpPK+sUnJRLeMIKcLVnl2uurqPhYin4JGNfsfKIS9s2oSnbtJnvTdmY+9DoAO/MKmTf2Lfrcfhm+9JSKmzMBouJiKMreN2MqysknKi7gmN3lHrNjJwLu2L4wlT43DsaXtu9+VHxSIuENIsjZvL1ce311bFISW3/ed4dl65YtHJO8/+/yok8X8OLYp3lrziwaNGhQ1l6Qn8/Vgy/ltvv+TIfTOtdIzUfjcC6gzwDGAj2BJgHtDwGfquoA94N/QcCy4oDvSw/z/SpbP1JV80WkUETSDzGLQlWzRGQZzuysCdBJRDa4dTQVkQWq2vMwa6pxvtRk8rbnkL8jl6i4GNZ+/T29rv592fKIyIYMf/rWstczn55Il4t740tNIn9HLtHxsYSEhlDwSx65234hpol3b4rWJF/LZPK2/UL+9hyiEmJYu2QVvW4cVLY8olFDho+7q+z1zIf+RZeh51o4VYEvLZm87dnkZ+USFe8es9dcVLY8olFDhj8zquz1zKfeoMvg3vjSksnPyiU6oeIxG1fzO+FB7TueyoZ169i0YSPHJCcx691pPPvPcs+A8d2Kldx78yj+NW0qib59D6Xs3r2bPw69kgGXDuH8i35fcdOedDiB8RqQq6rfikjPgPbG7HtoYkQVt1UiIuGqWgIsBF53700JMAC44hDrPwa8KCJD3MCKBgaq6oTATiLSCOgAPKmqU4C/u+1pwKzaEE4AIaEhdBtyLnOem4zf76f1GaeQkOwjc8ZnJKYmkXbKge+PbFuziRXzljjX8EXoflkfGkY3qsHqvSskNJRuI85nzuNv4PcrrXt2IOG4pmRO/YTE9GTSOh780eY3//QMJTuLKd1Tysalqzn/7iv2ewKwvgoJDaHb5X2Y8+xb+NVP626nkJDiI/N995jNOMQxO2exc8yGCN2HnkfDGDtmAcLCwhjz1JMMH3Ax/tJSBl8xlBNOOolnHn6Udqd24Ozz+/LYX+6nqKiIG4ePBCD5uON45e03mT3tPb7+92Jys7N5903nMYKn/v4iJ7dvF8xdOihRPfgVLBEpVNXoCm09gdtVtb+IdAXG4zyM8AEwTFXTRGQE0ElVb3TXmQWMVdUFIvIEcCGwzL0PNQq4yt38q6r6bECItHXXvx2IVtUxIiLAHcDVOI+YlwBPq+rESh4zf6PiI+0Vt30wvtQkHTj6qkN1M4crNuLQfcyRKSoJdgV10t1Dbgt2CXVWemzCUlXtVLH9kAFV31lAVRMLqOpjAVUtLKCqz4ECytuPcBhjjKm3LKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJYcEuoFaQYBdQB+3xB7uCOuvey+4Mdgl10sOTngh2CfWOzaCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ4UFuwBTNZtWrWXxlI9Qv3Ji91PIOO+MSvutW7aa+S9PY8A9I/GlJbF9/RYWTpwNgAId+/egRYfWNVi5t21auYbFE+ehfj8nntWBjAu6V9pv3dc/MP/5qQwYcw2+9GR+XrWWr6Z8QumeUkLDQjn90rNJOblFDVfvXQs+ms+YO++itLSUS4dfyf/dNqrc8leef4G3Xp9AWFgYCYlNGPv3FzmueXO+W7mSP98yioL8AkJDQ7nxjtu48OJBQdoLb6pPnwWHDCgRUWCSqg5zX4cBW4EvVbX/kbypiIxW1UePZN2AbdwOXAPsAkqA51V1gogsAJKAnUAD4BlVHeeuM9ddFgYsBP5PVUuPpo6a4Pf7WfTWPPrdchlR8bFMf+xfpLZvRXyyr1y/3buKWfXx1zRtkVzWlpDiY8DoqwgJDeHXvELeeehVUtu3IiTUJs9+v59FE+bQ785hRCXEMv3+V0k9tTXxKRXGdWcxqz78kqYtU8raGkY3os+tlxIVH0P2z9uZ/dQkhv3t1preBU8qLS3l3lG3MWnGeySlpHDBmb/jnPPP54STTizr06Z9ez5YuIDIRo1445VXefTe+3hpwutERjbimXEv0+L4lmzbupV+3c/irLN70zguLmj74yX17bOgKpUVAW1FJNJ9fQ6w+Sjfd/TRrCwif3TrOE1VM4DegAR0Geq2dwOeEJEIt/0SVT0FaAv4gMFHU0dNyVq/hcZN44n1xRMaFkrLTiezYcVP+/XLfP9zMs7rSmj4vvOOsIjwsgNwT8mecoNU32Wt3eyMa1N3XLu0YcOyH/frl/nuAjL6nVFuXBPTkoiKjwEgPsVH6e4SSkv21FjtXrY8cylp6emktmhBREQEF1w8kA8/+KBcnzPOOpPIRo0A6HBaZ7Zu2QJAeqvjaXF8SwCOTUoi0ecje8cvNbsDHlbfPguqGp2zgX7u95cBb+1dICKnicgSEflGRBaLSGu3fYSITBORuSLyk4g86bY/DkSKyHIRmeS2jRKRVe7XLW5bmoj8ICKviMh3IvJhQEiOBq5X1XwAVc1X1fGV1B2NE7Cle/u57WFABM5M1/OKcguIio8tex0VH0NRbkG5Pjv+u43CnHyatzt+v/W3r9/M1DHjeOfBV+g+tK+nz5hqUlFOAVFNGpe9jkqIpSinwrhu2Ephdh7NM0444HbWf/0DialJ5T4M6rNtW7aQfNy+2WZSSgr/27L1gP3fHv8GvzvnnP3al2cupWT3blLT7dLpXvXts6Cq1U0GLhWRhkB74MuAZauBHqraAbgPCLx0lwEMAdoBQ0SkmareDexU1QxVHSoiHYGRwOlAF+BaEengrt8KeFFV2wC5wCARiQViVHXdQeqdJCIrgR+BhwIv44nIPGA7UAC8U8X99zT1K0umzqfrxb0rXd60RQqDx1zHgHtGsnzuYvbYmX6VqF9Z8uaHdL3s3AP2yf55O19O+ZgeI/sdsI85sGmT32blN9/wh1v+VK79f9u2ccu11zH2Hy8REuLtD1EvqWufBVU65VPVlSKShjN7ml1hcWNgvIi0wpmRhAcs+1hV8wBE5HsgFdhUYf3uwHRVLXL7TQN6ADOA9aq63O23FEir0l45l/gyRcQHLBaRuaq60d2XPm7QTgJ6AR9VXFlErgOuA4hOiK24uMZFxcVQlJNf9roop4CouJiy1yXFxWRvzmLmXycBsDOvkHkvTaXPDYPxpSWV9YtPSiS8QQQ5m7PKtddXUfExFP2SV/a6KDu/7LIdQMmuYrJ/3s7Mx5zJ+c68QuY9O5k+t1yKLz2Zwux8PvrbFH533e+JPSahxuv3qmOTk9ny8767AFs3b+aY5P2Pt4WffsoLT45lytzZNGjQoKy9ID+fkYMGc8d9f+HU0zrXSM21RX37LDicaxIzgLFAT6BJQPtDwKeqOsANsQUBy4oDvi89zPerbP1IVc0XkUIRST/ELApVzRKRZTizs40B7btE5H3g91QSUO5DFeMAfKlJQb8M6EtLJm97Dvk7comKi2Ft5vf0uvr3ZcsjIhsy/K/7btDPfHoiXQb1xpeWRP6OXKLjYwkJDaHglzxyt/1CTGLjyt6m3vGlp5D3v2zys3KIio9l7Rff0ev6AWXLIxo1ZPhLd5S9nvnoeLpceg6+9GSKi3Yx9+m3OO2S3hx7QvNglO9Zp3Q8lfVr1/LfDRs4NjmZme9M47nXXi3XZ9WKFdzzp1t4Y/o0Epvuu8G/e/durr1sKAMvv4x+Ay6q4cq9r759FhxOYLwG5KrqtyLSM6C9MfsemhhRxW2ViEi4qpbgPE33untvSoABwBWHWP8x4EURGeIGVjQwUFUnBHYSkUZAB+BJt0+Mqm51n0Ts576354WEhtDt0nOZ87fJ+P1+Wnc7hYRkH5kzPiMxNYm0Uw58f2Tbmk2smLvEudYsQvfL+9AwulENVu9dIaEhdLuyL3OenIRfldZnZpBwXFMy3/2UxBbJpJ164Edwv5v/Ffn/y2bZ+5+z7P3PATj/zmFExkbVVPmeFRYWxkNPj+WKiwZSWlrKkCuG0frkk3j6oUdod2oHzu13Po/8+S/8WljE9VcMByC52XG8NmUys6ZN56t/LyY3O4d3Jr4JwNMvv0Sb9u2DuUueUd8+C0T14BMEESlU1egKbT2B21W1v4h0BcbjPIzwATBMVdNEZATQSVVvdNeZBYxV1QUi8gRwIbDMvQ81CrjK3fyrqvqsOxubpapt3fVvB6JVdYyICHAHcDXOI+YlwNOqOrGSx8zfUNVHReQYYJbbFgJ8Ctyqqge9COtLTdKBf77qYF3MkYi0Bwqqy70Dbwt2CXXSw5OeCHYJdda4Pzy6VFU7VWw/ZEDVdxZQ1cQCqtpYQFUPC6jqc6CAssdjjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCeFBbsAz1OgxB/sKuqcJ0b8Jdgl1Fl3TRgT7BLqpij7uKxpNoMyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBVUts+n4dbz84jslj/sHyD5ccsN+6b1Yz7sbHydq4tVx7YXYer416mhXzv6zuUmuV+fPm0bltG0496USeeerJ/Zb/e+FCzjq9M4mNGvL+tHfLLbvvnrvpmnEKp7dvx1233oKq1lTZtcKmVWt5+96/M3n0Syyfs/iA/dYtXc24ax8ha8MWALav38y7D7zCuw+8wjsPvML6ZatrquRaYdOKNbx9+4tMHvU8y2csOmC/dV/9wLhhD5K1zhnXXQW/MvOR8bx29WMsGj+npso9KocMKBFREZkY8DpMRLJEZNaRvqmIjD7SdQO2cbuIrBaR5SLytYhc6bYvEJEf3fYfROQ6t72RiHzgrvOdiDx+tDXUFL/fz6IpH9L3hksYfO+1rFn6PTlbd+zXb/euYlYtyKRpWvJ+y5ZM+4RmbdJrotxao7S0lDtu/hNTZ8zkixUrefftyaz+4ftyfZo1a8aLr/6Tiy+9tFz7l0sW8+WSxSxauozF3yznm6WZ/Pvzz2uyfE/z+/0senMufW++lMEP/oE1X31Hzpas/frt3lXMqo+/ommLfcdsQnJTBtx7NYPuv5bzb76UhRPn4C/112T5nuX3+1k0fg5977ycwU/ewJovviNncyXjurOYVfO+pGnLlLK20PAwOl/8O7pcfk5NlnxUqjKDKgLaikik+/ocYPNRvu9RBZSI/NGt4zRVzQB6AxLQZajb3g14QkQi3Paxqnoi0AHoJiJ9j6aOmpK1YSuNE+OJTYwjNCyUlqeezIaVP+3XL3PWQjLO6UJoWGi59g0r/kNMk8bEH5tYUyXXCku//or0li1JS08nIiKCgZcMYfbMmeX6NE9Lo2279oSElP9VERGKd+1i9+7dFBcXU1JSgq9p05os39Oy1m+hsS+BWF+8c8x2PpkNy/+zX7/M9z4j47yuhIaHlbWFNQgnJNQZ7z0lpUi5X+36LWvtZhofE09sU3dcu7Rhw9If9+uX+c4CMvqfUW5cwxtGcGzr5uXavK6ql/hmA/3c7y8D3tq7QEROE5ElIvKNiCwWkdZu+wgRmSYic0XkJxF50m1/HIh0ZziT3LZRIrLK/brFbUtzZ0CvuDOeDwNCcjRwvarmA6hqvqqOr6TuaJyALVXVX1X1U7f/bmAZcFwV9z+oivIKiIqPKXsdFR9DUV5BuT47Nm2jMCef5m2PL9deUryb5R99Qcfzu9dIrbXJ1i1bSGm27xBITklh6+aqnXud1qUrPc7qyYmpzTgxtRm9zjmX1iedVF2l1jpFuQVEJQQes7EU5VY4ZjdudY7Z9q32W3/7us1Mve9l3nlgHN2HnVcWWPVdUU4BUQmNy15HJcRSlFNhXNdvpTA7j+YdTqjp8n5zVf2pTwYuFZGGQHsg8EbGaqCHqnYA7gMeDViWAQwB2gFDRKSZqt4N7FTVDFUdKiIdgZHA6UAX4FoR6eCu3wp4UVXbALnAIBGJBWJUdd1B6p0kIiuBH4GHVLU0cKGIxAEXAB9Xcf89Tf3Kknc/puvAXvstW/rBItr16kx4g4hK1jRHat2aNfy4ejXfrdvA9+s3snDBpyxedOD7AaY89StLpsyn6+CzK13eND2FwQ/+gQF/vorlcxazp2RPDVdYO6lfWTLpQ7pefm6wS/lNVGmup6orRSQNZ/Y0u8LixsB4EWkFKBAesOxjVc0DEJHvgVRgU4X1uwPTVbXI7TcN6AHMANar6nK331IgrUp75VziyxQRH7BYROaq6kZ3+2E4M8DnDhRy7n2r6wCi42Or+JbVJ6pxTLmzpKKcAqIa7zs7LSkuJnvrDmb+7U0AduYXMe/ld+nzh0Fs37iFdctX8+V7n7J7ZzEiQmh4GG3P6ljj++E1ScnJbN70c9nrLZs3k5SScpA19pn1/vt0Pv10oqOjATi7z3l8/cUXnNHdZqoAUXExFGUHHrP5RMUFHLO7isneksXMsc7t7Z15hcx7YSp9bhyML+AeanxSIuENIsjZvL1ce30VFR9DUXZe2eui7PxyV1dKdhWT/fN2Zj7iXFDamVfIvL9Ops+oS/Gl177xO5yLkTOAsUBPoElA+0PAp6o6wA2xBQHLigO+Lz3M96ts/UhVzReRQhFJP8QsClXNEpFlOLOzjW7zOOAnVX32IOuNc/vha54U9EezfKlJ5GVlk78jl6i4GNYu+55eIy4sWx4R2ZDhT9xc9nrms5PoMqAXvtQkLrx1WFl75gcLCW8QYeHkOrVTZ9auWcPG9etJSklh2pS3eWXCG1Va97jmzZjwz39y6513oar8+/PP+eNNf6rmimsPX1oyeduzyc/KJSo+hrVff0+vay4qWx7RqCHDnxlV9nrmU2/QZXBvfGnJ5GflEp0QS0hoCAW/5JG77RdimsTV/E54kC89hbxt2eRvzyEqIZa1X3xHrxsGlC2PaNSQ4f+4o+z1zIfH0+Xyc2plOMHhBcZrQK6qfisiPQPaG7PvoYkRVdxWiYiEq2oJsBB43b03JcAA4IpDrP8Y8KKIDHEDKxoYqKoTAjuJSCOcByL23v962K33mirW6QkhoSF0u+Rc5rz4Nn5VWndpT0KSj8xZn5PYPIm0Sq7hm0MLCwvjyWf/xqD+/SgtLWXoiBGcdHIbHn1gDBmnduT8Cy5gWebXXHHJYHJzcpj7wQc8/uCDLFm+gt8PHMTnn35Kt1M7ICL0Pvdc+vbvH+xd8oyQ0BC6Xd6HOc++hV/9tO52CgkpPjLf/4zE1CTSMg58f2Tbmk2smLPYue8UInQfeh4NYxrVYPXeFRIaQrfhfZnz5CT8fqX1WRkkHNeUzHc+JbFFMmkdWx90/Tdv+RslO4sp3VPKxszVnH/3MOJTfDVU/eGTQ/3bDREpVNXoCm09gdtVtb+IdAXG4zyM8AEwTFXTRGQE0ElVb3TXmYXzFN0CEXkCuBBY5t6HGgVc5W7+VVV91p2NzVLVtu76twPRqjpGRAS4A7gaKHG/nlbViSKyAEgCdgINgDdU9VEROQ7n8uJq9s3MXlDVVw+2/77mSTrwrhEHHSNz+J649qFgl1Bn3TVhTLBLqJsiQw/dxxyRccMeXKqqnSq2HzKg6jsLqOphAVV9LKCqiQVUtTlQQNmzm8YYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTwoJdgNc1b5rCC9c/Euwy6pwbn7sn2CXUXZH2a10tsouDXUG9YzMoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40lhwS7AHNq8uXMZdeut+EtLGXn11dx5113lli/8/HNuGzWKb1euZOKbbzLo4ovLljUMD6dtu3YANGvWjOnvv1+jtXvdph/WsXj6x6gqJ57enoyzu1Tab92KH5n/+vsMuPUKfM2TKMjOY8rj/yTOlwBA09QkelzSpyZL97xN361l8ZSPUL9yYrdTyDjvjEr7rVu2mvnjpjHgnpH4UpPYvn4LCyfNBkAVOvbvQYsOrWuydE/b9P06Fk/7GPX7ObHrKWScc4BjdvmPzH/tPQbcfiW+5kll7YXZ+Ux59FU69u3GKb1Pr6myj8ghA0pEFJikqsPc12HAVuBLVe1/JG8qIqNV9dEjWTdgG7cD1wC7gBLgeVWdICILgCRgJ9AAeEZVx7nrPAJcCcSravTRvH9NKS0t5eabbmL2vHkcd9xxdD39dPpfcAEnn3xyWZ9mzZvz6muv8czTT++3fmRkJJnLltVkybWG3+9n0bvz6ffHS4iKi2H6MxNIbXs88ccmluu3e1cxqz5fStPUpHLtsU3iGHTHiBqsuPbw+/0semse/W6+jKj4WKY/9i9S27ciPtlXrt/uXcWs+uRrmrZILmtLSPEx4J6rCAkN4de8Qt55+FVS27ciJNQu+Pj9fhZN/Yh+/zfEOWbHjneO2aRKjtnPMvc7ZgGWTP+YZien11TJR6UqP/EioK2IRLqvzwE2H+X7jj6alUXkj24dp6lqBtAbkIAuQ932bsATIhLhts8ETjua965pX3/1FS1btiQ9PZ2IiAguGTKEmTNmlOuTlpZG+/btCQmxX+DDkfXfrTROjCM2MY7QsFBadjiJDavW7Ncvc84iMnqdTmiYXXCoqqwNW2jcNJ5YX7wztp1PZsPKn/brlznjczL6dC03tmER4WVhtKdkT7lf7Poua+NWGvsCjtlTT2LDt5WM6wcLyTi7C6Hh5Y/ZDSv/Q0yTuP1Owryqqp9os4F+7veXAW/tXSAip4nIEhH5RkQWi0hrt32EiEwTkbki8pOIPOm2Pw5EishyEZnkto0SkVXu1y1uW5qI/CAir4jIdyLyYUBIjgauV9V8AFXNV9XxldQdjROwpW6/L1R1a5VHxwM2b97Mcc2alb1OSUlhy+aqnx/s2rWLLqedRvczzuD9996rhgprr6LcQqLiYspeRzWOoSivoFyfHZu2UZhbQPM2LfdbvyA7j3fHvs7MF95k69pN1V5vbVKUU0BUfGzZ66i4GIpyKoztf7dRmJNP83bH77f+9vWbmfrAON556BW6X97XZk+uotwCouIqjGteYbk+BzpmS4p3s3z+l3Ts261Gav0tVPWUcDJwn4jMAtoDrwE93GWrgR6qukdEzgYeBQa5yzKADkAx8KOIPK+qd4vIje4MBxHpCIwETseZBX0pIp8BOUAr4DJVvVZEpgCDRGQGEKOq6w5S7yQRKXbXv0VVS6u4n7g1XQdcB9C8efPDWdVz1qxfT0pKCuvWraPP2WfTtl07Wrbc/8PW7E/9ypL3P6Xn5efvt6xRbBSX3/dHGkZFkrVpGx++Np3Bd11FRMMGQai09lG/smTqfHoOr/wuQdMWKQy+/zpytu5gweszada2JWHhNoM9FPUrS6Z/Qs+h/fZbtnTOItr17ER4g4hK1vSmKv3EVXWliKThzJ5mV1jcGBgvIq0ABcIDln2sqnkAIvI9kApUPNXsDkxX1SK33zSc8JsBrFfV5W6/pUBalfbKucSXKSI+YLGIzFXVjVVcF/ee1TiAjp06aVXXqw4pKSn8vGnfkG3evJnklJTDWh8gPT2dM886i+XffGMB5YqKi6Yod99ZfVFeAVGN982oSop3k71tBzNfcC4Y7CwoYt4/p9Hn6oH4mieVXZbyNTuW2CZx5G3PLnczuj6Lio+hKCe/7HVRbgFR8YFjW0z2lixm/nUSADvzC5n30lT63DAYX8B9k/ikRMIbRpCzJatce30VFRdDUW6FcW2873Z6SfFusrfuYObzbwKwM7+IeeOm0ee6gWzfsJV1y3/kyxkL2L2zGBEhNDyMtmd2rPH9qKrDOSWZAYwFegJNAtofAj5V1QFuiC0IWFYc8H3pYb5fZetHqmq+iBSKSPohZlGoapaILMOZnVU5oLykU+fOrFmzhvXuTGjK228zYeLEKq2bk5NDo0aNaNCgATt27GDJ4sXcdscd1Vxx7eFrlkReVg75v+QS1TiGtd/8QK9hF5Qtj4hswPCHbyp7PfOFt+hyYU98zZPYWfgrDRo1JCQkhPwdueTtyCGmSVwQ9sKbfKnJ5G3PIX9HLlFxMaz9+nt6Xf37suURkQ0Z/vStZa9nPj2RLhf3xpeaRP6OXKLjYwkJDaHglzxyt/1CTJPGwdgNz/E1r3DMLvuBXsMrHLOP/ans9czn3qTLRb/D1zyJC28ZWtaeOXsR4Q3CPR1OcHiB8RqQq6rfikjPgPbG7HtoYkQVt1UiIuGqWgIsBF53700JMAC44hDrPwa8KCJD3MCKBgaq6oTATiLSCOcS45NVrMtzwsLCePa55+jXty/+0lKGjxxJmzZtGHP//XTs2JELLryQzK+/ZvCgQeTk5PDBrFk8+MADrPj2W1b/8AM3XH89ISEh+P1+7rjzznJP/9V3IaEhdBt0NnNenorfr7Q+vR0JSYlkzllIYrNjSWvb6oDrbl27iaVzFhESGgoCPS4+l4ZRkQfsX9+EhIbQbci5zHluMn6/n9ZnnEJCso/MGZ+RmJpE2iknHHDdbWs2sWLeEue+kwjdL+tDw+hGNVi9d4WEhtDt4nOY89IU55jt0o6EJB+ZHywksfmxpLU78DFbG4nqwa9giUhhxUey3YC6XVX7i0hXYDzOwwgfAMNUNU1ERgCdVPVGd51ZwFhVXSAiTwAXAstUdaiIjAKucjf/qqo+687GZqlqW3f924FoVR0jIgLcAVyN84h5CfC0qk6s5DHzN/Y+0u4+qHE5kAxscd9rzMH2v2OnTvrFV18ddIzM4bvxuXuCXULdFWn3aqrF7sO6lW0Ow7g/PbFUVTtVbD9kQNV3FlDVwwKqGllAVQ8LqGpzoICyZzeNMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ4mqBrsGTxORLGBjsOuookRgR7CLqINsXKuPjW31qG3jmqqqvoqNFlB1iIhkqmqnYNdR19i4Vh8b2+pRV8bVLvEZY4zxJAsoY4wxnmQBVbeMC3YBdZSNa/Wxsa0edWJc7R6UMcYYT7IZlDHGGE+ygDLGGONJFlDGGGM8yQKqHhCRaBGJDnYddZmI2O/Sb0hEmohIfLDrqItEpLGIxAa7jqqwX6o6TkT6A28B80XkcrdNgltV3SAivUTkLgBV9VtI/TbcY/ZtYKaIXBXseuoSEbkYmAp8KCJXi0i7YNd0MGHBLsBUHxE5F3gcuBZIAB4VkW2q+klwK6v9RKQ38C7wtYjEqeo9e0NKVf3Brq+2EpF+wMPAdUAj4K8i8r6q/hLcymo/EUkB7gOGA42Bi4CWIhKvqp8Hs7YDsTO+OkpEIoDTgYdVdYmqfoDzbyO6B7eyOiMeeAy4EUgRkcfBZlJHwx23tsBdqvoV8COgwMMiMlREWga1wNqvIVAErFTVBcBLwK9AHxE5MZiFHYj9O6g6SEREVVVEmgObAQFKgauBs1X1UrdfqKqWBrHUWqfiDElEwoAM4CZgm6re5bZHq2phcKqsfQKO2RA35GNwZqiZwMfAUGCVqv41qIXWQnvH1v3+aSAXeEpVd4lIa+AWYImqTghelZWzM706JvBgBH4HJKvqHrdtNVDg9hsKDLL7UVXnjq3f/f5KEWmhqnuA5cALwLEicpeIjACuEpHQ4FVbe1Q4ZoeJSKqqFgDXqOpoVf0YeB04X0SiglZoLVQhnK4EfgJigctFJFJVfwSmAVeKSKMgllopC6g6JuBgHAxcgzNz2qsA+NU9UO/FOSO1KXQVVTK2u932PcA3wGjgcuAZYIHNTqumknHd47b/N6DbsUAhYPf3DkPA2F4CXAHMwzmhagPcLSLhOJerf8WDY2sPSdRBIpKBc8lpqqpucS9D+YEI4AacA3SAqq4OWpG1VMDYvqOqm92xLVXVPSLSB+dhlG6q+n0w66xtKhnXssvPInINcD0wQlV3BrHMWskd2xuB6aq6XkT+B2wDLsa5fBoFXKuqu4JXZeVsBlUHVHKZbiewCucSXoZ7ic8PrAVm4PyiWzhVQRXHdu8sdBfQ18Lp0KowrqXiOAE4ERiuqt/WeKG10EHG9iIROVVVf1XVj1X1epwZf29VXVbjhVaBPSRRy1W4xtwTZ1a8N3yuAI4D/qGq37r3RCLt5n3VHM7YBqXAWqqK4/qyqq50P2wjVLU4GLXWNocztkEp8DDZDKqWCzgY/w/nRv1g4EOca8yzgS3AHSLSTlVLLZyq7jDGtk3QiqyFqjiut4tIW3VYOFXRYYxtrThmLaBqOfcySBpwCc7lpT8AY4A73S4TgJVAVlAKrMUOY2ztH5EehsMY1x1BKbAWq2vHrAVULRR4jdk9Y9rifrUQkXBVnQJ8AvxJVTcBz6nqtuBUW7vY2FYPG9fqU5fH1gKqlqlwjTlNRJJVdTewCTgL5xozOP9ANx/AXW4Owca2eti4Vp+6Prb2kEQtUuFgHIXzlyHWAZ/hXG/+O86j5AqcjD35VGU2ttXDxrX61IextX8HVYsEHIynA6cC/YEGwCQgVFVHikgH4ATgXlXdEKxaaxsb2+ph41p96sPY2gyqFnGvNbcH3sD5G2XXuv9eJBXnT+gvUdWbg1ljbWVjWz1sXKtPfRhbuwflcRVvgKrqCmAs0Aro4t4E3QgMAdqLyDGV/EM9Uwkb2+ph41p96tvY2gyqlhDnj7u2ArYDE4F+wFXAg8AXqloiImHq/F04cxhsbKuHjWv1qS9jazMojxKRZBGJdL+/CefvlOUArXH+4OM8YDzwFNAJyv5oqTkEG9vqYeNaferr2NpDEh4kzv/58m5glYhMANKAm1X1S3f5aOBJVb1GRBrjPEJqqsDGtnrYuFaf+jy2NoPypi3AUpwp/FCcP1NyVsDyWbg/O1V9Ucv/bwnMwdnYVg8b1+pTb8fWZlAes/ffNojzv78+2f1aBtwkItmq+irQDkgTkTggb+/jpubgbGyrh41r9anvY2sPSXiQewP0dmAkzj++2wHEAYNwzpbOBIao6nfBqrG2srGtHjau1ac+j60FlAeJyINAgao+JSJ7/yeDXXGm+VOAQlW1P6R5BGxsq4eNa/Wpz2Nr96C8aRnQTUTaqOpuVX0WSAeaANl19WCsITa21cPGtfrU27G1e1DetADoDFwuIp8AkUAezl8hzg9mYXXAAmxsq8MCbFyrywLq6djaJT6PEpFkYKD7tQe4XWvJ/wXT62xsq4eNa/Wpr2NrAeVxIhKF83Oy/xPub8zGtnrYuFaf+ja2FlDGGGM8yR6SMMYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxnjS/wMYtJuIipCZdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "vireoSNP.plot.heat_matrix(res['matched_GPb_diff'], \n",
    "                          res['matched_donors1'] , \n",
    "                          res['matched_donors2'] )\n",
    "plt.title(\"Geno Prob Delta: %d SNPs\" %(res['matched_n_var']))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Option 2: using element functions for customised analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can use the build-in functions to perform each step, including [vireoSNP.load_VCF](https://github.com/single-cell-genetics/vireo/blob/master/vireoSNP/utils/vcf_utils.py#L68) to load VCF file and [vireoSNP.optimal_match](https://github.com/single-cell-genetics/vireo/blob/master/vireoSNP/utils/vireo_base.py#L143) to align donors that give minimal genotype differences by [Hungarian algorithm](https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.linear_sum_assignment.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare the data for this example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure your run vireo on the example data already\n",
    "\n",
    "```bash\n",
    "cd ../\n",
    "mkdir data/outs\n",
    "vireo -c data/cellSNP_mat -N 4 -o data/outs/cellSNP_noGT --randSeed 2\n",
    "```\n",
    "so you will have the estimated donor genotype `../data/outs/cellSNP_noGT/GT_donors.vireo.vcf.gz`\n",
    "\n",
    "Now, we will align the estimated donor genotype to another VCF estimated from bulk RNA-seq: `../data/donors.cellSNP.vcf.gz`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Often the known donor genotype VCF file can be very big, make sure you filter out unwanted variants first, as vireoSNP.load_VCF is not as efficient as [bcftools]()\n",
    "\n",
    "```bash\n",
    "bcftools view donor.vcf.gz -R cellSNP.cells.vcf.gz -Oz -o sub.vcf.gz\n",
    "```\n",
    "\n",
    "Also, add -s or -S for subsetting samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load genotypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load genotype from external omics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3784, 4, 3)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GT_tag0 = 'PL' # common ones: GT, GP, PL\n",
    "vcf_dat0 = vireoSNP.vcf.load_VCF(\"../data/donors.cellSNP.vcf.gz\",\n",
    "                                 biallelic_only=True, sparse=False, \n",
    "                                 format_list=[GT_tag0])\n",
    "\n",
    "GPb0_var_ids = np.array(vcf_dat0['variants'])\n",
    "GPb0_donor_ids = np.array(vcf_dat0['samples'])\n",
    "GPb0_tensor = vireoSNP.vcf.parse_donor_GPb(vcf_dat0['GenoINFO'][GT_tag0], GT_tag0)\n",
    "GPb0_tensor.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['chr1_1065797_G_C' 'chr1X_1217251_C_A' 'chr1_1230695_G_A'\n",
      " 'chr1_1722625_A_T']\n",
      "['MantonCB1' 'MantonCB2' 'MantonCB3' 'MantonCB4']\n"
     ]
    }
   ],
   "source": [
    "print(GPb0_var_ids[:4])\n",
    "print(GPb0_donor_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load donor genotype from vireo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3784, 4, 3)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GT_tag1 = 'PL' # common ones: GT, GP, PL\n",
    "vcf_dat1 = vireoSNP.vcf.load_VCF(\"../data/outs/cellSNP_noGT/GT_donors.vireo.vcf.gz\",\n",
    "                                 biallelic_only=True, sparse=False, \n",
    "                                 format_list=[GT_tag1])\n",
    "GPb1_var_ids = np.array(vcf_dat1['variants'])\n",
    "GPb1_donor_ids = np.array(vcf_dat1['samples'])\n",
    "GPb1_tensor = vireoSNP.vcf.parse_donor_GPb(vcf_dat1['GenoINFO'][GT_tag1], GT_tag1)\n",
    "GPb1_tensor.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['1_1065797_G_C' '1_1217251_C_A' '1_1230695_G_A' '1_1722625_A_T']\n",
      "['donor0' 'donor1' 'donor2' 'donor3']\n"
     ]
    }
   ],
   "source": [
    "print(GPb1_var_ids[:4])\n",
    "print(GPb1_donor_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Match variant ids in two VCF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# match variants\n",
    "mm_idx = vireoSNP.vcf.match_SNPs(GPb1_var_ids, GPb0_var_ids)\n",
    "idx1 = np.where(mm_idx != None)[0] #remove None for unmatched\n",
    "idx0 = mm_idx[idx1].astype(int)\n",
    "\n",
    "GPb1_var_ids_use = GPb1_var_ids[idx1]\n",
    "GPb0_var_ids_use = GPb0_var_ids[idx0]\n",
    "\n",
    "GPb1_tensor_use = GPb1_tensor[idx1]\n",
    "GPb0_tensor_use = GPb0_tensor[idx0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# match donors\n",
    "idx0, idx1, GPb_diff = vireoSNP.base.optimal_match(GPb0_tensor_use, GPb1_tensor_use, \n",
    "                                                   axis=1, return_delta=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "aligned donors:\n",
      "['MantonCB1' 'MantonCB2' 'MantonCB3' 'MantonCB4']\n",
      "['donor2' 'donor1' 'donor3' 'donor0']\n"
     ]
    }
   ],
   "source": [
    "print(\"aligned donors:\")\n",
    "print(GPb0_donor_ids[idx0])\n",
    "print(GPb1_donor_ids[idx1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.14587468, 0.18310327, 0.21593441, 0.22427656]), 0.7691889176247009)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GPb_diff[idx0, idx1], GPb_diff[idx0, idx1].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.43964819, 0.44643109, 0.14587468, 0.43192166],\n",
       "       [0.41247473, 0.18310327, 0.43770159, 0.42725593],\n",
       "       [0.42989822, 0.41561379, 0.42926244, 0.21593441],\n",
       "       [0.22427656, 0.40927258, 0.43158556, 0.42760225]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GPb_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the genotype probability difference\n",
    "\n",
    "Note, in this example data set the genotype estimation is not perfect as it is only based on ~250 cells in 10x data, namely not a decent coverage. Nevertheless, it is clear enough to find the donor identity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3gUlEQVR4nO3deXgV5dnH8e+dDUIWkpCDJhESgojKYhBQEFAKKiJoBURUUMCt1Ver4k6t4r5htW6taK0gKIKCArIoKgoFl4CAqFhZiywlmD1CCDn3+8cM4SQECGByJsn9ua5c5jzzzJx7nkzOb56ZIYqqYowxxnhNSLALMMYYYypjAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsqYaiIir4vIwzX4fmNEZGJNvZ8x1c0CyniCiFwqIl+KSJGIbHe/v0FEpIbrGCMiJSJSKCK5IrJYRLrWwPumiYi671soIv8TkVkics5Rbi/sN6zxUxHJEpF8EVkhIr8PWDY6oPZCEdkpIn4RSXSXJ4jI2yLyi4jsEJFJIhJblW1XUkeciLwmIttEpEBE/iMidwcsVxH5VkRCAtoeFpHXK4zN3lo3BK5vvMMCygSdiNwG/A14CjgWOAb4I9ANiAhCSW+rajTgAxYB0yoLShEJrYb3jnPf+xTgI2C6iIyohvc5EjcDSaoaC1wHTBSRJABVfVRVo/d+AU8AC1R1h7vuw0A80AJoifMzHlOVbVfiGSAaOAloDFwIrKnQJxm49BD7s3esLwPuE5HzDtHf1DALKBNUItIYeBC4QVXfUdUCdXyjqkNVtdjt10BExorIf93ZxT9EJNJd1lNEfhaR29zZ11YRGRn4HiIywT1D3ygi9waeXR+IqpYA43FCs4l7ye7vIjJbRIqA34nISSKywJ1tfSciF1bYTKKIfOSe6X8mIqlVGRdV3aaqf8P5EH9ib70ikiwi77r7sl5E/nSATXzu/jfXnSV0FZGWIvJJhVlMXFXqcWtaqap79r4EwoFmFfu5YX4lztjt1QJ4T1XzVTUPmA60OdxtuzoDb6pqjqr6VXW1qr5Toc+TwANVmUGq6hLgO6CtiCS6M9dcEckWkYVVOVZM9bCBN8HWFWgAvH+Ifo8DJwAZwPFACnBfwPJjcc6mU4CrgRdFJN5d9ry7LB04C+fDcySHICINgBHApoCZwOXAI0AM8CUwE/gQaArcBEwSkdYBmxkKPAQkAsuBSYd63wqmudtu7X5QzgRWuPvZG7hFRPpUst6Z7n/j3FnNEkCAx3BmFyfhBMCYgP19SUReOlgx7of3Lpx9XwBkVtKth1vzuwFtLwL9RSTe/bkMAuYcwbYBvgAeEZGRItLqAH2mAfk4P7+D7Y+ISDecsPwGuA34GWf2fAwwGicwTTCoqn3ZV9C+gGHAtgpti4FcYCfOB60ARUDLgD5dgfXu9z3dvmEBy7cDXYBQYDdwcsCyP+BcfqqsnjFu/1x3G58AHd1lrwMTAvr2ALYBIQFtbwFjAvpPDlgWDZQCzSp53zScD8KwCu0N3fZuwOnAfyssvwf4V0DtEw+2vQrrXgR8cwQ/s3CgLzDqAMv/CbxeoS0ZmA/43a+PgIjD3bbbJxInOJYCJTiX9/oGLFeck5jzgY04l4kf3ltTwNjkAjnAD8Cf3GUP4pwsHR/s3w37UptBmaD7BecyWNmlGFU9Q1Xj3GUhOGezjYCl7qWXXGCu2162Hd13iQjgV5xASMT50NsYsGwjzgzkQKaoapyqNlXVXqq6NGDZpoDvk3FmV/6DbLusv6oWAtnuelW1d1vZQCqQvHcM3HEYjXOmf0gicoyITBaRzSKSD0zEGZ/DoqolqjoHOLfiJU0RaQQMpvzlPYApwH9wZp6xwFr3/au87YA+O9W559URaOJue6qIJFToNxtnNvSHA+xKoqrGq+pJqvqc2/YUTuB9KCLr7OGJ4LKAMsG2BCgGDvjUFrADZ4bUxg2OOFVtrM4N7kPZgXOWHXjvpzmw+QjrDbzcswVoVuEeRcVtl91HEZFoIMFdr6oG4MzkfsQJu/UBYxCnqjGqev4h6tzrUbe9nToPIwzDmZ0eqTCcBx4q1puNc4kuUAbwsqoWuUH9D5wZzuFsez+qmo+zX1E497kq+jNOiDc61Lbc7RWo6m2qmo7z8MUoEeldlXXNb88CygSVquYCDwAvicjFIhIjIiEikoHzoYM7Q3kFeEZEmgKISMoB7r1U3H4pzhn2I+62U4FRVHL2fgS+xJmp3Ski4SLSE7gAmBzQ53wR6S4iETj3or5Q1U37bakCd7ZzI3A/cI87Bl8BBSJyl4hEikioiLQVkc6VbCIL51JaekBbDFAI5IlICnBHVXdURE4Ukb7u+4aLyDCcy6+fVeg6HOcyaMWA/Bq4xl0/EudJvZWHue29tfxFRDqLSISINMR5AjAXJ8TLUdUFwCq3rqrsZ38ROd590CMP55Ks/xCrmWpiAWWCTlWfxAmNO4H/uV8vA3fh3I/C/X4N8IV7eWo+0Hr/rVXqJpx7WOtwHht/E3jtN6h7N04g9cWZqb0EXKmqqwO6vYkTMtlAR5xZy8Hkuk8Ifoszwxisqq+571cK9MeZjax33/NVnAdAKtb2K87DHP92Lwd2wTkROBXng/cDnAcJyojzZOQ/DlCX4Nzj2o4TfjcDQ1R1WcD6KUAvYEIl61+Fc+/nZ5wZZjr7QuOQ2664e8C/3P3fApwD9HNnZpW5F2fmWhWtcI6tQpzZ/Uuq+mkV1zW/Mdn/RMcYY4wJPptBGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ40m/2l47rqobRjTSmyX4PSZmjFVqjf6S8fvHbg0/V4diEY4NdQp216pvlO1TVV7HdAuoQYpo0ZuDoq4JdRt0TG4w/Ul5PFJUEu4I66e4htwW7hDorPTZhY2XtdonPGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeFJYsAswVbPpu7UsnvIR6ldO7HYKGeedUWm/dctWM3/cNAbcMxJfahLb129h4aTZAKhCx/49aNGhdU2W7mmbVvzE4glzUb+fE393KhkX9qi037qvvmf+s1MY8PC1+NJT2FXwKx/9bQpZazdzwpkZdB/Zr4Yr975Nq9ayePKHzjHbI4OMvgc4ZpeuZv4/3mXAn0fiS0tm+/rNLJzgHrNAxwt60OLUE2uwcm/77KP5PHjXaPylpVwy/AquH3VLueWvvvAiU8a/QWhYGAmJiTz54vOkNG/G9yu/5S+33kZhQQEhoaH83+2j6D9oYHB2oooOGVAiosAkVR3mvg4DtgJfqmr/I3lTERmtqo8eyboB27gduAbYBZQAz6vqBBFZACQBO4EGwDOqOq7CujOAdFVtezQ11BS/38+it+bR7+bLiIqPZfpj/yK1fSvik33l+u3eVcyqT76maYvksraEFB8D7rmKkNAQfs0r5J2HXyW1fStCQm3y7Pf7WfSv2fS75wqimsQy/d5XSD21NfHHNS3Xb/fOYlbN/YKmx6eUtYWGh9H54t+R/fN2sjdtr+nSPc/v97Pozbn0u/Vy55h95DVSTznAMfvxV+WP2eSmDLj3aueYzS3gnQdfJfWUE+yYBUpLS7n/tjuZ8P40jk1J5qKevTn7/PNodeK+AG/Tvj3vf/YJkY0aMfHV13j8vvt5/vXXaBgZydiX/06L41vyv61bufDMXpzZuzexcY2DuEcHV5WfeBHQVkQi3dfnAJuP8n1HH83KIvJHt47TVDUD6A1IQJehbns34AkRiQhYdyBQeDTvX9OyNmyhcdN4Yn3xhIaF0rLzyWxY+dN+/TJnfE5Gn66Ehu077wiLCC/7xd5TsqfcINV3WWs20/iYBGKPSSA0LIyWXduyYemP+/XLnPoJGRd0JzR837iGN4zg2BNTy7WZfbLWb6GxL6H8Mbv8P/v1y3zvMzLO61puHMMaBB6zpYgdtWVWZC4lNb0FzVukERERQf9BA/nogznl+nQ9sweRjRoB0KFzJ7Zt3gJAeqvjaXF8SwCOSUqiiS+RX3bsqNkdOExVPSWZDey9hnEZ8NbeBSJymogsEZFvRGSxiLR220eIyDQRmSsiP4nIk27740CkiCwXkUlu2ygRWeV+3eK2pYnIDyLyioh8JyIfBoTkaOB6Vc0HUNV8VR1fSd3ROAFb6m4zGhgFPFzF/faEopwCouJjy15HxcVQlFNQrs+O/26jMCef5u2O32/97es3M/WBcbzz0Ct0v7yvnYm6inLyiWoSMK4JsRRl55frs2P9Fgp/yad5hxNqurxarSi3gKiEmLLXUfGxFOVWOGY3bnWO2fat9lt/+7rNTL3vZd55YBzdh51nx6xr29atJB23byaflJzM/7ZsPWD/KRMmctY5Z+/XviJzKSW7d5Oa3qJa6vytVPWnPhm4VEQaAu2BLwOWrQZ6qGoH4D4g8NJdBjAEaAcMEZFmqno3sFNVM1R1qIh0BEYCpwNdgGtFpIO7fivgRVVtA+QCg0QkFohR1XUHqXeSiKwEfgQeUtVSt/0h4Gng14PtrIhcJyKZIpK5q/CgXT1B/cqSqfPpOqh3pcubtkhh8P3XMeDukSyfu5g9JXtquMLaSf1+lkycR9dh5wa7lDpH/cqSKfPpOnj/D0+ApukpDH7wDwz481Usn2PH7JF4b/IUvv3mG669+aZy7du3bWPUddfz5EsvEBLi7eCv0vUJVV0pImk4s6fZFRY3BsaLSCuce5rhAcs+VtU8ABH5HkgFNlVYvzswXVWL3H7TgB7ADGC9qi53+y0F0qq0V84lvkwR8QGLRWQuEA+0VNVb3X052P6OA8YB+FKTtIrvWW2i4mMoytl3Zl+UW0BU/L6z05LiYrK3ZDHzr5MA2JlfyLyXptLnhsH4UpPK+sUnJRLeMIKcLVnl2uurqPhYin4JGNfsfKIS9s2oSnbtJnvTdmY+9DoAO/MKmTf2Lfrcfhm+9JSKmzMBouJiKMreN2MqysknKi7gmN3lHrNjJwLu2L4wlT43DsaXtu9+VHxSIuENIsjZvL1ce311bFISW3/ed4dl65YtHJO8/+/yok8X8OLYp3lrziwaNGhQ1l6Qn8/Vgy/ltvv+TIfTOtdIzUfjcC6gzwDGAj2BJgHtDwGfquoA94N/QcCy4oDvSw/z/SpbP1JV80WkUETSDzGLQlWzRGQZzuysCdBJRDa4dTQVkQWq2vMwa6pxvtRk8rbnkL8jl6i4GNZ+/T29rv592fKIyIYMf/rWstczn55Il4t740tNIn9HLtHxsYSEhlDwSx65234hpol3b4rWJF/LZPK2/UL+9hyiEmJYu2QVvW4cVLY8olFDho+7q+z1zIf+RZeh51o4VYEvLZm87dnkZ+USFe8es9dcVLY8olFDhj8zquz1zKfeoMvg3vjSksnPyiU6oeIxG1fzO+FB7TueyoZ169i0YSPHJCcx691pPPvPcs+A8d2Kldx78yj+NW0qib59D6Xs3r2bPw69kgGXDuH8i35fcdOedDiB8RqQq6rfikjPgPbG7HtoYkQVt1UiIuGqWgIsBF53700JMAC44hDrPwa8KCJD3MCKBgaq6oTATiLSCOgAPKmqU4C/u+1pwKzaEE4AIaEhdBtyLnOem4zf76f1GaeQkOwjc8ZnJKYmkXbKge+PbFuziRXzljjX8EXoflkfGkY3qsHqvSskNJRuI85nzuNv4PcrrXt2IOG4pmRO/YTE9GTSOh780eY3//QMJTuLKd1Tysalqzn/7iv2ewKwvgoJDaHb5X2Y8+xb+NVP626nkJDiI/N995jNOMQxO2exc8yGCN2HnkfDGDtmAcLCwhjz1JMMH3Ax/tJSBl8xlBNOOolnHn6Udqd24Ozz+/LYX+6nqKiIG4ePBCD5uON45e03mT3tPb7+92Jys7N5903nMYKn/v4iJ7dvF8xdOihRPfgVLBEpVNXoCm09gdtVtb+IdAXG4zyM8AEwTFXTRGQE0ElVb3TXmQWMVdUFIvIEcCGwzL0PNQq4yt38q6r6bECItHXXvx2IVtUxIiLAHcDVOI+YlwBPq+rESh4zf6PiI+0Vt30wvtQkHTj6qkN1M4crNuLQfcyRKSoJdgV10t1Dbgt2CXVWemzCUlXtVLH9kAFV31lAVRMLqOpjAVUtLKCqz4ECytuPcBhjjKm3LKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJYcEuoFaQYBdQB+3xB7uCOuvey+4Mdgl10sOTngh2CfWOzaCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ4UFuwBTNZtWrWXxlI9Qv3Ji91PIOO+MSvutW7aa+S9PY8A9I/GlJbF9/RYWTpwNgAId+/egRYfWNVi5t21auYbFE+ehfj8nntWBjAu6V9pv3dc/MP/5qQwYcw2+9GR+XrWWr6Z8QumeUkLDQjn90rNJOblFDVfvXQs+ms+YO++itLSUS4dfyf/dNqrc8leef4G3Xp9AWFgYCYlNGPv3FzmueXO+W7mSP98yioL8AkJDQ7nxjtu48OJBQdoLb6pPnwWHDCgRUWCSqg5zX4cBW4EvVbX/kbypiIxW1UePZN2AbdwOXAPsAkqA51V1gogsAJKAnUAD4BlVHeeuM9ddFgYsBP5PVUuPpo6a4Pf7WfTWPPrdchlR8bFMf+xfpLZvRXyyr1y/3buKWfXx1zRtkVzWlpDiY8DoqwgJDeHXvELeeehVUtu3IiTUJs9+v59FE+bQ785hRCXEMv3+V0k9tTXxKRXGdWcxqz78kqYtU8raGkY3os+tlxIVH0P2z9uZ/dQkhv3t1preBU8qLS3l3lG3MWnGeySlpHDBmb/jnPPP54STTizr06Z9ez5YuIDIRo1445VXefTe+3hpwutERjbimXEv0+L4lmzbupV+3c/irLN70zguLmj74yX17bOgKpUVAW1FJNJ9fQ6w+Sjfd/TRrCwif3TrOE1VM4DegAR0Geq2dwOeEJEIt/0SVT0FaAv4gMFHU0dNyVq/hcZN44n1xRMaFkrLTiezYcVP+/XLfP9zMs7rSmj4vvOOsIjwsgNwT8mecoNU32Wt3eyMa1N3XLu0YcOyH/frl/nuAjL6nVFuXBPTkoiKjwEgPsVH6e4SSkv21FjtXrY8cylp6emktmhBREQEF1w8kA8/+KBcnzPOOpPIRo0A6HBaZ7Zu2QJAeqvjaXF8SwCOTUoi0ecje8cvNbsDHlbfPguqGp2zgX7u95cBb+1dICKnicgSEflGRBaLSGu3fYSITBORuSLyk4g86bY/DkSKyHIRmeS2jRKRVe7XLW5bmoj8ICKviMh3IvJhQEiOBq5X1XwAVc1X1fGV1B2NE7Cle/u57WFABM5M1/OKcguIio8tex0VH0NRbkG5Pjv+u43CnHyatzt+v/W3r9/M1DHjeOfBV+g+tK+nz5hqUlFOAVFNGpe9jkqIpSinwrhu2Ephdh7NM0444HbWf/0DialJ5T4M6rNtW7aQfNy+2WZSSgr/27L1gP3fHv8GvzvnnP3al2cupWT3blLT7dLpXvXts6Cq1U0GLhWRhkB74MuAZauBHqraAbgPCLx0lwEMAdoBQ0SkmareDexU1QxVHSoiHYGRwOlAF+BaEengrt8KeFFV2wC5wCARiQViVHXdQeqdJCIrgR+BhwIv44nIPGA7UAC8U8X99zT1K0umzqfrxb0rXd60RQqDx1zHgHtGsnzuYvbYmX6VqF9Z8uaHdL3s3AP2yf55O19O+ZgeI/sdsI85sGmT32blN9/wh1v+VK79f9u2ccu11zH2Hy8REuLtD1EvqWufBVU65VPVlSKShjN7ml1hcWNgvIi0wpmRhAcs+1hV8wBE5HsgFdhUYf3uwHRVLXL7TQN6ADOA9aq63O23FEir0l45l/gyRcQHLBaRuaq60d2XPm7QTgJ6AR9VXFlErgOuA4hOiK24uMZFxcVQlJNf9roop4CouJiy1yXFxWRvzmLmXycBsDOvkHkvTaXPDYPxpSWV9YtPSiS8QQQ5m7PKtddXUfExFP2SV/a6KDu/7LIdQMmuYrJ/3s7Mx5zJ+c68QuY9O5k+t1yKLz2Zwux8PvrbFH533e+JPSahxuv3qmOTk9ny8767AFs3b+aY5P2Pt4WffsoLT45lytzZNGjQoKy9ID+fkYMGc8d9f+HU0zrXSM21RX37LDicaxIzgLFAT6BJQPtDwKeqOsANsQUBy4oDvi89zPerbP1IVc0XkUIRST/ELApVzRKRZTizs40B7btE5H3g91QSUO5DFeMAfKlJQb8M6EtLJm97Dvk7comKi2Ft5vf0uvr3ZcsjIhsy/K/7btDPfHoiXQb1xpeWRP6OXKLjYwkJDaHglzxyt/1CTGLjyt6m3vGlp5D3v2zys3KIio9l7Rff0ev6AWXLIxo1ZPhLd5S9nvnoeLpceg6+9GSKi3Yx9+m3OO2S3hx7QvNglO9Zp3Q8lfVr1/LfDRs4NjmZme9M47nXXi3XZ9WKFdzzp1t4Y/o0Epvuu8G/e/durr1sKAMvv4x+Ay6q4cq9r759FhxOYLwG5KrqtyLSM6C9MfsemhhRxW2ViEi4qpbgPE33untvSoABwBWHWP8x4EURGeIGVjQwUFUnBHYSkUZAB+BJt0+Mqm51n0Ts576354WEhtDt0nOZ87fJ+P1+Wnc7hYRkH5kzPiMxNYm0Uw58f2Tbmk2smLvEudYsQvfL+9AwulENVu9dIaEhdLuyL3OenIRfldZnZpBwXFMy3/2UxBbJpJ164Edwv5v/Ffn/y2bZ+5+z7P3PATj/zmFExkbVVPmeFRYWxkNPj+WKiwZSWlrKkCuG0frkk3j6oUdod2oHzu13Po/8+S/8WljE9VcMByC52XG8NmUys6ZN56t/LyY3O4d3Jr4JwNMvv0Sb9u2DuUueUd8+C0T14BMEESlU1egKbT2B21W1v4h0BcbjPIzwATBMVdNEZATQSVVvdNeZBYxV1QUi8gRwIbDMvQ81CrjK3fyrqvqsOxubpapt3fVvB6JVdYyICHAHcDXOI+YlwNOqOrGSx8zfUNVHReQYYJbbFgJ8Ctyqqge9COtLTdKBf77qYF3MkYi0Bwqqy70Dbwt2CXXSw5OeCHYJdda4Pzy6VFU7VWw/ZEDVdxZQ1cQCqtpYQFUPC6jqc6CAssdjjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCeFBbsAz1OgxB/sKuqcJ0b8Jdgl1Fl3TRgT7BLqpij7uKxpNoMyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBVUts+n4dbz84jslj/sHyD5ccsN+6b1Yz7sbHydq4tVx7YXYer416mhXzv6zuUmuV+fPm0bltG0496USeeerJ/Zb/e+FCzjq9M4mNGvL+tHfLLbvvnrvpmnEKp7dvx1233oKq1lTZtcKmVWt5+96/M3n0Syyfs/iA/dYtXc24ax8ha8MWALav38y7D7zCuw+8wjsPvML6ZatrquRaYdOKNbx9+4tMHvU8y2csOmC/dV/9wLhhD5K1zhnXXQW/MvOR8bx29WMsGj+npso9KocMKBFREZkY8DpMRLJEZNaRvqmIjD7SdQO2cbuIrBaR5SLytYhc6bYvEJEf3fYfROQ6t72RiHzgrvOdiDx+tDXUFL/fz6IpH9L3hksYfO+1rFn6PTlbd+zXb/euYlYtyKRpWvJ+y5ZM+4RmbdJrotxao7S0lDtu/hNTZ8zkixUrefftyaz+4ftyfZo1a8aLr/6Tiy+9tFz7l0sW8+WSxSxauozF3yznm6WZ/Pvzz2uyfE/z+/0senMufW++lMEP/oE1X31Hzpas/frt3lXMqo+/ommLfcdsQnJTBtx7NYPuv5bzb76UhRPn4C/112T5nuX3+1k0fg5977ycwU/ewJovviNncyXjurOYVfO+pGnLlLK20PAwOl/8O7pcfk5NlnxUqjKDKgLaikik+/ocYPNRvu9RBZSI/NGt4zRVzQB6AxLQZajb3g14QkQi3Paxqnoi0AHoJiJ9j6aOmpK1YSuNE+OJTYwjNCyUlqeezIaVP+3XL3PWQjLO6UJoWGi59g0r/kNMk8bEH5tYUyXXCku//or0li1JS08nIiKCgZcMYfbMmeX6NE9Lo2279oSElP9VERGKd+1i9+7dFBcXU1JSgq9p05os39Oy1m+hsS+BWF+8c8x2PpkNy/+zX7/M9z4j47yuhIaHlbWFNQgnJNQZ7z0lpUi5X+36LWvtZhofE09sU3dcu7Rhw9If9+uX+c4CMvqfUW5cwxtGcGzr5uXavK6ql/hmA/3c7y8D3tq7QEROE5ElIvKNiCwWkdZu+wgRmSYic0XkJxF50m1/HIh0ZziT3LZRIrLK/brFbUtzZ0CvuDOeDwNCcjRwvarmA6hqvqqOr6TuaJyALVXVX1X1U7f/bmAZcFwV9z+oivIKiIqPKXsdFR9DUV5BuT47Nm2jMCef5m2PL9deUryb5R99Qcfzu9dIrbXJ1i1bSGm27xBITklh6+aqnXud1qUrPc7qyYmpzTgxtRm9zjmX1iedVF2l1jpFuQVEJQQes7EU5VY4ZjdudY7Z9q32W3/7us1Mve9l3nlgHN2HnVcWWPVdUU4BUQmNy15HJcRSlFNhXNdvpTA7j+YdTqjp8n5zVf2pTwYuFZGGQHsg8EbGaqCHqnYA7gMeDViWAQwB2gFDRKSZqt4N7FTVDFUdKiIdgZHA6UAX4FoR6eCu3wp4UVXbALnAIBGJBWJUdd1B6p0kIiuBH4GHVLU0cKGIxAEXAB9Xcf89Tf3Kknc/puvAXvstW/rBItr16kx4g4hK1jRHat2aNfy4ejXfrdvA9+s3snDBpyxedOD7AaY89StLpsyn6+CzK13eND2FwQ/+gQF/vorlcxazp2RPDVdYO6lfWTLpQ7pefm6wS/lNVGmup6orRSQNZ/Y0u8LixsB4EWkFKBAesOxjVc0DEJHvgVRgU4X1uwPTVbXI7TcN6AHMANar6nK331IgrUp75VziyxQRH7BYROaq6kZ3+2E4M8DnDhRy7n2r6wCi42Or+JbVJ6pxTLmzpKKcAqIa7zs7LSkuJnvrDmb+7U0AduYXMe/ld+nzh0Fs37iFdctX8+V7n7J7ZzEiQmh4GG3P6ljj++E1ScnJbN70c9nrLZs3k5SScpA19pn1/vt0Pv10oqOjATi7z3l8/cUXnNHdZqoAUXExFGUHHrP5RMUFHLO7isneksXMsc7t7Z15hcx7YSp9bhyML+AeanxSIuENIsjZvL1ce30VFR9DUXZe2eui7PxyV1dKdhWT/fN2Zj7iXFDamVfIvL9Ops+oS/Gl177xO5yLkTOAsUBPoElA+0PAp6o6wA2xBQHLigO+Lz3M96ts/UhVzReRQhFJP8QsClXNEpFlOLOzjW7zOOAnVX32IOuNc/vha54U9EezfKlJ5GVlk78jl6i4GNYu+55eIy4sWx4R2ZDhT9xc9nrms5PoMqAXvtQkLrx1WFl75gcLCW8QYeHkOrVTZ9auWcPG9etJSklh2pS3eWXCG1Va97jmzZjwz39y6513oar8+/PP+eNNf6rmimsPX1oyeduzyc/KJSo+hrVff0+vay4qWx7RqCHDnxlV9nrmU2/QZXBvfGnJ5GflEp0QS0hoCAW/5JG77RdimsTV/E54kC89hbxt2eRvzyEqIZa1X3xHrxsGlC2PaNSQ4f+4o+z1zIfH0+Xyc2plOMHhBcZrQK6qfisiPQPaG7PvoYkRVdxWiYiEq2oJsBB43b03JcAA4IpDrP8Y8KKIDHEDKxoYqKoTAjuJSCOcByL23v962K33mirW6QkhoSF0u+Rc5rz4Nn5VWndpT0KSj8xZn5PYPIm0Sq7hm0MLCwvjyWf/xqD+/SgtLWXoiBGcdHIbHn1gDBmnduT8Cy5gWebXXHHJYHJzcpj7wQc8/uCDLFm+gt8PHMTnn35Kt1M7ICL0Pvdc+vbvH+xd8oyQ0BC6Xd6HOc++hV/9tO52CgkpPjLf/4zE1CTSMg58f2Tbmk2smLPYue8UInQfeh4NYxrVYPXeFRIaQrfhfZnz5CT8fqX1WRkkHNeUzHc+JbFFMmkdWx90/Tdv+RslO4sp3VPKxszVnH/3MOJTfDVU/eGTQ/3bDREpVNXoCm09gdtVtb+IdAXG4zyM8AEwTFXTRGQE0ElVb3TXmYXzFN0CEXkCuBBY5t6HGgVc5W7+VVV91p2NzVLVtu76twPRqjpGRAS4A7gaKHG/nlbViSKyAEgCdgINgDdU9VEROQ7n8uJq9s3MXlDVVw+2/77mSTrwrhEHHSNz+J649qFgl1Bn3TVhTLBLqJsiQw/dxxyRccMeXKqqnSq2HzKg6jsLqOphAVV9LKCqiQVUtTlQQNmzm8YYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTwoJdgNc1b5rCC9c/Euwy6pwbn7sn2CXUXZH2a10tsouDXUG9YzMoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ1lAGWOM8SQLKGOMMZ5kAWWMMcaTLKCMMcZ4kgWUMcYYT7KAMsYY40lhwS7AHNq8uXMZdeut+EtLGXn11dx5113lli/8/HNuGzWKb1euZOKbbzLo4ovLljUMD6dtu3YANGvWjOnvv1+jtXvdph/WsXj6x6gqJ57enoyzu1Tab92KH5n/+vsMuPUKfM2TKMjOY8rj/yTOlwBA09QkelzSpyZL97xN361l8ZSPUL9yYrdTyDjvjEr7rVu2mvnjpjHgnpH4UpPYvn4LCyfNBkAVOvbvQYsOrWuydE/b9P06Fk/7GPX7ObHrKWScc4BjdvmPzH/tPQbcfiW+5kll7YXZ+Ux59FU69u3GKb1Pr6myj8ghA0pEFJikqsPc12HAVuBLVe1/JG8qIqNV9dEjWTdgG7cD1wC7gBLgeVWdICILgCRgJ9AAeEZVx7nrPAJcCcSravTRvH9NKS0t5eabbmL2vHkcd9xxdD39dPpfcAEnn3xyWZ9mzZvz6muv8czTT++3fmRkJJnLltVkybWG3+9n0bvz6ffHS4iKi2H6MxNIbXs88ccmluu3e1cxqz5fStPUpHLtsU3iGHTHiBqsuPbw+/0semse/W6+jKj4WKY/9i9S27ciPtlXrt/uXcWs+uRrmrZILmtLSPEx4J6rCAkN4de8Qt55+FVS27ciJNQu+Pj9fhZN/Yh+/zfEOWbHjneO2aRKjtnPMvc7ZgGWTP+YZien11TJR6UqP/EioK2IRLqvzwE2H+X7jj6alUXkj24dp6lqBtAbkIAuQ932bsATIhLhts8ETjua965pX3/1FS1btiQ9PZ2IiAguGTKEmTNmlOuTlpZG+/btCQmxX+DDkfXfrTROjCM2MY7QsFBadjiJDavW7Ncvc84iMnqdTmiYXXCoqqwNW2jcNJ5YX7wztp1PZsPKn/brlznjczL6dC03tmER4WVhtKdkT7lf7Poua+NWGvsCjtlTT2LDt5WM6wcLyTi7C6Hh5Y/ZDSv/Q0yTuP1Owryqqp9os4F+7veXAW/tXSAip4nIEhH5RkQWi0hrt32EiEwTkbki8pOIPOm2Pw5EishyEZnkto0SkVXu1y1uW5qI/CAir4jIdyLyYUBIjgauV9V8AFXNV9XxldQdjROwpW6/L1R1a5VHxwM2b97Mcc2alb1OSUlhy+aqnx/s2rWLLqedRvczzuD9996rhgprr6LcQqLiYspeRzWOoSivoFyfHZu2UZhbQPM2LfdbvyA7j3fHvs7MF95k69pN1V5vbVKUU0BUfGzZ66i4GIpyKoztf7dRmJNP83bH77f+9vWbmfrAON556BW6X97XZk+uotwCouIqjGteYbk+BzpmS4p3s3z+l3Ts261Gav0tVPWUcDJwn4jMAtoDrwE93GWrgR6qukdEzgYeBQa5yzKADkAx8KOIPK+qd4vIje4MBxHpCIwETseZBX0pIp8BOUAr4DJVvVZEpgCDRGQGEKOq6w5S7yQRKXbXv0VVS6u4n7g1XQdcB9C8efPDWdVz1qxfT0pKCuvWraPP2WfTtl07Wrbc/8PW7E/9ypL3P6Xn5efvt6xRbBSX3/dHGkZFkrVpGx++Np3Bd11FRMMGQai09lG/smTqfHoOr/wuQdMWKQy+/zpytu5gweszada2JWHhNoM9FPUrS6Z/Qs+h/fZbtnTOItr17ER4g4hK1vSmKv3EVXWliKThzJ5mV1jcGBgvIq0ABcIDln2sqnkAIvI9kApUPNXsDkxX1SK33zSc8JsBrFfV5W6/pUBalfbKucSXKSI+YLGIzFXVjVVcF/ee1TiAjp06aVXXqw4pKSn8vGnfkG3evJnklJTDWh8gPT2dM886i+XffGMB5YqKi6Yod99ZfVFeAVGN982oSop3k71tBzNfcC4Y7CwoYt4/p9Hn6oH4mieVXZbyNTuW2CZx5G3PLnczuj6Lio+hKCe/7HVRbgFR8YFjW0z2lixm/nUSADvzC5n30lT63DAYX8B9k/ikRMIbRpCzJatce30VFRdDUW6FcW2873Z6SfFusrfuYObzbwKwM7+IeeOm0ee6gWzfsJV1y3/kyxkL2L2zGBEhNDyMtmd2rPH9qKrDOSWZAYwFegJNAtofAj5V1QFuiC0IWFYc8H3pYb5fZetHqmq+iBSKSPohZlGoapaILMOZnVU5oLykU+fOrFmzhvXuTGjK228zYeLEKq2bk5NDo0aNaNCgATt27GDJ4sXcdscd1Vxx7eFrlkReVg75v+QS1TiGtd/8QK9hF5Qtj4hswPCHbyp7PfOFt+hyYU98zZPYWfgrDRo1JCQkhPwdueTtyCGmSVwQ9sKbfKnJ5G3PIX9HLlFxMaz9+nt6Xf37suURkQ0Z/vStZa9nPj2RLhf3xpeaRP6OXKLjYwkJDaHglzxyt/1CTJPGwdgNz/E1r3DMLvuBXsMrHLOP/ans9czn3qTLRb/D1zyJC28ZWtaeOXsR4Q3CPR1OcHiB8RqQq6rfikjPgPbG7HtoYkQVt1UiIuGqWgIsBF53700JMAC44hDrPwa8KCJD3MCKBgaq6oTATiLSCOcS45NVrMtzwsLCePa55+jXty/+0lKGjxxJmzZtGHP//XTs2JELLryQzK+/ZvCgQeTk5PDBrFk8+MADrPj2W1b/8AM3XH89ISEh+P1+7rjzznJP/9V3IaEhdBt0NnNenorfr7Q+vR0JSYlkzllIYrNjSWvb6oDrbl27iaVzFhESGgoCPS4+l4ZRkQfsX9+EhIbQbci5zHluMn6/n9ZnnEJCso/MGZ+RmJpE2iknHHDdbWs2sWLeEue+kwjdL+tDw+hGNVi9d4WEhtDt4nOY89IU55jt0o6EJB+ZHywksfmxpLU78DFbG4nqwa9giUhhxUey3YC6XVX7i0hXYDzOwwgfAMNUNU1ERgCdVPVGd51ZwFhVXSAiTwAXAstUdaiIjAKucjf/qqo+687GZqlqW3f924FoVR0jIgLcAVyN84h5CfC0qk6s5DHzN/Y+0u4+qHE5kAxscd9rzMH2v2OnTvrFV18ddIzM4bvxuXuCXULdFWn3aqrF7sO6lW0Ow7g/PbFUVTtVbD9kQNV3FlDVwwKqGllAVQ8LqGpzoICyZzeNMcZ4kgWUMcYYT7KAMsYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxniSBZQxxhhPsoAyxhjjSRZQxhhjPMkCyhhjjCdZQBljjPEkCyhjjDGeZAFljDHGkyygjDHGeJIFlDHGGE+ygDLGGONJFlDGGGM8yQLKGGOMJ4mqBrsGTxORLGBjsOuookRgR7CLqINsXKuPjW31qG3jmqqqvoqNFlB1iIhkqmqnYNdR19i4Vh8b2+pRV8bVLvEZY4zxJAsoY4wxnmQBVbeMC3YBdZSNa/Wxsa0edWJc7R6UMcYYT7IZlDHGGE+ygDLGGONJFlDGGGM8yQKqHhCRaBGJDnYddZmI2O/Sb0hEmohIfLDrqItEpLGIxAa7jqqwX6o6TkT6A28B80XkcrdNgltV3SAivUTkLgBV9VtI/TbcY/ZtYKaIXBXseuoSEbkYmAp8KCJXi0i7YNd0MGHBLsBUHxE5F3gcuBZIAB4VkW2q+klwK6v9RKQ38C7wtYjEqeo9e0NKVf3Brq+2EpF+wMPAdUAj4K8i8r6q/hLcymo/EUkB7gOGA42Bi4CWIhKvqp8Hs7YDsTO+OkpEIoDTgYdVdYmqfoDzbyO6B7eyOiMeeAy4EUgRkcfBZlJHwx23tsBdqvoV8COgwMMiMlREWga1wNqvIVAErFTVBcBLwK9AHxE5MZiFHYj9O6g6SEREVVVEmgObAQFKgauBs1X1UrdfqKqWBrHUWqfiDElEwoAM4CZgm6re5bZHq2phcKqsfQKO2RA35GNwZqiZwMfAUGCVqv41qIXWQnvH1v3+aSAXeEpVd4lIa+AWYImqTghelZWzM706JvBgBH4HJKvqHrdtNVDg9hsKDLL7UVXnjq3f/f5KEWmhqnuA5cALwLEicpeIjACuEpHQ4FVbe1Q4ZoeJSKqqFgDXqOpoVf0YeB04X0SiglZoLVQhnK4EfgJigctFJFJVfwSmAVeKSKMgllopC6g6JuBgHAxcgzNz2qsA+NU9UO/FOSO1KXQVVTK2u932PcA3wGjgcuAZYIHNTqumknHd47b/N6DbsUAhYPf3DkPA2F4CXAHMwzmhagPcLSLhOJerf8WDY2sPSdRBIpKBc8lpqqpucS9D+YEI4AacA3SAqq4OWpG1VMDYvqOqm92xLVXVPSLSB+dhlG6q+n0w66xtKhnXssvPInINcD0wQlV3BrHMWskd2xuB6aq6XkT+B2wDLsa5fBoFXKuqu4JXZeVsBlUHVHKZbiewCucSXoZ7ic8PrAVm4PyiWzhVQRXHdu8sdBfQ18Lp0KowrqXiOAE4ERiuqt/WeKG10EHG9iIROVVVf1XVj1X1epwZf29VXVbjhVaBPSRRy1W4xtwTZ1a8N3yuAI4D/qGq37r3RCLt5n3VHM7YBqXAWqqK4/qyqq50P2wjVLU4GLXWNocztkEp8DDZDKqWCzgY/w/nRv1g4EOca8yzgS3AHSLSTlVLLZyq7jDGtk3QiqyFqjiut4tIW3VYOFXRYYxtrThmLaBqOfcySBpwCc7lpT8AY4A73S4TgJVAVlAKrMUOY2ztH5EehsMY1x1BKbAWq2vHrAVULRR4jdk9Y9rifrUQkXBVnQJ8AvxJVTcBz6nqtuBUW7vY2FYPG9fqU5fH1gKqlqlwjTlNRJJVdTewCTgL5xozOP9ANx/AXW4Owca2eti4Vp+6Prb2kEQtUuFgHIXzlyHWAZ/hXG/+O86j5AqcjD35VGU2ttXDxrX61IextX8HVYsEHIynA6cC/YEGwCQgVFVHikgH4ATgXlXdEKxaaxsb2+ph41p96sPY2gyqFnGvNbcH3sD5G2XXuv9eJBXnT+gvUdWbg1ljbWVjWz1sXKtPfRhbuwflcRVvgKrqCmAs0Aro4t4E3QgMAdqLyDGV/EM9Uwkb2+ph41p96tvY2gyqlhDnj7u2ArYDE4F+wFXAg8AXqloiImHq/F04cxhsbKuHjWv1qS9jazMojxKRZBGJdL+/CefvlOUArXH+4OM8YDzwFNAJyv5oqTkEG9vqYeNaferr2NpDEh4kzv/58m5glYhMANKAm1X1S3f5aOBJVb1GRBrjPEJqqsDGtnrYuFaf+jy2NoPypi3AUpwp/FCcP1NyVsDyWbg/O1V9Ucv/bwnMwdnYVg8b1+pTb8fWZlAes/ffNojzv78+2f1aBtwkItmq+irQDkgTkTggb+/jpubgbGyrh41r9anvY2sPSXiQewP0dmAkzj++2wHEAYNwzpbOBIao6nfBqrG2srGtHjau1ac+j60FlAeJyINAgao+JSJ7/yeDXXGm+VOAQlW1P6R5BGxsq4eNa/Wpz2Nr96C8aRnQTUTaqOpuVX0WSAeaANl19WCsITa21cPGtfrU27G1e1DetADoDFwuIp8AkUAezl8hzg9mYXXAAmxsq8MCbFyrywLq6djaJT6PEpFkYKD7tQe4XWvJ/wXT62xsq4eNa/Wpr2NrAeVxIhKF83Oy/xPub8zGtnrYuFaf+ja2FlDGGGM8yR6SMMYY40kWUMYYYzzJAsoYY4wnWUAZY4zxJAsoY4wxnmQBZYwxxpMsoIwxxnjS/wMYtJuIipCZdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "vireoSNP.plot.heat_matrix(GPb_diff[idx0, :][:, idx1], \n",
    "                          GPb0_donor_ids[idx0], \n",
    "                          GPb1_donor_ids[idx1])\n",
    "plt.title(\"Geno Prob Delta: %d SNPs\" %(len(GPb0_var_ids_use)))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "TFProb",
   "language": "python",
   "name": "tfprob"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}