(B) Create Gene Categories.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook takes all of the genes from H37Rv and places them into different categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "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",
    "\n",
    "import Bio\n",
    "from Bio.Alphabet import IUPAC\n",
    "from Bio.Blast.Applications import NcbiblastnCommandline\n",
    "from Bio.Blast import NCBIXML\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.SeqFeature import SeqFeature, FeatureLocation\n",
    "from Bio import pairwise2\n",
    "from Bio import SeqIO\n",
    "from Bio.Graphics import GenomeDiagram\n",
    "from Bio.SeqUtils import GC\n",
    "\n",
    "from Bio.Align.Applications import MuscleCommandline\n",
    "from StringIO import StringIO\n",
    "from Bio import AlignIO\n",
    "from Bio.Align import AlignInfo\n",
    "from Bio.Seq import MutableSeq\n",
    "\n",
    "import networkx as nx\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "############ Filter Out tRNAs from H37Rv annotation ############\n",
    "\n",
    "#H37Rv REFERENCE ANNOTATION\n",
    "H37Rv_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt', '\\t').set_index('name')\n",
    "\n",
    "#filter out annotation for t-RNAs and other misc RNAs\n",
    "RNAs_filter = []\n",
    "for CDS_i in H37Rv_annotation.index:\n",
    "    if ('tRNA' in H37Rv_annotation.loc[CDS_i , 'description'].split(' ')) or ('Rvn' in CDS_i):\n",
    "        RNAs_filter.append(True)\n",
    "        \n",
    "    else:\n",
    "        RNAs_filter.append(False)\n",
    "\n",
    "non_RNAs_filter = [not RNA for RNA in RNAs_filter]\n",
    "\n",
    "H37Rv_annotation_no_RNA = H37Rv_annotation[non_RNAs_filter]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Categories from Comas et. al. 2010 (use to get essential, family-proteins & mobile genetic elements)\n",
    "gene_categories_from_comas_2010 = pd.read_csv( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories_from_comas_2010.csv' , sep = ',' , names = ['gene_category'])\n",
    "\n",
    "############ Family Proteins ############\n",
    "family_protein_genes = []\n",
    "#Family Proteins from Comas et. al. 2010\n",
    "for i in gene_categories_from_comas_2010.index:\n",
    "    \n",
    "    gene_category_list = gene_categories_from_comas_2010.loc[i].values[0].split(' ')\n",
    "    gene = gene_category_list[0]\n",
    "        \n",
    "    #look at PE, PE-PGRS & PPE family proteins\n",
    "    if ('PE' in gene_category_list) or ('PE-PGRS' in gene_category_list) or ('PPE' in gene_category_list):\n",
    "        family_protein_genes.append(gene)\n",
    "        \n",
    "############ Antibiotic Resistance Genes (from Farhat website & Maha & Table with AR associated regions) ############\n",
    "antibiotic_resistance_genes = ['Rv0005','Rv0006','Rv0341','Rv0342','Rv0343','Rv0667','Rv0668','Rv0682','Rv1207','Rv1484','Rv1630','Rv1694','Rv1908c','Rv2043c','Rv2245','Rv2428','Rv2447c','Rv2754c','Rv2763c','Rv2764c','Rv3423c','Rv3608c','Rv3793','Rv3794','Rv3795','Rv3806c','Rv3854c','Rv3919c']\n",
    "        \n",
    "############ Antigens (Not Including Family Proteins or Antibiotic Resistance Genes) ############\n",
    "#Antigens from IEDB epitopes\n",
    "epitope_peptide_seqs_genomic_coords = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/epitope_peptide_sequences/iedb_filtered_epitope_peptides_seqs_and_genomic_coords_05_23_18.csv' ,  sep = ',').set_index('Unnamed: 0')\n",
    "antigen_genes = list( set(epitope_peptide_seqs_genomic_coords.loc[:, 'H37Rv Homolog']) - set(family_protein_genes) - set(antibiotic_resistance_genes))\n",
    "\n",
    "############ Essential & Mobile Genetic Elements ############\n",
    "essential_genes = []\n",
    "mobile_genetic_element_genes = []\n",
    "\n",
    "#iterate through gene categories from comas et. al. 2010\n",
    "for i in gene_categories_from_comas_2010.index:\n",
    "    \n",
    "    gene_category_list = gene_categories_from_comas_2010.loc[i].values[0].split(' ')\n",
    "    gene = gene_category_list[0]\n",
    "    \n",
    "    if ('essential' in gene_category_list) and (gene not in antigen_genes) and (gene not in antibiotic_resistance_genes):\n",
    "        essential_genes.append(gene)\n",
    "    \n",
    "    if ('transposase' in gene_category_list) or ('integrase' in gene_category_list) or ('phage' in gene_category_list) or ('prophage' in gene_category_list) or ('IS' in gene_category_list) or ('IS1533' in gene_category_list) or ('prophage' in gene_category_list) or ('phiRV1' in gene_category_list) or ('phiRv1' in gene_category_list) or ('phiRv2' in gene_category_list) or ('putative' in gene_category_list) and (gene not in antibiotic_resistance_genes):\n",
    "        mobile_genetic_element_genes.append(gene)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of Family-Protein genes also flagged as Antigen genes (kept in seperate categories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "89"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len( set(epitope_peptide_seqs_genomic_coords.loc[:, 'H37Rv Homolog']).intersection(set(family_protein_genes)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "############ Construct into DataFrame ############\n",
    "gene_categories = pd.DataFrame(index = H37Rv_annotation_no_RNA.index , columns = ['gene_id', 'Gene_Category'])\n",
    "\n",
    "for gene_id in gene_categories.index:\n",
    "    \n",
    "    if gene_id in mobile_genetic_element_genes:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'Mobile Genetic Element']\n",
    "    \n",
    "    elif gene_id in family_protein_genes:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'PE/PPE']\n",
    "    \n",
    "    elif gene_id in antigen_genes:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'Antigen']\n",
    "        \n",
    "    elif gene_id in antibiotic_resistance_genes:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'Antibiotic Resistance']\n",
    "    \n",
    "    elif gene_id in essential_genes:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'Essential']\n",
    "        \n",
    "    else:\n",
    "        gene_categories.loc[gene_id,:] = [gene_id , 'Non-Essential']"
   ]
  },
  {
   "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>gene_id</th>\n",
       "      <th>Gene_Category</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Rv0001</th>\n",
       "      <td>Rv0001</td>\n",
       "      <td>Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rv0002</th>\n",
       "      <td>Rv0002</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rv0003</th>\n",
       "      <td>Rv0003</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rv0004</th>\n",
       "      <td>Rv0004</td>\n",
       "      <td>Non-Essential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rv0005</th>\n",
       "      <td>Rv0005</td>\n",
       "      <td>Antibiotic Resistance</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       gene_id          Gene_Category\n",
       "name                                 \n",
       "Rv0001  Rv0001              Essential\n",
       "Rv0002  Rv0002          Non-Essential\n",
       "Rv0003  Rv0003          Non-Essential\n",
       "Rv0004  Rv0004          Non-Essential\n",
       "Rv0005  Rv0005  Antibiotic Resistance"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_categories.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3994, 2)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(gene_categories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "682"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'Essential')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2752"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'Non-Essential')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "257"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'Antigen')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "167"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'PE/PPE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'Antibiotic Resistance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "108"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(gene_categories.Gene_Category == 'Mobile Genetic Element') #all of these CDS regions will be dropped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "############ Save as CSV file ############\n",
    "gene_categories.to_csv( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories.csv' , sep = ',')"
   ]
  }
 ],
 "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
}