iPythonNotebooks/From_functional_roles_to_gap-filling.ipynb

{
 "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"
  },
  "name": "",
  "signature": "sha256:bfffc1cc29ba4bf4fe871b04c45c069e7bfc687d1ad1eb87751bf215c79e8a5b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# How to gap-fill a genome scale metabolic model\n",
      "\n",
      "## Getting started\n",
      "\n",
      "### Installing libraries\n",
      "\n",
      "Before you start, you will need to install a couple of libraries:\n",
      "   \n",
      "The [ModelSeedDatabase](https://github.com/ModelSEED/ModelSEEDDatabase) has all the biochemistry we'll need. You can install that with `git clone`.\n",
      "   \n",
      "The [PyFBA](http://linsalrob.github.io/PyFBA) library has detailed [installation instructions](http://linsalrob.github.io/PyFBA/installation.html). Don't be scared, its mostly just `pip install`.\n",
      "\n",
      "(Optional) Also, get the [SEED Servers](https://github.com/linsalrob/SEED_Servers_Python) as you can get a lot of information from them. You can install the git python repo from github.  Make sure that the SEED_Servers_Python is in your PYTHONPATH.\n",
      "\n",
      "We start with importing some modules that we are going to use. \n",
      "\n",
      "We import *sys* so that we can use standard out and standard error if we have some error messages.<br>\n",
      "We import *copy* so that we can make a deep copy of data structures for later comparisons.<br>\n",
      "We import *print_function* so that we can use the print using Python 3.X syntax.<br>\n",
      "Then we import the *PyFBA* module to get started."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "import copy\n",
      "import PyFBA\n",
      "from __future__ import print_function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Running a basic model\n",
      "\n",
      "The SBML model is great if you have built a model elsewhere, but what if you want to build a model from a genome?\n",
      "\n",
      "We typically start with an *assigned_functions* file from RAST. The easiest way to find that is in the RAST directory by choosing `Genome Directory` from the `Downloads` menu on the job details page.\n",
      "\n",
      "For this example, [here is an *assigned_functions* file](https://raw.githubusercontent.com/linsalrob/PyFBA/master/example_data/Citrobacter/ungapfilled_model/citrobacter.assigned_functions) from our *Citrobacter* model that you can download to the same directory as this iPython notebook. Notice that it has two columns: the first column is the protein ID (using SEED standard IDs that start with `fig|`, and then have the taxonomy ID and version number of the genome, and then `peg` to indicate *protein encoding gene*, `rna` to indicate *RNA*, `crispr_spacer` to indicate *crispr spacers* or other acronym, followed by the feature number. After the tab is the *functional role* of that feature. Download that file to use in this test. \n",
      "\n",
      "We start by converting this *assigned_functions* file to a list of reactions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "assigned_functions = PyFBA.parse.read_assigned_functions(\"citrobacter.assigned_functions\")\n",
      "roles = set([i[0] for i in [list(j) for j in assigned_functions.values()]])\n",
      "print(\"There are {} unique roles in this genome\".format(len(roles)))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": false
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "There are 3594 unique roles in this genome\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Convert those roles to reactions. We start with a dict of roles and reactions, but we only need a list of unique reactions, so we convert the keys to a set."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "roles_to_reactions = PyFBA.filters.roles_to_reactions(roles)\n",
      "reactions_to_run = set()\n",
      "for role in roles_to_reactions:\n",
      "    reactions_to_run.update(roles_to_reactions[role])\n",
      "print(\"There are {}\".format(len(reactions_to_run)),\n",
      "      \"unique reactions associated with this genome\".format(len(reactions_to_run)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "There are 1392 unique reactions associated with this genome\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Read all the reactions and compounds in our database\n",
      "\n",
      "We read all the reactions, compounds, and enzymes in the [ModelSEEDDatabase](https://github.com/ModelSEED/ModelSEEDDatabase) into three data structures. Each one is a dictionary with a string representation of the object as the key and the PyFBA object as the value.\n",
      "\n",
      "We modify the reactions specifically for Gram negative models (there are also options for Gram positive models, Mycobacterial models, general microbial models, and plant models)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "compounds, reactions, enzymes = \\\n",
      "    PyFBA.parse.model_seed.compounds_reactions_enzymes('gramnegative')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Update reactions to run, making sure that all reactions are in the list!\n",
      "\n",
      "There are some reactions that come from functional roles that do not appear in the reactions list. We're working on tracking these down, but for now we just check that all reaction IDs in *reactions_to_run* are in *reactions*, too."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tempset = set()\n",
      "for r in reactions_to_run:\n",
      "    if r in reactions:\n",
      "        tempset.add(r)\n",
      "    else:\n",
      "        print(\"Reaction ID {}\".format(r),\n",
      "              \"is not in our reactions list. Skipped\",\n",
      "              file=sys.stderr)\n",
      "reactions_to_run = tempset"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "Reaction ID rxn37218 is not in our reactions list. Skipped\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Test whether these reactions grow on ArgonneLB media\n",
      "\n",
      "We can test whether this set of reactions grows on ArgonneLB media. The media is the same one we used above, and you can download the [ArgonneLB.txt](https://raw.githubusercontent.com/linsalrob/PyFBA/master/media/ArgonneLB.txt) and text file and put it in the same directory as this iPython notebook to run it.\n",
      "\n",
      "(Note: we don't need to convert the media components, because the media and compounds come from the same source.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "media = PyFBA.parse.read_media_file('ArgonneLB.txt')\n",
      "print(\"Our media has {} components\".format(len(media)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Our media has 65 components\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Define a biomass equation\n",
      "\n",
      "The biomass equation is the part that says whether the model will grow! This is a [metabolism.reaction.Reaction](https://github.com/linsalrob/PyFBA/blob/master/PyFBA/metabolism/reaction.py) object."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "biomass_equation = PyFBA.metabolism.biomass_equation('gramnegative')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Run the FBA\n",
      "\n",
      "With the reactions, compounds, reactions_to_run, media, and biomass model, we can test whether the model grows on this media."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"Initial run has a biomass flux value of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": false
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Initial run has a biomass flux value of 4.06345328793e-14 --> Growth: False\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Gap-fill the model\n",
      "\n",
      "Since the model does not grow on ArgonneLB we need to gap-fill it to ensure growth. There are several ways that we can gap-fill, and we will work through them until we get growth.\n",
      "\n",
      "As you will see, we update the *reactions_to_run list* each time, and keep the media and everything else consistent. Then we just need to run the FBA like we have done above and see if we get growth.\n",
      "\n",
      "We also keep a copy of the original *reactions_to_run*, and a list with all the reactions that we are adding, so once we are done we can go back and bisect the reactions that are added."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "added_reactions = []\n",
      "original_reactions_to_run = copy.copy(reactions_to_run)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Media import reactions\n",
      "\n",
      "We need to make sure that the cell can import everything that is in the media... otherwise it won't be able to grow. Be sure to only do this step if you are certain that the cell can grow on the media you are testing."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "media_reactions = PyFBA.gapfill.suggest_from_media(compounds, reactions,\n",
      "                                                   reactions_to_run, media)\n",
      "added_reactions.append((\"media\", media_reactions))\n",
      "reactions_to_run.update(media_reactions)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The biomass reaction has a flux of 1.4583713765e-15 --> Growth: False\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Reactions from closely related organisms\n",
      "\n",
      "We also gap-fill on closely related organisms. We assume that an organism is most likely to have reactions in its genome that are similar to those in closely related organisms. \n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reactions_from_other_orgs = PyFBA.gapfill.suggest_from_roles(\"closest.genomes.roles\",\n",
      "                                                             reactions)\n",
      "added_reactions.append((\"close genomes\", reactions_from_other_orgs))\n",
      "reactions_to_run.update(reactions_from_other_orgs)"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The biomass reaction has a flux of -3.23181490595e-14 --> Growth: False\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Essential reactions\n",
      "\n",
      "There are ~100 reactions that are in every model we have tested, and we construe these to be essential for all models, so we typically add these next!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "essential_reactions = PyFBA.gapfill.suggest_essential_reactions()\n",
      "added_reactions.append((\"essential\", essential_reactions))\n",
      "reactions_to_run.update(essential_reactions)"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": false
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The biomass reaction has a flux of 3.38622145944e-14 --> Growth: False\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Subsystems\n",
      "\n",
      "The reactions connect us to subsystems (see [Overbeek et al. 2014](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3965101/)), and this test ensures that all the subsystems are complete. We add reactions required to complete the subsystem."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subsystem_reactions = \\\n",
      "    PyFBA.gapfill.suggest_reactions_from_subsystems(reactions,\n",
      "                                                    reactions_to_run,\n",
      "                                                    threshold=0.5)\n",
      "added_reactions.append((\"subsystems\", subsystem_reactions))\n",
      "reactions_to_run.update(subsystem_reactions)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The biomass reaction has a flux of 4.13676770622e-13 --> Growth: False\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Orphan compounds\n",
      "\n",
      "Orphan compounds are those compounds which are only associated with one reaction. They are either produced, or trying to be consumed. We need to add reaction(s) that complete the network of those compounds.\n",
      "\n",
      "You can change the maximum number of reactions that a compound is in to be considered an orphan (try increasing it to 2 or 3)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "orphan_reactions = PyFBA.gapfill.suggest_by_compound(compounds, reactions,\n",
      "                                                     reactions_to_run,\n",
      "                                                     max_reactions=1)\n",
      "added_reactions.append((\"orphans\", orphan_reactions))\n",
      "reactions_to_run.update(orphan_reactions)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, reactions_to_run,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {
      "scrolled": false
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The biomass reaction has a flux of 1000.0 --> Growth: True\n"
       ]
      }
     ],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Trimming the model\n",
      "Now that the model has been shown to grow on ArgonneLB media after several gap-fill iterations, we should trim down the reactions to only the required reactions necessary to observe growth."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "reqd_additional = set()\n",
      "\n",
      "# Begin loop through all gap-filled reactions\n",
      "while added_reactions:\n",
      "    ori = copy.copy(original_reactions_to_run)\n",
      "    ori.update(reqd_additional)\n",
      "    # Test next set of gap-filled reactions\n",
      "    # Each set is based on a method described above\n",
      "    how, new = added_reactions.pop()\n",
      "    sys.stderr.write(\"Testing reactions from {}\".format(how))\n",
      "    \n",
      "    # Get all the other gap-filled reactions we need to add\n",
      "    for tple in added_reactions:\n",
      "        ori.update(tple[1])\n",
      "    \n",
      "    # Use minimization function to determine the minimal\n",
      "    # set of gap-filled reactions from the current method\n",
      "    new_essential = PyFBA.gapfill.minimize_additional_reactions(ori, new, compounds,\n",
      "                                                                reactions, media,\n",
      "                                                                biomass_equation)\n",
      "    sys.stderr.write(\"Saved {} reactions from {}\"})\n",
      "    # Record the method used to determine\n",
      "    # how the reaction was gap-filled\n",
      "    for new_r in new_essential:\n",
      "        reactions[new_r].is_gapfilled = True\n",
      "        reactions[new_r].gapfill_method = how\n",
      "    reqd_additional.update(new_essential)\n",
      "\n",
      "# Combine old and new reactions\n",
      "all_reactions = original_reactions_to_run.union(reqd_additional)"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "status, value, growth = PyFBA.fba.run_fba(compounds, reactions, all_reactions,\n",
      "                                          media, biomass_equation)\n",
      "print(\"The biomass reaction has a flux of {} --> Growth: {}\".format(value, growth))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Other gap-filling techniques\n",
      "\n",
      "Besides those methods we have described above, listed here are other methods that can be used to gap-fill your model. This list will continue to grow over time as we create new techniques to identify reactions and compounds that should be added to your model."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "collapsed": true
     },
     "source": [
      "###Probable reactions\n",
      "Probable reactions are those reactions whose probability is based on whether there is a protein associated with the reaction and if the reaction's compounds are currently present in the model. Above a certain probability threshold, those reactions will be added to the model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "probable_reactions = PyFBA.gapfill.compound_probability(reactions, reactions_to_run,\n",
      "                                                        cutoff=0, rxn_with_protiens=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    }
   ],
   "metadata": {}
  }
 ]
}