8_custom_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Using a custom model data structure"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This notebooks shows how to utilize the generic accessors and modification\n",
    "functions in COBREXA.jl to run the analysis on any custom model type. We will\n",
    "create a simple dictionary-style structure that describes the model, allow\n",
    "COBREXA to run a FVA on it, and create a simple reaction-removing\n",
    "modification."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "First, let's define a very simple stoichiometry-only structure for the model:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using COBREXA\n",
    "\n",
    "mutable struct MyReaction\n",
    "    max_rate::Float64 # maximum absolute conversion rate\n",
    "    stoi::Dict{String,Float64} # stoichimetry of the reaction\n",
    "\n",
    "    MyReaction() = new(0.0, Dict{String,Float64}())\n",
    "end\n",
    "\n",
    "mutable struct MyModel <: MetabolicModel\n",
    "    optimization_target::String # the \"objective\" reaction name\n",
    "    reactions::Dict{String,MyReaction} # dictionary of reactions\n",
    "\n",
    "    MyModel() = new(\"\", Dict{String,MyReaction}())\n",
    "    MyModel(o, r) = new(o, r)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "With this, we can start defining the accessors:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "COBREXA.n_reactions(m::MyModel) = length(m.reactions)\n",
    "COBREXA.reactions(m::MyModel) = sort(collect(keys(m.reactions)))"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Metabolites are defined only very implicitly, so let's just make a function\n",
    "that collects all names. `n_metabolites` can be left at the default\n",
    "definition that just measures the output of `metabolites`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function COBREXA.metabolites(m::MyModel)\n",
    "    mets = Set{String}()\n",
    "    for (_, r) in m.reactions\n",
    "        for (m, _) in r.stoi\n",
    "            push!(mets, m)\n",
    "        end\n",
    "    end\n",
    "    return sort(collect(mets))\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, the extraction of the linear model. Remember the order of element in the\n",
    "vectors must match the order in the output of `reactions` and\n",
    "`metabolites`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SparseArrays\n",
    "\n",
    "function COBREXA.bounds(m::MyModel)\n",
    "    max_rates = [m.reactions[r].max_rate for r in reactions(m)]\n",
    "    (sparse(-max_rates), sparse(max_rates))\n",
    "end\n",
    "\n",
    "function COBREXA.objective(m::MyModel)\n",
    "    if m.optimization_target in keys(m.reactions)\n",
    "        c = spzeros(n_reactions(m))\n",
    "        c[first(indexin([m.optimization_target], reactions(m)))] = 1.0\n",
    "        c\n",
    "    else\n",
    "        throw(\n",
    "            DomainError(\n",
    "                m.optimization_target,\n",
    "                \"The target reaction for flux optimization not found\",\n",
    "            ),\n",
    "        )\n",
    "    end\n",
    "end\n",
    "\n",
    "function COBREXA.stoichiometry(m::MyModel)\n",
    "    sparse([\n",
    "        get(m.reactions[rxn].stoi, met, 0.0) for met in metabolites(m), rxn in reactions(m)\n",
    "    ])\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now the model is complete! We can generate a random one and see how it\n",
    "performs"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "import Random\n",
    "Random.seed!(123)\n",
    "\n",
    "rxn_names = [\"Reaction $i\" for i = 'A':'Z'];\n",
    "metabolite_names = [\"Metabolite $i\" for i = 1:20];\n",
    "\n",
    "m = MyModel();\n",
    "for i in rxn_names\n",
    "    m.reactions[i] = MyReaction()\n",
    "end\n",
    "\n",
    "for i = 1:50\n",
    "    rxn = rand(rxn_names)\n",
    "    met = rand(metabolite_names)\n",
    "    m.reactions[rxn].stoi[met] = rand([-3, -2, -1, 1, 2, 3])\n",
    "    m.reactions[rxn].max_rate = rand()\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's see what the model looks like now:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Metabolic model of type Main.##273.MyModel\n\n⠐⠀⠀⡈⠀⠀⠈⠁⠂⠀⠠⠀⠀\n⠀⠄⠀⡀⠀⠀⠁⠒⡀⠀⠐⠈⠀\n⠂⠀⠐⠀⠆⠠⠡⢀⠠⠐⡀⠀⢀\n⠀⠨⠁⡂⠀⠘⠀⠃⠀⠀⠀⠀⠀\n⠡⠉⠐⠀⢀⠀⠈⠀⢄⠈⠀⠀⡄\nNumber of reactions: 26\nNumber of metabolites: 20\n"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "m"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can run most of the standard function on the model data right away:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Dict{String, Float64} with 26 entries:\n  \"Reaction E\" => 0.000322931\n  \"Reaction F\" => -0.0\n  \"Reaction T\" => 0.174544\n  \"Reaction G\" => 0.0\n  \"Reaction B\" => -4.55624e-12\n  \"Reaction J\" => -0.102165\n  \"Reaction K\" => -0.0\n  \"Reaction P\" => 0.00193759\n  \"Reaction R\" => -0.102165\n  \"Reaction U\" => 0.000443018\n  \"Reaction Y\" => -1.849e-11\n  \"Reaction N\" => 0.56429\n  \"Reaction Z\" => 0.000350518\n  \"Reaction D\" => 0.0\n  \"Reaction W\" => -0.0\n  \"Reaction M\" => -0.0\n  \"Reaction X\" => 0.0\n  \"Reaction C\" => 0.0\n  \"Reaction A\" => 0.30446\n  ⋮            => ⋮"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "using Tulip\n",
    "m.optimization_target = \"Reaction A\"\n",
    "flux_balance_analysis_dict(m, Tulip.Optimizer)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "To be able to use the model conveniently in functions such as\n",
    "`screen`, you usually want to be able to easily specify the\n",
    "modifications. In this example, we enable use of\n",
    "`with_removed_reactions` by overloading the internal\n",
    "`remove_reactions` for this specific model type:\n",
    "\n",
    "We need to make an as-shallow-as-possible copy of the model that allows us to\n",
    "remove the reactions without breaking the original model."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function COBREXA.remove_reactions(m::MyModel, rxns::AbstractVector{String})\n",
    "    m = MyModel(m.optimization_target, copy(m.reactions))\n",
    "    delete!.(Ref(m.reactions), rxns)\n",
    "    return m\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "The screening is ready now!"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "25-element Vector{Pair{String, Float64}}:\n \"Reaction B\" => 0.3044596278261141\n \"Reaction C\" => 0.3044596423819445\n \"Reaction D\" => 0.3044596423819445\n \"Reaction E\" => 0.303975242543509\n \"Reaction F\" => 0.3044596423819445\n \"Reaction G\" => 0.3044596424143211\n \"Reaction H\" => 0.04135032969093828\n \"Reaction I\" => 0.03084320823855454\n \"Reaction J\" => 0.26359370990717723\n \"Reaction K\" => 0.3044596423819445\n              ⋮\n \"Reaction R\" => 0.263593709330887\n \"Reaction S\" => 0.3044596423819445\n \"Reaction T\" => 0.10281069417903423\n \"Reaction U\" => 0.3044596424137387\n \"Reaction V\" => 0.3044596423819445\n \"Reaction W\" => 0.3044596423819445\n \"Reaction X\" => 0.3044596423819445\n \"Reaction Y\" => 0.3044596260552588\n \"Reaction Z\" => 0.3044596424137387"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "reactions_to_remove = (\"Reaction $i\" for i = 'B':'Z')\n",
    "\n",
    "reactions_to_remove .=> screen_variants(\n",
    "    m,\n",
    "    [[with_removed_reactions([r])] for r in reactions_to_remove],\n",
    "    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)[\"Reaction A\"],\n",
    ")"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "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.6.0"
  },
  "kernelspec": {
   "name": "julia-1.6",
   "display_name": "Julia 1.6.0",
   "language": "julia"
  }
 },
 "nbformat": 4
}