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.##290.MyModel\nsparse([1, 4, 8, 2, 5, 4, 16, 5, 15, 11 … 10, 8, 12, 6, 8, 13, 18, 18, 5, 15], [1, 1, 1, 2, 2, 3, 3, 4, 6, 7 … 19, 20, 20, 21, 22, 22, 23, 24, 26, 26], [1.0, -3.0, 1.0, -1.0, 2.0, 2.0, -2.0, -3.0, 1.0, 3.0 … 3.0, 3.0, 1.0, 2.0, 2.0, -1.0, 2.0, -1.0, 1.0, -3.0], 19, 26)\nNumber of reactions: 26\nNumber of metabolites: 19\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.0\n \"Reaction F\" => -0.422872\n \"Reaction T\" => -0.0013608\n \"Reaction G\" => 0.0\n \"Reaction B\" => -0.0\n \"Reaction J\" => 0.000486989\n \"Reaction K\" => 0.000486989\n \"Reaction P\" => -3.52357e-11\n \"Reaction R\" => -0.0\n \"Reaction U\" => 1.28067e-9\n \"Reaction Y\" => -0.0\n \"Reaction N\" => -0.0\n \"Reaction Z\" => -0.14112\n \"Reaction D\" => -0.0470399\n \"Reaction W\" => 0.000120942\n \"Reaction M\" => -0.0\n \"Reaction X\" => 0.000241884\n \"Reaction C\" => 1.10241e-10\n \"Reaction A\" => 1.38911e-11\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\" => 1.3891081053642988e-11\n \"Reaction C\" => -6.053049545510096e-11\n \"Reaction D\" => 1.2248420082882448e-11\n \"Reaction E\" => 1.3891081053642988e-11\n \"Reaction F\" => 1.3891081053642988e-11\n \"Reaction G\" => 1.3891081053642988e-11\n \"Reaction H\" => 1.6597522344141004e-11\n \"Reaction I\" => 1.3891081053642988e-11\n \"Reaction J\" => 1.7764900680738363e-11\n \"Reaction K\" => 1.3247412533182346e-11\n ⋮\n \"Reaction R\" => 1.3891081053642988e-11\n \"Reaction S\" => 5.991945842008422e-10\n \"Reaction T\" => 1.3686861323282557e-11\n \"Reaction U\" => 1.2424891707562945e-9\n \"Reaction V\" => 1.6597522344141004e-11\n \"Reaction W\" => 1.3595610516296605e-11\n \"Reaction X\" => 8.980140092800575e-12\n \"Reaction Y\" => 1.3891081053642988e-11\n \"Reaction Z\" => 1.3891081053642988e-11"},"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.1.0"},"kernelspec":{"name":"julia-1.1","display_name":"Julia 1.1.0","language":"julia"}},"nbformat":4}