1 - Ground-truth generation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shutil\n",
    "import subprocess\n",
    "import os\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import biolqm\n",
    "\n",
    "from colomoto.minibn import BooleanNetwork\n",
    "from merrin.bnet2flexflux import bn2flexflux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models = [\"covert\", \"covert_no_reg\"]\n",
    "pipeline = \"nbi\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metabolic_network = \"data/{}/metabolic_network.xml\"\n",
    "regulatory_network = \"data/{}/regulatory_network.{}\"\n",
    "\n",
    "FlexFlux_output = \"data/{}/result/pipeline_{}/FlexFlux\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_sbmlqual(bnet_file, sbmlqual_file, regulators):\n",
    "    bnet = BooleanNetwork.load(bnet_file)\n",
    "    bn2flexflux(bnet, sbmlqual_file, regulators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulators = {'RPO2', 'RPb', 'RPcl', 'RPh'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for model in models:\n",
    "    bnet_file = regulatory_network.format(model, 'bnet')\n",
    "    reg_file = regulatory_network.format(model, 'sbml')\n",
    "    write_sbmlqual(bnet_file, reg_file, regulators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lqm = biolqm.load(regulatory_network.format(models[0], 'sbml'))\n",
    "bn = biolqm.to_minibn(lqm)\n",
    "bn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constraints = \"\"\"\\\n",
    "obj : MAX(Growth)\n",
    "\n",
    "\n",
    "Tc1 0 10.5\n",
    "Tc2 0 10.5\n",
    "Td 0 12\n",
    "Te 0 12\n",
    "Tf 0 5\n",
    "Th 0 5\n",
    "To2 0 15\n",
    "\"\"\"\n",
    "\n",
    "initParams = {\n",
    "    \"biomassReaction\": \"Growth\",\n",
    "    \"timeStep\": \"0.01\",\n",
    "    \"nSteps\": 300,\n",
    "    \"initBiomass\": \"0.01\"\n",
    "}\n",
    "\n",
    "initConcentrations = {\n",
    "    \"Carbon1\": 10,\n",
    "    \"Carbon2\": 10,\n",
    "    \"Hext\": 0,\n",
    "    \"Dext\": 0,\n",
    "    \"Eext\": 0,\n",
    "    \"Fext\": 0,\n",
    "    \"Oxygen\": 100,\n",
    "    \"Biomass\": 0\n",
    "}\n",
    "\n",
    "default = {\n",
    "    \"constraints\": constraints,\n",
    "    \"initParams\": initParams,\n",
    "    \"initConcentrations\": initConcentrations,\n",
    "}\n",
    "\n",
    "simulations = {\n",
    "    \"fig5\": default,\n",
    "    \"fig6\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Oxygen\": 5,\n",
    "    }},\n",
    "    \"fig7\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Hext\": 2,\n",
    "    }},\n",
    "    \"fig8\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Carbon2\": 5,\n",
    "        \"Hext\": 10,\n",
    "    }},\n",
    "    \"fig9\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 1,\n",
    "        \"Hext\": 5,\n",
    "        \"Fext\": 0.1,\n",
    "    }},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatically set initial value of regulatory nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def auto_initReg(x0):\n",
    "    return {i:v for i,v in bn(x0).items() if v == 1 and i in regulators}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cfg in simulations.values():\n",
    "    initReg = auto_initReg(cfg[\"initConcentrations\"])\n",
    "    if initReg:\n",
    "        cfg[\"initRegStates\"] = initReg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def met4j_dynamicRFBA(mn, rn, cons, init, output):\n",
    "    subprocess.run([\"java\", \"-Dlog4j.configuration=\",\n",
    "        \"-jar\", \"/opt/met4j-dynamicRFBA.jar\",\n",
    "        \"-sol\", \"GLPK\", \n",
    "        \"-s\", mn, \"-reg\", rn,\n",
    "        \"-cons\", cons, \"-init\", init, \"-out\", output], check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_constraints(cfg, filename):\n",
    "    with open(filename, \"w\") as fp:\n",
    "        fp.write(cfg[\"constraints\"])\n",
    "\n",
    "def _write_state(fp, sect, state):\n",
    "    fp.write(f\"{sect}\\n\")\n",
    "    for node, value in state.items():\n",
    "        fp.write(f\"{node}\\t{value}\\n\")\n",
    "        \n",
    "def write_initial(cfg, filename):\n",
    "    with open(filename, \"w\") as fp:\n",
    "        for sect, value in cfg[\"initParams\"].items():\n",
    "            fp.write(f\"{sect}\\n{value}\\n\\n\")\n",
    "        sect = \"initConcentrations\"\n",
    "        _write_state(fp, sect, cfg[sect])\n",
    "        sect = \"initRegStates\"\n",
    "        if sect in cfg:\n",
    "            fp.write(\"\\n\")\n",
    "            _write_state(fp, sect, cfg[sect])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_simulation(out):\n",
    "    df = pd.read_csv(out, sep=\"\\t\")\n",
    "    df.set_index(\"Time\", inplace=True)\n",
    "    return df\n",
    "\n",
    "def run_simulations(mn_file, reg_file, output_dir, output_tmpl):\n",
    "    output = {}\n",
    "\n",
    "    for name, cfg in tqdm(simulations.items()):\n",
    "        out = os.path.join(output_dir, output_tmpl.format(name))\n",
    "        cons = os.path.join(output_dir, f\"constraints-{name}.txt\")\n",
    "        init = os.path.join(output_dir, f\"initialParameters-{name}.txt\")\n",
    "        write_constraints(cfg, cons)\n",
    "        write_initial(cfg, init)\n",
    "        met4j_dynamicRFBA(mn_file, reg_file, cons, init, out)\n",
    "        output[name] = read_simulation(out)\n",
    "        \n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_tmpl = \"out_times-{}.csv\"\n",
    "model_outputs = {}\n",
    "\n",
    "for model in models:\n",
    "    output_dir = FlexFlux_output.format(model, pipeline)\n",
    "    if os.path.isdir(output_dir):\n",
    "        shutil.rmtree(output_dir)\n",
    "    os.makedirs(output_dir)\n",
    "\n",
    "    mn_file = metabolic_network.format(model)\n",
    "    reg_file = regulatory_network.format(model, 'sbml')\n",
    "\n",
    "    model_outputs[model] = run_simulations(mn_file, reg_file, output_dir, output_tmpl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}