Using a custom model data structure

This notebooks shows how to utilize the generic accessors and modification functions in COBREXA.jl to run the analysis on any custom model type. We will create a simple dictionary-style structure that describes the model, allow COBREXA to run a FVA on it, and create a simple reaction-removing modification.

First, let's define a very simple stoichiometry-only structure for the model:

using COBREXA

mutable struct MyReaction
    max_rate::Float64 # maximum absolute conversion rate
    stoi::Dict{String,Float64} # stoichimetry of the reaction

    MyReaction() = new(0.0, Dict{String,Float64}())
end

mutable struct MyModel <: MetabolicModel
    optimization_target::String # the "objective" reaction name
    reactions::Dict{String,MyReaction} # dictionary of reactions

    MyModel() = new("", Dict{String,MyReaction}())
    MyModel(o, r) = new(o, r)
end

With this, we can start defining the accessors:

COBREXA.n_reactions(m::MyModel) = length(m.reactions)
COBREXA.reactions(m::MyModel) = sort(collect(keys(m.reactions)))

Metabolites are defined only very implicitly, so let's just make a function that collects all names. n_metabolites can be left at the default definition that just measures the output of metabolites.

function COBREXA.metabolites(m::MyModel)
    mets = Set{String}()
    for (_, r) in m.reactions
        for (m, _) in r.stoi
            push!(mets, m)
        end
    end
    return sort(collect(mets))
end

Now, the extraction of the linear model. Remember the order of element in the vectors must match the order in the output of reactions and metabolites.

using SparseArrays

function COBREXA.bounds(m::MyModel)
    max_rates = [m.reactions[r].max_rate for r in reactions(m)]
    (sparse(-max_rates), sparse(max_rates))
end

function COBREXA.objective(m::MyModel)
    if m.optimization_target in keys(m.reactions)
        c = spzeros(n_reactions(m))
        c[first(indexin([m.optimization_target], reactions(m)))] = 1.0
        c
    else
        throw(
            DomainError(
                m.optimization_target,
                "The target reaction for flux optimization not found",
            ),
        )
    end
end

function COBREXA.stoichiometry(m::MyModel)
    sparse([
        get(m.reactions[rxn].stoi, met, 0.0) for met in metabolites(m), rxn in reactions(m)
    ])
end

Now the model is complete! We can generate a random one and see how it performs

import Random
Random.seed!(123)

rxn_names = ["Reaction $i" for i = 'A':'Z'];
metabolite_names = ["Metabolite $i" for i = 1:20];

m = MyModel();
for i in rxn_names
    m.reactions[i] = MyReaction()
end

for i = 1:50
    rxn = rand(rxn_names)
    met = rand(metabolite_names)
    m.reactions[rxn].stoi[met] = rand([-3, -2, -1, 1, 2, 3])
    m.reactions[rxn].max_rate = rand()
end

Let's see what the model looks like now:

Metabolic model of type Main.ex-8_custom_model.MyModel
sparse([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)
Number of reactions: 26
Number of metabolites: 19

We can run most of the standard function on the model data right away:

using Tulip
m.optimization_target = "Reaction A"
flux_balance_analysis_dict(m, Tulip.Optimizer)

Dict{String, Float64} with 26 entries:
  "Reaction E" => -0.0
  "Reaction F" => -0.422872
  "Reaction T" => -0.0013608
  "Reaction G" => 0.0
  "Reaction B" => -0.0
  "Reaction J" => 0.000486989
  "Reaction K" => 0.000486989
  "Reaction P" => -3.52357e-11
  "Reaction R" => -0.0
  "Reaction U" => 1.28067e-9
  "Reaction Y" => -0.0
  "Reaction N" => -0.0
  "Reaction Z" => -0.14112
  "Reaction D" => -0.0470399
  "Reaction W" => 0.000120942
  "Reaction M" => -0.0
  "Reaction X" => 0.000241884
  "Reaction C" => 1.10241e-10
  "Reaction A" => 1.38911e-11
  ⋮            => ⋮

To be able to use the model conveniently in functions such as screen, you usually want to be able to easily specify the modifications. In this example, we enable use of with_removed_reactions by overloading the internal remove_reactions for this specific model type:

We need to make an as-shallow-as-possible copy of the model that allows us to remove the reactions without breaking the original model.

function COBREXA.remove_reactions(m::MyModel, rxns::AbstractVector{String})
    m = MyModel(m.optimization_target, copy(m.reactions))
    delete!.(Ref(m.reactions), rxns)
    return m
end

The screening is ready now!

reactions_to_remove = ("Reaction $i" for i = 'B':'Z')

reactions_to_remove .=> screen_variants(
    m,
    [[with_removed_reactions([r])] for r in reactions_to_remove],
    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["Reaction A"],
)

25-element Vector{Pair{String, Float64}}:
 "Reaction B" => 1.3891081053642988e-11
 "Reaction C" => -6.053049545510096e-11
 "Reaction D" => 1.2248420082882448e-11
 "Reaction E" => 1.3891081053642988e-11
 "Reaction F" => 1.3891081053642988e-11
 "Reaction G" => 1.3891081053642988e-11
 "Reaction H" => 1.6597522344141004e-11
 "Reaction I" => 1.3891081053642988e-11
 "Reaction J" => 1.7764900680738363e-11
 "Reaction K" => 1.3247412533182346e-11
              ⋮
 "Reaction R" => 1.3891081053642988e-11
 "Reaction S" => 5.991945842008422e-10
 "Reaction T" => 1.3686861323282557e-11
 "Reaction U" => 1.2424891707562945e-9
 "Reaction V" => 1.6597522344141004e-11
 "Reaction W" => 1.3595610516296605e-11
 "Reaction X" => 8.980140092800575e-12
 "Reaction Y" => 1.3891081053642988e-11
 "Reaction Z" => 1.3891081053642988e-11

This page was generated using Literate.jl.