(D) whiB7 and whiB7 promoter analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook was created to 1. collect (SNP & INDEL) variants in *whiB7* and the *whiB7* promoter and 2. detect instances of co-occurence of LoF variants in *whiB7* and the *whiB7* promoter"
   ]
  },
  {
   "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",
    "from Bio import Phylo\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\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *whiB7*\n",
    "\n",
    "promoter: 3,568,680 - 3,569,082\n",
    "\n",
    "coding: 3,568,401 - 3,568,679\n",
    "\n",
    "H37Rv region: 3,568,401 - 3,569,082"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "H37Rv_region_start = 3568401\n",
    "H37Rv_region_end = 3569082\n",
    "regulator_id = 'Rv3197A_Rv3198c'\n",
    "regulon_id = 'Rv3197A'\n",
    "gene_symbol = 'whiB7'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Function* to convert lineage-calls to new/simpler lineage call scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_lineage_calls(isolate_annotation_DF):\n",
    "\n",
    "    lineage_hierarchincal_to_lineage_map_df = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/lineage_hierachical_to_lineage_map.csv')\n",
    "    lineage_hierarchincal_to_lineage_map_series = pd.Series(list(lineage_hierarchincal_to_lineage_map_df.loc[:, 'lineage']), index = lineage_hierarchincal_to_lineage_map_df.lineage_hierarchical)\n",
    "\n",
    "    lineage_call_update_list = []\n",
    "\n",
    "    for hierarchical_lineage_call in isolate_annotation_DF.lineage_call:\n",
    "\n",
    "        # if lineage call has an updated shorter call\n",
    "        if hierarchical_lineage_call in lineage_hierarchincal_to_lineage_map_series.index:\n",
    "            lineage_call_update_list.append(lineage_hierarchincal_to_lineage_map_series[hierarchical_lineage_call])\n",
    "\n",
    "        # if lineage call doesn't have an updated shorter call\n",
    "        else:\n",
    "            lineage_call_update_list.append(hierarchical_lineage_call)\n",
    "\n",
    "    # replace the hierarchical calls with the shorter ones\n",
    "    isolate_annotation_DF.loc[:,'lineage_call'] = lineage_call_update_list\n",
    "    \n",
    "    return isolate_annotation_DF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] INDELs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [1.1] Load INDEL genotype matrix and Annotation Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 31,428 isolates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "INDEL_isolate_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')\n",
    "INDEL_isolate_annotation_DF = INDEL_isolate_annotation_DF.loc[:, ['isolate_ID','lineage_call']] #drop columns\n",
    "INDEL_isolate_annotation_DF = convert_lineage_calls(INDEL_isolate_annotation_DF)\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/homoplasy_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_INDEL_functional_annotation.pkl')\n",
    "INDEL_annotation_DF.reset_index(inplace = True, drop = False)\n",
    "\n",
    "#load Genotypes Matrix\n",
    "INDEL_genotypes_array =  np.load('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape_indels/Genotypes_Filtered_2/genotypes_matrix.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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_call</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SAMN13051687</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SAMN09100245</td>\n",
       "      <td>4.3.i4.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SAMN08732238</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SAMN07658260</td>\n",
       "      <td>3.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SAMN03648003</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     isolate_ID  lineage_call\n",
       "0  SAMN13051687  2.2.1.1.1.i3\n",
       "1  SAMN09100245      4.3.i4.1\n",
       "2  SAMN08732238     2.2.1.1.1\n",
       "3  SAMN07658260         3.1.1\n",
       "4  SAMN03648003     2.2.1.1.1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_isolate_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "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_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",
       "  </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>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>1</th>\n",
       "      <td>TC_1549_T</td>\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>2</th>\n",
       "      <td>T_1552_TAA</td>\n",
       "      <td>1552</td>\n",
       "      <td>T</td>\n",
       "      <td>TAA</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",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TAA_1552_T</td>\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>4</th>\n",
       "      <td>T_1552_TA</td>\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": [
       "                key   pos         ref  alt gene_category gene_name  \\\n",
       "0  ACCGACGAAG_313_A   313  ACCGACGAAG    A     Essential      dnaA   \n",
       "1         TC_1549_T  1549          TC    T          None      None   \n",
       "2        T_1552_TAA  1552           T  TAA          None      None   \n",
       "3        TAA_1552_T  1552         TAA    T          None      None   \n",
       "4         T_1552_TA  1552           T   TA          None      None   \n",
       "\n",
       "         gene_id gene_pos ins_del  INDEL_type codon_pos  \n",
       "0         Rv0001      313     del     inframe       105  \n",
       "1  Rv0001_Rv0002     None     del  frameshift      None  \n",
       "2  Rv0001_Rv0002     None     ins  frameshift      None  \n",
       "3  Rv0001_Rv0002     None     del  frameshift      None  \n",
       "4  Rv0001_Rv0002     None     ins  frameshift      None  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50260, 11)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50260, 31428)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### additional 12 eis C-14T mutants with AG MICs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_DF_extra_strains = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/eis_promoter_mutant_strains_from_Milan/Genotypes/Genotypes_Filtered/genotypes_isolate_annotation.pkl')\n",
    "isolate_annotation_DF_extra_strains = isolate_annotation_DF_extra_strains.loc[:, ['isolate_ID','lineage_call']]\n",
    "isolate_annotation_DF_extra_strains = convert_lineage_calls(isolate_annotation_DF_extra_strains)\n",
    "\n",
    "#load Genotypes Matrix\n",
    "genotypes_array_extra_strains =  np.load('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/eis_promoter_mutant_strains_from_Milan/Genotypes_indels/Genotypes_Filtered/genotypes_matrix.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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_call</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>168-19</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>622-19</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>655-19</td>\n",
       "      <td>4.3.i3.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>IT1070</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>IT123</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  isolate_ID  lineage_call\n",
       "0     168-19     2.2.1.1.1\n",
       "1     622-19  2.2.1.1.1.i3\n",
       "2     655-19      4.3.i3.1\n",
       "3     IT1070  2.2.1.1.1.i3\n",
       "4      IT123     2.2.1.1.1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_DF_extra_strains.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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, ..., 9, 0, 0],\n",
       "       [0, 0, 0, ..., 0, 0, 0]], dtype=int8)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array_extra_strains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50260, 12)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_array_extra_strains)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge *isolate annotation files* and *genotypes matrices* together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge isolate annotation dataframes\n",
    "INDEL_isolate_annotation_DF = INDEL_isolate_annotation_DF.append(isolate_annotation_DF_extra_strains)\n",
    "INDEL_isolate_annotation_DF.reset_index(inplace = True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31440, 2)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge genotypes matrices\n",
    "INDEL_genotypes_array = np.hstack((INDEL_genotypes_array, genotypes_array_extra_strains))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50260, 31440)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [1.2] Subset to targeted chomosomal region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep only indels that occur within chromosomal region\n",
    "indels_in_targeted_region = list(INDEL_annotation_DF[(INDEL_annotation_DF.pos >= H37Rv_region_start) & (INDEL_annotation_DF.pos <= H37Rv_region_end)].key)\n",
    "\n",
    "# construct boolean filter\n",
    "indels_to_keep_filter = [indel_i in indels_in_targeted_region for indel_i in INDEL_annotation_DF.key]\n",
    "\n",
    "# filter Genotypes Matrix\n",
    "INDEL_genotypes_array = INDEL_genotypes_array[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": 21,
   "metadata": {
    "scrolled": true
   },
   "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_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",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C_3568429_CCG</td>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C_3568434_CAA</td>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>82</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>GCTT_3568439_G</td>\n",
       "      <td>3568439</td>\n",
       "      <td>GCTT</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>241</td>\n",
       "      <td>del</td>\n",
       "      <td>inframe</td>\n",
       "      <td>81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C_3568487_CG</td>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GC_3568488_G</td>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              key      pos   ref  alt  gene_category gene_name  gene_id  \\\n",
       "0   C_3568429_CCG  3568429     C  CCG  Non-Essential     whiB7  Rv3197A   \n",
       "1   C_3568434_CAA  3568434     C  CAA  Non-Essential     whiB7  Rv3197A   \n",
       "2  GCTT_3568439_G  3568439  GCTT    G  Non-Essential     whiB7  Rv3197A   \n",
       "3    C_3568487_CG  3568487     C   CG  Non-Essential     whiB7  Rv3197A   \n",
       "4    GC_3568488_G  3568488    GC    G  Non-Essential     whiB7  Rv3197A   \n",
       "\n",
       "  gene_pos ins_del  INDEL_type codon_pos  \n",
       "0      251     ins  frameshift        84  \n",
       "1      246     ins  frameshift        82  \n",
       "2      241     del     inframe        81  \n",
       "3      193     ins  frameshift        65  \n",
       "4      192     del  frameshift        64  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26, 11)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26, 31440)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2] SNPs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 31,428 isolates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [2.1] Load SNP genotype matrix and Annotation Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "SNP_isolate_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')\n",
    "SNP_isolate_annotation_DF = SNP_isolate_annotation_DF.loc[:, ['isolate_ID','lineage_call']] #drop columns\n",
    "SNP_isolate_annotation_DF = convert_lineage_calls(SNP_isolate_annotation_DF)\n",
    "\n",
    "#load INDEL annotation file (rows of Genotype Matrix) with gene annotation information\n",
    "SNP_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_SNP_annotation.pkl')\n",
    "SNP_annotation_DF = SNP_annotation_DF.loc[:, ['pos','ref','alt']] # drop columns\n",
    "SNP_annotation_DF.reset_index(inplace = True, drop = False)\n",
    "\n",
    "#load Genotypes Matrix\n",
    "SNP_genotypes_array =  np.load('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_matrix.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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_call</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SAMEA3558733</td>\n",
       "      <td>4.3.i4.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SAMN03648641</td>\n",
       "      <td>4.4.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SAMN03647419</td>\n",
       "      <td>3.1.1.i1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SAMEA3671418</td>\n",
       "      <td>4.3.i2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SAMN07659096</td>\n",
       "      <td>1.1.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     isolate_ID lineage_call\n",
       "0  SAMEA3558733     4.3.i4.2\n",
       "1  SAMN03648641      4.4.1.1\n",
       "2  SAMN03647419     3.1.1.i1\n",
       "3  SAMEA3671418       4.3.i2\n",
       "4  SAMN07659096        1.1.3"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_isolate_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 2)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "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>index</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>0</td>\n",
       "      <td>48</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>64</td>\n",
       "      <td>G</td>\n",
       "      <td>[C]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>67</td>\n",
       "      <td>G</td>\n",
       "      <td>[A, T]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>69</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>71</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   index  pos ref     alt\n",
       "0      0   48   C     [T]\n",
       "1      1   64   G     [C]\n",
       "2      2   67   G  [A, T]\n",
       "3      3   69   C     [T]\n",
       "4      4   71   C     [T]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 4)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1, 1, ..., 1, 1, 1],\n",
       "       [2, 2, 2, ..., 2, 2, 2],\n",
       "       [2, 2, 2, ..., 2, 2, 2],\n",
       "       ...,\n",
       "       [1, 9, 9, ..., 1, 1, 1],\n",
       "       [1, 9, 1, ..., 1, 1, 1],\n",
       "       [1, 9, 1, ..., 1, 1, 1]], dtype=int8)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 31428)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### additional 12 eis C-14T mutants with AG MICs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_DF_extra_strains = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/eis_promoter_mutant_strains_from_Milan/Genotypes/Genotypes_Filtered/genotypes_isolate_annotation.pkl')\n",
    "isolate_annotation_DF_extra_strains = isolate_annotation_DF_extra_strains.loc[:, ['isolate_ID','lineage_call']]\n",
    "isolate_annotation_DF_extra_strains = convert_lineage_calls(isolate_annotation_DF_extra_strains)\n",
    "\n",
    "#load Genotypes Matrix\n",
    "genotypes_array_extra_strains =  np.load('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/eis_promoter_mutant_strains_from_Milan/Genotypes/Genotypes_Filtered/genotypes_matrix.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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_call</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>168-19</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>622-19</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>655-19</td>\n",
       "      <td>4.3.i3.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>IT1070</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>IT123</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  isolate_ID  lineage_call\n",
       "0     168-19     2.2.1.1.1\n",
       "1     622-19  2.2.1.1.1.i3\n",
       "2     655-19      4.3.i3.1\n",
       "3     IT1070  2.2.1.1.1.i3\n",
       "4      IT123     2.2.1.1.1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_DF_extra_strains.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1, 1, ..., 1, 1, 1],\n",
       "       [2, 2, 2, ..., 2, 2, 2],\n",
       "       [2, 2, 2, ..., 2, 2, 2],\n",
       "       ...,\n",
       "       [1, 1, 1, ..., 1, 1, 1],\n",
       "       [1, 1, 1, ..., 1, 1, 1],\n",
       "       [1, 1, 1, ..., 1, 1, 1]], dtype=int8)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array_extra_strains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 12)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(genotypes_array_extra_strains)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge *isolate annotation files* and *genotypes matrices* together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge isolate annotation dataframes\n",
    "SNP_isolate_annotation_DF = SNP_isolate_annotation_DF.append(isolate_annotation_DF_extra_strains)\n",
    "SNP_isolate_annotation_DF.reset_index(inplace = True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31440, 2)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge genotypes matrices\n",
    "SNP_genotypes_array = np.hstack((SNP_genotypes_array, genotypes_array_extra_strains))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 31440)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [2.2] Subset to targeted chomosomal region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep only SNPs that occur within chromosomal region\n",
    "SNPs_in_targeted_region = list(SNP_annotation_DF[(SNP_annotation_DF.pos >= H37Rv_region_start) & (SNP_annotation_DF.pos <= H37Rv_region_end)].pos)\n",
    "\n",
    "# construct boolean filter\n",
    "SNPs_to_keep_filter = [SNP_i in SNPs_in_targeted_region for SNP_i in SNP_annotation_DF.pos]\n",
    "\n",
    "# filter Genotypes Matrix\n",
    "SNP_genotypes_array = SNP_genotypes_array[SNPs_to_keep_filter , :]\n",
    "\n",
    "# filter SNP annotation file\n",
    "SNP_annotation_DF = SNP_annotation_DF[SNPs_to_keep_filter]\n",
    "SNP_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)\n",
    "SNP_annotation_DF.loc[:,'array_index'] = SNP_annotation_DF.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "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>index</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>array_index</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>631520</td>\n",
       "      <td>3568406</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>631521</td>\n",
       "      <td>3568410</td>\n",
       "      <td>A</td>\n",
       "      <td>[C]</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>631522</td>\n",
       "      <td>3568417</td>\n",
       "      <td>T</td>\n",
       "      <td>[G]</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>631523</td>\n",
       "      <td>3568427</td>\n",
       "      <td>G</td>\n",
       "      <td>[C]</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>631524</td>\n",
       "      <td>3568435</td>\n",
       "      <td>G</td>\n",
       "      <td>[A]</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    index      pos ref  alt  array_index\n",
       "0  631520  3568406   C  [T]            0\n",
       "1  631521  3568410   A  [C]            1\n",
       "2  631522  3568417   T  [G]            2\n",
       "3  631523  3568427   G  [C]            3\n",
       "4  631524  3568435   G  [A]            4"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(150, 5)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Re-format SNP annotation DataFrame so that each alternate allele is its own row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNP_annotation_DF_mod = pd.DataFrame(columns = SNP_annotation_DF.columns)\n",
    "new_SNP_index = 0\n",
    "for SNP_i_row in SNP_annotation_DF.index:\n",
    "    \n",
    "    SNP_i_info = SNP_annotation_DF.loc[SNP_i_row, :]\n",
    "    \n",
    "    #iterate through alternate alleles, store a row for each one\n",
    "    for SNP_i_alt_allele in SNP_i_info.alt:\n",
    "        \n",
    "        SNP_annotation_DF_mod.loc[new_SNP_index, :] = list(SNP_i_info.loc[['index','pos','ref']]) + [SNP_i_alt_allele] + list(SNP_i_info.loc[['array_index']])\n",
    "        new_SNP_index += 1\n",
    "\n",
    "SNP_annotation_DF = SNP_annotation_DF_mod\n",
    "#create new column to merge func annotation on\n",
    "SNP_annotation_DF.loc[:, 'key'] = [str(SNP_i_pos) + '_' + SNP_i_alt for SNP_i_pos, SNP_i_alt in zip(SNP_annotation_DF.pos, SNP_annotation_DF.alt)]"
   ]
  },
  {
   "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>index</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>array_index</th>\n",
       "      <th>key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>631520</td>\n",
       "      <td>3568406</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>0</td>\n",
       "      <td>3568406_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>631521</td>\n",
       "      <td>3568410</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>1</td>\n",
       "      <td>3568410_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>631522</td>\n",
       "      <td>3568417</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>2</td>\n",
       "      <td>3568417_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>631523</td>\n",
       "      <td>3568427</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>3</td>\n",
       "      <td>3568427_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>631524</td>\n",
       "      <td>3568435</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>4</td>\n",
       "      <td>3568435_A</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    index      pos ref alt array_index        key\n",
       "0  631520  3568406   C   T           0  3568406_T\n",
       "1  631521  3568410   A   C           1  3568410_C\n",
       "2  631522  3568417   T   G           2  3568417_G\n",
       "3  631523  3568427   G   C           3  3568427_C\n",
       "4  631524  3568435   G   A           4  3568435_A"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(160, 6)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Merge functional annotation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNP_func_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_SNP_functional_annotation.pkl')\n",
    "SNP_func_annotation_DF.loc[:, 'key'] = SNP_func_annotation_DF.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "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>SNP_type</th>\n",
       "      <th>AA_change</th>\n",
       "      <th>key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>48_T</th>\n",
       "      <td>48</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>48</td>\n",
       "      <td>S</td>\n",
       "      <td>V16V</td>\n",
       "      <td>48_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64_C</th>\n",
       "      <td>64</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>64</td>\n",
       "      <td>N</td>\n",
       "      <td>G22R</td>\n",
       "      <td>64_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67_A</th>\n",
       "      <td>67</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>67</td>\n",
       "      <td>N</td>\n",
       "      <td>D23N</td>\n",
       "      <td>67_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67_T</th>\n",
       "      <td>67</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>67</td>\n",
       "      <td>N</td>\n",
       "      <td>D23Y</td>\n",
       "      <td>67_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69_T</th>\n",
       "      <td>69</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>69</td>\n",
       "      <td>S</td>\n",
       "      <td>D23D</td>\n",
       "      <td>69_T</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      pos ref alt gene_category gene_name gene_id gene_pos SNP_type AA_change  \\\n",
       "48_T   48   C   T     Essential      dnaA  Rv0001       48        S      V16V   \n",
       "64_C   64   G   C     Essential      dnaA  Rv0001       64        N      G22R   \n",
       "67_A   67   G   A     Essential      dnaA  Rv0001       67        N      D23N   \n",
       "67_T   67   G   T     Essential      dnaA  Rv0001       67        N      D23Y   \n",
       "69_T   69   C   T     Essential      dnaA  Rv0001       69        S      D23D   \n",
       "\n",
       "       key  \n",
       "48_T  48_T  \n",
       "64_C  64_C  \n",
       "67_A  67_A  \n",
       "67_T  67_T  \n",
       "69_T  69_T  "
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_func_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNP_annotation_DF = SNP_annotation_DF.merge(SNP_func_annotation_DF.loc[: , ['gene_category','gene_name','gene_id','gene_pos','SNP_type','AA_change','key']], how = 'left', on = 'key')\n",
    "SNP_annotation_DF = SNP_annotation_DF.loc[:, ['key','pos','ref','alt','gene_category','gene_name','gene_id','gene_pos','SNP_type','AA_change','array_index']] #re-arrange columns\n",
    "SNP_annotation_DF.rename(columns={'SNP_type':'mut_type', 'AA_change':'codon_pos'}, inplace = True) #change column names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "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_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>array_index</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3568406_T</td>\n",
       "      <td>3568406</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>274</td>\n",
       "      <td>N</td>\n",
       "      <td>A92T</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3568410_C</td>\n",
       "      <td>3568410</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>270</td>\n",
       "      <td>S</td>\n",
       "      <td>A90A</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3568417_G</td>\n",
       "      <td>3568417</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>263</td>\n",
       "      <td>N</td>\n",
       "      <td>K88T</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3568427_C</td>\n",
       "      <td>3568427</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>253</td>\n",
       "      <td>N</td>\n",
       "      <td>R85G</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3568435_A</td>\n",
       "      <td>3568435</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>245</td>\n",
       "      <td>N</td>\n",
       "      <td>P82L</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         key      pos ref alt  gene_category gene_name  gene_id gene_pos  \\\n",
       "0  3568406_T  3568406   C   T  Non-Essential     whiB7  Rv3197A      274   \n",
       "1  3568410_C  3568410   A   C  Non-Essential     whiB7  Rv3197A      270   \n",
       "2  3568417_G  3568417   T   G  Non-Essential     whiB7  Rv3197A      263   \n",
       "3  3568427_C  3568427   G   C  Non-Essential     whiB7  Rv3197A      253   \n",
       "4  3568435_A  3568435   G   A  Non-Essential     whiB7  Rv3197A      245   \n",
       "\n",
       "  mut_type codon_pos array_index  \n",
       "0        N      A92T           0  \n",
       "1        S      A90A           1  \n",
       "2        N      K88T           2  \n",
       "3        N      R85G           3  \n",
       "4        N      P82L           4  "
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(160, 11)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1, 1, ..., 1, 1, 1],\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [3, 3, 3, ..., 3, 3, 3],\n",
       "       ...,\n",
       "       [0, 0, 0, ..., 0, 0, 0],\n",
       "       [1, 1, 1, ..., 1, 1, 1],\n",
       "       [1, 1, 1, ..., 1, 1, 1]], dtype=int8)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(150, 31440)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3] Pull info for isolates with mutant alleles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [3.1] Get genotypes for SNPs of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Function* to get the genotypes for a specifc SNP and get attributes for isolates that support SNP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_genotypes_for_SNP(SNP_i_pos, SNP_i_alt_allele):\n",
    "    \n",
    "    #get the number code that corresponds to the alt allele, 9 indicated bad quality call\n",
    "    base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "    SNP_i_alt_allele_code = base_code_dict[SNP_i_alt_allele]\n",
    "    SNP_i_genotypes = SNP_genotypes_array[SNP_annotation_DF[SNP_annotation_DF.pos == SNP_i_pos].array_index.values[0] , :] \n",
    "\n",
    "    #count the number of isolates that support this SNP call\n",
    "    num_isolates_SNP_i = np.sum(SNP_i_genotypes == SNP_i_alt_allele_code)\n",
    "\n",
    "    #create a boolean filter and extract information for the isolates that support this call\n",
    "    SNP_i_isolate_filter = SNP_i_genotypes == SNP_i_alt_allele_code\n",
    "    isolate_annotation_SNP_i = SNP_isolate_annotation_DF[SNP_i_isolate_filter]\n",
    "\n",
    "    #find number of sub-lineages that have SNP support in at least 1 isolate\n",
    "    num_sublineages_with_SNP_i = len(set(list(isolate_annotation_SNP_i.lineage_call.values)))\n",
    "    \n",
    "    #get list of sublineages w/ at least 1 isolates that supports SNP call\n",
    "    sublineages_with_SNP_count_dict = Counter(isolate_annotation_SNP_i.lineage_call.values)\n",
    "    sublineages_with_SNP_list = []\n",
    "    for sublineage_i in sublineages_with_SNP_count_dict.keys():\n",
    "\n",
    "        sublineage_i_with_SNP = sublineage_i + '({0})'.format(str(sublineages_with_SNP_count_dict[sublineage_i]))\n",
    "        sublineages_with_SNP_list.append(sublineage_i_with_SNP)\n",
    "\n",
    "    sublineages_with_SNP = ' - '.join(sublineages_with_SNP_list)\n",
    "\n",
    "    return [num_isolates_SNP_i , num_sublineages_with_SNP_i , sublineages_with_SNP]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "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_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>array_index</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3568406_T</td>\n",
       "      <td>3568406</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>274</td>\n",
       "      <td>N</td>\n",
       "      <td>A92T</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3568410_C</td>\n",
       "      <td>3568410</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>270</td>\n",
       "      <td>S</td>\n",
       "      <td>A90A</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3568417_G</td>\n",
       "      <td>3568417</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>263</td>\n",
       "      <td>N</td>\n",
       "      <td>K88T</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3568427_C</td>\n",
       "      <td>3568427</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>253</td>\n",
       "      <td>N</td>\n",
       "      <td>R85G</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3568435_A</td>\n",
       "      <td>3568435</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>245</td>\n",
       "      <td>N</td>\n",
       "      <td>P82L</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         key      pos ref alt  gene_category gene_name  gene_id gene_pos  \\\n",
       "0  3568406_T  3568406   C   T  Non-Essential     whiB7  Rv3197A      274   \n",
       "1  3568410_C  3568410   A   C  Non-Essential     whiB7  Rv3197A      270   \n",
       "2  3568417_G  3568417   T   G  Non-Essential     whiB7  Rv3197A      263   \n",
       "3  3568427_C  3568427   G   C  Non-Essential     whiB7  Rv3197A      253   \n",
       "4  3568435_A  3568435   G   A  Non-Essential     whiB7  Rv3197A      245   \n",
       "\n",
       "  mut_type codon_pos array_index  \n",
       "0        N      A92T           0  \n",
       "1        S      A90A           1  \n",
       "2        N      K88T           2  \n",
       "3        N      R85G           3  \n",
       "4        N      P82L           4  "
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(160, 11)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_isolates_with_SNP = []\n",
    "num_sublineages_with_SNP = []\n",
    "sublineages_with_SNP = []\n",
    "\n",
    "for SNP_i_pos, SNP_i_alt_allele in zip(SNP_annotation_DF.pos, SNP_annotation_DF.alt):\n",
    "    \n",
    "    num_isolates_with_SNP_i , num_sublineages_with_SNP_i , sublineages_with_SNP_i = get_genotypes_for_SNP(SNP_i_pos, SNP_i_alt_allele)\n",
    "    \n",
    "    num_isolates_with_SNP.append(num_isolates_with_SNP_i)\n",
    "    num_sublineages_with_SNP.append(num_sublineages_with_SNP_i)\n",
    "    sublineages_with_SNP.append(sublineages_with_SNP_i)\n",
    "    \n",
    "SNP_annotation_DF.loc[: , 'num_isolates'] = num_isolates_with_SNP\n",
    "SNP_annotation_DF.loc[: , 'num_sublineages'] = num_sublineages_with_SNP\n",
    "SNP_annotation_DF.loc[: , 'sublineages'] = sublineages_with_SNP\n",
    "\n",
    "#drop SNPs present in 0 isolates\n",
    "SNP_annotation_DF = SNP_annotation_DF[SNP_annotation_DF.num_isolates > 0]\n",
    "SNP_annotation_DF.reset_index(inplace = True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "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_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>array_index</th>\n",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3568406_T</td>\n",
       "      <td>3568406</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>274</td>\n",
       "      <td>N</td>\n",
       "      <td>A92T</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i4.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3568410_C</td>\n",
       "      <td>3568410</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>270</td>\n",
       "      <td>S</td>\n",
       "      <td>A90A</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3568417_G</td>\n",
       "      <td>3568417</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>263</td>\n",
       "      <td>N</td>\n",
       "      <td>K88T</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3568427_C</td>\n",
       "      <td>3568427</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>253</td>\n",
       "      <td>N</td>\n",
       "      <td>R85G</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3568435_A</td>\n",
       "      <td>3568435</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>245</td>\n",
       "      <td>N</td>\n",
       "      <td>P82L</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1.i2(1)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         key      pos ref alt  gene_category gene_name  gene_id gene_pos  \\\n",
       "0  3568406_T  3568406   C   T  Non-Essential     whiB7  Rv3197A      274   \n",
       "1  3568410_C  3568410   A   C  Non-Essential     whiB7  Rv3197A      270   \n",
       "2  3568417_G  3568417   T   G  Non-Essential     whiB7  Rv3197A      263   \n",
       "3  3568427_C  3568427   G   C  Non-Essential     whiB7  Rv3197A      253   \n",
       "4  3568435_A  3568435   G   A  Non-Essential     whiB7  Rv3197A      245   \n",
       "\n",
       "  mut_type codon_pos array_index  num_isolates  num_sublineages   sublineages  \n",
       "0        N      A92T           0             1                1   4.3.i4.2(1)  \n",
       "1        S      A90A           1             2                1      1.1.2(2)  \n",
       "2        N      K88T           2             2                1    4.10.i1(2)  \n",
       "3        N      R85G           3             2                1  2.2.1.1.1(2)  \n",
       "4        N      P82L           4             1                1   3.1.1.i2(1)  "
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(160, 14)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [3.2] Get genotypes for indels of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Function* to get the genotypes for a specifc indel and get attributes for isolates that support indel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_genotypes_for_indel(indel_i_key):\n",
    "    \n",
    "    #0 supports Ref, 1 supports Alt, 9 indicated bad quality call\n",
    "    indel_i_genotypes = INDEL_genotypes_array[INDEL_annotation_DF[INDEL_annotation_DF.key == indel_i_key].index[0] , :] \n",
    "\n",
    "    #count the number of isolates that support this indel call\n",
    "    num_isolates_indel_i = np.sum(indel_i_genotypes == 1)\n",
    "\n",
    "    #create a boolean filter and extract information for the isolates taht support this call\n",
    "    indel_i_isolate_filter = indel_i_genotypes == 1\n",
    "    isolate_annotation_indel_i = INDEL_isolate_annotation_DF[indel_i_isolate_filter]\n",
    "\n",
    "    #find number of sub-lineages that have indel support in at least 1 isolate\n",
    "    num_sublineages_with_indel_i = len(set(list(isolate_annotation_indel_i.lineage_call.values)))\n",
    "    \n",
    "    #get list of sublineages w/ at least 1 isolates that supports indel call\n",
    "    sublineage_with_indel_count_dict = Counter(isolate_annotation_indel_i.lineage_call.values)\n",
    "    sublineages_with_indel_list = []\n",
    "    for sublineage_i in sublineage_with_indel_count_dict.keys():\n",
    "\n",
    "        sublineage_i_with_indel = sublineage_i + '({0})'.format(str(sublineage_with_indel_count_dict[sublineage_i]))\n",
    "        sublineages_with_indel_list.append(sublineage_i_with_indel)\n",
    "\n",
    "    sublineages_with_indel = ' - '.join(sublineages_with_indel_list)\n",
    "    \n",
    "    return [num_isolates_indel_i , num_sublineages_with_indel_i , sublineages_with_indel]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "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_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",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C_3568429_CCG</td>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C_3568434_CAA</td>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>82</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>GCTT_3568439_G</td>\n",
       "      <td>3568439</td>\n",
       "      <td>GCTT</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>241</td>\n",
       "      <td>del</td>\n",
       "      <td>inframe</td>\n",
       "      <td>81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C_3568487_CG</td>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GC_3568488_G</td>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              key      pos   ref  alt  gene_category gene_name  gene_id  \\\n",
       "0   C_3568429_CCG  3568429     C  CCG  Non-Essential     whiB7  Rv3197A   \n",
       "1   C_3568434_CAA  3568434     C  CAA  Non-Essential     whiB7  Rv3197A   \n",
       "2  GCTT_3568439_G  3568439  GCTT    G  Non-Essential     whiB7  Rv3197A   \n",
       "3    C_3568487_CG  3568487     C   CG  Non-Essential     whiB7  Rv3197A   \n",
       "4    GC_3568488_G  3568488    GC    G  Non-Essential     whiB7  Rv3197A   \n",
       "\n",
       "  gene_pos ins_del  INDEL_type codon_pos  \n",
       "0      251     ins  frameshift        84  \n",
       "1      246     ins  frameshift        82  \n",
       "2      241     del     inframe        81  \n",
       "3      193     ins  frameshift        65  \n",
       "4      192     del  frameshift        64  "
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26, 11)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_isolates_with_indel = []\n",
    "num_sublineages_with_indel = []\n",
    "sublineages_with_indel = []\n",
    "\n",
    "for indel_i_key in INDEL_annotation_DF.key:\n",
    "    \n",
    "    num_isolates_with_indel_i , num_sublineages_with_indel_i , sublineages_with_indel_i = get_genotypes_for_indel(indel_i_key)\n",
    "    \n",
    "    num_isolates_with_indel.append(num_isolates_with_indel_i)\n",
    "    num_sublineages_with_indel.append(num_sublineages_with_indel_i)\n",
    "    sublineages_with_indel.append(sublineages_with_indel_i)\n",
    "    \n",
    "INDEL_annotation_DF.loc[: , 'num_isolates'] = num_isolates_with_indel\n",
    "INDEL_annotation_DF.loc[: , 'num_sublineages'] = num_sublineages_with_indel\n",
    "INDEL_annotation_DF.loc[: , 'sublineages'] = sublineages_with_indel\n",
    "\n",
    "#drop indels present in 0 isolates\n",
    "INDEL_annotation_DF = INDEL_annotation_DF[INDEL_annotation_DF.num_isolates > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "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_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",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C_3568429_CCG</td>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>84</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C_3568434_CAA</td>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>82</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>GCTT_3568439_G</td>\n",
       "      <td>3568439</td>\n",
       "      <td>GCTT</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>241</td>\n",
       "      <td>del</td>\n",
       "      <td>inframe</td>\n",
       "      <td>81</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.1.i1.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C_3568487_CG</td>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>65</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GC_3568488_G</td>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>64</td>\n",
       "      <td>573</td>\n",
       "      <td>3</td>\n",
       "      <td>1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              key      pos   ref  alt  gene_category gene_name  gene_id  \\\n",
       "0   C_3568429_CCG  3568429     C  CCG  Non-Essential     whiB7  Rv3197A   \n",
       "1   C_3568434_CAA  3568434     C  CAA  Non-Essential     whiB7  Rv3197A   \n",
       "2  GCTT_3568439_G  3568439  GCTT    G  Non-Essential     whiB7  Rv3197A   \n",
       "3    C_3568487_CG  3568487     C   CG  Non-Essential     whiB7  Rv3197A   \n",
       "4    GC_3568488_G  3568488    GC    G  Non-Essential     whiB7  Rv3197A   \n",
       "\n",
       "  gene_pos ins_del  INDEL_type codon_pos  num_isolates  num_sublineages  \\\n",
       "0      251     ins  frameshift        84             1                1   \n",
       "1      246     ins  frameshift        82             2                1   \n",
       "2      241     del     inframe        81             1                1   \n",
       "3      193     ins  frameshift        65             3                3   \n",
       "4      192     del  frameshift        64           573                3   \n",
       "\n",
       "                                   sublineages  \n",
       "0                                  4.3.i3.1(1)  \n",
       "1                                 2.2.1.1.1(2)  \n",
       "2                              4.1.i1.1.1.1(1)  \n",
       "3  4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)  \n",
       "4  1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)  "
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26, 14)"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(INDEL_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Check co-occurrence between *regulator* SNPs/INDELs & *regulon* SNPs/INDELs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SNPs in regulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_SNPs_DF = SNP_annotation_DF[SNP_annotation_DF.gene_id == regulator_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "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_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>array_index</th>\n",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>3568686_C</td>\n",
       "      <td>3568686</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>56</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>3568686_G</td>\n",
       "      <td>3568686</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>56</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>3568692_T</td>\n",
       "      <td>3568692</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>57</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>3568701_A</td>\n",
       "      <td>3568701</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>58</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>63</th>\n",
       "      <td>3568702_T</td>\n",
       "      <td>3568702</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>59</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1(1)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          key      pos ref alt gene_category gene_name          gene_id  \\\n",
       "59  3568686_C  3568686   T   C          None      None  Rv3197A_Rv3198c   \n",
       "60  3568686_G  3568686   T   G          None      None  Rv3197A_Rv3198c   \n",
       "61  3568692_T  3568692   C   T          None      None  Rv3197A_Rv3198c   \n",
       "62  3568701_A  3568701   G   A          None      None  Rv3197A_Rv3198c   \n",
       "63  3568702_T  3568702   G   T          None      None  Rv3197A_Rv3198c   \n",
       "\n",
       "   gene_pos mut_type codon_pos array_index  num_isolates  num_sublineages  \\\n",
       "59     None        I      None          56             1                1   \n",
       "60     None        I      None          56             1                1   \n",
       "61     None        I      None          57             1                1   \n",
       "62     None        I      None          58             1                1   \n",
       "63     None        I      None          59             1                1   \n",
       "\n",
       "     sublineages  \n",
       "59    1.2.2.1(1)  \n",
       "60  2.2.1.1.1(1)  \n",
       "61    1.2.2.1(1)  \n",
       "62    4.10.i1(1)  \n",
       "63      3.1.1(1)  "
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulator_SNPs_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(101, 14)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(regulator_SNPs_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### INDELs in regulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_INDELs_DF = INDEL_annotation_DF[INDEL_annotation_DF.gene_id == regulator_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "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_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",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>C_3568730_CG</td>\n",
       "      <td>3568730</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>CA_3568735_C</td>\n",
       "      <td>3568735</td>\n",
       "      <td>CA</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>T_3568795_TC</td>\n",
       "      <td>3568795</td>\n",
       "      <td>T</td>\n",
       "      <td>TC</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.1.i1.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>AC_3568854_A</td>\n",
       "      <td>3568854</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>C_3568881_CGG</td>\n",
       "      <td>3568881</td>\n",
       "      <td>C</td>\n",
       "      <td>CGG</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              key      pos ref  alt gene_category gene_name          gene_id  \\\n",
       "11   C_3568730_CG  3568730   C   CG          None      None  Rv3197A_Rv3198c   \n",
       "12   CA_3568735_C  3568735  CA    C          None      None  Rv3197A_Rv3198c   \n",
       "13   T_3568795_TC  3568795   T   TC          None      None  Rv3197A_Rv3198c   \n",
       "14   AC_3568854_A  3568854  AC    A          None      None  Rv3197A_Rv3198c   \n",
       "15  C_3568881_CGG  3568881   C  CGG          None      None  Rv3197A_Rv3198c   \n",
       "\n",
       "   gene_pos ins_del  INDEL_type codon_pos  num_isolates  num_sublineages  \\\n",
       "11     None     ins  frameshift      None             1                1   \n",
       "12     None     del  frameshift      None             1                1   \n",
       "13     None     ins  frameshift      None             1                1   \n",
       "14     None     del  frameshift      None             1                1   \n",
       "15     None     ins  frameshift      None             1                1   \n",
       "\n",
       "        sublineages  \n",
       "11       4.10.i1(1)  \n",
       "12         3.1.1(1)  \n",
       "13  4.1.i1.1.1.1(1)  \n",
       "14     2.2.1.1.1(1)  \n",
       "15      4.3.i3.1(1)  "
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulator_INDELs_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 14)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(regulator_INDELs_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop **INDEL_type** column and rename **ins_del** column to fit with SNPs DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rv76/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n",
      "/home/rv76/anaconda2/lib/python2.7/site-packages/pandas/core/frame.py:2746: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  **kwargs)\n"
     ]
    }
   ],
   "source": [
    "regulator_INDELs_DF.drop(['INDEL_type'], axis=1, inplace = True)\n",
    "regulator_INDELs_DF.rename(columns={'ins_del':'mut_type'}, inplace = True) #change column names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SNPs in regulon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to include SNPs in a gene that either \n",
    "\n",
    "- cause a **premature STOP codon** (nonsense mutation)\n",
    "- **synonymous mutations at the first codon that abolish the start codon** (i.e. eis, whiB7 and mmpR which start with a **valine** as this is not possible with ahpC, mmpS5 or mmpL5, which start with a **methionine**)\n",
    "- **non-synonymous mutations at the first codon that abolish the start codon** (i.e. that change the start to an AA other than **valine** and **methionine**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulon_SNPs_DF = SNP_annotation_DF[SNP_annotation_DF.gene_id == regulon_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "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_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>array_index</th>\n",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [key, pos, ref, alt, gene_category, gene_name, gene_id, gene_pos, mut_type, codon_pos, array_index, num_isolates, num_sublineages, sublineages]\n",
       "Index: []"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulon_SNPs_DF[[(AA_change[-1]=='*') | ((len(AA_change)==3) and (AA_change[1]=='1')) for AA_change in regulon_SNPs_DF.codon_pos]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulon_SNPs_DF = regulon_SNPs_DF[[(AA_change[-1]=='*') | ((len(AA_change)==3) and (AA_change[1]=='1')) for AA_change in regulon_SNPs_DF.codon_pos]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save SNPs that occur in the first codon seperately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulon_SNPs_DF[[((len(AA_change)==3) and (AA_change[1]=='1')) for AA_change in regulon_SNPs_DF.codon_pos]].to_excel('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/START_codon_SNVs/whiB7_START_codon_SNVs.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### INDELs in regulon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### keep only *frameshift* indels for gene body"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "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_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",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C_3568429_CCG</td>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>84</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C_3568434_CAA</td>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>82</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C_3568487_CG</td>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>65</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GC_3568488_G</td>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>64</td>\n",
       "      <td>573</td>\n",
       "      <td>3</td>\n",
       "      <td>1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>GC_3568501_G</td>\n",
       "      <td>3568501</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>179</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>60</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(8)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             key      pos ref  alt  gene_category gene_name  gene_id gene_pos  \\\n",
       "0  C_3568429_CCG  3568429   C  CCG  Non-Essential     whiB7  Rv3197A      251   \n",
       "1  C_3568434_CAA  3568434   C  CAA  Non-Essential     whiB7  Rv3197A      246   \n",
       "3   C_3568487_CG  3568487   C   CG  Non-Essential     whiB7  Rv3197A      193   \n",
       "4   GC_3568488_G  3568488  GC    G  Non-Essential     whiB7  Rv3197A      192   \n",
       "5   GC_3568501_G  3568501  GC    G  Non-Essential     whiB7  Rv3197A      179   \n",
       "\n",
       "  ins_del  INDEL_type codon_pos  num_isolates  num_sublineages  \\\n",
       "0     ins  frameshift        84             1                1   \n",
       "1     ins  frameshift        82             2                1   \n",
       "3     ins  frameshift        65             3                3   \n",
       "4     del  frameshift        64           573                3   \n",
       "5     del  frameshift        60             8                1   \n",
       "\n",
       "                                   sublineages  \n",
       "0                                  4.3.i3.1(1)  \n",
       "1                                 2.2.1.1.1(2)  \n",
       "3  4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)  \n",
       "4  1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)  \n",
       "5                                   4.10.i1(8)  "
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_annotation_DF[(INDEL_annotation_DF.gene_id == regulon_id) & (INDEL_annotation_DF.INDEL_type == 'frameshift')].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulon_INDELs_DF = INDEL_annotation_DF[(INDEL_annotation_DF.gene_id == regulon_id) & (INDEL_annotation_DF.INDEL_type == 'frameshift')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "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_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",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C_3568429_CCG</td>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>84</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C_3568434_CAA</td>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>82</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C_3568487_CG</td>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>65</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GC_3568488_G</td>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>64</td>\n",
       "      <td>573</td>\n",
       "      <td>3</td>\n",
       "      <td>1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>GC_3568501_G</td>\n",
       "      <td>3568501</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>179</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>60</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(8)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             key      pos ref  alt  gene_category gene_name  gene_id gene_pos  \\\n",
       "0  C_3568429_CCG  3568429   C  CCG  Non-Essential     whiB7  Rv3197A      251   \n",
       "1  C_3568434_CAA  3568434   C  CAA  Non-Essential     whiB7  Rv3197A      246   \n",
       "3   C_3568487_CG  3568487   C   CG  Non-Essential     whiB7  Rv3197A      193   \n",
       "4   GC_3568488_G  3568488  GC    G  Non-Essential     whiB7  Rv3197A      192   \n",
       "5   GC_3568501_G  3568501  GC    G  Non-Essential     whiB7  Rv3197A      179   \n",
       "\n",
       "  ins_del  INDEL_type codon_pos  num_isolates  num_sublineages  \\\n",
       "0     ins  frameshift        84             1                1   \n",
       "1     ins  frameshift        82             2                1   \n",
       "3     ins  frameshift        65             3                3   \n",
       "4     del  frameshift        64           573                3   \n",
       "5     del  frameshift        60             8                1   \n",
       "\n",
       "                                   sublineages  \n",
       "0                                  4.3.i3.1(1)  \n",
       "1                                 2.2.1.1.1(2)  \n",
       "3  4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)  \n",
       "4  1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)  \n",
       "5                                   4.10.i1(8)  "
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulon_INDELs_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 14)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(regulon_INDELs_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop **INDEL_type** column (since all indels are frameshifts) and rename **ins_del** column to fit with SNPs DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rv76/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "regulon_INDELs_DF.drop(['INDEL_type'], axis=1, inplace = True)\n",
    "regulon_INDELs_DF.rename(columns={'ins_del':'mut_type'}, inplace = True) #change column names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save DataFrame for *regulator* & *regulon* mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_SNPs_DF_for_export = regulator_SNPs_DF.drop(['array_index'], axis=1)\n",
    "regulon_SNPs_DF_for_export = regulon_SNPs_DF.drop(['array_index'], axis=1)\n",
    "regulator_regulon_mutations_df = regulator_SNPs_DF_for_export.append(regulon_SNPs_DF_for_export.append(regulator_INDELs_DF.append(regulon_INDELs_DF)))\n",
    "\n",
    "regulator_regulon_mutations_df.drop(['key'], axis=1, inplace = True) #drop col\n",
    "regulator_regulon_mutations_df.sort_values(by = 'pos', inplace = True) #sort by reference position\n",
    "regulator_regulon_mutations_df.reset_index(inplace = True, drop = True) #reset index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(126, 12)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(regulator_regulon_mutations_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "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>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3568429</td>\n",
       "      <td>C</td>\n",
       "      <td>CCG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>251</td>\n",
       "      <td>ins</td>\n",
       "      <td>84</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3568434</td>\n",
       "      <td>C</td>\n",
       "      <td>CAA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>246</td>\n",
       "      <td>ins</td>\n",
       "      <td>82</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>65</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>64</td>\n",
       "      <td>573</td>\n",
       "      <td>3</td>\n",
       "      <td>1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3568501</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>179</td>\n",
       "      <td>del</td>\n",
       "      <td>60</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(8)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3568547</td>\n",
       "      <td>TCA</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>133</td>\n",
       "      <td>del</td>\n",
       "      <td>45</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>4.3.i4.2(1) - 2.2.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3568626</td>\n",
       "      <td>CA</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>54</td>\n",
       "      <td>del</td>\n",
       "      <td>18</td>\n",
       "      <td>61</td>\n",
       "      <td>2</td>\n",
       "      <td>2.2.1.1(1) - 2.2.1.1.1(60)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3568646</td>\n",
       "      <td>TG</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>34</td>\n",
       "      <td>del</td>\n",
       "      <td>12</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>3568646</td>\n",
       "      <td>T</td>\n",
       "      <td>TG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>34</td>\n",
       "      <td>ins</td>\n",
       "      <td>12</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(3)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3568648</td>\n",
       "      <td>G</td>\n",
       "      <td>GA</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>32</td>\n",
       "      <td>ins</td>\n",
       "      <td>11</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>4.6.2.1.2(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>3568686</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>3568686</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>3568692</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>3568701</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>3568702</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>3568703</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>4.11(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>3568704</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>3568705</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>3568708</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>3568712</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>4.1.i1.1.1.1(9)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>3568714</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>3568716</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>3568717</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>3568727</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>3568730</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>3568730</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>3568731</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i4.2(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>3568731</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.3(5)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>3568732</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(9)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>3568734</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>3568936</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>3568938</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>3568943</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.3(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>3568944</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>3568950</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.2.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>101</th>\n",
       "      <td>3568951</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>102</th>\n",
       "      <td>3568957</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>103</th>\n",
       "      <td>3568965</td>\n",
       "      <td>GCC</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.2(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>104</th>\n",
       "      <td>3568966</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>6(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>105</th>\n",
       "      <td>3568969</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>106</th>\n",
       "      <td>3568974</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>4.6.2.2(3)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107</th>\n",
       "      <td>3568978</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.1.1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108</th>\n",
       "      <td>3568985</td>\n",
       "      <td>G</td>\n",
       "      <td>GCC</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>None</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1.i1(11)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>109</th>\n",
       "      <td>3568995</td>\n",
       "      <td>G</td>\n",
       "      <td>GC</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>ins</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1.1.1.1(1) - 1.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>3569014</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>111</th>\n",
       "      <td>3569015</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1.2.2.1(3)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>112</th>\n",
       "      <td>3569022</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4.3.i3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>113</th>\n",
       "      <td>3569026</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(2)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>114</th>\n",
       "      <td>3569029</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>249</td>\n",
       "      <td>3</td>\n",
       "      <td>4.2.1.1.1.2.2.1.1.1(5) - 4.6.1.1.1.1.2(173) - ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>115</th>\n",
       "      <td>3569030</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>4.1.i1.2.1(1) - 3.1.1(4)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>116</th>\n",
       "      <td>3569032</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>4.1.i1.1.1.1(8)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117</th>\n",
       "      <td>3569035</td>\n",
       "      <td>CTGGAGG</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>del</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td>3569037</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.2.1.1.1.i3(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>119</th>\n",
       "      <td>3569040</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>4.1.i1.2.1(1) - 3.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>120</th>\n",
       "      <td>3569047</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1.2.2.1(1) - 1.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>121</th>\n",
       "      <td>3569048</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>122</th>\n",
       "      <td>3569058</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>4.10.i1(4)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>123</th>\n",
       "      <td>3569065</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>66</td>\n",
       "      <td>4</td>\n",
       "      <td>4.3.i1(49) - 3.1.1(1) - 4.2.1.1.1.2.2.1(15) - ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>124</th>\n",
       "      <td>3569069</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td>3569081</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1.1.2(2)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>126 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         pos      ref  alt  gene_category gene_name          gene_id gene_pos  \\\n",
       "0    3568429        C  CCG  Non-Essential     whiB7          Rv3197A      251   \n",
       "1    3568434        C  CAA  Non-Essential     whiB7          Rv3197A      246   \n",
       "2    3568487        C   CG  Non-Essential     whiB7          Rv3197A      193   \n",
       "3    3568488       GC    G  Non-Essential     whiB7          Rv3197A      192   \n",
       "4    3568501       GC    G  Non-Essential     whiB7          Rv3197A      179   \n",
       "5    3568547      TCA    T  Non-Essential     whiB7          Rv3197A      133   \n",
       "6    3568626       CA    C  Non-Essential     whiB7          Rv3197A       54   \n",
       "7    3568646       TG    T  Non-Essential     whiB7          Rv3197A       34   \n",
       "8    3568646        T   TG  Non-Essential     whiB7          Rv3197A       34   \n",
       "9    3568648        G   GA  Non-Essential     whiB7          Rv3197A       32   \n",
       "10   3568686        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "11   3568686        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "12   3568692        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "13   3568701        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "14   3568702        G    T           None      None  Rv3197A_Rv3198c     None   \n",
       "15   3568703        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "16   3568704        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "17   3568705        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "18   3568708        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "19   3568712        G    C           None      None  Rv3197A_Rv3198c     None   \n",
       "20   3568714        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "21   3568716        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "22   3568717        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "23   3568727        A    G           None      None  Rv3197A_Rv3198c     None   \n",
       "24   3568730        C   CG           None      None  Rv3197A_Rv3198c     None   \n",
       "25   3568730        C    G           None      None  Rv3197A_Rv3198c     None   \n",
       "26   3568731        G    C           None      None  Rv3197A_Rv3198c     None   \n",
       "27   3568731        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "28   3568732        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "29   3568734        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "..       ...      ...  ...            ...       ...              ...      ...   \n",
       "96   3568936        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "97   3568938        A    G           None      None  Rv3197A_Rv3198c     None   \n",
       "98   3568943        G    T           None      None  Rv3197A_Rv3198c     None   \n",
       "99   3568944        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "100  3568950        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "101  3568951        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "102  3568957        C    G           None      None  Rv3197A_Rv3198c     None   \n",
       "103  3568965      GCC    G           None      None  Rv3197A_Rv3198c     None   \n",
       "104  3568966        C    G           None      None  Rv3197A_Rv3198c     None   \n",
       "105  3568969        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "106  3568974        T    G           None      None  Rv3197A_Rv3198c     None   \n",
       "107  3568978        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "108  3568985        G  GCC           None      None  Rv3197A_Rv3198c     None   \n",
       "109  3568995        G   GC           None      None  Rv3197A_Rv3198c     None   \n",
       "110  3569014        A    G           None      None  Rv3197A_Rv3198c     None   \n",
       "111  3569015        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "112  3569022        A    C           None      None  Rv3197A_Rv3198c     None   \n",
       "113  3569026        G    C           None      None  Rv3197A_Rv3198c     None   \n",
       "114  3569029        T    C           None      None  Rv3197A_Rv3198c     None   \n",
       "115  3569030        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "116  3569032        G    T           None      None  Rv3197A_Rv3198c     None   \n",
       "117  3569035  CTGGAGG    C           None      None  Rv3197A_Rv3198c     None   \n",
       "118  3569037        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "119  3569040        G    C           None      None  Rv3197A_Rv3198c     None   \n",
       "120  3569047        C    T           None      None  Rv3197A_Rv3198c     None   \n",
       "121  3569048        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "122  3569058        G    A           None      None  Rv3197A_Rv3198c     None   \n",
       "123  3569065        A    G           None      None  Rv3197A_Rv3198c     None   \n",
       "124  3569069        C    A           None      None  Rv3197A_Rv3198c     None   \n",
       "125  3569081        C    A           None      None  Rv3197A_Rv3198c     None   \n",
       "\n",
       "    mut_type codon_pos  num_isolates  num_sublineages  \\\n",
       "0        ins        84             1                1   \n",
       "1        ins        82             2                1   \n",
       "2        ins        65             3                3   \n",
       "3        del        64           573                3   \n",
       "4        del        60             8                1   \n",
       "5        del        45             2                2   \n",
       "6        del        18            61                2   \n",
       "7        del        12             2                1   \n",
       "8        ins        12             3                1   \n",
       "9        ins        11             2                1   \n",
       "10         I      None             1                1   \n",
       "11         I      None             1                1   \n",
       "12         I      None             1                1   \n",
       "13         I      None             1                1   \n",
       "14         I      None             1                1   \n",
       "15         I      None             2                1   \n",
       "16         I      None             1                1   \n",
       "17         I      None             1                1   \n",
       "18         I      None             1                1   \n",
       "19         I      None             9                1   \n",
       "20         I      None             1                1   \n",
       "21         I      None             1                1   \n",
       "22         I      None             1                1   \n",
       "23         I      None             2                1   \n",
       "24       ins      None             1                1   \n",
       "25         I      None             1                1   \n",
       "26         I      None             2                1   \n",
       "27         I      None             5                1   \n",
       "28         I      None             9                1   \n",
       "29         I      None             1                1   \n",
       "..       ...       ...           ...              ...   \n",
       "96         I      None             1                1   \n",
       "97         I      None             1                1   \n",
       "98         I      None             1                1   \n",
       "99         I      None             1                1   \n",
       "100        I      None             1                1   \n",
       "101        I      None             2                1   \n",
       "102        I      None             1                1   \n",
       "103      del      None             2                1   \n",
       "104        I      None             2                1   \n",
       "105        I      None             1                1   \n",
       "106        I      None             3                1   \n",
       "107        I      None             2                1   \n",
       "108      ins      None            11                1   \n",
       "109      ins      None             2                2   \n",
       "110        I      None             1                1   \n",
       "111        I      None             3                1   \n",
       "112        I      None             1                1   \n",
       "113        I      None             2                1   \n",
       "114        I      None           249                3   \n",
       "115        I      None             5                2   \n",
       "116        I      None             8                1   \n",
       "117      del      None             1                1   \n",
       "118        I      None             1                1   \n",
       "119        I      None             2                2   \n",
       "120        I      None             2                2   \n",
       "121        I      None             1                1   \n",
       "122        I      None             4                1   \n",
       "123        I      None            66                4   \n",
       "124        I      None             1                1   \n",
       "125        I      None             2                1   \n",
       "\n",
       "                                           sublineages  \n",
       "0                                          4.3.i3.1(1)  \n",
       "1                                         2.2.1.1.1(2)  \n",
       "2          4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)  \n",
       "3          1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)  \n",
       "4                                           4.10.i1(8)  \n",
       "5                               4.3.i4.2(1) - 2.2.2(1)  \n",
       "6                           2.2.1.1(1) - 2.2.1.1.1(60)  \n",
       "7                                                 5(2)  \n",
       "8                                         2.2.1.1.1(3)  \n",
       "9                                         4.6.2.1.2(2)  \n",
       "10                                          1.2.2.1(1)  \n",
       "11                                        2.2.1.1.1(1)  \n",
       "12                                          1.2.2.1(1)  \n",
       "13                                          4.10.i1(1)  \n",
       "14                                            3.1.1(1)  \n",
       "15                                             4.11(2)  \n",
       "16                                              4.1(1)  \n",
       "17                                          1.2.2.1(1)  \n",
       "18                                          4.10.i2(1)  \n",
       "19                                     4.1.i1.1.1.1(9)  \n",
       "20                                            1.1.2(1)  \n",
       "21                                           4.3.i2(1)  \n",
       "22                                        2.2.1.1.1(1)  \n",
       "23                                            1.1.2(2)  \n",
       "24                                          4.10.i1(1)  \n",
       "25                                         4.3.i3.1(1)  \n",
       "26                                         4.3.i4.2(2)  \n",
       "27                                            1.1.3(5)  \n",
       "28                                            1.1.2(9)  \n",
       "29                                         4.3.i3.1(1)  \n",
       "..                                                 ...  \n",
       "96                                            3.2.1(1)  \n",
       "97                                         4.3.i3.1(1)  \n",
       "98                                            1.1.3(1)  \n",
       "99                                        2.2.1.1.1(1)  \n",
       "100                                         4.2.1.1(1)  \n",
       "101                                           3.1.1(2)  \n",
       "102                                        4.3.i3.1(1)  \n",
       "103                                           2.2.2(2)  \n",
       "104                                               6(2)  \n",
       "105                                         1.1.1.2(1)  \n",
       "106                                         4.6.2.2(3)  \n",
       "107                                         1.1.1.1(2)  \n",
       "108                                       3.1.1.i1(11)  \n",
       "109                          1.1.1.1(1) - 1.2.1.1.1(1)  \n",
       "110                                         1.2.2.1(1)  \n",
       "111                                         1.2.2.1(3)  \n",
       "112                                        4.3.i3.1(1)  \n",
       "113                                         4.10.i1(2)  \n",
       "114  4.2.1.1.1.2.2.1.1.1(5) - 4.6.1.1.1.1.2(173) - ...  \n",
       "115                           4.1.i1.2.1(1) - 3.1.1(4)  \n",
       "116                                    4.1.i1.1.1.1(8)  \n",
       "117                                         1.1.1.1(1)  \n",
       "118                                    2.2.1.1.1.i3(1)  \n",
       "119                           4.1.i1.2.1(1) - 3.1.1(1)  \n",
       "120                            1.2.2.1(1) - 1.1.1.1(1)  \n",
       "121                                           1.1.2(1)  \n",
       "122                                         4.10.i1(4)  \n",
       "123  4.3.i1(49) - 3.1.1(1) - 4.2.1.1.1.2.2.1(15) - ...  \n",
       "124                                           3.1.1(1)  \n",
       "125                                           1.1.2(2)  \n",
       "\n",
       "[126 rows x 12 columns]"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulator_regulon_mutations_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_regulon_mutations_df.to_excel('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/mutations_detected_within_isolates/promoter_whiB7_gene_body_mutations_in_31440_isolates.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save DataFrame for *regulator* & *regulon* mutations for *Table 1*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**REGULATOR** Subset to the 4 most commonly occurring mutations in the regulator regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_mutations_df_subset = regulator_regulon_mutations_df[regulator_regulon_mutations_df.gene_id == regulator_id].sort_values(by = 'num_isolates', ascending = False).head(n=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**REGULON** Subset to mutations that occured in at least 2 Mtb sub-lineages in the regulon regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulon_mutations_df_subset = regulator_regulon_mutations_df[regulator_regulon_mutations_df.gene_id == regulon_id]\n",
    "regulon_mutations_df_subset = regulon_mutations_df_subset[regulon_mutations_df_subset.num_sublineages >= 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Export as excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_regulon_subset_mutations_df = regulator_mutations_df_subset.append(regulon_mutations_df_subset)\n",
    "regulator_regulon_subset_mutations_df.sort_values(by = 'pos', ascending = True, inplace = True)\n",
    "regulator_regulon_subset_mutations_df.reset_index(inplace = True, drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "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>mut_type</th>\n",
       "      <th>codon_pos</th>\n",
       "      <th>num_isolates</th>\n",
       "      <th>num_sublineages</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3568487</td>\n",
       "      <td>C</td>\n",
       "      <td>CG</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>193</td>\n",
       "      <td>ins</td>\n",
       "      <td>65</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3568488</td>\n",
       "      <td>GC</td>\n",
       "      <td>G</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>192</td>\n",
       "      <td>del</td>\n",
       "      <td>64</td>\n",
       "      <td>573</td>\n",
       "      <td>3</td>\n",
       "      <td>1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3568547</td>\n",
       "      <td>TCA</td>\n",
       "      <td>T</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>133</td>\n",
       "      <td>del</td>\n",
       "      <td>45</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>4.3.i4.2(1) - 2.2.2(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3568626</td>\n",
       "      <td>CA</td>\n",
       "      <td>C</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>54</td>\n",
       "      <td>del</td>\n",
       "      <td>18</td>\n",
       "      <td>61</td>\n",
       "      <td>2</td>\n",
       "      <td>2.2.1.1(1) - 2.2.1.1.1(60)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3568779</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>256</td>\n",
       "      <td>2</td>\n",
       "      <td>4.4.1.1(244) - 4.2.1.2.2.1(12)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3568857</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>73</td>\n",
       "      <td>2</td>\n",
       "      <td>3.1.1(72) - 3.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3568921</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>117</td>\n",
       "      <td>2</td>\n",
       "      <td>4.1.i1.1.2(110) - 4.1.i1.1(7)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3569029</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>249</td>\n",
       "      <td>3</td>\n",
       "      <td>4.2.1.1.1.2.2.1.1.1(5) - 4.6.1.1.1.1.2(173) - ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pos  ref alt  gene_category gene_name          gene_id gene_pos  \\\n",
       "0  3568487    C  CG  Non-Essential     whiB7          Rv3197A      193   \n",
       "1  3568488   GC   G  Non-Essential     whiB7          Rv3197A      192   \n",
       "2  3568547  TCA   T  Non-Essential     whiB7          Rv3197A      133   \n",
       "3  3568626   CA   C  Non-Essential     whiB7          Rv3197A       54   \n",
       "4  3568779    A   G           None      None  Rv3197A_Rv3198c     None   \n",
       "5  3568857    G   A           None      None  Rv3197A_Rv3198c     None   \n",
       "6  3568921    C   G           None      None  Rv3197A_Rv3198c     None   \n",
       "7  3569029    T   C           None      None  Rv3197A_Rv3198c     None   \n",
       "\n",
       "  mut_type codon_pos  num_isolates  num_sublineages  \\\n",
       "0      ins        65             3                3   \n",
       "1      del        64           573                3   \n",
       "2      del        45             2                2   \n",
       "3      del        18            61                2   \n",
       "4        I      None           256                2   \n",
       "5        I      None            73                2   \n",
       "6        I      None           117                2   \n",
       "7        I      None           249                3   \n",
       "\n",
       "                                         sublineages  \n",
       "0        4.10.i1(1) - 4.1.i1.1.1.1(1) - 2.2.1.1.2(1)  \n",
       "1        1.2.1.1.2(74) - 1.2.1.1(5) - 1.2.1.1.1(494)  \n",
       "2                             4.3.i4.2(1) - 2.2.2(1)  \n",
       "3                         2.2.1.1(1) - 2.2.1.1.1(60)  \n",
       "4                     4.4.1.1(244) - 4.2.1.2.2.1(12)  \n",
       "5                                 3.1.1(72) - 3.1(1)  \n",
       "6                      4.1.i1.1.2(110) - 4.1.i1.1(7)  \n",
       "7  4.2.1.1.1.2.2.1.1.1(5) - 4.6.1.1.1.1.2(173) - ...  "
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulator_regulon_subset_mutations_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulator_regulon_subset_mutations_df.to_excel('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/mutations_detected_within_isolates/promoter_whiB7_gene_body_mutations_in_31440_isolates_for_Table_1.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct DataFrame (which displays each pair of co-occuring mutations as a seperate row with info for each pair of mutations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "variant_mutation_pair_df = pd.DataFrame(columns=['A_type', 'A_pos', 'A_ref', 'A_alt', 'A_mut_type', 'A_codon_pos', 'A_gene_id', 'A_gene_name', 'A_gene_pos' ,'A_num_isolates', 'B_type', 'B_pos', 'B_ref', 'B_alt', 'B_mut_type', 'B_codon_pos', 'B_gene_id', 'B_gene_name', 'B_gene_pos', 'B_num_isolates', 'num_isolates_co_occur', 'isolate_IDs', 'sublineages'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the ordering of the INDELs genotypes columns in terms of the SNPs genotypes columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNP_isolate_annotation_DF['isolate_order'] = SNP_isolate_annotation_DF.index\n",
    "SNP_isolate_annotation_DF.set_index('isolate_ID', drop = True, inplace = True)\n",
    "INDEL_isolate_ordering_for_SNP_genotypes_array = SNP_isolate_annotation_DF.loc[INDEL_isolate_annotation_DF.isolate_ID, :].isolate_order.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  188,  2733, 24409, ..., 31437, 31438, 31439])"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "INDEL_isolate_ordering_for_SNP_genotypes_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill in dataframe by quering the genotypes matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SNPs or INDELs in regulator x SNPs or INDELs in regulon\n",
    "regulator_mut_keys = list(regulator_SNPs_DF.key) + list(regulator_INDELs_DF.key)\n",
    "regulon_mut_keys = list(regulon_SNPs_DF.key) + list(regulon_INDELs_DF.key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#index for DataFrame\n",
    "variant_mutation_pair_df_index = 0\n",
    "\n",
    "for variant_i in regulator_mut_keys: #iterate through regulator SNPs & INDELs\n",
    "    \n",
    "    if variant_i.count('_') == 1: # if variant in regulator is a SNP\n",
    "\n",
    "        SNP_i_pos = int(variant_i.split('_')[0])\n",
    "        SNP_i_alt_allele = variant_i.split('_')[1]\n",
    "\n",
    "        #get the number code that corresponds to the alt allele, 9 indicated bad quality call\n",
    "        base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "        SNP_i_alt_allele_code = base_code_dict[SNP_i_alt_allele]\n",
    "\n",
    "        #genotypes for each SNP site\n",
    "        genotypes_for_SNP_i = SNP_genotypes_array[SNP_annotation_DF[SNP_annotation_DF.pos == SNP_i_pos].array_index.values[0] , :]\n",
    "\n",
    "        #re-order isolates according to the order of the columns in the INDEL genotypes matrix\n",
    "        genotypes_for_SNP_i = genotypes_for_SNP_i[INDEL_isolate_ordering_for_SNP_genotypes_array]\n",
    "\n",
    "        for variant_j in regulon_mut_keys: #iterate through mutations in regulon\n",
    "\n",
    "            if variant_j.count('_') == 1: #if variant in regulon is a SNP\n",
    "\n",
    "                SNP_j_pos = int(variant_j.split('_')[0])\n",
    "                SNP_j_alt_allele = variant_j.split('_')[1]\n",
    "\n",
    "                #get the number code that corresponds to the alt allele, 9 indicated bad quality call\n",
    "                SNP_j_alt_allele_code = base_code_dict[SNP_j_alt_allele]\n",
    "\n",
    "                #genotypes for each SNP site\n",
    "                genotypes_for_SNP_j = SNP_genotypes_array[SNP_annotation_DF[SNP_annotation_DF.pos == SNP_j_pos].array_index.values[0] , :]\n",
    "\n",
    "                #re-order isolates according to the order of the columns in the INDEL genotypes matrix\n",
    "                genotypes_for_SNP_j = genotypes_for_SNP_j[INDEL_isolate_ordering_for_SNP_genotypes_array]\n",
    "\n",
    "                #boolean filter of isolates with both variants\n",
    "                isolates_with_both_variants_filter = [(SNP_i_called and SNP_j_called) for SNP_i_called , SNP_j_called in zip(list(genotypes_for_SNP_i == SNP_i_alt_allele_code) , list(genotypes_for_SNP_j == SNP_j_alt_allele_code))]\n",
    "\n",
    "                #number of isolates that have both variants\n",
    "                num_isolates_with_both_variants = np.sum(isolates_with_both_variants_filter)\n",
    "\n",
    "                #check to see if the number of isolates carrying both variants was >= 1\n",
    "                if num_isolates_with_both_variants >= 1:\n",
    "\n",
    "                    # get info for mutation in regulator\n",
    "                    regulator_SNP_info = list(regulator_SNPs_DF[(regulator_SNPs_DF.pos == SNP_i_pos) & (regulator_SNPs_DF.alt == SNP_i_alt_allele)].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    # get info for mutation in regulon\n",
    "                    regulon_SNP_info = list(regulon_SNPs_DF[(regulon_SNPs_DF.pos == SNP_j_pos) & (regulon_SNPs_DF.alt == SNP_j_alt_allele)].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    strain_names = [' - '.join(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].isolate_ID))]\n",
    "\n",
    "                    # sublineage calls\n",
    "                    strain_lineage_calls_count_dict = Counter(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].lineage_call))\n",
    "                    strain_lineage_calls_list = []\n",
    "                    for sublineage_i in strain_lineage_calls_count_dict.keys():\n",
    "                        sublineage_i_with_mutations = sublineage_i + '({0})'.format(str(strain_lineage_calls_count_dict[sublineage_i]))\n",
    "                        strain_lineage_calls_list.append(sublineage_i_with_mutations)\n",
    "                    strain_lineage_calls = ' - '.join(strain_lineage_calls_list)\n",
    "\n",
    "                    co_occurring_SNP_i_and_SNP_or_indel_j_info = ['SNP'] + regulator_SNP_info + ['SNP'] + regulon_SNP_info + [num_isolates_with_both_variants] + strain_names + [strain_lineage_calls]\n",
    "                    variant_mutation_pair_df.loc[variant_mutation_pair_df_index, :] = co_occurring_SNP_i_and_SNP_or_indel_j_info\n",
    "                    variant_mutation_pair_df_index += 1\n",
    "\n",
    "            if variant_j.count('_') == 2: #if variant in regulon is an INDEL\n",
    "\n",
    "                #genotypes for each INDEL\n",
    "                genotypes_for_indel_j = INDEL_genotypes_array[INDEL_annotation_DF[INDEL_annotation_DF.key == variant_j].index[0] , :]\n",
    "\n",
    "                #boolean filter of isolates with both variants\n",
    "                isolates_with_both_variants_filter = [(SNP_i_called and indel_j_called) for SNP_i_called , indel_j_called in zip(list(genotypes_for_SNP_i == SNP_i_alt_allele_code) , list(genotypes_for_indel_j == 1))]\n",
    "\n",
    "                #number of isolates that have both variants\n",
    "                num_isolates_with_both_variants = np.sum(isolates_with_both_variants_filter)\n",
    "\n",
    "                #check to see if the number of isolates carrying both variants was >= 1\n",
    "                if num_isolates_with_both_variants >= 1:\n",
    "\n",
    "                    # get info for mutation in regulator\n",
    "                    regulator_SNP_info = list(regulator_SNPs_DF[(regulator_SNPs_DF.pos == SNP_i_pos) & (regulator_SNPs_DF.alt == SNP_i_alt_allele)].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    # get info for mutation in regulon\n",
    "                    regulon_INDEL_info = list(regulon_INDELs_DF[regulon_INDELs_DF.key == variant_j].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    strain_names = [' - '.join(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].isolate_ID))]\n",
    "\n",
    "                    # sublineage calls\n",
    "                    strain_lineage_calls_count_dict = Counter(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].lineage_call))\n",
    "                    strain_lineage_calls_list = []\n",
    "                    for sublineage_i in strain_lineage_calls_count_dict.keys():\n",
    "                        sublineage_i_with_mutations = sublineage_i + '({0})'.format(str(strain_lineage_calls_count_dict[sublineage_i]))\n",
    "                        strain_lineage_calls_list.append(sublineage_i_with_mutations)\n",
    "                    strain_lineage_calls = ' - '.join(strain_lineage_calls_list)\n",
    "\n",
    "                    co_occurring_SNP_i_and_SNP_or_indel_j_info = ['SNP'] + regulator_SNP_info + ['INDEL'] + regulon_INDEL_info + [num_isolates_with_both_variants] + strain_names + [strain_lineage_calls]\n",
    "                    variant_mutation_pair_df.loc[variant_mutation_pair_df_index, :] = co_occurring_SNP_i_and_SNP_or_indel_j_info\n",
    "                    variant_mutation_pair_df_index += 1\n",
    "    \n",
    "    elif variant_i.count('_') == 2: # if variant in regulator is an INDEL\n",
    "        \n",
    "        #genotypes for each INDEL\n",
    "        genotypes_for_indel_i = INDEL_genotypes_array[INDEL_annotation_DF[INDEL_annotation_DF.key == variant_i].index[0] , :]\n",
    "\n",
    "        for variant_j in regulon_mut_keys: #iterate through mutations in regulon\n",
    "\n",
    "            if variant_j.count('_') == 1: #if variant in regulon is a SNP\n",
    "\n",
    "                SNP_j_pos = int(variant_j.split('_')[0])\n",
    "                SNP_j_alt_allele = variant_j.split('_')[1]\n",
    "\n",
    "                #get the number code that corresponds to the alt allele, 9 indicated bad quality call\n",
    "                SNP_j_alt_allele_code = base_code_dict[SNP_j_alt_allele]\n",
    "\n",
    "                #genotypes for each SNP site\n",
    "                genotypes_for_SNP_j = SNP_genotypes_array[SNP_annotation_DF[SNP_annotation_DF.pos == SNP_j_pos].array_index.values[0] , :]\n",
    "\n",
    "                #re-order isolates according to the order of the columns in the INDEL genotypes matrix\n",
    "                genotypes_for_SNP_j = genotypes_for_SNP_j[INDEL_isolate_ordering_for_SNP_genotypes_array]\n",
    "\n",
    "                #boolean filter of isolates with both variants\n",
    "                isolates_with_both_variants_filter = [(indel_i_called and SNP_j_called) for indel_i_called , SNP_j_called in zip(list(genotypes_for_indel_i == 1) , list(genotypes_for_SNP_j == SNP_j_alt_allele_code))]\n",
    "\n",
    "                #number of isolates that have both variants\n",
    "                num_isolates_with_both_variants = np.sum(isolates_with_both_variants_filter)\n",
    "\n",
    "                #check to see if the number of isolates carrying both variants was >= 1\n",
    "                if num_isolates_with_both_variants >= 1:\n",
    "\n",
    "                    # get info for mutation in regulator\n",
    "                    regulator_INDEL_info = list(regulator_INDELs_DF[regulator_INDELs_DF.key == variant_i].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    # get info for mutation in regulon\n",
    "                    regulon_SNP_info = list(regulon_SNPs_DF[(regulon_SNPs_DF.pos == SNP_j_pos) & (regulon_SNPs_DF.alt == SNP_j_alt_allele)].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    strain_names = [' - '.join(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].isolate_ID))]\n",
    "\n",
    "                    # sublineage calls\n",
    "                    strain_lineage_calls_count_dict = Counter(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].lineage_call))\n",
    "                    strain_lineage_calls_list = []\n",
    "                    for sublineage_i in strain_lineage_calls_count_dict.keys():\n",
    "                        sublineage_i_with_mutations = sublineage_i + '({0})'.format(str(strain_lineage_calls_count_dict[sublineage_i]))\n",
    "                        strain_lineage_calls_list.append(sublineage_i_with_mutations)\n",
    "                    strain_lineage_calls = ' - '.join(strain_lineage_calls_list)\n",
    "\n",
    "                    co_occurring_SNP_i_and_SNP_or_indel_j_info = ['INDEL'] + regulator_INDEL_info + ['SNP'] + regulon_SNP_info + [num_isolates_with_both_variants] + strain_names + [strain_lineage_calls]\n",
    "                    variant_mutation_pair_df.loc[variant_mutation_pair_df_index, :] = co_occurring_SNP_i_and_SNP_or_indel_j_info\n",
    "                    variant_mutation_pair_df_index += 1\n",
    "\n",
    "            if variant_j.count('_') == 2: # if variant in regulon is an INDEL\n",
    "\n",
    "                #genotypes for each INDEL\n",
    "                genotypes_for_indel_j = INDEL_genotypes_array[INDEL_annotation_DF[INDEL_annotation_DF.key == variant_j].index[0] , :]\n",
    "\n",
    "                #boolean filter of isolates with both variants\n",
    "                isolates_with_both_variants_filter = [(indel_i_called and indel_j_called) for indel_i_called , indel_j_called in zip(list(genotypes_for_indel_i == 1) , list(genotypes_for_indel_j == 1))]\n",
    "\n",
    "                #number of isolates that have both variants\n",
    "                num_isolates_with_both_variants = np.sum(isolates_with_both_variants_filter)\n",
    "\n",
    "                #check to see if the number of isolates carrying both variants was >= 1\n",
    "                if num_isolates_with_both_variants >= 1:\n",
    "\n",
    "                    # get info for mutation in regulator\n",
    "                    regulator_INDEL_info = list(regulator_INDELs_DF[regulator_INDELs_DF.key == variant_i].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    # get info for mutation in regulon\n",
    "                    regulon_INDEL_info = list(regulon_INDELs_DF[regulon_INDELs_DF.key == variant_j].loc[:, ['pos','ref','alt','mut_type','codon_pos','gene_id','gene_name','gene_pos','num_isolates']].values[0])\n",
    "\n",
    "                    strain_names = [' - '.join(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].isolate_ID))]\n",
    "\n",
    "                    # sublineage calls\n",
    "                    strain_lineage_calls_count_dict = Counter(list(INDEL_isolate_annotation_DF[isolates_with_both_variants_filter].lineage_call))\n",
    "                    strain_lineage_calls_list = []\n",
    "                    for sublineage_i in strain_lineage_calls_count_dict.keys():\n",
    "                        sublineage_i_with_mutations = sublineage_i + '({0})'.format(str(strain_lineage_calls_count_dict[sublineage_i]))\n",
    "                        strain_lineage_calls_list.append(sublineage_i_with_mutations)\n",
    "                    strain_lineage_calls = ' - '.join(strain_lineage_calls_list)\n",
    "\n",
    "                    co_occurring_SNP_i_and_SNP_or_indel_j_info = ['INDEL'] + regulator_INDEL_info + ['INDEL'] + regulon_INDEL_info + [num_isolates_with_both_variants] + strain_names + [strain_lineage_calls]\n",
    "                    variant_mutation_pair_df.loc[variant_mutation_pair_df_index, :] = co_occurring_SNP_i_and_SNP_or_indel_j_info\n",
    "                    variant_mutation_pair_df_index += 1"
   ]
  },
  {
   "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>A_type</th>\n",
       "      <th>A_pos</th>\n",
       "      <th>A_ref</th>\n",
       "      <th>A_alt</th>\n",
       "      <th>A_mut_type</th>\n",
       "      <th>A_codon_pos</th>\n",
       "      <th>A_gene_id</th>\n",
       "      <th>A_gene_name</th>\n",
       "      <th>A_gene_pos</th>\n",
       "      <th>A_num_isolates</th>\n",
       "      <th>...</th>\n",
       "      <th>B_alt</th>\n",
       "      <th>B_mut_type</th>\n",
       "      <th>B_codon_pos</th>\n",
       "      <th>B_gene_id</th>\n",
       "      <th>B_gene_name</th>\n",
       "      <th>B_gene_pos</th>\n",
       "      <th>B_num_isolates</th>\n",
       "      <th>num_isolates_co_occur</th>\n",
       "      <th>isolate_IDs</th>\n",
       "      <th>sublineages</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SNP</td>\n",
       "      <td>3568826</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>G</td>\n",
       "      <td>del</td>\n",
       "      <td>64</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>192</td>\n",
       "      <td>573</td>\n",
       "      <td>1</td>\n",
       "      <td>SAMN07766333</td>\n",
       "      <td>1.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>INDEL</td>\n",
       "      <td>3568893</td>\n",
       "      <td>TG</td>\n",
       "      <td>T</td>\n",
       "      <td>del</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>9</td>\n",
       "      <td>...</td>\n",
       "      <td>G</td>\n",
       "      <td>del</td>\n",
       "      <td>64</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>192</td>\n",
       "      <td>573</td>\n",
       "      <td>1</td>\n",
       "      <td>SAMN08436033</td>\n",
       "      <td>1.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>INDEL</td>\n",
       "      <td>3568995</td>\n",
       "      <td>G</td>\n",
       "      <td>GC</td>\n",
       "      <td>ins</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv3197A_Rv3198c</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>G</td>\n",
       "      <td>del</td>\n",
       "      <td>64</td>\n",
       "      <td>Rv3197A</td>\n",
       "      <td>whiB7</td>\n",
       "      <td>192</td>\n",
       "      <td>573</td>\n",
       "      <td>1</td>\n",
       "      <td>SAMN07766299</td>\n",
       "      <td>1.2.1.1.1(1)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 23 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  A_type    A_pos A_ref A_alt A_mut_type A_codon_pos        A_gene_id  \\\n",
       "0    SNP  3568826     A     G          I        None  Rv3197A_Rv3198c   \n",
       "1  INDEL  3568893    TG     T        del        None  Rv3197A_Rv3198c   \n",
       "2  INDEL  3568995     G    GC        ins        None  Rv3197A_Rv3198c   \n",
       "\n",
       "  A_gene_name A_gene_pos A_num_isolates      ...      B_alt B_mut_type  \\\n",
       "0        None       None              1      ...          G        del   \n",
       "1        None       None              9      ...          G        del   \n",
       "2        None       None              2      ...          G        del   \n",
       "\n",
       "  B_codon_pos B_gene_id B_gene_name B_gene_pos B_num_isolates  \\\n",
       "0          64   Rv3197A       whiB7        192            573   \n",
       "1          64   Rv3197A       whiB7        192            573   \n",
       "2          64   Rv3197A       whiB7        192            573   \n",
       "\n",
       "  num_isolates_co_occur   isolate_IDs   sublineages  \n",
       "0                     1  SAMN07766333  1.2.1.1.1(1)  \n",
       "1                     1  SAMN08436033  1.2.1.1.1(1)  \n",
       "2                     1  SAMN07766299  1.2.1.1.1(1)  \n",
       "\n",
       "[3 rows x 23 columns]"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variant_mutation_pair_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 23)"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(variant_mutation_pair_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save dataframe as CSV file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "variant_mutation_pair_df.to_excel('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/co_occurring_mutations_within_isolates/whiB7_promoter_gene_body_co_occurence_mutations_in_31440_isolates.xlsx')"
   ]
  }
 ],
 "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
}