Newer
Older
{"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.5","display_name":"Julia 1.5","language":"julia"}},"nbformat":4}