calculateResistome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculations of the resistome datasets from nanopore, Sanger and PacBio\n",
    "This notebook takes the Sanger and PacBio dataset and BLAST it against the CARD database. Next all the nanopore, Sanger and PacBio datasets are blast'ed against the Sanger data to get a sense of the sequence identity\n",
    "\n",
    "* Requires [blast command line tool](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:1011: UserWarning: Duplicate key in file \"//anaconda/lib/python2.7/site-packages/matplotlib/mpl-data/stylelib/dark_background.mplstyle\", line #9\n",
      "  (fname, cnt))\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from collections import OrderedDict\n",
    "import pandas as pd\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from matplotlib import pyplot as plt\n",
    "import os\n",
    "from Bio import SeqIO\n",
    "from types import *\n",
    "import numpy as np\n",
    "from subprocess import Popen, PIPE\n",
    "import os\n",
    "import logging\n",
    "\n",
    "logger = logging.getLogger('CARDBlast')\n",
    "logging.basicConfig()\n",
    "logger.setLevel(logging.INFO)\n",
    "\n",
    "CORES = 4\n",
    "\n",
    "CARDdb = 'inputFiles/n.fasta.protein.homolog.fasta' #Card database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## run poreFUME\n",
    "The following block will execute poreFUME, be aware that this can take a long time depending on the amount of cores avialable (set with the ```--cores``` flag). The reads are split up by barcode over the cores. Also be sure to have run ```install.sh``` prior to running this notebook. In the ```poreFUME/``` directory you can run ```nosetests``` to check poreFUME is working correctly. Furthermore poreFUME relies on avialable paths, so if you did not add the paths to you ```~./profile``` run ```source env.sh``` to load them in the current shell.  \n",
    "** The data in the ```output/``` directory will be overwritten by running poreFUME with the current ```--overwriteDemux --overwriteNanocorrect --overwriteCARD --overwriteNanopolish``` flags. If you just wish to analyze the data head over to the ```analyzeResistome.ipynb``` which uses the already generated data in ```output/``` **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "HERE = os.getcwd()\n",
    "nanocorrectDir = os.path.join(HERE,'poreFUME','nanocorrect')\n",
    "\n",
    "#run poreFUME (see integrationTest.sh for an short example run on 75 raw reads)\n",
    "!python poreFUME/poreFUME.py inputData/testSet75.fasta inputData/pb_39.fasta --PacBioLegacyBarcode --cores 8 --pathCARD=inputData/n.fasta.protein.homolog.fasta --pathNanocorrect=---abspahgoeshere--/nanocorrect/ --pathRawreads=---abspahgoeshere--/test/data/testSet75 --overwriteNanocorrect --pathNanopolish=---abspahgoeshere--/nanopolish/ --overwriteNanopolish --overwriteDemux --overwriteCARD\n",
    "\n",
    "!python poreFUME/poreFUME.py inputFiles/LejlaControl.2D.min500bp.fasta poreFUME/inputData/pb_39b.fasta  --barcodeEdge=60 --pathNanocorrect={nanocorrectDir} --PacBioLegacyBarcode --cores 8 --pathCARD=inputData/n.fasta.protein.homolog.fasta --pathNanocorrect=---abspahgoeshere--/nanocorrect/ --pathRawreads=---abspahgoeshereToENAFiles- --overwriteNanocorrect --pathNanopolish=---abspahgoeshere--/nanopolish/ --overwriteNanopolish --overwriteDemux --overwriteCARD\n",
    "!python poreFUME/poreFUME.py inputFiles/poreCamp.2D.min500.fasta poreFUME/inputData/pb_39b.fasta  --pathNanocorrect={nanocorrectDir}  --barcodeEdge=120 --PacBioLegacyBarcode --cores 8 --pathCARD=inputData/n.fasta.protein.homolog.fasta --pathNanocorrect=---abspahgoeshere--/nanocorrect/ --pathRawreads=---abspahgoeshereToENAFiles- --overwriteNanocorrect --pathNanopolish=---abspahgoeshere--/nanopolish/ --overwriteNanopolish --overwriteDemux --overwriteCARD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cmdExists(cmd):\n",
    "    return any(\n",
    "        os.access(os.path.join(path, cmd), os.X_OK)\n",
    "        for path in os.environ[\"PATH\"].split(os.pathsep)\n",
    "    )\n",
    "if not cmdExists('blastn'):\n",
    "    raise  IOError('blastn not found in PATH, install from for example https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download')\n",
    "if not cmdExists('makeblastdb'):\n",
    "    raise  IOError('makeblastdb not found in PATH, install from for example https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load in the paths to all the files. \n",
    "In the python file ```paths.py``` a ```OrderedDict``` keeps track of the paths to the FASTA files. It is loaded in using the ```%load``` magic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load paths.py\n",
    "qPathOrdered = OrderedDict([\n",
    "        \n",
    "    ('S_raw','inputFiles/sanger.total.aftertrim.removeCT.min500bp.fasta'),\n",
    "    ('S_ass','inputFiles/SangerAssembledAndNonAssembled.min500bp.withoutCandT.fasta'),\n",
    "\n",
    "\n",
    "    ('PB_raw','inputFiles/PB.Cell1and2.raw.fasta'),\n",
    "    ('PB_ass','inputFiles/PacBio.Correct.C1.fasta'),\n",
    "\n",
    "    ('MI_raw','inputFiles/LejlaControl.2D.min500bp.fasta'),\n",
    "    ('MI_demux'   , 'output/LejlaControl.2D.min500bp.afterBC.fasta'),\n",
    "    ('MI_corr1'   , 'output/LejlaControl.2D.min500bp.afterNC1.fasta'),\n",
    "    ('MI_corr2'   , 'output/LejlaControl.2D.min500bp.afterNC2.fasta'), \n",
    "    ('MI_polish'   , 'output/LejlaControl.2D.min500bp.afterNP.fasta'),\n",
    "    ('MI_polishNoFilter'   , 'output/LejlaControl.2D.min500bp.afterNPnofilter.fasta'), #This file will be generated by calculateResistome.ipynb\n",
    "\n",
    "\n",
    "    ('MIpc_raw', 'inputFiles/poreCamp.2D.min500.fasta'), \n",
    "    ('MIpc_demux' , 'output/poreCamp.2D.min500.afterBC.fasta'),  \n",
    "    ('MIpc_corr1' , 'output/poreCamp.2D.min500.afterNC1.fasta'), \n",
    "    ('MIpc_corr2' , 'output/poreCamp.2D.min500.afterNC2.fasta'),\n",
    "    #('MIpc_polish'   , 'output/poreCamp.2D.min500.afterNP.fasta'),\n",
    "    #('MIpc_polishNoFilter'   , 'output/poreCamp.2D.min500.afterNPnofilter.fasta') #This file will be generated by calculateResistome.ipynb\n",
    "\n",
    "        \n",
    "      \n",
    "    \n",
    "    ])\n",
    "\n",
    "annotationPathOrdered = OrderedDict([\n",
    "        \n",
    "    ('S_raw','output/annotation/sanger.total.aftertrim.removeCT.min500bp.fasta.annotated.csv'), \n",
    "    ('S_ass','output/annotation/SangerAssembledAndNonAssembled.min500bp.withoutCandT.fasta.annotated.csv'),\n",
    "\n",
    "\n",
    "    ('PB_raw','output/annotation/PB.Cell1and2.raw.fasta.annotated.csv'), \n",
    "    ('PB_ass','output/annotation/PacBio.Correct.C1.fasta.annotated.csv'),\n",
    "\n",
    "    ('MI_raw','output/annotation/LejlaControl.2D.min500bp/LejlaControl.2D.min500bp.annotated.csv'), \n",
    "    ('MI_demux'   , 'output/annotation/LejlaControl.2D.min500bp/LejlaControl.2D.min500bp.afterBC.annotated.csv'), \n",
    "    ('MI_corr1'   , 'output/annotation/LejlaControl.2D.min500bp/LejlaControl.2D.min500bp.afterNC1.annotated.csv'), \n",
    "    ('MI_corr2'   , 'output/annotation/LejlaControl.2D.min500bp/LejlaControl.2D.min500bp.afterNC2.annotated.csv'), \n",
    "    ('MI_polishNoFilter'   , 'output/annotation/LejlaControl.2D.min500bp.afterNPnofilter.fasta.annotated.csv'), \n",
    "\n",
    "\n",
    "    ('MIpc_raw', 'output/annotation/poreCamp.2D.min500/poreCamp.2D.min500.annotated.csv'), \n",
    "    ('MIpc_demux' , 'output/annotation/poreCamp.2D.min500/poreCamp.2D.min500.afterBC.annotated.csv'),  \n",
    "    ('MIpc_corr1' , 'output/annotation/poreCamp.2D.min500/poreCamp.2D.min500.afterNC1.annotated.csv'),  \n",
    "    ('MIpc_corr2' , 'output/annotation/poreCamp.2D.min500/poreCamp.2D.min500.afterNC2.annotated.csv'), \n",
    "    #('MIpc_polishNoFilter'   , 'output/annotation/poreCamp.2D.min500.afterNPnofilter.fasta.annotated.csv'), \n",
    "    \n",
    "\n",
    "    \n",
    "    ])\n",
    "\n",
    "\n",
    "readLabels = OrderedDict([\n",
    "        \n",
    "    ('S_raw','Sanger raw'),\n",
    "    ('S_ass','Sanger assembled'),\n",
    "\n",
    "\n",
    "    ('PB_raw','PacBio raw'), \n",
    "    ('PB_ass','PacBio assembled'),\n",
    "\n",
    "\n",
    "    ('MIpc_raw', 'nanopore lib. B:\\n2D reads'), \n",
    "    ('MIpc_demux' , 'nanopore lib. B:\\nafter debarcoding'),\n",
    "    ('MIpc_corr1' , 'nanopore lib. B:\\nafter nanocorrect round 1'), \n",
    "    ('MIpc_corr2' , 'nanopore lib. B:\\nafter nanocorrect round 2'), \n",
    "    #('MIpc_polishNoFilter'   , 'nanopore lib. B:\\nafter nanopolish')  ,\n",
    "\n",
    "    ('MI_raw','nanopore lib. A: 2D reads'), \n",
    "    ('MI_demux'   , 'nanopore lib. A:\\nafter debarcoding'),\n",
    "    ('MI_corr1'   , 'nanopore lib. A:\\nafter nanocorrect round 1'), \n",
    "    ('MI_corr2'   , 'nanopore lib. A:\\nafter nanocorrect round 2')  ,\n",
    "    ('MI_polishNoFilter'   , 'nanopore lib. A:\\nafter nanopolish ')  \n",
    "\n",
    "    \n",
    "    ])\n",
    "\n",
    "readLabelsSample = OrderedDict([\n",
    "        \n",
    "    ('S_raw','Sanger'), \n",
    "    ('S_ass','Sanger'),\n",
    "\n",
    "\n",
    "    ('PB_raw','PacBio'), \n",
    "    ('PB_ass','PacBio'),\n",
    "\n",
    "\n",
    "    ('MIpc_raw'   , 'nanopore lib. B'),\n",
    "    ('MIpc_demux' , 'nanopore lib. B'), \n",
    "    ('MIpc_corr1' , 'nanopore lib. B'), \n",
    "    ('MIpc_corr2' , 'nanopore lib. B'), \n",
    "    #('MIpc_polishNoFilter'   , 'nanopore lib. B') ,\n",
    "    \n",
    "    ('MI_raw'     , 'nanopore lib. A'),\n",
    "    ('MI_demux'   , 'nanopore lib. A'),\n",
    "    ('MI_corr1'   , 'nanopore lib. A'), \n",
    "    ('MI_corr2'   , 'nanopore lib. A'),\n",
    "    ('MI_polishNoFilter'   , 'nanopore lib. A') \n",
    "\n",
    "    \n",
    "    ])\n",
    "\n",
    "readLabelsWorkflow = OrderedDict([\n",
    "        \n",
    "    ('S_raw','raw'), \n",
    "    ('S_ass','assembled'),\n",
    "\n",
    "\n",
    "    ('PB_raw','raw'), \n",
    "    ('PB_ass','assembled'),\n",
    "\n",
    "\n",
    "    ('MIpc_raw', 'raw 2D'), \n",
    "    ('MIpc_demux' , 'after debarcoding'),\n",
    "    ('MIpc_corr1' , 'after nanocorrect round 1'), \n",
    "    ('MIpc_corr2' , 'after nanocorrect round 2'),\n",
    "    #('MIpc_polishNoFilter'   , 'after nanopolish')  ,\n",
    "\n",
    "    ('MI_raw','raw 2D'), \n",
    "    ('MI_demux'   , 'after debarcoding'),\n",
    "    ('MI_corr1'   , 'after nanocorrect round 1'), \n",
    "    ('MI_corr2'   , 'after nanocorrect round 2') ,\n",
    "    ('MI_polishNoFilter'   , 'after nanopolish')  \n",
    "    \n",
    "    ])\n",
    "\n",
    "pacbioMuxData = 'inputFiles/pacbioMux.p'\n",
    "\n",
    "\n",
    "colorSet = {'MI': '#4daf4a',\n",
    "        'PB': '#377eb8',\n",
    "        'MIpc':'#ff7f00',\n",
    "        'S':'#984ea3'}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine nanocorrect and nanopolish data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:CARDBlast:Start merging the nanocorrect and nanopolish data\n",
      "ERROR:CARDBlast:The output/poreCamp.2D.min500.afterNPnofilter.fasta excist already, we are going to overwrite this!\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     582\n",
      "    1810\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:CARDBlast:Created a combined set of \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    2392\r\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "def createFullPolish(polishedPath ,correctedPath,nofilterPath):\n",
    "#def createFullPolish(polishedPath (=qPathOrdered['MI_polish']),correctedPath(=qPathOrdered['MI_corr2']),nofilterPath(=qPathOrdered['MI_polishNoFilter'])):\n",
    "    \"\"\"\n",
    "    We ran nanopore in such as fashion that only high coverage reads are polished (>30).\n",
    "    Reads that do not match this criteria (but can still have a high quality) are dropped from nanopolishing.\n",
    "    Consequently we lose sequence diversity, therefore we combine all the non-polished reads (thus just taking the reads after two rounds of nanocorrect) with the polished-reads. \n",
    "    But in order to this we first need to 'calculate' which reads were not polished (for the next version we can just let poreFUME take care of this)\n",
    "    \"\"\"\n",
    "    logger.info(\"Start merging the nanocorrect and nanopolish data\")\n",
    "    \n",
    "    if not os.path.isfile(polishedPath):\n",
    "        raise IOError('thisFilePath does not excist',polishedPath)\n",
    "\n",
    "    if not os.path.isfile(correctedPath):\n",
    "        raise IOError('thisFilePath does not excist',correctedPath)\n",
    "\n",
    "    if os.path.isfile(nofilterPath):\n",
    "        logger.error('The %s excist already, we are going to overwrite this!',nofilterPath)\n",
    "        \n",
    "    polishList = []\n",
    "    for record in SeqIO.parse(polishedPath, \"fasta\"): #Read in the records of the nanopolished data\n",
    "        polishList.append((record.id)) #Store all the name in a list\n",
    "\n",
    "    #Quick sanity check if list is unique\n",
    "    assert len(polishList) == len(set(polishList)), \"The polishList record names are not unique! poreFUME gives unique names, is something wrong?\"\n",
    "\n",
    "\n",
    "    polishSet = set(polishList)\n",
    "\n",
    "    nonPolished = []\n",
    "\n",
    "    for record in SeqIO.parse(correctedPath, \"fasta\"): #Read in the records of the nanopolished data\n",
    "        if record.id not in polishSet:\n",
    "            #print record.id\n",
    "            nonPolished.append(record) #Wrire record to a 'non-polished' list\n",
    "\n",
    "\n",
    "    output_handle = open(\"nonpolished.fasta.tmp\", \"w\")\n",
    "    SeqIO.write(nonPolished, output_handle, \"fasta\")\n",
    "    output_handle.close()\n",
    "\n",
    "    #!cat {polishedPath} | grep \">\" | wc -l\n",
    "    !cat nonpolished.fasta.tmp {polishedPath} > {nofilterPath} \n",
    "    #!cat nonpolished.fasta.tmp | grep \">\" | wc -l\n",
    "    !rm nonpolished.fasta.tmp\n",
    "    \n",
    "    logger.info(\"Created a combined set of \")\n",
    "    !cat {nofilterPath} | grep \">\" | wc -l\n",
    "\n",
    "#createFullPolish(qPathOrdered['MI_polish'],qPathOrdered['MI_corr2'],qPathOrdered['MI_polishNoFilter'])\n",
    "createFullPolish(qPathOrdered['MIpc_polish'],qPathOrdered['MIpc_corr2'],qPathOrdered['MIpc_polishNoFilter'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BLAST against the CARD database\n",
    "The [CARD database](https://card.mcmaster.ca/) is a _\"bioinformatic database of resistance genes, their products and associated phenotypes\"_. The reads from the various sequencing platforms are annotated using BLAST against the nucleotide CARD database. Please cite [McArthur et al. 2013. The Comprehensive Antibiotic Resistance Database. Antimicrobial Agents and Chemotherapy, 57, 3348-3357.](http://www.ncbi.nlm.nih.gov/pubmed/23650175) when using. \n",
    "\n",
    "The ```calcSegments()``` function will prune the BLAST report so only the highest scoring (by bitscore) segments are reported back. This is more robust than ORF calling (data not shown).\n",
    "_Note that this function takes quite a while on larger datasets (raw PacBio), even on multicore since calcSegements() is not optimized_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:CARDBlast:Initialization done\n",
      "INFO:CARDBlast:Starting CARD analysis loop with MI_polishNoFilter stored at output/LejlaControl.2D.min500bp.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:Start builing blast database\n",
      "INFO:CARDBlast:\n",
      "\n",
      "Building a new DB, current time: 11/19/2016 13:40:03\n",
      "New DB name:   /Users/evand/Downloads/nanorevision/production_setup/poreFUME_paper/inputFiles/n.fasta.protein.homolog.fasta\n",
      "New DB title:  inputFiles/n.fasta.protein.homolog.fasta\n",
      "Sequence type: Nucleotide\n",
      "Keep Linkouts: T\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 2181 sequences in 0.149022 seconds.\n",
      "\n",
      "INFO:CARDBlast:Start BLASTing for subsample output/LejlaControl.2D.min500bp.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "INFO:CARDBlast:Finished BLASTing for subsample output/LejlaControl.2D.min500bp.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:Parse BLAST data\n",
      "INFO:CARDBlast:Identified BLAST hits: 307663\n",
      "INFO:CARDBlast:Start identifing relevant segments\n",
      "INFO:CARDBlast:Finished identifing relevant segments\n",
      "INFO:CARDBlast:Identified single hits: 10807\n",
      "INFO:CARDBlast:\n",
      "ERROR:CARDBlast:mkdir: output/annotation: File exists\n",
      "\n",
      "INFO:CARDBlast:Saved annotationed in: output/annotation/LejlaControl.2D.min500bp.afterNPnofilter.fasta.annotated.csv\n",
      "INFO:CARDBlast:Starting CARD analysis loop with MIpc_polishNoFilter stored at output/poreCamp.2D.min500.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:Start builing blast database\n",
      "INFO:CARDBlast:\n",
      "\n",
      "Building a new DB, current time: 11/19/2016 13:45:52\n",
      "New DB name:   /Users/evand/Downloads/nanorevision/production_setup/poreFUME_paper/inputFiles/n.fasta.protein.homolog.fasta\n",
      "New DB title:  inputFiles/n.fasta.protein.homolog.fasta\n",
      "Sequence type: Nucleotide\n",
      "Keep Linkouts: T\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 2181 sequences in 0.155289 seconds.\n",
      "\n",
      "INFO:CARDBlast:Start BLASTing for subsample output/poreCamp.2D.min500.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "INFO:CARDBlast:Finished BLASTing for subsample output/poreCamp.2D.min500.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:Parse BLAST data\n",
      "INFO:CARDBlast:Identified BLAST hits: 88342\n",
      "INFO:CARDBlast:Start identifing relevant segments\n",
      "INFO:CARDBlast:Finished identifing relevant segments\n",
      "INFO:CARDBlast:Identified single hits: 2130\n",
      "INFO:CARDBlast:\n",
      "ERROR:CARDBlast:mkdir: output/annotation: File exists\n",
      "\n",
      "INFO:CARDBlast:Saved annotationed in: output/annotation/poreCamp.2D.min500.afterNPnofilter.fasta.annotated.csv\n"
     ]
    }
   ],
   "source": [
    "logger.info(\"Initialization done\")\n",
    "\n",
    "#########\n",
    "#Functions\n",
    "#########\n",
    "\n",
    "def blastDatabase(queryFile,dbFile):\n",
    "    \"\"\"\n",
    "    Returns a dataframe of a BLAST operation of the queryFile on the dbFile\n",
    "    \n",
    "    Parameters:\n",
    "    queryFile = path to file that will be used as query in the BLAST process\n",
    "    dbFile = path to file that will be used as database in the BLAST process\n",
    "    \"\"\"\n",
    "    #Create a blast DB http://www.ncbi.nlm.nih.gov/books/NBK279688/\n",
    "    \n",
    "    \n",
    "    assert type(queryFile) is StringType, \"queryFile is not a string: %r\" % queryFile \n",
    "    assert type(dbFile) is StringType, \"databaseFile is not a string: %r\" % dbFile \n",
    "        \n",
    "    #Test inputs\n",
    "    if not os.path.isfile(queryFile):\n",
    "        raise IOError('queryFile does not excist',queryFile)\n",
    "    \n",
    "    if not os.path.isfile(dbFile):\n",
    "        raise IOError('queryFile does not excist',dbFile)\n",
    "    \n",
    "    \n",
    "    #Create dbase\n",
    "    logger.info(\"Start builing blast database\")\n",
    "    process = Popen(['makeblastdb', '-in', str(dbFile), '-dbtype', 'nucl'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #wait till finished\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "\n",
    "    #Search with blast\n",
    "    logger.info(\"Start BLASTing for subsample %s\",queryFile)\n",
    "    process = Popen(['blastn','-db',str(dbFile),'-query',str(queryFile),'-max_hsps', '1', '-max_target_seqs', '1000', '-num_threads', str(CORES), '-outfmt' ,'10','-out','blastn.tmp.output'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #Wait till finished\n",
    "\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "    logger.info(\"Finished BLASTing for subsample %s\",queryFile)\n",
    "\n",
    "    if not os.path.isfile('blastn.tmp.output'):\n",
    "        raise RuntimeError('BLAST did not produce an alignment, one can implement an exception here. But for now stop')\n",
    "\n",
    "    #Read in blast results\n",
    "    dfBlast = None\n",
    "    dfBlast = pd.read_csv('blastn.tmp.output',names=['qseqid' ,'sseqid' ,'pident', 'length', 'mismatch', 'gapopen', 'qstart', 'qend' ,'sstart' ,'send','evalue', 'bitscore']) #load demuxed and corrected Minion dbase blasted by Sanger\n",
    "\n",
    "\n",
    "\n",
    "    os.remove(str(dbFile) +'.nhr') #cleanup subsample file blastDB\n",
    "    os.remove(str(dbFile) +'.nin') #cleanup subsample file blastDB\n",
    "    os.remove(str(dbFile) +'.nsq') #cleanup subsample file blastDB\n",
    "\n",
    "    os.remove('blastn.tmp.output') #cleanup blast output\n",
    "    \n",
    "    return dfBlast\n",
    "\n",
    "\n",
    "\n",
    "def calcGeneLength(row):\n",
    "    \"\"\"\n",
    "    The CARD database contains the position of the subject gene in the header name, this is originally stored in this script in \n",
    "    id2 but split into a subjectGeneSTart and subjectGeneEnd. Based on this we can calculate the original gene length\n",
    "\n",
    "    \"\"\"\n",
    "   \n",
    "    return abs(int(row['subjectGeneStart'])- int(row['subjectGeneEnd']))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def calcCoverage(row):\n",
    "    \"\"\"\n",
    "    We can calculate the coverage of the alignment by dividing the length of the alignemnt over the subjectGeneLength\n",
    "    There is a glitch in the length of the header and the real length of the DNA in the card database, so limit of at 100%\n",
    "    \"\"\"\n",
    "    #return (int(row['length'])/ float(row['subjectGeneLength']))*100\n",
    "\n",
    "    return (int(row['length'])/ float(row['subjectGeneLength']))*100 if (int(row['length'])/ float(row['subjectGeneLength']))*100 < 100 else 100\n",
    "\n",
    "\n",
    "def calcSegments(thisDF):\n",
    "    \"\"\"\n",
    "    Calculates the most relevant hit for each segment on the read\n",
    "    \n",
    "    Paramters:\n",
    "    thisDF = dataframe with the BLAST result for an INDIVIDUAL query. So don't pass the full BLAST table, but only from 1 qseqid, ie. thisDF[thisDF.qseqid == thisSeqid], where thisSeqid is the current seqid of interest\n",
    "    \"\"\"\n",
    "    #from PIL import Image, ImageDraw #For viz purposes\n",
    "    \n",
    "    #im = Image.new('RGBA', (6000, 1000), (0, 0, 0, 0))  #initialize a debug drawing screen\n",
    "    #draw = ImageDraw.Draw(im) \n",
    "    thisDF.reset_index(drop=True,inplace=True)\n",
    "\n",
    "    dy = 50 #off set for debug drawing\n",
    "    highscore = [] #Store the coordinates of the visited positions in a [start,end] format\n",
    "    highindex = [] #Store the index number (=blast hit) of each visited position\n",
    "    #for row in thisDF.itertuples(): #go through all the rows of the dataframe, not the most effient way! When this scales up write a new implemetation based on sorting\n",
    "\n",
    "   \n",
    "    for index,qstart,qend in zip(thisDF.index.values, thisDF.qstart, thisDF.qend): #this speeds up 5x compared to using iterrows(). However it is still slow. \n",
    "        #print row\n",
    "        if len(highscore) == 0: #First hit\n",
    "            highscore.append([qstart,qend]) #Add position to the position list\n",
    "            highindex.append(index)  #also keep track of the index\n",
    "            #draw.line((row['qstart'],dy,row['qend'],dy), fill=(255,255,255)) #draw bright\n",
    "        else:\n",
    "            doesOverlap = False\n",
    "            for thisScore in highscore: #Go through all the set positions\n",
    "                if calcOverlap(thisScore[0],thisScore[1],qstart,qend) > 0: #If there is an overlap \n",
    "                    #print 'Overlapping, skip'\n",
    "                    doesOverlap = True\n",
    "                    break #Exit loop, no need to continue\n",
    "            if doesOverlap == False: #We found a new segment\n",
    "                highscore.append([qstart,qend]) #Add position to list\n",
    "                #draw.line((row['qstart'],dy,row['qend'],dy), fill=(255,255,255))  #Draw bright\n",
    "                highindex.append(index)\n",
    "\n",
    "            else: #Old segment\n",
    "                pass\n",
    "                #draw.line((row['qstart'],dy,row['qend'],dy), fill=(55,55,55))     #Draw darker for debugging\n",
    "\n",
    "\n",
    "\n",
    "        #print row['qstart'],row['qend']\n",
    "\n",
    "\n",
    "        dy = dy + 2\n",
    "    #im.save('out.png',\"PNG\")\n",
    "    #im.show() #show a debug figure\n",
    " \n",
    "    return thisDF.iloc[highindex] #Return a dataframe with the most relevant hit on each segment\n",
    "\n",
    "def calcOverlap(a,b,c,d):\n",
    "    \"\"\"\n",
    "    Returns the overlap between two segments (ab) vs (cd)\n",
    "    \n",
    "    Parameters:\n",
    "    a = start position of segment AB\n",
    "    b = end position of segment AB\n",
    "    c = start position of segment CD\n",
    "    d = end position of segment CD\n",
    "    \n",
    "    Can be in any direction and any order\n",
    "    \"\"\"\n",
    "    return min([max([a,b]), max([c,d])]) - max([min([c,d]), min([a,b])]) #calculate overlap between two segments\n",
    "\n",
    "\n",
    "\n",
    "###########\n",
    "## Production loop\n",
    "###########\n",
    "\n",
    "#for thisQueryFile in qpath:\n",
    "for thisRunName, thisQueryFile in qPathOrdered.items():    \n",
    "    \n",
    "    if (thisRunName.split('_')[0] == 'S') | (thisRunName.split('_')[0] == 'PB'): #MI and #MIpc are already done by poreFUME (except the nofilter)\n",
    "        pass\n",
    "    else:\n",
    "        if (thisRunName.split('_')[1] == 'polishNoFilter'):  #the MI_polishNoFilter and MIpc_polishNoFilter we just generated, so they need to be annotated\n",
    "            pass\n",
    "        else:\n",
    "            continue\n",
    "\n",
    "    logger.info('Starting CARD analysis loop with %s stored at %s',thisRunName,thisQueryFile)\n",
    "    dfBlast =  blastDatabase(thisQueryFile,CARDdb) \n",
    "\n",
    "    #Parse the CARD database identifiers\n",
    "    logger.info(\"Parse BLAST data\")\n",
    "    dfBlast[['gb','id1','id2','ARO','GeneName']]= dfBlast.sseqid.str.split('|',expand=True) #Split the CARD header out\n",
    "    dfBlast[['subjectGeneStart','subjectGeneEnd']]= dfBlast.id2.str.split('-',expand=True) #Split the CARD gene length\n",
    "    dfBlast = dfBlast.drop('id2',1) #remove id2 redundand\n",
    "    dfBlast = dfBlast.drop('gb',1) #remove gb, does not contain information\n",
    "\n",
    "    dfBlast['subjectGeneLength'] = dfBlast.apply(calcGeneLength, axis=1)  \n",
    "\n",
    "    dfBlast['coverageOfSubjectGene'] = dfBlast.apply(calcCoverage, axis=1)\n",
    "\n",
    "\n",
    "    dfBlastSort = dfBlast.sort_values('bitscore',ascending=False) #sort by bitscore (in case it is not)\n",
    "\n",
    "    dfBlastSort = dfBlastSort.reset_index(drop=True) #Reindex based on bitscore\n",
    "    logger.info(\"Identified BLAST hits: \" + str(dfBlastSort.shape[0]))\n",
    "\n",
    "    logger.info(\"Start identifing relevant segments\")\n",
    "\n",
    "    logger.info(\"Finished identifing relevant segments\")\n",
    "\n",
    "    #In order to calculate each the most relevant BLAST hits we go through each query sequence and calculate the relevant hits\n",
    "    dfTotal = pd.DataFrame() #Initialze the final table\n",
    "\n",
    "    for thisSeqid in dfBlastSort.qseqid.unique(): #go through each query sequence individually\n",
    "        \n",
    "        dfTotal = dfTotal.append(calcSegments(dfBlastSort[dfBlastSort.qseqid == thisSeqid]), ignore_index=True) #Make a subset of the BLAST table containing only the current query seqience, next calculate the relevant blast hits and add them to the dfTotal table\n",
    "\n",
    "    logger.info(\"Identified single hits: \" + str(dfTotal.shape[0]))\n",
    "\n",
    "    ###Save results in querySuffix + '.annotated.csv'\n",
    "\n",
    "\n",
    "    dfTotal =  dfTotal.add_prefix('CARD:') \n",
    "\n",
    "    queryPostfix, querySuffix = os.path.split(thisQueryFile) #Parse filename\n",
    "\n",
    "    #Make sure there a directory to store the result\n",
    "    process = Popen(['mkdir', 'output/annotation'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #wait till finished\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "\n",
    "\n",
    "    dfTotal.to_csv('output/annotation/' + querySuffix + '.annotated.csv')\n",
    "\n",
    "    logger.info(\"Saved annotationed in: output/annotation/\" + querySuffix + \".annotated.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BLAST against Sanger dataset\n",
    "In order to investigate the sequence identity of the obtained reads a Sanger dataset is avialble whith contains a subset of the total colonies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def blastDatabaseSanger(queryFile,dbFile,fillUp = False,maxHits = 1):\n",
    "    \"\"\"\n",
    "    Returns a dataframe of a BLAST operation of the queryFile on the dbFile. Optionally the non-matched queries can be added with pident, length and bitscore = 0\n",
    "    \n",
    "    Parameters:\n",
    "    queryFile = path to file that will be used as query in the BLAST process\n",
    "    dbFile = path to file that will be used as database in the BLAST process\n",
    "    fillUp = False, turn on to add non-matched query sequences\n",
    "    maxHits = 1, the max number of hits per query.\n",
    "    \"\"\"\n",
    "    #Create a blast DB http://www.ncbi.nlm.nih.gov/books/NBK279688/\n",
    "    \n",
    "    \n",
    "    assert type(queryFile) is StringType, \"queryFile is not a string: %r\" % queryFile \n",
    "    assert type(dbFile) is StringType, \"databaseFile is not a string: %r\" % dbFile \n",
    "        \n",
    "    #Test inputs\n",
    "    if not os.path.isfile(queryFile):\n",
    "        raise IOError('queryFile does not exist',queryFile)\n",
    "    \n",
    "    if not os.path.isfile(dbFile):\n",
    "        raise IOError('queryFile does not exist',dbFile)\n",
    "    \n",
    "    \n",
    "    #Create dbase\n",
    "    logger.info(\"Start builing blast database\")\n",
    "    process = Popen(['makeblastdb', '-in', str(dbFile), '-dbtype', 'nucl'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #wait till finished\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "\n",
    "    #Search with blast\n",
    "    logger.info(\"Start BLASTing for subsample %s\",queryFile)\n",
    "    process = Popen(['blastn','-db',str(dbFile),'-query',str(queryFile),'-max_hsps', '1', '-max_target_seqs', str(maxHits), '-num_threads', str(CORES), '-outfmt' ,'10','-out','blastn.tmp.output'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #Wait till finished\n",
    "\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "    logger.info(\"Finished BLASTing for subsample %s\",queryFile)\n",
    "\n",
    "    if not os.path.isfile('blastn.tmp.output'):\n",
    "        raise RuntimeError('BLAST did not produce an alignment, one can implement an exception here. But for now stop')\n",
    "\n",
    "    #Read in blast results\n",
    "    dfBlast = None\n",
    "    dfBlast = pd.read_csv('blastn.tmp.output',names=['qseqid' ,'sseqid' ,'pident', 'length', 'mismatch', 'gapopen', 'qstart', 'qend' ,'sstart' ,'send','evalue', 'bitscore']) #load demuxed and corrected Minion dbase blasted by Sanger\n",
    "\n",
    "    if fillUp:\n",
    "        queryList = []\n",
    "        queryLength = []\n",
    "        for record in SeqIO.parse(queryFile, \"fasta\"):\n",
    "            queryList.append(record.id)\n",
    "            queryLength.append( len(record.seq))\n",
    "\n",
    "        blastList = dfBlast.qseqid.tolist()\n",
    "\n",
    "\n",
    "        nohitList = set(queryList) - set(blastList)\n",
    "        for thisHit in nohitList:\n",
    "\n",
    "            thisAddition = {\n",
    "            'qseqid': thisHit,\n",
    "            'length' : 0,\n",
    "            'pident': 0,\n",
    "            'bitscore': 0\n",
    "            }\n",
    "\n",
    "            dfBlast = dfBlast.append(thisAddition,ignore_index=True)\n",
    "\n",
    "        #Sanity check\n",
    "        if len(queryList) == dfBlast.shape[0]:\n",
    "            logger.info('OK queryList and blast result are same length. Continue')\n",
    "        else:\n",
    "            raise ValueError('the resulting dataframe and the initial query list are not the same length!')\n",
    "        \n",
    "\n",
    "\n",
    "    os.remove(str(dbFile) +'.nhr') #cleanup subsample file blastDB\n",
    "    os.remove(str(dbFile) +'.nin') #cleanup subsample file blastDB\n",
    "    os.remove(str(dbFile) +'.nsq') #cleanup subsample file blastDB\n",
    "\n",
    "    os.remove('blastn.tmp.output') #cleanup blast output\n",
    "    \n",
    "    return dfBlast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:CARDBlast:Start builing blast database\n",
      "INFO:CARDBlast:\n",
      "\n",
      "Building a new DB, current time: 11/19/2016 13:47:48\n",
      "New DB name:   /Users/evand/Downloads/nanorevision/production_setup/poreFUME_paper/inputFiles/sanger.total.aftertrim.removeCT.min500bp.fasta\n",
      "New DB title:  inputFiles/sanger.total.aftertrim.removeCT.min500bp.fasta\n",
      "Sequence type: Nucleotide\n",
      "Keep Linkouts: T\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 1568 sequences in 0.102869 seconds.\n",
      "\n",
      "INFO:CARDBlast:Start BLASTing for subsample output/LejlaControl.2D.min500bp.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "INFO:CARDBlast:Finished BLASTing for subsample output/LejlaControl.2D.min500bp.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "ERROR:CARDBlast:mkdir: output/blastAgainstSanger: File exists\n",
      "\n",
      "INFO:CARDBlast:Saved annotationed in: output/blastAgainstSanger/LejlaControl.2D.min500bp.afterNPnofilter.fasta.againstSanger.csv\n",
      "INFO:CARDBlast:Start builing blast database\n",
      "INFO:CARDBlast:\n",
      "\n",
      "Building a new DB, current time: 11/19/2016 13:48:26\n",
      "New DB name:   /Users/evand/Downloads/nanorevision/production_setup/poreFUME_paper/inputFiles/sanger.total.aftertrim.removeCT.min500bp.fasta\n",
      "New DB title:  inputFiles/sanger.total.aftertrim.removeCT.min500bp.fasta\n",
      "Sequence type: Nucleotide\n",
      "Keep Linkouts: T\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 1568 sequences in 0.1087 seconds.\n",
      "\n",
      "INFO:CARDBlast:Start BLASTing for subsample output/poreCamp.2D.min500.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "INFO:CARDBlast:Finished BLASTing for subsample output/poreCamp.2D.min500.afterNPnofilter.fasta\n",
      "INFO:CARDBlast:\n",
      "ERROR:CARDBlast:mkdir: output/blastAgainstSanger: File exists\n",
      "\n",
      "INFO:CARDBlast:Saved annotationed in: output/blastAgainstSanger/poreCamp.2D.min500.afterNPnofilter.fasta.againstSanger.csv\n"
     ]
    }
   ],
   "source": [
    "for thisRunName, thisQueryFile in qPathOrdered.items():\n",
    "    if thisRunName == 'S_raw':\n",
    "        continue\n",
    "        \n",
    "    thisDF = blastDatabaseSanger(thisQueryFile,qPathOrdered['S_raw'],fillUp=False)\n",
    "    process = Popen(['mkdir', 'output/blastAgainstSanger'], stdout=PIPE, stderr=PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    process.wait() #wait till finished\n",
    "    logger.info(str(stdout))\n",
    "    if stderr:\n",
    "        logger.error(str(stderr))\n",
    "\n",
    "    queryPostfix, querySuffix = os.path.split(thisQueryFile) #Parse filename\n",
    "\n",
    "    thisDF.to_csv('output/blastAgainstSanger/' + querySuffix + '.againstSanger.csv')\n",
    "\n",
    "    \n",
    "    logger.info(\"Saved annotationed in: output/blastAgainstSanger/\" + querySuffix + \".againstSanger.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}