build based on 9497ad2

dev/.documenter-siteinfo.json

@@ -1 +1 @@
-{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2024-09-20T13:29:08","documenter_version":"1.7.0"}}
+{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2024-09-24T19:23:48","documenter_version":"1.7.0"}}

dev/examples/05g-gapfilling.ipynb

@@ -0,0 +1,356 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Gap filling\n",
+    "\n",
+    "Gapfilling is used to find easiest additions to the models that would make\n",
+    "them feasible and capable of growth.\n",
+    "\n",
+    "Typically, an infeasible model (\"with gaps\") is used together with an\n",
+    "universal model (which contains \"everything\"), and the algorithm attempts to\n",
+    "find the minimal amount of reactions from the universal model that make the\n",
+    "gappy model happy. In turn, the gapfilling optimization problem becomes a\n",
+    "MILP.\n",
+    "\n",
+    "Gapfilling is sometimes used to produce \"viable\" genome-scale\n",
+    "reconstructions from partial ones, but without additional manual intervention\n",
+    "the gapfilling results are almost never biologically valid. A good use of\n",
+    "gapfilling is to find problems in a model that cause infeasibility as\n",
+    "follows: First the modeller makes a set of (unrealistic) universal reactions\n",
+    "that supply or remove metabolites, and after gapfilling, metabolites that had\n",
+    "to be supplied or removed to make the model feasible mark possible problems,\n",
+    "thus guiding the modeller towards correct solution."
+   ],
+   "metadata": {}
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "We will use a partially crippled *E. coli* toy model and see the minimal\n",
+    "amount of reactions that may save it."
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[ Info: using cached `e_coli_core.json'\n"
+     ]
+    },
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)"
+     },
+     "metadata": {},
+     "execution_count": 1
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "using COBREXA\n",
+    "\n",
+    "download_model(\n",
+    "    \"http://bigg.ucsd.edu/static/models/e_coli_core.json\",\n",
+    "    \"e_coli_core.json\",\n",
+    "    \"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8\",\n",
+    ")\n",
+    "\n",
+    "import JSONFBCModels, HiGHS\n",
+    "model = load_model(\"e_coli_core.json\")"
+   ],
+   "metadata": {},
+   "execution_count": 1
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "First, let's produce an infeasible model:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [],
+   "cell_type": "code",
+   "source": [
+    "import AbstractFBCModels.CanonicalModel as CM\n",
+    "infeasible_model = convert(CM.Model, model)\n",
+    "\n",
+    "for rxn in [\"TALA\", \"PDH\", \"PGI\", \"PYK\"]\n",
+    "    infeasible_model.reactions[rxn].lower_bound = 0.0\n",
+    "    infeasible_model.reactions[rxn].upper_bound = 0.0\n",
+    "end"
+   ],
+   "metadata": {},
+   "execution_count": 2
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "After removing the above reactions, the model will fail to solve:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "nothing\n"
+     ]
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "flux_balance_analysis(infeasible_model, optimizer = HiGHS.Optimizer) |> println"
+   ],
+   "metadata": {},
+   "execution_count": 3
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "To avoid very subtle semantic issues, we are going to remove the ATP\n",
+    "maintenance pseudoreaction from the universal model:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "Dict{String, AbstractFBCModels.CanonicalModel.Reaction} with 94 entries:\n  \"ACALD\"       => Reaction(\"Acetaldehyde dehydrogenase (acetylating)\", -1000.0…\n  \"PTAr\"        => Reaction(\"Phosphotransacetylase\", -1000.0, 1000.0, Dict(\"coa…\n  \"ALCD2x\"      => Reaction(\"Alcohol dehydrogenase (ethanol)\", -1000.0, 1000.0,…\n  \"PDH\"         => Reaction(\"Pyruvate dehydrogenase\", 0.0, 1000.0, Dict(\"coa_c\"…\n  \"PYK\"         => Reaction(\"Pyruvate kinase\", 0.0, 1000.0, Dict(\"adp_c\"=>-1.0,…\n  \"CO2t\"        => Reaction(\"CO2 transporter via diffusion\", -1000.0, 1000.0, D…\n  \"EX_nh4_e\"    => Reaction(\"Ammonia exchange\", -1000.0, 1000.0, Dict(\"nh4_e\"=>…\n  \"MALt2_2\"     => Reaction(\"Malate transport via proton symport (2 H)\", 0.0, 1…\n  \"CS\"          => Reaction(\"Citrate synthase\", 0.0, 1000.0, Dict(\"coa_c\"=>1.0,…\n  \"PGM\"         => Reaction(\"Phosphoglycerate mutase\", -1000.0, 1000.0, Dict(\"3…\n  \"TKT1\"        => Reaction(\"Transketolase\", -1000.0, 1000.0, Dict(\"g3p_c\"=>1.0…\n  \"EX_mal__L_e\" => Reaction(\"L-Malate exchange\", 0.0, 1000.0, Dict(\"mal__L_e\"=>…\n  \"ACONTa\"      => Reaction(\"Aconitase (half-reaction A, Citrate hydro-lyase)\",…\n  \"EX_pi_e\"     => Reaction(\"Phosphate exchange\", -1000.0, 1000.0, Dict(\"pi_e\"=…\n  \"GLNS\"        => Reaction(\"Glutamine synthetase\", 0.0, 1000.0, Dict(\"adp_c\"=>…\n  \"ICL\"         => Reaction(\"Isocitrate lyase\", 0.0, 1000.0, Dict(\"succ_c\"=>1.0…\n  \"EX_o2_e\"     => Reaction(\"O2 exchange\", -1000.0, 1000.0, Dict(\"o2_e\"=>-1.0),…\n  \"FBA\"         => Reaction(\"Fructose-bisphosphate aldolase\", -1000.0, 1000.0, …\n  \"EX_gln__L_e\" => Reaction(\"L-Glutamine exchange\", 0.0, 1000.0, Dict(\"gln__L_e…\n  ⋮             => ⋮"
+     },
+     "metadata": {},
+     "execution_count": 4
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "universal_model = convert(CM.Model, model)\n",
+    "delete!(universal_model.reactions, \"ATPM\")"
+   ],
+   "metadata": {},
+   "execution_count": 4
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Making the model feasible with a minimal set of reactions"
+   ],
+   "metadata": {}
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Which of the reactions we have to fill back to get the model working again?\n",
+    "First, let's run `gap_filling_analysis` to get a solution for a\n",
+    "system that implements the reaction patching:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "ConstraintTrees.Tree{Float64} with 6 elements:\n  :fill_flags            => ConstraintTrees.Tree{Float64}(#= 94 elements =#)\n  :n_filled              => 1.0\n  :stoichiometry         => ConstraintTrees.Tree{Float64}(#= 72 elements =#)\n  :system                => ConstraintTrees.Tree{Float64}(#= 3 elements =#)\n  :universal_flux_bounds => ConstraintTrees.Tree{Float64}(#= 94 elements =#)\n  :universal_fluxes      => ConstraintTrees.Tree{Float64}(#= 94 elements =#)"
+     },
+     "metadata": {},
+     "execution_count": 5
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "x = gap_filling_analysis(\n",
+    "    infeasible_model,\n",
+    "    universal_model,\n",
+    "    0.05,\n",
+    "    optimizer = HiGHS.Optimizer,\n",
+    ")"
+   ],
+   "metadata": {},
+   "execution_count": 5
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "The reactions that had to be re-added can be found from the `fill_flags`:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "1-element Vector{Symbol}:\n :TALA"
+     },
+     "metadata": {},
+     "execution_count": 6
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "filled_reactions = [k for (k, v) in x.fill_flags if v != 0]"
+   ],
+   "metadata": {},
+   "execution_count": 6
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "If we want to try to generate another solution, we have to explicitly ask the\n",
+    "system to avoid the previous solution. That is done via setting the argument\n",
+    "`known_fill`. We can also set the `max_cost` to avoid finding too benevolent\n",
+    "solutions:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "1-element Vector{Symbol}:\n :PGI"
+     },
+     "metadata": {},
+     "execution_count": 7
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "x2 = gap_filling_analysis(\n",
+    "    infeasible_model,\n",
+    "    universal_model,\n",
+    "    0.05,\n",
+    "    max_cost = 2.0,\n",
+    "    known_fills = [x.fill_flags],\n",
+    "    optimizer = HiGHS.Optimizer,\n",
+    ")\n",
+    "\n",
+    "other_filled_reactions = [k for (k, v) in x2.fill_flags if v != 0]"
+   ],
+   "metadata": {},
+   "execution_count": 7
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "## Model debugging: which metabolite is missing?\n",
+    "\n",
+    "Gap-filling is great for detecting various broken links and imbalances in\n",
+    "metabolic models. We show how to find the metabolites are causing the\n",
+    "imbalance for our \"broken\" E. coli model.\n",
+    "\n",
+    "First, we construct a few completely unnatural reactions that create/remove\n",
+    "the metabolites from/to nowhere:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [],
+   "cell_type": "code",
+   "source": [
+    "magic_model = convert(CM.Model, model)\n",
+    "empty!(magic_model.genes)\n",
+    "empty!(magic_model.reactions)\n",
+    "\n",
+    "for mid in keys(magic_model.metabolites)\n",
+    "    magic_model.reactions[mid] = CM.Reaction(\n",
+    "        lower_bound = -100.0,\n",
+    "        upper_bound = 100.0,\n",
+    "        stoichiometry = Dict(mid => 1.0),\n",
+    "    )\n",
+    "end"
+   ],
+   "metadata": {},
+   "execution_count": 8
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Gapfilling now points to the metabolites that need to be somehow taken care\n",
+    "of by the modeller in order for the model to become feasible:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "1-element Vector{Symbol}:\n :e4p_c"
+     },
+     "metadata": {},
+     "execution_count": 9
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "xm = gap_filling_analysis(infeasible_model, magic_model, 0.05, optimizer = HiGHS.Optimizer)\n",
+    "\n",
+    "blocking_metabolites = [k for (k, v) in xm.fill_flags if v != 0]"
+   ],
+   "metadata": {},
+   "execution_count": 9
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "We can also have a look at how much of a given metabolite was used to make\n",
+    "the model feasible again:"
+   ],
+   "metadata": {}
+  },
+  {
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": "0.6866985638297876"
+     },
+     "metadata": {},
+     "execution_count": 10
+    }
+   ],
+   "cell_type": "code",
+   "source": [
+    "xm.universal_fluxes[first(blocking_metabolites)]"
+   ],
+   "metadata": {},
+   "execution_count": 10
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "---\n",
+    "\n",
+    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
+   ],
+   "metadata": {}
+  }
+ ],
+ "nbformat_minor": 3,
+ "metadata": {
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.4"
+  },
+  "kernelspec": {
+   "name": "julia-1.9",
+   "display_name": "Julia 1.9.4",
+   "language": "julia"
+  }
+ },
+ "nbformat": 4
+}