Nucleotide diversity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.spatial.distance\n",
    "import seaborn as sns\n",
    "import os.path\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "sns.set_style(\"white\")\n",
    "from IPython.core.pylabtools import figsize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load haplotypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Catalog ID</th>\n",
       "      <th>Cnt</th>\n",
       "      <th>CMHAMM10_0002</th>\n",
       "      <th>CMHAMM10_0005</th>\n",
       "      <th>CMHAMM10_0008</th>\n",
       "      <th>CMHAMM10_0011</th>\n",
       "      <th>CMHAMM10_0012</th>\n",
       "      <th>CMHAMM10_0014</th>\n",
       "      <th>CMHAMM10_0015</th>\n",
       "      <th>CMHAMM10_0016</th>\n",
       "      <th>...</th>\n",
       "      <th>CMUW10_0003</th>\n",
       "      <th>CMUW10_0004</th>\n",
       "      <th>CMUW10_0007</th>\n",
       "      <th>CMUW10_0008</th>\n",
       "      <th>CMUW10_0009</th>\n",
       "      <th>CMUW10_0010</th>\n",
       "      <th>CMUW10_0011</th>\n",
       "      <th>CMUW10_0012</th>\n",
       "      <th>CMUW10_0015</th>\n",
       "      <th>CMUW10_0016</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>194</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>185</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>-</td>\n",
       "      <td>-</td>\n",
       "      <td>A/G</td>\n",
       "      <td>...</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4</td>\n",
       "      <td>186</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>-</td>\n",
       "      <td>consensus</td>\n",
       "      <td>-</td>\n",
       "      <td>...</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>195</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>190</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>...</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 202 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Catalog ID  Cnt CMHAMM10_0002 CMHAMM10_0005 CMHAMM10_0008 CMHAMM10_0011  \\\n",
       "0           1  194         GG/TA         GG/TA         GG/TA         GG/TA   \n",
       "1           2  185           A/G             G             G             G   \n",
       "2           4  186     consensus     consensus     consensus     consensus   \n",
       "3           7  195   GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   \n",
       "4           9  190           A/C             C           A/C             C   \n",
       "\n",
       "  CMHAMM10_0012 CMHAMM10_0014 CMHAMM10_0015 CMHAMM10_0016     ...      \\\n",
       "0         GG/TA         GG/TA         GG/TA         GG/TA     ...       \n",
       "1             G             -             -           A/G     ...       \n",
       "2     consensus             -     consensus             -     ...       \n",
       "3         GA/TA         GA/TA      GA/GG/TA   GA/GG/TA/TG     ...       \n",
       "4           A/C             C             A             C     ...       \n",
       "\n",
       "  CMUW10_0003 CMUW10_0004 CMUW10_0007 CMUW10_0008 CMUW10_0009 CMUW10_0010  \\\n",
       "0       GG/TA          TA       GG/TA       GG/TA          TA          TA   \n",
       "1         A/G           G           G         A/G           G           G   \n",
       "2   consensus   consensus   consensus   consensus   consensus   consensus   \n",
       "3       GG/TA       GA/TA       GA/TA       GA/TA       GA/TA    GA/GG/TA   \n",
       "4           C           A           A         A/C           A           A   \n",
       "\n",
       "  CMUW10_0011 CMUW10_0012  CMUW10_0015 CMUW10_0016  \n",
       "0          TA       GG/TA        GG/TA       GG/TA  \n",
       "1           G         A/G            G           G  \n",
       "2   consensus   consensus    consensus   consensus  \n",
       "3       GG/TA       GA/TA  GA/GG/TA/TG       GA/TA  \n",
       "4           C         A/C            A           -  \n",
       "\n",
       "[5 rows x 202 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "haplotypes = pd.read_csv('/media/Shared/Data/chum/populations/pstacks/batch_4/batch_4.haplotypes.tsv', sep = '\\t')\n",
    "haplotypes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CMHAMM10_0002</th>\n",
       "      <th>CMHAMM10_0005</th>\n",
       "      <th>CMHAMM10_0008</th>\n",
       "      <th>CMHAMM10_0011</th>\n",
       "      <th>CMHAMM10_0012</th>\n",
       "      <th>CMHAMM10_0014</th>\n",
       "      <th>CMHAMM10_0015</th>\n",
       "      <th>CMHAMM10_0016</th>\n",
       "      <th>CMHAMM10_0017</th>\n",
       "      <th>CMHAMM10_0018</th>\n",
       "      <th>...</th>\n",
       "      <th>CMUW10_0003</th>\n",
       "      <th>CMUW10_0004</th>\n",
       "      <th>CMUW10_0007</th>\n",
       "      <th>CMUW10_0008</th>\n",
       "      <th>CMUW10_0009</th>\n",
       "      <th>CMUW10_0010</th>\n",
       "      <th>CMUW10_0011</th>\n",
       "      <th>CMUW10_0012</th>\n",
       "      <th>CMUW10_0015</th>\n",
       "      <th>CMUW10_0016</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Catalog ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>-</td>\n",
       "      <td>-</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>...</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>-</td>\n",
       "      <td>consensus</td>\n",
       "      <td>-</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>...</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "      <td>consensus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>...</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 200 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           CMHAMM10_0002 CMHAMM10_0005 CMHAMM10_0008 CMHAMM10_0011  \\\n",
       "Catalog ID                                                           \n",
       "1                  GG/TA         GG/TA         GG/TA         GG/TA   \n",
       "2                    A/G             G             G             G   \n",
       "4              consensus     consensus     consensus     consensus   \n",
       "7            GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   \n",
       "9                    A/C             C           A/C             C   \n",
       "\n",
       "           CMHAMM10_0012 CMHAMM10_0014 CMHAMM10_0015 CMHAMM10_0016  \\\n",
       "Catalog ID                                                           \n",
       "1                  GG/TA         GG/TA         GG/TA         GG/TA   \n",
       "2                      G             -             -           A/G   \n",
       "4              consensus             -     consensus             -   \n",
       "7                  GA/TA         GA/TA      GA/GG/TA   GA/GG/TA/TG   \n",
       "9                    A/C             C             A             C   \n",
       "\n",
       "           CMHAMM10_0017 CMHAMM10_0018     ...     CMUW10_0003 CMUW10_0004  \\\n",
       "Catalog ID                                 ...                               \n",
       "1                  GG/TA         GG/TA     ...           GG/TA          TA   \n",
       "2                      G             G     ...             A/G           G   \n",
       "4              consensus     consensus     ...       consensus   consensus   \n",
       "7                  GA/TA         GA/TA     ...           GG/TA       GA/TA   \n",
       "9                      C             C     ...               C           A   \n",
       "\n",
       "           CMUW10_0007 CMUW10_0008 CMUW10_0009 CMUW10_0010 CMUW10_0011  \\\n",
       "Catalog ID                                                               \n",
       "1                GG/TA       GG/TA          TA          TA          TA   \n",
       "2                    G         A/G           G           G           G   \n",
       "4            consensus   consensus   consensus   consensus   consensus   \n",
       "7                GA/TA       GA/TA       GA/TA    GA/GG/TA       GG/TA   \n",
       "9                    A         A/C           A           A           C   \n",
       "\n",
       "           CMUW10_0012  CMUW10_0015 CMUW10_0016  \n",
       "Catalog ID                                       \n",
       "1                GG/TA        GG/TA       GG/TA  \n",
       "2                  A/G            G           G  \n",
       "4            consensus    consensus   consensus  \n",
       "7                GA/TA  GA/GG/TA/TG       GA/TA  \n",
       "9                  A/C            A           -  \n",
       "\n",
       "[5 rows x 200 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "haplotypes = haplotypes.set_index('Catalog ID')\n",
    "haplotypes = haplotypes.drop('Cnt', axis =1 )\n",
    "haplotypes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exclude individuals excluded during genotype filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['CMHAMM10_0002',\n",
       " 'CMHAMM10_0005',\n",
       " 'CMHAMM10_0008',\n",
       " 'CMHAMM10_0011',\n",
       " 'CMHAMM10_0012',\n",
       " 'CMHAMM10_0014',\n",
       " 'CMHAMM10_0015',\n",
       " 'CMHAMM10_0016']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indlist = pd.read_csv(os.path.join('data','batch_4','pop_genotypes','complete.nosex'), sep = '\\t', header = None)\n",
    "indlist = indlist.drop(0, axis =1)\n",
    "inds_to_keep = indlist[1].values.tolist()\n",
    "inds_to_keep[:8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "haplotypes = haplotypes[inds_to_keep]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Replace:\n",
    "'-' is missing data\n",
    "\n",
    "'consensus' is scored as invariant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CMHAMM10_0002</th>\n",
       "      <th>CMHAMM10_0005</th>\n",
       "      <th>CMHAMM10_0008</th>\n",
       "      <th>CMHAMM10_0011</th>\n",
       "      <th>CMHAMM10_0012</th>\n",
       "      <th>CMHAMM10_0014</th>\n",
       "      <th>CMHAMM10_0015</th>\n",
       "      <th>CMHAMM10_0016</th>\n",
       "      <th>CMHAMM10_0017</th>\n",
       "      <th>CMHAMM10_0018</th>\n",
       "      <th>...</th>\n",
       "      <th>CMSTILL10_0067</th>\n",
       "      <th>CMSTILL10_0080</th>\n",
       "      <th>CMUW10_0001</th>\n",
       "      <th>CMUW10_0002</th>\n",
       "      <th>CMUW10_0003</th>\n",
       "      <th>CMUW10_0004</th>\n",
       "      <th>CMUW10_0007</th>\n",
       "      <th>CMUW10_0008</th>\n",
       "      <th>CMUW10_0009</th>\n",
       "      <th>CMUW10_0011</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Catalog ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>TA</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>...</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "      <td>A/G</td>\n",
       "      <td>G</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>...</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GA/GG/TA/TG</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>...</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GA/TA</td>\n",
       "      <td>GG/TA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>C</td>\n",
       "      <td>...</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A/C</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>A</td>\n",
       "      <td>A/C</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 174 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           CMHAMM10_0002 CMHAMM10_0005 CMHAMM10_0008 CMHAMM10_0011  \\\n",
       "Catalog ID                                                           \n",
       "1                  GG/TA         GG/TA         GG/TA         GG/TA   \n",
       "2                    A/G             G             G             G   \n",
       "4                      A             A             A             A   \n",
       "7            GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   GA/GG/TA/TG   \n",
       "9                    A/C             C           A/C             C   \n",
       "\n",
       "           CMHAMM10_0012 CMHAMM10_0014 CMHAMM10_0015 CMHAMM10_0016  \\\n",
       "Catalog ID                                                           \n",
       "1                  GG/TA         GG/TA         GG/TA         GG/TA   \n",
       "2                      G           NaN           NaN           A/G   \n",
       "4                      A           NaN             A           NaN   \n",
       "7                  GA/TA         GA/TA      GA/GG/TA   GA/GG/TA/TG   \n",
       "9                    A/C             C             A             C   \n",
       "\n",
       "           CMHAMM10_0017 CMHAMM10_0018     ...     CMSTILL10_0067  \\\n",
       "Catalog ID                                 ...                      \n",
       "1                  GG/TA         GG/TA     ...              GG/TA   \n",
       "2                      G             G     ...                A/G   \n",
       "4                      A             A     ...                  A   \n",
       "7                  GA/TA         GA/TA     ...              GG/TA   \n",
       "9                      C             C     ...                A/C   \n",
       "\n",
       "           CMSTILL10_0080 CMUW10_0001 CMUW10_0002 CMUW10_0003 CMUW10_0004  \\\n",
       "Catalog ID                                                                  \n",
       "1                   GG/TA          TA       GG/TA       GG/TA          TA   \n",
       "2                       G           G           G         A/G           G   \n",
       "4                       A           A           A           A           A   \n",
       "7                   GA/TA    GA/GG/TA       GA/TA       GG/TA       GA/TA   \n",
       "9                     A/C         A/C         A/C           C           A   \n",
       "\n",
       "           CMUW10_0007 CMUW10_0008 CMUW10_0009 CMUW10_0011  \n",
       "Catalog ID                                                  \n",
       "1                GG/TA       GG/TA          TA          TA  \n",
       "2                    G         A/G           G           G  \n",
       "4                    A           A           A           A  \n",
       "7                GA/TA       GA/TA       GA/TA       GG/TA  \n",
       "9                    A         A/C           A           C  \n",
       "\n",
       "[5 rows x 174 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "haplotypes = haplotypes.replace({'consensus':'A', '-':np.NaN})\n",
    "haplotypes = haplotypes.dropna(how = 'all')\n",
    "haplotypes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29131, 174)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "haplotypes.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## from each sample, generate two sequences per locus\n",
    "    - 1,2,3,4 are variable nuceotides (substituted in below)\n",
    "\n",
    "    - 9's are invariate positions\n",
    "\n",
    "    - 0's are missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "chromosomes = []\n",
    "for catid, haps in haplotypes.iterrows():\n",
    "    chrs = []\n",
    "    for hap_call in haps:\n",
    "        # get chr1       \n",
    "        try:\n",
    "            hap1 = hap_call.split('/')[0]\n",
    "            var_sites = len(hap1)\n",
    "            chrs.append(hap1 + '9'*(94-var_sites))\n",
    "        except AttributeError:\n",
    "            chrs.append('0'*94)\n",
    "        # get chr2\n",
    "        try:\n",
    "            hap2 = hap_call.split('/')[1]\n",
    "            chrs.append(hap2 + '9'*(94-var_sites))\n",
    "        except AttributeError:\n",
    "            chrs.append('0'*94)\n",
    "        except IndexError: # homozygous\n",
    "            hap2 = hap_call.split('/')[0]\n",
    "            chrs.append(hap2 + '9'*(94-var_sites))\n",
    "    chromosomes.append(chrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### maintain locus IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Catalog ID</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>4</th>\n",
       "      <th>7</th>\n",
       "      <th>9</th>\n",
       "      <th>13</th>\n",
       "      <th>14</th>\n",
       "      <th>15</th>\n",
       "      <th>16</th>\n",
       "      <th>21</th>\n",
       "      <th>...</th>\n",
       "      <th>86838</th>\n",
       "      <th>86853</th>\n",
       "      <th>86996</th>\n",
       "      <th>87091</th>\n",
       "      <th>87172</th>\n",
       "      <th>87322</th>\n",
       "      <th>87493</th>\n",
       "      <th>87988</th>\n",
       "      <th>87992</th>\n",
       "      <th>88006</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>GG99999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>GA99999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>CAGGCC9999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>T999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>ACA9999999999999999999999999999999999999999999...</td>\n",
       "      <td>...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TA99999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>GG99999999999999999999999999999999999999999999...</td>\n",
       "      <td>C999999999999999999999999999999999999999999999...</td>\n",
       "      <td>TGTCGA9999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>T999999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>ACA9999999999999999999999999999999999999999999...</td>\n",
       "      <td>...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>GG99999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>GA99999999999999999999999999999999999999999999...</td>\n",
       "      <td>C999999999999999999999999999999999999999999999...</td>\n",
       "      <td>TGTCGA9999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>T999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>ACA9999999999999999999999999999999999999999999...</td>\n",
       "      <td>...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA99999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>GG99999999999999999999999999999999999999999999...</td>\n",
       "      <td>C999999999999999999999999999999999999999999999...</td>\n",
       "      <td>TGTCGA9999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>T999999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>TGG9999999999999999999999999999999999999999999...</td>\n",
       "      <td>...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GG99999999999999999999999999999999999999999999...</td>\n",
       "      <td>G999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>GA99999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>CAGGCC9999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>T999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>ACA9999999999999999999999999999999999999999999...</td>\n",
       "      <td>...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "      <td>0000000000000000000000000000000000000000000000...</td>\n",
       "      <td>A999999999999999999999999999999999999999999999...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 29131 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "Catalog ID                                              1      \\\n",
       "0           GG99999999999999999999999999999999999999999999...   \n",
       "1           TA99999999999999999999999999999999999999999999...   \n",
       "2           GG99999999999999999999999999999999999999999999...   \n",
       "3           TA99999999999999999999999999999999999999999999...   \n",
       "4           GG99999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              2      \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           G999999999999999999999999999999999999999999999...   \n",
       "2           G999999999999999999999999999999999999999999999...   \n",
       "3           G999999999999999999999999999999999999999999999...   \n",
       "4           G999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              4      \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              7      \\\n",
       "0           GA99999999999999999999999999999999999999999999...   \n",
       "1           GG99999999999999999999999999999999999999999999...   \n",
       "2           GA99999999999999999999999999999999999999999999...   \n",
       "3           GG99999999999999999999999999999999999999999999...   \n",
       "4           GA99999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              9      \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           C999999999999999999999999999999999999999999999...   \n",
       "2           C999999999999999999999999999999999999999999999...   \n",
       "3           C999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              13     \\\n",
       "0           CAGGCC9999999999999999999999999999999999999999...   \n",
       "1           TGTCGA9999999999999999999999999999999999999999...   \n",
       "2           TGTCGA9999999999999999999999999999999999999999...   \n",
       "3           TGTCGA9999999999999999999999999999999999999999...   \n",
       "4           CAGGCC9999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              14     \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              15     \\\n",
       "0           T999999999999999999999999999999999999999999999...   \n",
       "1           T999999999999999999999999999999999999999999999...   \n",
       "2           T999999999999999999999999999999999999999999999...   \n",
       "3           T999999999999999999999999999999999999999999999...   \n",
       "4           T999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              16     \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           G999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           G999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              21     \\\n",
       "0           ACA9999999999999999999999999999999999999999999...   \n",
       "1           ACA9999999999999999999999999999999999999999999...   \n",
       "2           ACA9999999999999999999999999999999999999999999...   \n",
       "3           TGG9999999999999999999999999999999999999999999...   \n",
       "4           ACA9999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                        ...                          \\\n",
       "0                                 ...                           \n",
       "1                                 ...                           \n",
       "2                                 ...                           \n",
       "3                                 ...                           \n",
       "4                                 ...                           \n",
       "\n",
       "Catalog ID                                              86838  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           0000000000000000000000000000000000000000000000...   \n",
       "\n",
       "Catalog ID                                              86853  \\\n",
       "0           0000000000000000000000000000000000000000000000...   \n",
       "1           0000000000000000000000000000000000000000000000...   \n",
       "2           0000000000000000000000000000000000000000000000...   \n",
       "3           0000000000000000000000000000000000000000000000...   \n",
       "4           0000000000000000000000000000000000000000000000...   \n",
       "\n",
       "Catalog ID                                              86996  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              87091  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              87172  \\\n",
       "0           0000000000000000000000000000000000000000000000...   \n",
       "1           0000000000000000000000000000000000000000000000...   \n",
       "2           0000000000000000000000000000000000000000000000...   \n",
       "3           0000000000000000000000000000000000000000000000...   \n",
       "4           0000000000000000000000000000000000000000000000...   \n",
       "\n",
       "Catalog ID                                              87322  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           0000000000000000000000000000000000000000000000...   \n",
       "3           0000000000000000000000000000000000000000000000...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              87493  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           0000000000000000000000000000000000000000000000...   \n",
       "3           0000000000000000000000000000000000000000000000...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              87988  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           A999999999999999999999999999999999999999999999...   \n",
       "3           A999999999999999999999999999999999999999999999...   \n",
       "4           A999999999999999999999999999999999999999999999...   \n",
       "\n",
       "Catalog ID                                              87992  \\\n",
       "0           A999999999999999999999999999999999999999999999...   \n",
       "1           A999999999999999999999999999999999999999999999...   \n",
       "2           0000000000000000000000000000000000000000000000...   \n",
       "3           0000000000000000000000000000000000000000000000...   \n",
       "4           0000000000000000000000000000000000000000000000...   \n",
       "\n",
       "Catalog ID                                              88006  \n",
       "0           A999999999999999999999999999999999999999999999...  \n",
       "1           A999999999999999999999999999999999999999999999...  \n",
       "2           0000000000000000000000000000000000000000000000...  \n",
       "3           0000000000000000000000000000000000000000000000...  \n",
       "4           A999999999999999999999999999999999999999999999...  \n",
       "\n",
       "[5 rows x 29131 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hap_df = pd.DataFrame(chromosomes).T\n",
    "hap_df.columns = haplotypes.index\n",
    "hap_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate pi\n",
    "as the jaccard similary of sequences:\n",
    "\n",
    "\"Jaccard distance between the points. Given two vectors, u and v, the Jaccard distance is the proportion of those elements u[i] and v[i] that disagree where at least one of them is non-zero.\"\n",
    "\n",
    "[http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.spatial.distance.pdist.html]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Per locus I produce a distance matrix of pairwise nucleotide difference between chromosomes, some within population and some between.\n",
    "\n",
    "#### Different subsets of this distance matrix show within and between population comparisons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Still need one more dimension - break up each locus into sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Pi_of_locus = dict()\n",
    "pdist_of_locus = dict()\n",
    "for locus in hap_df.columns:\n",
    "    myseqs = []\n",
    "    for x in hap_df[locus].values.tolist():\n",
    "        try:\n",
    "            myseqs.append(list(x))\n",
    "        except TypeError:\n",
    "            myseqs.append([np.NaN]*94)\n",
    "    myseqs = pd.DataFrame(myseqs)\n",
    "    myseqs = myseqs.replace({'A':1, 'C':2, 'G':3, 'T':4})\n",
    "    pdist_mat = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(myseqs, 'jaccard'))\n",
    "    ## found values of 1 to be from comparing missing values to anything non-missing - remove\n",
    "    pdist_mat = pd.DataFrame(pdist_mat).replace({1:np.NaN})\n",
    "    pdist_of_locus[locus] = pdist_mat    \n",
    "    ## remove chomosomes with mostly missing values \n",
    "    cleaned_pdist_mat = pdist_mat.dropna(thresh = 50, axis = 0)\n",
    "    cleaned_pdist_mat = cleaned_pdist_mat.dropna(thresh = 50, axis = 1)\n",
    "    # mean calulated without zeros as self-compairisons    \n",
    "    Pi_of_locus[locus] = np.nanmean(scipy.spatial.distance.squareform(cleaned_pdist_mat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "for locus in [7]:\n",
    "    myseqs = []\n",
    "    for x in hap_df[locus].values.tolist():\n",
    "        try:\n",
    "            myseqs.append(list(x))\n",
    "        except TypeError:\n",
    "            myseqs.append([np.NaN]*94)\n",
    "    myseqs = pd.DataFrame(myseqs)\n",
    "    myseqs = myseqs.replace({'A':1, 'C':2, 'G':3, 'T':4})\n",
    "    pdist_mat = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(myseqs, 'jaccard'))\n",
    "    ## found values of 1 to be from comparing missing values to anything non-missing - remove\n",
    "    pdist_mat = pd.DataFrame(pdist_mat).replace({1:np.NaN})\n",
    "    pdist_of_locus[locus] = pdist_mat    \n",
    "    ## remove chomosomes with mostly missing values \n",
    "    cleaned_pdist_mat = pdist_mat.dropna(thresh = 50, axis = 0)\n",
    "    cleaned_pdist_mat = cleaned_pdist_mat.dropna(thresh = 50, axis = 1)\n",
    "    # mean calulated without zeros as self-compairisons    \n",
    "    Pi_of_locus[locus] = np.nanmean(scipy.spatial.distance.squareform(cleaned_pdist_mat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### merge with individual IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0008</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0012</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0              1\n",
       "0  1  CMHAMM10_0002\n",
       "1  1  CMHAMM10_0005\n",
       "2  1  CMHAMM10_0008\n",
       "3  1  CMHAMM10_0011\n",
       "4  1  CMHAMM10_0012"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indlist = pd.read_csv(os.path.join('data','batch_4','pop_genotypes','complete.nosex'), sep = '\\t', header = None)\n",
    "indlist.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pop</th>\n",
       "      <th>ind</th>\n",
       "      <th>chr</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0002</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0002</td>\n",
       "      <td>b</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0005</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0005</td>\n",
       "      <td>b</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>CMHAMM10_0008</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pop            ind chr\n",
       "0    1  CMHAMM10_0002   a\n",
       "1    1  CMHAMM10_0002   b\n",
       "2    1  CMHAMM10_0005   a\n",
       "3    1  CMHAMM10_0005   b\n",
       "4    1  CMHAMM10_0008   a"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chr_index = pd.concat([indlist, indlist]).sort_values([0,1])\n",
    "chr_index.columns = ['pop', 'ind']\n",
    "chr_index['chr'] = ['a','b']*174\n",
    "chr_index.shape\n",
    "chr_index = chr_index.reset_index(drop = True)\n",
    "chr_index.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29131"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(pdist_of_locus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### calculate pairwise nucleotide diversity between all chromosomes within each population\n",
    "chromosomes are not compared to themselves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_within_pi(pdist_df):\n",
    "    within_pi_of_pop = dict()\n",
    "    for xx in range(1,11):\n",
    "        within_pdist_vals = pdist_df.loc[xx,xx].values\n",
    "        np.fill_diagonal(within_pdist_vals,np.NaN)\n",
    "        within_pi_of_pop[xx] = np.nanmean(within_pdist_vals)\n",
    "    return(within_pi_of_pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate d_xy, the mean pairwise necelotide diversity between chromosomes from different populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_dxy(pdist_df):\n",
    "    between_pi_of_pops = dict()\n",
    "    for xx in range(1,11):\n",
    "        for yy in range(1,xx):\n",
    "            between_pdist_vals = pdist_df.loc[yy,xx]\n",
    "            between_pi_of_pops[yy,xx] = np.nanmean(between_pdist_vals)\n",
    "    return(between_pi_of_pops)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## apply to each locus' pairwise pi matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 192\n",
      "200 398\n",
      "300 619\n",
      "400 858\n",
      "500 1087\n",
      "600 1277\n",
      "700 1481\n",
      "800 1667\n",
      "900 1887\n",
      "1000 2095\n",
      "1100 2296\n",
      "1200 2492\n",
      "1300 2720\n",
      "1400 2921\n",
      "1500 3164\n",
      "1600 3346\n",
      "1700 3559\n",
      "1800 3804\n",
      "1900 4022\n",
      "2000 4222\n",
      "2100 4421\n",
      "2200 4623\n",
      "2300 4838\n",
      "2400 5062\n",
      "2500 5287\n",
      "2600 5506\n",
      "2700 5725\n",
      "2800 5928\n",
      "2900 6152\n",
      "3000 6377\n",
      "3100 6589\n",
      "3200 6787\n",
      "3300 6992\n",
      "3400 7219\n",
      "3500 7418\n",
      "3600 7626\n",
      "3700 7842\n",
      "3800 8090\n",
      "3900 8328\n",
      "4000 8537\n",
      "4100 8770\n",
      "4200 9002\n",
      "4300 9254\n",
      "4400 9480\n",
      "4500 9685\n",
      "4600 9893\n",
      "4700 10110\n",
      "4800 10364\n",
      "4900 10608\n",
      "5000 10821\n",
      "5100 11078\n",
      "5200 11293\n",
      "5300 11528\n",
      "5400 11768\n",
      "5500 11982\n",
      "5600 12225\n",
      "5700 12454\n",
      "5800 12653\n",
      "5900 12892\n",
      "6000 13115\n",
      "6100 13320\n",
      "6200 13529\n",
      "6300 13773\n",
      "6400 14014\n",
      "6500 14231\n",
      "6600 14451\n",
      "6700 14680\n",
      "6800 14913\n",
      "6900 15137\n",
      "7000 15351\n",
      "7100 15571\n",
      "7200 15783\n",
      "7300 16027\n",
      "7400 16280\n",
      "7500 16490\n",
      "7600 16701\n",
      "7700 16950\n",
      "7800 17197\n",
      "7900 17429\n",
      "8000 17725\n",
      "8100 17966\n",
      "8200 18209\n",
      "8300 18423\n",
      "8400 18636\n",
      "8500 18866\n",
      "8600 19091\n",
      "8700 19337\n",
      "8800 19585\n",
      "8900 19812\n",
      "9000 20040\n",
      "9100 20284\n",
      "9200 20514\n",
      "9300 20771\n",
      "9400 21017\n",
      "9500 21246\n",
      "9600 21470\n",
      "9700 21716\n",
      "9800 21950\n",
      "9900 22203\n",
      "10000 22426\n",
      "10100 22628\n",
      "10200 22842\n",
      "10300 23056\n",
      "10400 23282\n",
      "10500 23510\n",
      "10600 23725\n",
      "10700 24007\n",
      "10800 24241\n",
      "10900 24467\n",
      "11000 24688\n",
      "11100 24923\n",
      "11200 25154\n",
      "11300 25405\n",
      "11400 25656\n",
      "11500 25847\n",
      "11600 26080\n",
      "11700 26295\n",
      "11800 26506\n",
      "11900 26728\n",
      "12000 26944\n",
      "12100 27161\n",
      "12200 27377\n",
      "12300 27574\n",
      "12400 27801\n",
      "12500 28026\n",
      "12600 28252\n",
      "12700 28466\n",
      "12800 28669\n",
      "12900 28853\n",
      "13000 29055\n",
      "13100 29282\n",
      "13200 29526\n",
      "13300 29746\n",
      "13400 29951\n",
      "13500 30159\n",
      "13600 30373\n",
      "13700 30577\n",
      "13800 30791\n",
      "13900 31011\n",
      "14000 31204\n",
      "14100 31411\n",
      "14200 31625\n",
      "14300 31860\n",
      "14400 32090\n",
      "14500 32325\n",
      "14600 32538\n",
      "14700 32729\n",
      "14800 32923\n",
      "14900 33127\n",
      "15000 33306\n",
      "15100 33518\n",
      "15200 33741\n",
      "15300 33941\n",
      "15400 34157\n",
      "15500 34360\n",
      "15600 34571\n",
      "15700 34765\n",
      "15800 34977\n",
      "15900 35165\n",
      "16000 35339\n",
      "16100 35556\n",
      "16200 35771\n",
      "16300 35975\n",
      "16400 36189\n",
      "16500 36405\n",
      "16600 36622\n",
      "16700 36830\n",
      "16800 37035\n",
      "16900 37237\n",
      "17000 37429\n",
      "17100 37665\n",
      "17200 37846\n",
      "17300 38036\n",
      "17400 38243\n",
      "17500 38448\n",
      "17600 38645\n",
      "17700 38839\n",
      "17800 39054\n",
      "17900 39247\n",
      "18000 39458\n",
      "18100 39678\n",
      "18200 39888\n",
      "18300 40078\n",
      "18400 40273\n",
      "18500 40463\n",
      "18600 40666\n",
      "18700 40851\n",
      "18800 41036\n",
      "18900 41253\n",
      "19000 41466\n",
      "19100 41693\n",
      "19200 41919\n",
      "19300 42128\n",
      "19400 42334\n",
      "19500 42550\n",
      "19600 42742\n",
      "19700 42929\n",
      "19800 43121\n",
      "19900 43359\n",
      "20000 43573\n",
      "20100 43754\n",
      "20200 43966\n",
      "20300 44173\n",
      "20400 44399\n",
      "20500 44601\n",
      "20600 44817\n",
      "20700 45015\n",
      "20800 45230\n",
      "20900 45455\n",
      "21000 45645\n",
      "21100 45840\n",
      "21200 46054\n",
      "21300 46299\n",
      "21400 46492\n",
      "21500 46695\n",
      "21600 46874\n",
      "21700 47068\n",
      "21800 47276\n",
      "21900 47487\n",
      "22000 47670\n",
      "22100 47889\n",
      "22200 48086\n",
      "22300 48310\n",
      "22400 48519\n",
      "22500 48726\n",
      "22600 48958\n",
      "22700 49170\n",
      "22800 49395\n",
      "22900 49595\n",
      "23000 49808\n",
      "23100 50060\n",
      "23200 50281\n",
      "23300 50514\n",
      "23400 50733\n",
      "23500 50973\n",
      "23600 51191\n",
      "23700 51412\n",
      "23800 51647\n",
      "23900 51859\n",
      "24000 52084\n",
      "24100 52291\n",
      "24200 52491\n",
      "24300 52698\n",
      "24400 52928\n",
      "24500 53143\n",
      "24600 53365\n",
      "24700 53588\n",
      "24800 53814\n",
      "24900 54048\n",
      "25000 54245\n",
      "25100 54465\n",
      "25200 54697\n",
      "25300 54916\n",
      "25400 55136\n",
      "25500 55350\n",
      "25600 55584\n",
      "25700 55818\n",
      "25800 56033\n",
      "25900 56241\n",
      "26000 56471\n",
      "26100 56679\n",
      "26200 56903\n",
      "26300 57128\n",
      "26400 57375\n",
      "26500 57604\n",
      "26600 57807\n",
      "26700 58046\n",
      "26800 58269\n",
      "26900 58635\n",
      "27000 59559\n",
      "27100 60401\n",
      "27200 61165\n",
      "27300 62004\n",
      "27400 62776\n",
      "27500 63537\n",
      "27600 64234\n",
      "27700 65123\n",
      "27800 65926\n",
      "27900 66523\n",
      "28000 66878\n",
      "28100 67279\n",
      "28200 67702\n",
      "28300 68066\n",
      "28400 68429\n",
      "28500 69072\n",
      "28600 71150\n",
      "28700 73487\n",
      "28800 75625\n",
      "28900 77946\n",
      "29000 80300\n",
      "29100 84528\n"
     ]
    }
   ],
   "source": [
    "count = 0\n",
    "with open(os.path.join('results','batch_4','pop_analysis','within_pi.txt'), 'w') as WITHIN:\n",
    "    with open(os.path.join('results','batch_4','pop_analysis','dxy.txt'), 'w') as DXY:\n",
    "        for locus,pdist in pdist_of_locus.items():\n",
    "            count +=1 \n",
    "            if count %100 == 0:\n",
    "                print count, locus\n",
    "            pdist = pd.concat([chr_index, pdist], axis = 1)\n",
    "            pdist = pdist.set_index(['pop', 'ind', 'chr'])\n",
    "            pdist.columns = pdist.index\n",
    "            for pop,within_pi in get_within_pi(pdist).items():\n",
    "                WITHIN.write('{}\\t{}\\t{}\\n'.format(locus,pop,within_pi))\n",
    "            for pops,dxy in get_dxy(pdist).items():\n",
    "                pop1 = pops[0]\n",
    "                pop2 = pops[1]\n",
    "                DXY.write('{}\\t{}\\t{}\\t{}\\n'.format(locus,pop1, pop2, dxy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assert False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 308,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "locus\n",
       "1        0.010642\n",
       "2        0.002187\n",
       "4        0.000000\n",
       "7        0.008766\n",
       "9        0.005354\n",
       "13       0.031836\n",
       "14       0.000000\n",
       "15       0.004284\n",
       "16       0.004655\n",
       "21       0.016047\n",
       "22       0.020448\n",
       "23       0.000000\n",
       "24       0.000000\n",
       "25       0.010155\n",
       "26       0.007089\n",
       "28       0.003045\n",
       "30       0.004830\n",
       "31       0.004274\n",
       "33       0.000000\n",
       "36       0.009776\n",
       "37       0.026216\n",
       "38       0.000000\n",
       "39       0.003515\n",
       "43       0.000000\n",
       "44       0.000000\n",
       "45       0.003883\n",
       "46       0.003524\n",
       "47       0.000000\n",
       "48       0.026571\n",
       "50       0.018493\n",
       "           ...   \n",
       "84780    0.000000\n",
       "85288    0.000000\n",
       "85541    0.000000\n",
       "85567    0.000000\n",
       "85607    0.000000\n",
       "85611    0.000000\n",
       "85699    0.000000\n",
       "85703    0.000000\n",
       "85814    0.000000\n",
       "85820    0.000000\n",
       "85959    0.000000\n",
       "86062    0.000000\n",
       "86111    0.000000\n",
       "86147    0.000000\n",
       "86331    0.000000\n",
       "86425    0.000000\n",
       "86494    0.000000\n",
       "86564    0.000000\n",
       "86653    0.000000\n",
       "86723    0.000000\n",
       "86838    0.000000\n",
       "86853    0.000000\n",
       "86996    0.000000\n",
       "87091    0.000000\n",
       "87172    0.000000\n",
       "87322    0.000000\n",
       "87493    0.000000\n",
       "87988    0.000000\n",
       "87992    0.000000\n",
       "88006    0.000000\n",
       "Name: dxy, dtype: float64"
      ]
     },
     "execution_count": 308,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "within_pi = pd.read_csv('/home/ipseg/Desktop/waples/temp/within_pi.txt', sep = 't', header = None)\n",
    "dxy = pd.read_csv('/home/ipseg/Desktop/waples/temp/dxy.txt', sep = '\\t', header = None)\n",
    "dxy.columns = ['locus', 'pop1', 'pop2', 'dxy']\n",
    "dxy.groupby('locus').mean()['dxy']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f0062542b50>"
      ]
     },
     "execution_count": 310,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEDCAYAAADZUdTgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X901PWd7/HnfGcmkx+TH0zIBI0UXWtBATHSywp2ETjY\nHLyuljYiK8Rel97tsYWr57BqKWp7enp7bVoXt3qpVsutWu6FU1AP3XaB2gWrawRqrCAWLdEiBjKZ\nISGEJPP7e/8YZvgRkvlBkskwr8c5aOY7me+8ZwKv+eT9/Xw/X4tpmiYiIpIXjGwXICIiI0ehLyKS\nRxT6IiJ5RKEvIpJHFPoiInlEoS8ikkdsg925a9cu7rvvPq666ioAJk6cyNe+9jUeeOABotEoVVVV\nNDY2UlBQwJYtW3jhhRcwDINFixZRX18/Ii9ARERSN2joA/zt3/4t//qv/5q4vWrVKpYuXUpdXR1r\n1qxh8+bN3H777axdu5ZNmzZht9upr6/n5ptvpry8fFiLFxGR9CRt75x77tbu3buZN28eAHPnzqWp\nqYm9e/cydepUnE4nDoeD2tpampubh6diERHJ2KAjfYvFQktLC/feey9dXV1885vfpK+vD7vdDoDL\n5aK9vR2fz4fL5Uo8rrKyEq/XO7yVi4hI2gYN/QkTJrB8+XIWLFjA4cOHaWhoIBKJJO4faAWHZCs7\n+P1+3nvvPaqqqrBarRmULSKSfyKRCF6vlylTplBYWJjRPgYN/erqahYsWADA+PHjGTt2LPv37ycY\nDFJQUIDH48HtduN2u/H5fInHeTweamtrB9zve++9x5IlSzIqWEQk361fv57Pf/7zGT120ND/9a9/\njdfr5R//8R/xer10dHTw5S9/ma1bt3Lbbbexfft2Zs+ezbRp03j44Yfp7u7GMAyam5tZvXr1gPut\nqqpKFD5u3LiMChcRyTdtbW0sWbIkkaGZGDT0582bx8qVK/n9739PKBTiu9/9LldffTUPPfQQGzdu\npKamhoULF2K1Wlm5ciXLli3DYrGwYsUKnE7ngPuNt3TGjRvHZZddlnHxIiL56ELa4oOGfklJCU8/\n/XS/7evWreu3ra6ujrq6uowLERGR4aczckVE8ohCX0Qkjyj0RUTyiEJfRCSPKPRFRPKIQl9EJI8o\n9EVE8ohCX0Qkjyj0RUTySNKLqMj5vfqHPVisBYnbkUiYGdd+looKXThGREYvhX6GTvSaFJ9xZTB/\noJdgMJjFikREklN7R0Qkjyj0RUTyiEJfRCSPKPRFRPKIQl9EJI8o9EVE8ohCX0Qkjyj0RUTyiEJf\nRCSPKPRFRPKIQl9EJI8o9EVE8ohCX0Qkjyj0RUTyiEJfRCSPKPRFRPKIQl9EJI8o9EVE8ohCX0Qk\njyj0RUTyiEJfRCSPKPRFRPKIQl9EJI8o9EVE8ohCX0Qkj6QU+n6/n/nz5/Pyyy9z9OhRGhoaWLJk\nCffffz/BYBCALVu2UF9fz6JFi9i0adOwFi0iIplJKfR/+tOfMmbMGAB+8pOfsHTpUtavX8+ECRPY\nvHkzvb29rF27ll/84he8+OKLPP/883R1dQ1r4SIikr6kod/S0kJLSws33XQTALt372bevHkAzJ07\nl6amJvbu3cvUqVNxOp04HA5qa2tpbm4e3spFRCRtSUO/sbGRVatWJW739fVht9sBcLlctLe34/P5\ncLlcie+prKzE6/UOQ7kiInIhBg39V155heuuu46amhoATNM86/5zbyfbLiIi2WUb7M7XXnuNw4cP\ns3PnTtra2igoKKCkpIRAIIDD4cDj8eB2u3G73fh8vsTjPB4PtbW1w168iIikZ9DQX7NmTeLrp556\nipqaGt555x22bdvGbbfdxvbt25k9ezbTpk3j4Ycfpru7G8MwaG5uZvXq1cNevIiIpGfQ0D+XxWJh\nxYoVPPTQQ2zcuJGamhoWLlyI1Wpl5cqVLFu2LPE9TqdzuGoWEZEMpRz6y5cvT3y9bt26fvfX1dVR\nV1c3NFWJiMiw0Bm5IiJ5RKEvIpJHFPoiInlEoS8ikkcU+iIieUShLyKSRxT6IiJ5RKEvIpJHFPoi\nInlEoS8ikkcU+iIieUShLyKSRxT6IiJ5RKEvIpJHFPoiInlEoS8ikkcU+iIieUShLyKSRxT6Gdrz\n4Qne/Ys322WIiKRFoZ+BcCTKe4d6eP/jjmyXIiKSFoV+BnzH+4BY+IuI5BKFfga8p0I/otAXkRyj\n0M+At7MXgHDUzHIlIiLpUehnwNupkb6I5CaFfgbaO+M9fRPT1GhfRHKHQj8D7afaOwBRtXhEJIco\n9DMQb+9AbLQvIpIrFPppMk0zMXsHNG1TRHKLQj9NJ3qCBEORxG2FvojkEoV+ms7s5wNE1NMXkRyi\n0E9TvJ9vscRua6QvIrlEoZ+m+HTNsiIroNAXkdyi0E+T93isvVNWYgMgotk7IpJDFPppird3yotj\noa+RvojkEoV+mrydvRTYDEoK4+0djfRFJHfYBruzr6+Pb33rW3R0dBAIBPjGN77BxIkTefDBB4lG\no1RVVdHY2EhBQQFbtmzhhRdewDAMFi1aRH19/Ui9hhHV3tlH1ZgibNbYkdxIVCN9Eckdg4b+jh07\nuPbaa1m2bBlHjhzhnnvu4frrr2fp0qXU1dWxZs0aNm/ezO23387atWvZtGkTdrud+vp6br75ZsrL\ny0fqdYyYk30hLh1bgtWIjfA10heRXDJoe+eWW25h2bJlABw5coRx48axe/du5s2bB8DcuXNpampi\n7969TJ06FafTicPhoLa2lubm5uGvfoRFoibRqEmB3Yr11DunlTZFJJcMOtKPW7x4Me3t7fz0pz/l\nnnvuwW63A+ByuWhvb8fn8+FyuRLfX1lZidd78V0/NhSOnYlrsxmJ9o4O5IpILkkp9Dds2MCBAwf4\n53/+57O2D7Ss8MW63HAoHAt4u9XAasRD/+J8rSJycRq0vbN//37a2toAmDRpEpFIhJKSEgKBAAAe\njwe3243b7cbn8yUe5/F4qK6uHsaysyMe+rH2zqkDuRrpi0gOGTT09+zZw7p16wDw+Xz09fUxc+ZM\ntm3bBsD27duZPXs206ZNY9++fXR3d9PT00NzczPTp08f/upHWGKkr/aOiOSoQds7//AP/8C3v/1t\nlixZgt/v5zvf+Q6TJ0/moYceYuPGjdTU1LBw4UKsVisrV65k2bJlWCwWVqxYgdPpHKnXMGLiPX27\n7XR7RwuuiUguGTT0HQ4Hjz/+eL/t8dH/merq6qirqxu6ykah8/f0NdIXkdyR0oFciUmEvt2K7dQA\nXwdyRSSXaBmGNJzZ09dIX0RykUb6aTirp2/GZ+9opC8iuUMj/TQEz+jpG4YFw6KRvojkFoV+Gk73\n9GNvm9VqKPRFJKco9NNwuqcfW1bZZjU0ZVNEcopCPw3heE//1GprNqtFI30RySkK/TQEE8swnG7v\n6ECuiOQShX4azpyyCRrpi0juUeinoV9P3zC04JqI5BSFfhrOXIYBYu2dqAlRHcwVkRyh0E9D4uQs\n++n2DmiuvojkDoV+Gvr39GP/17RNEckVCv00nNvTt54KfY30RSRXKPTTcL7ZO6DQF5HcodBPQ/CM\nBdfg9Ehfc/VFJFco9NPQb6Sv5ZVFJMco9NMQHrCnr5G+iOQGhX4a1NMXkVyn0E9DsN+Ca/Epmwp9\nEckNCv00hMJRbFYLxqlevk3tHRHJMQr9NITC0URrB8BqjV8yUSN9EckNCv00xEb61sRtjfRFJNco\n9NMQCkcSa+lDbJVN0IFcEckdCv00qL0jIrlOoZ+Gc0Nf7R0RyTUK/TSEwlHsZ/X0NdIXkdyi0E9D\nKBxJrKUPZ5yRq6WVRSRHKPRTFI2ahCPmOe0dnZErIrlFoZ+ieLDHz8aFM1fZVOiLSG5Q6Kcovu5O\ngV3z9EUkdyn0UxRfd8d2ZntHSyuLSI5R6Kfo3BU2gcQaPLqIiojkCoV+ihJr6Z/R07dYLNisBmGt\nsikiOUKhn6Lz9fQhNoNH7R0RyRUK/RSde33cOJvVUHtHRHKGLZVvamxspLm5mXA4zNe//nWmTJnC\ngw8+SDQapaqqisbGRgoKCtiyZQsvvPAChmGwaNEi6uvrh7v+EXO+nj7EQj8QimSjJBGRtCUN/bfe\neouWlhY2bNjA8ePH+dKXvsTMmTNZunQpdXV1rFmzhs2bN3P77bezdu1aNm3ahN1up76+nptvvpny\n8vKReB3DLnSenj6AzWahx6/2jojkhqTtnRkzZvDEE08AUFpaSl9fH3v27GHevHkAzJ07l6amJvbu\n3cvUqVNxOp04HA5qa2tpbm4e3upHUCL0z+3pG4ZOzhKRnJE09A3DoKioCIBNmzZx00030dvbi91u\nB8DlctHe3o7P58PlciUeV1lZidfrHaayR96A7R2bQdSMLdMgIjLapXwg99VXX+Wll17i0UcfPWu7\naZ4/7AbanqtCgxzIBZ2gJSK5IaXQf/3113nmmWd49tlncTqdFBcXEwwGAfB4PLjdbtxuNz6fL/EY\nj8dDdXX18FSdBQP29BOLrl1cH3IicnFKGvrd3d00NjbyzDPPUFZWBsCsWbPYunUrANu3b2f27NlM\nmzaNffv20d3dTU9PD83NzUyfPn14qx9BA/b0tbyyiOSQpLN3fvvb33L8+HHuu+8+IHYW6mOPPcbD\nDz/Mxo0bqampYeHChVitVlauXMmyZcuwWCysWLECp9M57C9gpAzW0we1d0QkNyQN/TvvvJM777yz\n3/Z169b121ZXV0ddXd3QVDbKDNjTN+LLK2ukLyKjn87ITdHA8/Q10heR3KHQT9Fga++ADuSKSG5Q\n6KdosGUYQO0dEckNCv0UDbbgGqi9IyK5QaGfoqSzdzRlU0RygEI/RadDf6Cevkb6IjL6KfRTFFZP\nX0QuAgr9FA3U07cnevoKfREZ/RT6KRqop2/VgVwRySEK/RQN3NPXSF9EcodCP0WhcBTDsGA1LGdt\nt9tityNRjfRFZPRT6KcoFI70a+3Ame0djfRFZPRT6KcoFI5ScJ7Q14FcEcklCv0UhcLRQUf6uk6u\niOQChX6KQpEotnMO4gJYDQuGRSN9EckNCv0UhULRfssqx9msBmEdyBWRHKDQT1EoHKHAPkDo2wyd\nkSsiOUGhn6KBevpwaqSv0BeRHKDQT4FpmoQi0X4nZsXFQl/tHREZ/RT6KQiGo5gmOAoGCn2LRvoi\nkhMU+inwB8IAFA4Y+gaRqIlpKvhFZHRT6KcgEIytsOmwDxz6cHp9HhGR0Uqhn4JAKBb6hQW2894f\nv3pWUKEvIqOcQj8F/mCsvTNYTx8gGFLoi8joptBPQaK9M0hPHzTSF5HRT6GfAn+KPX2N9EVktFPo\npyBpTz8R+pERq0lEJBMK/RQEgsmnbILaOyIy+in0U5C8p3/qQK5CX0RGOYV+CuI9/aRTNtXTF5FR\nTqGfgnhPP+mBXI30RWSUU+inIL4Mg8ORbJ6+DuSKyOim0E9BqsswqL0jIqOdQj8FqU7Z1No7IjLa\nKfRT4NcZuSJykVDopyCQmL0zQOifmr0TUHtHREa5lEL/ww8/ZP78+axfvx6Ao0eP0tDQwJIlS7j/\n/vsJBoMAbNmyhfr6ehYtWsSmTZuGr+oRllhwbYCevl3tHRHJEUlDv7e3l+9///vMmjUrse0nP/kJ\nS5cuZf369UyYMIHNmzfT29vL2rVr+cUvfsGLL77I888/T1dX17AWP1ICoQg2q4HVev63y6rZOyKS\nI5KGvsPh4Nlnn8Xtdie27d69m3nz5gEwd+5cmpqa2Lt3L1OnTsXpdOJwOKitraW5uXn4Kh9BgWBk\nwH4+nB7pq6cvIqPd+aejnMFqtWK1nh14fX192O12AFwuF+3t7fh8PlwuV+J7Kisr8Xq9Q1xudgSC\nkQH7+aAzckUkd1zwgdyBrgt7MV0v1h8MD9jPB7AaWntHRHJDRqFfXFycOHjr8Xhwu9243W58Pl/i\nezweD9XV1UNTZZYFQpEB5+gDWCwWrIZlSA7kth3r4Q/vfMrJ3uAF70tE5Fwph/6ZI/dZs2axdetW\nALZv387s2bOZNm0a+/bto7u7m56eHpqbm5k+ffrQVzzCTNPEn6SnD7GlGOJTOy/EExve4Ue/fJuG\n727l8fVvE41ePL8xiUj2Je3p79+/n8cee4zW1lbsdjvbtm3jRz/6EatWrWLjxo3U1NSwcOFCrFYr\nK1euZNmyZVgsFlasWIHT6RyJ1zCswpEo0aiZNPQLC6yc6A1d0HN1nQzw54+P4XYVY5omO5s/5a66\nSVwytuSC9isiEpc09CdPnsyLL77Yb/u6dev6baurq6Ourm5oKhslkp2YFVdWYudwey/dvUFKiwsy\neq63D7QTNWHBzMsJR6Ks33qAVu9Jhb6IDBmdkZvE6evjDv75WFYSC/oj3pMZP9fu99sA+C/XVFNT\nFfstqfUC9icici6FfhKJxdYGWFY5rqw4NoX1iK8no+cJhaM0H2in2lXMZ6pLT4d+u0JfRIaOQj+J\nxFr6g0zZBCh3xkf6mYX++x8doy8QZsbkcVgsFi491dLRSF9EhlLSnn6+S1w1K1lPPzHSTz+ko9Eo\nr//pEADXfMZJZ2cnAJVlDg57TiSWsygtLcUw9DktIplT6CeRbFnlOGeRDZvVklF7p/N4F6+9c5QC\nm8ER7wn2HziIYbVRYLdw9FiAf3/zY0LBPm6bcw3l5eUZvQ4REVB7J6lAkouix1ksFtwVhRz1nkz7\nbOR3/nIMfzDKxMtdlJWVU1RcTFFRCZUVsRZPyCyguFgzeETkwin0kwgkWVb5TG5XET3+MCd60jub\ndmfzUQAmX+E6a3uF0wHA8ZOBtPYnIjIQhX4Spy+VmDz0q8cUAukdzG3v6GXvwQ6qyh1Ulheddd+Y\n0lOh363QF5GhodBPItWePpwR+mkczN2+6xAm8Lnxpf3uqzgV+p0KfREZIgr9JBJXzUrS0weodsVG\n6qkezN3/0TE27ziIs8jG5eP69+ydxQVYDYvaOyIyZBT6SQQSZ+QOPtKPRqMU22IfEIeOdNLV1UVX\nVxfR6PlX3mz1nuR//p9dmKbJ8q9Mxm7r/6MwLBbKnQ6OdwcuqqWqRSR7NGUziVTX3unr7WHPe13Y\nrBb2HPBx74/ewMSkrMSB21XC9EnVXHFJGf5gmKa9n/KHd9sIhaP897+fyPhKg6Pe84d6RamDjhN+\n+gK6FKOIXDiFfhKnl2FI/lYVlziZfEUJLa1dRE2TSNTkRE+QI75e/vTh2VcRcxbZmDHJRSgU5N/f\nOEyJsxxn/7Z+YgZPV8+FreApIgIK/aT8gdTaO3FfuK6GL1xXA0C7p5VAIEhhcTmtvj56/GH8fT1U\nVji5duJ4jFNX3Orp6R5wf/EZPBe6bLOICCj0kwqEYn36VKZsnk9RUQljx7qoGhu73e5pxTDsicBP\nJj6D54RG+iIyBHQgN4l0pmwOh3joq70jIkNBoZ9EIBjBMCzYrNl5qwoLbLGrcin0RWQIKPSTCAQj\nOOxWLJbU2jHDoaLUQXdfmHDkwi+8LiL5TaGfRCAUzrifP1QqSh2YJrR3+rNah4jkPoV+Ev5gJGv9\n/LgxztjyDkeP9Wa1DhHJfQr9JPzBSNJllYdb/GCuQl9ELpRCP4l4Tz+b4qHfdqwvq3WISO5T6A8i\nEIoQjkQpKszuSL+8pAALGumLyIVT6A/Ccyy2Wma1qzirdVitBs5im0JfRC6YQn8QbadCdlxl9i9V\nWFZs50RPKO2rcomInEmhP4ijp0b64yqzO9KH0xdo2b3/aJYrEZFcptAfRFsi9LM/0r/iklgNO5s/\nzXIlIpLLFPqDGE3tndJiO1ddVsbegz6OdWkWj4hkRqE/iLZjPZQW23EW2bNdCgCzplZjmvD6n1ov\neF/RaDRxda9kV/kSkYuHQn8A0aiJp6OX6lEwyo+7YXIVVsPCzuZPL/jyiX/91Mvj6//Ii/9+gK1N\nH7Nl5/t0dw+8rr+IXBy0nv4AOk74CYWjXDKKQr+0uIDrJ7nZ876H7z77Fv/t1mu4/JKytBaDC4Yi\n/NsbH/N/t/2ZQCgKLSex2wzm1bqHsXIRGS0U+gMYTTN34HQ7ZsnNV9DnD9L8QTvNH7QzpszBpAku\nJk1wcdX4CiorCnEWFRAKRwiEIgSCEfyBCMdP+vnL4eNs33WIrpNBSovtTP2bCkIRg30tPl7f62Xh\nnCspL8/2KxWR4aTQH0Cbb/TM3AHo6+thW1MnLlcl068qp7qigAOHOjnpj9C07yhN+1KbyllSaOPW\nWeOZPbWC/X/txllaTqHDxq79baz7tw95eFllVpeRFpHhpdAfQFtHbObOaGrvFBWVUOIsA2BSaTmu\nkgh+fwBHcQXe4346u4N0dp0kFLVQVOjAZliwWi2EAn0UOqxUucq5rKoIm9XgD39sSVyM/fpJbj4+\n0snuP3vZ2vRXFsy64oJrjUaj5z1GUFpaimHoUJJItij0BxAf6VePkvbOQIqLnYytclFdFbsdvwbv\n2KrTPfrzbTvzYuyGxcLsa6vYtsfDMy/vY3x1KVOuHHtBdXV3d/Pyf+znaKdJy5FuSovtfKbKzj1/\nfy0VFRUXtG8RydyQh/4PfvAD9u7dC8Dq1auZOnXqUD/FiGjr6MFmNagsL8p2KSPCWWTnf9RP5rFf\nvsv/en4PDyydzrSrqjJq9ZimSdN7Hn67p4O+QOwaw57OAAdboatvPw/ePYPiwqGbBvv+x8fY8fan\nNB/w4Cwu4Ku3XMP1k4b+wHQ0GuV4VxfWc35T0W8vkkuGNPR3797NJ598woYNG2hpaWH16tVs2LBh\nKJ9ixBz19VLtKsJq5Ed/OxqNcukYC19dcBXrfvMhjzzTxGeqS/hsTRmV5YWUFNkodtgoctiocpVR\nUVZIudNBSaEt8cHQFwjz9gEP2946xJ8+9GI1LNR+roprPzuWzu4Au947wtsf+Hjgydd5qOHzfGZc\n2QXV6z3Wyf979SN+/8cjABQX2vAe7+M7zzYxfWIlX13wOVxljgsO5UjUpGnfEX716gd8dKSbkkIb\nleUFTLuygiJbmNvmXEP5EBwBj7fEgqEIHScCuMcUYRgWfajIkBrS0H/rrbeYP38+AFdeeSVdXV30\n9PRQUjJ6+uKpME2TqGlyZU3+tCHOPFB868xL2ffRcQ55evjE0zPo46yGBWexnVAoQu+pUT3ApM+U\ncvWEMi6tjrWJnMUFVBSNo+14iO27W1n+4x383XU1fGHapYyrLKGkyI7damC1Gtisltj/DQuGYTnr\nt41gKMKn7Sd5dddHbN99mEDIpMJpZ8bVldgiXfQEHRxojfD2B8fY27KLq8cXs3TB1Uz6m0vS/q3l\nRE+QN95t5ZWdLRw91oMFqCwroC8Y5RNPL594evlsjZPrr+m5oNA3TZNPPN1s/c+D7Gw+wkl/7H0s\nLDCoqXTw1f96NbVXj894/yJnGtLQ9/l8TJ48OXHb5XLh9XpzLvQtFgv/+4G5lAxhCyIXxA8Ulzhh\nQk0Vra2H6Qsa2ItKCZ6a/tnR2UlfIIzFWoA/EKUvGKHXH8RmWBjnKsQ9ppAJ1cVE/Z0UGOee4Wty\n+0w3k68Yw0uv/ZU/vNPKH94Z/OxiC2C1WrAasQ+CPn+Y+GlpBTaDz1/t5vOT3FitBu2eMKWldiZP\nrOL9jzv4z71H2PvxSR5cu4fCAivOIhslRXZKCm047AZRM3YSXuTUn2j8j2nS1ROi40QAALvVwtzr\nL+HvplRwyNOHs7Scw55u3nj3CAdbT/Ktp/dwZU0pl19SiruiEJvVoKioEKs1NjoPhaMEgmFO9vQR\nCkcJhqOEQhGC4Si9/ggfHe2m62Rs9VS7zcKlY2Mfgp+2n6TlaB+PPtfMVeM/4rrPVXH5JWWUFNmx\nWY1TfyyJrwGiZuw1mOapr02TYCiCPxihzx+i60QPgVO3A6EowVCEMmcxleWFuMoKGVNWSJHDhvXU\nh61hWDAsFgyDAT80L/REwUF//ud5zn5bLOfePM9j0vi8j5omJ3qCBIIRrqwpT/wcLxbDeiDXNM3z\n/tAikdhIpq2tbTif/oL1DXKCavvRQ9iPeRO3T3R1EsROKHB6zfuODh+GYcXfeyKtbZk+bries9AY\nQxFQZIOo9Rgup5WKijHE/rXZ6Og4jmFYqagoAkIQ6uJ4Vwcnu7v67e/D94OUlpVzlcvE4e+kN2Qh\nQgHhiIlpmgSCIQAMqz0RXOFQOPa1YcVptVBSZKXQCHCZ20ml8yRHPz3Z7zWUW2Hu5AIOHvJy7ESY\nQMCgs9fEE44y0GoTllP/sVjAagGnA8ZWFPGZKgc2/1G2/n4/RcVOKirGYAFuvMrgw0PdfOoLceAv\nHRz4y/n3m0yh3aDSaaXMEaKmykmlqxgIMcldwGFPHyf6TP78wUe8/8FHmT2BZOxrt0/hhimXZLuM\nhHhmxjM0E0Ma+m63G5/Pl7jd3t5OVVVVv+/zemNhuWTJkqF8ehGRIbX6P7Jdwfl5vV4mTJiQ0WOH\nNPRvvPFGnnzySe688072799PdXU1xcX9pzxOmTKF9evXU1VVhdWa3evPiojkikgkgtfrZcqUKRnv\nw2IOcUPu8ccfZ8+ePVitVh599FEmTpw4lLsXEZELMOShLyIio9fFdVhaREQGpdAXEckjCn0RkTwy\nrKEfCoVYuXIld911Fw0NDRw+fLjf92zZsoX6+noWLVrEpk2bEtt3797NrFmz2Llz53CWOKAf/OAH\nLF68mMWLF7Nv376z7nvzzTe54447WLx4MWvXrk3pMSMtk/o//PBD5s+fz/r160e63LNkUntjYyOL\nFy+mvr6e3/3udyNd8lnSrb+vr4/77ruPhoYGFi1alLW/83GZvP8Afr+f+fPn8/LLL49kuf2kW/+u\nXbu44YYbaGhooKGhge9///vZKDshk/d/y5Yt3H777Xz5y1/mtddeG/wJzGH00ksvmd/73vdM0zTN\nN954w7z//vvPur+np8esq6szu7u7Tb/fb956663m8ePHzUOHDpn33nuvuXz5cnPHjh3DWeJ57dq1\ny/z6179umqZpHjx40LzzzjvPuv+WW24x29razGg0at51113mwYMHkz5mJGVSf29vr9nQ0GA+8sgj\n5i9/+cv/LBwyAAAFG0lEQVRslG2aZma1NzU1mf/0T/9kmqZpdnZ2mnPmzBnxuuMyqf83v/mN+dxz\nz5mmaZqtra3mF7/4xRGvOy6T+uP+5V/+xfzKV75ivvzyyyNa85ky/be7YsWKbJTbTyb1d3R0mF/8\n4hfNnp4es7293XzkkUcGfY5hHemfuRbPzJkzaW5uPuv+d999l6lTp+J0OnE4HNTW1tLc3Izb7eap\np57C6XQOZ3kDGmgNIYDDhw9TXl5OdXU1FouFm266iaampkEfkwv1FxQU8Oyzz+J2Z/eyiZnUPmPG\nDJ544gkgtuJlb2/vsC4NMNT133LLLSxbtgyAI0eOcMkl2TsDNJP6AVpaWmhpaWHOnDlZe+8vpP7R\nIpP6m5qamDVrFsXFxVRVVfG9731v0OcY1tD3+Xy4XK7YExkGFouFcDicuP/YsWOJ+wEqKyvxer0U\nFhZmdVVBn8/HmDFjErddLlfiTGOv13tWzfH1hc73mPiZxyMtk/qtVisOh2PEaz1XJrUbhkFRUWwJ\n7E2bNjFnzpysXf0rk/rjFi9ezAMPPMCqVatGruBzZFp/Y2NjVuuOy7T+lpYW7r33Xu666y7efPPN\nkS36DJnU39rait/v595772XJkiVJP8iG7IzcX/3qV2f15CE2kj9TshFANkcIgxmsroHuMwdYdygb\nMql/tEin9ldffZXNmzezbt264S4rZenUv2HDBg4cOMADDzzAli1bhru0lKRS/yuvvMJ1111HTU3N\nqPv7lEr9EyZMYPny5SxYsIDDhw9z991387vf/Q6bLfvXmEqlftM06ezsZO3atbS2tnL33XezY8eO\nAR83ZK/qjjvu4I477jhr26pVq/B6vUycOJFQKIRpmme9keeu1ePxeKitrT1rH9kIzsHWEKquru5X\ns9vtxm63p7Tu0EjIpP7RItPaX3/9dX72s5/x3HPPZa0tCJnVv3//fiorKxk3bhyTJk0iEonQ0dFx\n1qhuNNf/2muvcfjwYXbu3ElbWxsFBQWMGzeOmTNn5kT91dXVLFiwAIDx48czduxYPB4PNTU1I1s8\nmdVfVFTE9ddfj2EYjB8/npKSkkH//gxrD+XGG29k69atAOzYsYMbbrjhrPuvvfZa9u3bR3d3Nz09\nPTQ3NzN9+vTE/aZpZmXkcOONN7Jt2zaAfmsI1dTUcPLkSVpbWwmHw+zcuZMvfOELgz4mF+qPy/ZI\nLZPau7u7aWxs5Omnn6asLPMLswyFTOrfs2dP4rcTn89Hb29vVgI/0/rXrFnDpk2b2LhxI3fccQff\n/OY3sxL4mdb/61//OvH+e71ejh07RnV1dc7Uf+ONN/LWW28lRvzJ/v4M6zIM0WiU1atXc+jQIRwO\nB4899hjV1dX87Gc/Y8aMGVx33XVs27aNn//851gsFhoaGrj11lt57bXXeO655/j4448ZM2YMbreb\nn//858NV5nmdu4bQ+++/T2lpKfPnz+ePf/wjP/7xjwGoq6vjnnvuOe9jsrnuULr1v/fee/zwhz+k\ntbUVu91OdXU1Tz755JBcEWq4a9+4cSNPPfUUl19+eWIfjY2NWTsgmm79gUCAb3/727S1teH3+1mx\nYgVz5szJSu2Z1H+mp556issuu4wvfelL2SgdSL/+np4eVq5cSXd3N6FQiOXLlzN79uycqR9g48aN\nifb6N77xDebOnTvg/rX2johIHtEZuSIieUShLyKSRxT6IiJ5RKEvIpJHFPoiInlEoS8ikkcU+iIi\neUShLyKSR/4/G0vkI8WZA9sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7effe8c079d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(dxy.groupby('locus').mean()['dxy'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>pop</th>\n",
       "      <th colspan=\"10\" halign=\"left\">1</th>\n",
       "      <th>...</th>\n",
       "      <th colspan=\"10\" halign=\"left\">10</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>ind</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMHAMM10_0002</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMHAMM10_0005</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMHAMM10_0008</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMHAMM10_0011</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMHAMM10_0012</th>\n",
       "      <th>...</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMUW10_0004</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMUW10_0007</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMUW10_0008</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMUW10_0009</th>\n",
       "      <th colspan=\"2\" halign=\"left\">CMUW10_0011</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>...</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pop</th>\n",
       "      <th>ind</th>\n",
       "      <th>chr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">1</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">CMHAMM10_0002</th>\n",
       "      <th>a</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b</th>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>...</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">CMHAMM10_0005</th>\n",
       "      <th>a</th>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>...</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b</th>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>...</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMHAMM10_0008</th>\n",
       "      <th>a</th>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>...</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.010638</td>\n",
       "      <td>0.021277</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010638</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 348 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "pop                              1                                     \\\n",
       "ind                   CMHAMM10_0002           CMHAMM10_0005             \n",
       "chr                               a         b             a         b   \n",
       "pop ind           chr                                                   \n",
       "1   CMHAMM10_0002 a        0.000000  0.010638      0.010638  0.010638   \n",
       "                  b        0.010638  0.000000      0.021277  0.021277   \n",
       "    CMHAMM10_0005 a        0.010638  0.021277      0.000000  0.000000   \n",
       "                  b        0.010638  0.021277      0.000000  0.000000   \n",
       "    CMHAMM10_0008 a        0.010638  0.000000      0.021277  0.021277   \n",
       "\n",
       "pop                                                                    \\\n",
       "ind                   CMHAMM10_0008           CMHAMM10_0011             \n",
       "chr                               a         b             a         b   \n",
       "pop ind           chr                                                   \n",
       "1   CMHAMM10_0002 a        0.010638  0.010638      0.000000  0.021277   \n",
       "                  b        0.000000  0.021277      0.010638  0.010638   \n",
       "    CMHAMM10_0005 a        0.021277  0.000000      0.010638  0.010638   \n",
       "                  b        0.021277  0.000000      0.010638  0.010638   \n",
       "    CMHAMM10_0008 a        0.000000  0.021277      0.010638  0.010638   \n",
       "\n",
       "pop                                              ...             10            \\\n",
       "ind                   CMHAMM10_0012              ...    CMUW10_0004             \n",
       "chr                               a         b    ...              a         b   \n",
       "pop ind           chr                            ...                            \n",
       "1   CMHAMM10_0002 a        0.000000  0.021277    ...       0.000000  0.010638   \n",
       "                  b        0.010638  0.010638    ...       0.010638  0.000000   \n",
       "    CMHAMM10_0005 a        0.010638  0.010638    ...       0.010638  0.021277   \n",
       "                  b        0.010638  0.010638    ...       0.010638  0.021277   \n",
       "    CMHAMM10_0008 a        0.010638  0.010638    ...       0.010638  0.000000   \n",
       "\n",
       "pop                                                                            \\\n",
       "ind                   CMUW10_0007           CMUW10_0008           CMUW10_0009   \n",
       "chr                             a         b           a         b           a   \n",
       "pop ind           chr                                                           \n",
       "1   CMHAMM10_0002 a      0.000000  0.010638    0.010638  0.000000    0.010638   \n",
       "                  b      0.010638  0.021277    0.021277  0.010638    0.021277   \n",
       "    CMHAMM10_0005 a      0.010638  0.000000    0.021277  0.010638    0.021277   \n",
       "                  b      0.010638  0.000000    0.021277  0.010638    0.021277   \n",
       "    CMHAMM10_0008 a      0.010638  0.021277    0.021277  0.010638    0.021277   \n",
       "\n",
       "pop                                                    \n",
       "ind                             CMUW10_0011            \n",
       "chr                           b           a         b  \n",
       "pop ind           chr                                  \n",
       "1   CMHAMM10_0002 a    0.010638    0.010638  0.021277  \n",
       "                  b    0.000000    0.000000  0.010638  \n",
       "    CMHAMM10_0005 a    0.021277    0.021277  0.010638  \n",
       "                  b    0.021277    0.021277  0.010638  \n",
       "    CMHAMM10_0008 a    0.000000    0.000000  0.010638  \n",
       "\n",
       "[5 rows x 348 columns]"
      ]
     },
     "execution_count": 266,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_pdist = pdist_of_locus[193]\n",
    "test_pdist = pd.concat([chr_index, test_pdist], axis = 1)\n",
    "test_pdist = test_pdist.set_index(['pop', 'ind', 'chr'])\n",
    "test_pdist.columns = test_pdist.index\n",
    "test_pdist.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the above functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: 0.0096699399890889243,\n",
       " 2: 0.0066560473318921379,\n",
       " 3: 0.0084424440807419534,\n",
       " 4: 0.010606636271529886,\n",
       " 5: 0.0095828345625619434,\n",
       " 6: 0.0085659021829234596,\n",
       " 7: 0.0070359112912304397,\n",
       " 8: 0.0081560283687943259,\n",
       " 9: 0.0099181669394435342,\n",
       " 10: 0.0098404255319148922}"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_within_pi(test_pdist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0088474336550121487"
      ]
     },
     "execution_count": 232,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(get_within_pi(test_pdist).values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(1, 2): 0.0083541927409261571,\n",
       " (1, 3): 0.0091888297872340426,\n",
       " (1, 4): 0.010355718085106384,\n",
       " (1, 5): 0.0097028613352898013,\n",
       " (1, 6): 0.0093568665377175998,\n",
       " (1, 7): 0.0085486322188449843,\n",
       " (1, 8): 0.0091090425531914904,\n",
       " (1, 9): 0.009942716857610473,\n",
       " (1, 10): 0.0096742021276595737,\n",
       " (2, 3): 0.0074155193992490607,\n",
       " (2, 4): 0.0088196182728410497,\n",
       " (2, 5): 0.0081243796124465927,\n",
       " (2, 6): 0.0075378313801342583,\n",
       " (2, 7): 0.0066601108528517774,\n",
       " (2, 8): 0.0076658322903629525,\n",
       " (2, 9): 0.0083517858862039076,\n",
       " (2, 10): 0.0079787234042553185,\n",
       " (3, 4): 0.0095827792553191491,\n",
       " (3, 5): 0.0089600146735143056,\n",
       " (3, 6): 0.0082688588007736926,\n",
       " (3, 7): 0.0076367781155015193,\n",
       " (3, 8): 0.0087101063829787221,\n",
       " (3, 9): 0.0092062193126022911,\n",
       " (3, 10): 0.0087765957446808516,\n",
       " (4, 5): 0.0099676724137931043,\n",
       " (4, 6): 0.0097769584139264995,\n",
       " (4, 7): 0.0088692059270516696,\n",
       " (4, 8): 0.0094331781914893595,\n",
       " (4, 9): 0.010024549918166937,\n",
       " (4, 10): 0.0099110704787234022,\n",
       " (5, 6): 0.0091542719935970105,\n",
       " (5, 7): 0.0082145477413269035,\n",
       " (5, 8): 0.0088499633162142342,\n",
       " (5, 9): 0.0095236751509678863,\n",
       " (5, 10): 0.009354365370506235,\n",
       " (6, 7): 0.0078405636916275209,\n",
       " (6, 8): 0.0091876208897485497,\n",
       " (6, 9): 0.0094851956554084221,\n",
       " (6, 10): 0.0088853965183752423,\n",
       " (7, 8): 0.0075987841945288747,\n",
       " (7, 9): 0.0083586626139817623,\n",
       " (7, 10): 0.0081212006079027354,\n",
       " (8, 9): 0.0087970540098199664,\n",
       " (8, 10): 0.0089760638297872331,\n",
       " (9, 10): 0.0095130932896890352}"
      ]
     },
     "execution_count": 242,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_dxy(test_pdist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Distance matrix 'X' must be symmetric.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-169-80a7124faa25>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;31m# fails with missing data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mspatial\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdistance\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msquareform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtest_pdist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mspatial\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdistance\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msquareform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtest_pdist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m/home/ipseg/anaconda/lib/python2.7/site-packages/scipy/spatial/distance.pyc\u001b[0m in \u001b[0;36msquareform\u001b[1;34m(X, force, checks)\u001b[0m\n\u001b[0;32m   1484\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'The matrix argument must be square.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1485\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mchecks\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1486\u001b[1;33m             \u001b[0mis_valid_dm\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthrow\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'X'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1487\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1488\u001b[0m         \u001b[1;31m# One-side of the dimensions is set here.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/ipseg/anaconda/lib/python2.7/site-packages/scipy/spatial/distance.pyc\u001b[0m in \u001b[0;36mis_valid_dm\u001b[1;34m(D, tol, throw, name, warning)\u001b[0m\n\u001b[0;32m   1567\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1568\u001b[0m                     raise ValueError(('Distance matrix \\'%s\\' must be '\n\u001b[1;32m-> 1569\u001b[1;33m                                      'symmetric.') % name)\n\u001b[0m\u001b[0;32m   1570\u001b[0m                 \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1571\u001b[0m                     \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Distance matrix must be symmetric.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Distance matrix 'X' must be symmetric."
     ]
    }
   ],
   "source": [
    "test_pdist.loc[3,3]\n",
    "# fails with missing data\n",
    "scipy.spatial.distance.squareform(test_pdist.loc[3,3]).mean()\n",
    "scipy.spatial.distance.squareform(test_pdist.loc[1,1]).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f008af47310>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEDCAYAAADZUdTgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9w1PW97/Hn/s6PTUISsgkNHPxZ7AGEiMcKdhC4aC5e\nR8RGpEI610OnXlsYnaHWUmzPmY63Y3OOpcd6qVbL9cfhXJiC9tI5LT+0QO0VARtrAKXWqBgD2eyS\nHyzJJvvre/9Idk2AZH+QZLPu6zHD4H6/u999bwZf+8n7+/l+vibDMAxERCQrmNNdgIiIjB2FvohI\nFlHoi4hkEYW+iEgWUeiLiGQRhb6ISBaxDrfz0KFDPPjgg1x99dUATJs2jW984xs8/PDDRCIRysrK\nqKurw263s3PnTl588UXMZjPLly+npqZmTD6AiIgkbtjQB/jyl7/Mv/3bv8Uer1+/nlWrVlFdXc3G\njRvZsWMHS5cuZdOmTWzfvh2bzUZNTQ233HILRUVFo1q8iIgkJ2575/xrtw4fPsyiRYsAWLhwIQcP\nHqShoYGZM2fidDpxOBxUVVVRX18/OhWLiEjKhh3pm0wmGhsbeeCBB+js7OTb3/42fr8fm80GQElJ\nCa2trXi9XkpKSmKvKy0txePxjG7lIiKStGFDf+rUqaxZs4YlS5bQ1NREbW0t4XA4tn+oFRzirezQ\n09PDsWPHKCsrw2KxpFC2iEj2CYfDeDweZsyYQU5OTkrHGDb0y8vLWbJkCQBTpkxh4sSJHD9+nEAg\ngN1ux+1243K5cLlceL3e2OvcbjdVVVVDHvfYsWOsXLkypYJFRLLdli1buP7661N67bCh/9vf/haP\nx8M//uM/4vF4aGtr46677mLXrl3ccccd7Nmzh/nz5zNr1iweffRRfD4fZrOZ+vp6NmzYMORxy8rK\nYoVXVFSkVLiISLZpaWlh5cqVsQxNxbChv2jRItatW8drr71GMBjkn//5n/nSl77EI488wrZt26is\nrGTZsmVYLBbWrVvH6tWrMZlMrF27FqfTOeRxoy2diooKJk+enHLxIiLZ6FLa4sOGfn5+Pk8//fQF\n2zdv3nzBturqaqqrq1MuRERERp+uyBURySIKfRGRLKLQFxHJIgp9EZEsotAXEckiCn0RkSyi0BcR\nySIKfRGRLKLQFxHJInFvoiIXt+sPb2KxfbbKXTAQ4Cv/cA2FhYVprEpEZHgK/RR1B8zk5U6IPe7p\nPUcwGEpjRSIi8am9IyKSRRT6IiJZRKEvIpJFFPoiIllEoS8ikkUU+iIiWUShLyKSRRT6IiJZRKEv\nIpJFFPoiIllEoS8ikkUU+iIiWUShLyKSRRT6IiJZRKEvIpJFFPoiIllEoS8ikkUU+iIiWUShLyKS\nRRT6IiJZRKEvIpJFFPoiIllEoS8ikkUU+iIiWUShLyKSRRIK/Z6eHhYvXswrr7zC6dOnqa2tZeXK\nlTz00EMEAgEAdu7cSU1NDcuXL2f79u2jWrSIiKQmodD/xS9+QXFxMQBPPvkkq1atYsuWLUydOpUd\nO3bQ3d3Npk2beP7553nppZd44YUX6OzsHNXCRUQkeXFDv7GxkcbGRm6++WYADh8+zKJFiwBYuHAh\nBw8epKGhgZkzZ+J0OnE4HFRVVVFfXz+6lYuISNLihn5dXR3r16+PPfb7/dhsNgBKSkpobW3F6/VS\nUlISe05paSkej2cUyhURkUsxbOj/5je/Yfbs2VRWVgJgGMag/ec/jrddRETSyzrczgMHDtDU1MT+\n/ftpaWnBbreTn59Pb28vDocDt9uNy+XC5XLh9Xpjr3O73VRVVY168SIikpxhQ3/jxo2x/37qqaeo\nrKzk7bffZvfu3dxxxx3s2bOH+fPnM2vWLB599FF8Ph9ms5n6+no2bNgw6sWLiEhyhg3985lMJtau\nXcsjjzzCtm3bqKysZNmyZVgsFtatW8fq1atjz3E6naNVs4iIpCjh0F+zZk3svzdv3nzB/urqaqqr\nq0emKhERGRW6IldEJIso9EVEsohCX0Qkiyj0RUSyiEJfRCSLKPRFRLKIQl9EJIso9EVEsohCX0Qk\niyj0RUSyiEJfRCSLKPRFRLKIQl9EJIso9EVEsohCX0Qkiyj0RUSyiEJfRCSLKPRT1OYL0nmuN91l\niIgkRaGfAsMw2PXnM/zhraZ0lyIikhSFfgrOdgXoDRr4e0PpLkVEJCkK/RR42v0AhMJGmisREUmO\nQj8Fno5uAMKRSJorERFJjkI/BZ+N9BX6IpJZFPop8HSovSMimUmhn4LoSD8SMTAMBb+IZA6Ffgqi\nPX3QaF9EMotCPwXRkT5AWH19EckgCv0kBUNh2n2fXYkbimikLyKZQ6GfJG9Hz6DHmsEjIplEoZ+k\ngf18UHtHRDKLQj9J0X6+3WoCdCJXRDKLQj9J0Tn6RflWQCN9EcksCv0kRUf6hXl9oa+evohkEoV+\nkjztfT39ojwLoPaOiGQWhX6SPB1+CvLsOGx9PzotuiYimcQ63E6/38/3vvc92tra6O3t5Vvf+hbT\npk3ju9/9LpFIhLKyMurq6rDb7ezcuZMXX3wRs9nM8uXLqampGavPMKbOdPqpKM3HYtGJXBHJPMOG\n/r59+7j22mtZvXo1p06d4r777uO6665j1apVVFdXs3HjRnbs2MHSpUvZtGkT27dvx2azUVNTwy23\n3EJRUdFYfY4xYRgGPYEwuQ4rFnPfCF89fRHJJMO2d2677TZWr14NwKlTp6ioqODw4cMsWrQIgIUL\nF3Lw4EEaGhqYOXMmTqcTh8NBVVUV9fX1o1/9GAtHDAwDbFYzVnPfSD+skb6IZJBhR/pRK1asoLW1\nlV/84hfcd9992Gw2AEpKSmhtbcXr9VJSUhJ7fmlpKR6PZ3QqTqNgqG9Ub7NasFg00heRzJNQ6G/d\nupUTJ07wne98Z9D2oZYV/rwuN/xZ6JuxmKM9fYW+iGSOYds7x48fp6WlBYBrrrmGcDhMfn4+vb19\nC4653W5cLhculwuv1xt7ndvtpry8fBTLTo9gKAxE2zt929TeEZFMMmzoHzlyhM2bNwPg9Xrx+/3M\nnTuX3bt3A7Bnzx7mz5/PrFmzOHr0KD6fj66uLurr65kzZ87oVz/GBo30o7N3NGVTRDLIsO2dr33t\na3z/+99n5cqV9PT08E//9E9Mnz6dRx55hG3btlFZWcmyZcuwWCysW7eO1atXYzKZWLt2LU6nc6w+\nw5gZ1NNXe0dEMtCwoe9wOHjiiScu2B4d/Q9UXV1NdXX1yFU2Dg0c6VvR7B0RyTwJnciVPtGevt1q\nxmJopC8imUfLMCQhOtK3Dpi9o5G+iGQSjfSTEBh4Irc/6zXSF5FMopF+EkLR0LdYMJlMWMwmrb0j\nIhlFoZ+EgSdyAawWs1bZFJGMotBPQuxEbv+yyhaLSe0dEckoCv0kXHSkr/aOiGQQhX4SAgN6+qCR\nvohkHoV+EgZO2YS+kb5O5IpIJlHoJ2HggmsAVrOZcDjyuV1VVEQ+fxT6SYhO2Rx4ItcAIgp9EckQ\nCv0kBMODe/pWS//N0dXiEZEModBPQiB4/uwdrb8jIplFoZ+E83v6lv6Rvk7mikimUOgn4WLz9AHC\nGumLSIZQ6Cch1tO39vf0dSMVEckwCv0kBM/r6au9IyKZRqGfhAvbO/1r6mvRNRHJEAr9JOhErohk\nOoV+EqI9/egJXE3ZFJFMo9BPQjAUwWoxY+4/gauLs0Qk0yj0kxAMRmKtHRjY3tFIX0Qyg0I/CcFw\neFDoq70jIplGoZ+EYOi8kb5Z7R0RySwK/SScH/rRnn5IUzZFJEMo9JPQF/qW2OPYPH21d0QkQyj0\nkxAMhYc4kav2johkBoV+Ei5s7+hErohkFoV+giIRg1DYuGhPXydyRSRTKPQTFIrdNWvg7B2N9EUk\nsyj0ExSM3R934Inc/pF+RCN9EckMCv0EBfoXW7NebMqmRvoikiEU+gk6f1llALPZhMmk0BeRzKHQ\nT1AodGFPH/pG+5qyKSKZQqGfoIuN9KEv9HVxlohkCoV+gi52Ihf65uqrvSMimcKayJPq6uqor68n\nFApx//33M2PGDL773e8SiUQoKyujrq4Ou93Ozp07efHFFzGbzSxfvpyamprRrn/MBM67a1aU1WKm\nJxBOR0kiIkmLG/pvvvkmjY2NbN26lY6ODu68807mzp3LqlWrqK6uZuPGjezYsYOlS5eyadMmtm/f\njs1mo6amhltuuYWioqKx+ByjLjhsTz+YjpJERJIWt71zww038LOf/QyAgoIC/H4/R44cYdGiRQAs\nXLiQgwcP0tDQwMyZM3E6nTgcDqqqqqivrx/d6sdQNPStFxnph0IRDEMnc0Vk/Isb+mazmdzcXAC2\nb9/OzTffTHd3NzabDYCSkhJaW1vxer2UlJTEXldaWorH4xmlssfeZydyz+vpW00YgFZXFpFMkPCJ\n3FdffZWXX36ZH/7wh4O2DzXC/byNfIP9PX277cKRPmhNfRHJDAmF/uuvv84zzzzDs88+i9PpJC8v\nj0AgAIDb7cblcuFyufB6vbHXuN1uysvLR6fqNBiupw9adE1EMkPc0Pf5fNTV1fHMM89QWFgIwLx5\n89i1axcAe/bsYf78+cyaNYujR4/i8/no6uqivr6eOXPmjG71Y2i4efqgNfVFJDPEnb3zu9/9jo6O\nDh588EEATCYTjz/+OI8++ijbtm2jsrKSZcuWYbFYWLduHatXr8ZkMrF27VqcTueof4CxMmRPP7qm\nvto7IpIB4ob+Pffcwz333HPB9s2bN1+wrbq6murq6pGpbJyJN9JXe0dEMoGuyE1QcJiLs0DtHRHJ\nDAr9BA050rdqeWURyRwK/QSpvSMinwcK/QQFw/FO5Cr0RWT8U+gnKP6UTbV3RGT8U+gnKBAc4kSu\nVe0dEckcCv0EhYYa6ZujyzAo9EVk/FPoJ2jInr61v6ev9o6IZACFfoKG6unbNHtHRDKIQj9BQ12c\nZdHFWSKSQRT6CQoEI5hNYDGbBm3X0soikkkU+gkKhiNYrRZMpvNDv++x2jsikgkU+gkKhSIXtHZg\n4DIMCn0RGf8U+gkKhsIXD31dnCUiGUShn6DgECN9s8mE2WxSe0dEMoJCP0GBUAT7RUIf+qZt6uIs\nEckECv0E9Y30LRfdZ7WY1NMXkYyg0E9QMBSJnbQ9n8ViJqyevohkAIV+AgzDIBgKD93esaq9IyKZ\nQaGfgEAogmFAjv3itxS2mM06kSsiGUGhn4DeQN8SDA77xXv6NquJcMQgYij4RWR8U+gnoCcQAoYO\n/ehc/WBQfX0RGd8U+gmIjfRtw4d+b/+NVkRExiuFfgKioT9UTz820g9ppC8i45tCPwHREfyQ7R1r\ndKSv0BeR8U2hn4BoTz9nqNDvX245oJG+iIxzCv0ExO3pW6MnctXTF5HxTaGfgJ7YlM3he/q9GumL\nyDin0E9A3J6+pmyKSIZQ6CegN15Pv//uWZqyKSLjnUI/AYnO09eJXBEZ7xT6CeiJN0+//0RuQO0d\nERnnFPoJSLSnr5G+iIx3Cv0E9PQm1tMPqKcvIuOcQj8BGumLyOeFQj8BvQnO01dPX0TGu4RC//33\n32fx4sVs2bIFgNOnT1NbW8vKlSt56KGHCAQCAOzcuZOamhqWL1/O9u3bR6/qMZbw7B21d0RknIsb\n+t3d3Tz22GPMmzcvtu3JJ59k1apVbNmyhalTp7Jjxw66u7vZtGkTzz//PC+99BIvvPACnZ2do1r8\nWOkNhrGYTdiGuF1ibPaO2jsiMs7FDX2Hw8Gzzz6Ly+WKbTt8+DCLFi0CYOHChRw8eJCGhgZmzpyJ\n0+nE4XBQVVVFfX396FU+hnoCoSH7+TDwRK5CX0TGt4s3qQewWCxYLIMDz+/3Y7PZACgpKaG1tRWv\n10tJSUnsOaWlpXg8nhEuNz16A+EhZ+6A2jsikjku+USuMcR9YYfanol6AmEctqG/Hy1aWllEMkRK\noZ+Xlxc7eet2u3G5XLhcLrxeb+w5breb8vLykakyzXqD4WHbOyaTCavFNCLtnVA4QndP8JKPIyJy\nMQmH/sCR+7x589i1axcAe/bsYf78+cyaNYujR4/i8/no6uqivr6eOXPmjHzFadAbp6cPYLOaOee/\n9LB+bPMhvvbo7/jOk39k76GTl3w8EZGB4vb0jx8/zuOPP05zczM2m43du3fzL//yL6xfv55t27ZR\nWVnJsmXLsFgsrFu3jtWrV2MymVi7di1Op3MsPsOoCoUjhMLGkNM1owrzbLR29BAMRYac5RNPh6+X\n+r+2kmO38remDj5o6mDBnCkpH09E5HxxQ3/69Om89NJLF2zfvHnzBduqq6uprq4emcrGiXg3RY8q\nzLPhbu+h5UwXU8oLUnqvt95rwTBgxS3TaHL7ePXIJ5d0PBGR82kIGUe8JRiiCvP7ZjOd9nal/F6H\n33UD8OUZFVS6+n5LavacS/l4IiLnU+jHEe+m6FGFeX2hf8qbWkgHgmHq/9pKZZmTyjInk/tD/9NW\nhb6IjJy47Z1sF28JhqjoSP+UJ7WRfsMHXnoDYW6YXkEkEqEot+/E+UfNbbErmwsKCjCb9T0tIqlT\ngsTx2WJrw4d+QYoj/UgkQmdnJ6+/3TdTZ/pUJ59++imH3vkIkwlOfNzOq4dPsnP/u/h8vhQ+gYjI\nZzTSjyPeCptRNquZCU47p5Ls6ft8Pl75w3H+X8MZHDYzn5zu4C/HWsl3FlGYb+dsd4h8Z2HK9YuI\nDKSRfhyJ9vQBykty8Xb4k16OwXMWeoMRvnR5KQWFReTm5QFQXJBDTyCMv/8mLiIil0qhH0eis3eg\nL/QNA06fSW60f6LpLADTrygdtH1CgQPom78vIjISFPpxfHZT9PihX1GcAyR3Mrep9Ryt7b1McTmZ\n4HQM2lfcH/rtvp6EjyciMhyFfhyfzd4ZvqcfiUQo6Mt8PvzUS2dnJ52dnUQiw6/H8+pbpwCYceXE\nC/ZFvwQ00heRkaITuXFEe/rx2jt+fxfuM32za+pPeMixGnR3d3HHgr+nqKjooq95/S/N/OGtU+Tn\nWLls0oUna2PtnXO9QOYvaSEi6aeRfhzJ9PRLJ/Qtl+DpDNDWZcLuyB3yuW+95+an//FnHHYLi65z\nYe5fnnmgXIcVh81Cu0b6IjJCNNKPozeJnr7FYqKsOBdPu5+dr3+IxWziLx+cZc7fT2L21WXk5Vpp\nbOrgd298yLEP27FaTHzz9svp7Lr47ByTycSEAgee9m4ikc/P/QlEJH0U+nEkekVu1F0LrqLZc45P\n3ec42dLBsY/aOfZROy+c97xJpTlUXVXMhydPke8swjnEmmrFBQ7cbd34/Jq2KSKXTqEfx2fz9BP7\nUVktZqZWFDK1opCrKwzOnuvBH8nl9Bk/4bCBwxxgcnkhX7xiMgCt7uHX4I/29c926cYqInLpFPpx\nJNPTv5gJhQVcVeZi5tV9j1vdzZjNtsRf3x/6nV2BlN5fRGQgnciNoyfBtXdGS3H/PNBOjfRFZAQo\n9OOI9vTt1vSEflG+HRMKfREZGQr9OKL3x73YlMqxYLGYKci3q6cvIiNCoR9HbzCc8Myd0VJc4KAn\nEBmRG6+LSHZT6MfREwgnNEd/NE3o7+uf9nantQ4RyXwK/Th6A+G0ncSNii68dvqMQl9ELo1CP46e\nQDjuDVRG24RY6PvTWoeIZD6F/jCCoTCBYJg8x3gJfY30ReTSKPSH4W7rC1lXcV5a68hzWLFZTZxS\nT19ELpFCfxgt/SPrionpDX2TyURRvh13m5/uHs3gEZHUKfSH0dJ/28OKkvw0VwKTy3IJRwwOHW9J\ndykiksEU+sOI3ut20sT0h/7lFX01/PHt5jRXIiKZTKE/DHd/e6e8JL3tHYAip52pFU7e/msrZ7X4\nmoikSKE/jJYzXeQ6rBTm29NdCgBzZ7gIRwzeaDh1SceJRCKxe/gO/BPvfr4ikvm0tPIQDMOgpa2b\nyolOTKb0rLtzvhunu9j66ofs+3MTt355asrrAfl8Pv7j90fp9Jtw2M0UO+1YCAx7P18R+XxQ6A+h\nw9dLbyBMeWn6WztRE4tyuPaqiTR84OV7/+tP3L9sJldUFiX1pRSJGOw90sye+jZC4b5bMJpNcOs/\nVIxW2SIyjij0hxCdrjmpNP0nceGzlsz9d3yRF35v4vB7Hh7aeIA8h4XLJhVw2aQCJpflUZhvZ2Jx\nIZhMGIaBAWCAgcHHp86y+9BJPmnx4bCZuWF6BaGQwZF3W3i9wcNdC65CA32RzzeF/hCiM3cqxslI\n3+/vYvfBdkpKSvn7qU4Kcs2895GXTr/Bux938O7HHQkdx2I2cf20Yi6flI9rYikAEcPgrffcPP/7\n91n/3+eO5scQkTRT6A/BHQv98THSB8jNzSffWQjANc5CSvLDmM02CieU4unw03muF++ZNnp6Q+Tm\n5oLJhAno7j6HyWSmuKiAyyflc67Ti9nIiR33+i+V8/GpDt442sreQye55ctTL7nWSCSCz+cbtK2g\noACzWXMHRNJJoT+Elv4lGMZT6A/FbrNQWeakssxJa14PZrONiWWu2P7ofXmj28LBwcs5WMwmbp5V\nxu8Pt/D0yw1cNWUCl3/h0vo8Pp+P/7vvOB6fiZPuLvLssPTmq7jx2qnj5sS4SDYa8dD/8Y9/TEND\nAwAbNmxg5syZI/0WY+K0twuz2URZcW66SxkTBXk2/sfSa/jptmP8z/99mLV3z+baqyemHNDvftzO\nvoazeDt7Y9uOv/gO86u8rF0+m5wRXLnU3dbNK/v/xut/aebvyp185dpybpzuwmoxj+hvF9HfXsKR\nCJb+Y+q3F8k0Ixr6hw8f5pNPPmHr1q00NjayYcMGtm7dOpJvMWbcbV2UTcjFasmO/6EjkQhXVthY\nNn8qr/zxJI8+8wZTXPlMmphHSYGDCQV2igscFDvtVFaUMHFCLnk5tkHH6O4J8s7fvOw5dJK33nMD\ncNXkIq6b5uJM+1k+afXzx7eb+dR9joe+VnXJv030BsNs3fNXXtn/AeGIgc1q4tiH7Rz7sJ1tr33I\nDV90Uvvfrr3kaaiGYfDnE61s2/seH506S2/QoKzIwfSpuXzjztkjNs01+qUSCkcwDLBZR/ZLSwRG\nOPTffPNNFi9eDMCVV15JZ2cnXV1d5OeP/xbJ+SIRuKIye6ayDDxRfPvcL9DQ2EFTaxdNrV1DvsZh\nM5Ofa8NmNePvDeHrDmL0zQLl6slOvji5gMsm97WU8mwhapdM4/+8dpLdb57koZ/u57/8w99x3TUu\nKsuc5Nit2KxmbFYzdpsFu82C5bzrEAzD4Jw/yCnPOf78npvX3vqE1vYeyibksOSGcgwDIuZc6v/a\nynsft7G3vo2Q6X2W3vxFLptUmPRvLT29Id44eorfv/ExJ062A1CYZ6PIaaW13c/+hl7auhr4xp3X\nXvIXmKfdz+v1H7HrzZO0dgQJRwxy7GZmXlHC7fOvZvbVZWm7T7N8voxo6Hu9XqZPnx57XFJSgsfj\nycjQf3LdAnLSvI7+WIueKM53wtTKMtwtnxIIWXDkF9HlD9LlD9LqbaOrJ0TYsOIP9N1voKu7F6vV\njGuCg4qSXCrLcjEFOnA6Pjt2JBKh65yP2lsvZ9YVRWzZ+wF7D3/C3sOfDFmPzWrGYTNjt5oJhg16\nekME+68tgL7rC6ZfVsjsq4ppO+Mh31nExLIJLLp+Cld8oYh9f/6E1946xWtvnSIvx8qkklzyc619\n9zw29X2xRwyDSMQgHDEwDINIBMKRCO2+AG1newlH+t7vui+W8l+vL6P5TA/OgiLcbd386S9NNDS2\n8eAT+7myspDLv+CkIM+GxWzCYjZjNvdNmw2FI4TDBmEDwuEI4YhBOBwhFDFo9wU45e3G3fbZDXKK\nCxzk5Vjp8PVw5ISXIye8lBTmMOcaF1f/XTETnHZyHVYsFjNWsxmLxYTV0vc3QDAUIRSKEOz/Ewr3\n/d0bDHHW14W/N4y/N0RPIIxhQOmEfAqdDgpy7RTk2yjIs5ObY8WEiej3pKl/UgAmBm0fimEMv/9i\nzj9mIt/R53+RX/CSixzDdN7GgYeIGAZnuwL0BsJcWVmE5XP4m/6oppphGBcdXYXDYQBaWsb3ipFD\nj3Gh9fRJbGc8scedne2EsRPwn4tta2vzYjZb6Ok+O+S2RJ4z0tuSfd2EcBc2YIIVIvY2XDkWJkwo\nJrqKR1ubr+95E3KBIISDtHW2cc7XOeg93383QEFh34j4mjKD/GA7/qCJkMlOJGwQMaC3N0DEAJPF\nRjhkEAwa+INhTGawW6047SZy7GZyzb1MmuhkYmmYc51eOjvO4DvbEXs/K1BV6aO1LUCb34KvO8Jf\n28MJh5HNAjk2mFycw+SJDvKsZ3j1wN/IzXP2f3a4akIbhUaA5k4zJz5o472/JXbsC9/LxASnlQJ7\nkIqJ+XzB1TdI6u6OMLm8iHc+PMexjzz8bv+lLb8hyVm9dAZzZ0xKdxmDRDMzmqGpGNHQd7lceL3e\n2OPW1lbKysoueJ7H0xeWK1euHMm3FxEZMY/+Id0VDM3j8TB1ampTq0c09G+66SZ+/vOfc88993D8\n+HHKy8vJy7vw4qYZM2awZcsWysrKsFjSe9NxEZFMEQ6H8Xg8zJgxI+VjmAwjle7b0J544gmOHDmC\nxWLhhz/8IdOmTRvJw4uIyCUY8dAXEZHx6/N3alpERIak0BcRySIKfRGRLDKqoR8MBlm3bh333nsv\ntbW1NDU1XfCcnTt3UlNTw/Lly9m+fXts++HDh5k3bx779+8fzRKH9OMf/5gVK1awYsUKjh49Omjf\nG2+8wd13382KFSvYtGlTQq8Za6nU//7777N48WK2bNky1uUOkkrtdXV1rFixgpqaGvbu3TvWJQ+S\nbP1+v58HH3yQ2tpali9fnrZ/81Gp/PwBenp6WLx4Ma+88spYlnuBZOs/dOgQN954I7W1tdTW1vLY\nY4+lo+yYVH7+O3fuZOnSpdx1110cOHBg+DcwRtHLL79s/OhHPzIMwzD+9Kc/GQ899NCg/V1dXUZ1\ndbXh8/mMnp4e4/bbbzc6OjqMkydPGg888ICxZs0aY9++faNZ4kUdOnTIuP/++w3DMIwPPvjAuOee\newbtv+1S5x9TAAAFQ0lEQVS224yWlhYjEokY9957r/HBBx/Efc1YSqX+7u5uo7a21vjBD35g/Pu/\n/3s6yjYMI7XaDx48aHzzm980DMMw2tvbjQULFox53VGp1P+f//mfxnPPPWcYhmE0Nzcbt95665jX\nHZVK/VE//elPja9+9avGK6+8MqY1D5Tq/7tr165NR7kXSKX+trY249ZbbzW6urqM1tZW4wc/+MGw\n7zGqI/2Ba/HMnTuX+vr6QfvfeecdZs6cidPpxOFwUFVVRX19PS6Xi6eeegqn0zma5Q1pqDWEAJqa\nmigqKqK8vByTycTNN9/MwYMHh31NJtRvt9t59tlncblcwx161KVS+w033MDPfvYzoG/Vy+7ubow0\nTUpLpf7bbruN1atXA3Dq1CkmTUrfVaCp1A/Q2NhIY2MjCxYsSNvP/lLqHy9Sqf/gwYPMmzePvLw8\nysrK+NGPfjTse4xq6Hu9XkpKSvreyGzGZDIRCoVi+8+cORPbD1BaWorH4yEnJyetKwt6vV6Ki4tj\nj0tKSmJXGns8nkE1R9cXuthrolcej7VU6rdYLDgcjguONdZSqd1sNvfdNAbYvn07CxYsSNua/anU\nH7VixQoefvhh1q9fP3YFnyfV+uvq6tJad1Sq9Tc2NvLAAw9w77338sYbb4xt0QOkUn9zczM9PT08\n8MADrFy5Mu4X2YhdkfvrX/96UE8e+kbyA8UbAaRzhDCc4eoaap8xxLpD6ZBK/eNFMrW/+uqr7Nix\ng82bN492WQlLpv6tW7dy4sQJHn74YXbu3DnapSUkkfp/85vfMHv2bCorK8fdv6dE6p86dSpr1qxh\nyZIlNDU18fWvf529e/ditaZ/wcVE6jcMg/b2djZt2kRzczNf//rX2bdv35CvG7FPdffdd3P33XcP\n2rZ+/Xo8Hg/Tpk0jGAxiGMagH+T5a/W43W6qqqoGHSMdwTncGkLl5eUX1OxyubDZbAmtOzQWUql/\nvEi19tdff51f/vKXPPfcc2lrC0Jq9R8/fpzS0lIqKiq45pprCIfDtLW1DRrVjef6Dxw4QFNTE/v3\n76elpQW73U5FRQVz5479/ZZTqb+8vJwlS5YAMGXKFCZOnIjb7aaysnJsiye1+nNzc7nuuuswm81M\nmTKF/Pz8Yf/9jGoP5aabbmLXrl0A7Nu3jxtvvHHQ/muvvZajR4/i8/no6uqivr6eOXPmxPYbhpGW\nkcNNN93E7t27AS5YQ6iyspJz587R3NxMKBRi//79fOUrXxn2NZlQf1S6R2qp1O7z+airq+Ppp5+m\nsLAwneWnVP+RI0div514vV66u7vTEvip1r9x40a2b9/Otm3buPvuu/n2t7+dlsBPtf7f/va3sZ+/\nx+PhzJkzlJeXZ0z9N910E2+++WZsxB/v38+oLsMQiUTYsGEDJ0+exOFw8Pjjj1NeXs4vf/lLbrjh\nBmbPns3u3bv51a9+hclkora2lttvv50DBw7w3HPP8dFHH1FcXIzL5eJXv/rVaJV5UeevIfTuu+9S\nUFDA4sWLeeutt/jXf/1XAKqrq7nvvvsu+pp0rjuUbP3Hjh3jJz/5Cc3NzdhsNsrLy/n5z38+YneF\nGs3at23bxlNPPcVll10WO0ZdXV3aTogmW39vby/f//73aWlpoaenh7Vr17JgwYK01J5K/QM99dRT\nTJ48mTvvvDMdpQPJ19/V1cW6devw+XwEg0HWrFnD/PnzM6Z+gG3btsXa69/61rdYuHDhkMfX2jsi\nIllEV+SKiGQRhb6ISBZR6IuIZBGFvohIFlHoi4hkEYW+iEgWUeiLiGQRhb6ISBb5/2Hp6zPXZ4Vh\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0073c4c950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(Pi_of_locus.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assert False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myseqs = []\n",
    "for x in hap_df[1].values.tolist():\n",
    "    try:\n",
    "        myseqs.append(list(x))\n",
    "    except TypeError:\n",
    "        myseqs.append([np.NaN]*94)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myseqs = pd.DataFrame(myseqs)\n",
    "myseqs = myseqs.replace({'A':1, 'C':2, 'G':3, 'T':4})\n",
    "myseqs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## jaccard distance will all missing values for one array is 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "scipy.spatial.distance.jaccard(myseqs.ix[22].values, myseqs.ix[20].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myseqs.ix[22].values, myseqs.ix[21].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pdist_mat = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(myseqs, 'jaccard'))\n",
    "pdist_mat = pd.DataFrame(pdist_mat).replace({1:np.NaN})\n",
    "pdist_mat = pdist_mat.dropna(thresh = 50, axis = 0)\n",
    "pdist_mat = pdist_mat.dropna(thresh = 50, axis = 1)\n",
    "pdist_mat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print 'mean', np.nanmean(scipy.spatial.distance.squareform(pdist_mat)) # mean calulated without zeros as self-compairisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "scipy.spatial.distance.squareform(pdist_mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figsize(12,12)\n",
    "sns.clustermap(pdist_mat, square=True, linewidths=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[x for x in 'GGXXXXXXXXXXXXXXXXX']\n",
    "[x for x in np.NaN]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Generate random sequences\n",
    "NUCLEOTIDES = '11111111111111111111111111111111111111111111111111111234'\n",
    "NUCLEOTIDES = np.array([x for x in NUCLEOTIDES])\n",
    "SEQ_LEN = 100\n",
    "NUM_SEQ = 10\n",
    "\n",
    "randseqs = np.random.choice(NUCLEOTIDES, size = (NUM_SEQ, SEQ_LEN), replace = True)\n",
    "pd.DataFrame(randseqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "normalized Hamming distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print 'mean', scipy.spatial.distance.pdist(randseqs, 'hamming').mean() # mean calulated without zeros as self-compairisons\n",
    "scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(randseqs, 'hamming')).shape"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}