HP_gem_refinement.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": "\n(function(global) {\n  function now() {\n    return new Date();\n  }\n\n  var force = \"1\";\n\n  if (typeof (window._bokeh_onload_callbacks) === \"undefined\" || force !== \"\") {\n    window._bokeh_onload_callbacks = [];\n    window._bokeh_is_loading = undefined;\n  }\n\n  function run_callbacks() {\n    window._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n    delete window._bokeh_onload_callbacks\n    console.info(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(js_urls, callback) {\n    window._bokeh_onload_callbacks.push(callback);\n    if (window._bokeh_is_loading > 0) {\n      console.log(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    window._bokeh_is_loading = js_urls.length;\n    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      var s = document.createElement('script');\n      s.src = url;\n      s.async = false;\n      s.onreadystatechange = s.onload = function() {\n        window._bokeh_is_loading--;\n        if (window._bokeh_is_loading === 0) {\n          console.log(\"Bokeh: all BokehJS libraries loaded\");\n          run_callbacks()\n        }\n      };\n      s.onerror = function() {\n        console.warn(\"failed to load library \" + url);\n      };\n      console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.getElementsByTagName(\"head\")[0].appendChild(s);\n    }\n  };\n\n  var js_urls = ['https://cdn.pydata.org/bokeh/release/bokeh-0.12.1.min.js', 'https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.1.min.js', 'https://cdn.pydata.org/bokeh/release/bokeh-compiler-0.12.1.min.js'];\n\n  var inline_js = [\n    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\n    \n    function(Bokeh) {\n      \n      Bokeh.$(\"#307b97a9-6384-4ab9-9ece-ec737b084d78\").text(\"BokehJS successfully loaded\");\n      var kernel = Jupyter.notebook.kernel\n      if (kernel.execute !== undefined) {\n          kernel.execute(\"import bokeh.io; bokeh.io._nb_loaded = True\");\n      }\n\n    },\n    function(Bokeh) {\n      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.12.1.min.css\");\n      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.12.1.min.css\");\n      console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.1.min.css\");\n      Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.1.min.css\");\n    }\n  ];\n\n  function run_inline_js() {\n    for (var i = 0; i < inline_js.length; i++) {\n      inline_js[i](window.Bokeh);\n    }\n  }\n\n  if (window._bokeh_is_loading === 0) {\n    console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(js_urls, function() {\n      console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(this));"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cobra.test\n",
    "import cobra\n",
    "from cobra import Model, Reaction, Metabolite\n",
    "import os\n",
    "import pandas as pd\n",
    "from Bio import SeqIO\n",
    "import pickle\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "HP_model = cobra.io.read_sbml_model(\"draft_model/HP_draft_model.xml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Add missing ATP reaction </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nsummary = model_coeli.metabolites.atp_c.summary()\\nsummary = summary.to_frame()\\nrxns = summary.index.get_level_values(1).to_list()\\n\\ngene_ids = []\\n\\nfor rxn in rxns:\\n    print(\"Genes\")\\n    print(\"-----\")\\n    for x in model_coeli.genes:\\n        associated_ids = [i.id for i in x.reactions if i.id == rxn]\\n        if len(associated_ids)>0:\\n            gene_ids.append([x.id, associated_ids[0]])\\n            print(\"%s is associated with reactions: %s\" %\\n                  (x.id, associated_ids))\\n    print(\"\")\\n\\nprint(gene_ids)\\n'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_coeli = pickle.load(open('template_model/sco_model.p','rb'))\n",
    "'''\n",
    "summary = model_coeli.metabolites.atp_c.summary()\n",
    "summary = summary.to_frame()\n",
    "rxns = summary.index.get_level_values(1).to_list()\n",
    "\n",
    "gene_ids = []\n",
    "\n",
    "for rxn in rxns:\n",
    "    print(\"Genes\")\n",
    "    print(\"-----\")\n",
    "    for x in model_coeli.genes:\n",
    "        associated_ids = [i.id for i in x.reactions if i.id == rxn]\n",
    "        if len(associated_ids)>0:\n",
    "            gene_ids.append([x.id, associated_ids[0]])\n",
    "            print(\"%s is associated with reactions: %s\" %\n",
    "                  (x.id, associated_ids))\n",
    "    print(\"\")\n",
    "\n",
    "print(gene_ids)\n",
    "'''"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "##### check sco genes' EC number\n",
    "sco = SeqIO.read('template_model/NC_003888.3.gb,'genbank')\n",
    "ec_list = []\n",
    "for gene,rxn in gene_ids:\n",
    "    for seq_feature in sco.features:\n",
    "        if seq_feature.type=='CDS':\n",
    "            if seq_feature.qualifiers.get('locus_tag'):\n",
    "                locus_tag = seq_feature.qualifiers.get('locus_tag')[0]\n",
    "                if locus_tag == gene and seq_feature.qualifiers.get('EC_number'):\n",
    "                    ec_number = seq_feature.qualifiers.get('EC_number')[0]\n",
    "                    if ec_number not in ec_list:\n",
    "                        ec_list.append(ec_number)\n",
    "                    print(locus_tag, seq_feature.qualifiers.get('EC_number'))\n",
    "print(ec_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "reaction_atps4r = model_coeli.reactions.ATPS4r\n",
    "HP_model.add_reactions([reaction_atps4r])\n",
    "\n",
    "#check corresponding sco genes from template_model/temp_target_BBH_dict.txt\n",
    "\n",
    "HP_model.reactions.ATPS4r.gene_reaction_rule = 'NJMEEBPL_03421 and NJMEEBPL_03423 and NJMEEBPL_03416 and NJMEEBPL_03417 and NJMEEBPL_03418 and NJMEEBPL_03419 and NJMEEBPL_03420'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "reaction_ACCOAC = model_coeli.reactions.ACCOAC\n",
    "HP_model.add_reactions([reaction_ACCOAC])\n",
    "HP_model.reactions.ACCOAC.gene_reaction_rule = '( NJMEEBPL_07611 or ( NJMEEBPL_03221 and NJMEEBPL_03220 ) ) and ( NJMEEBPL_06272 or NJMEEBPL_05660 or NJMEEBPL_09895 )'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Solution 0.074 at 0x7f49225be610>\n",
      "0.07407495361456679\n"
     ]
    }
   ],
   "source": [
    "solution = HP_model.optimize()\n",
    "print(solution)\n",
    "print(solution.objective_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"2\" halign=\"left\">IN_FLUXES</th>\n",
       "      <th colspan=\"2\" halign=\"left\">OUT_FLUXES</th>\n",
       "      <th colspan=\"2\" halign=\"left\">OBJECTIVES</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "      <th>ID</th>\n",
       "      <th>FLUX</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>o2_e</td>\n",
       "      <td>1.794722</td>\n",
       "      <td>h2o_e</td>\n",
       "      <td>3.625673</td>\n",
       "      <td>BIOMASS_SCO</td>\n",
       "      <td>0.074075</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>glc__D_e</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>co2_e</td>\n",
       "      <td>1.872247</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>nh4_e</td>\n",
       "      <td>0.658388</td>\n",
       "      <td>h_e</td>\n",
       "      <td>0.606916</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>pi_e</td>\n",
       "      <td>0.075283</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<cobra.core.summary.model_summary.ModelSummary at 0x7f49225beb50>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HP_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>PERCENT</th>\n",
       "      <th>FLUX</th>\n",
       "      <th>REACTION_STRING</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RXN_STAT</th>\n",
       "      <th>ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">PRODUCING</th>\n",
       "      <th>ATPS4r</th>\n",
       "      <td>78.584889</td>\n",
       "      <td>7.264846</td>\n",
       "      <td>adp_c + 4.0 h_e + pi_c &lt;=&gt; atp_c + h2o_c + 3.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PGK</th>\n",
       "      <td>13.963353</td>\n",
       "      <td>1.290854</td>\n",
       "      <td>3pg_c + atp_c &lt;=&gt; 13dpg_c + adp_c</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SUCOAS</th>\n",
       "      <td>5.538858</td>\n",
       "      <td>0.512044</td>\n",
       "      <td>atp_c + coa_c + succ_c &lt;=&gt; adp_c + pi_c + succ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NDPK9</th>\n",
       "      <td>1.912900</td>\n",
       "      <td>0.176840</td>\n",
       "      <td>atp_c + idp_c &lt;=&gt; adp_c + itp_c</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">CONSUMING</th>\n",
       "      <th>BIOMASS_SCO</th>\n",
       "      <td>60.745356</td>\n",
       "      <td>5.615656</td>\n",
       "      <td>0.004309 10fthf_c + 0.019494 15dap_c + 0.01169...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ATPM</th>\n",
       "      <td>28.665432</td>\n",
       "      <td>2.650000</td>\n",
       "      <td>atp_c + h2o_c --&gt; adp_c + h_c + pi_c</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GLNS</th>\n",
       "      <td>1.534612</td>\n",
       "      <td>0.141868</td>\n",
       "      <td>atp_c + glu__L_c + nh4_c --&gt; adp_c + gln__L_c ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NDPK2</th>\n",
       "      <td>0.834255</td>\n",
       "      <td>0.077123</td>\n",
       "      <td>atp_c + udp_c &lt;=&gt; adp_c + utp_c</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CBPS</th>\n",
       "      <td>0.792308</td>\n",
       "      <td>0.073246</td>\n",
       "      <td>2.0 atp_c + gln__L_c + h2o_c + hco3_c --&gt; 2.0 ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<cobra.core.summary.metabolite_summary.MetaboliteSummary at 0x7f4922565a90>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HP_model.metabolites.atp_c.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Add rapamycin producing pathway </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_metabolites(io, sheet_name, model):\n",
    "    df_metab = pd.read_excel(io, sheet_name)\n",
    "    for i in range(df_metab.shape[0]):\n",
    "        metab = Metabolite(df_metab.ID[i], \n",
    "                           formula = df_metab.Formula[i],\n",
    "                           name = df_metab.Name[i],\n",
    "                           compartment = df_metab.Compartment[i])\n",
    "        model.add_metabolites([metab])\n",
    "        print(metab.id,' is added to the model.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_reactions(io, sheet_name, model):\n",
    "    df_reac = pd.read_excel(io, sheet_name)\n",
    "    for i in range(df_reac.shape[0]):\n",
    "        reac = Reaction(df_reac.ID[i])\n",
    "        reac.name = df_reac.Name[i]\n",
    "        model.add_reactions([reac])\n",
    "        reac.build_reaction_from_string(df_reac.Equation[i], fwd_arrow = '=>', reversible_arrow = '<=>')\n",
    "        \n",
    "        if df_reac.Reversible[i] == 'Y':\n",
    "            reac.lower_bound = -1000.0\n",
    "            reac.upper_bound = 1000.0\n",
    "        else:\n",
    "            reac.lower_bound = 0.0\n",
    "            reac.upper_bound = 1000.0\n",
    "        \n",
    "        \n",
    "        print(reac.id, ' is added to the model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dhchc_c  is added to the model.\n",
      "pipecol__L_c  is added to the model.\n",
      "rap_c  is added to the model.\n"
     ]
    }
   ],
   "source": [
    "add_metabolites('inputs/metabolite_reaction.xlsx', 0, HP_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DHCHCS  is added to the model\n",
      "MNXR113634  is added to the model\n",
      "Rapamycin  is added to the model\n"
     ]
    }
   ],
   "source": [
    "add_reactions('inputs/metabolite_reaction.xlsx', 1, HP_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong><em>Optimal</em> solution with objective value 0.142</strong><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fluxes</th>\n",
       "      <th>reduced_costs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ACTS10</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-2.170680e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACTS13</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-2.170680e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACTS14</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACTS15</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACTS16</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACCOAC</th>\n",
       "      <td>0.995594</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DHCHCS</th>\n",
       "      <td>0.142228</td>\n",
       "      <td>-4.857226e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MNXR113634</th>\n",
       "      <td>0.142228</td>\n",
       "      <td>-2.428613e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Rapamycin</th>\n",
       "      <td>0.142228</td>\n",
       "      <td>1.172673e-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DM_rap_c</th>\n",
       "      <td>0.142228</td>\n",
       "      <td>4.033232e-16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2040 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "<Solution 0.142 at 0x7f48ff822890>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HP_model.reactions.EX_glc__D_e.lower_bound = -10.0\n",
    "#HP_model.reactions.EX_glyc_e.lower_bound = -5.0\n",
    "#HP_model.reactions.BIOMASS_SCO.lower_bound = 0.75*1.309\n",
    "DM_rap = HP_model.add_boundary(HP_model.metabolites.rap_c, type='demand')\n",
    "HP_model.objective = HP_model.reactions.Rapamycin\n",
    "HP_model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "cobra.io.write_sbml_model(HP_model,\"final_model/HP_model_final.xml\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gmsm",
   "language": "python",
   "name": "gmsm"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}