docs/Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Workflow and web application for annotation of NCBI BioProject transcriptome data \n",
    "\n",
    "Roberto Vera Alvarez$^{1}$, Newton Medeiros Vidal$^{1}$, Gina A Garzón-Martínez$^{2}$, Luz S Barrero$^{2}$, Richa Agarwala$^{3}$, David Landsman$^{1}$ and Leonardo Mariño-Ramírez$^{1,*}$\n",
    "\n",
    "$^{1}$Computational Biology Branch, National Center for Biotechnology Information, National Library of Medicine, National Institutes of Health, Bethesda, Maryland, USA\n",
    "\n",
    "$^{2}$Colombian Corporation for Agricultural Research (CORPOICA), Bogota, Colombia\n",
    "\n",
    "$^{3}$Information Engineering Branch, National Center for Biotechnology Information, National Library of Medicine, National Institutes of Health, Bethesda, Maryland, USA\n",
    "\n",
    "$^{*}$Corresponding author: Email: marino@ncbi.nlm.nih.gov Address: Building 38A, Room 6S614-M, 8600 Rockville Pike MSC 6075. Bethesda, MD 20894-6075 Tel: 1-301-402-3708\n",
    "\n",
    "## Abstract\n",
    "\n",
    "Transcriptome data is growing exponentially due to the latest improvements of experimental devices. In response to that, large central resources, such as those of the National Center for Biotechnology Information (NCBI), are continually adapting their computational infrastructure to accommodate this data influx. New specialized databases have been created like Transcriptome Shotgun Assembly Sequence Database (TSA) and Sequence Read Archive (SRA) aiding the creation of centralized repositories. Although these databases are in a constant development, they do not include automatic pipelines to increase the annotation of the deposited data. Therefore, third party applications are required to achieve that aim. Here we present an automatic workflow and web application for the annotation of transcriptome data. The workflow creates secondary data like sequencing read and BLAST alignments which are available through the web application. It is based on freely available bioinformatics tools and a set of in-house developed scripts. The interactive web application provides a search engine and different browse utilities. Additionally, graphical views of the transcript alignments are available through SeqViewer, a friendly embedded tool developed by NCBI for viewing biological sequence data. The web application is tightly integrated to the other NCBI web applications and tools extending the possibility of data processing and interconnectivity. A case study is presented for the species Physalis peruviana with data generated from the BioProject with ID 67621.\n",
    "\n",
    "PhysalisDB URL: http://www.ncbi.nlm.nih.gov/projects/physalis/\n",
    "\n",
    "Workflow Github repository: https://github.com/r78v10a07/trans-annot-notebook\n",
    "\n",
    "Jupyter Notebook: https://github.com/r78v10a07/trans-annot-notebook/blob/master/docs/Notebook.ipynb \n",
    "\n",
    "\n",
    "## Notebook\n",
    "\n",
    "This is a Jupyter Notebook which illustrate the workflow used for the transcript annotation process. This Github project can be clone and executed in a compatible environment. Please, note, this is just a demo version which will annotate only 10 transcripts and can be executed in a personal laptop.  \n",
    "\n",
    "## Requirements\n",
    "\n",
    "### Basic programs\n",
    "\n",
    "1. Python3\n",
    "2. virtualenv3\n",
    "3. C++\n",
    "\n",
    "### Bioinformatics programs\n",
    "\n",
    "These programs have to be installed and included in the user PATH variable.\n",
    "\n",
    "1. EDirect (https://www.ncbi.nlm.nih.gov/books/NBK179288/)\n",
    "2. SRA Toolkit (https://www.ncbi.nlm.nih.gov/books/NBK158900/)\n",
    "3. Bowtie 2 (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml)\n",
    "4. SAMTools (http://www.htslib.org/)\n",
    "5. BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi)\n",
    "\n",
    "## Instalation of the Notebook\n",
    "\n",
    "Execute this commands in a BASH terminal just the first time to get the notebook dependencies installed.\n",
    "\n",
    "```\n",
    "$ git clone https://github.com/r78v10a07/trans-annot-notebook\n",
    "$ cd trans-annot-notebook\n",
    "\n",
    "# Compiling the myBlast2GO C++ tool\n",
    "$ cd src/myblast2go/\n",
    "$ make\n",
    "$ cd ../../\n",
    "\n",
    "# Creating a virtual environment\n",
    "$ virtualenv3 venv\n",
    "$ source venv/bin/activate\n",
    "$ pip install -r requirements/jupyter.txt\n",
    "$ chmod a+x ./bin/*\n",
    "```\n",
    "\n",
    "Execute this script to start the jupyter server.\n",
    "\n",
    "```\n",
    "$ ./bin/runNotebook.sh\n",
    "```\n",
    "\n",
    "After this open a web browser with this URL: http://localhost:9096\n",
    "\n",
    "Click on docs and then, click on Notebook.ipynb to open the notebook\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Global imports\n",
    "import os\n",
    "import re\n",
    "import csv\n",
    "import sqlite3\n",
    "import xmltodict\n",
    "from Bio import Entrez\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Email to be used by Entrez\n",
    "Entrez.email = \"veraalva@ncbi.nlm.nih.gov\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Environment variables\n",
    "os.environ['WORKDIR'] = os.path.abspath('../')\n",
    "os.environ['CONFIG'] = os.environ['WORKDIR'] + '/config'\n",
    "os.environ['DATA'] = os.environ['WORKDIR'] + '/data'\n",
    "os.environ['BIN'] = os.environ['WORKDIR'] + '/bin'\n",
    "os.environ['RESULTS'] = os.environ['WORKDIR'] + '/results'\n",
    "os.environ['DOC'] = os.environ['WORKDIR'] + '/docs'\n",
    "os.environ['SRC'] = os.environ['WORKDIR'] + '/src'\n",
    "os.environ['SQLITEDBSCHEMA'] = os.environ['CONFIG'] + '/sqlite_database_schema.sql'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# BioProject, SRA Id and SQLite database\n",
    "os.environ['BIOPROJECT'] = 'PRJNA67621'\n",
    "os.environ['SRA'] = 'SRP005904'\n",
    "os.environ['SQLITEDB'] = os.environ['RESULTS'] + '/' + os.environ['BIOPROJECT'] + '.sqlite'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "memory\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "# Initialize SQLite database\n",
    "if [ ! -e \"${SQLITEDB}\" ]\n",
    "then\n",
    "    cat $SQLITEDBSCHEMA | sqlite3 $SQLITEDB\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The data is ready\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "# Downloading the SRR data using fastq-dump\n",
    "cd $DATA\n",
    "if [ ! -e \"${SRA}\" ]\n",
    "then\n",
    "    mkdir ${SRA}\n",
    "    sh $BIN/download_SRA.sh -b ${BIOPROJECT} -w ./${SRA}\n",
    "else\n",
    "    echo \"The data is ready\"\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Download the UniProt ID mapping to link Gi with Protein ID\n",
    "cd $DATA\n",
    "if [ ! -e \"idmapping_selected.tab\" ]\n",
    "then\n",
    "    wget -o wget_idmapping.log ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/idmapping_selected.tab.gz\n",
    "    gunzip idmapping_selected.tab.gz\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Downloading ec2go map\n",
    "cd $DATA\n",
    "if [ ! -e \"ec2go\" ]\n",
    "then\n",
    "    wget -o wget_ec2go.log http://geneontology.org/external2go/ec2go\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Get the SRA IDs related to the BioProject from NCBI\n",
    "handle = Entrez.esearch(db=\"sra\", term=os.environ['BIOPROJECT'], retmax=50)\n",
    "sra_id = Entrez.read(handle)\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Retrieving the first 10 transcripts from the BioProject \n",
    "# Modify the parameter retmax in order to increase the number of transcripts to be processed\n",
    "handle = Entrez.esearch(db=\"nucleotide\", retmax=11, term=os.environ['BIOPROJECT'] + \"[BioProject]\")\n",
    "record = Entrez.read(handle)\n",
    "handle.close()\n",
    "with open(os.environ['DATA'] + '/bioproject-' + os.environ['BIOPROJECT'] + '.fsa', 'w') as fout:\n",
    "    for id in record['IdList']:\n",
    "        handle = Entrez.efetch(db=\"nucleotide\", id=id, rettype=\"fasta\", retmode=\"text\")\n",
    "        for rec in SeqIO.parse(handle, \"fasta\"):\n",
    "            if rec.seq:\n",
    "                SeqIO.write(rec, fout, \"fasta\")\n",
    "        handle.close"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Retrieving the bioproject from NCBI\n",
    "handle = Entrez.esearch(db=\"bioproject\", term=os.environ['BIOPROJECT'])\n",
    "bioproject_ids = Entrez.read(handle)\n",
    "handle.close()\n",
    "os.environ['BIOPROJECTID'] = bioproject_ids['IdList'][0]\n",
    "handle = Entrez.efetch(db=\"bioproject\", id=os.environ['BIOPROJECTID'])\n",
    "bioproject = xmltodict.parse(handle.read())\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Inserting the BioProject data into the SQLite DB\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "c.execute(\"INSERT INTO 'bioprojects' ('id', 'path') VALUES (\" + os.environ['BIOPROJECTID'] + \", '\" + os.environ['WORKDIR'] + \"')\")\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Retirve the first experiment related to the SRA IDs from NCBI\n",
    "sra = {}\n",
    "for id in sra_id['IdList']:\n",
    "    handle = Entrez.efetch(db=\"sra\", id=id, rettype=\"full\", retmode=\"full\")\n",
    "    sra[id] = xmltodict.parse(handle.read())\n",
    "    handle.close()\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Insert the experiments and runs in the database\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "for key in sra:\n",
    "    c.execute(\"INSERT INTO 'experiments' ('accession', 'title', 'bioproject_id') VALUES ('\" + key + \"', '\" + sra[key]['EXPERIMENT_PACKAGE_SET']['EXPERIMENT_PACKAGE']['EXPERIMENT']['TITLE'] + \"', '\" + os.environ['BIOPROJECTID'] + \"')\")\n",
    "    conn.commit()\n",
    "    c.execute(\"SELECT id FROM 'experiments' where accession = '\" + key + \"'\")\n",
    "    experiment_id = c.fetchone()\n",
    "    if type(sra[key]['EXPERIMENT_PACKAGE_SET']['EXPERIMENT_PACKAGE']['RUN_SET']['RUN']) is list: \n",
    "        for run in sra[key]['EXPERIMENT_PACKAGE_SET']['EXPERIMENT_PACKAGE']['RUN_SET']['RUN']:\n",
    "            c.execute(\"INSERT INTO 'runs' ('accession', 'experiment_id') VALUES ('\" + run['@accession'] + \"', \" + str(experiment_id[0]) + \")\")\n",
    "            conn.commit()  \n",
    "    else:\n",
    "        c.execute(\"INSERT INTO 'runs' ('accession', 'experiment_id') VALUES ('\" + sra[key]['EXPERIMENT_PACKAGE_SET']['EXPERIMENT_PACKAGE']['RUN_SET']['RUN']['@accession'] + \"', \" + str(experiment_id[0]) + \")\")\n",
    "        conn.commit()        \n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Parsinf the FASTA file and insert the IDs into the database\n",
    "# A folder for each transcript is created into the results folder\n",
    "# Bowtie2 and blast are executed\n",
    "\n",
    "if not os.path.isdir(os.environ['RESULTS'] + '/transcripts'):\n",
    "    os.makedirs(os.environ['RESULTS'] + '/transcripts')\n",
    "    \n",
    "os.chdir(os.environ['RESULTS'] + '/transcripts')\n",
    "\n",
    "rpoint = re.compile('[\\|\\.]')\n",
    "rPipe = re.compile('\\|')\n",
    "\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open(os.environ['DATA'] + '/bioproject-' + os.environ['BIOPROJECT'] + '.fsa', \"r\") as handle:\n",
    "    for record in SeqIO.parse(handle, \"fasta\"):\n",
    "        acc = record.id.split('.')\n",
    "        title = record.description.strip()\n",
    "        \n",
    "\n",
    "        if not os.path.exists(acc[0]):\n",
    "            os.makedirs(acc[0])\n",
    "\n",
    "        os.chdir(acc[0])\n",
    "        \n",
    "        toInsert = (acc[0], acc[1], title, len(record.seq), os.path.abspath('./'), os.environ['BIOPROJECTID'])\n",
    "        c.execute(\"INSERT INTO nucleotides ('accession', 'version', 'title', 'len', 'path','bioproject_id')  VALUES (?, ?, ?, ?, ?, ?)\", toInsert)\n",
    "        \n",
    "        if not os.path.exists(acc[0] + \".fasta\"):  \n",
    "            output = open(acc[0] + \".fasta\", \"w\")\n",
    "            SeqIO.FastaIO.FastaWriter(output, wrap=80).write_file([record])\n",
    "            output.close()\n",
    "        \n",
    "        if not os.path.exists(acc[0] + \"_sorted.bam\"):\n",
    "            os.environ['ID'] = acc[0]\n",
    "            !sh $BIN/processTranscript.sh -i $ID -s $DATA/$SRA\n",
    "        \n",
    "        os.chdir(\"../\")\n",
    "\n",
    "\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Retrieving taxonomies, iteration \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "# Proecessing Blast results to be inserted into the database\n",
    "cd $RESULTS\n",
    "sh $BIN/processBlastResults.sh -t ./transcripts -b $BIN/buildLine.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loading the taxonomies to the SQLite DB\n",
    "os.chdir(os.environ['RESULTS'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('taxonomies_to_load.tsv') as handle:\n",
    "    rows = csv.reader(handle, delimiter='\\t')\n",
    "    c.executemany(\"INSERT INTO taxonomies (id, name, division) VALUES (?, ?, ?);\", rows)\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loading the proteins to the SQLite DB\n",
    "os.chdir(os.environ['RESULTS'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('proteins_to_load.tsv') as handle:\n",
    "    rows = csv.reader(handle, delimiter='\\t')\n",
    "    c.executemany(\"INSERT INTO proteins (taxonomy_id, gi, accession, title) VALUES (?, ?, ?, ?);\", rows)\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loading the Blast results to the SQLite DB\n",
    "os.chdir(os.environ['RESULTS'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('blastx_to_load.tsv') as handle:\n",
    "    rows = csv.reader(handle, delimiter='\\t')\n",
    "    c.executemany(\"INSERT INTO blastx_temp (qacc, sgi, slen, qstart, qend, sstart, send, evalue, bitscore, score, length, pident, nident, mismatch, positive, gapopen, gaps, ppos, qcovs, qcovhsp) VALUES (?, ?, ?, ?, ?, ?, ?, ?,?, ?, ?, ?,?, ?, ?, ?,?, ?, ?, ?);\", rows)\n",
    "c.execute(\"insert into blastx (slen,qstart,qend,sstart,send,evalue,bitscore,score,length,pident,nident,mismatch,positive,gapopen,gaps,ppos,qcovs,qcovhsp,protein_id, nucleotide_id) select bt.slen,bt.qstart, bt.qend, bt.sstart, bt.send, bt.evalue, bt.bitscore, bt.score, bt.length, bt.pident, bt.nident, bt.mismatch, bt.positive, bt.gapopen, bt.gaps, bt.ppos, bt.qcovs, bt.qcovhsp, p.id, n.id from blastx_temp bt inner join nucleotides n on n.accession || '.' || n.version = bt.qacc inner join proteins p on bt.sgi = p.gi\")    \n",
    "c.execute(\"delete from blastx_temp\")\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Running our Blast to GO annotation program\n",
    "cd $RESULTS\n",
    "if [ ! -e \"blastGI2GO.tsv\" ]\n",
    "then\n",
    "    !$BIN/myBlast2GO -d $DATA/idmapping_selected.tab -g 2 -i ./transcripts/ -l 9 -c 50 -o blastGI2GO.tsv -e blastGI2GO.err -r -s _blastx.tsv\n",
    "fi    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loading the GI2GO results to the SQLite DB\n",
    "os.chdir(os.environ['RESULTS'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('blastGI2GO.tsv') as handle:\n",
    "    rows = csv.reader(handle, delimiter='\\t')\n",
    "    c.executemany(\"INSERT INTO gi2go (gi, go) VALUES (?, ?);\", rows)\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Loading the GO data from a previusly created tsv file from JBioWH to the SQLite DB\n",
    "os.chdir(os.environ['DATA'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('go_JBioWH.tsv') as handle:\n",
    "    rows = csv.reader(handle, delimiter='\\t')\n",
    "    c.executemany(\"INSERT INTO go (accession, name, description, name_space) VALUES (?, ?, ?, ?);\", rows)\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Creating links between nucleotides and GO terms\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "c.execute(\"insert OR IGNORE into nucleotides_go (nucleotides_id, go_id) select n.id,go.id from nucleotides n inner join blastx b on n.id = b.nucleotide_id inner join proteins p on p.id = b.protein_id inner join gi2go g on g.gi = p.gi inner join go on go.accession = g.go;\")\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Creates the GO statistics for the 20 more relevant terms by GO name space\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "c.execute(\"insert into go_statistics (go_id,transcripts) select g.id,count(n.nucleotides_id) as c from go g inner join nucleotides_go n on n.go_id = g.id where g.name_space = 'biological_process' group by g.id order by c desc limit 20 ;\")\n",
    "c.execute(\"insert into go_statistics (go_id,transcripts) select g.id,count(n.nucleotides_id) as c from go g inner join nucleotides_go n on n.go_id = g.id where g.name_space = 'molecular_function' group by g.id order by c desc limit 20 ;\")\n",
    "c.execute(\"insert into go_statistics (go_id,transcripts) select g.id,count(n.nucleotides_id) as c from go g inner join nucleotides_go n on n.go_id = g.id where g.name_space = 'cellular_component' group by g.id order by c desc limit 20 ;\")\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Inserting EC2GO data and create cross-references with the transcripts\n",
    "os.chdir(os.environ['DATA'])\n",
    "conn = sqlite3.connect(os.environ['SQLITEDB'])\n",
    "c = conn.cursor()\n",
    "with open('ec2go') as handle:\n",
    "    rows = csv.reader(handle, delimiter=';')\n",
    "    for r in rows:\n",
    "        if len(r) == 2:\n",
    "            ec_rows = r[0].split('>')\n",
    "            c.execute(\"insert into ec2go (ec, go) VALUES ('\" + ec_rows[0].strip() + \"', '\" + r[1].strip() + \"')\")\n",
    "conn.commit()\n",
    "c.execute(\"insert into enzymes ('accession') select ec from ec2go group by ec\")\n",
    "c.execute(\"insert into nucleotides_enzymes ('nucleotides_id', 'enzymes_id') select n.nucleotides_id,e.id  from nucleotides_go n inner join go g on g.id = n.go_id inner join ec2go eg on eg.go = g.accession inner join enzymes e on e.accession = eg.ec group by n.nucleotides_id,e.id\")\n",
    "conn.commit()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}