(E) Wrangle Genotypes Matrix for Insertions and Deletions with Mixed Allele Freqs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook was created to wrangle all of the insertion/deletion variants in the 31,428 isolates and store the allele frequencies for variants with an alternate allele frequency $\\ge 10\\%$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vcf\n",
    "\n",
    "%matplotlib inline\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.ticker as ticker\n",
    "from itertools import compress\n",
    "from pylab import MaxNLocator\n",
    "import seaborn as sns; sns.set()\n",
    "from matplotlib.colors import LogNorm\n",
    "from matplotlib import gridspec\n",
    "import ast\n",
    "import itertools\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "import fastcluster\n",
    "from sklearn import cluster, datasets\n",
    "import scipy.cluster.hierarchy as hier\n",
    "from sklearn.cluster import KMeans\n",
    "import time\n",
    "import sys\n",
    "import pickle\n",
    "\n",
    "import Bio\n",
    "from Bio.Alphabet import IUPAC\n",
    "from Bio.Blast.Applications import NcbiblastnCommandline\n",
    "from Bio.Blast import NCBIXML\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.SeqFeature import SeqFeature, FeatureLocation\n",
    "from Bio import pairwise2\n",
    "from Bio import SeqIO\n",
    "from Bio.Graphics import GenomeDiagram\n",
    "from Bio.SeqUtils import GC\n",
    "\n",
    "from Bio.Align.Applications import MuscleCommandline\n",
    "from StringIO import StringIO\n",
    "from Bio import AlignIO\n",
    "from Bio.Align import AlignInfo\n",
    "from Bio.Seq import MutableSeq\n",
    "import itertools\n",
    "import gzip\n",
    "\n",
    "import networkx as nx\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] This section of the notebook loads all of the Lineage Calls that were made on the isolates in rollingDB as well as the mean coverage / isolates for all of the isolates in rollingDB. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We runs through the pilon-generated VCF for each isolate (N = 33,875) and scrapes the High Quality INDEL Variant calls made against H37Rv. The DataFrame containing these high quality INDEL calls (with the Reference Position, Reference Allele and Alterante Allele) for each isolate is then stored in a directory containing all of the DFs for all isolates. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Side-Note: Some isolate have ~20x coverage and we require 20 reads at a given position to make an INDEL call"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scrape Luca's DB\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0500152951973\n",
      "0.100030590395\n",
      "0.150045885592\n",
      "0.200061180789\n",
      "0.250076475987\n",
      "0.300091771184\n",
      "0.350107066381\n",
      "0.400122361578\n",
      "0.450137656776\n",
      "0.500152951973\n",
      "0.55016824717\n",
      "0.600183542368\n",
      "0.650198837565\n"
     ]
    }
   ],
   "source": [
    "## genomic data directory\n",
    "rolling_DB_dir = '/n/data1/hms/dbmi/farhat/rollingDB/genomic_data/'\n",
    "\n",
    "#get all folders (each folder corresponds to a different sequenced isolate)\n",
    "isolate_directories = os.listdir(rolling_DB_dir)\n",
    "\n",
    "#iterate through each sequenced isolate\n",
    "isolate_i = 0\n",
    "\n",
    "for isolate_ID in isolate_directories:\n",
    "\n",
    "    #directory that stores files for each sequenced isolate\n",
    "    directory_for_sequenced_isolate = rolling_DB_dir + isolate_ID\n",
    "\n",
    "    #check to see if megapipe successfully ran on sequenced isolate\n",
    "    try:\n",
    "        #existence of a PILON and QUALIMAP directories and corresponding VCF file\n",
    "        if ( 'pilon' in os.listdir(directory_for_sequenced_isolate) ):\n",
    "            \n",
    "            #existence of a REDUCED VCF and GENOME-QUALITY files and existence of FULL VCF file in relevent directories \n",
    "            if ( 'vcf' in list( itertools.chain.from_iterable( [filename.split('.') for filename in os.listdir(directory_for_sequenced_isolate + '/pilon/')] ) ) ) and ( isolate_ID + '_full' in list( itertools.chain.from_iterable( [filename.split('.') for filename in os.listdir(directory_for_sequenced_isolate + '/pilon/')] ) ) ):\n",
    "                \n",
    "                #we have a valid VCF file so megapipe ran successfully, let's keep the variant call information for this sequenced isolate\n",
    "                \n",
    "                #PILON-generated VCF DATA\n",
    "                ########################################################################################################################\n",
    "                \n",
    "                #dictionaries to construct DataFrame from later that hold information about each high-quality SNP found for a specific isolate from rollingDB\n",
    "                ref_dict = {} #reference allele\n",
    "                alt_dict = {} #alternate allele\n",
    "                \n",
    "                #load in VCF file if one was found\n",
    "                vcf_reader = vcf.Reader( open(directory_for_sequenced_isolate + '/pilon/' + isolate_ID + '.vcf' , 'r') )\n",
    "\n",
    "                #iterate through each Variant Call \n",
    "                for record in vcf_reader:\n",
    "                    \n",
    "                    #check to see if the call is a PASS by Pilon (PASS by Pilon ; there's exactly 1 alternate allele ; call supports the alternate allele)\n",
    "                    if (record.FILTER == []) and (len(record.ALT) == 1) and (record.ALT != [None]):\n",
    "                        \n",
    "                        ref_pos = int( record.POS )\n",
    "                        ref_allele = str( record.REF )  #allele on H37Rv\n",
    "                        alt_allele = str( record.ALT[0] ) #alternate allele supported by reads\n",
    "                        \n",
    "                        #check to see if variant is INDEL (either Refeference or Alternate Alleles are > 1 base but not both ; 1 allele is 1bp, the other allele is > 1bp)\n",
    "                        if ( (len(ref_allele) > 1) or (len(alt_allele) > 1) ) and ( (len(ref_allele) == 1) or (len(alt_allele) == 1) ):\n",
    "                        \n",
    "                            #keep only INDELs <= 10bp long (1 allele is 1bp ; the other allele is <= 11bp)\n",
    "                            if (len(ref_allele) <= 11) and (len(alt_allele) <= 11):\n",
    "                                \n",
    "                                #analyze only INDELs with Depth metrics (larger structural variants don't have this)\n",
    "                                quality_metrics = record.INFO.keys()\n",
    "                                if ('MQ' in quality_metrics) and ('DC' in quality_metrics) and ('IC' in quality_metrics) and ('TD' in quality_metrics) and ('DP' in quality_metrics):\n",
    "                            \n",
    "                                    ##### Retrieve Relevant information for filtering quality of Base Call #####\n",
    "                                    # Mean Mapping Quality @ locus\n",
    "                                    MQ = record.INFO['MQ']\n",
    "                                    # Number of Reads w/ Deletion \n",
    "                                    DC = record.INFO['DC']\n",
    "                                    # Number of Reads w/ Insertion\n",
    "                                    IC = record.INFO['IC']\n",
    "                                    # Depth of Valid Reads in Pileup\n",
    "                                    VD = record.INFO['DP']\n",
    "                                    # Total Depth in Pileup\n",
    "                                    TD = record.INFO['TD']\n",
    "\n",
    "                                    ### Filtering Criteria\n",
    "\n",
    "                                    #---> Mean Mapping Quality > 30\n",
    "                                    #---> Number of High Quality Reads >= 20\n",
    "                                    if (MQ > 30) and (VD >= 20): #INDEL passed filtering criteria!\n",
    "\n",
    "                                        # Calculate Alternate Allele Frequency\n",
    "                                        ########################################################################\n",
    "\n",
    "                                        #calculate INDEL (Alternate Allele) Frequency from Deletion/Insertion Count & Total Depth (all metrics include ALL reads)\n",
    "                                        alt_allele_frequency = float(max(DC , IC)) / float(TD)\n",
    "\n",
    "                                        ########################################################################\n",
    "\n",
    "                                        # use INDEL calls that pass the following filters\n",
    "                                        #---> alternate (INDEL) allele frequency >= 10%\n",
    "                                        #---> Reads supporting INDEL variant >= 10 reads\n",
    "                                        if (alt_allele_frequency >= 0.10) and (float(max(DC , IC)) >= 10):\n",
    "\n",
    "                                            # Store good-quality SNV in SNV dictionaries\n",
    "                                            ref_dict[ref_pos] = ref_allele\n",
    "                                            alt_dict[ref_pos] = alt_allele\n",
    "\n",
    "                #construct DataFrame from dictionaries\n",
    "                INDEL_DF = pd.DataFrame( index = pd.Series( ref_dict.keys() ) )\n",
    "                INDEL_DF['ref'] = pd.Series(ref_dict)\n",
    "                INDEL_DF['alt'] = pd.Series(alt_dict)\n",
    "\n",
    "                #make reference-position a column and re-index, sort according to reference position\n",
    "                INDEL_DF['pos'] = INDEL_DF.index\n",
    "                INDEL_DF.sort_values(by = 'pos' , ascending = True , inplace = True)\n",
    "                INDEL_DF.reset_index(inplace = True , drop = True)\n",
    "\n",
    "                #pickle DataFrame holding High Quality SNVs from isolate\n",
    "                INDEL_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB/' + isolate_ID + '_INDELs.pkl')\n",
    "                ########################################################################################################################\n",
    "                \n",
    "                #keep track of progress   \n",
    "                isolate_i += 1\n",
    "                if isolate_i % np.ceil(0.05*len(isolate_directories)) == 0:\n",
    "                    print float(isolate_i) / float(len(isolate_directories))\n",
    "        \n",
    "    except OSError: #hit some file that is not another directory with genomic data \n",
    "        continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2]  The purpose of this section of the notebook is generate the master catalogue of (high-quality) INDEL calls made in rollingDB by iterating over all of the HQ SNV DataFrames for each isolate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vcf\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from pylab import MaxNLocator\n",
    "import ast\n",
    "import itertools\n",
    "import time\n",
    "import sys\n",
    "import itertools\n",
    "import gzip\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list of pickles DataFrame files for INDELs for each isolate in rollingDB\n",
    "INDEL_DF_list = os.listdir('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def store_INDEL_info_from_INDEL(INDEL_DF_row):\n",
    "\n",
    "    ref_allele , alt_allele , ref_pos = list(INDEL_DF_row)\n",
    "    \n",
    "    INDEL_i_key = ref_allele + '_' + str(ref_pos) + '_' + alt_allele\n",
    "\n",
    "    #INDEL is not documented at this Reference Position\n",
    "    if INDEL_i_key not in ref_pos_dict.keys():\n",
    "\n",
    "        #store a new row for INDEL at this reference position; ref position, ref allele, alt allele \n",
    "        ref_pos_dict[INDEL_i_key] = ref_pos\n",
    "        ref_allele_dict[INDEL_i_key] = ref_allele\n",
    "        alt_allele_dict[INDEL_i_key] = alt_allele\n",
    "        num_isolates_dict[INDEL_i_key] = 1\n",
    "\n",
    "    #a INDEL is already documented at this Reference Position\n",
    "    else:\n",
    "\n",
    "        #increment # of isolates that have an indel here\n",
    "        num_isolates_dict[INDEL_i_key] = num_isolates_dict[INDEL_i_key] + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0500056028687\n",
      "0.100011205737\n",
      "0.150016808606\n",
      "0.200022411475\n",
      "0.250028014343\n",
      "0.300033617212\n",
      "0.350039220081\n",
      "0.400044822949\n",
      "0.450050425818\n",
      "0.500056028687\n",
      "0.550061631555\n",
      "0.600067234424\n",
      "0.650072837293\n",
      "0.700078440161\n",
      "0.75008404303\n",
      "0.800089645899\n",
      "0.850095248767\n",
      "0.900100851636\n",
      "0.950106454505\n"
     ]
    }
   ],
   "source": [
    "#create dictionaries to store INDEL info for each isolate\n",
    "ref_pos_dict = {}\n",
    "ref_allele_dict = {}\n",
    "alt_allele_dict = {}\n",
    "num_isolates_dict = {}\n",
    "\n",
    "#iterate through each sequenced isolate\n",
    "isolate_i = 0\n",
    "\n",
    "#iterate through all DataFrames for each isolate in rollingDB\n",
    "for pickled_INDEL_DF in INDEL_DF_list:\n",
    "    \n",
    "    #load SNV DF for isolates from pickles object\n",
    "    INDEL_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB/' + pickled_INDEL_DF)\n",
    "\n",
    "    #extract important info from SNV df\n",
    "    _ = INDEL_DF.apply(lambda row: store_INDEL_info_from_INDEL(row) , axis = 1 )  \n",
    "    \n",
    "    #keep track of progress   \n",
    "    isolate_i += 1\n",
    "    if isolate_i % np.ceil(0.05*len(INDEL_DF_list)) == 0:\n",
    "        print float(isolate_i) / float(len(INDEL_DF_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#construct DataFrame from dictionaries\n",
    "INDEL_DF = pd.DataFrame( index = pd.Series( ref_pos_dict.keys() ) )\n",
    "INDEL_DF['indel_id'] = INDEL_DF.index\n",
    "INDEL_DF['ref_pos'] = pd.Series(ref_pos_dict)\n",
    "INDEL_DF['ref_allele'] = pd.Series(ref_allele_dict)\n",
    "INDEL_DF['alt_allele'] = pd.Series(alt_allele_dict)\n",
    "INDEL_DF['num_isolates'] = pd.Series(num_isolates_dict)\n",
    "\n",
    "#sort according to reference position\n",
    "INDEL_DF.sort_values(by = 'ref_pos' , ascending = True , inplace = True)\n",
    "\n",
    "INDEL_DF.reset_index(inplace = True , drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>indel_id</th>\n",
       "      <th>ref_pos</th>\n",
       "      <th>ref_allele</th>\n",
       "      <th>alt_allele</th>\n",
       "      <th>num_isolates</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>31</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           indel_id  ref_pos  ref_allele alt_allele  num_isolates\n",
       "0  ACCGACGAAG_313_A      313  ACCGACGAAG          A             3\n",
       "1         TC_1549_T     1549          TC          T             3\n",
       "2        TAA_1552_T     1552         TAA          T             1\n",
       "3         TA_1552_T     1552          TA          T           479\n",
       "4         T_1552_TA     1552           T         TA            31"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>indel_id</th>\n",
       "      <th>ref_pos</th>\n",
       "      <th>ref_allele</th>\n",
       "      <th>alt_allele</th>\n",
       "      <th>num_isolates</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>59502</th>\n",
       "      <td>C_4411205_CA</td>\n",
       "      <td>4411205</td>\n",
       "      <td>C</td>\n",
       "      <td>CA</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59503</th>\n",
       "      <td>AG_4411208_A</td>\n",
       "      <td>4411208</td>\n",
       "      <td>AG</td>\n",
       "      <td>A</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59504</th>\n",
       "      <td>C_4411423_CTG</td>\n",
       "      <td>4411423</td>\n",
       "      <td>C</td>\n",
       "      <td>CTG</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59505</th>\n",
       "      <td>CCTG_4411459_C</td>\n",
       "      <td>4411459</td>\n",
       "      <td>CCTG</td>\n",
       "      <td>C</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59506</th>\n",
       "      <td>AC_4411486_A</td>\n",
       "      <td>4411486</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             indel_id  ref_pos ref_allele alt_allele  num_isolates\n",
       "59502    C_4411205_CA  4411205          C         CA             1\n",
       "59503    AG_4411208_A  4411208         AG          A             1\n",
       "59504   C_4411423_CTG  4411423          C        CTG             3\n",
       "59505  CCTG_4411459_C  4411459       CCTG          C             4\n",
       "59506    AC_4411486_A  4411486         AC          A             1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_DF.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507, 5)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output Pickle files for downstream analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_DF.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB_catalogue/all_HighQuality_INDELs_found_in_rollingDB_catalogue.pkl' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3] The purpose of this section of the notebook is to analyze the catalogue of (high-quality) INDEL calls made in rollingDB. For population genetic analyses, we will retain only INDELs that are present in at least 1 isolate. We then construct a skeleton genotype DF (rows = Reference Position , columns = isolates) for all of the remaining INDELs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct a series for the isolate annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "HQ_INDEL_per_isolate_files = os.listdir('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB/')\n",
    "isolate_ids_list = [isolate_id.split('_INDELs')[0] for isolate_id in HQ_INDEL_per_isolate_files]\n",
    "isolate_id_series = pd.Series(isolate_ids_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_id_series.head() #index = column indices of genotype matrix, values isolate IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35696,)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_id_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remove 2 bad isolates from isolate annotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remove:** Delete ID for isolate **IDR1400002533** since Full GZipped VCF looks corrupt and remove from list accordingly!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series = isolate_id_series[isolate_id_series != 'IDR1400002533']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remove:** Delete ID for **Peru4994** since F2 calculations keeps failing for this isolate (lineage-defining SNP depths df is empty) **Think this is because sample has very low coverage**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series = isolate_id_series[isolate_id_series != 'Peru4994']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series.reset_index(inplace = True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_id_series.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35694"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_id_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_isolate_annotation.pkl' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load in catologue of INDELs detected within rollingDB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#INDEL data for all isoaltes in rollingDB \n",
    "INDEL_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/all_HighQuality_INDELs_found_in_rollingDB_catalogue/all_HighQuality_INDELs_found_in_rollingDB_catalogue.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>indel_id</th>\n",
       "      <th>ref_pos</th>\n",
       "      <th>ref_allele</th>\n",
       "      <th>alt_allele</th>\n",
       "      <th>num_isolates</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>31</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           indel_id  ref_pos  ref_allele alt_allele  num_isolates\n",
       "0  ACCGACGAAG_313_A      313  ACCGACGAAG          A             3\n",
       "1         TC_1549_T     1549          TC          T             3\n",
       "2        TAA_1552_T     1552         TAA          T             1\n",
       "3         TA_1552_T     1552          TA          T           479\n",
       "4         T_1552_TA     1552           T         TA            31"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507, 5)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')\n",
    "plt.rcParams['lines.linewidth']=3\n",
    "plt.rcParams['axes.facecolor']='1.0'\n",
    "plt.rcParams['xtick.color']='black'\n",
    "plt.rcParams['axes.grid']=True\n",
    "plt.rcParams['axes.edgecolor']='black'\n",
    "plt.rcParams['grid.color']= '1.0'\n",
    "plt.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAFxCAYAAAAmtUHWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XtcFGXfP/DPQmiyYUGUyUoIy0kwQ9RVK1NTUVD0NvBQpmmUeSoRtDKr2ycT0lBvn8dOWuJjdyolEplI2MEzeUjkRghEVwWVTBATFlZlmd8f/pjHFViWgT3h5/168Xq51+zMfC9mi89ec82MTBAEAUREREQS2Fm6ACIiIrJdDBJEREQkGYMEERERScYgQURERJIxSBAREZFkDBJEREQkGYME2aS+fftCoVDg4MGDli4FOp0O0dHR6NatGxQKBd5//32T7zM6OhoKhQIrVqww+b6shTUdc1Npi8f14MGDUCgU6NatW6PvuRuObVvGIEGS1f3HHxwcjOrqagDAiRMnoFAooFAoLFyd+ezYsQPffvst7O3t8dJLL6Ffv3713mPM/0xNqW7/ffv2Nfk+7vwZNmyYyfbZlLoaTpw4AeD//lB36dJFbAMghsC6P2QrVqwQ1/Xy8sJjjz2GsWPHYsOGDaitrRXXu/19t/+89957AJp/3AcOHIioqCgEBwc3u6/FxcUN1mKpz1xzTJw4EVFRUejcubOlSyEJ7rF0AWT7Ll26hI0bN+LVV1+1dCktcvPmTTg4ODR7PbVaDQB45plnsGTJktYuy+a0a9cOkydPFl9b4x8HQRCQkJCADRs2GHyfh4cHBg0aBLVajf379+Pw4cPYv38/1q5dCzs7O733DR06VHzdv39/SXWNHTsWY8eOlbTu7SZPnox27doBADp06NDi7ZnavHnzLF0CtQBHJKjFZDIZPvnkE3FU4k5134yKi4sB/N+3uOjoaABAUlISFAoFhg4disWLF8PHxweDBg3CiRMnsHz5cvj7+6N///7Ys2dPvW3n5ORg2LBh8PX1RVRUFK5cuSIuO3z4MCIjIxEQEIDg4GDExMSIy2//9rZx40YEBwfjueeea7D+4uJiTJ8+HT179kRAQAAiIyNx7NgxsS8fffQRACA5ORkKhQJJSUlN/s4EQUB8fDx69+4NT09PBAUF4fnnnxfrKysrw/z589GnTx/4+flh1KhR+PXXXxvd3t69exESEgJ/f394eHhApVIhISEBwK1vxePGjQMAnD9/Xm/EqLq6GnFxcXjyySfh7e2N4cOHIz09XW+7I0aMgLe3N/z9/TF8+HCkpaUZ7Nu9996L999/X/yZOXMmgFtBbeLEiQgKCkLXrl3RrVs3TJ06FRcuXBDXvXjxIubOnQuVSgUvLy8MHDgQx48f19v+iRMnEBoaCm9vb0yePBlXr15t8vd9J5lMhl27diE7O9vg+/z9/REXF4ctW7YgMTERALBz505s37693vtu73NoaGizawLqn9rIycnB2LFj4efnBx8fHzzzzDP43//93ya3s3DhQrGWhQsXiu2GPstA/VMMdf9tRkZGAtAf2frv//5v9OjRAz169MCnn34qbuPatWuYMWMG/Pz8MHToUOTk5DRZ7537jYyMhEKhQHx8PJ599lkolUqMGTMG58+fN+K3SObGIEEtNmrUKJSWljb57a4p+fn5yMrKgo+PDwoLCzFu3Djs2LEDwcHBKCoqQkxMTL11Vq5cicceewwuLi5IT0/Hm2++KW5rwoQJ+M9//oNBgwbBy8sLSUlJePXVV3HnXeGXLVuGwYMHo3fv3vW2X1VVhfHjx2PHjh3w8vLCk08+iczMTIwfPx5nz55FcHAwevbsCQDw8fFBVFQUfH19m+zrvn37sGbNGtjb22PixIno27cv8vPzodFoUFtbi6lTp2Lz5s1wcXFBSEgIcnJyMGXKFBw5cqTB7ZWUlMDFxQWjR49GREQENBoNVq1ahdTUVHTu3BlhYWEAgPvuuw9RUVGIiooCAMTGxuLjjz+Gk5MTwsLCcPHiRbz88svi/9DnzZuHvLw8hIWFISwsDHZ2digoKDDYN61Wi/fee0/8SUlJAQDU1tbir7/+wsCBA/H888/j0Ucfxa5du/DGG28AuBVqxo8fj61bt6J9+/aIiIjA/fffj0uXLultf/ny5fD19UX79u3xyy+/YO3atU3+vu80fPhwODg4iGHLGMOGDROP9a5du/SW5efn6/U5Kyur2TU15N1338Xhw4fx9NNP4x//+Afuv/9+o/4wx8fHi7XUBaCmPsvNcf78eWzbtg0qlQplZWVYunSpODL37rvvYvv27XByckKPHj2watWqZve7zieffAI3Nze4uLjg6NGjWLZsmeRtkenw1Aa12OjRo1FQUIBPP/0UvXr1krwdR0dHbNmyBVlZWRg3bhyuXbuG1NRUuLm5wc/PD3/++SfKysrw4IMPiussWLAAL7/8MnJzcxESEoK0tDRoNBps3LgRN27cQGBgIFxdXeHq6opjx47h4MGDOH36NNq3by9u4/PPP8dTTz3VYE0//fQTioqK4OHhga1bt0ImkyEqKgrp6enYvHkzFi5ciGPHjiErKwtBQUFGT7SsqakBAHTt2hXh4eHw9fXFgw8+CEEQkJ2djWPHjkEulyMlJQWOjo5wcXHBF198gQ0bNqBPnz71tjdu3Di4uroiJycH5eXl8PDwwNWrV3HgwAGMGTMG06ZNQ1paGh544AGxxrKyMqSmpsLOzg69e/eGnZ0d/Pz8kJmZia+++gpPPPEEampq0L59e3G0w8vLq14Qu9ONGzfw5Zdf6tU2duxYtG/fHl988QV27dqFy5cvw9/fHydOnEBmZiZqa2vx888/48yZM+jUqRMyMjLEIfmbN2/qbT82NhYzZ85EQkICVq1apTfXwVgKhQITJkzAv//9b71v5E3p0qULsrKyUFpaqtd+7tw5vT4HBgaKoaMl6j4nQ4YMQVBQEJRKpd4plcZ89dVX4r/79++PadOmGfVZNpa9vT2++eYbPPzww1CpVLhw4QLy8vLg4eGB77//HgDw8ccfo2/fvujevTvefffdZvb8lsmTJyMuLg5JSUmIiYlBbm6upO2QaTFIUIvZ2dlh3rx5mDlzpvjtxxCdTtdgu7u7Ozp06ICOHTuKbUqlEvb29uLrqqoqvSDh4+MDAPD29hbbSkpKxNMoWVlZ9b4dnjlzBv7+/uLrhv4w16kbSlUqlZDJZHr7un1IvrkGDhyIF198EcnJyeJph6CgIKxfv16s3c3NDY6Ojnr7bGxo96233sLXX39dr72srKzRGur2U1tbW++41X1DXbZsGZYsWSLOf3F2dsbSpUsxZsyYRrfbsWNH/PHHH/XaDx06hHHjxtU7/tevX0dlZaVYj7+/v955/TvnrQQGBor7AW59JqR4/fXX8e233zbrCom637+rq6te+/Dhw7F+/XpJdRjyz3/+EwsXLsT8+fMhCALkcjnmz5+P6dOnG1wvLy8P999/v16blM/y7RNLb/fQQw/h4YcfBnDrOFy4cAEajQZXrlzBjRs3xP0AgJeXlzFdbVD37t3FfQCARqORvC0yHZ7aoFYRHh6Obt261Tt3DPzfZK+KigoAaHRovKFvWreHiIYUFhYCAE6dOiW2de7cGe7u7gCA6dOn48KFC+JPZmZmvasIbh+duFOXLl0A3JpQWfdN/PTp0wDQoitTdDodli5dij/++AMHDhxAZGQkjh8/jk2bNom1X7x4UZx3UrfPunruVPd7X7VqFYqLizFlyhQAEGuu+z3e/oehbj/t2rVDTk6O+Ds6e/YsvvjiCwDA4MGDceDAAeTk5GDt2rUoLy+XPLy8Y8cO6HQ6DB48GKdOncIPP/wgLhMEQawnPz9fb75N3bfyOvfcc+v7T90fQ6kUCgWee+457N69W/xsGrJr1y4xlJrrSpQePXrgp59+Ql5eHrZu3YqamhrEx8fX+50Yw5jPcl1wraysBHDrWDSk7hgA+sfBxcVFnORZt+26Ux5S1H1uW3qsybQ4IkGtQiaTYd68eQ1+U+revTuOHDmCd955B0qlEj/++GOr7fejjz5CXl6eeE4/NDQUcrkckyZNwqZNm/Dll1+iqKgILi4uKCwsxNGjR5s1YWvo0KFwd3fH2bNnMW7cODg7O2Pnzp249957MXHiRMl1Hz16FPPmzUOvXr3wwAMP4OjRowCA+++/H48//jh69uyJrKwsjB07Fr6+vkhNTYVMJsOLL77Y4PZcXV1x7do1rF+/Hnv37sXOnTv1lru5uQG4NVozf/58eHp6Yvbs2QgPD8f27dsxatQoDBgwAOXl5Th8+DAmT56M2NhYDB8+HO7u7nBzc8PFixcBQG/EqCF1cyTqdOjQAQsXLsRDDz0E4NYo0TvvvIPffvtNb70hQ4bA09MTZ86cwfDhw9GvXz+cPn0a06dPx/Dhw5vx2zXenDlzsGXLFmi12gaX5+fnY9GiRVCr1di3bx+AW5+x8PDweu+7vc/dunXTm7yr0WgwatQo8fV9992HLVu2NFnf1KlTodPp0LVrV1y7dg3Xr1+Hs7NzkwG7IcZ8lrt3746TJ09i+fLl2L9/v94pEmPY29sjPDwcycnJmD17Np5++mns2LGj2bWSbeGIBLWasLAwcdj5dh988AG6deuG3NxclJSUYMKECa22z5iYGOTk5KCsrAwhISFYvnw5gFvD31u2bEHfvn3x22+/ITU1FZWVlZg9e3aztu/o6IikpCSEhYXh1KlT2L9/P/r374+kpCR4enpKrvuRRx6Bp6cn9u/fj02bNqG6uhqTJ0/GpEmTYGdnhw0bNmDChAkoLS1Feno6unfvjsTERKhUqga3l5CQAG9vbxQWFqKyshIvvPCC3nJ3d3fMmDEDHTt2xObNm5GcnAzg1lUnc+bMgUwmw7fffoujR4+iV69eGDRoEABgwIABOH36NLZu3YojR46gf//+TU5QrJsjUfezceNGAMC0adMwYsQIXL9+Hb/99htef/11vfU6dOiApKQkREREoLq6Glu3bkVpaSk6deok5VdslM6dO2PSpEmNLj937hw2bdqEEydOoE+fPoiLi6t36Wfd+27v852TMXU6nXiaLSsrq8mrRer0798fly5dwrZt2/Dzzz8jKCgIn376qaRv6MZ8lt944w306tUL586dQ05ODqZOndrs/SxZsgQjR47EtWvXcPz4ccyZM6fZ2yDbIhOamjllYcXFxQgLC4Ofnx+AWxPjbj9HTkRERJZjE6c2+vXrh3Xr1lm6DCIiIrqDTZzaOHr0KMaOHYv4+PgmLz0jIiIi8zFbkEhMTERoaCg8PT3FOxrWKS8vR1RUFLy9vaFSqcQb2ADAww8/jAMHDmDbtm0oKytr8q56REREZD5mCxKdOnXC3LlzG5xot2jRIjg4OCA7Oxtr1qzBwoULxUsE27dvD0dHR8hkMoSGhiIvL89cJRMREVETzBYkwsLCMGLECDg7O+u1V1VVIS0tDQsWLIBcLodKpcKwYcPEWeW3X999+PBhdO3a1VwlExERURMsPkdCrVbD3t5evAsacOvSvboRiczMTIwYMQJjx45FSUlJs5+MxxuZEBERmY7Fr9rQaDRwcnLSa3NychJvhRoSEoKQkBBLlEZERERNsPiIhFwur3d72oqKCsjlcgtVRERERMayeJDw8vKCTqfTux97Xl6eeAMqIiIisl5mCxI1NTXQarWora2FTqeDVqtFTU0NHB0dERoaioSEBFRVVeHIkSPIyMhARESEuUojIiIiicx2i+wVK1Zg5cqVem0xMTGIjY1FeXk5YmNjsXfvXjg7O+Ptt99u9qTKxshkMt7EioiIyESs/lkbLcUgQUREZDoWnyNBREREtotBgoiIiCRjkCAiIiLJGCSIiIhIMgYJIiIikoxBgoiIiCRjkJBAp9NJWkZERNTWWPyhXbbI3t4er7zyit5tvYFbt/tet26dhaoiIiIyPwYJidRqNfLz8y1dBhERkUXx1AYRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQMEkRERCQZgwQRERFJxiBBREREkjFIEBERkWQ2EyS+++47PPbYY5Yug4iIiG5jE0GitrYWO3bsgJubm6VLISIiotvYRJBISUnByJEjYWdnE+USERHdNcz2lzkxMRGhoaHw9PREdHS03rLy8nJERUXB29sbKpUKKSkp4jKdToft27dj9OjR5iqViIiIjHSPuXbUqVMnzJ07F7t374ZWq9VbtmjRIjg4OCA7Oxu5ubmYMmUKAgIC4Ofnh+TkZISHh3M0goiIyAqZ7a9zWFgYRowYAWdnZ732qqoqpKWlYcGCBZDL5VCpVBg2bBiSk5MBAIWFhdi6dSsmTZqEM2fO4N133zVXyURERNQEs41INEatVsPe3h5KpVJsCwwMRGZmJoBboxV1QkNDsWTJErPXSERERA2z+PkCjUYDJycnvTYnJydoNJp67925c6e5yiIiIiIjWDxIyOVyVFRU6LVVVFRALpdbqCIiIiIylsWDhJeXF3Q6HdRqtdiWl5cHPz8/C1ZFRERExjBbkKipqYFWq0VtbS10Oh20Wi1qamrg6OiI0NBQJCQkoKqqCkeOHEFGRgYiIiLMVRoRERFJZLYgsXr1aiiVSqxZswbbtm2DUqnE6tWrAQBxcXHQarXo0aMHZs2ahfj4eI5IEBER2QCZIAiCpYswJZlMBlN0cciQIcjPz9dr8/f3x88//9zq+yIiIrJWFp8jQURERLaLQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkMypIXL16FRcuXAAAXL58GatXr8aGDRtw8+ZNkxZHRERE1u0eY940e/ZsFBQU4OjRo5g2bRqOHz8OmUyGwsJCLF261NQ1EhERkZUyakQiLy8Pffr0QVlZGY4fP47IyEg8+uijSE9PN3V9REREZMWMChJ///03XFxccPLkSchkMsTExGDQoEEoLy83dX1ERERkxYw6tfHggw/il19+gVqtxn333YdHH30UV69eRceOHU1dHxEREVkxo0Yknn32WRQXF2Pfvn2IiIgAAPz+++8ICAgwaXFERERk3YwakVi4cCGeeuop1NTUYODAgbh58yZWrVqFzp07m7o+IiIismJGBQkAGDBggPhvOzs7/Prrrzhz5gzWrVtnksKIiIjI+hkdJO6UmZmJ48ePt2YtREREZGN4Z0siIiKSjEGCiIiIJDN4auPll19udJlarW71YoiIiMi2GAwSTd25UiaTtWoxREREZFsMBomYmBhz1UFEREQ2iEGCiIiIJDM42TIjIwMFBQUAbj1vo+5R4gBw6NAhfPnll6atjoiIiKyawSDx0ksvYcuWLQCAlStXol+/fuKytLQ0LF682KTFERERkXVr8vJPQRDMUQcRERHZIMl3tjSXy5cvIyoqCg4ODrCzs8OaNWvQqVMnS5dFREREMCJInDhxAl9++SVyc3MBAOvXr4cgCOJrU3NxccF3330HOzs7JCUlYfPmzYiOjjbLvomIiMiwJoPEoUOHcOjQIQC3TnP885//FP9tjvtI2Nvbi//WaDTw8/Mz+T6JiIjIOAaDxO2TK1sqMTER33zzDfLz8zFmzBj861//EpeVl5dj/vz52LNnD1xcXLBw4UKMHTtWXH7ixAm89dZb+Pvvv7F58+ZWq4mIiIhaxmCQ2Lp1a6vtqFOnTpg7dy52794NrVart2zRokVwcHBAdnY2cnNzMWXKFAQEBIijD927d8cPP/yA77//Hv/zP/+DZcuWtVpdREREJJ1Rky1LSkqwe/dulJaWwtXVFYMGDULnzp2btaOwsDAAQHZ2NkpKSsT2qqoqpKWl4eeff4ZcLodKpcKwYcOQnJyMt99+G9evX0f79u0BAB07dkSHDh2atV8iIiIynSaDxNq1a/Hhhx/i5s2bYpuDgwMWLVqEqKioFhegVqthb28PpVIptgUGBiIzMxMAkJOTg7i4ONjZ2eHee+/FihUrWrxPIiIiah0Gg8SePXvw/vvvAwDs7Ozg4uKCK1eu4MaNG1i8eDF8fX0xYMCAFhWg0Wjg5OSk1+bk5ASNRgMA6N27N7Zt29aifRAREZFpGLwhVWJiIgBg5syZKCgowPHjx1FQUIAZM2ZAEARxeUvI5XJUVFTotVVUVEAul7d420RERGRaBoPEiRMn0L17dyxatEicm9ChQwe88847CAwMRE5OTosL8PLygk6ng1qtFtvy8vJ4mScREZENMBgkrl69ih49ejS47PHHH0d5ebnRO6qpqYFWq0VtbS10Oh20Wi1qamrg6OiI0NBQJCQkoKqqCkeOHEFGRgYiIiKa1xMiIiIyO4NzJLRaLdRqdYNP+VSr1bh+/brRO1q9ejVWrlwpvt62bRtiYmIQGxuLuLg4xMbGokePHnB2dkZ8fDxHJIiIiGyATDDwVK4uXbo0evfKujtbFhcXm6y41iCTyUzy4LEhQ4YgPz9fr83f3x8///xzq++LiIjIWpntzpZERETU9pjtzpZERETU9hicbElERERkiMERCXd3d4Mry2QyFBUVtWpBREREZDsMBglTTFIkIiKitsNgkPjxxx/1XldVVSExMRHbt2+HIAgIDAw0aXFERERk3QwGibqgoNVqsWHDBnz22WcoLS2Fv78/YmJixCd6EhER0d3JYJC4fv06Nm7ciE8++QSlpaXw8fHBkiVLEB4ebq76iIiIyIoZDBJPPPEE/vrrL8hkMoSHh2P06NGws7NDRkaG+J6QkBCTF0lERETWyWCQuHTpknhnyO3bt2P79u16y3nVBhER0d3NYJBQKBSN3iKbiIiIyGCQOHTokLnqICIiIhvEO1sSERGRZEZd/mlIbm5uqxVDREREtsVgkPj7778Nrsz5E0RERHc3g0Hi22+/NVcdREREZIMMBon+/fubqw4iIiKyQZxsSURERJIxSBAREZFkDBJEREQkWaNB4oMPPsCuXbsA3Jp0mZ2dbbaiiIiIyDY0GiQ+++wz7N+/HwAwb948fPfdd2YrioiIiGxDo1dtdOzYESkpKeK9JPbs2YOYmBi998hkMqxYscK0FRIREZHVajRIhISEYOvWrdi6dStkMhlOnjyJkydP6r2HQYKIiOju1miQWLlyJYYMGYLTp08jISEBPXv2xODBg81ZGxEREVm5RoOEnZ0dwsPDAQC1tbUIDg7GoEGDzFUXERER2QCDd7asExMTg6KiIixbtgznz5+Hu7s7JkyYAA8PD1PXR0RERFbMqCCRlZWFCRMmoLq6GoIgAAC++OILJCUloWfPniYtkIiIiKyXUTekiouLQ1VVFYYPH47Zs2djxIgRqKqqQnx8vKnrIyIiIitm1IhEXl4ewsLCsHbtWrFt+vTpOHDggMkKIyIiIutn1IhE+/btceXKFb22K1euoH379iYpioiIiGyDUSMSvXv3xs6dO/Hkk09CqVRCrVbj3LlzCAsLM3V9REREZMWMChLvvfcecnJycO7cOZw7dw4A4O7ujnfeecekxREREZF1MypIdOnSBb/88gt++uknnD9/Hl26dMHQoUPRoUMHU9dHREREVsyoIAEAHTp0EG9QRURERAQYOdmSiIiIqCEMEkRERCQZgwQRERFJ1mSQqKmpwfjx47F8+XJz1ENEREQ2pMkgcc899+D06dO4fPmyOeohIiIiG2LUqY3o6Gikp6cjMzMTN2/eNHVNREREZCOMuvxz4cKFkMlkGD9+vF67TCZDUVGRSQojIiIi62f0fSTqHh/eVBsRERHdPYwKEr/99pup62gTXF1dodPpYG9v3+ByQ8uIiIhskdG3yAaA06dP4+zZsxgyZIhJi7JVTk5OsLe3xyuvvAK1Wq23zMvLC+vWrbNQZURERKZhVJAoLy/HjBkzcPDgQchkMhw4cABPPvkk5syZgzfeeMPUNeLo0aP4r//6L7Rr1w6dOnXC6tWr4eDgYPL9SqVWq5Gfn2/pMoiIiEzOqKs2PvjgAxw4cADt2rWDIAhwd3dHcHAwfvrpJ1PXB+DWiMg333yD5ORkPProo/jxxx/Nsl8iIiIyzKggsXv3bvTv3x8vvPCC2Obr62u2KzYeeeQR8Umj99xzD+zseENOIiIia2DUX2StVotOnTrptV25cqXZpxcSExMRGhoKT09PREdH6y0rLy9HVFQUvL29oVKpkJKSUm/9oqIi/Prrrxg6dGiz9ktERESmYdQciW7duuGnn36Cr68vAGDJkiXYtWsX+vXr16ydderUCXPnzsXu3buh1Wr1li1atAgODg7Izs5Gbm4upkyZgoCAAPj5+QEAKioqEB0djdWrV6Ndu3bN2i8RERGZhlEjEm+88QZu3LiBY8eOQRAEfP7557Czs0NsbGyzdhYWFoYRI0bA2dlZr72qqgppaWlYsGAB5HI5VCoVhg0bhuTkZAC3nvcxa9YsxMTEwNvbu1n7JCIiItMxKkioVCqkpaVhypQpeOaZZ/Diiy9ix44dUKlUrVKEWq2Gvb09lEql2BYYGIiCggIAwHfffYesrCz861//QmRkJFJTU1tlv0RERNQyRt/Z0t/fH++++y5KSkrQuXNncfJja9BoNHByctJrc3JygkajAQBERkYiMjKy1fZHRERErcOoEYmrV69i+vTp8PX1xcCBA+Hr64vp06ejvLy8VYqQy+WoqKjQa6uoqIBcLm+V7RMREZFpGBU/GPgeAAAa+ElEQVQk5s+fj7S0NAiCIP7s3LkTCxYsaJUivLy8oNPp9O4GmZeXJ060JCIiIutkVJDYt28fHn30UWRkZECtVuPHH39Ely5dsG/fvmbtrKamBlqtFrW1tdDpdNBqtaipqYGjoyNCQ0ORkJCAqqoqHDlyBBkZGYiIiJDUKSIiIjIPo4KEh4cHnnjiCQQEBKBdu3YIDAzEE088ga5duzZrZ6tXr4ZSqcSaNWuwbds2KJVKrF69GgAQFxcHrVaLHj16YNasWYiPj+eIBBERkZVrdLLl7U/8jIyMxKpVqxAYGAilUolTp04hLS0NMTExzdpZbGxso5eMOjs7Y/369c3aHhEREVlWo0EiMjISMplMr+29994T/y0IApYsWYJXXnnFdNURERGRVWs0SCgUinpBgoiIiOh2jQaJQ4cOmbMOIiIiskFG35AKuPWgrurqar02hULRqgURERGR7TAqSOzduxexsbH4888/9dplMpnZHiVORERE1seoIPHWW2+hpKSkXrsgCK1eEBEREdkOo4JEeXk5Bg0ahM8//5y3rbYyOp0O9vb2zV5GRETUGowKEtOmTcPOnTtRUlICpVLJqzmsiL29PV555RW924sDt247vm7dOgtVRUREdwujgkRYWBg2bNiAwYMH67VzjoR1UKvVyM/Pt3QZRER0FzIqSLz22mu4du1avXbOkSAiIrq7GRUkLly4gKCgICxatAgdO3Y0dU1ERERkI4wKEpMmTcLx48fRq1cvtGvXztQ1ERERkY0wKkjs378fBQUF6NGjB7p06SJeCSCTyZCenm7SAomIiMh6GRUk/vjjDwBAZWWl3qQ+Xr1hPFdXV16qSUREbY5RQWLlypWmrqPNc3Jy4qWaRETU5hgVJMaPH2/qOu4avFSTiIjaEqOCRExMTIPtMpkMK1asaNWCiIiIyHYYFSS++eYbyGQy8b4Rdf9mkCAiIrq7NXtEQqfTIT8/HxkZGZg4caLJCiMiIiLrJ/nUxhtvvNHgE0HJOvAqESIiMgej72x5u4qKCpw+fRp5eXkmKYpajleJEBGRORgVJPr169dgu7+/f6sWQ62PV4kQEZEp2RnzJkEQ9H7uvfde9O7dG6tWrTJ1fURERGTFjBqROH/+vKnrICIiIhtk1IgEkTXS6XSSlhERUesxOCLh7u5ucGWZTIaioqJWLYjIWJxMSkRkeQaDRN0NqIisFSeTEhFZlsEg8eOPP+q9rqqqQmJiIrZv3w5BEBAYGGjS4oiIiMi6GQwSdUFBq9Viw4YN+Oyzz1BaWgp/f3/ExMQgLCzMLEUSERGRdTIYJK5fv46NGzfik08+QWlpKXx8fLBkyRKEh4ebqz4iIiKyYgaDxBNPPIG//voLMpkM4eHhGD16NOzs7JCRkSG+JyQkxORFEhERkXUyGCQuXbokPulz+/bt2L59u95yXrVBRER0dzMYJBQKBWQymblqoUbwAVtERGStDAaJQ4cOmasOMqCx+yUMGDAAixcvtkxRREREMPIW2WRaTT3yG2j4fgmenp6mLo2IiMggBgkrYOiR3xx1ICIia8YgYUU46kBERLaGD+0iIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpLM6oPEtWvXMHLkSPj4+NS7xwIRERFZltUHiQ4dOmDjxo0YOXKkpUshIiKiO1h9kHBwcMCDDz5o6TKIiIioAWYLEomJiQgNDYWnpyeio6P1lpWXlyMqKgre3t5QqVRISUkxV1lERERWS6fTSVpmTmZ71kanTp0wd+5c7N69G1qtVm/ZokWL4ODggOzsbOTm5mLKlCkICAiAn5+fucojIiKyOo090NHLywvr1q2zUFX6zBYkwsLCAADZ2dkoKSkR26uqqpCWloaff/4ZcrkcKpUKw4YNQ3JyMt5++21zlUdERGSVGnqgozWx+BwJtVoNe3t7KJVKsS0wMBAFBQXi68mTJ2Pv3r1YsGABkpKSLFEmGcEWhuCIiKh1Wfwx4hqNBk5OTnptTk5O0Gg04uuvvvrK3GWRBLYwBEdERK3L4kFCLpejoqJCr62iogJyudxCFVFLWPsQHDWfTqeDvb19s5cR0d3B4kHCy8sLOp0OarUaXl5eAIC8vDxOtCSyEhxpIiJDzBYkampqUFNTg9raWuh0Omi1Wtxzzz1wdHREaGgoEhISkJCQgNzcXGRkZCA1NdVcpRFREzjSRESNMdtky9WrV0OpVGLNmjXYtm0blEolVq9eDQCIi4uDVqtFjx49MGvWLMTHx3NEwoKsadIkJ2kSEVk3s41IxMbGIjY2tsFlzs7OWL9+vblKoSZY01B2Y7UMGDAAixcvNmstRERUn8XnSJB1sqah7IZq8fT0tFA1RER0O4vfR4LaPldXV6s6XUJEZE1s/f+BHJEgk3NycrKq0yVERNbE1k/hMkiQ2VjT6RIiImtiy6dweWqDiIiIJGOQIKM1NdeBbBePKxFJxVMbZDRDcx1s5VweNYzHlYikYpCgZrPlc3nUOB5XIpKCpzaIiIhIMgYJIiIikoxB4i5kTZMmrakWIiJqPs6RuAtZ06RJa6qFiIiaj0HiLmZNk+usqRYiIjIeT20QERGRZAwSREREJBmDBJEF8amoRHeHtvzfM+dIEFkQn4pKdHdoy5PKGSSILIxPRSW6O7TVSeU8tUFERESSMUgQERGRZAwSREREJBmDBBEREUnGIEF0l2jLl5+R9eAlzXcfXrVBdJdoy5efkfXgJc13HwYJortIW738jKwLL2m+u/DUBhEREUnGIEFERESSMUgQERGRZAwSREREJBmDBBEREUnGIEFERESSMUgQERGRZAwSRERkcq6urpLvemmKu2UaWq+mpsastdg63pCKiIhMzsnJSfJdL01xt8ym7vTamnW29bvHMkgQEZHZSL3rpSnulmnoTq+tWWdbv3ssT20QERGRZAwSREREJBmDBBEREUnGIEFERESSMUgQERGRZAwSREREJBmDBBEREUnGIEFERESSMUgQERGRZDYRJD788EOMHTsWr7zyCqqrqy1dDhEREf1/Vh8k8vPzce7cOaSkpGDAgAHYsmWLpUsiIiKi/8/qg8ShQ4cwePBgAMAzzzyDw4cPW7giIiIiqmO2IJGYmIjQ0FB4enoiOjpab1l5eTmioqLg7e0NlUqFlJQUcdnff/+Njh07Arj19LirV6+aq2QiIiJqgtme/tmpUyfMnTsXu3fvhlar1Vu2aNEiODg4IDs7G7m5uZgyZQoCAgLg5+eH+++/H9euXQMAVFRU4IEHHjBXyURERNQEs41IhIWFYcSIEXB2dtZrr6qqQlpaGhYsWAC5XA6VSoVhw4YhOTkZAKBSqbBnzx4AwO7du9GnTx9zlUxERERNsPgcCbVaDXt7eyiVSrEtMDAQBQUFAIBu3bqhS5cuGDt2LPbs2YOJEydaqlSyEa6urtDpdI0uN7TMWrSFPkit39z9NrTNmpoasy6T2j+pfbCWz1FTn3ep6xnqO7Ues53aaIxGo4GTk5Nem5OTEzQajfh64cKF5i6LbJiTkxPs7e3xyiuvQK1W6y3z8vLCunXrLFSZ8dpCH6TWb+5+N7a/AQMGYPHixWZb1pL+SemDNX2ODH3e6/ogdb3mbpOaz+JBQi6Xo6KiQq+toqICcrncQhVRW6FWq5Gfn2/pMlrE1vsgtX5z97uh/Xl6epp1WUuZe3+mYKgPUteTuk0ynsVPbXh5eUGn0+klxry8PPj5+VmwKiIiIjKG2YJETU0NtFotamtrodPpoNVqUVNTA0dHR4SGhiIhIQFVVVU4cuQIMjIyEBERYa7SiIiISCKzBYnVq1dDqVRizZo12LZtG5RKJVavXg0AiIuLg1arRY8ePTBr1izEx8dzRIKIiMgGmG2ORGxsLGJjYxtc5uzsjPXr15urFCIiImolFp8jQURERLaLQYKIiIgkY5AgIiIiyRgkiIiISDIGCSIiIpKMQYKIiIgkkwmCIFi6CFOSyWSWLoGIiMhmNDcWWPxZG6bWxnMSERGRRfHUBhEREUnGIEFERESSMUgQERGRZAwSREREJBmDBBEREUnGIEFERESSMUgQERGRZAwSREREJBmDBBEREUnGINEM5eXliIqKgre3N1QqFVJSUixdUrNERkbCy8sLPj4+8PHxwYABA8RlKSkpUKlU8Pb2xksvvYTy8nJxWVP9NrSuOSQmJiI0NBSenp6Ijo7WW7Zv3z48/fTTUCqViIyMxPnz58Vl169fR0xMDPz8/BAUFITPP/+81dY1pcb6W1xcDIVCIR5fHx8frFq1yuiara2/169fR2xsLFQqFXx9fRESEoJffvmlVeq1pb62teMKAK+99hp69uwJPz8/PPXUU9i0aVOr1GtLfW1Tx1Ugo82cOVN49dVXhcrKSuHQoUOCn5+fkJ+fb+myjBYRESF8/fXX9drz8/MFHx8fITMzU6isrBRmzZolzJgxQ1xuqN9NrWsOO3bsEHbu3Cm8+eabwty5c8X2srIywc/PT/j++++F6upq4f333xdGjhwpLo+LixP+8Y9/COXl5cLJkyeFxx9/XPjll19avK6l+ltUVCS4ubkJN2/ebHA9W+uvRqMREhIShKKiIkGn0wkZGRmCj4+PUFRU1OaOraG+trXjKgi3/r+h1WoFQRCEwsJC4fHHHxeys7Pb3HE11Ne2dFwZJIyk0WgEDw8P4dSpU2LbnDlzhKVLl1qwquZpLEjExcUJs2bNEl+fOXNG8PDwECoqKprst6F1ze3DDz/U+8P61VdfCeHh4eJrjUYjeHl5CYWFhYIgCEJwcLCwe/ducfmyZcvEENSSdc3lzv429T8mW++vIAjCkCFDhB9++KHNH1tB+L++tvXjWlhYKAQFBQmpqalt/rje3te2dFx5asNIarUa9vb2UCqVYltgYCAKCgosWFXzxcfHo3v37hgzZgwOHjwIADh58iQCAgLE93Tt2hUODg5Qq9VN9tvQupZ2Z22Ojo7w8PBAQUEBrl69ij///FNveUBAAE6ePNnidS2tb9++6NWrF+bNm4crV64AQJvo7+XLl6FWq+Hn59fmj+3tfa3T1o7rwoULoVQqMXDgQDz88MMYMmRImz2uDfW1Tls4rgwSRtJoNHByctJrc3JygkajsVBFzff2228jMzMTv//+OyZNmoSpU6fi7NmzDfatY8eOqKysbLLfhta1NI1Gg44dO+q1OTk5if2qe13n9rpbsq6luLi4IC0tDYcOHUJ6ejoqKysxZ84cALD5/t68eRNz5sxBZGQkvL292/SxvbOvbfW4xsfH4+TJk0hJSUFoaCjatWvXZo9rQ31tS8eVQcJIcrkcFRUVem0VFRWQy+UWqqj5goODcd9996F9+/YYP348+vTpg19++QVyubzeh6yiogL33Xdfk/02tK6lNVR7ZWWl2K+613Vur7sl61qKXC7H448/jnvuuQcPPfQQli5dij179ugdL1vsb21tLV5//XW0a9cOS5cubXG9ttjXtnhcAcDe3h4qlQolJSXYuHFjmz2uQMN9bSvHlUHCSF5eXtDpdHpD9nl5eXpDj7ZGJpNBEAT4+voiLy9PbD937hxu3LgBLy+vJvttaF1Lu7O2qqoqnD17Fn5+fnjggQfQqVMnveV5eXnw9fVt8brWQiaTAQAEQbDZ/gqCgNjYWFy+fBlr166Fg4NDi+u1tb7eqS0c1zvpdDqcO3euTR7XO9X19U42fVxNNvuiDZoxY4Ywc+ZMQaPRCIcPH7apqzauXr0q/Prrr0J1dbVw8+ZNITk5WVAqlUJhYaGQn58v+Pr6Cr/99pug0WiEOXPm6E3MMdTvptY1h5s3bwrV1dVCXFycMGfOHLGPpaWlgp+fn/DDDz8I1dXVwgcffKA3s3np0qXCs88+K5SXl4uToOpmNrdkXUv19/fffxcKCwsFnU4nlJWVCa+++qoQERFh0/194403hJEjRwqVlZV67W3x2DbW17Z2XC9fvix89913QmVlpVBTUyP8+uuvglKpFNLT09vccTXU17Z0XBkkmuHKlSvCtGnTBKVSKfTu3VvYtm2bpUsyWmlpqRAaGir4+PgI/v7+wsiRI4U9e/aIy7dt2yb07t1bUCqVwtSpU4UrV66Iy5rqt6F1zSEhIUFwc3PT+0lISBAEQRD27NkjDBgwQPDy8hIiIiKEoqIicT2tVivMmzdP8PX1FXr06CF89tlnetttybqm1Fh/U1JShL59+wpKpVIICgoSXnvtNeHSpUs229/i4mLBzc1N8PT0FLy9vcWf5OTkFtdrS31ta8e1tLRUePbZZwV/f3/B19dXeOaZZ4R///vfrVKvLfW1LR1XmSAIgunGO4iIiKgt4xwJIiIikoxBgoiIiCRjkCAiIiLJGCSIiIhIMgYJIiIikoxBgoiIiCRjkCC6C/Xt2xcKhUJ8cJsl6XQ6REdHo1u3blAoFHj//fdNvs/o6GgoFAqsWLHC5PsiausYJIgspO6PeXBwMKqrqwEAJ06cgEKhgEKhsHB15rNjxw58++23sLe3x0svvYR+/frVe8/BgwehUCjQrVs3C1T4f/vv27evRfZPZM0YJIgs7NKlS9i4caOly2ixmzdvSlqv7jkuzzzzDJYsWYKQkJDWLIuITIxBgsjCZDIZPvnkE3FU4k51IxTFxcUAgBUrVkChUCA6OhoAkJSUBIVCgaFDh2Lx4sXw8fHBoEGDcOLECSxfvhz+/v7o378/9uzZU2/bOTk5GDZsGHx9fREVFYUrV66Iyw4fPozIyEgEBAQgODgYMTEx4vLi4mKxro0bNyI4OBjPPfdcg/UXFxdj+vTp6NmzJwICAhAZGYljx46Jffnoo48AAMnJyVAoFEhKSmrydyYIAuLj49G7d294enoiKCgIzz//vFhfWVkZ5s+fjz59+sDPzw+jRo3Cr7/+2uj29u7di5CQEPj7+8PDwwMqlQoJCQkAbo1GjBs3DgBw/vx5vRGj6upqxMXF4cknn4S3tzeGDx+O9PR0ve2OGDEC3t7e8Pf3x/Dhw5GWltZk/4hsCYMEkYWNGjUKpaWl2LBhQ4u2k5+fj6ysLPj4+KCwsBDjxo3Djh07EBwcjKKiIsTExNRbZ+XKlXjsscfg4uKC9PR0vPnmm+K2JkyYgP/85z8YNGgQvLy8kJSUhFdffRV33lV/2bJlGDx4MHr37l1v+1VVVRg/fjx27NgBLy8vPPnkk8jMzMT48eNx9uxZBAcHo2fPngAAHx8fREVFGfWUwn379mHNmjWwt7fHxIkT0bdvX+Tn50Oj0aC2thZTp07F5s2b4eLigpCQEOTk5GDKlCk4cuRIg9srKSmBi4sLRo8ejYiICGg0GqxatQqpqano3LkzwsLCAAD33XcfoqKiEBUVBQCIjY3Fxx9/DCcnJ4SFheHixYt4+eWXxbkn8+bNQ15eHsLCwhAWFgY7OzsUFBQ02T8iW3KPpQsgutuNHj0aBQUF+PTTT9GrVy/J23F0dMSWLVuQlZWFcePG4dq1a0hNTYWbmxv8/Pzw559/oqysDA8++KC4zoIFC/Dyyy8jNzcXISEhSEtLg0ajwcaNG3Hjxg0EBgbC1dUVrq6uOHbsGA4ePIjTp0+jffv24jY+//xzPPXUUw3W9NNPP6GoqAgeHh7YunUrZDIZoqKikJ6ejs2bN2PhwoU4duwYsrKyEBQUZPREy5qaGgBA165dER4eDl9fXzz44IMQBAHZ2dk4duwY5HI5UlJS4OjoCBcXF3zxxRfYsGED+vTpU29748aNg6urK3JyclBeXg4PDw9cvXoVBw4cwJgxYzBt2jSkpaXhgQceEGssKytDamoq7Ozs0Lt3b9jZ2cHPzw+ZmZn46quv8MQTT6Cmpgbt27cXRzu8vLzqBTEiW8cgQWRhdnZ2mDdvHmbOnInExMQm36/T6Rpsd3d3R4cOHdCxY0exTalUwt7eXnxdVVWlFyR8fHwAAN7e3mJbSUmJeBolKysLWVlZevs5c+YM/P39xdcN/WGuc/78ebEOmUymt68LFy4Y6KVhAwcOxIsvvojk5GTxtENQUBDWr18v1u7m5gZHR0e9fdbVc6e33noLX3/9db32srKyRmuo209tbW2943b27FkAt0ZrlixZgldffRUA4OzsjKVLl2LMmDHGdpXI6vHUBpEVCA8PR7du3bB9+/Z6yzp06AAAqKioAIBGh8bt7Or/53x7iGhIYWEhAODUqVNiW+fOneHu7g4AmD59Oi5cuCD+ZGZmYtiwYXrbuH104k5dunQBcGtCZd038dOnTwNAi65M0el0WLp0Kf744w8cOHAAkZGROH78ODZt2iTWfvHiRXHeSd0+6+q5U93vfdWqVSguLsaUKVMAQKy57vdYW1srrlO3n3bt2iEnJ0f8HZ09exZffPEFAGDw4ME4cOAAcnJysHbtWpSXl2PZsmWS+01kjTgiQWQFZDIZ5s2bh+nTp9db1r17dxw5cgTvvPMOlEolfvzxx1bb70cffYS8vDzxnH5oaCjkcjkmTZqETZs24csvv0RRURFcXFxQWFiIo0ePNvqtviFDhw6Fu7s7zp49i3HjxsHZ2Rk7d+7Evffei4kTJ0qu++jRo5g3bx569eqFBx54AEePHgUA3H///Xj88cfRs2dPZGVlYezYsfD19UVqaipkMhlefPHFBrfn6uqKa9euYf369di7dy927typt9zNzQ3ArdGa+fPnw9PTE7Nnz0Z4eDi2b9+OUaNGYcCAASgvL8fhw4cxefJkxMbGYvjw4XB3d4ebmxsuXrwIAHojRkRtAUckiKxEWFgYAgMD67V/8MEH6NatG3Jzc1FSUoIJEya02j5jYmKQk5ODsrIyhISEYPny5QCAwMBAbNmyBX379sVvv/2G1NRUVFZWYvbs2c3avqOjI5KSkhAWFoZTp05h//796N+/P5KSkuDp6Sm57kceeQSenp7Yv38/Nm3ahOrqakyePBmTJk2CnZ0dNmzYgAkTJqC0tBTp6eno3r07EhMToVKpGtxeQkICvL29UVhYiMrKSrzwwgt6y93d3TFjxgx07NgRmzdvRnJyMoBbV53MmTMHMpkM3377LY4ePYpevXph0KBBAIABAwbg9OnT2Lp1K44cOYL+/fuLV4MQtRUygTN/iIiISCKOSBAREZFkDBJEREQkGYMEERERScYgQURERJIxSBAREZFkDBJEREQkGYMEERERScYgQURERJIxSBAREZFk/w/z8vv0A0cyxwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 540x378 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig , ax = plt.subplots()\n",
    "\n",
    "n , bins , patches = plt.hist(INDEL_DF.num_isolates.astype(float) , bins = 60 , color = 'k' , rwidth = 1.0 , edgecolor='white', linewidth=1.35)\n",
    "\n",
    "plt.title('Number of Isolates Each INDEL is Found in', fontweight = 'bold' , fontsize = 12, color = 'k')\n",
    "plt.xlabel('Number of Isolates', fontweight = 'bold' , fontsize = 12, color = 'k')\n",
    "plt.ylabel('Number of INDELs', fontweight = 'bold' , fontsize = 12, color = 'k')\n",
    "\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "\n",
    "ax.tick_params(labelcolor = 'k')\n",
    "ax.set_yscale(\"log\", nonposy='clip')\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(7.5, 5.25)\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.tick_params(axis='y', which='major', labelsize=12 , labelcolor = 'k')\n",
    "plt.tick_params(axis='x', which='major', labelsize=12 , labelcolor = 'k')\n",
    "\n",
    "#file_name = '/home/rv76/Farhat_Lab/Figures/Final_5/Scraping_rollingDB/Mean_Coverage_Distribution_across_all_isolates.pdf'\n",
    "#plt.savefig(file_name, bbox_inches='tight', dpi = 300 , transparent = True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Re-Structure DataFrame to include the INDEL id, Reference Position, Reference Allele and Alternate Allele(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF = pd.DataFrame(index = range( 0 , np.shape(INDEL_DF)[0] ) ) #index for each INDEL\n",
    "INDEL_annotation_DF['key'] = list(INDEL_DF.indel_id) #key for corresponding INDEL\n",
    "INDEL_annotation_DF['pos'] = list(INDEL_DF.ref_pos) #H37Rv Reference Position\n",
    "INDEL_annotation_DF['ref'] = list(INDEL_DF.ref_allele) #reference allele\n",
    "INDEL_annotation_DF['alt'] = list(INDEL_DF.alt_allele) #alterante allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507, 4)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This is the DataFrame which we'll use to extract INDELs of interest when iterating over all of the VCF files in rollingDB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head() #index = row indices of genotype matrix , 0 supports ref, 1 supports alt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>59502</th>\n",
       "      <td>C_4411205_CA</td>\n",
       "      <td>4411205</td>\n",
       "      <td>C</td>\n",
       "      <td>CA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59503</th>\n",
       "      <td>AG_4411208_A</td>\n",
       "      <td>4411208</td>\n",
       "      <td>AG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59504</th>\n",
       "      <td>C_4411423_CTG</td>\n",
       "      <td>4411423</td>\n",
       "      <td>C</td>\n",
       "      <td>CTG</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59505</th>\n",
       "      <td>CCTG_4411459_C</td>\n",
       "      <td>4411459</td>\n",
       "      <td>CCTG</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59506</th>\n",
       "      <td>AC_4411486_A</td>\n",
       "      <td>4411486</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  key      pos   ref  alt\n",
       "59502    C_4411205_CA  4411205     C   CA\n",
       "59503    AG_4411208_A  4411208    AG    A\n",
       "59504   C_4411423_CTG  4411423     C  CTG\n",
       "59505  CCTG_4411459_C  4411459  CCTG    C\n",
       "59506    AC_4411486_A  4411486    AC    A"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output Annotation for SNVs of interest in Genotypes folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_INDEL_annotation.pkl' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a text file with a single column of all Reference Positions that we're interested in based off of DF above. We will use this text file to subset the *full* VCF files to only the positions that Reference Positions that have an INDEL in at least 1 isolate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get number of unique Reference Positions (with at least 1 INDEL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_ref_pos_with_indels = list( np.sort( list( set( INDEL_annotation_DF.pos ) ) ) )\n",
    "unique_ref_pos_with_indels = [str(ref_pos) for ref_pos in unique_ref_pos_with_indels]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "56184"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(unique_ref_pos_with_indels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['313', '1549', '1552', '1560', '1561']"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unique_ref_pos_with_indels[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['4411205', '4411208', '4411423', '4411459', '4411486']"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unique_ref_pos_with_indels[-5:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_ref_positions = ['NC_000962.3\\t' + ref_pos for ref_pos in unique_ref_pos_with_indels ]\n",
    "\n",
    "#save each element as a new row in a text file\n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/ref_positions_to_select.txt', 'w') as f:\n",
    "    for item in INDEL_ref_positions:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] This section of the notebook prepares the script that iterates over all of the VCF files in rollingDB and outputs a genotype vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vcf\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_ID = 'Peru4988' #INPUT - isolate ID"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load INDEL-Index annotation (for indexing genotypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_INDEL_annotation.pkl' )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert INDEL annotation into dictionaries with keys = ROW indices of genotype array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "#keys: INDEL identifier , values: Genotype Row Index\n",
    "INDEL_dict = dict( zip( list(INDEL_annotation_DF.key) , list(INDEL_annotation_DF.index) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create set of H37Rv Reference Positions that are INDEL sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "H37Rv_ref_positions_of_interest = set(INDEL_annotation_DF.pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construct 'skeleton' Numpy Array that will store the supported base for all INDELs of interest for a single isolate in rollingDB\n",
    "\n",
    "- **0** INDEL not called in isolate (but good quality call otherwise)\n",
    "- **10 - 100** INDEL called in isolate (Reference & Alternate Alleles match), % alternate allele frequency if >= 10%\n",
    "- **-9** for Bad Quality calls\n",
    "\n",
    "#### In constructing the Genotype Vector, we're going to iterate through all Reference Positions of interest and find the supported variant call from the VCF file for each position, assume bad quality Base Call/Missing Data unless evidence from VCF file supports a high-quality variant call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "#number of INDELs we're interested in (where High Quality INDEL calls were made)\n",
    "number_of_INDELs = len(INDEL_dict.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "59507"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "number_of_INDELs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create an empty numpy array to store genotypes\n",
    "genotypes_array = np.full(number_of_INDELs , -9 , dtype = int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-9, -9, -9, ..., -9, -9, -9])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507,)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iterate through the subset VCF file corresponding to the isolate in rollingDB and fill in genotypes vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#genomic data directory\n",
    "VCF_genotypes_dir = '/n/data1/hms/dbmi/farhat/rollingDB/genomic_data/'\n",
    "\n",
    "#directory that stores files for each sequenced isolate\n",
    "directory_for_sequenced_isolate = VCF_genotypes_dir + isolate_ID + '/pilon/'\n",
    "\n",
    "#load in VCF file if one was found\n",
    "vcf_reader = vcf.Reader( open(directory_for_sequenced_isolate  + isolate_ID + '_full.vcf.gz' , 'r') )\n",
    "    \n",
    "#iterate through each Variant Call \n",
    "for record in vcf_reader:\n",
    "    \n",
    "    #check to see if H37Rv Reference Position is one of the INDEL sites of interest\n",
    "    if record.POS in H37Rv_ref_positions_of_interest:\n",
    "\n",
    "        ref_pos = int( record.POS )\n",
    "        ref_allele = str( record.REF )  #allele on H37Rv\n",
    "        alt_allele = str( record.ALT[0] ) #alternate allele supported by reads\n",
    "        quality_metrics = record.INFO.keys()\n",
    "        variant_i_key = ref_allele + '_' + str(ref_pos) + '_' + alt_allele #variant key, unique identifier for a specific variant\n",
    "\n",
    "        #get the all relevant INDEL keys for possible INDELs @ this reference position\n",
    "        INDEL_keys_at_ref_pos = list( INDEL_annotation_DF[INDEL_annotation_DF.pos == ref_pos].key )\n",
    "\n",
    "         #BOOLEAN for whether INDEL was found\n",
    "        INDEL_FOUND = False\n",
    "\n",
    "        ###################################################################################################\n",
    "        #The variant call at this reference position was high quality and was an INDEL of interest\n",
    "        ###################################################################################################\n",
    "        #check to see if the call is a PASS by Pilon (PASS by Pilon ; there's exactly 1 alternate allele ; call supports the alternate allele)\n",
    "        if (record.FILTER == []) and (len(record.ALT) == 1) and (record.ALT != [None]):\n",
    "\n",
    "            #check to see if variant is INDEL (either Refeference or Alternate Alleles are > 1 base but not both ; 1 allele is 1bp, the other allele is > 1bp)\n",
    "            if ( (len(ref_allele) > 1) or (len(alt_allele) > 1) ) and ( (len(ref_allele) == 1) or (len(alt_allele) == 1) ):\n",
    "\n",
    "                #keep only INDELs <= 10bp long (1 allele is 1bp ; the other allele is <= 11bp)\n",
    "                if (len(ref_allele) <= 11) and (len(alt_allele) <= 11):\n",
    "\n",
    "                    #analyze only INDELs with Depth metrics (larger structural variants don't have this)\n",
    "                    if ('MQ' in quality_metrics) and ('DC' in quality_metrics) and ('IC' in quality_metrics) and ('TD' in quality_metrics) and ('DP' in quality_metrics):\n",
    "\n",
    "                        ##### Retrieve Relevant information for filtering quality of Base Call #####\n",
    "                        # Mean Mapping Quality @ locus\n",
    "                        MQ = record.INFO['MQ']\n",
    "                        # Number of Reads w/ Deletion \n",
    "                        DC = record.INFO['DC']\n",
    "                        # Number of Reads w/ Insertion\n",
    "                        IC = record.INFO['IC']\n",
    "                        # Depth of Valid Reads in Pileup\n",
    "                        VD = record.INFO['DP']\n",
    "                        # Total Depth in Pileup\n",
    "                        TD = record.INFO['TD']\n",
    "\n",
    "                        ### Filtering Criteria\n",
    "\n",
    "                        #---> Mean Mapping Quality > 30\n",
    "                        #---> Number of High Quality Reads >= 20\n",
    "                        if (MQ > 30) and (VD >= 20): #INDEL passed filtering criteria!\n",
    "                            \n",
    "                            # Calculate Alternate Allele Frequency\n",
    "                            ########################################################################\n",
    "\n",
    "                            #calculate INDEL (Alternate Allele) Frequency from Deletion/Insertion Count & Total Depth (all metrics include ALL reads)\n",
    "                            alt_allele_frequency = float(max(DC , IC)) / float(TD)\n",
    "\n",
    "                            ########################################################################\n",
    "\n",
    "                            # use INDEL calls that pass the following filters\n",
    "                            #---> alternate (INDEL) allele frequency >= 10%\n",
    "                            #---> Reads supporting INDEL variant >= 10 reads\n",
    "                            if (alt_allele_frequency >= 0.10) and (float(max(DC , IC)) >= 10):\n",
    "\n",
    "                                #INDEL passed filters, find relevant INDEL index and store '1' and store '0' in for all other possible INDELs at this reference position\n",
    "                                if variant_i_key in INDEL_dict.keys():\n",
    "\n",
    "                                    #BOOLEAN for whether INDEL was found\n",
    "                                    INDEL_FOUND = True\n",
    "\n",
    "                                    #store a '10' - '100' in Genotype Vector (INDEL index) referencing that the INDEL is supported in this isolate and alt AF >= 10%\n",
    "                                    genotypes_array[INDEL_dict[variant_i_key]] = int(alt_allele_frequency * 100)\n",
    "\n",
    "                                    #store a '0' in Genotype Vector (INDEL index) for all other INDELs @ this reference position (if there were multiple indels here)\n",
    "                                    INDELs_not_detected = INDEL_keys_at_ref_pos.remove(variant_i_key)\n",
    "\n",
    "                                    #if list of not empty, then there are some INDEL @ this ref pos that were not detected in this isolate\n",
    "                                    if INDELs_not_detected != None: \n",
    "\n",
    "                                        for INDEL_i_key in INDELs_not_detected:\n",
    "\n",
    "                                            genotypes_array[INDEL_dict[INDEL_i_key]] = 0\n",
    "\n",
    "        ###################################################################################################\n",
    "        #The variant call at this reference position was high quality but was not an INDEL of interest\n",
    "        ###################################################################################################\n",
    "        if (INDEL_FOUND == False) and (record.FILTER == []) and ('MQ' in quality_metrics) and ('DP' in quality_metrics):\n",
    "\n",
    "            ##### Retrieve Relevant information for filtering quality of Base Call #####\n",
    "            # Mean Mapping Quality @ locus\n",
    "            MQ = record.INFO['MQ']\n",
    "            # Depth of Valid Reads in Pileup\n",
    "            VD = record.INFO['DP']\n",
    "\n",
    "            ### Filtering Criteria\n",
    "\n",
    "            #---> Mean Mapping Quality > 30\n",
    "            #---> Number of High Quality Reads >= 20\n",
    "            if (MQ > 30) and (VD >= 20): #variant passed filtering criteria!\n",
    "\n",
    "                #store a '0' in Genotype Vector (INDEL index) for all INDELs @ this reference position\n",
    "                for INDEL_i_key in INDEL_keys_at_ref_pos:\n",
    "\n",
    "                    genotypes_array[INDEL_dict[INDEL_i_key]] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pickle Genotypes Array for downstream analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, ..., 0, 0, 0])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/Genotypes_per_isolate/' + isolate_ID + '/' + isolate_ID + '_genotypes' , genotypes_array , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] This section of the notebook launches the script that was created in the previous section of the notebook and submits jobs for batches of isolates to make genotypes arrays for each isolate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "from slurmpy import Slurm\n",
    "import numpy as np\n",
    "import itertools\n",
    "import shutil"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every isolate in the pool of **35,694** isolates in rollingDB had a reduced VCF file and a corresponding Full VCF file. We will iterate through each of these Full VCF files and extract the base calls for only the **56,184** reference positions that we're interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_isolate_annotation.pkl' )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_id_series.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35694"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_id_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create directory for each isolate to store output for genotyping (only run this once)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "for tag in isolate_id_series.values:\n",
    "\n",
    "    #create directory to store subset VCF and genotypes for each isolate\n",
    "    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/Genotypes_per_isolate/' + tag\n",
    "    if os.path.exists(isolate_genotype_dir):\n",
    "        shutil.rmtree(isolate_genotype_dir)\n",
    "        os.makedirs(isolate_genotype_dir)\n",
    "    elif not os.path.exists(isolate_genotype_dir):\n",
    "        os.makedirs(isolate_genotype_dir)\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Check Progress of Jobs:__ Isolates that need to have Full VCF >> Genotypes Array Reference Positions of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_ids_to_run_VCF_extraction = []\n",
    "isolate_ids_completed = []\n",
    "\n",
    "for tag in isolate_id_series.values:\n",
    "\n",
    "    #create directory to store subset VCF and genotypes for each isolate\n",
    "    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/Genotypes_per_isolate/' + tag\n",
    "    \n",
    "    #genotypes array file\n",
    "    genotypes_array_path = isolate_genotype_dir + '/' + tag + '_genotypes.npy'\n",
    "    \n",
    "    #check for the existence of a genotypes file, if TRUE job has either started or completed\n",
    "    if os.path.exists(genotypes_array_path):\n",
    "            \n",
    "        isolate_ids_completed.append(tag)\n",
    "            \n",
    "    #genotypes file doesn't exist, task never started\n",
    "    else:\n",
    "        \n",
    "        isolate_ids_to_run_VCF_extraction.append(tag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35370"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_ids_completed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "324"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_ids_to_run_VCF_extraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construct genotypes array for reference positions of interest for isolates that had Full VCF files "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series = pd.Series(isolate_ids_to_run_VCF_extraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract genotypes for each isolate from VCF file containing a subset of Reference Positions from the Full VCF files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create BASH scripts to extract genotypes of N isolates per job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_isolates_per_job = 120\n",
    "\n",
    "njobs = int( np.ceil( float( len( isolate_id_series ) ) / float(N_isolates_per_job) ) ) #number of jobs required if we split for every N isolates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33\n"
     ]
    }
   ],
   "source": [
    "print njobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "bash_scripts = []\n",
    "\n",
    "#split jobs up into tasks of N\n",
    "for first_isolate_i_for_job in np.arange(0 , (njobs)*N_isolates_per_job , N_isolates_per_job):\n",
    "\n",
    "    last_isolate_i_for_job = first_isolate_i_for_job + N_isolates_per_job\n",
    "    \n",
    "    commands_list = []\n",
    "\n",
    "    for isolate_i in range(first_isolate_i_for_job , last_isolate_i_for_job): \n",
    "\n",
    "        try:\n",
    "            tag = isolate_id_series[isolate_i]\n",
    "\n",
    "            commands_list.append( 'python /home/rv76/Farhat_Lab/Python_Scripts/homoplasy_project/mmpR_BDQ_mutant_project/rollingDB_extract_HQ_mixed_INDEL_genotypes.py {}'.format(tag) )\n",
    "            \n",
    "        except KeyError: #applies only to the last job if it has < N isolates\n",
    "            continue\n",
    "        \n",
    "    bash_scripts.append(commands_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Submit each job to O2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704688\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno1 : 18704688\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704689\n",
      "submitted: Submitted batch job 18704690\n",
      "submitted: Submitted batch job 18704691\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno2 : 18704689\n",
      "Geno3 : 18704690\n",
      "Geno4 : 18704691\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704692\n",
      "submitted: Submitted batch job 18704693\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno5 : 18704692\n",
      "Geno6 : 18704693\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704694\n",
      "submitted: Submitted batch job 18704695\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno7 : 18704694\n",
      "Geno8 : 18704695\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704696\n",
      "submitted: Submitted batch job 18704697\n",
      "submitted: Submitted batch job 18704698\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno9 : 18704696\n",
      "Geno10 : 18704697\n",
      "Geno11 : 18704698\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704699\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno12 : 18704699\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704700\n",
      "submitted: Submitted batch job 18704701\n",
      "submitted: Submitted batch job 18704702\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno13 : 18704700\n",
      "Geno14 : 18704701\n",
      "Geno15 : 18704702\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704703\n",
      "submitted: Submitted batch job 18704704\n",
      "submitted: Submitted batch job 18704705\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno16 : 18704703\n",
      "Geno17 : 18704704\n",
      "Geno18 : 18704705\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704706\n",
      "submitted: Submitted batch job 18704707\n",
      "submitted: Submitted batch job 18704708\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno19 : 18704706\n",
      "Geno20 : 18704707\n",
      "Geno21 : 18704708\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704709\n",
      "submitted: Submitted batch job 18704710\n",
      "submitted: Submitted batch job 18704711\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno22 : 18704709\n",
      "Geno23 : 18704710\n",
      "Geno24 : 18704711\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704712\n",
      "submitted: Submitted batch job 18704713\n",
      "submitted: Submitted batch job 18704714\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno25 : 18704712\n",
      "Geno26 : 18704713\n",
      "Geno27 : 18704714\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704715\n",
      "submitted: Submitted batch job 18704716\n",
      "submitted: Submitted batch job 18704717\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno28 : 18704715\n",
      "Geno29 : 18704716\n",
      "Geno30 : 18704717\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704718\n",
      "submitted: Submitted batch job 18704719\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geno31 : 18704718\n",
      "Geno32 : 18704719\n",
      "Geno33 : 18704720\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 18704720\n"
     ]
    }
   ],
   "source": [
    "job_num = 1\n",
    "for job_i_commands_list in bash_scripts:\n",
    "    \n",
    "    #append all commands in a single string to be submitted as a job\n",
    "    genotypes_extract_job = ''\n",
    "    for command_i in job_i_commands_list:\n",
    "        genotypes_extract_job  = genotypes_extract_job  + '\\n' + command_i\n",
    "    \n",
    "    #directory where you want output + error files\n",
    "    os.chdir('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/O2_SLURM_logs/make_genotypes_arrays_from_VCF/')\n",
    "\n",
    "    job_name = 'Geno' + str(job_num)\n",
    "\n",
    "    s = Slurm(job_name , {'partition':'short' , 'N':'1' , 't':'0-12:00:00' , 'mem':'2G' , 'mail-type':'FAIL' , 'mail-user':'roger_vargas@g.harvard.edu'})\n",
    "\n",
    "    #submits the job\n",
    "    job_id = s.run(genotypes_extract_job)\n",
    "\n",
    "    print job_name  + ' : ' +  str(job_id)\n",
    "    job_num += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [6] Inspect the distribution of missing values after Genotype Matrix is filled in and filter out INDELs (Reference Positions) where values are missing in many isolates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.ticker as ticker\n",
    "from itertools import compress\n",
    "from matplotlib.colors import LogNorm\n",
    "from matplotlib import gridspec\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make sure all isolates have a corresponding _genotypes_ array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_isolate_annotation.pkl' )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35694,)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_id_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "__Check for Genotypes Arrays:__ Get isolate IDs that have a corresponding genotypes array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_ids_with_genotypes = []\n",
    "isolate_ids_without_genotypes = []\n",
    "\n",
    "for tag in isolate_id_series.values:\n",
    "\n",
    "    #create directory to store subset VCF and genotypes for each isolate\n",
    "    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/Genotypes_per_isolate/' + tag\n",
    "    \n",
    "    #genotypes_array\n",
    "    genotypes_array_pickle_file = isolate_genotype_dir + '/' + tag + '_genotypes.npy'\n",
    "    \n",
    "    #check for the existence of a genotypes file for this isolate\n",
    "    if os.path.exists(genotypes_array_pickle_file):\n",
    "            \n",
    "        isolate_ids_with_genotypes.append(tag)\n",
    "            \n",
    "    else:\n",
    "        \n",
    "        isolate_ids_without_genotypes.append(tag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35694"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_ids_with_genotypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(isolate_ids_without_genotypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "del isolate_ids_with_genotypes #clear variable from memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concatenate genotypes arrays into a genotypes matrix for all isolates in rollingDB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Store the __genotypes array__ for each isolate into a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.0280159130386\n",
      "0.0560318260772\n",
      "0.0840477391158\n",
      "0.112063652154\n",
      "0.140079565193\n",
      "0.168095478232\n",
      "0.19611139127\n",
      "0.224127304309\n",
      "0.252143217347\n",
      "0.280159130386\n",
      "0.308175043425\n",
      "0.336190956463\n",
      "0.364206869502\n",
      "0.39222278254\n",
      "0.420238695579\n",
      "0.448254608618\n",
      "0.476270521656\n",
      "0.504286434695\n",
      "0.532302347734\n",
      "0.560318260772\n",
      "0.588334173811\n",
      "0.616350086849\n",
      "0.644365999888\n",
      "0.672381912927\n",
      "0.700397825965\n",
      "0.728413739004\n",
      "0.756429652042\n",
      "0.784445565081\n",
      "0.81246147812\n",
      "0.840477391158\n",
      "0.868493304197\n",
      "0.896509217235\n",
      "0.924525130274\n",
      "0.952541043313\n",
      "0.980556956351\n"
     ]
    }
   ],
   "source": [
    "genotypes_directory = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/Genotypes_per_isolate/'\n",
    "genotype_arrays_per_isolate_list = []\n",
    "\n",
    "#iterate through the genotypes arrays for each isolate\n",
    "isolate_i = 0\n",
    "for isolate_ID in isolate_id_series:\n",
    "    \n",
    "    #load genotypes array\n",
    "    genotypes_array = np.load(genotypes_directory + isolate_ID + '/' + isolate_ID + '_genotypes.npy')\n",
    "    \n",
    "    #convert genotypes array from dtype 'int64' to 'int8' to save lots of memory\n",
    "    genotypes_array = genotypes_array.astype('int8')\n",
    "    \n",
    "    #store genotypes array for this isolate in a list with all other genotypes arrays\n",
    "    genotype_arrays_per_isolate_list.append(genotypes_array)\n",
    "    \n",
    "    if isolate_i % 1000 == 0:\n",
    "        print float(isolate_i) / float(len(isolate_id_series))\n",
    "    isolate_i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stack genotype arrays for all isolates to create one __genotypes matrix__ for all isolates "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "genotypes_matrix = np.stack(genotype_arrays_per_isolate_list , axis = -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507, 35694)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       ...,\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0]], dtype=int8)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2124042858 bytes\n"
     ]
    }
   ],
   "source": [
    "print(\"%d bytes\" % (genotypes_matrix.size * genotypes_matrix.itemsize))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genotypes Matrix is ~2.12 GB large"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __Genotypes Matrix__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load __INDEL annotation file__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59507, 4)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "del genotype_arrays_per_isolate_list #clear variable from memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [6.1] Drop INDEL sites that *[1]* have low Empirical Base Pair Recall (EBPR) Scores, *[2]* are located in Mobile Genetic Element (MGE) regions or *[3]* are located in a region with overlapping genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Filter 1)  *FILTER* out INDELs with low EBR scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load list of low EBR score Reference positions to Drop (IL mapping/variant calling: sites with empirical score < 0.90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes/H37Rv_sites_with_low_EBR_to_drop.pkl', 'rb') as f:\n",
    "    H37Rv_positions_to_drop = pickle.load(f)\n",
    "    \n",
    "#convert to set, faster to check\n",
    "H37Rv_positions_to_drop = set(H37Rv_positions_to_drop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter out sites with low EBR scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop unique INDELs from the __GENOTYPE MATRIX__ and __INDEL ANNOTATION FILE__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "#construct boolean filter\n",
    "INDELs_to_keep_filter = []\n",
    "for H37Rv_pos_i in INDEL_annotation_DF.pos:\n",
    "    \n",
    "    if H37Rv_pos_i not in H37Rv_positions_to_drop:\n",
    "        \n",
    "        INDELs_to_keep_filter.append(True)\n",
    "        \n",
    "    elif H37Rv_pos_i in H37Rv_positions_to_drop:\n",
    "        \n",
    "        INDELs_to_keep_filter.append(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter Genotypes Matrix\n",
    "genotypes_matrix = genotypes_matrix[INDELs_to_keep_filter , :]\n",
    "\n",
    "#filter SNP annotation file\n",
    "INDEL_annotation_DF = INDEL_annotation_DF[INDELs_to_keep_filter]\n",
    "INDEL_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered INDEL annotation DF (so new index matches indexing of genotypes matrix rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3246"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(INDELs_to_keep_filter) - sum(INDELs_to_keep_filter) #number of sites dropped by this filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56261, 35694)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56261, 4)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF) #Annotation for Rows of Genotypes Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __Genotypes Matrix__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __INDEL annotation file__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Filter 2) *FILTER* out INDELs in MGE regions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Collect list of MGE Reference Positions to Drop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "#csv files with category for each gene\n",
    "gene_categories_df = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>Gene_Category</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Rv0001</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Rv0002</td>\n",
       "      <td>Rv0002</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Rv0003</td>\n",
       "      <td>Rv0003</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Rv0004</td>\n",
       "      <td>Rv0004</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Rv0005</td>\n",
       "      <td>Rv0005</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     name gene_id          Gene_Category\n",
       "0  Rv0001  Rv0001              Essential\n",
       "1  Rv0002  Rv0002          Non-Essential\n",
       "2  Rv0003  Rv0003          Non-Essential\n",
       "3  Rv0004  Rv0004          Non-Essential\n",
       "4  Rv0005  Rv0005  Antibiotic Resistance"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_categories_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reference genome annotation\n",
    "reference_genome_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt', '\\t').set_index('name')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collect positions for **MGEs**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "H37Rv_positions_MGE_genes = []\n",
    "\n",
    "for gene_id in gene_categories_df[gene_categories_df.Gene_Category == 'Mobile Genetic Element'].gene_id:\n",
    "    \n",
    "    #get reference positions tht corresond to each gene\n",
    "    gene_id_info = reference_genome_annotation.loc[gene_id, :]\n",
    "    gene_id = gene_id_info.name\n",
    "\n",
    "    chrom_start = gene_id_info.chromStart\n",
    "    chrom_end = gene_id_info.chromEnd\n",
    "\n",
    "    #find the position of the first base relative to H37Rv in 5' -> 3'\n",
    "    H37Rv_start = min(chrom_start , chrom_end)\n",
    "    H37Rv_end = max(chrom_start , chrom_end)\n",
    "    \n",
    "    #store all corresponding H37Rv Reference Positions to gene_id in list of reference positions to scan\n",
    "    H37Rv_positions_MGE_genes.append([H37Rv_start+1 , H37Rv_end+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "108"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(H37Rv_positions_MGE_genes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FILTER out sites in MGE regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "#construct boolean filter\n",
    "INDELs_to_keep_filter = []\n",
    "for H37Rv_pos_i in INDEL_annotation_DF.pos:\n",
    "\n",
    "    bad_pos_i = False\n",
    "\n",
    "    for problem_region in H37Rv_positions_MGE_genes:\n",
    "\n",
    "        if (H37Rv_pos_i > problem_region[0]) and (H37Rv_pos_i < problem_region[1]):\n",
    "            bad_pos_i = True\n",
    "            break\n",
    "\n",
    "    if bad_pos_i == False:\n",
    "        INDELs_to_keep_filter.append(True)\n",
    "    else:\n",
    "        INDELs_to_keep_filter.append(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter Genotypes Matrix\n",
    "genotypes_matrix = genotypes_matrix[INDELs_to_keep_filter , :]\n",
    "\n",
    "#filter INDEL annotation file\n",
    "INDEL_annotation_DF = INDEL_annotation_DF[INDELs_to_keep_filter]\n",
    "INDEL_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered INDEL annotation DF (so new index matches indexing of genotypes matrix rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "762"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(INDELs_to_keep_filter) - sum(INDELs_to_keep_filter) #number of sites dropped by this filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55499, 35694)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55499, 4)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF) #Annotation for Rows of Genotypes Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __Genotypes Matrix__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __INDEL annotation file__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Filter 3) *FILTER* out INDELs that located in H37Rv Reference Positions with overlapping genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load dictionary that contains the mapping between each H37Rv Reference Position and H37Rv locus tag, from the file: */n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/pickled_files/dicts_for_SNP_annotation/H37Rv_coord_gene_mapping.pickle', 'rb') as handle:\n",
    "    ReferencePosition_Gene_mapping = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Rv0001']"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ReferencePosition_Gene_mapping[1] #one gene mapped to position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Rv2035', 'Rv2036']"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ReferencePosition_Gene_mapping[2282101] #two genes mapped to position (overlapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ReferencePosition_Gene_mapping[2281250] #no genes mapped to position (intergenic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FILTER out sites that map to 2 (or more) H37Rv locus tags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDELs_to_keep_filter = np.array( [(len(ReferencePosition_Gene_mapping[ref_pos_i]) < 2) for ref_pos_i in INDEL_annotation_DF.pos] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "55268"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(INDELs_to_keep_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter Genotypes Matrix\n",
    "genotypes_matrix = genotypes_matrix[INDELs_to_keep_filter , :]\n",
    "\n",
    "#filter INDEL annotation file\n",
    "INDEL_annotation_DF = INDEL_annotation_DF[INDELs_to_keep_filter]\n",
    "INDEL_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered INDEL annotation DF (so new index matches indexing of genotypes matrix rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "231"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(INDELs_to_keep_filter) - sum(INDELs_to_keep_filter) #number of sites dropped by this filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 35694)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 4)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF) #Annotation for Rows of Genotypes Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __Genotypes Matrix__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save __INDEL annotation file__ in case kernel crashes when filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [6.2] Drop low-quality isolates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Filter 4) *FILTER* to retain the same isolates that form the *SNPs* Genotypes Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep things consistent with the SNP genotypes matrix, we will drop the isolates that were dropped as part of constructing the (**SNPs**) Genotypes Matrix in which SNPs were dropped that were poorly called in > 10% of isolates and Isolates were dropped with > 10% missing calls. These are also the isolates that were used in constructing the phylogenies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_annotation_df_from_SNP_analysis = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lineage_1</th>\n",
       "      <th>lineage_2</th>\n",
       "      <th>lineage_3</th>\n",
       "      <th>lineage_4</th>\n",
       "      <th>lineage_5</th>\n",
       "      <th>lineage_6</th>\n",
       "      <th>lineage_7</th>\n",
       "      <th>lineage_8</th>\n",
       "      <th>lineage_9</th>\n",
       "      <th>lineage_10</th>\n",
       "      <th>lineage_11</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3558733</td>\n",
       "      <td>4.2.1.2.1.1.i3.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03648641</td>\n",
       "      <td>4.2.1.2.2.1.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03647419</td>\n",
       "      <td>3.1.1.i1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3671418</td>\n",
       "      <td>4.2.1.2.1.1.i1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN07659096</td>\n",
       "      <td>1.1.1.2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  lineage_1 lineage_2 lineage_3 lineage_4 lineage_5 lineage_6 lineage_7  \\\n",
       "0         4         2         1         2         1         1        i3   \n",
       "1         4         2         1         2         2         1         1   \n",
       "2         3         1         1        i1       NaN       NaN       NaN   \n",
       "3         4         2         1         2         1         1        i1   \n",
       "4         1         1         1         2       NaN       NaN       NaN   \n",
       "\n",
       "  lineage_8 lineage_9 lineage_10 lineage_11    isolate_ID      lineage_call  \\\n",
       "0         1       NaN        NaN        NaN  SAMEA3558733  4.2.1.2.1.1.i3.1   \n",
       "1       NaN       NaN        NaN        NaN  SAMN03648641     4.2.1.2.2.1.1   \n",
       "2       NaN       NaN        NaN        NaN  SAMN03647419          3.1.1.i1   \n",
       "3       NaN       NaN        NaN        NaN  SAMEA3671418    4.2.1.2.1.1.i1   \n",
       "4       NaN       NaN        NaN        NaN  SAMN07659096           1.1.1.2   \n",
       "\n",
       "  group  \n",
       "0    4B  \n",
       "1    4B  \n",
       "2     3  \n",
       "3    4B  \n",
       "4     1  "
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_df_from_SNP_analysis.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 14)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_df_from_SNP_analysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 35694)"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop isolates from the __GENOTYPE MATRIX__ that were dropped in SNP genotypes analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35694"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len( np.array( [isolate_id in list(isolate_annotation_df_from_SNP_analysis.isolate_ID) for isolate_id in list(isolate_id_series)] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolates_to_keep_filter = np.array( [isolate_id in list(isolate_annotation_df_from_SNP_analysis.isolate_ID) for isolate_id in list(isolate_id_series)] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31428"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(isolates_to_keep_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "genotypes_matrix = genotypes_matrix[: , isolates_to_keep_filter]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 31428)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_matrix) #Genotype Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop isolates with a high abundance of missing values from __ISOLATE ANNOTATION FILE__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_annotation_series = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes/genotypes_isolate_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_series.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(35694,)"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_series) #Annotation for Columns of Genotypes Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_annotation_series = isolate_annotation_series[isolates_to_keep_filter]\n",
    "\n",
    "#re-index new filtered isolate annotation series (so new index matches indexing of genotypes matrix columns)\n",
    "isolate_annotation_series.reset_index(drop = True , inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_series.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428,)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_series) #Annotation for Columns of Genotypes Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save annotation files and genotypes matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_series.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_isolate_annotation.pkl')\n",
    "\n",
    "#save SNP annotation file (rows of Genotype Matrix)\n",
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation.pkl')\n",
    "\n",
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [7] Annotate INDELs w/ H37Rv Locus Info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Retrieve information for H37Rv Reference Genome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rv76/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:18: DeprecationWarning: \n",
      ".ix is deprecated. Please use\n",
      ".loc for label based indexing or\n",
      ".iloc for positional indexing\n",
      "\n",
      "See the documentation here:\n",
      "http://pandas.pydata.org/pandas-docs/stable/indexing.html#ix-indexer-is-deprecated\n"
     ]
    }
   ],
   "source": [
    "####################################################################################################################################################################################\n",
    "#load reference genome and reference annotation\n",
    "reference_genome = '/home/rv76/Farhat_Lab/Reference_Seqs/H37Rv/h37rv.fasta'\n",
    "for reference_genome in SeqIO.parse(reference_genome, \"fasta\"):\n",
    "    reference_genome.seq.alphabet = IUPAC.unambiguous_dna\n",
    "\n",
    "reference_genome_annotation = pd.read_csv('/home/rv76/Farhat_Lab/CSV_files/h37rv_genome_summary.txt', '\\t').set_index('name')\n",
    "\n",
    "####### Make a dictionary of all H37Rv reference positions and corresponding genes #######\n",
    "ReferencePosition_Gene_mapping = {} #keys: H37Rv Reference Positions , values: gene_ids (may be multiple genes for some Reference Positions)\n",
    "\n",
    "#store a list corresponding to every Reference Position (to store all genes that map to the Reference Position)\n",
    "for H37Rv_RefPos in range(0 , len(reference_genome.seq) + 1):\n",
    "    ReferencePosition_Gene_mapping[H37Rv_RefPos] = []\n",
    "\n",
    "for gene_id_index in range(0 , len(reference_genome_annotation.index)):\n",
    "\n",
    "    gene_id_info = reference_genome_annotation.ix[gene_id_index , :]\n",
    "    gene_id = gene_id_info.name\n",
    "\n",
    "    chrom_start = gene_id_info.chromStart\n",
    "    chrom_end = gene_id_info.chromEnd\n",
    "\n",
    "    #find the position of the first base relative to H37Rv in 5' -> 3'\n",
    "    H37Rv_start = min(chrom_start , chrom_end)\n",
    "    H37Rv_end = max(chrom_start , chrom_end)\n",
    "\n",
    "    #store all corresponding H37Rv Reference Positions to gene_id in dictionary\n",
    "    for H37Rv_RefPos in range(H37Rv_start+1 , H37Rv_end+1):\n",
    "\n",
    "        ReferencePosition_Gene_mapping[H37Rv_RefPos].append(gene_id) #append gene_id to list already in dict\n",
    "        \n",
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Annotate INDELs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To annotate INDELs, we get gene_ids that the INDELs 'land in' and the gene-relative coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A\n",
       "1         TC_1549_T  1549          TC   T\n",
       "2        TAA_1552_T  1552         TAA   T\n",
       "3         TA_1552_T  1552          TA   T\n",
       "4         T_1552_TA  1552           T  TA"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#load INDEL annotation file (rows of Genotype Matrix)\n",
    "INDEL_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation.pkl')\n",
    "\n",
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 4)"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gene Categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get Gene Categories\n",
    "gene_categories = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories.csv').set_index('name')\n",
    "\n",
    "#create dictionary, keys: H37Rv locus, values: gene category\n",
    "gene_categories_map = dict( zip( list(gene_categories.gene_id) , list(gene_categories.Gene_Category) ) ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gene Names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse through Mycobrowser (tuberculist) to collect gene names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "mycobrowser_table = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/Mycobrowser_Release3/Mycobacterium_tuberculosis_H37Rv_txt_v3.txt' , sep = '\\t')\n",
    "\n",
    "#retreive columns for locus tag & gene name\n",
    "mycobrowser_table = mycobrowser_table.loc[: , ['Locus' , 'Name']]\n",
    "\n",
    "#create dictionary, keys: H37Rv loucs, values: gene category\n",
    "gene_names_map = dict( zip( list(mycobrowser_table.Locus) , list(mycobrowser_table.Name) ) ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Iterate through each INDEL to annotate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some INDELs are present in more than 1 gene, have to account for this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.200007237461\n",
      "0.400014474922\n",
      "0.600021712383\n",
      "0.800028949844\n"
     ]
    }
   ],
   "source": [
    "#create dictionaries for each columsn of information\n",
    "locus_tag_dict = {} #H37Rv id\n",
    "genomic_coord_dict = {} #coordinate on gene, 5' -> 3'\n",
    "gene_categ_dict = {} #gene Category\n",
    "gene_name_dict = {} #gene Name\n",
    "\n",
    "for INDEL_index in INDEL_annotation_DF.index:\n",
    "\n",
    "    #reference position, ref allele and alt allele for INDEL\n",
    "    INDEL_info = INDEL_annotation_DF.loc[INDEL_index , :]\n",
    "\n",
    "    #find the gene that INDEL occurs on\n",
    "    ref_seq_position = INDEL_info.pos\n",
    "\n",
    "    #check list corresponding to H37Rv coordinate to see if there are any genes associated with RefPosition\n",
    "    if len(ReferencePosition_Gene_mapping[ref_seq_position]) == 1:\n",
    "        \n",
    "        #get gene that ReferencePosition is mapped to (i.e. INDEL filtering ensures that INDEL corresponds to only 1 gene)\n",
    "        gene_intergenic_id = ReferencePosition_Gene_mapping[ref_seq_position][0]\n",
    "\n",
    "        #find genomic coordinate of INDEL relative to gene (subtract 1 since reference seq starts counting at 1)\n",
    "        gene_relative_coord = (ref_seq_position - 1) - min( reference_genome_annotation.loc[gene_intergenic_id , 'chromStart'] , reference_genome_annotation.loc[gene_intergenic_id , 'chromEnd'] )\n",
    "\n",
    "        #store all information for INDEL in relevant dictionaries\n",
    "        locus_tag_dict[INDEL_index] = gene_intergenic_id\n",
    "        genomic_coord_dict[INDEL_index] = gene_relative_coord\n",
    "\n",
    "        #use gene id to get category and names from mappings\n",
    "        if gene_intergenic_id in gene_categories_map.keys():\n",
    "            gene_categ_dict[INDEL_index] = gene_categories_map[gene_intergenic_id]\n",
    "        else:\n",
    "            gene_categ_dict[INDEL_index] = 'None'\n",
    "\n",
    "        if gene_intergenic_id in gene_names_map.keys():\n",
    "            gene_name_dict[INDEL_index] = gene_names_map[gene_intergenic_id] \n",
    "        else:\n",
    "            gene_name_dict[INDEL_index] = 'None'\n",
    "            \n",
    "    #H37Rv RefSeq position was not found in dictionary so Reference position must be intergenic    \n",
    "    elif len(ReferencePosition_Gene_mapping[ref_seq_position]) == 0:\n",
    "\n",
    "        #store all information for INDEL in relevant dictionaries\n",
    "        locus_tag_dict[INDEL_index] = 'intergenic'\n",
    "        genomic_coord_dict[INDEL_index] = np.nan\n",
    "        gene_categ_dict[INDEL_index] = 'None'\n",
    "        gene_name_dict[INDEL_index] = 'None'\n",
    "\n",
    "    #keep track of progress of loop\n",
    "    num_INDELs = float( np.shape(INDEL_annotation_DF)[0] )\n",
    "    if INDEL_index % np.ceil(0.20 * num_INDELs) == 0:\n",
    "        print float(INDEL_index) / num_INDELs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Update INDEL annotation DF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "#append columns for info from annotation\n",
    "INDEL_annotation_DF['gene_id'] = pd.Series(locus_tag_dict)\n",
    "INDEL_annotation_DF['gene_pos'] = pd.Series(genomic_coord_dict)\n",
    "INDEL_annotation_DF['gene_category'] = pd.Series(gene_categ_dict)\n",
    "INDEL_annotation_DF['gene_name'] = pd.Series(gene_name_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 8)"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape( INDEL_annotation_DF )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>312.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt     gene_id  gene_pos gene_category  \\\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A      Rv0001     312.0     Essential   \n",
       "1         TC_1549_T  1549          TC   T  intergenic       NaN          None   \n",
       "2        TAA_1552_T  1552         TAA   T  intergenic       NaN          None   \n",
       "3         TA_1552_T  1552          TA   T  intergenic       NaN          None   \n",
       "4         T_1552_TA  1552           T  TA  intergenic       NaN          None   \n",
       "\n",
       "  gene_name  \n",
       "0      dnaA  \n",
       "1      None  \n",
       "2      None  \n",
       "3      None  \n",
       "4      None  "
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>55263</th>\n",
       "      <td>C_4411205_CA</td>\n",
       "      <td>4411205</td>\n",
       "      <td>C</td>\n",
       "      <td>CA</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55264</th>\n",
       "      <td>AG_4411208_A</td>\n",
       "      <td>4411208</td>\n",
       "      <td>AG</td>\n",
       "      <td>A</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55265</th>\n",
       "      <td>C_4411423_CTG</td>\n",
       "      <td>4411423</td>\n",
       "      <td>C</td>\n",
       "      <td>CTG</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55266</th>\n",
       "      <td>CCTG_4411459_C</td>\n",
       "      <td>4411459</td>\n",
       "      <td>CCTG</td>\n",
       "      <td>C</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55267</th>\n",
       "      <td>AC_4411486_A</td>\n",
       "      <td>4411486</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  key      pos   ref  alt     gene_id  gene_pos gene_category  \\\n",
       "55263    C_4411205_CA  4411205     C   CA  intergenic       NaN          None   \n",
       "55264    AG_4411208_A  4411208    AG    A  intergenic       NaN          None   \n",
       "55265   C_4411423_CTG  4411423     C  CTG  intergenic       NaN          None   \n",
       "55266  CCTG_4411459_C  4411459  CCTG    C  intergenic       NaN          None   \n",
       "55267    AC_4411486_A  4411486    AC    A  intergenic       NaN          None   \n",
       "\n",
       "      gene_name  \n",
       "55263      None  \n",
       "55264      None  \n",
       "55265      None  \n",
       "55266      None  \n",
       "55267      None  "
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "#update INDEL annotation file (rows of Genotype Matrix) with extra INDEL information\n",
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation_with_gene_info.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to investigate SNPs that appeared within a particular gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>10057</th>\n",
       "      <td>G_760090_GGTCGAT</td>\n",
       "      <td>760090</td>\n",
       "      <td>G</td>\n",
       "      <td>GGTCGAT</td>\n",
       "      <td>Rv0667</td>\n",
       "      <td>283.0</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "      <td>rpoB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10058</th>\n",
       "      <td>A_760095_ATGTCGT</td>\n",
       "      <td>760095</td>\n",
       "      <td>A</td>\n",
       "      <td>ATGTCGT</td>\n",
       "      <td>Rv0667</td>\n",
       "      <td>288.0</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "      <td>rpoB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10059</th>\n",
       "      <td>TCTTCGG_761077_T</td>\n",
       "      <td>761077</td>\n",
       "      <td>TCTTCGG</td>\n",
       "      <td>T</td>\n",
       "      <td>Rv0667</td>\n",
       "      <td>1270.0</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "      <td>rpoB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10060</th>\n",
       "      <td>TCGGCAC_761080_T</td>\n",
       "      <td>761080</td>\n",
       "      <td>TCGGCAC</td>\n",
       "      <td>T</td>\n",
       "      <td>Rv0667</td>\n",
       "      <td>1273.0</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "      <td>rpoB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10061</th>\n",
       "      <td>GGCACCA_761082_G</td>\n",
       "      <td>761082</td>\n",
       "      <td>GGCACCA</td>\n",
       "      <td>G</td>\n",
       "      <td>Rv0667</td>\n",
       "      <td>1275.0</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "      <td>rpoB</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    key     pos      ref      alt gene_id  gene_pos  \\\n",
       "10057  G_760090_GGTCGAT  760090        G  GGTCGAT  Rv0667     283.0   \n",
       "10058  A_760095_ATGTCGT  760095        A  ATGTCGT  Rv0667     288.0   \n",
       "10059  TCTTCGG_761077_T  761077  TCTTCGG        T  Rv0667    1270.0   \n",
       "10060  TCGGCAC_761080_T  761080  TCGGCAC        T  Rv0667    1273.0   \n",
       "10061  GGCACCA_761082_G  761082  GGCACCA        G  Rv0667    1275.0   \n",
       "\n",
       "               gene_category gene_name  \n",
       "10057  Antibiotic Resistance      rpoB  \n",
       "10058  Antibiotic Resistance      rpoB  \n",
       "10059  Antibiotic Resistance      rpoB  \n",
       "10060  Antibiotic Resistance      rpoB  \n",
       "10061  Antibiotic Resistance      rpoB  "
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF[['Rv0667' in gene_list for gene_list in INDEL_annotation_DF.gene_id]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [8] Merge isolate annotation from INDEL genotypes matrix & Save annotation files and genotypes matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load INDEL genotype matrix and Annotation Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_series = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_isolate_annotation.pkl')\n",
    "\n",
    "#load INDEL annotation file (rows of Genotype Matrix) with gene annotation information\n",
    "INDEL_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_INDEL_annotation_with_gene_info.pkl')\n",
    "\n",
    "#load Genotypes Matrix\n",
    "genotypes_array =  np.load('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered/genotypes_matrix.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428,)"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 8)"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 31428)"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Isolate Annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    SAMN13051687\n",
       "1    SAMN09100245\n",
       "2    SAMN08732238\n",
       "3    SAMN07658260\n",
       "4    SAMN03648003\n",
       "dtype: object"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_series.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load file with isolate annotation DF from the SNP analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_annotation_df_from_SNP_analysis = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lineage_1</th>\n",
       "      <th>lineage_2</th>\n",
       "      <th>lineage_3</th>\n",
       "      <th>lineage_4</th>\n",
       "      <th>lineage_5</th>\n",
       "      <th>lineage_6</th>\n",
       "      <th>lineage_7</th>\n",
       "      <th>lineage_8</th>\n",
       "      <th>lineage_9</th>\n",
       "      <th>lineage_10</th>\n",
       "      <th>lineage_11</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3558733</td>\n",
       "      <td>4.2.1.2.1.1.i3.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03648641</td>\n",
       "      <td>4.2.1.2.2.1.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03647419</td>\n",
       "      <td>3.1.1.i1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3671418</td>\n",
       "      <td>4.2.1.2.1.1.i1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN07659096</td>\n",
       "      <td>1.1.1.2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  lineage_1 lineage_2 lineage_3 lineage_4 lineage_5 lineage_6 lineage_7  \\\n",
       "0         4         2         1         2         1         1        i3   \n",
       "1         4         2         1         2         2         1         1   \n",
       "2         3         1         1        i1       NaN       NaN       NaN   \n",
       "3         4         2         1         2         1         1        i1   \n",
       "4         1         1         1         2       NaN       NaN       NaN   \n",
       "\n",
       "  lineage_8 lineage_9 lineage_10 lineage_11    isolate_ID      lineage_call  \\\n",
       "0         1       NaN        NaN        NaN  SAMEA3558733  4.2.1.2.1.1.i3.1   \n",
       "1       NaN       NaN        NaN        NaN  SAMN03648641     4.2.1.2.2.1.1   \n",
       "2       NaN       NaN        NaN        NaN  SAMN03647419          3.1.1.i1   \n",
       "3       NaN       NaN        NaN        NaN  SAMEA3671418    4.2.1.2.1.1.i1   \n",
       "4       NaN       NaN        NaN        NaN  SAMN07659096           1.1.1.2   \n",
       "\n",
       "  group  \n",
       "0    4B  \n",
       "1    4B  \n",
       "2     3  \n",
       "3    4B  \n",
       "4     1  "
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_df_from_SNP_analysis.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 14)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_df_from_SNP_analysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge isolate annotation series (from INDEL genotypes construction) with Lineage/Sub-Lineage calls annotation (from SNP genotypes construction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_annotation_DF = pd.DataFrame()\n",
    "isolate_annotation_DF['isolate_ID'] = isolate_annotation_series\n",
    "\n",
    "#merge Lineage/Sub-Lineage calls\n",
    "isolate_annotation_DF = isolate_annotation_DF.merge(isolate_annotation_df_from_SNP_analysis , how = 'inner' , on = 'isolate_ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_1</th>\n",
       "      <th>lineage_2</th>\n",
       "      <th>lineage_3</th>\n",
       "      <th>lineage_4</th>\n",
       "      <th>lineage_5</th>\n",
       "      <th>lineage_6</th>\n",
       "      <th>lineage_7</th>\n",
       "      <th>lineage_8</th>\n",
       "      <th>lineage_9</th>\n",
       "      <th>lineage_10</th>\n",
       "      <th>lineage_11</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SAMN13051687</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SAMN09100245</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>4.2.1.2.1.1.i3.2</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SAMN08732238</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SAMN07658260</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>3.1.1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SAMN03648003</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     isolate_ID lineage_1 lineage_2 lineage_3 lineage_4 lineage_5 lineage_6  \\\n",
       "0  SAMN13051687         2         2         1         1         1        i3   \n",
       "1  SAMN09100245         4         2         1         2         1         1   \n",
       "2  SAMN08732238         2         2         1         1         1       NaN   \n",
       "3  SAMN07658260         3         1         1       NaN       NaN       NaN   \n",
       "4  SAMN03648003         2         2         1         1         1       NaN   \n",
       "\n",
       "  lineage_7 lineage_8 lineage_9 lineage_10 lineage_11      lineage_call group  \n",
       "0       NaN       NaN       NaN        NaN        NaN      2.2.1.1.1.i3     2  \n",
       "1        i3         2       NaN        NaN        NaN  4.2.1.2.1.1.i3.2    4B  \n",
       "2       NaN       NaN       NaN        NaN        NaN         2.2.1.1.1     2  \n",
       "3       NaN       NaN       NaN        NaN        NaN             3.1.1     3  \n",
       "4       NaN       NaN       NaN        NaN        NaN         2.2.1.1.1     2  "
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 14)"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### INDEL Annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>312.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt     gene_id  gene_pos gene_category  \\\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A      Rv0001     312.0     Essential   \n",
       "1         TC_1549_T  1549          TC   T  intergenic       NaN          None   \n",
       "2        TAA_1552_T  1552         TAA   T  intergenic       NaN          None   \n",
       "3         TA_1552_T  1552          TA   T  intergenic       NaN          None   \n",
       "4         T_1552_TA  1552           T  TA  intergenic       NaN          None   \n",
       "\n",
       "  gene_name  \n",
       "0      dnaA  \n",
       "1      None  \n",
       "2      None  \n",
       "3      None  \n",
       "4      None  "
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 8)"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save SNP annotation file (rows of Genotype Matrix)\n",
    "INDEL_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Genotype Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       ...,\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0]], dtype=int8)"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 31428)"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save Genotypes Matrix\n",
    "np.save('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_matrix' , genotypes_array , allow_pickle = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [9] Create DataFrame with the functional annotation for INDELs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This block of code will determine the functional consequence of each INDEL in our genotypes array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load SNP annotation file (rows of Genotype Matrix)\n",
    "INDEL_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_INDEL_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACCGACGAAG_313_A</td>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>312.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TC_1549_T</td>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>TAA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TA_1552_T</td>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T_1552_TA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                key   pos         ref alt     gene_id  gene_pos gene_category  \\\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG   A      Rv0001     312.0     Essential   \n",
       "1         TC_1549_T  1549          TC   T  intergenic       NaN          None   \n",
       "2        TAA_1552_T  1552         TAA   T  intergenic       NaN          None   \n",
       "3         TA_1552_T  1552          TA   T  intergenic       NaN          None   \n",
       "4         T_1552_TA  1552           T  TA  intergenic       NaN          None   \n",
       "\n",
       "  gene_name  \n",
       "0      dnaA  \n",
       "1      None  \n",
       "2      None  \n",
       "3      None  \n",
       "4      None  "
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 8)"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a dataframe for each unique _reference_ allele - _alternate_ allele pairing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_functional_annotation_DF = INDEL_annotation_DF.loc[: , ['key' , 'pos' , 'ref' , 'alt' , 'gene_category' , 'gene_name']]\n",
    "INDEL_functional_annotation_DF.set_index('key' , inplace = True , drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>key</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>ACCGACGAAG_313_A</th>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TC_1549_T</th>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TAA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_1552_TA</th>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   pos         ref alt gene_category gene_name\n",
       "key                                                           \n",
       "ACCGACGAAG_313_A   313  ACCGACGAAG   A     Essential      dnaA\n",
       "TC_1549_T         1549          TC   T          None      None\n",
       "TAA_1552_T        1552         TAA   T          None      None\n",
       "TA_1552_T         1552          TA   T          None      None\n",
       "T_1552_TA         1552           T  TA          None      None"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_functional_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 5)"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_functional_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert **gene_category** and **gene_name** columns so that the elements are lists only when there is more than 1 label in a cell (when a Reference Position corresponds to two overlapping genes, these positions have been filtered out in the most recent version of analysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_category_list = []\n",
    "gene_name_list = []\n",
    "\n",
    "for INDEL_i in INDEL_functional_annotation_DF.index:\n",
    "    \n",
    "    gene_category_INDEL_i = INDEL_functional_annotation_DF.loc[INDEL_i , 'gene_category']\n",
    "    gene_name_INDEL_i = INDEL_functional_annotation_DF.loc[INDEL_i , 'gene_name']\n",
    "    \n",
    "    if len(gene_category_INDEL_i) == 1:\n",
    "        \n",
    "        gene_category_list.append(gene_category_INDEL_i[0])\n",
    "        gene_name_list.append(gene_name_INDEL_i[0])\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        gene_category_list.append(gene_category_INDEL_i)\n",
    "        gene_name_list.append(gene_name_INDEL_i)\n",
    "        \n",
    "INDEL_functional_annotation_DF['gene_category'] = gene_category_list\n",
    "INDEL_functional_annotation_DF['gene_name'] = gene_name_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>key</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>ACCGACGAAG_313_A</th>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TC_1549_T</th>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TAA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_1552_TA</th>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   pos         ref alt gene_category gene_name\n",
       "key                                                           \n",
       "ACCGACGAAG_313_A   313  ACCGACGAAG   A     Essential      dnaA\n",
       "TC_1549_T         1549          TC   T          None      None\n",
       "TAA_1552_T        1552         TAA   T          None      None\n",
       "TA_1552_T         1552          TA   T          None      None\n",
       "T_1552_TA         1552           T  TA          None      None"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_functional_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 5)"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_functional_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load in work-horse for functional INDEL annotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relevant Information for H37Rv sequence and INDEL functional annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "####### Collect all DNA and Amino Acid sequences corresponding to genes on H37Rv #######\n",
    "#load reference genome and reference annotation\n",
    "reference_genome = '/n/data1/hms/dbmi/farhat/bin/work-horse/bin/h37rv.fasta'\n",
    "for reference_genome in SeqIO.parse(reference_genome, \"fasta\"):\n",
    "    reference_genome.seq.alphabet = IUPAC.unambiguous_dna\n",
    "\n",
    "reference_genome_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt', '\\t').set_index('name')\n",
    "\n",
    "####### Load in dictionaries for SNP annotation #######\n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/pickled_files/dicts_for_SNP_annotation/H37Rv_gene_seq_records.pickle', 'rb') as handle:\n",
    "    ref_gene_sequences_records = pickle.load(handle)\n",
    "    \n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/pickled_files/dicts_for_SNP_annotation/H37Rv_protein_seq_records.pickle', 'rb') as handle:\n",
    "    ref_protein_sequences_records = pickle.load(handle)\n",
    "    \n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/pickled_files/dicts_for_SNP_annotation/H37Rv_coord_gene_mapping.pickle', 'rb') as handle:\n",
    "    ReferencePosition_Gene_mapping = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to annotate Intergenic INDELs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_flanking_genes_for_intergenic_region(intergenic_ref_pos): \n",
    "\n",
    "    #this function finds the genes flagging an intergenic region given a reference position\n",
    "\n",
    "    #find gene immediately in the 5' direction\n",
    "    for i in range(0 , 100000):\n",
    "\n",
    "        #move toward 5' direction\n",
    "        if ReferencePosition_Gene_mapping[intergenic_ref_pos - i] != []:\n",
    "\n",
    "            gene_to_left = ReferencePosition_Gene_mapping[intergenic_ref_pos - i][0]\n",
    "            break\n",
    "\n",
    "    #find gene immediately in the 3' direction       \n",
    "    for i in range(0 , 100000):\n",
    "\n",
    "        #move toward 3' direction\n",
    "        try:\n",
    "            if ReferencePosition_Gene_mapping[intergenic_ref_pos + i] != []:\n",
    "\n",
    "                gene_to_right = ReferencePosition_Gene_mapping[intergenic_ref_pos + i][0]\n",
    "                break\n",
    "        \n",
    "        #KeyError means we have hit the 'end' of the chromosome, the intergenic region at then end of H37Rv in 5' > 3' orientation \n",
    "        #since TB chromosome is circular the gene to the 'right' is Rv0001    \n",
    "        except KeyError:\n",
    "            \n",
    "            gene_to_right = 'Rv0001'\n",
    "            break\n",
    "    \n",
    "    return gene_to_left + '_' + gene_to_right"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to determine whether INDELs are *Inframe* or *Frameshift*; Return codon position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "def INDEL_type(ref_allele_i , ref_seq_position , alt_allele_i):\n",
    "    \n",
    "    '''\n",
    "    This function takes as input a reference position on H37Rv located within a \n",
    "    gene, a reference allelel and an alternate allele and returns whether the indel\n",
    "    would correspond to an inframe or frameshift indel.\n",
    "    \n",
    "    '''\n",
    "    gene_intergenic_id_list = []\n",
    "    genomic_coord_list = []\n",
    "    IF_FS_list = []\n",
    "    INS_DEL_list = []\n",
    "    codon_coord_list = []\n",
    "    \n",
    "    #find the gene that INDEL occurs on\n",
    "    \n",
    "    #check list corresponding to H37Rv coordinate to see if there are any genes associated with RefPosition\n",
    "    if len(ReferencePosition_Gene_mapping[ref_seq_position]) > 0:\n",
    "\n",
    "        #iterate through all genes that ReferencePosition is mapped to (i.e. SNP might correspond to 2 genes)\n",
    "        for gene_intergenic_id in ReferencePosition_Gene_mapping[ref_seq_position]:\n",
    "\n",
    "            #find genomic coordinate of SNP relative to gene (subtract 1 since reference seq starts counting at 1)\n",
    "            gene_relative_coord = (ref_seq_position - 1) - min( reference_genome_annotation.loc[gene_intergenic_id , 'chromStart'] , reference_genome_annotation.loc[gene_intergenic_id , 'chromEnd'] )\n",
    "\n",
    "            #find out if indel is inframe or frameshift\n",
    "            indel_length = abs(len(ref_allele_i) - len(alt_allele_i))\n",
    "            \n",
    "            if indel_length % 3 == 0:\n",
    "                INDEL_type = 'inframe'\n",
    "                \n",
    "            elif indel_length % 3 != 0:\n",
    "                INDEL_type = 'frameshift'\n",
    "                \n",
    "            #find out whether insertion or deletion\n",
    "            if len(ref_allele_i) > len(alt_allele_i):\n",
    "                INS_or_DEL = 'del'\n",
    "                \n",
    "            elif len(ref_allele_i) < len(alt_allele_i):\n",
    "                INS_or_DEL = 'ins'\n",
    "            \n",
    "            #get the codon number where the SNP occurs within\n",
    "            ## take the genomic coordinate (relative to the gene, in the 5' to 3' direction), divide by 3, then take the ceiling of this number (will be fraction if SNP occurs in 1st or 2nd position on codon)\n",
    "            strand = reference_genome_annotation.loc[gene_intergenic_id, 'strand']\n",
    "            if strand == '+':\n",
    "                 genomic_5_to_3_coord = (ref_seq_position) - reference_genome_annotation.loc[gene_intergenic_id , 'chromStart']\n",
    "\n",
    "            elif strand == '-':\n",
    "                 genomic_5_to_3_coord = (reference_genome_annotation.loc[gene_intergenic_id , 'chromEnd']) - (ref_seq_position-1)\n",
    "\n",
    "            codon_coord = int(np.ceil( float( genomic_5_to_3_coord) / 3.0 ))\n",
    "            \n",
    "            #store relevant info in lists    \n",
    "            gene_intergenic_id_list.append(gene_intergenic_id)\n",
    "            genomic_coord_list.append(genomic_5_to_3_coord)\n",
    "            IF_FS_list.append(INDEL_type)\n",
    "            INS_DEL_list.append(INS_or_DEL)\n",
    "            codon_coord_list.append(codon_coord)\n",
    "    \n",
    "    #if no gene in H37Rv corresponds to the Reference Position for INDEL, then INDEL must be intergenic\n",
    "    else:\n",
    "        \n",
    "        gene_intergenic_id = find_flanking_genes_for_intergenic_region(ref_seq_position)\n",
    "        \n",
    "        #find out if indel is inframe or frameshift\n",
    "        indel_length = abs(len(ref_allele_i) - len(alt_allele_i))\n",
    "\n",
    "        if indel_length % 3 == 0:\n",
    "            INDEL_type = 'inframe'\n",
    "\n",
    "        elif indel_length % 3 != 0:\n",
    "            INDEL_type = 'frameshift'\n",
    "            \n",
    "        #find out whether insertion or deletion\n",
    "        if len(ref_allele_i) > len(alt_allele_i):\n",
    "            INS_or_DEL = 'del'\n",
    "\n",
    "        elif len(ref_allele_i) < len(alt_allele_i):\n",
    "            INS_or_DEL = 'ins'\n",
    "                \n",
    "        genomic_5_to_3_coord = 'None'\n",
    "        codon_coord = 'None'\n",
    "        \n",
    "        #store relevant info in lists    \n",
    "        gene_intergenic_id_list.append(gene_intergenic_id)\n",
    "        genomic_coord_list.append('None')\n",
    "        IF_FS_list.append('I')\n",
    "        INS_DEL_list.append(INS_or_DEL)\n",
    "        codon_coord_list.append('None')\n",
    "    \n",
    "    #if there is only a single gene associated with this INDEL, just return the individual elememts\n",
    "    if len(gene_intergenic_id_list) == 1:\n",
    "        return [gene_intergenic_id , genomic_5_to_3_coord , INS_or_DEL , INDEL_type , codon_coord]\n",
    "    \n",
    "    #else if there are two genes associated with this SNP, return elements for each SNP annotation in a list\n",
    "    elif len(gene_intergenic_id_list) > 1:\n",
    "        return [gene_intergenic_id_list , genomic_coord_list , INS_DEL_list , IF_FS_list , codon_coord_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Iterate through each unique Reference Allele - Alternate Allele pair and get relevant annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create lists for each new column of the INDEL functional annotation DF\n",
    "gene_intergenic_id_list = []\n",
    "gene_pos_list = []\n",
    "INS_DEL_list = []\n",
    "INDEL_type_list = []\n",
    "codon_coord_list = []\n",
    "\n",
    "for INDEL_i in INDEL_functional_annotation_DF.index:\n",
    "    \n",
    "    #get the reference allele , reference position & alternate allele from index key\n",
    "    ref_allele_i , ref_pos_i , alt_allele_i = INDEL_i.split('_')\n",
    "    \n",
    "    #get functional annotation for INDEL\n",
    "    gene_intergenic_id_i , gene_pos_i , INS_or_DEL_i , INDEL_type_i , codon_coord_i = INDEL_type(ref_allele_i , int(ref_pos_i) , alt_allele_i)\n",
    "    \n",
    "    #store functional annotation info in lists\n",
    "    gene_intergenic_id_list.append(gene_intergenic_id_i)\n",
    "    gene_pos_list.append(gene_pos_i)\n",
    "    INS_DEL_list.append(INS_or_DEL_i)\n",
    "    INDEL_type_list.append(INDEL_type_i)\n",
    "    codon_coord_list.append(codon_coord_i)\n",
    "    \n",
    "#append new columns to INDEL functional annotation\n",
    "INDEL_functional_annotation_DF['gene_id'] = gene_intergenic_id_list\n",
    "INDEL_functional_annotation_DF['gene_pos'] = gene_pos_list\n",
    "INDEL_functional_annotation_DF['ins_del'] = INS_DEL_list\n",
    "INDEL_functional_annotation_DF['INDEL_type'] = INDEL_type_list\n",
    "INDEL_functional_annotation_DF['codon_pos'] = codon_coord_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55268, 10)"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_functional_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>ins_del</th>\n",
       "      <th>INDEL_type</th>\n",
       "      <th>codon_pos</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>key</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>ACCGACGAAG_313_A</th>\n",
       "      <td>313</td>\n",
       "      <td>ACCGACGAAG</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>313</td>\n",
       "      <td>del</td>\n",
       "      <td>inframe</td>\n",
       "      <td>105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TC_1549_T</th>\n",
       "      <td>1549</td>\n",
       "      <td>TC</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv0001_Rv0002</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TAA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TAA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv0001_Rv0002</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TA_1552_T</th>\n",
       "      <td>1552</td>\n",
       "      <td>TA</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv0001_Rv0002</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_1552_TA</th>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TA</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv0001_Rv0002</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   pos         ref alt gene_category gene_name        gene_id  \\\n",
       "key                                                                             \n",
       "ACCGACGAAG_313_A   313  ACCGACGAAG   A     Essential      dnaA         Rv0001   \n",
       "TC_1549_T         1549          TC   T          None      None  Rv0001_Rv0002   \n",
       "TAA_1552_T        1552         TAA   T          None      None  Rv0001_Rv0002   \n",
       "TA_1552_T         1552          TA   T          None      None  Rv0001_Rv0002   \n",
       "T_1552_TA         1552           T  TA          None      None  Rv0001_Rv0002   \n",
       "\n",
       "                 gene_pos ins_del  INDEL_type codon_pos  \n",
       "key                                                      \n",
       "ACCGACGAAG_313_A      313     del     inframe       105  \n",
       "TC_1549_T            None     del  frameshift      None  \n",
       "TAA_1552_T           None     del  frameshift      None  \n",
       "TA_1552_T            None     del  frameshift      None  \n",
       "T_1552_TA            None     ins  frameshift      None  "
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_functional_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pickle functional INDEL annotation for downstream analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Note__: \n",
    "\n",
    "The gene position/coordinate for the original INDEL annotation file (whose rows map to the rows of the Genotypes Matrix) starts counting from 0, this is to accomodate the functional INDEL annotation. \n",
    "\n",
    "The gene position/coordinate for the functional INDEL annotation (above) starts counting from 1 which is more intuitive for reporting purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEL_functional_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_INDEL_functional_annotation.pkl')"
   ]
  }
 ],
 "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}