(C) Phylogeny Construction for sub-lineage 2.2.1.1.1.i3 cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook was created to prepare the data and construct a phylogeny for the sub-lineage 2.2.1.1.1.i3 isolate cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import itertools\n",
    "import gzip\n",
    "import vcf\n",
    "from slurmpy import Slurm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### t-SNE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import t-SNE embeddings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_SNE_coords = np.loadtxt('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/pairwise_distance_matrix/tb_output/tSNE_embeddings_from_pairwise_SNP_dist.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 15.77039433, -31.88985825],\n",
       "       [ 47.15923691, -36.95353317],\n",
       "       [ 12.14623451,  71.26150513],\n",
       "       ...,\n",
       "       [-15.47960186,  38.40433884],\n",
       "       [ 16.71702766, -43.72228622],\n",
       "       [ 68.94006348,  35.38974762]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_SNE_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(t_SNE_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##################################################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [1] Interactive Session to prepare data for IQTree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##################################################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Different isolates groups: '1', '2', '3', '4A', '4B', '4C', '5', '6'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_lineage = '2.2.1.1.1.i3_cluster'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load in the isolate annotation DF for Genotypes Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load isolate annotation file (columns of Genotype Matrix)\n",
    "isolate_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/homoplasy_project/rolling_DB_scrape/Genotypes_Filtered_2/genotypes_isolate_annotation.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lineage_1</th>\n",
       "      <th>lineage_2</th>\n",
       "      <th>lineage_3</th>\n",
       "      <th>lineage_4</th>\n",
       "      <th>lineage_5</th>\n",
       "      <th>lineage_6</th>\n",
       "      <th>lineage_7</th>\n",
       "      <th>lineage_8</th>\n",
       "      <th>lineage_9</th>\n",
       "      <th>lineage_10</th>\n",
       "      <th>lineage_11</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3558733</td>\n",
       "      <td>4.2.1.2.1.1.i3.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03648641</td>\n",
       "      <td>4.2.1.2.2.1.1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN03647419</td>\n",
       "      <td>3.1.1.i1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3671418</td>\n",
       "      <td>4.2.1.2.1.1.i1</td>\n",
       "      <td>4B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN07659096</td>\n",
       "      <td>1.1.1.2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  lineage_1 lineage_2 lineage_3 lineage_4 lineage_5 lineage_6 lineage_7  \\\n",
       "0         4         2         1         2         1         1        i3   \n",
       "1         4         2         1         2         2         1         1   \n",
       "2         3         1         1        i1       NaN       NaN       NaN   \n",
       "3         4         2         1         2         1         1        i1   \n",
       "4         1         1         1         2       NaN       NaN       NaN   \n",
       "\n",
       "  lineage_8 lineage_9 lineage_10 lineage_11    isolate_ID      lineage_call  \\\n",
       "0         1       NaN        NaN        NaN  SAMEA3558733  4.2.1.2.1.1.i3.1   \n",
       "1       NaN       NaN        NaN        NaN  SAMN03648641     4.2.1.2.2.1.1   \n",
       "2       NaN       NaN        NaN        NaN  SAMN03647419          3.1.1.i1   \n",
       "3       NaN       NaN        NaN        NaN  SAMEA3671418    4.2.1.2.1.1.i1   \n",
       "4       NaN       NaN        NaN        NaN  SAMN07659096           1.1.1.2   \n",
       "\n",
       "  group  \n",
       "0    4B  \n",
       "1    4B  \n",
       "2     3  \n",
       "3    4B  \n",
       "4     1  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolate_annotation_DF.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31428, 14)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolate_annotation_DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a text file with a list of all of the isolates that belong to a particular lineage & filter based off of t-SNE coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lineage 2.2.1.1.1.i3 (SUBSET this sublineage to cluster of isolates with -2 <= t-SNE_1 <= 8 , 26 <= t-SNE_2 <= 38)\n",
    "isolate_sublineage_filter = [isolate_i_lineage_call[0:12] == '2.2.1.1.1.i3' for isolate_i_lineage_call in isolate_annotation_DF.lineage_call] # n = 1887\n",
    "isolate_cluster_filter = [( -2 <= tSNE_1 <= 8 ) and ( 26 <= tSNE_2 <= 38 ) for tSNE_1, tSNE_2 in zip(t_SNE_coords[: , 0], t_SNE_coords[: , 1])] #n = 438\n",
    "\n",
    "#create filter based off of isolates belonging to cluster within sublineage\n",
    "isolate_sublineage_cluster_filter = [(isolate_in_sublineage and isolate_in_cluster) for isolate_in_sublineage, isolate_in_cluster in zip(isolate_sublineage_filter, isolate_cluster_filter)] #n = 438"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolates_belonging_to_sub_lineage = isolate_annotation_DF[isolate_sublineage_cluster_filter]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lineage_1</th>\n",
       "      <th>lineage_2</th>\n",
       "      <th>lineage_3</th>\n",
       "      <th>lineage_4</th>\n",
       "      <th>lineage_5</th>\n",
       "      <th>lineage_6</th>\n",
       "      <th>lineage_7</th>\n",
       "      <th>lineage_8</th>\n",
       "      <th>lineage_9</th>\n",
       "      <th>lineage_10</th>\n",
       "      <th>lineage_11</th>\n",
       "      <th>isolate_ID</th>\n",
       "      <th>lineage_call</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>132</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA1015996</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>134</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA1404001</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>219</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMN08708725</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>269</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA3443144</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>338</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>i3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SAMEA1403810</td>\n",
       "      <td>2.2.1.1.1.i3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    lineage_1 lineage_2 lineage_3 lineage_4 lineage_5 lineage_6 lineage_7  \\\n",
       "132         2         2         1         1         1        i3       NaN   \n",
       "134         2         2         1         1         1        i3       NaN   \n",
       "219         2         2         1         1         1        i3       NaN   \n",
       "269         2         2         1         1         1        i3       NaN   \n",
       "338         2         2         1         1         1        i3       NaN   \n",
       "\n",
       "    lineage_8 lineage_9 lineage_10 lineage_11    isolate_ID  lineage_call  \\\n",
       "132       NaN       NaN        NaN        NaN  SAMEA1015996  2.2.1.1.1.i3   \n",
       "134       NaN       NaN        NaN        NaN  SAMEA1404001  2.2.1.1.1.i3   \n",
       "219       NaN       NaN        NaN        NaN  SAMN08708725  2.2.1.1.1.i3   \n",
       "269       NaN       NaN        NaN        NaN  SAMEA3443144  2.2.1.1.1.i3   \n",
       "338       NaN       NaN        NaN        NaN  SAMEA1403810  2.2.1.1.1.i3   \n",
       "\n",
       "    group  \n",
       "132     2  \n",
       "134     2  \n",
       "219     2  \n",
       "269     2  \n",
       "338     2  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isolates_belonging_to_sub_lineage.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(438, 14)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(isolates_belonging_to_sub_lineage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_IDs_for_sub_lineage = [isolate_ID for isolate_ID in list(isolates_belonging_to_sub_lineage.isolate_ID)]\n",
    "\n",
    "#save each element as a new row in a text file\n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_' + sub_lineage + '.txt', 'w') as f:\n",
    "    for item in isolate_IDs_for_sub_lineage:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a text file with a list of the additional 12 eis C-14T mutants with AG MICs that belong to a particular lineage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolate_IDs_for_sub_lineage_extra_strains = ['622-19','IT1070','IT233','IT634','IT77','IT947']\n",
    "\n",
    "#save each element as a new row in a text file\n",
    "with open('/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_' + '2.2.1.1.1.i3_cluster_extra_strains' + '.txt', 'w') as f:\n",
    "    for item in isolate_IDs_for_sub_lineage_extra_strains:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### [interactive session] BASH function to prepare data for iqtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prep_data_for_iqtree(){\n",
    "    \n",
    "# Load BCFtools\n",
    "module load bcftools/1.9\n",
    "    \n",
    "## The default ulimit is 1024. I should increase it. I think on o2 the hard limit is\n",
    "ulimit -n 50000  #'ulimit -Hn' to check the hard limit, do this for lineage 4 and find login node with high ulimit capacity\n",
    "    \n",
    "# I initialize the variables\n",
    "local OUTDIR=$1\n",
    "local TAG=$2\n",
    "local LIST_ISOL=$3\n",
    "local LIST_ISOL_EXTRA=$4\n",
    "\n",
    "echo \"OUTDIR: ${1}\"\n",
    "echo \"TAG: ${2}\"\n",
    "echo \"LIST_ISOL: ${3}\"\n",
    "echo \"LIST_ISOL_EXTRA: ${4}\"\n",
    "    \n",
    "# I create the output directory\n",
    "rm -rf ${OUTDIR}\n",
    "mkdir -p ${OUTDIR}\n",
    "    \n",
    "# Change the current directory to the output directory\n",
    "cd ${OUTDIR}\n",
    "\n",
    "# I clean the directory with the bcf data\n",
    "rm -rf bcf/\n",
    "mkdir -p bcf/\n",
    "\n",
    "# I generate the bcf files, change the headers and index them\n",
    "for i in `cat ${LIST_ISOL}`;do\n",
    "bcftools view /n/data1/hms/dbmi/farhat/rollingDB/genomic_data/${i}/pilon/${i}.vcf --types snps -O b -o ${i}_snps.bcf\n",
    "bcftools index ${i}_snps.bcf\n",
    "bcftools view ${i}_snps.bcf --regions-file /home/lf61/lf61/mic_assemblies/40-full-analysis/lin-sp-var-10k/results/data_cleaning/regions_to_include_final.bed -o ${i}.bcf -O b;\n",
    "bcftools reheader -s <(echo \"${i}\") ${i}.bcf -o bcf/${i}_renamed.bcf;\n",
    "bcftools index bcf/${i}_renamed.bcf;\n",
    "done\n",
    "    \n",
    "# I generate the bcf files, change the headers and index them\n",
    "for i in `cat ${LIST_ISOL_EXTRA}`;do\n",
    "bcftools view /n/data1/hms/dbmi/farhat/lfreschi/repos/megapipe/megapipe_snakemake/results_mmpR_eis/${i}/pilon/${i}.vcf --types snps -O b -o ${i}_snps.bcf\n",
    "bcftools index ${i}_snps.bcf\n",
    "bcftools view ${i}_snps.bcf --regions-file /home/lf61/lf61/mic_assemblies/40-full-analysis/lin-sp-var-10k/results/data_cleaning/regions_to_include_final.bed -o ${i}.bcf -O b;\n",
    "bcftools reheader -s <(echo \"${i}\") ${i}.bcf -o bcf/${i}_renamed.bcf;\n",
    "bcftools index bcf/${i}_renamed.bcf;\n",
    "done\n",
    "\n",
    "# I add canettii\n",
    "bcftools view /n/data1/hms/dbmi/farhat/lfreschi/mic_assemblies/40-full-analysis/lin-sp-var-10k/results/vcf_canettii/canettii.vcf --types snps -O b -o canettii_snps.bcf\n",
    "bcftools index canettii_snps.bcf\n",
    "bcftools view canettii_snps.bcf --regions-file /home/lf61/lf61/mic_assemblies/40-full-analysis/lin-sp-var-10k/results/data_cleaning/regions_to_include_final.bed -o canettii.bcf -O b;\n",
    "bcftools reheader -s <(echo \"canettii\") canettii.bcf -o bcf/canettii_renamed.bcf;\n",
    "bcftools index bcf/canettii_renamed.bcf;\n",
    "\n",
    "# cleaning\n",
    "rm -rf *.bcf\n",
    "rm -rf *.bcf.csi\n",
    "\n",
    "# I merge the bcf files\n",
    "ls bcf/*|grep \"bcf$\" > list_bcfs.txt\n",
    "bcftools merge -l list_bcfs.txt -o ${OUTDIR}/${TAG}.vcf -O v -0\n",
    "\n",
    "## I use vcf2phylip (from luca's home directory)\n",
    "/home/lf61/sw/vcf2phylip/1.5/vcf2phylip.py -i ${OUTDIR}/${TAG}.vcf -f\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get variables for data preperation function & for iqtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUTDIR = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/tree_output_files/phylogeny_lineage_' + sub_lineage +  '/'\n",
    "TAG = 'lineage_' + sub_lineage\n",
    "LIST_ISOL = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_' + sub_lineage + '.txt'\n",
    "LIST_ISOL_EXTRA = '/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_2.2.1.1.1.i3_cluster_extra_strains.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/tree_output_files/phylogeny_lineage_2.2.1.1.1.i3_cluster/'"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "OUTDIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'lineage_2.2.1.1.1.i3_cluster'"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TAG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_2.2.1.1.1.i3_cluster.txt'"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LIST_ISOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_2.2.1.1.1.i3_cluster_extra_strains.txt'"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LIST_ISOL_EXTRA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### [interactive session] Call function above in an interactive session to prepare the data for iqtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prep_data_for_iqtree /n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/tree_output_files/phylogeny_lineage_2.2.1.1.1.i3_cluster/ lineage_2.2.1.1.1.i3_cluster /n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_2.2.1.1.1.i3_cluster.txt /n/data1/hms/dbmi/farhat/Roger/mmpR_BDQ_mutant_project/phylogenies/sublineage_isolate_lists/isolate_list_lineage_2.2.1.1.1.i3_cluster_extra_strains.txt\n"
     ]
    }
   ],
   "source": [
    "print 'prep_data_for_iqtree' + ' ' + OUTDIR + ' ' + TAG + ' ' + LIST_ISOL + ' ' + LIST_ISOL_EXTRA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##################################################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [2] Submit job to O2 to create trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##################################################################################################################################################################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use code below to submit a iqtree job (uses ModelFinder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iq_L2.2.1.1.1.i3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "submitted: Submitted batch job 26950900\n"
     ]
    }
   ],
   "source": [
    "#iqtree_job_p1 = 'VIRTUAL_ENV_DISABLE_PROMPT=true source activate iqtree_virtualenv'\n",
    "iqtree_job_p1 = 'set +eu\\nVIRTUAL_ENV_DISABLE_PROMPT=true source activate iqtree_virtualenv\\nset -eu' # to fix \"/home/rv76/anaconda2/etc/profile.d/conda.sh: line 55: PS1: unbound variable\" ERROR\n",
    "iqtree_job_p2 = 'time iqtree -s {0}{1}.min4.fasta -mem 7G -pre {0}tree_{1}_iqtree_FINAL -m MFP -mset GTR -bb 1000 -alrt 1000 -nt AUTO -ntmax 8 -redo'.format(OUTDIR, TAG)\n",
    "\n",
    "iqtree_job = iqtree_job_p1 + '\\n' + iqtree_job_p2\n",
    "\n",
    "#directory where you want output + error files\n",
    "os.chdir(OUTDIR)\n",
    "\n",
    "job_name = 'iq_L' + TAG.split('_')[1]  \n",
    "\n",
    "s = Slurm(job_name , {'partition':'priority' , 'N':'1' , 'c':'8' , 't':'0-12:00:00' , 'mem-per-cpu':'8G' , 'mail-type':'ALL' , 'mail-user':'roger_vargas@g.harvard.edu' , 'o':'out_run_iqtree_{}.txt'.format(TAG)})\n",
    "\n",
    "#submits the job\n",
    "job_id = s.run(iqtree_job)\n",
    "\n",
    "print job_name"
   ]
  }
 ],
 "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
}