{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a Model\n",
    "\n",
    "This simple example demonstrates how to create a model, create a reaction, and then add the reaction to the model.\n",
    "\n",
    "We'll use the '3OAS140' reaction from the STM_1.0 model:\n",
    "\n",
    "1.0 malACP[c] + 1.0 h[c] + 1.0 ddcaACP[c] $\\rightarrow$ 1.0 co2[c] + 1.0 ACP[c] + 1.0 3omrsACP[c]\n",
    "\n",
    "First, create the model and reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from cobra import Model, Reaction, Metabolite\n",
    "# Best practise: SBML compliant IDs\n",
    "cobra_model = Model('example_cobra_model')\n",
    "\n",
    "reaction = Reaction('3OAS140')\n",
    "reaction.name = '3 oxoacyl acyl carrier protein synthase n C140 '\n",
    "reaction.subsystem = 'Cell Envelope Biosynthesis'\n",
    "reaction.lower_bound = 0.  # This is the default\n",
    "reaction.upper_bound = 1000.  # This is the default\n",
    "reaction.objective_coefficient = 0. # this is the default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to create metabolites as well. If we were using an existing model, we could use get_by_id to get the apporpriate Metabolite objects instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ACP_c = Metabolite(\n",
    "    'ACP_c',\n",
    "    formula='C11H21N2O7PRS',\n",
    "    name='acyl-carrier-protein',\n",
    "    compartment='c')\n",
    "omrsACP_c = Metabolite(\n",
    "    '3omrsACP_c',\n",
    "    formula='C25H45N2O9PRS',\n",
    "    name='3-Oxotetradecanoyl-acyl-carrier-protein',\n",
    "    compartment='c')\n",
    "co2_c = Metabolite(\n",
    "    'co2_c',\n",
    "    formula='CO2',\n",
    "    name='CO2',\n",
    "    compartment='c')\n",
    "malACP_c = Metabolite(\n",
    "    'malACP_c',\n",
    "    formula='C14H22N2O10PRS',\n",
    "    name='Malonyl-acyl-carrier-protein',\n",
    "    compartment='c')\n",
    "h_c = Metabolite(\n",
    "    'h_c',\n",
    "    formula='H',\n",
    "    name='H',\n",
    "    compartment='c')\n",
    "ddcaACP_c = Metabolite(\n",
    "    'ddcaACP_c',\n",
    "    formula='C23H43N2O8PRS',\n",
    "    name='Dodecanoyl-ACP-n-C120ACP',\n",
    "    compartment='c')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding metabolites to a reaction requires using a dictionary of the metabolites and their stoichiometric coefficients. A group of metabolites can be added all at once, or they can be added one at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ddcaACP_c + h_c + malACP_c --> 3omrsACP_c + ACP_c + co2_c'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction.add_metabolites({malACP_c: -1.0,\n",
    "                          h_c: -1.0,\n",
    "                          ddcaACP_c: -1.0,\n",
    "                          co2_c: 1.0,\n",
    "                          ACP_c: 1.0,\n",
    "                          omrsACP_c: 1.0})\n",
    "\n",
    "\n",
    "reaction.reaction  # This gives a string representation of the reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gene_reaction_rule is a boolean representation of the gene requirements for this reaction to be active as described in [Schellenberger et al 2011 Nature Protocols 6(9):1290-307](http://dx.doi.org/doi:10.1038/nprot.2011.308). We will assign the gene reaction rule string, which will automatically create the corresponding gene objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({<Gene STM2378 at 0x7fada4592908>, <Gene STM1197 at 0x7fada45927f0>})"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reaction.gene_reaction_rule = '( STM2378 or STM1197 )'\n",
    "reaction.genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point in time, the model is still empty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 reactions initially\n",
      "0 metabolites initially\n",
      "0 genes initially\n"
     ]
    }
   ],
   "source": [
    "print('%i reactions initially' % len(cobra_model.reactions))\n",
    "print('%i metabolites initially' % len(cobra_model.metabolites))\n",
    "print('%i genes initially' % len(cobra_model.genes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will add the reaction to the model, which will also add all associated metabolites and genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 reaction\n",
      "6 metabolites\n",
      "2 genes\n"
     ]
    }
   ],
   "source": [
    "cobra_model.add_reaction(reaction)\n",
    "\n",
    "# Now there are things in the model\n",
    "print('%i reaction' % len(cobra_model.reactions))\n",
    "print('%i metabolites' % len(cobra_model.metabolites))\n",
    "print('%i genes' % len(cobra_model.genes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can iterate through the model objects to observe the contents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reactions\n",
      "---------\n",
      "3OAS140 : ddcaACP_c + h_c + malACP_c --> 3omrsACP_c + ACP_c + co2_c\n",
      "\n",
      "Metabolites\n",
      "-----------\n",
      "3omrsACP_c : C25H45N2O9PRS\n",
      "ddcaACP_c : C23H43N2O8PRS\n",
      "    ACP_c : C11H21N2O7PRS\n",
      "    co2_c : CO2\n",
      " malACP_c : C14H22N2O10PRS\n",
      "      h_c : H\n",
      "\n",
      "Genes\n",
      "-----\n",
      "STM2378 is associated with reactions: {3OAS140}\n",
      "STM1197 is associated with reactions: {3OAS140}\n"
     ]
    }
   ],
   "source": [
    "# Iterate through the the objects in the model\n",
    "print(\"Reactions\")\n",
    "print(\"---------\")\n",
    "for x in cobra_model.reactions:\n",
    "    print(\"%s : %s\" % (x.id, x.reaction))\n",
    "\n",
    "print(\"\")\n",
    "print(\"Metabolites\")\n",
    "print(\"-----------\")\n",
    "for x in cobra_model.metabolites:\n",
    "    print('%9s : %s' % (x.id, x.formula))\n",
    "\n",
    "print(\"\")\n",
    "print(\"Genes\")\n",
    "print(\"-----\")\n",
    "for x in cobra_model.genes:\n",
    "    associated_ids = (i.id for i in x.reactions)\n",
    "    print(\"%s is associated with reactions: %s\" %\n",
    "          (x.id, \"{\" + \", \".join(associated_ids) + \"}\"))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}