Dion_al_2020_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sqlite3\n",
    "import os\n",
    "from collections import Counter, defaultdict\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import axes\n",
    "import matplotlib.colors\n",
    "import numpy as np\n",
    "from ete3 import NCBITaxa\n",
    "ncbi = NCBITaxa()\n",
    "import pandas as pd\n",
    "import statistics\n",
    "import random\n",
    "import scipy.stats as stats\n",
    "import plotly.graph_objs as go\n",
    "from plotly.offline import iplot, init_notebook_mode\n",
    "init_notebook_mode()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Streamlining CRISPR spacer-based bacterial host predictions to decipher the viral dark matter (Dion et al. 2020)</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goal: optimize phage host predictions using CRISPR spacers, determine sensitivity and specificity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Build spacers database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Download all bacterial genomes on NCBI\n",
    "\n",
    "#### Each bacterial genome has a ftp address which corresponds to the location of the file on the NCBI server. The list of all ftp addresses for all bacterial genomes is found in this summary document in the column ftp_path:\n",
    "\n",
    "#### ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt\n",
    "\n",
    "#### As of March 23rd 2020, there were 580 384 bacterial genomes deposited on NCBI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Identify CRISPR loci with CRISPRDetect\n",
    "\n",
    "#### We used CRISPRDetect to identify CRISPR loci as it is the most sensitive and accurate software. More information on CRISPRDetect is available in the related publication and Github page of the authors:\n",
    "\n",
    "#### https://pubmed.ncbi.nlm.nih.gov/27184979/\n",
    "\n",
    "#### https://github.com/ambarishbiswas/CRISPRDetect_2.2\n",
    "\n",
    "#### We added the option '-array_quality_score_cutoff 3' as recommended by the authors because we used fasta files. We also added the option '-tmp_dir $SLURM_TMPDIR' because CRISPRDetect creates temporary files and does not delete them which increases very quickly the number of files. Three output files are generated by CRISPRDetect and we used the gff file to parse the information. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Composition and diversity of the spacer database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SpacersDB import CrisprOpenDB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn = sqlite3.connect('SpacersDB/CrisprOpenDB.sqlite')\n",
    "c = conn.cursor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of crispr positive bacteria\n",
    "c.execute('''select DISTINCT ORG.ASSEMBLY_ACCESSION \\\n",
    "from SPACER_TABLE ST, ORGANISM ORG \\\n",
    "where ST.GENEBANK_ID=ORG.GENEBANK_ID''')\n",
    "len(c.fetchall())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of crispr positive bacteria\n",
    "c.execute('''select DISTINCT ORG.ASSEMBLY_ACCESSION \\\n",
    "from SPACER_TABLE ST, ORGANISM ORG \\\n",
    "where ST.GENEBANK_ID=ORG.GENEBANK_ID''')\n",
    "accession = [x[0] for x in c.fetchall()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of spacers\n",
    "c.execute('''select COUNT(*) from SPACER_TABLE''')\n",
    "c.fetchall()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shorter spacer\n",
    "c.execute('''select MIN(SPACER_LENGTH) from SPACER_TABLE''')\n",
    "c.fetchall()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#longer spacer\n",
    "c.execute('''select MAX(SPACER_LENGTH) from SPACER_TABLE''')\n",
    "c.fetchall()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of spacers between 28 and 43 nt\n",
    "c.execute('''select COUNT(*) from SPACER_TABLE where SPACER_LENGTH BETWEEN 28 AND 43''')\n",
    "c.fetchall()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "11674395/11767782 * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db_explorer = CrisprOpenDB.CrisprOpenDB(os.path.join(\"SpacersDB\", \"CrisprOpenDB.sqlite\"))\n",
    "\n",
    "df = pd.read_sql_query(\"select ST.SPACER, ORG.GENUS \\\n",
    "from SPACER_TABLE ST, ORGANISM ORG \\\n",
    "where ST.GENEBANK_ID=ORG.GENEBANK_ID\", db_explorer._connection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spacers_diversity = df[df.GENUS != 'Unknown']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spacers_diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_spacers = Counter(spacers_diversity.GENUS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_spacers.most_common(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "8325687/588364"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_total = ['$\\it{Salmonella}$',\n",
    "          '$\\it{Listeria}$',\n",
    "          '$\\it{Escherichia}$',\n",
    "          '$\\it{Clostridioides}$',\n",
    "          '$\\it{Mycobacterium}$',\n",
    "          '$\\it{Campylobacter}$',\n",
    "          '$\\it{Pseudomonas}$',\n",
    "          '$\\it{Klebsiella}$',\n",
    "          '$\\it{Acinetobacter}$',\n",
    "          '$\\it{Streptococcus}$']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total number of spacers for Salmonella is 14 times greater than the second most abundant\n",
    "# genus, so we will change the figure a little bit to have more readable numbers and add\n",
    "# an axis break manually on Illustrator.\n",
    "x_total = [800000, 588364, 368069, 271784, 199856, 147421, 121125, 79980, 78313, 58885]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6, 5.5))\n",
    "\n",
    "y_total_pos = np.arange(len(y_total))\n",
    "\n",
    "barlist = ax.barh(y_total_pos, x_total, align = 'center')\n",
    "barlist[0].set_color('#d55e00')\n",
    "barlist[1].set_color('#949494')\n",
    "barlist[2].set_color('#949494')\n",
    "barlist[3].set_color('#949494')\n",
    "barlist[4].set_color('#949494')\n",
    "barlist[5].set_color('#cc78bc')\n",
    "barlist[6].set_color('#0173b2')\n",
    "barlist[7].set_color('#949494')\n",
    "barlist[8].set_color('#029e73')\n",
    "barlist[9].set_color('#de8f05')\n",
    "ax.set_yticks(y_total_pos)\n",
    "ax.set_yticklabels(y_total, size = 20)\n",
    "ax.set_xticklabels(['0', '100', '200', '300', '400', '500','600', '', '8000'], size = 18)\n",
    "\n",
    "ax.invert_yaxis()\n",
    "ax.set_xlabel('Total spacer sequences (thousands)', size = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_spacers = Counter(spacers_diversity.drop_duplicates(subset=['SPACER','GENUS']).GENUS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(spacers_diversity.SPACER.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_spacers.most_common(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_unique = ['$\\it{Salmonella}$',\n",
    "          '$\\it{Clostridium}$',\n",
    "          '$\\it{Streptomyces}$',\n",
    "          '$\\it{Lactobacillus}$',\n",
    "          '$\\it{Acinetobacter}$',\n",
    "          '$\\it{Streptococcus}$',\n",
    "          '$\\it{Bifidobacterium}$',\n",
    "          '$\\it{Pseudomonas}$',\n",
    "          '$\\it{Campylobacter}$',\n",
    "          '$\\it{Corynebacterium}$']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6, 5.5))\n",
    "\n",
    "y_unique_pos = np.arange(len(y_unique))\n",
    "x_unique = [x[1] for x in unique_spacers.most_common(10)]\n",
    "\n",
    "barlist = ax.barh(y_unique_pos, x_unique, align = 'center')\n",
    "barlist[0].set_color('#d55e00')\n",
    "barlist[1].set_color('#949494')\n",
    "barlist[2].set_color('#949494')\n",
    "barlist[3].set_color('#949494')\n",
    "barlist[4].set_color('#029e73')\n",
    "barlist[5].set_color('#de8f05')\n",
    "barlist[6].set_color('#949494')\n",
    "barlist[7].set_color('#0173b2')\n",
    "barlist[8].set_color('#cc78bc')\n",
    "barlist[9].set_color('#949494')\n",
    "ax.set_yticks(y_unique_pos)\n",
    "ax.set_yticklabels(y_unique, size = 20)\n",
    "ax.set_xticklabels(['0', '10', '20', '30', '40', '50','60', '70'], size = 18)\n",
    "\n",
    "\n",
    "ax.invert_yaxis()\n",
    "ax.set_xlabel('Unique spacer sequences (thousands)', size = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "under95_x = []\n",
    "under95_y = []\n",
    "above95_x = []\n",
    "above95_y = []\n",
    "ratio_unique_total_distribution = []\n",
    "\n",
    "for genus, unique in unique_spacers.items() :\n",
    "    ratio = unique/total_spacers[genus]\n",
    "    ratio_unique_total_distribution.append(ratio)\n",
    "    if ratio >= 0.95 :\n",
    "        above95_x.append(total_spacers[genus])\n",
    "        above95_y.append(unique)\n",
    "    else :\n",
    "        under95_x.append(total_spacers[genus])\n",
    "        under95_y.append(unique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (8,9))\n",
    "\n",
    "plt.scatter(x = under95_x, \n",
    "            y = under95_y, \n",
    "            alpha = 0.5, s = 12, label = '< 0.95', c = 'gainsboro')\n",
    "\n",
    "plt.scatter(x = above95_x, \n",
    "            y = above95_y, \n",
    "            alpha = 0.5, s = 12, label = '≥ 0.95', c = 'mediumvioletred')\n",
    "\n",
    "legend = plt.legend(title = 'unique/total\\n   spacers', fontsize = 22, numpoints = 1, bbox_to_anchor=(0.4,1))\n",
    "legend.legendHandles[0]._sizes = [80]\n",
    "legend.legendHandles[1]._sizes = [80]\n",
    "plt.setp(legend.get_title(),fontsize=22)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Total number of spacers', size = 22)\n",
    "plt.ylabel('Number of unique spacers', size = 22)\n",
    "plt.ylim(1,10000000)\n",
    "plt.xticks(size = 22)\n",
    "plt.yticks(size = 22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(ratio_unique_total_distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_spacers['Salmonella']/total_spacers['Salmonella']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(i >= 0.95 for i in ratio_unique_total_distribution)/len(ratio_unique_total_distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (8,3))\n",
    "plt.hist(ratio_unique_total_distribution, bins = 20, color = 'k')\n",
    "plt.xlabel('Ratio unique/total spacers', size = 20)\n",
    "plt.ylabel('Frequency', size = 20)\n",
    "plt.xticks(size = 20)\n",
    "plt.yticks(size = 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Evaluate recall and precision using phages with known hosts\n",
    "\n",
    "### 2.1 Download phage genomes\n",
    "\n",
    "#### We used the NCBIVirus webtool to obtain the list of all phages and download their genomes.\n",
    "\n",
    "#### https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide\n",
    "\n",
    "#### In the left column, we selected Virus > Bacteriophages, all taxids and Nucleotide Completeness > complete.  \n",
    "\n",
    "#### We then separeted the results by Sequence Type (Genbank or RefSeq) because almost all RefSeq are duplicates from Genbank and we wanted to avoid duplicates as much as possible. We downloaded the two resulting tables in csv format, keeping \"Accession\", \"Species\", \"Genus\", \"Family\", \"Length\", \"Host\" and \"GenBank_Title\". We also downloaded all the genomes in FASTA format.\n",
    "\n",
    "### 2.2 Clean up dataset\n",
    "\n",
    "#### There are some phages in this list for which the bacterial host is unknown. This is the case mostly for phages sequenced from metagenomics. These phages will be identified as \"organism from which the phage was sequenced\" + phage (e.g. bees, lynx, chimpanzees, etc.). These phages need to be removed from our dataset because we need to know the host of each phage that will be used to test our prediction method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phage_table = pd.read_csv('genbank_20200423.csv')\n",
    "\n",
    "# For some phages, the column 'host' is empty. For the majority of them, \n",
    "# it's because the host is really unknown, but for some others it looks\n",
    "# like a mistake because in the phage name the host name is clearly there.\n",
    "\n",
    "phage_table = phage_table.fillna(0)\n",
    "\n",
    "for index, row in phage_table.iterrows():\n",
    "    if row['Host'] == 0 :\n",
    "        deduced_host = row['Species'].split(' ')[0]\n",
    "        phage_table.loc[index:index, 'Host'] = deduced_host\n",
    "\n",
    "all_host = phage_table.Host.unique().tolist()\n",
    "host_taxid = ncbi.get_name_translator(all_host)\n",
    "bacteria_host = {}\n",
    "\n",
    "for i in all_host :\n",
    "    if i not in host_taxid.keys() :\n",
    "        bacteria_host[i] = 0\n",
    "    else :\n",
    "        if 2 in ncbi.get_lineage(host_taxid[i][0]) :\n",
    "            bacteria_host[i] = 1\n",
    "        else :\n",
    "            bacteria_host[i] = 0\n",
    "\n",
    "for index, row in phage_table.iterrows():\n",
    "    if bacteria_host[row['Host']] == 0 :\n",
    "        phage_table = phage_table.drop([index])\n",
    "        \n",
    "# Next thing to do is remove \"biological\" duplicates, i.e. phages that have a different\n",
    "# Genbank_Title but are the result of laboratory evolution experiments. We need to remove\n",
    "# these duplicates because they would cause an unjustified overrepresentation of some sequences.\n",
    "# Most of these \"biological\" duplicates have either \"mutant\" or \"clone\" in their name, so\n",
    "# we'll simply filter these out.\n",
    "\n",
    "duplicates = ['mutant', 'Mutant', 'clone','isolate','evolved']\n",
    "\n",
    "for i in duplicates :\n",
    "    phage_table = phage_table[~phage_table.GenBank_Title.str.contains(i)]\n",
    "            \n",
    "# Last step is to add phages categorized in RefSeq and remove duplicates. We'll do the same\n",
    "# cleaning up steps then compare their \"Genbank_Title\" and keep only RefSeq phages that\n",
    "# are NOT already in our table (this should be a very small number)\n",
    "\n",
    "refseq_table = pd.read_csv('refseq_20200423.csv')\n",
    "\n",
    "refseq_table = refseq_table.fillna(0)\n",
    "\n",
    "for index, row in refseq_table.iterrows():\n",
    "    if row['Host'] == 0 :\n",
    "        deduced_host = row['Species'].split(' ')[0]\n",
    "        refseq_table.loc[index:index, 'Host'] = deduced_host\n",
    "\n",
    "all_host = refseq_table.Host.unique().tolist()\n",
    "host_taxid = ncbi.get_name_translator(all_host)\n",
    "bacteria_host = {}\n",
    "\n",
    "for i in all_host :\n",
    "    if i not in host_taxid.keys() :\n",
    "        bacteria_host[i] = 0\n",
    "    else :\n",
    "        if 2 in ncbi.get_lineage(host_taxid[i][0]) :\n",
    "            bacteria_host[i] = 1\n",
    "        else :\n",
    "            bacteria_host[i] = 0\n",
    "\n",
    "for index, row in refseq_table.iterrows():\n",
    "    if bacteria_host[row['Host']] == 0 :\n",
    "        refseq_table = refseq_table.drop([index])\n",
    "        \n",
    "\n",
    "title_genbank = phage_table.GenBank_Title.tolist()\n",
    "\n",
    "for index, row in refseq_table.iterrows():\n",
    "    if row['GenBank_Title'] in title_genbank :\n",
    "        refseq_table = refseq_table.drop([index])\n",
    "        \n",
    "phage_table = pd.concat([phage_table, refseq_table])\n",
    "phage_table.set_index('Accession', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Search for homology between phage genomes and the spacers database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Perform predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_spacer_table():\n",
    "    db_explorer = CrisprOpenDB.CrisprOpenDB(os.path.join(\"SpacersDB\", \"CrisprOpenDB.sqlite\"))\n",
    "    df = pd.read_sql_query(\"select ST.SPACER_ID, ST.GENEBANK_ID, ORG.GENUS, ORG.FAMILY, ORG.TORDER, ST.SPACER_LENGTH, SAL.COUNT_SPACER, ST.POSITION_INSIDE_LOCUS  \\\n",
    "        from ORGANISM ORG, SPACER_TABLE ST, SPACER_ARRAY_LENGTH SAL \\\n",
    "        where ST.GENEBANK_ID=ORG.GENEBANK_ID and ST.GENEBANK_ID=SAL.GENEBANK_ID and ST.NUMERO_LOCUS=SAL.NUMERO_LOCUS\", db_explorer._connection)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spacer_table = load_spacer_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_blastn_results(blastn_file):\n",
    "    columns = [\"Query\", \"SPACER_ID\", \"identity\", \"alignement_length\", \n",
    "    \"mismatch\", \"gap\", \"q_start\", \"q_end\", \"s_start\", \"s_end\", \"e_value\", \"score\"]\n",
    "    blastn_result_table = pd.read_csv(blastn_file, names=columns, sep=\"\\t\")\n",
    "    return blastn_result_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# blastn alignment results between phage genomes and spacer database\n",
    "# file available upon request\n",
    "# enter path to file here\n",
    "path_to_alignment_result = ''\n",
    "alignement_results = load_blastn_results(path_to_alignment_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_table = pd.merge(alignement_results, spacer_table, on=\"SPACER_ID\", how=\"left\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findHost(summary_table, criteria, limit):\n",
    "    \n",
    "    result_table = summary_table.copy()\n",
    "    result_table = result_table.dropna()\n",
    "    result_table = result_table[(result_table.gap <= 0) & (result_table.GENUS != 'Unknown')]\n",
    "\n",
    "    if criteria == 'mismatch' :\n",
    "        result_table.loc[:,\"true_num_mismatch\"] = (result_table.loc[:,\"SPACER_LENGTH\"] - result_table.loc[:,\"alignement_length\"]) + result_table.loc[:,\"mismatch\"]\n",
    "        result_table  = result_table[result_table.true_num_mismatch <= limit]\n",
    "    \n",
    "    elif criteria == 'e_value' :\n",
    "        result_table  = result_table[result_table.e_value <= limit]\n",
    "        \n",
    "    else :\n",
    "        print('wrong criteria')\n",
    "        return\n",
    "    \n",
    "    if result_table.empty:\n",
    "        return\n",
    "\n",
    "    result_table[\"mean_position\"] = np.array((result_table[\"q_start\"] + result_table[\"q_end\"]) / 2, dtype=int)\n",
    "\n",
    "\n",
    "    #Criteria 1: If only one genus, it is the host.\n",
    "\n",
    "    if len(set(result_table[\"GENUS\"])) == 1:\n",
    "        prediction = {\"PREDICTION\": result_table[\"GENUS\"].iloc[0]}\n",
    "        return(prediction)\n",
    "    \n",
    "    else:\n",
    "        no_pred = {\"level_1\": result_table[\"GENUS\"].unique().tolist()}\n",
    "\n",
    "\n",
    "    #Criteria 2: Number of different position on phage genome (use mean(start, end) to check positions) (ex: MF153391)\n",
    "\n",
    "    genus = np.array(result_table[\"GENUS\"])\n",
    "    position = np.array(result_table[\"mean_position\"])\n",
    "    sets_to_count = list(set([(j, position[i]) for i,j in enumerate(genus)]))\n",
    "    counted_genus = Counter([i[0] for i in sets_to_count])\n",
    "\n",
    "    most_commons_genus = counted_genus.most_common()\n",
    "\n",
    "    if most_commons_genus[0][1] != most_commons_genus[1][1]: # If count is not equal, we found host\n",
    "        prediction = {\"PREDICTION\": most_commons_genus[0][0]}\n",
    "        return({**no_pred, **prediction})\n",
    "    genuses_to_keep = []\n",
    "    for i, g in enumerate(most_commons_genus):\n",
    "        if most_commons_genus[i][1] == most_commons_genus[0][1]:\n",
    "            genuses_to_keep.append(g[0])\n",
    "\n",
    "    result_table = result_table[result_table[\"GENUS\"].isin(genuses_to_keep)]\n",
    "    no_pred['level_2'] = genuses_to_keep\n",
    "\n",
    "    #Criteria 3: If 2 is equal, find relative position most in 5' (MF158036)\n",
    "    \n",
    "    five_prime_relative_position = (np.array(result_table[\"POSITION_INSIDE_LOCUS\"], dtype=float) -1) / (np.array(result_table[\"COUNT_SPACER\"]) - 1)\n",
    "    result_table[\"five_prime_relative_position\"] = five_prime_relative_position\n",
    "    result_table.sort_values(by=\"five_prime_relative_position\", inplace=True)\n",
    "    result_table.reset_index(inplace=True, drop=True)\n",
    "\n",
    "    sub_section = result_table[result_table[\"five_prime_relative_position\"] == result_table[\"five_prime_relative_position\"][0]]\n",
    "\n",
    "    if len(sub_section[\"GENUS\"].unique().tolist()) == 1:\n",
    "        prediction = {\"PREDICTION\": sub_section[\"GENUS\"].iloc[0]}\n",
    "        return({**no_pred, **prediction})\n",
    "    else:\n",
    "        result_table = sub_section\n",
    "        no_pred['level_3'] = sub_section[\"GENUS\"].unique().tolist()\n",
    "\n",
    "    #Criteria 4: Last common ancester (does not return a genus)\n",
    "\n",
    "    if len(set(result_table[\"FAMILY\"])) == 1:\n",
    "        prediction = {\"PREDICTION\": result_table[\"FAMILY\"].iloc[0]}\n",
    "        return({**no_pred, **prediction})\n",
    "    elif len(set(result_table[\"TORDER\"])) == 1:\n",
    "        prediction = {\"PREDICTION\": result_table[\"TORDER\"].iloc[0]}\n",
    "        return({**no_pred, **prediction})\n",
    "    else:\n",
    "        return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(criteria, value) :\n",
    "    pred_dict = {}\n",
    "    for phage in phage_table.index:\n",
    "        df = summary_table[summary_table[\"Query\"] == phage]\n",
    "        pred = findHost(df, criteria, value)\n",
    "        if pred != None :\n",
    "            pred_dict[phage] = pred\n",
    "    return pred_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_lca_rank(host, pred_host):\n",
    "    \n",
    "    translator = ncbi.get_name_translator([host, pred_host])\n",
    "    \n",
    "    lineage_host = ncbi.get_lineage(translator[host][0])\n",
    "    lineage_pred_host = ncbi.get_lineage(translator[pred_host][0])\n",
    "    \n",
    "    lca = [i for i, j in zip(lineage_host, lineage_pred_host) if i == j]\n",
    "    rank = ncbi.get_rank([lca[-1]])[lca[-1]]\n",
    "    return(rank)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def specificity_sensitivity(pred_host_dict, total_number_of_phages):\n",
    "    \n",
    "    specificity_by_rank = defaultdict(int)\n",
    "\n",
    "    true_positive = 0\n",
    "    false_positive = 0\n",
    "    \n",
    "    for phage in pred_host_dict.keys():\n",
    "        prediction = pred_host_dict[phage]['PREDICTION']\n",
    "        real_host = phage_table.loc[phage]['Host']              \n",
    "        lca_rank = find_lca_rank(real_host, prediction)\n",
    "        specificity_by_rank[lca_rank] += 1\n",
    "\n",
    "    number_of_predictions = len(pred_host_dict)\n",
    "    sensitivity = number_of_predictions/total_number_of_phages * 100\n",
    "    \n",
    "    specificity_by_rank['family'] = specificity_by_rank['family'] + specificity_by_rank['genus']\n",
    "    specificity_by_rank['order'] = specificity_by_rank['order'] + specificity_by_rank['family']\n",
    "    \n",
    "    specificity_by_rank['genus'] = specificity_by_rank['genus'] / number_of_predictions * 100\n",
    "    specificity_by_rank['family'] = specificity_by_rank['family'] / number_of_predictions * 100\n",
    "    specificity_by_rank['order'] = specificity_by_rank['order'] / number_of_predictions * 100\n",
    "\n",
    "    \n",
    "    return specificity_by_rank, sensitivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_0_mismatch = predict('mismatch',0)\n",
    "pred_1_mismatch = predict('mismatch',1)\n",
    "pred_2_mismatch = predict('mismatch',2)\n",
    "pred_3_mismatch = predict('mismatch',3)\n",
    "pred_4_mismatch = predict('mismatch',4)\n",
    "pred_5_mismatch = predict('mismatch',5)\n",
    "pred_6_mismatch = predict('mismatch',6)\n",
    "pred_7_mismatch = predict('mismatch',7)\n",
    "pred_8_mismatch = predict('mismatch',8)\n",
    "pred_9_mismatch = predict('mismatch',9)\n",
    "pred_10_mismatch = predict('mismatch',10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specificity_by_rank_0 = specificity_sensitivity(pred_0_mismatch, 9484)\n",
    "specificity_by_rank_1 = specificity_sensitivity(pred_1_mismatch, 9484)\n",
    "specificity_by_rank_2 = specificity_sensitivity(pred_2_mismatch, 9484)\n",
    "specificity_by_rank_3 = specificity_sensitivity(pred_3_mismatch, 9484)\n",
    "specificity_by_rank_4 = specificity_sensitivity(pred_4_mismatch, 9484)\n",
    "specificity_by_rank_5 = specificity_sensitivity(pred_5_mismatch, 9484)\n",
    "specificity_by_rank_6 = specificity_sensitivity(pred_6_mismatch, 9484)\n",
    "specificity_by_rank_7 = specificity_sensitivity(pred_7_mismatch, 9484)\n",
    "specificity_by_rank_8 = specificity_sensitivity(pred_8_mismatch, 9484)\n",
    "specificity_by_rank_9 = specificity_sensitivity(pred_9_mismatch, 9484)\n",
    "specificity_by_rank_10 = specificity_sensitivity(pred_10_mismatch, 9484)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-deep')\n",
    "genus = []\n",
    "family = []\n",
    "order = []\n",
    "\n",
    "sensitivity_recap = []\n",
    "\n",
    "for i in [specificity_by_rank_0, specificity_by_rank_1, specificity_by_rank_2,\n",
    "         specificity_by_rank_3, specificity_by_rank_4, specificity_by_rank_5,\n",
    "         specificity_by_rank_6, specificity_by_rank_7, specificity_by_rank_8,\n",
    "         specificity_by_rank_9, specificity_by_rank_10] :\n",
    "    genus.append(i[0]['genus'])\n",
    "    family.append(i[0]['family'])\n",
    "    order.append(i[0]['order'])\n",
    "\n",
    "    sensitivity_recap.append(i[1])\n",
    "    \n",
    "x_axis = [0,1,2,3,4,5,6,7,8,9,10]\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "\n",
    "plt.plot(x_axis, sensitivity_recap, linestyle = '--', color = 'black', marker = 'o')\n",
    "plt.plot(x_axis, sensitivity_recap, alpha = 0)\n",
    "plt.plot(x_axis, order, marker='v')\n",
    "plt.plot(x_axis, family, marker='o')\n",
    "plt.plot(x_axis, genus, marker='D')\n",
    "\n",
    "plt.xticks(x_axis, fontsize=20)\n",
    "plt.xlabel(\"Maximal number of mismatches\", fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.ylim([0,100])\n",
    "plt.ylabel(\"%\", fontsize=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_000000001_evalue = predict('e_value', 0.000000001)\n",
    "pred_00000001_evalue = predict('e_value', 0.00000001)\n",
    "pred_0000001_evalue = predict('e_value', 0.0000001)\n",
    "pred_000001_evalue = predict('e_value', 0.000001)\n",
    "pred_00001_evalue = predict('e_value', 0.00001)\n",
    "pred_0001_evalue = predict('e_value', 0.0001)\n",
    "pred_001_evalue = predict('e_value', 0.001)\n",
    "pred_01_evalue = predict('e_value', 0.01)\n",
    "pred_1_evalue = predict('e_value', 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specificity_by_rank_evalue_000000001 = specificity_sensitivity(pred_000000001_evalue, 9484)\n",
    "specificity_by_rank_evalue_00000001 = specificity_sensitivity(pred_00000001_evalue, 9484)\n",
    "specificity_by_rank_evalue_0000001 = specificity_sensitivity(pred_0000001_evalue, 9484)\n",
    "specificity_by_rank_evalue_000001 = specificity_sensitivity(pred_000001_evalue, 9484)\n",
    "specificity_by_rank_evalue_00001 = specificity_sensitivity(pred_00001_evalue, 9484)\n",
    "specificity_by_rank_evalue_0001 = specificity_sensitivity(pred_0001_evalue, 9484)\n",
    "specificity_by_rank_evalue_001 = specificity_sensitivity(pred_001_evalue, 9484)\n",
    "specificity_by_rank_evalue_01 = specificity_sensitivity(pred_01_evalue, 9484)\n",
    "specificity_by_rank_evalue_1 = specificity_sensitivity(pred_1_evalue, 9484)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-deep')\n",
    "genus = []\n",
    "family = []\n",
    "order = []\n",
    "\n",
    "sensitivity_recap = []\n",
    "\n",
    "for i in [specificity_by_rank_evalue_000000001, specificity_by_rank_evalue_00000001, specificity_by_rank_evalue_0000001,\n",
    "         specificity_by_rank_evalue_000001, specificity_by_rank_evalue_00001, specificity_by_rank_evalue_0001,\n",
    "         specificity_by_rank_evalue_001, specificity_by_rank_evalue_01, specificity_by_rank_evalue_1] :\n",
    "    genus.append(i[0]['genus'])\n",
    "    family.append(i[0]['family'])\n",
    "    order.append(i[0]['order'])\n",
    "\n",
    "    sensitivity_recap.append(i[1])\n",
    "    \n",
    "x_axis = [0.000000001, 0.00000001, 0.0000001, 0.000001, 0.00001,\n",
    "          0.0001, 0.001, 0.01, 0.1]\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(x_axis, sensitivity_recap, linestyle = '--', color = 'black', marker = 'o')\n",
    "plt.plot(x_axis, sensitivity_recap, alpha = 0)\n",
    "plt.plot(x_axis, order, marker='v')\n",
    "plt.plot(x_axis, family, marker='o')\n",
    "plt.plot(x_axis, genus, marker='D')\n",
    "\n",
    "plt.xticks(x_axis, fontsize=20)\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"Maximum e-value\", fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.ylim([0,100])\n",
    "plt.legend([\"Recall\", 'Precision:',\"Order\",\"Family\",\"Genus\"], bbox_to_anchor=(1.05, 1), prop={'size': 20})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Data exploration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_pred_2_mismatch = []\n",
    "\n",
    "for phage in pred_2_mismatch.keys():\n",
    "    final_pred_2_mismatch.append(pred_2_mismatch[phage]['PREDICTION'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Contribution of each filter in identifying one final predicted host\n",
    "\n",
    "no_filter = []\n",
    "crit_1 = []\n",
    "crit_2 = []\n",
    "crit_3 = []\n",
    "crit_4 = []\n",
    "\n",
    "for phage in pred_2_mismatch.keys():\n",
    "    len_no_filter = len(summary_table[summary_table[\"Query\"] == phage]['GENUS'].unique().tolist())\n",
    "    if len_no_filter > 10 :\n",
    "        len_no_filter = 10\n",
    "    no_filter.append(len_no_filter)\n",
    "\n",
    "    if len(pred_2_mismatch[phage].keys()) == 1 :\n",
    "        crit_1.append(1)\n",
    "        crit_2.append(1)\n",
    "        crit_3.append(1)\n",
    "        crit_4.append(1)\n",
    "    if len(pred_2_mismatch[phage].keys()) == 2 :\n",
    "        len_1 = len(pred_2_mismatch[phage]['level_1'])\n",
    "        if len_1 > 10 :\n",
    "            len_1 = 10\n",
    "        crit_1.append(len_1)\n",
    "        crit_2.append(1)\n",
    "        crit_3.append(1)\n",
    "        crit_4.append(1)\n",
    "    if len(pred_2_mismatch[phage].keys()) == 3 :\n",
    "        len_1 = len(pred_2_mismatch[phage]['level_1'])\n",
    "        if len_1 > 10 :\n",
    "            len_1 = 10\n",
    "        len_2 = len(pred_2_mismatch[phage]['level_2'])\n",
    "        if len_2 > 10 :\n",
    "            len_2 = 10\n",
    "        crit_1.append(len_1)\n",
    "        crit_2.append(len_2)\n",
    "        crit_3.append(1)\n",
    "        crit_4.append(1)\n",
    "    if len(pred_2_mismatch[phage].keys()) == 4 :\n",
    "        len_1 = len(pred_2_mismatch[phage]['level_1'])\n",
    "        if len_1 > 10 :\n",
    "            len_1 = 10\n",
    "        len_2 = len(pred_2_mismatch[phage]['level_2'])\n",
    "        if len_2 > 10 :\n",
    "            len_2 = 10\n",
    "        len_3 = len(pred_2_mismatch[phage]['level_3'])\n",
    "        if len_3 > 10 :\n",
    "            len_3 = 10\n",
    "        crit_1.append(len_1)\n",
    "        crit_2.append(len_2)\n",
    "        crit_3.append(len_3)\n",
    "        crit_4.append(1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crit_2.count(1)/len(pred_2_mismatch)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(crit_3.count(1)-crit_2.count(1))/len(pred_2_mismatch)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(crit_4.count(1)-crit_3.count(1))/len(pred_2_mismatch)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-deep')\n",
    "plt.figure(figsize = (8,6))\n",
    "plt.hist([no_filter,crit_1,crit_2, crit_3, crit_4], \n",
    "         bins = range(1, 11 + 1, 1))\n",
    "plt.legend(labels = ['No filter',\n",
    "                     'After 2 mismatches',\n",
    "                     'After highest number of targets',\n",
    "                     'After 5\\' end neighbouring spacers',\n",
    "                    'After last common ancestor'],\n",
    "          prop={'size': 14})\n",
    "plt.xlabel('Number of predicted hosts (genus level)', size = 20)\n",
    "plt.ylabel('Number of phages', size = 20)\n",
    "plt.xticks(np.arange(1+1/2, 11, 1), [\"1\",\"2\",\"3\",\"4\",\"5\",\"6\",\"7\",\"8\",\"9\",\"10+\"], size = 20)\n",
    "plt.yticks(size = 20)\n",
    "#plt.savefig(\"contributution_per_filter_20200812.png\", dpi =300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is there a difference in the accuracy of prediction depending on the filter used to identify the predicted host?\n",
    "accuracy_by_filter = defaultdict(list)\n",
    "\n",
    "for phage, pred in pred_2_mismatch.items() :\n",
    "    level = len(pred)\n",
    "    pred_host = pred_2_mismatch[phage]['PREDICTION']\n",
    "    real_host = phage_table.loc[phage]['Host']\n",
    "    rank = find_lca_rank(real_host, pred_host)\n",
    "    if rank == 'genus' :\n",
    "        accuracy_by_filter[level].append('Y')\n",
    "    else :\n",
    "        accuracy_by_filter[level].append('N')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_by_filter[1].count('Y')/len(accuracy_by_filter[1])*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_by_filter[2].count('Y')/len(accuracy_by_filter[2])*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_by_filter[3].count('Y')/len(accuracy_by_filter[3])*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-deep')\n",
    "plt.figure(figsize=(8,6))\n",
    "N = len(accuracy_by_filter)\n",
    "yes = [accuracy_by_filter[1].count('Y'), \n",
    "       accuracy_by_filter[2].count('Y'), \n",
    "       accuracy_by_filter[3].count('Y'), \n",
    "       accuracy_by_filter[4].count('Y')]\n",
    "no = [accuracy_by_filter[1].count('N'), \n",
    "      accuracy_by_filter[2].count('N'), \n",
    "      accuracy_by_filter[3].count('N'), \n",
    "      accuracy_by_filter[4].count('N')]\n",
    "ind = np.arange(N)\n",
    "width = 0.35\n",
    "\n",
    "p1 = plt.bar(ind, yes, width, color = 'mediumseagreen')\n",
    "p2 = plt.bar(ind, no, width, bottom=yes, color = 'indianred')\n",
    "\n",
    "plt.ylabel('Number of phages', size = 20)\n",
    "plt.yticks(size = 20)\n",
    "plt.xticks(ind, ('Filter 1', 'Filter 2', 'Filter 3', 'Filter 4'), size = 20)\n",
    "legend = plt.legend(title = 'Prediction', labels = ['Accurate', 'Wrong'], prop={'size': 14})\n",
    "legend.get_title().set_fontsize('14')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_accuracy = defaultdict(list)\n",
    "for phage in pred_2_mismatch.keys():\n",
    "        prediction = pred_2_mismatch[phage]['PREDICTION']\n",
    "        real_host = phage_table.loc[phage]['Host'].split(' ')[0]      \n",
    "        lca_rank = find_lca_rank(real_host, prediction)\n",
    "        if lca_rank in ['genus','family'] :\n",
    "            if lca_rank == 'genus' :\n",
    "                accuracy_genus = 'Y'\n",
    "                accuracy_family = 'Y'\n",
    "            else :\n",
    "                accuracy_genus = 'N'\n",
    "                accuracy_family = 'Y'\n",
    "        else :\n",
    "            accuracy_genus = 'N'\n",
    "            accuracy_family = 'N'\n",
    "            \n",
    "        pred_accuracy[phage].extend((prediction, real_host, accuracy_genus, accuracy_family))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_accuracy_table = pd.DataFrame.from_dict(pred_accuracy, orient='index', columns = ['prediction','real_host','accuracy_genus','accuracy_family'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_accuracy_table = pred_accuracy_table.merge(phage_table[['Length','Family']], right_index=True, left_index=True, how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_accuracy_table = pred_accuracy_table[pred_accuracy_table.Family != 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_accuracy_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(set(phage_table.Host))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(set(pred_accuracy_table.real_host))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_by_real_host = defaultdict(list)\n",
    "host_total = Counter([x.split(' ')[0] for x in phage_table.Host.tolist()])\n",
    "\n",
    "for i in host_total.keys() :\n",
    "    number_of_pred = len(pred_accuracy_table[pred_accuracy_table.real_host == i].accuracy_genus)\n",
    "    sensitivity = round(number_of_pred / host_total[i] * 100,2)\n",
    "    \n",
    "    if sensitivity == 0 :\n",
    "        specificity_genus = 0\n",
    "        specificity_family = 0\n",
    "    \n",
    "    else :  \n",
    "        try :\n",
    "            yes_genus = pred_accuracy_table[pred_accuracy_table.real_host == i].accuracy_genus.tolist().count('Y')\n",
    "            specificity_genus = round(yes_genus/number_of_pred*100,2)\n",
    "        except :\n",
    "            specificity_genus = 0\n",
    "\n",
    "        try :\n",
    "            yes_family = pred_accuracy_table[pred_accuracy_table.real_host == i].accuracy_family.tolist().count('Y')\n",
    "            specificity_family = round(yes_family/number_of_pred*100,2)\n",
    "        except :\n",
    "            yes_family = 0\n",
    "    \n",
    "    accuracy_by_real_host[i].extend((number_of_pred,\n",
    "                                     host_total[i],\n",
    "                                     sensitivity,\n",
    "                                     specificity_genus, \n",
    "                                     specificity_family))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(accuracy_by_real_host)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_by_host_table = pd.DataFrame.from_dict(accuracy_by_real_host, orient = 'index', \n",
    "                       columns=['number_pred','number_phages', 'sensitivity','accuracy_genus','accuracy_family'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_unique_spacers = []\n",
    "\n",
    "for index, row in accuracy_by_host_table.iterrows() :\n",
    "    host = index\n",
    "    number_unique_spacers.append(unique_spacers[host])\n",
    "    \n",
    "accuracy_by_host_table.loc[:,'number_unique_spacers'] = number_unique_spacers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,3))\n",
    "plt.axvline(x = 50, color = 'grey', linewidth=1)\n",
    "plt.axhline(y = 50, color = 'grey',linewidth=1)\n",
    "\n",
    "plt.scatter(x = accuracy_by_host_table[(accuracy_by_host_table.number_phages > 10) ].sensitivity,\n",
    "            y = accuracy_by_host_table[(accuracy_by_host_table.number_phages > 10)].accuracy_genus,\n",
    "           c = accuracy_by_host_table[(accuracy_by_host_table.number_phages > 10) ].number_unique_spacers,\n",
    "           norm = matplotlib.colors.LogNorm(),\n",
    "           s = 10)\n",
    "plt.colorbar(label = 'Number of host unique spacers')\n",
    "plt.xlabel('Recall (%)', size = 10)\n",
    "plt.ylabel('Precision (%)', size = 10)\n",
    "plt.xticks(size = 10)\n",
    "plt.yticks(size = 10)\n",
    "\n",
    "\n",
    "plt.annotate('3', (95,55), size = 10, color = 'k')\n",
    "plt.annotate('2', (1,55), size = 10, color = 'k')\n",
    "plt.annotate('1', (1,40), size = 10, color = 'k')\n",
    "plt.annotate('4', (95,40), size = 10, color = 'k')\n",
    "\n",
    "plt.annotate('$\\it{Synechococcus}$\\n    $\\it{Prochlorococcus}$', (-4,12), size = 5, color = 'k')\n",
    "plt.plot([0,-2],[2,14], color='k', linewidth=0.5)\n",
    "plt.plot([0,2],[2,10], color='k', linewidth=0.5)\n",
    "\n",
    "plt.annotate('$\\it{Cellulophaga}$', (17,5), size = 5, color = 'k')\n",
    "plt.plot([19.5,19.5],[2,4], color='k', linewidth=0.5)\n",
    "\n",
    "plt.annotate('$\\it{Clostridioides}$', (78,92), size = 5, color = 'k')\n",
    "plt.annotate('$\\it{Prevotella}$', (90,88.5), size = 5, color = 'k')\n",
    "plt.plot([100,102],[98,91], color='k', linewidth=0.5)\n",
    "plt.plot([100,98],[98,94.5], color='k', linewidth=0.5)\n",
    "\n",
    "plt.annotate('$\\it{Enterococcus}$', (62,98), size = 5, color = 'k')\n",
    "plt.plot([92.5,84],[97.5,98.5], color='k', linewidth=0.5)\n",
    "\n",
    "\n",
    "#plt.savefig(\"sensitivity_specificity_by_host_20201007.png\", dpi =300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.spearmanr(accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].sensitivity,\n",
    "            accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].accuracy_genus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.spearmanr(accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].sensitivity,\n",
    "            accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].number_unique_spacers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.spearmanr(accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].accuracy_genus,\n",
    "            accuracy_by_host_table[accuracy_by_host_table.number_phages > 10].number_unique_spacers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn-deep')\n",
    "N = 6\n",
    "sensitivity = [accuracy_by_real_host['Synechococcus'][2]+0.3, \n",
    "               accuracy_by_real_host['Prochlorococcus'][2]+0.3, \n",
    "               accuracy_by_real_host['Cellulophaga'][2],\n",
    "               accuracy_by_real_host['Clostridioides'][2], \n",
    "               accuracy_by_real_host['Enterococcus'][2], \n",
    "               accuracy_by_real_host['Prevotella'][2]]\n",
    "\n",
    "specificity = [accuracy_by_real_host['Synechococcus'][3]+0.3, \n",
    "               accuracy_by_real_host['Prochlorococcus'][3]+0.3, \n",
    "               accuracy_by_real_host['Cellulophaga'][3]+0.3,\n",
    "               accuracy_by_real_host['Clostridioides'][3], \n",
    "               accuracy_by_real_host['Enterococcus'][3], \n",
    "               accuracy_by_real_host['Prevotella'][3]]\n",
    "\n",
    "ind = np.arange(N)\n",
    "width = 0.35\n",
    "\n",
    "plt.figure(figsize= (8,6))\n",
    "plt.bar(ind, sensitivity, width, color = 'palevioletred')\n",
    "plt.bar(ind+width, specificity, width, color = 'cornflowerblue')\n",
    "\n",
    "plt.xticks(ind,('$\\it{Synechococcus}$', '$\\it{Prochlorococcus}$', '$\\it{Cellulophaga}$',\n",
    "                '$\\it{Clostridioides}$', '$\\it{Enterococcus}$', '$\\it{Prevotella}$'),\n",
    "           rotation = 20, size = 16)\n",
    "\n",
    "plt.yticks(size = 20)\n",
    "plt.ylim(0,119)\n",
    "\n",
    "legend = plt.legend(labels = ['Recall', 'Precision'], prop={'size': 14},\n",
    "                   loc='upper left', bbox_to_anchor=(0.01,0.9))\n",
    "legend.get_title().set_fontsize('14')\n",
    "\n",
    "plt.axvline(x = 2.7, color = 'grey')\n",
    "plt.annotate('Category 1', (0.5, 110), size = 16, color = 'grey')\n",
    "plt.annotate('Category 3', (3.5, 110), size = 16, color = 'grey')\n",
    "plt.ylabel('%', size = 20)\n",
    "\n",
    "#plt.savefig(\"category_1_and_3_20201009.png\", dpi =300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Statistical analyses\n",
    "\n",
    "#### We will determine if our prediction approach is better than picking a random hit from the alignment result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phage_all_hit_dict = {}\n",
    "result_table = summary_table.dropna()\n",
    "\n",
    "for phage in set(np.array(result_table[\"Query\"])) :\n",
    "    phage_all_hit_dict[phage] = result_table[result_table[\"Query\"] == phage]['GENUS'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Statistical analysis: null model to verify if prediction approach better than random pick among hits\n",
    "\n",
    "random.seed(12345)\n",
    "null_prediction_dict = defaultdict(list)\n",
    "for i in range(1,1001,1):\n",
    "    for phage, hit in phage_all_hit_dict.items():\n",
    "        random_host = random.choice(hit)\n",
    "        real_host = phage_table.loc[phage]['Host']\n",
    "        rank = find_lca_rank(real_host, random_host)\n",
    "        if rank == 'genus' :\n",
    "            null_prediction_dict[i].append('Y')\n",
    "        else :\n",
    "            null_prediction_dict[i].append('N')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "success_rate_null_model = {}\n",
    "for i, success_list in null_prediction_dict.items():\n",
    "    yes = success_list.count('Y')\n",
    "    success_rate_null_model[i] = yes/9884*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(success_rate_null_model.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize= (8,6))\n",
    "plt.hist(success_rate_null_model.values())\n",
    "plt.xlabel('Precision (%)', size = 16)\n",
    "plt.ylabel('Frequency', size = 16)\n",
    "plt.xticks(size = 16)\n",
    "plt.yticks(size = 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Compare performance with other methods\n",
    "\n",
    "### 3.1 Select existing software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WIsH phage list\n",
    "# https://github.com/soedinglab/WIsH/blob/master/benchmark/PhageSequences.lst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Confirm the phage list used for the two other methods are present in our dataset\n",
    "\n",
    "#### WIsH used only phages from RefSeq. In the previous step (2.1) of cleaning up the dataset, we removed most of the RefSeq phage genomes as they are duplicates of Genbank phage genomes. We return to the RefSeq table to retrieve the name for each accession number and search for that name in the final phage dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import Entrez\n",
    "\n",
    "def NCBIEntrez(accession):\n",
    "    Entrez.email = \"moira.dion.1@ulaval.ca\"\n",
    "    handle = Entrez.efetch(db=\"nucleotide\", id=accession, rettype=\"gb\", retmode=\"text\")\n",
    "    for line in handle:\n",
    "        if line.startswith('  ORGANISM'):\n",
    "            line = line.strip().split()\n",
    "            species = ' '.join(line[1:])\n",
    "            return species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refseq_table = pd.read_csv('refseq_20200423.csv')\n",
    "refseq_table.set_index('Accession', inplace=True)\n",
    "refseq_accession = refseq_table.index.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wish_dataset = []\n",
    "path_to_wish_list = ''\n",
    "with open(path_to_wish_list,'r') as f:\n",
    "    for line in f:\n",
    "        wish_dataset.append(line.strip())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in wish_dataset :\n",
    "    if i not in refseq_accession and i not in refseq_archaeal_virus_accession :\n",
    "        print(i + '\\t' + NCBIEntrez(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following lines were manually verified\n",
    "\n",
    "# NC_003525 is Stx2 converting phage I and host is unknown (removed)\n",
    "# NC_005808 refers to AY526908\n",
    "# NC_005809 refers to AY526909\n",
    "# NC_005949 (Vibrio phage Vf12) no longer exists, but shares exactly the same sequence\n",
    "## with NC_005948 (Vibrio phage Vf33). They were deposited at the same time from the\n",
    "## same authors so we will use Vf33 as substitute.\n",
    "# NC_007457 refers to DQ222851\n",
    "# NC_008355 refers to NC_002194\n",
    "# NC_010762 refers to EU568876 \n",
    "# NC_010763 refers to EU676000\n",
    "# NC_011105 refers to EU056923\n",
    "# NC_011284 refers to FJ174694\n",
    "# NC_011286 refers to FJ168660\n",
    "# NC_020836 (Synechococcus phage KBS-M-1A) no longer exists. We tried finding a similar cyanophage \n",
    "## under a different name but we did not find anything. No other phage in our database has\n",
    "## the same genome size (removed)\n",
    "# NC_021296 refers to KC691254\n",
    "# NC_021303 refers to KC661274\n",
    "# NC_021778 refers to JN255163\n",
    "# NC_021859 refers to KC787108\n",
    "# NC_022326 refers to KF279416\n",
    "# NC_022327 refers to KF279413\n",
    "# NC_023596 refers to KF981875\n",
    "# NC_023718 refers to JN175269\n",
    "# NC_027338 refers to KJ578763\n",
    "# NC_027343 refers to KJ578775\n",
    "# NC_027980 refers to JQ267518"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wish_dataset_clean = []\n",
    "\n",
    "for i in wish_dataset :\n",
    "    if i in refseq_accession or i in refseq_archaeal_virus_accession :\n",
    "        wish_dataset_clean.append(i)\n",
    "        \n",
    "wish_dataset_clean.extend(['AY526908','AY526909','NC_005948','DQ222851','NC_002194',\n",
    "                          'EU568876','EU676000','EU056923','FJ174694','FJ168660',\n",
    "                          'KC691254','KC661274','JN255163','KC787108','KF279416',\n",
    "                          'KF279413','KF981875','JN175269','KJ578763','KJ578775','JQ267518'])\n",
    "\n",
    "print(len(wish_dataset), len(wish_dataset_clean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_wish(criteria, value) :\n",
    "    pred_dict = {}\n",
    "    for query in wish_dataset_clean:\n",
    "        df = summary_table[summary_table[\"Query\"] == query]\n",
    "        pred = findHost(df, criteria, value)\n",
    "        if pred != None :\n",
    "            pred_dict[query] = pred\n",
    "    return pred_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wish_pred_2_mismatch = predict_wish('mismatch',2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wish_performance = specificity_sensitivity(wish_pred_2_mismatch, len(wish_dataset_clean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(wish_pred_2_mismatch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wish_performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Proof of concept on already published virome\n",
    "\n",
    "### We are using the virome dataset by Shkoporov et al (2019) Cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downoladed here: \n",
    "# https://figshare.com/articles/The_human_gut_virome_is_highly_diverse_stable_and_individual-specific_/9248864/1\n",
    "\n",
    "# Retrieved from this folder: Dataset_S1/selected_contigs.fa\n",
    "# For a total of 57 721 contigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# blast alignment results between shkoporov virome and spacer database\n",
    "# file available on request\n",
    "# enter path to file here\n",
    "path_to_shkoporov_alignment = ''\n",
    "shkoporov_results = load_blastn_results(path_to_shkoporov_alignment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shkoporov_summary_table = pd.merge(shkoporov_results, spacer_table, on=\"SPACER_ID\", how=\"left\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shkoporov_summary_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shkoporov_vibrant_list = []\n",
    "with open('shkoporov_virome.phages_combined.txt','r') as f:\n",
    "    for line in f:\n",
    "        line = line.strip()\n",
    "        shkoporov_vibrant_list.append(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_shkoporov(criteria, value) :\n",
    "    pred_dict = {}\n",
    "    for phage in shkoporov_vibrant_list:\n",
    "        df = shkoporov_summary_table[shkoporov_summary_table[\"Query\"] == phage]\n",
    "        pred = findHost(df, criteria, value)\n",
    "        if pred != None :\n",
    "            pred_dict[phage] = pred\n",
    "    return pred_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shkoporov_pred_2_mismatch = predict_shkoporov('mismatch', 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shkoporov_final_pred_2_mismatch = {}\n",
    "\n",
    "for phage, pred in shkoporov_pred_2_mismatch.items():\n",
    "    shkoporov_final_pred_2_mismatch[phage] = pred['PREDICTION']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Supplementary table S2 from Shkoporov et al 2019 paper\n",
    "# enter path to file here\n",
    "path_to_table_S2 = ''\n",
    "virome_description = pd.read_excel(path_to_table_S2, skiprows=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "virome_description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "virome_for_comparison = virome_description[(virome_description['Contig name'].isin(shkoporov_final_pred_2_mismatch)) &\n",
    "                  (virome_description['CRISPR-based host'].notna())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_shkop_spacer_pred = defaultdict(list)\n",
    "\n",
    "for index, row in virome_for_comparison.iterrows() :\n",
    "    contig = row['Contig name']\n",
    "    \n",
    "    shkop_pred = row['CRISPR-based host']\n",
    "    compare_shkop_spacer_pred[contig].append(shkop_pred)\n",
    "    \n",
    "    spacer_pred = shkoporov_final_pred_2_mismatch[contig]\n",
    "    compare_shkop_spacer_pred[contig].append(spacer_pred)\n",
    "    \n",
    "    try :\n",
    "        lca = find_lca_rank(shkop_pred, spacer_pred)\n",
    "    except :\n",
    "        lca = 'unknown'\n",
    "    \n",
    "    compare_shkop_spacer_pred[contig].append(lca)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Counter([x[2] for x in compare_shkop_spacer_pred.values()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(shkoporov_vibrant_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(shkoporov_final_pred_2_mismatch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(compare_shkop_spacer_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "virome_description[(virome_description['Contig name'].isin(shkoporov_vibrant_list)) &\n",
    "                  (virome_description['CRISPR-based host'].notna())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(shkoporov_final_pred_2_mismatch)/len(shkoporov_vibrant_list)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1392/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "160+109+105+75+65+64+60+50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_eight = Counter(shkoporov_final_pred_2_mismatch.values()).most_common(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1392-688"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ['<i>' + i[0] for i in top_eight]\n",
    "labels.append('Other')\n",
    "values = [i[1] for i in top_eight]\n",
    "values.append(704)\n",
    "trace_pie_org_names = go.Pie(labels=labels, values=values, sort = False)\n",
    "layout = go.Layout(title=\"Most common hosts species\")\n",
    "fig = go.Figure(data=[trace_pie_org_names])\n",
    "fig.update_layout(legend=go.layout.Legend(x = 1, y = 1, font=dict(size=20)))\n",
    "fig.update_traces(textinfo='value', textfont_size=20)\n",
    "iplot(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clostridium_contigs = []\n",
    "for phage, pred in shkoporov_final_pred_2_mismatch.items() :\n",
    "    if pred == 'Clostridium' :\n",
    "        clostridium_contigs.append(phage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clostridium_table = shkoporov_summary_table[shkoporov_summary_table['Query'].isin(clostridium_contigs)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clostridium_table.loc[:,'mean_position'] = np.array((clostridium_table.loc[:,\"q_start\"] + clostridium_table.loc[:,\"q_end\"]) / 2, dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clostridium_table.loc[:,\"true_num_mismatch\"] = (clostridium_table.loc[:,\"SPACER_LENGTH\"] - clostridium_table.loc[:,\"alignement_length\"]) + clostridium_table.loc[:,\"mismatch\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clostridium_targets = {}\n",
    "\n",
    "for i in clostridium_contigs :\n",
    "    clostridium_targets[i] = Counter(clostridium_table[(clostridium_table['Query'] == i) &\n",
    "                  (clostridium_table['true_num_mismatch'] <= 2) & \n",
    "                  (clostridium_table['GENUS'] == 'Clostridium')]['mean_position'].tolist()).most_common()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_targets(vOTU, genome_size):\n",
    "\n",
    "    x_coord = []\n",
    "    y_coord = []\n",
    "    for target in clostridium_targets[vOTU] :\n",
    "        i = 1\n",
    "        while i <= target[1] :\n",
    "            x_coord.append(target[0])\n",
    "            y_coord.append(0)\n",
    "            i += 1\n",
    "    width = round(genome_size * 0.0009)\n",
    "    size = 20\n",
    "    plt.figure(figsize=(width, 1))\n",
    "    plt.scatter(x_coord, y_coord, s = 100, color = '#00CC96')\n",
    "    plt.xlim(0,genome_size+500)\n",
    "    plt.xticks(size = size)\n",
    "    plt.ylim(-0.1,0.5)\n",
    "    plt.yticks([])\n",
    "    label = 'Position on putative phage contig (bp)'\n",
    "    plt.xlabel(label, size = size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_targets('metaspades_NG-13376_921T3_lib202033_5478_NODE_37_length_22057_cov_6.42255',22057)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}