(H) Check for atpE, rrs mutations in specific strains and analyze co-occurrence of eis promoter mutations with rrs AG resistance mutations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook was created to 1. check for mutations in *atpE* in *mmpR* / *mmpL5* double mutants, 2. check for mutations in *rrs* in *eis* promoter / *eis* double mutants and 3. check for co-occurence of 4 most common *eis* promoter mutations with 3 *rrs* aminoglycoside resistance conferring mutations"
   ]
  },
  {
   "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": 3,
   "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": [
    "### *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": [
    "### *Function* to retrieve isolates with INDELs in a specific genomic loci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_strains_with_indels_in_gene(gene_intergenic_id, strains_of_interest):\n",
    "\n",
    "    # filter STRAINS\n",
    "    strains_to_keep = [strain_i in strains_of_interest for strain_i in INDEL_isolate_annotation_DF.isolate_ID] # construct boolean filter for strains to keep\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    INDEL_genotypes_array_filtered = INDEL_genotypes_array[: , strains_to_keep]\n",
    "\n",
    "    # filter Isolate annotation file\n",
    "    INDEL_isolate_annotation_DF_filtered = INDEL_isolate_annotation_DF[strains_to_keep]\n",
    "    INDEL_isolate_annotation_DF_filtered.reset_index(drop = True , inplace = True) #re-index new filtered INDEL annotation DF (so new index matches indexing of genotypes matrix rows)\n",
    "\n",
    "    # filter VARIANTS\n",
    "    indels_in_targeted_region = list(INDEL_annotation_DF[INDEL_annotation_DF.gene_id == gene_intergenic_id].key) # keep only indels that occur within this gene/intergenic region\n",
    "    indels_to_keep_filter = [indel_i in indels_in_targeted_region for indel_i in INDEL_annotation_DF.key] # construct boolean filter\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    INDEL_genotypes_array_filtered = INDEL_genotypes_array_filtered[indels_to_keep_filter , :]\n",
    "\n",
    "    # filter INDEL annotation file\n",
    "    INDEL_annotation_DF_filtered = INDEL_annotation_DF[indels_to_keep_filter]\n",
    "    INDEL_annotation_DF_filtered.reset_index(drop = True , inplace = True) #re-index new filtered INDEL annotation DF (so new index matches indexing of genotypes matrix rows)\n",
    "\n",
    "    isolate_with_variant_df = pd.DataFrame(columns = ['key','pos','ref','alt','gene_category','gene_name','gene_id','gene_pos','ins_del','INDEL_type','codon_pos','isolate_ID','lineage_call','group'])\n",
    "    isolate_with_variant_df_index = 0\n",
    "    # iterate through each variant (row) and store a row into a dataframe for each isolate with an INDEL\n",
    "    for row_i in INDEL_annotation_DF_filtered.index:\n",
    "\n",
    "        #iterate through each strain for this variant row\n",
    "        for strain_j in INDEL_isolate_annotation_DF_filtered.index:\n",
    "\n",
    "            genotypes_row_i_strain_j = INDEL_genotypes_array_filtered[row_i, strain_j]\n",
    "\n",
    "            #check if strain_i has variant\n",
    "            if genotypes_row_i_strain_j == 1:\n",
    "\n",
    "                #store a row indicating that this strain had this variant\n",
    "                isolate_j_with_variant_i_row = list(INDEL_annotation_DF_filtered.loc[row_i]) + list(INDEL_isolate_annotation_DF_filtered.loc[strain_j, ['isolate_ID','lineage_call','group']])\n",
    "                isolate_with_variant_df.loc[isolate_with_variant_df_index, :] = isolate_j_with_variant_i_row\n",
    "                isolate_with_variant_df_index += 1\n",
    "            \n",
    "    return isolate_with_variant_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AC_777875_A</td>\n",
       "      <td>777875</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>mmpL5</td>\n",
       "      <td>Rv0676c</td>\n",
       "      <td>606</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>202</td>\n",
       "      <td>99R893</td>\n",
       "      <td>4.11</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AC_777875_A</td>\n",
       "      <td>777875</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>mmpL5</td>\n",
       "      <td>Rv0676c</td>\n",
       "      <td>606</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>202</td>\n",
       "      <td>Peru4670</td>\n",
       "      <td>4.11</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AC_777875_A</td>\n",
       "      <td>777875</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>mmpL5</td>\n",
       "      <td>Rv0676c</td>\n",
       "      <td>606</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>202</td>\n",
       "      <td>SAMEA5569914</td>\n",
       "      <td>4.11</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AC_777875_A</td>\n",
       "      <td>777875</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>mmpL5</td>\n",
       "      <td>Rv0676c</td>\n",
       "      <td>606</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>202</td>\n",
       "      <td>Peru5003</td>\n",
       "      <td>4.11</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AC_777875_A</td>\n",
       "      <td>777875</td>\n",
       "      <td>AC</td>\n",
       "      <td>A</td>\n",
       "      <td>Non-Essential</td>\n",
       "      <td>mmpL5</td>\n",
       "      <td>Rv0676c</td>\n",
       "      <td>606</td>\n",
       "      <td>del</td>\n",
       "      <td>frameshift</td>\n",
       "      <td>202</td>\n",
       "      <td>Peru3367</td>\n",
       "      <td>4.11</td>\n",
       "      <td>NaN</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  AC_777875_A  777875  AC   A  Non-Essential     mmpL5  Rv0676c      606   \n",
       "1  AC_777875_A  777875  AC   A  Non-Essential     mmpL5  Rv0676c      606   \n",
       "2  AC_777875_A  777875  AC   A  Non-Essential     mmpL5  Rv0676c      606   \n",
       "3  AC_777875_A  777875  AC   A  Non-Essential     mmpL5  Rv0676c      606   \n",
       "4  AC_777875_A  777875  AC   A  Non-Essential     mmpL5  Rv0676c      606   \n",
       "\n",
       "  ins_del  INDEL_type codon_pos    isolate_ID lineage_call group  \n",
       "0     del  frameshift       202        99R893         4.11   NaN  \n",
       "1     del  frameshift       202      Peru4670         4.11   NaN  \n",
       "2     del  frameshift       202  SAMEA5569914         4.11   NaN  \n",
       "3     del  frameshift       202      Peru5003         4.11   NaN  \n",
       "4     del  frameshift       202      Peru3367         4.11   NaN  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_strains_with_indels_in_gene('Rv0676c', ['99R893','Peru4670','SAMEA5569914','Peru5003','Peru3367'])"
   ]
  },
  {
   "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": 22,
   "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 SNP 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_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[:, 'mut_key'] = SNP_func_annotation_DF.index\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": 23,
   "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": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_isolate_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 2)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>48</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>47.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>64</td>\n",
       "      <td>G</td>\n",
       "      <td>[C]</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>63.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>67</td>\n",
       "      <td>G</td>\n",
       "      <td>[A, T]</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>66.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>69</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>68.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>71</td>\n",
       "      <td>C</td>\n",
       "      <td>[T]</td>\n",
       "      <td>Rv0001</td>\n",
       "      <td>70.0</td>\n",
       "      <td>Essential</td>\n",
       "      <td>dnaA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pos ref     alt gene_id  gene_pos gene_category gene_name\n",
       "0   48   C     [T]  Rv0001      47.0     Essential      dnaA\n",
       "1   64   G     [C]  Rv0001      63.0     Essential      dnaA\n",
       "2   67   G  [A, T]  Rv0001      66.0     Essential      dnaA\n",
       "3   69   C     [T]  Rv0001      68.0     Essential      dnaA\n",
       "4   71   C     [T]  Rv0001      70.0     Essential      dnaA"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 7)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "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>mut_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",
       "     mut_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": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_func_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(844429, 10)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_func_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_genotypes_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 31428)"
      ]
     },
     "execution_count": 30,
     "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": 31,
   "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": 32,
   "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": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_DF_extra_strains.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "genotypes_array_extra_strains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 12)"
      ]
     },
     "execution_count": 34,
     "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": 35,
   "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": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31440, 2)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge genotypes matrices\n",
    "SNP_genotypes_array = np.hstack((SNP_genotypes_array, genotypes_array_extra_strains))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(782565, 31440)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_genotypes_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3] Look for resistance mutations for specific sets of strains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Function* to retrieve isolates with SNVs in a specific genomic loci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_strains_with_SNPs_in_gene(gene_intergenic_id, strains_of_interest):\n",
    "\n",
    "    #get the number code that corresponds to the minor allele, 9 indicated bad quality call\n",
    "    base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "\n",
    "    # filter STRAINS\n",
    "    strains_to_keep = [strain_i in strains_of_interest for strain_i in SNP_isolate_annotation_DF.isolate_ID] # construct boolean filter for strains to keep\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    SNP_genotypes_array_filtered = SNP_genotypes_array[: , strains_to_keep]\n",
    "\n",
    "    # filter Isolate annotation file\n",
    "    SNP_isolate_annotation_DF_filtered = SNP_isolate_annotation_DF[strains_to_keep]\n",
    "    SNP_isolate_annotation_DF_filtered.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)\n",
    "\n",
    "    # filter VARIANTS\n",
    "    SNPs_in_targeted_region = list(set(SNP_func_annotation_DF[SNP_func_annotation_DF.gene_id == gene_intergenic_id].pos)) # keep only SNPs that occur within this gene/intergenic region\n",
    "    SNPs_to_keep_filter = [SNP_i in SNPs_in_targeted_region for SNP_i in SNP_annotation_DF.pos] # construct boolean filter\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    SNP_genotypes_array_filtered = SNP_genotypes_array_filtered[SNPs_to_keep_filter , :]\n",
    "\n",
    "    # filter SNP annotation file\n",
    "    SNP_annotation_DF_filtered = SNP_annotation_DF[SNPs_to_keep_filter]\n",
    "    SNP_annotation_DF_filtered.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)\n",
    "\n",
    "    isolate_with_variant_df = pd.DataFrame(columns = ['mut_key','pos','ref','alt','isolate_ID','lineage_call'])\n",
    "    isolate_with_variant_df_index = 0\n",
    "    # iterate through each variant (row) and store a row into a dataframe for each isolate with an INDEL\n",
    "    for row_i in SNP_annotation_DF_filtered.index:\n",
    "\n",
    "        # get a list of alternate alleles for this SNP site (row\n",
    "        row_i_alt_allele_list = SNP_annotation_DF_filtered.loc[row_i, 'alt']\n",
    "\n",
    "        # iterate through each strain for this variant row\n",
    "        for strain_j in SNP_isolate_annotation_DF_filtered.index:\n",
    "\n",
    "            genotypes_row_i_strain_j = SNP_genotypes_array_filtered[row_i, strain_j]\n",
    "\n",
    "            # check if strain_i has variant for each alternate allele\n",
    "            for alt_allele_i in row_i_alt_allele_list:\n",
    "\n",
    "                if genotypes_row_i_strain_j == base_code_dict[alt_allele_i]:\n",
    "\n",
    "                    # store a row indicating that this strain had this variant\n",
    "                    isolate_j_with_variant_i_row = [str(SNP_annotation_DF_filtered.loc[row_i, 'pos']) + '_' + alt_allele_i] + list(SNP_annotation_DF_filtered.loc[row_i, ['pos','ref']]) + [alt_allele_i] + list(SNP_isolate_annotation_DF_filtered.loc[strain_j, ['isolate_ID','lineage_call']])\n",
    "                    isolate_with_variant_df.loc[isolate_with_variant_df_index, :] = isolate_j_with_variant_i_row\n",
    "                    isolate_with_variant_df_index += 1\n",
    "                    \n",
    "    # get functional annotation for df\n",
    "    isolate_with_variant_df = isolate_with_variant_df.merge(SNP_func_annotation_DF.loc[:,['gene_category','gene_name','gene_id','gene_pos','SNP_type','AA_change','mut_key']], how='left', on='mut_key')\n",
    "            \n",
    "    return isolate_with_variant_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Function* to retrieve genotypes for a SNP site given a set of strains, gene and gene coordinate for SNP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_genotypes_for_SNP(strains_of_interest, gene_id, gene_pos_SNP_i):\n",
    "\n",
    "    # get the number code that corresponds to the minor allele, 9 indicated bad quality call\n",
    "    base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "    base_code_dict_r = {0:'A', 1:'C', 2:'G', 3:'T', 9:'N'}\n",
    "\n",
    "    # filter STRAINS\n",
    "    strains_to_keep = [strain_i in strains_of_interest for strain_i in SNP_isolate_annotation_DF.isolate_ID] # construct boolean filter for strains to keep\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    SNP_genotypes_array_filtered = SNP_genotypes_array[: , strains_to_keep]\n",
    "\n",
    "    # filter for SNP position for VARIANT\n",
    "    H37Rv_ref_pos_for_SNP = SNP_func_annotation_DF[(SNP_func_annotation_DF.gene_id == gene_id) & (SNP_func_annotation_DF.gene_pos == gene_pos_SNP_i)].pos[0]\n",
    "    SNP_to_keep_filter = SNP_annotation_DF.pos == H37Rv_ref_pos_for_SNP # construct boolean filter\n",
    "\n",
    "    # filter Genotypes Matrix\n",
    "    SNP_genotypes_vector = SNP_genotypes_array_filtered[SNP_to_keep_filter , :]\n",
    "    SNP_genotypes_vector = [base_code_dict_r[allele_i] for allele_i in SNP_genotypes_vector[0]]\n",
    "\n",
    "    # func annotation for SNP\n",
    "    func_annot_for_SNP_list = list(SNP_func_annotation_DF[SNP_func_annotation_DF.pos == H37Rv_ref_pos_for_SNP].loc[: , ['pos','ref','gene_pos']].values[0]) + list(SNP_annotation_DF[SNP_annotation_DF.pos == H37Rv_ref_pos_for_SNP].alt)\n",
    "\n",
    "    # append the genotypes for the strains to the func annotation for the SNP\n",
    "    return func_annot_for_SNP_list + SNP_genotypes_vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Side-Analysis\n",
    "mutation: *mmpL5* Y300* 4.10.i1(1) - 1.1.1.1(292); check whether the single typing result for the L4 isolate is accurate? In Merker et al. 2020, it was found to be a marker for 1.1.1.1, which is in line with the remaining typing results.\n",
    "\n",
    "\n",
    "Some probing for this - there were 556 isolates typed as 1.1.1.1 in our dataset, 292/556 had this mutation. Furthermore, there were 2227 isolates typed as 4.10.i1 in our dataset, 1/2227 had this mutation and it was not one of the 14/31428 that were mis-typed according to our SNP barcoding scheme, so this isolate is certainly lineage 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1.1.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the number of isolates that are typed as sub-lineage **1.1.1.1**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "556"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(SNP_isolate_annotation_DF[SNP_isolate_annotation_DF.lineage_call == '1.1.1.1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the genotypes as H37Rv position **777,581** for the **1.1.1.1** isolates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos_777581_1_1_1_1_genotypes = get_genotypes_for_SNP(list(SNP_isolate_annotation_DF[SNP_isolate_annotation_DF.lineage_call == '1.1.1.1'].isolate_ID), 'Rv0676c', 900)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "556"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(pos_777581_1_1_1_1_genotypes[4:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many of these **1.1.1.1** isolates have the *mmpL5* Y300* allele?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "292"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(np.array(pos_777581_1_1_1_1_genotypes[4:]) == 'T')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.10.i1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the number of isolates that are typed as sub-lineage **4.10.i1**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2227"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(SNP_isolate_annotation_DF[SNP_isolate_annotation_DF.lineage_call == '4.10.i1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the genotypes as H37Rv position **777,581** for the **4.10.i1** isolates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos_777581_4_10_i1_genotypes = get_genotypes_for_SNP(list(SNP_isolate_annotation_DF[SNP_isolate_annotation_DF.lineage_call == '4.10.i1'].isolate_ID), 'Rv0676c', 900)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2227"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(pos_777581_4_10_i1_genotypes[4:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many of these **4.10.i1** isolates have the *mmpL5* Y300* allele?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(np.array(pos_777581_4_10_i1_genotypes[4:]) == 'T')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What was the breakdown of **alleles** for hte 2,227 **4.10.i1** isolates and what was the isolate ID of the single **4.10.i1** isolate that had the *mmpL5* Y300* allele?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "strains_of_interest = list(SNP_isolate_annotation_DF[SNP_isolate_annotation_DF.lineage_call == '4.10.i1'].isolate_ID)\n",
    "gene_id = 'Rv0676c'\n",
    "gene_pos_SNP_i = 900\n",
    "\n",
    "# get the number code that corresponds to the minor allele, 9 indicated bad quality call\n",
    "base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "base_code_dict_r = {0:'A', 1:'C', 2:'G', 3:'T', 9:'N'}\n",
    "\n",
    "# filter STRAINS\n",
    "strains_to_keep = [strain_i in strains_of_interest for strain_i in SNP_isolate_annotation_DF.isolate_ID] # construct boolean filter for strains to keep\n",
    "SNP_isolate_annotation_DF_filtered = SNP_isolate_annotation_DF[strains_to_keep]\n",
    "\n",
    "# filter Genotypes Matrix\n",
    "SNP_genotypes_array_filtered = SNP_genotypes_array[: , strains_to_keep]\n",
    "\n",
    "# filter for SNP position for VARIANT\n",
    "H37Rv_ref_pos_for_SNP = SNP_func_annotation_DF[(SNP_func_annotation_DF.gene_id == gene_id) & (SNP_func_annotation_DF.gene_pos == gene_pos_SNP_i)].pos[0]\n",
    "SNP_to_keep_filter = SNP_annotation_DF.pos == H37Rv_ref_pos_for_SNP # construct boolean filter\n",
    "\n",
    "# filter Genotypes Matrix\n",
    "SNP_genotypes_vector = SNP_genotypes_array_filtered[SNP_to_keep_filter , :]\n",
    "SNP_genotypes_vector = [base_code_dict_r[allele_i] for allele_i in SNP_genotypes_vector[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "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>4304</th>\n",
       "      <td>SAMEA1118012</td>\n",
       "      <td>4.10.i1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        isolate_ID lineage_call\n",
       "4304  SAMEA1118012      4.10.i1"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_isolate_annotation_DF_filtered[np.array(SNP_genotypes_vector) == 'T']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_isolate_annotation_DF_filtered[np.array(SNP_genotypes_vector) == 'T'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2200, 2)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape((SNP_isolate_annotation_DF_filtered[np.array(SNP_genotypes_vector) == 'G']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26, 2)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape((SNP_isolate_annotation_DF_filtered[np.array(SNP_genotypes_vector) == 'N']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] Check *atpE* for mmpR / mmpL5 double mutant isolates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [5.1] Check *atpE* for mmpR / mmpL5 double mutant isolates that were part of the Peruvian cluster (n=82)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get list of isolates that had both the *mmpR* insertion **and** the *mmpL5* deletion (from other notebook)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "mmpR_mmpL5_double_mutant_isolate_IDs = ['SAMEA5569929', 'SAMEA3402909', 'SAMEA5569553', 'Peru2918', 'Peru4515', 'SAMEA3392629', 'SAMEA5569985', 'SAMEA5569889', 'Peru4671', 'Peru4670', '01R1305', 'Peru2959', '02R0099', 'Peru4683', 'Peru4652', 'Peru4722', 'Peru4961', 'Peru4900', 'SAMEA1101329', 'Peru3418', 'SAMEA3414462', 'Peru2965', 'Peru3056', 'SAMEA2682981', 'Peru5003', '00R1547', 'SAMEA5570040', '02R0417', 'SAMEA5569638', '02R0948', 'Peru3376', 'Peru4521', 'Peru5012', '99R893', 'SAMEA5569636', 'SAMEA2683080', 'SAMEA3558289', 'Peru4988', 'SAMEA3401008', 'Peru4668', 'SAMEA2682695', 'Peru4685', 'SAMEA5569938', '01R0272', '1791', 'Peru4565', 'Peru5036', 'Peru4577', 'Peru4932', 'Peru4646', 'Peru4647', 'Peru4936', 'Peru4937', '02R1942', 'Peru3367', 'Peru4582', 'SAMEA5569790', 'Peru4549', 'Peru4057', 'SAMEA3558270', 'SAMEA2683134', '01R0451', 'SAMEA2682679', 'Peru4492', 'Peru3407', 'Peru5076', 'SAMN02414923', 'SAMEA5569666', 'Peru3342', 'Peru3389', 'Peru3408', 'Peru5445', 'Peru3405', 'Peru5115', 'Peru5114', 'Peru4707', 'SAMEA5569914', 'SAMEA5569642', 'SAMEA5570058', 'Peru3031', 'SAMEA3558288', 'Peru4919']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "82"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mmpR_mmpL5_double_mutant_isolate_IDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set of SNP sites in **atpE** in which strains are probed for the presence of an *alternate* allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "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>mut_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1461068_T</th>\n",
       "      <td>1461068</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>24</td>\n",
       "      <td>S</td>\n",
       "      <td>G8G</td>\n",
       "      <td>1461068_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461080_T</th>\n",
       "      <td>1461080</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>36</td>\n",
       "      <td>S</td>\n",
       "      <td>G12G</td>\n",
       "      <td>1461080_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461081_A</th>\n",
       "      <td>1461081</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>37</td>\n",
       "      <td>N</td>\n",
       "      <td>G13S</td>\n",
       "      <td>1461081_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461090_G</th>\n",
       "      <td>1461090</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>46</td>\n",
       "      <td>N</td>\n",
       "      <td>I16V</td>\n",
       "      <td>1461090_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461095_A</th>\n",
       "      <td>1461095</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>51</td>\n",
       "      <td>N</td>\n",
       "      <td>M17I</td>\n",
       "      <td>1461095_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461096_T</th>\n",
       "      <td>1461096</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>52</td>\n",
       "      <td>N</td>\n",
       "      <td>A18S</td>\n",
       "      <td>1461096_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461119_A</th>\n",
       "      <td>1461119</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>75</td>\n",
       "      <td>S</td>\n",
       "      <td>G25G</td>\n",
       "      <td>1461119_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461120_G</th>\n",
       "      <td>1461120</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>76</td>\n",
       "      <td>N</td>\n",
       "      <td>I26V</td>\n",
       "      <td>1461120_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461143_T</th>\n",
       "      <td>1461143</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>99</td>\n",
       "      <td>S</td>\n",
       "      <td>N33N</td>\n",
       "      <td>1461143_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461154_G</th>\n",
       "      <td>1461154</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>110</td>\n",
       "      <td>N</td>\n",
       "      <td>S37C</td>\n",
       "      <td>1461154_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461159_A</th>\n",
       "      <td>1461159</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>115</td>\n",
       "      <td>N</td>\n",
       "      <td>V39I</td>\n",
       "      <td>1461159_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461176_A</th>\n",
       "      <td>1461176</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>132</td>\n",
       "      <td>S</td>\n",
       "      <td>E44E</td>\n",
       "      <td>1461176_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461176_C</th>\n",
       "      <td>1461176</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>132</td>\n",
       "      <td>N</td>\n",
       "      <td>E44D</td>\n",
       "      <td>1461176_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461179_A</th>\n",
       "      <td>1461179</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>135</td>\n",
       "      <td>S</td>\n",
       "      <td>A45A</td>\n",
       "      <td>1461179_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461186_A</th>\n",
       "      <td>1461186</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>142</td>\n",
       "      <td>S</td>\n",
       "      <td>R48R</td>\n",
       "      <td>1461186_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461194_A</th>\n",
       "      <td>1461194</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>150</td>\n",
       "      <td>N</td>\n",
       "      <td>F50L</td>\n",
       "      <td>1461194_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461199_T</th>\n",
       "      <td>1461199</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>155</td>\n",
       "      <td>N</td>\n",
       "      <td>P52L</td>\n",
       "      <td>1461199_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461206_T</th>\n",
       "      <td>1461206</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>162</td>\n",
       "      <td>S</td>\n",
       "      <td>F54F</td>\n",
       "      <td>1461206_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461209_T</th>\n",
       "      <td>1461209</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>165</td>\n",
       "      <td>S</td>\n",
       "      <td>I55I</td>\n",
       "      <td>1461209_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461215_T</th>\n",
       "      <td>1461215</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>171</td>\n",
       "      <td>S</td>\n",
       "      <td>V57V</td>\n",
       "      <td>1461215_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461221_A</th>\n",
       "      <td>1461221</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>177</td>\n",
       "      <td>S</td>\n",
       "      <td>L59L</td>\n",
       "      <td>1461221_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461227_A</th>\n",
       "      <td>1461227</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>183</td>\n",
       "      <td>S</td>\n",
       "      <td>E61E</td>\n",
       "      <td>1461227_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461230_A</th>\n",
       "      <td>1461230</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>186</td>\n",
       "      <td>S</td>\n",
       "      <td>A62A</td>\n",
       "      <td>1461230_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461232_T</th>\n",
       "      <td>1461232</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>188</td>\n",
       "      <td>N</td>\n",
       "      <td>A63V</td>\n",
       "      <td>1461232_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461240_G</th>\n",
       "      <td>1461240</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>196</td>\n",
       "      <td>N</td>\n",
       "      <td>I66V</td>\n",
       "      <td>1461240_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461251_T</th>\n",
       "      <td>1461251</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>207</td>\n",
       "      <td>S</td>\n",
       "      <td>A69A</td>\n",
       "      <td>1461251_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461277_T</th>\n",
       "      <td>1461277</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>233</td>\n",
       "      <td>N</td>\n",
       "      <td>T78I</td>\n",
       "      <td>1461277_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461283_C</th>\n",
       "      <td>1461283</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>239</td>\n",
       "      <td>N</td>\n",
       "      <td>V80A</td>\n",
       "      <td>1461283_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461284_T</th>\n",
       "      <td>1461284</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>240</td>\n",
       "      <td>S</td>\n",
       "      <td>V80V</td>\n",
       "      <td>1461284_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  \\\n",
       "1461068_T  1461068   C   T     Essential      atpE  Rv1305       24        S   \n",
       "1461080_T  1461080   C   T     Essential      atpE  Rv1305       36        S   \n",
       "1461081_A  1461081   G   A     Essential      atpE  Rv1305       37        N   \n",
       "1461090_G  1461090   A   G     Essential      atpE  Rv1305       46        N   \n",
       "1461095_A  1461095   G   A     Essential      atpE  Rv1305       51        N   \n",
       "1461096_T  1461096   G   T     Essential      atpE  Rv1305       52        N   \n",
       "1461119_A  1461119   T   A     Essential      atpE  Rv1305       75        S   \n",
       "1461120_G  1461120   A   G     Essential      atpE  Rv1305       76        N   \n",
       "1461143_T  1461143   C   T     Essential      atpE  Rv1305       99        S   \n",
       "1461154_G  1461154   C   G     Essential      atpE  Rv1305      110        N   \n",
       "1461159_A  1461159   G   A     Essential      atpE  Rv1305      115        N   \n",
       "1461176_A  1461176   G   A     Essential      atpE  Rv1305      132        S   \n",
       "1461176_C  1461176   G   C     Essential      atpE  Rv1305      132        N   \n",
       "1461179_A  1461179   G   A     Essential      atpE  Rv1305      135        S   \n",
       "1461186_A  1461186   C   A     Essential      atpE  Rv1305      142        S   \n",
       "1461194_A  1461194   C   A     Essential      atpE  Rv1305      150        N   \n",
       "1461199_T  1461199   C   T     Essential      atpE  Rv1305      155        N   \n",
       "1461206_T  1461206   C   T     Essential      atpE  Rv1305      162        S   \n",
       "1461209_T  1461209   C   T     Essential      atpE  Rv1305      165        S   \n",
       "1461215_T  1461215   C   T     Essential      atpE  Rv1305      171        S   \n",
       "1461221_A  1461221   G   A     Essential      atpE  Rv1305      177        S   \n",
       "1461227_A  1461227   G   A     Essential      atpE  Rv1305      183        S   \n",
       "1461230_A  1461230   G   A     Essential      atpE  Rv1305      186        S   \n",
       "1461232_T  1461232   C   T     Essential      atpE  Rv1305      188        N   \n",
       "1461240_G  1461240   A   G     Essential      atpE  Rv1305      196        N   \n",
       "1461251_T  1461251   G   T     Essential      atpE  Rv1305      207        S   \n",
       "1461277_T  1461277   C   T     Essential      atpE  Rv1305      233        N   \n",
       "1461283_C  1461283   T   C     Essential      atpE  Rv1305      239        N   \n",
       "1461284_T  1461284   C   T     Essential      atpE  Rv1305      240        S   \n",
       "\n",
       "          AA_change    mut_key  \n",
       "1461068_T       G8G  1461068_T  \n",
       "1461080_T      G12G  1461080_T  \n",
       "1461081_A      G13S  1461081_A  \n",
       "1461090_G      I16V  1461090_G  \n",
       "1461095_A      M17I  1461095_A  \n",
       "1461096_T      A18S  1461096_T  \n",
       "1461119_A      G25G  1461119_A  \n",
       "1461120_G      I26V  1461120_G  \n",
       "1461143_T      N33N  1461143_T  \n",
       "1461154_G      S37C  1461154_G  \n",
       "1461159_A      V39I  1461159_A  \n",
       "1461176_A      E44E  1461176_A  \n",
       "1461176_C      E44D  1461176_C  \n",
       "1461179_A      A45A  1461179_A  \n",
       "1461186_A      R48R  1461186_A  \n",
       "1461194_A      F50L  1461194_A  \n",
       "1461199_T      P52L  1461199_T  \n",
       "1461206_T      F54F  1461206_T  \n",
       "1461209_T      I55I  1461209_T  \n",
       "1461215_T      V57V  1461215_T  \n",
       "1461221_A      L59L  1461221_A  \n",
       "1461227_A      E61E  1461227_A  \n",
       "1461230_A      A62A  1461230_A  \n",
       "1461232_T      A63V  1461232_T  \n",
       "1461240_G      I66V  1461240_G  \n",
       "1461251_T      A69A  1461251_T  \n",
       "1461277_T      T78I  1461277_T  \n",
       "1461283_C      V80A  1461283_C  \n",
       "1461284_T      V80V  1461284_T  "
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_func_annotation_DF[SNP_func_annotation_DF.gene_name == 'atpE']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29, 10)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_func_annotation_DF[SNP_func_annotation_DF.gene_name == 'atpE'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "atpE_muts_for_mmpR_mmpL5_double_mutants = get_strains_with_SNPs_in_gene('Rv1305', mmpR_mmpL5_double_mutant_isolate_IDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "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>isolate_ID</th>\n",
       "      <th>lineage_call</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>mut_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [pos, ref, alt, isolate_ID, lineage_call, gene_category, gene_name, gene_id, gene_pos, SNP_type, AA_change, mut_key]\n",
       "Index: []"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atpE_muts_for_mmpR_mmpL5_double_mutants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [5.2] Check *atpE* for mmpR / mmpL5 double mutant isolates that were note part of Peruvian cluster (n=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get list of isolates that had both a *mmpR* frameshift **and** a *mmpL5* insertion (from other notebook)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "mmpR_mmpL5_double_mutant_isolate_IDs = ['SAMN12551172','SAMD00117890']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mmpR_mmpL5_double_mutant_isolate_IDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set of SNP sites in **atpE** in which strains are probed for the presence of an *alternate* allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "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>mut_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1461068_T</th>\n",
       "      <td>1461068</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>24</td>\n",
       "      <td>S</td>\n",
       "      <td>G8G</td>\n",
       "      <td>1461068_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461080_T</th>\n",
       "      <td>1461080</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>36</td>\n",
       "      <td>S</td>\n",
       "      <td>G12G</td>\n",
       "      <td>1461080_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461081_A</th>\n",
       "      <td>1461081</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>37</td>\n",
       "      <td>N</td>\n",
       "      <td>G13S</td>\n",
       "      <td>1461081_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461090_G</th>\n",
       "      <td>1461090</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>46</td>\n",
       "      <td>N</td>\n",
       "      <td>I16V</td>\n",
       "      <td>1461090_G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1461095_A</th>\n",
       "      <td>1461095</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Essential</td>\n",
       "      <td>atpE</td>\n",
       "      <td>Rv1305</td>\n",
       "      <td>51</td>\n",
       "      <td>N</td>\n",
       "      <td>M17I</td>\n",
       "      <td>1461095_A</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               pos ref alt gene_category gene_name gene_id gene_pos SNP_type  \\\n",
       "1461068_T  1461068   C   T     Essential      atpE  Rv1305       24        S   \n",
       "1461080_T  1461080   C   T     Essential      atpE  Rv1305       36        S   \n",
       "1461081_A  1461081   G   A     Essential      atpE  Rv1305       37        N   \n",
       "1461090_G  1461090   A   G     Essential      atpE  Rv1305       46        N   \n",
       "1461095_A  1461095   G   A     Essential      atpE  Rv1305       51        N   \n",
       "\n",
       "          AA_change    mut_key  \n",
       "1461068_T       G8G  1461068_T  \n",
       "1461080_T      G12G  1461080_T  \n",
       "1461081_A      G13S  1461081_A  \n",
       "1461090_G      I16V  1461090_G  \n",
       "1461095_A      M17I  1461095_A  "
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_func_annotation_DF[SNP_func_annotation_DF.gene_name == 'atpE'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29, 10)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(SNP_func_annotation_DF[SNP_func_annotation_DF.gene_name == 'atpE'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "atpE_muts_for_mmpR_mmpL5_double_mutants = get_strains_with_SNPs_in_gene('Rv1305', mmpR_mmpL5_double_mutant_isolate_IDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "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>isolate_ID</th>\n",
       "      <th>lineage_call</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>mut_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [pos, ref, alt, isolate_ID, lineage_call, gene_category, gene_name, gene_id, gene_pos, SNP_type, AA_change, mut_key]\n",
       "Index: []"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atpE_muts_for_mmpR_mmpL5_double_mutants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['SAMN12551172', 'SAMD00117890']"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mmpR_mmpL5_double_mutant_isolate_IDs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [6] Check *rrs* for *eis* promoter / *eis* double mutant isolates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>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>mut_key</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1471848_A</th>\n",
       "      <td>1471848</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>3</td>\n",
       "      <td>N</td>\n",
       "      <td>F1L</td>\n",
       "      <td>1471848_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1471850_T</th>\n",
       "      <td>1471850</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>5</td>\n",
       "      <td>N</td>\n",
       "      <td>C2F</td>\n",
       "      <td>1471850_T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1471861_A</th>\n",
       "      <td>1471861</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>16</td>\n",
       "      <td>N</td>\n",
       "      <td>L6M</td>\n",
       "      <td>1471861_A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1471861_C</th>\n",
       "      <td>1471861</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>16</td>\n",
       "      <td>S</td>\n",
       "      <td>L6L</td>\n",
       "      <td>1471861_C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1471861_G</th>\n",
       "      <td>1471861</td>\n",
       "      <td>T</td>\n",
       "      <td>G</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>16</td>\n",
       "      <td>N</td>\n",
       "      <td>L6V</td>\n",
       "      <td>1471861_G</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               pos ref alt gene_category gene_name gene_id gene_pos SNP_type  \\\n",
       "1471848_A  1471848   T   A          None      None  Rvnr01        3        N   \n",
       "1471850_T  1471850   G   T          None      None  Rvnr01        5        N   \n",
       "1471861_A  1471861   T   A          None      None  Rvnr01       16        N   \n",
       "1471861_C  1471861   T   C          None      None  Rvnr01       16        S   \n",
       "1471861_G  1471861   T   G          None      None  Rvnr01       16        N   \n",
       "\n",
       "          AA_change    mut_key  \n",
       "1471848_A       F1L  1471848_A  \n",
       "1471850_T       C2F  1471850_T  \n",
       "1471861_A       L6M  1471861_A  \n",
       "1471861_C       L6L  1471861_C  \n",
       "1471861_G       L6V  1471861_G  "
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_func_annotation_DF[(SNP_func_annotation_DF.pos >= 1471846) & (SNP_func_annotation_DF.pos <= 1473382)].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get list of isolates that were additional 12 *eis* C-14T mutants with AG MICs and/or had both the *eis* promoter SNP **and** the *eis* LOF mutations\n",
    "\n",
    "Check additional 12 *eis* C-14T mutants with AG MICs Strains: \n",
    "1. IT952\n",
    "1. IT947\n",
    "1. IT77\n",
    "1. IT634\n",
    "1. IT524\n",
    "1. IT233\n",
    "1. IT184\n",
    "1. IT123\n",
    "1. IT1070\n",
    "1. 655-19\n",
    "1. 622-19\n",
    "1. 168-19\n",
    "\n",
    "Other strains that have a *eis* promoter mutation & *eis* LOF mutation:\n",
    "1. SAMN07956543\n",
    "1. SAMN02584676\n",
    "1. SAMEA1016073\n",
    "1. SAMN04633319\n",
    "1. SAMEA1403638\n",
    "1. SAMEA1403685\n",
    "1. SAMN02419586\n",
    "1. SAMN02419559\n",
    "1. SAMN02419535\n",
    "1. SAMN07236283\n",
    "1. SAMN02419543\n",
    "1. SAMN07236283\n",
    "1. Peru3354\n",
    "1. Peru2946\n",
    "1. SAMN08376196\n",
    "1. SAMN08709032\n",
    "1. SAMN06210015\n",
    "1. SAMN02584612"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check to see that strains that (additional 12 *eis* C-14T mutants with AG MICs) have the G-14A *eis* promoter mutation (@ position 2715346)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "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>mut_key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</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",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>168-19</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>622-19</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>655-19</td>\n",
       "      <td>4.3.i3.1</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT1070</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT123</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT184</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT233</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT524</td>\n",
       "      <td>4.10.i2</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT634</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT77</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT947</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>2715346_A</td>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>IT952</td>\n",
       "      <td>4.10.i2</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>Rv2416c_Rv2417c</td>\n",
       "      <td>None</td>\n",
       "      <td>I</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      mut_key      pos ref alt isolate_ID  lineage_call gene_category  \\\n",
       "0   2715346_A  2715346   G   A     168-19     2.2.1.1.1          None   \n",
       "1   2715346_A  2715346   G   A     622-19  2.2.1.1.1.i3          None   \n",
       "2   2715346_A  2715346   G   A     655-19      4.3.i3.1          None   \n",
       "3   2715346_A  2715346   G   A     IT1070  2.2.1.1.1.i3          None   \n",
       "4   2715346_A  2715346   G   A      IT123     2.2.1.1.1          None   \n",
       "5   2715346_A  2715346   G   A      IT184     2.2.1.1.1          None   \n",
       "6   2715346_A  2715346   G   A      IT233  2.2.1.1.1.i3          None   \n",
       "7   2715346_A  2715346   G   A      IT524       4.10.i2          None   \n",
       "8   2715346_A  2715346   G   A      IT634  2.2.1.1.1.i3          None   \n",
       "9   2715346_A  2715346   G   A       IT77  2.2.1.1.1.i3          None   \n",
       "10  2715346_A  2715346   G   A      IT947  2.2.1.1.1.i3          None   \n",
       "11  2715346_A  2715346   G   A      IT952       4.10.i2          None   \n",
       "\n",
       "   gene_name          gene_id gene_pos SNP_type AA_change  \n",
       "0       None  Rv2416c_Rv2417c     None        I      None  \n",
       "1       None  Rv2416c_Rv2417c     None        I      None  \n",
       "2       None  Rv2416c_Rv2417c     None        I      None  \n",
       "3       None  Rv2416c_Rv2417c     None        I      None  \n",
       "4       None  Rv2416c_Rv2417c     None        I      None  \n",
       "5       None  Rv2416c_Rv2417c     None        I      None  \n",
       "6       None  Rv2416c_Rv2417c     None        I      None  \n",
       "7       None  Rv2416c_Rv2417c     None        I      None  \n",
       "8       None  Rv2416c_Rv2417c     None        I      None  \n",
       "9       None  Rv2416c_Rv2417c     None        I      None  \n",
       "10      None  Rv2416c_Rv2417c     None        I      None  \n",
       "11      None  Rv2416c_Rv2417c     None        I      None  "
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_strains_with_SNPs_in_gene('Rv2416c_Rv2417c', ['IT952','IT947','IT77','IT634','IT524','IT233','IT184','IT123','IT1070','655-19','622-19','168-19'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30\n"
     ]
    }
   ],
   "source": [
    "strains_of_interest = ['IT952','IT947','IT77','IT634','IT524','IT233','IT184','IT123','IT1070','655-19','622-19','168-19','SAMN07956543','SAMN02584676','SAMEA1016073','SAMN04633319','SAMEA1403638','SAMEA1403685','SAMN02419586','SAMN02419559','SAMN02419535','SAMN07236283','SAMN02419543','SAMN07236283','Peru3354','Peru2946','SAMN08376196','SAMN08709032','SAMN06210015','SAMN02584612']\n",
    "print(len(strains_of_interest))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter Isolate annotation file\n",
    "strains_to_keep = [strain_i in strains_of_interest for strain_i in SNP_isolate_annotation_DF.isolate_ID] # construct boolean filter for strains to keep\n",
    "SNP_isolate_annotation_DF_filtered = SNP_isolate_annotation_DF[strains_to_keep]\n",
    "SNP_isolate_annotation_DF_filtered.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "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>SAMEA1016073</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SAMN08376196</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SAMEA1403685</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SAMEA1403638</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SAMN06210015</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>SAMN02419535</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>SAMN02419543</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SAMN07236283</td>\n",
       "      <td>1.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>SAMN02584676</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>SAMN02419586</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Peru2946</td>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>SAMN02419559</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>SAMN02584612</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>SAMN08709032</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>SAMN04633319</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Peru3354</td>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>SAMN07956543</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>168-19</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>622-19</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>655-19</td>\n",
       "      <td>4.3.i3.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>IT1070</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>IT123</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>IT184</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>IT233</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>IT524</td>\n",
       "      <td>4.10.i2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>IT634</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>IT77</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>IT947</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>IT952</td>\n",
       "      <td>4.10.i2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      isolate_ID  lineage_call\n",
       "0   SAMEA1016073  2.2.1.1.1.i3\n",
       "1   SAMN08376196     2.2.1.1.1\n",
       "2   SAMEA1403685  2.2.1.1.1.i3\n",
       "3   SAMEA1403638  2.2.1.1.1.i3\n",
       "4   SAMN06210015     2.2.1.1.1\n",
       "5   SAMN02419535  2.2.1.1.1.i3\n",
       "6   SAMN02419543  2.2.1.1.1.i3\n",
       "7   SAMN07236283       1.1.1.1\n",
       "8   SAMN02584676  2.2.1.1.1.i3\n",
       "9   SAMN02419586  2.2.1.1.1.i3\n",
       "10      Peru2946  4.1.i1.1.1.1\n",
       "11  SAMN02419559  2.2.1.1.1.i3\n",
       "12  SAMN02584612     2.2.1.1.1\n",
       "13  SAMN08709032     2.2.1.1.1\n",
       "14  SAMN04633319  2.2.1.1.1.i3\n",
       "15      Peru3354  4.1.i1.1.1.1\n",
       "16  SAMN07956543  2.2.1.1.1.i3\n",
       "17        168-19     2.2.1.1.1\n",
       "18        622-19  2.2.1.1.1.i3\n",
       "19        655-19      4.3.i3.1\n",
       "20        IT1070  2.2.1.1.1.i3\n",
       "21         IT123     2.2.1.1.1\n",
       "22         IT184     2.2.1.1.1\n",
       "23         IT233  2.2.1.1.1.i3\n",
       "24         IT524       4.10.i2\n",
       "25         IT634  2.2.1.1.1.i3\n",
       "26          IT77  2.2.1.1.1.i3\n",
       "27         IT947  2.2.1.1.1.i3\n",
       "28         IT952       4.10.i2"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_isolate_annotation_DF_filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrs_genotypes_for_strains_of_interest = pd.DataFrame(columns = ['pos','ref','gene_pos','alt'] + list(SNP_isolate_annotation_DF_filtered.isolate_ID))\n",
    "\n",
    "#add strain lineage calls\n",
    "rrs_genotypes_for_strains_of_interest.loc['lineage_call', :] = [np.nan, np.nan, np.nan, np.nan] + list(SNP_isolate_annotation_DF_filtered.lineage_call)\n",
    "\n",
    "# A1401G\n",
    "rrs_genotypes_for_strains_of_interest.loc['rrs_SNP_1', :] = get_genotypes_for_SNP(strains_of_interest, 'Rvnr01', 1401)\n",
    "\n",
    "# C1402T\n",
    "rrs_genotypes_for_strains_of_interest.loc['rrs_SNP_2', :] = get_genotypes_for_SNP(strains_of_interest, 'Rvnr01', 1402)\n",
    "\n",
    "# G1484T\n",
    "rrs_genotypes_for_strains_of_interest.loc['rrs_SNP_3', :] = get_genotypes_for_SNP(strains_of_interest, 'Rvnr01', 1484)"
   ]
  },
  {
   "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>lineage_call</th>\n",
       "      <th>rrs_SNP_1</th>\n",
       "      <th>rrs_SNP_2</th>\n",
       "      <th>rrs_SNP_3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>pos</th>\n",
       "      <td>NaN</td>\n",
       "      <td>1473246</td>\n",
       "      <td>1473247</td>\n",
       "      <td>1473329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ref</th>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_pos</th>\n",
       "      <td>NaN</td>\n",
       "      <td>1401</td>\n",
       "      <td>1402</td>\n",
       "      <td>1484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alt</th>\n",
       "      <td>NaN</td>\n",
       "      <td>[G]</td>\n",
       "      <td>[A, T]</td>\n",
       "      <td>[C, T]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMEA1016073</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN08376196</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMEA1403685</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMEA1403638</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN06210015</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02419535</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02419543</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN07236283</th>\n",
       "      <td>1.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02584676</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02419586</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Peru2946</th>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02419559</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN02584612</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN08709032</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN04633319</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Peru3354</th>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SAMN07956543</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>168-19</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>622-19</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>655-19</th>\n",
       "      <td>4.3.i3.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT1070</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT123</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT184</th>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT233</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT524</th>\n",
       "      <td>4.10.i2</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT634</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT77</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT947</th>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IT952</th>\n",
       "      <td>4.10.i2</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              lineage_call rrs_SNP_1 rrs_SNP_2 rrs_SNP_3\n",
       "pos                    NaN   1473246   1473247   1473329\n",
       "ref                    NaN         A         C         G\n",
       "gene_pos               NaN      1401      1402      1484\n",
       "alt                    NaN       [G]    [A, T]    [C, T]\n",
       "SAMEA1016073  2.2.1.1.1.i3         A         C         G\n",
       "SAMN08376196     2.2.1.1.1         A         C         G\n",
       "SAMEA1403685  2.2.1.1.1.i3         A         C         G\n",
       "SAMEA1403638  2.2.1.1.1.i3         A         C         G\n",
       "SAMN06210015     2.2.1.1.1         A         C         G\n",
       "SAMN02419535  2.2.1.1.1.i3         A         C         G\n",
       "SAMN02419543  2.2.1.1.1.i3         A         C         G\n",
       "SAMN07236283       1.1.1.1         A         C         G\n",
       "SAMN02584676  2.2.1.1.1.i3         A         C         G\n",
       "SAMN02419586  2.2.1.1.1.i3         A         C         G\n",
       "Peru2946      4.1.i1.1.1.1         A         C         G\n",
       "SAMN02419559  2.2.1.1.1.i3         A         C         G\n",
       "SAMN02584612     2.2.1.1.1         A         C         G\n",
       "SAMN08709032     2.2.1.1.1         A         C         G\n",
       "SAMN04633319  2.2.1.1.1.i3         A         C         G\n",
       "Peru3354      4.1.i1.1.1.1         A         C         G\n",
       "SAMN07956543  2.2.1.1.1.i3         A         C         G\n",
       "168-19           2.2.1.1.1         A         C         G\n",
       "622-19        2.2.1.1.1.i3         A         C         G\n",
       "655-19            4.3.i3.1         A         C         G\n",
       "IT1070        2.2.1.1.1.i3         A         C         G\n",
       "IT123            2.2.1.1.1         A         C         G\n",
       "IT184            2.2.1.1.1         A         C         G\n",
       "IT233         2.2.1.1.1.i3         A         C         N\n",
       "IT524              4.10.i2         N         N         N\n",
       "IT634         2.2.1.1.1.i3         A         C         G\n",
       "IT77          2.2.1.1.1.i3         N         N         N\n",
       "IT947         2.2.1.1.1.i3         A         C         G\n",
       "IT952              4.10.i2         A         C         G"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rrs_genotypes_for_strains_of_interest.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output as CSV file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrs_genotypes_for_strains_of_interest.T.to_csv('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/rrs_resistance_mut_check_extra_strains_and_eis_prom_eis_double_muts.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check that the *rrs* mutations using the other function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrs_muts_for_strains_of_interest = get_strains_with_SNPs_in_gene('Rvnr01', strains_of_interest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrs_muts_for_strains_of_interest.drop(['SNP_type','AA_change','gene_category'],axis=1,inplace=True)\n",
    "rrs_muts_for_strains_of_interest.loc[:, 'gene_name'] = ['rrs']*np.shape(rrs_muts_for_strains_of_interest)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "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>mut_key</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1472359_C</td>\n",
       "      <td>1472359</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>655-19</td>\n",
       "      <td>4.3.i3.1</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>514</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMEA1016073</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMEA1403685</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMEA1403638</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN02419535</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN02419543</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN07236283</td>\n",
       "      <td>1.1.1.1</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN02584676</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN02419586</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN02419559</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN04633319</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>SAMN07956543</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>168-19</td>\n",
       "      <td>2.2.1.1.1</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>IT1070</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>IT233</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>IT634</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1472362_T</td>\n",
       "      <td>1472362</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>IT947</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>1472751_G</td>\n",
       "      <td>1472751</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Peru2946</td>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1472751_G</td>\n",
       "      <td>1472751</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>Peru3354</td>\n",
       "      <td>4.1.i1.1.1.1</td>\n",
       "      <td>rrs</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>906</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      mut_key      pos ref alt    isolate_ID  lineage_call gene_name gene_id  \\\n",
       "0   1472359_C  1472359   A   C        655-19      4.3.i3.1       rrs  Rvnr01   \n",
       "1   1472362_T  1472362   C   T  SAMEA1016073  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "2   1472362_T  1472362   C   T  SAMEA1403685  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "3   1472362_T  1472362   C   T  SAMEA1403638  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "4   1472362_T  1472362   C   T  SAMN02419535  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "5   1472362_T  1472362   C   T  SAMN02419543  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "6   1472362_T  1472362   C   T  SAMN07236283       1.1.1.1       rrs  Rvnr01   \n",
       "7   1472362_T  1472362   C   T  SAMN02584676  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "8   1472362_T  1472362   C   T  SAMN02419586  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "9   1472362_T  1472362   C   T  SAMN02419559  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "10  1472362_T  1472362   C   T  SAMN04633319  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "11  1472362_T  1472362   C   T  SAMN07956543  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "12  1472362_T  1472362   C   T        168-19     2.2.1.1.1       rrs  Rvnr01   \n",
       "13  1472362_T  1472362   C   T        IT1070  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "14  1472362_T  1472362   C   T         IT233  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "15  1472362_T  1472362   C   T         IT634  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "16  1472362_T  1472362   C   T         IT947  2.2.1.1.1.i3       rrs  Rvnr01   \n",
       "17  1472751_G  1472751   A   G      Peru2946  4.1.i1.1.1.1       rrs  Rvnr01   \n",
       "18  1472751_G  1472751   A   G      Peru3354  4.1.i1.1.1.1       rrs  Rvnr01   \n",
       "\n",
       "   gene_pos  \n",
       "0       514  \n",
       "1       517  \n",
       "2       517  \n",
       "3       517  \n",
       "4       517  \n",
       "5       517  \n",
       "6       517  \n",
       "7       517  \n",
       "8       517  \n",
       "9       517  \n",
       "10      517  \n",
       "11      517  \n",
       "12      517  \n",
       "13      517  \n",
       "14      517  \n",
       "15      517  \n",
       "16      517  \n",
       "17      906  \n",
       "18      906  "
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rrs_muts_for_strains_of_interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [7] Check co-occurence of 4 most common *eis* promoter mutations with 3 *rrs* aminoglycoside resistance conferring mutations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Jnawali et al. illustrates why checking for *rrs* A1401G, C1402T and G1484T is important (i.e. it describes two *eis* -14 mutants with different LoF mutations that coincide with *rrs* 1401)\n",
    "\n",
    "- How often do each of the four classical eis promoter mutations (i.e. c-14t, c-37t, g-10a, and c-12t) **occur without** one of the three rrs mutations in out sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*eis* promoter mutations associated with KAN resistance\n",
    "\n",
    "- eis G-10A, 2715342, C>T\n",
    "- eis C-12T, 2715344, G>A\n",
    "- eis C-14T, 2715346, G>A\n",
    "- eis G-37T, 2715369, C>A\n",
    "\n",
    "*rrs* mutations associated with AG resistance\n",
    "- A1401G, 1473246, A>G\n",
    "- C1402T, 1473247, C>T\n",
    "- G1484T, 1473329, G>T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNP_sites_of_interest = [2715342, 2715344, 2715346, 2715369, 1473246, 1473247, 1473329]\n",
    "\n",
    "# get the number code that corresponds to the minor allele, 9 indicated bad quality call\n",
    "base_code_dict = {'A':0, 'C':1, 'G':2, 'T':3}\n",
    "base_code_dict_r = {0:'A', 1:'C', 2:'G', 3:'T', 9:'N'}\n",
    "\n",
    "# filter SNP sites\n",
    "SNP_sites_to_keep = [SNP_i in SNP_sites_of_interest for SNP_i in SNP_annotation_DF.pos] # construct boolean filter for SNP sites to keep\n",
    "SNP_annotation_DF_filtered = SNP_annotation_DF[SNP_sites_to_keep]\n",
    "SNP_annotation_DF_filtered.reset_index(inplace = True, drop = True)\n",
    "SNP_annotation_DF_filtered.loc[:, 'variant'] = ['rrs A1401G', 'rrs C1402T', 'rrs G1484T', 'eis G-10A', 'eis C-12T', 'eis C-14T', 'eis G-37T']\n",
    "SNP_annotation_DF_filtered.loc[:, 'resistance_allele'] = ['G', 'T', 'T', 'T', 'A', 'A', 'A']\n",
    "\n",
    "# filter Genotypes Matrix\n",
    "SNP_genotypes_array_filtered = SNP_genotypes_array[SNP_sites_to_keep , :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>gene_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>variant</th>\n",
       "      <th>resistance_allele</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1473246</td>\n",
       "      <td>A</td>\n",
       "      <td>[G]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1400.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs A1401G</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1473247</td>\n",
       "      <td>C</td>\n",
       "      <td>[A, T]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1401.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs C1402T</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1473329</td>\n",
       "      <td>G</td>\n",
       "      <td>[C, T]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1483.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs G1484T</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2715342</td>\n",
       "      <td>C</td>\n",
       "      <td>[T, G]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis G-10A</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2715344</td>\n",
       "      <td>G</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis C-12T</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis C-14T</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2715369</td>\n",
       "      <td>C</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis G-37T</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pos ref     alt     gene_id  gene_pos gene_category gene_name  \\\n",
       "0  1473246   A     [G]      Rvnr01    1400.0          None      None   \n",
       "1  1473247   C  [A, T]      Rvnr01    1401.0          None      None   \n",
       "2  1473329   G  [C, T]      Rvnr01    1483.0          None      None   \n",
       "3  2715342   C  [T, G]  intergenic       NaN          None      None   \n",
       "4  2715344   G     [A]  intergenic       NaN          None      None   \n",
       "5  2715346   G     [A]  intergenic       NaN          None      None   \n",
       "6  2715369   C     [A]  intergenic       NaN          None      None   \n",
       "\n",
       "      variant resistance_allele  \n",
       "0  rrs A1401G                 G  \n",
       "1  rrs C1402T                 T  \n",
       "2  rrs G1484T                 T  \n",
       "3   eis G-10A                 T  \n",
       "4   eis C-12T                 A  \n",
       "5   eis C-14T                 A  \n",
       "6   eis G-37T                 A  "
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF_filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 9, ..., 9, 0, 0],\n",
       "       [1, 1, 9, ..., 9, 1, 1],\n",
       "       [2, 2, 9, ..., 9, 2, 2],\n",
       "       ...,\n",
       "       [2, 2, 2, ..., 2, 2, 2],\n",
       "       [2, 2, 2, ..., 0, 0, 0],\n",
       "       [1, 1, 1, ..., 1, 1, 1]], dtype=int8)"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_genotypes_array_filtered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7, 31440)"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_genotypes_array_filtered.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the number of isolates with each SNV (just to double check)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_isolates_with_variant_allele = []\n",
    "for SNV_index, SNV_res_allele in zip(SNP_annotation_DF_filtered.index,SNP_annotation_DF_filtered.resistance_allele):\n",
    "    \n",
    "    num_isolates_with_variant_allele.append(np.sum(SNP_genotypes_array_filtered[SNV_index,:] == base_code_dict[SNV_res_allele]))\n",
    "    \n",
    "SNP_annotation_DF_filtered.loc[:,'num_isolates'] = num_isolates_with_variant_allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "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_id</th>\n",
       "      <th>gene_pos</th>\n",
       "      <th>gene_category</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>variant</th>\n",
       "      <th>resistance_allele</th>\n",
       "      <th>num_isolates</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1473246</td>\n",
       "      <td>A</td>\n",
       "      <td>[G]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1400.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs A1401G</td>\n",
       "      <td>G</td>\n",
       "      <td>1452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1473247</td>\n",
       "      <td>C</td>\n",
       "      <td>[A, T]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1401.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs C1402T</td>\n",
       "      <td>T</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1473329</td>\n",
       "      <td>G</td>\n",
       "      <td>[C, T]</td>\n",
       "      <td>Rvnr01</td>\n",
       "      <td>1483.0</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>rrs G1484T</td>\n",
       "      <td>T</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2715342</td>\n",
       "      <td>C</td>\n",
       "      <td>[T, G]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis G-10A</td>\n",
       "      <td>T</td>\n",
       "      <td>293</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2715344</td>\n",
       "      <td>G</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis C-12T</td>\n",
       "      <td>A</td>\n",
       "      <td>332</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2715346</td>\n",
       "      <td>G</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis C-14T</td>\n",
       "      <td>A</td>\n",
       "      <td>181</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2715369</td>\n",
       "      <td>C</td>\n",
       "      <td>[A]</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>eis G-37T</td>\n",
       "      <td>A</td>\n",
       "      <td>285</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pos ref     alt     gene_id  gene_pos gene_category gene_name  \\\n",
       "0  1473246   A     [G]      Rvnr01    1400.0          None      None   \n",
       "1  1473247   C  [A, T]      Rvnr01    1401.0          None      None   \n",
       "2  1473329   G  [C, T]      Rvnr01    1483.0          None      None   \n",
       "3  2715342   C  [T, G]  intergenic       NaN          None      None   \n",
       "4  2715344   G     [A]  intergenic       NaN          None      None   \n",
       "5  2715346   G     [A]  intergenic       NaN          None      None   \n",
       "6  2715369   C     [A]  intergenic       NaN          None      None   \n",
       "\n",
       "      variant resistance_allele  num_isolates  \n",
       "0  rrs A1401G                 G          1452  \n",
       "1  rrs C1402T                 T            18  \n",
       "2  rrs G1484T                 T            26  \n",
       "3   eis G-10A                 T           293  \n",
       "4   eis C-12T                 A           332  \n",
       "5   eis C-14T                 A           181  \n",
       "6   eis G-37T                 A           285  "
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNP_annotation_DF_filtered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the number of times an *eis* promoter KAN res mutation occurs with an *rrs* AG res mutation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "eis_prom_rrs_co_occur_count = pd.DataFrame(index = ['rrs A1401G', 'rrs C1402T', 'rrs G1484T'], columns = ['eis G-10A', 'eis C-12T', 'eis C-14T', 'eis G-37T'])\n",
    "\n",
    "for rrs_SNV_index in [0,1,2]:\n",
    "    \n",
    "    #boolean filter that is TRUE for isolates that have the rrs SNV\n",
    "    rrs_SNV_bool = SNP_genotypes_array_filtered[rrs_SNV_index,:] == base_code_dict[SNP_annotation_DF_filtered.loc[rrs_SNV_index,'resistance_allele']]\n",
    "    \n",
    "    for eis_prom_SNV_index in [3,4,5,6]:\n",
    "        \n",
    "        #boolean filter that is TRUE for isolates that have the eis prom SNV\n",
    "        eis_prom_SNV_bool = SNP_genotypes_array_filtered[eis_prom_SNV_index,:] == base_code_dict[SNP_annotation_DF_filtered.loc[eis_prom_SNV_index,'resistance_allele']]\n",
    "        \n",
    "        #count the number of isolates that have the eis prom SNV and DO NOT have the rrs SNV\n",
    "        num_isolates_both_SNVs = np.sum([(rrs_SNV_i and eis_prom_SNV_i) for rrs_SNV_i, eis_prom_SNV_i in zip(rrs_SNV_bool, eis_prom_SNV_bool)])\n",
    "        \n",
    "        #store in DataFrame\n",
    "        eis_prom_rrs_co_occur_count.loc[SNP_annotation_DF_filtered.loc[rrs_SNV_index,'variant'],SNP_annotation_DF_filtered.loc[eis_prom_SNV_index,'variant']] = num_isolates_both_SNVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "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>eis G-10A</th>\n",
       "      <th>eis C-12T</th>\n",
       "      <th>eis C-14T</th>\n",
       "      <th>eis G-37T</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>rrs A1401G</th>\n",
       "      <td>2</td>\n",
       "      <td>23</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs C1402T</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs G1484T</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           eis G-10A eis C-12T eis C-14T eis G-37T\n",
       "rrs A1401G         2        23         2         2\n",
       "rrs C1402T         0         0         0         0\n",
       "rrs G1484T         1         1         0         0"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eis_prom_rrs_co_occur_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add in row/column that indicates the total number of isolates have each allele"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "eis_prom_rrs_co_occur_count.loc['# isolates eis prom SNV',:] = [293,332,181,285] # eis promoter mutations\n",
    "eis_prom_rrs_co_occur_count.loc[:,'# isolates rrs SNV'] = [1452,18,26,0] # rrs mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "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>eis G-10A</th>\n",
       "      <th>eis C-12T</th>\n",
       "      <th>eis C-14T</th>\n",
       "      <th>eis G-37T</th>\n",
       "      <th># isolates rrs SNV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>rrs A1401G</th>\n",
       "      <td>2</td>\n",
       "      <td>23</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs C1402T</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs G1484T</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th># isolates eis prom SNV</th>\n",
       "      <td>293</td>\n",
       "      <td>332</td>\n",
       "      <td>181</td>\n",
       "      <td>285</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         eis G-10A  eis C-12T  eis C-14T  eis G-37T  \\\n",
       "rrs A1401G                       2         23          2          2   \n",
       "rrs C1402T                       0          0          0          0   \n",
       "rrs G1484T                       1          1          0          0   \n",
       "# isolates eis prom SNV        293        332        181        285   \n",
       "\n",
       "                         # isolates rrs SNV  \n",
       "rrs A1401G                             1452  \n",
       "rrs C1402T                               18  \n",
       "rrs G1484T                               26  \n",
       "# isolates eis prom SNV                   0  "
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eis_prom_rrs_co_occur_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add in row that indicates the total number of isolates with each of the four classical eis promoter mutations (i.e. c-14t, c-37t, g-10a, and c-12t) that **occur without** one of the three rrs mutations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [],
   "source": [
    "eis_prom_rrs_co_occur_count.loc['# isolates eis prom SNV and no rrs SNV',:] = [293-3,332-24,181-2,285-2,0] # eis promoter mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "eis_prom_rrs_co_occur_count = eis_prom_rrs_co_occur_count.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "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>eis G-10A</th>\n",
       "      <th>eis C-12T</th>\n",
       "      <th>eis C-14T</th>\n",
       "      <th>eis G-37T</th>\n",
       "      <th># isolates rrs SNV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>rrs A1401G</th>\n",
       "      <td>2</td>\n",
       "      <td>23</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs C1402T</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rrs G1484T</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th># isolates eis prom SNV</th>\n",
       "      <td>293</td>\n",
       "      <td>332</td>\n",
       "      <td>181</td>\n",
       "      <td>285</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th># isolates eis prom SNV and no rrs SNV</th>\n",
       "      <td>290</td>\n",
       "      <td>308</td>\n",
       "      <td>179</td>\n",
       "      <td>283</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                        eis G-10A  eis C-12T  eis C-14T  \\\n",
       "rrs A1401G                                      2         23          2   \n",
       "rrs C1402T                                      0          0          0   \n",
       "rrs G1484T                                      1          1          0   \n",
       "# isolates eis prom SNV                       293        332        181   \n",
       "# isolates eis prom SNV and no rrs SNV        290        308        179   \n",
       "\n",
       "                                        eis G-37T  # isolates rrs SNV  \n",
       "rrs A1401G                                      2                1452  \n",
       "rrs C1402T                                      0                  18  \n",
       "rrs G1484T                                      0                  26  \n",
       "# isolates eis prom SNV                       285                   0  \n",
       "# isolates eis prom SNV and no rrs SNV        283                   0  "
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eis_prom_rrs_co_occur_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "eis_prom_rrs_co_occur_count.to_excel('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/CSV files/eis_prom_SNV_rrs_SNV_co_occur_count.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
}