Re-organized reaction utilities and added tests. Also reorganized fluxes and renamed atom_fluxes

fix tests

fix test

docs/src/notebooks/3_basic_stdmodel_usage.jl

@@ -55,18 +55,18 @@ fluxes = flux_balance_analysis_dict(
 # However, deeper inspection of flux results is possible when using
 # `StandardModel`.
 
-# ## Inspecting the flux solution: `atom_exchange`
+# ## Inspecting the flux solution: `atom_fluxes`
 
 # It is sometimes interesting to keep track of the atoms entering and leaving
 # the system through boundary reactions. This can be inspected by calling
-# [`atom_exchange`](@ref). That gives you the flux of individual atoms getting
+# [`atom_fluxes`](@ref). That gives you the flux of individual atoms getting
 # consumed and produced by all reactions, based on `fluxes`. We erase the
 # reaction that consumes the atoms for creating biomass, to see how much mass
 # the "rest" of the reaction produces for it:
 
 fluxes_without_biomass = copy(fluxes);
 delete!(fluxes_without_biomass, "BIOMASS_Ecoli_core_w_GAM");
-atom_exchange(model, fluxes_without_biomass)
+atom_fluxes(model, fluxes_without_biomass)
 
 # ## Inspecting the flux solution: `metabolite_fluxes`
 
@@ -174,10 +174,16 @@ pgm_duplicate
 #
 check_duplicate_reaction(pgm_duplicate, model.reactions; only_metabolites = false) # can also just check if only the metabolites are the same but different stoichiometry is used
 
-# ## Checking the internals of `StandardModel`s: `is_mass_balanced`
+# ## Checking the internals of `StandardModel`s: `reaction_mass_balanced`
 
-# Finally, [`is_mass_balanced`](@ref) can be used to check if a reaction is mass
+# Finally, [`reaction_mass_balanced`](@ref) can be used to check if a reaction is mass
 # balanced based on the formulas of the reaction equation.
 
-pgm_duplicate.metabolites = Dict{String,Float64}("3pg_c" => 1, "2pg_c" => -1, "h2o_c" => 1) # not mass balanced now
-is_bal, extra_atoms = is_mass_balanced(model, pgm_duplicate) # extra_atoms shows which atoms are in excess/deficit
+rxn_dict = Dict{String,Float64}("3pg_c" => 1, "2pg_c" => -1, "h2o_c" => 1)
+reaction_mass_balanced(model, rxn_dict)
+
+# Now to determine which atoms are unbalanced, you can use `reaction_atom_balance`
+reaction_atom_balance(model, rxn_dict)
+
+# Note, since `pgm_duplicate` is not in the model, we cannot use the other variants of this
+# function because they find the reaction equation stored inside the `model`.

src/base/utils/Metabolite.jl

-"""
-    get_atoms(met::Metabolite)::MetaboliteFormula
-
-Simple wrapper for getting the atom dictionary count out of a `Metabolite`.
-
-See also: [`metabolite_formula`](@ref)
-"""
-get_atoms(met::Metabolite)::MetaboliteFormula = _maybemap(_parse_formula, met.formula)

src/base/utils/Reaction.jl

@@ -7,7 +7,7 @@ Otherwise, compares metabolite `id`s and the absolute value of their stoichiomet
 If `rxn` has the same reaction equation as another reaction in `rxns`, the return the `id`.
 Otherwise return `nothing`.
 
-See also: [`is_mass_balanced`](@ref)
+See also: [`reaction_mass_balanced`](@ref)
 """
 function check_duplicate_reaction(
     crxn::Reaction,
@@ -39,54 +39,64 @@ function check_duplicate_reaction(
 end
 
 """
-    is_boundary(rxn::Reaction)
+    is_boundary(reaction::Reaction)
 
 Return true if reaction is a boundary reaction, otherwise return false.
-Checks if boundary by inspecting number of metabolites in reaction equation.
+Checks if on boundary by inspecting the number of metabolites in reaction equation.
 Boundary reactions have only one metabolite, e.g. an exchange reaction, or a sink/demand reaction.
 """
 function is_boundary(rxn::Reaction)::Bool
     length(keys(rxn.metabolites)) == 1 ? true : false
 end
 
+is_boundary(model::StandardModel, rxn_id::String) = is_boundary(model.reactions[rxn_id])
+
+is_boundary(model::StandardModel, rxn::Reaction) = is_boundary(rxn) # for consistency with functions below
+
 """
-    is_mass_balanced(model::StandardModel, rxn::Reaction)
+    reaction_atom_balance(model::StandardModel, rxn)
 
-Checks if `rxn` is atom balanced. Returns a boolean for whether the reaction is balanced,
-and the associated balance of atoms for convenience (useful if not balanced).
+Returns a dictionary mapping the stoichiometry of atoms through a single reaction. Uses the
+metabolite information in `model` to determine the mass balance. Accepts a reaction
+dictionary, a reaction string id or a `Reaction` as an argument for `rxn`.
 
-See also: [`get_atoms`](@ref), [`check_duplicate_reaction`](@ref)
+See also: [`reaction_mass_balanced`](@ref)
 """
-function is_mass_balanced(model::StandardModel, rxn::Reaction)
-    atom_balances = Dict{String,Float64}() # float here because stoichiometry is not Int
-    for (met, stoich) in rxn.metabolites
-        atoms = metabolite_formula(model, met)
-        isnothing(atoms) &&
+function reaction_atom_balance(model::StandardModel, reaction_dict::Dict{String,Float64})
+    atom_balances = Dict{String,Float64}()
+    for (met, stoich_rxn) in reaction_dict
+        adict = metabolite_formula(model, met)
+        isnothing(adict) &&
             throw(ErrorException("Metabolite $met does not have a formula assigned to it."))
-        for (k, v) in atoms
-            atom_balances[k] = get(atom_balances, k, 0) + v * stoich
+        for (atom, stoich_molecule) in adict
+            atom_balances[atom] =
+                get(atom_balances, atom, 0.0) + stoich_rxn * stoich_molecule
         end
     end
+    return atom_balances
+end
 
-    return all(sum(values(atom_balances)) == 0), atom_balances
+function reaction_atom_balance(model::StandardModel, rxn_id::String)
+    reaction_atom_balance(model, model.reactions[rxn_id].metabolites)
 end
 
+reaction_atom_balance(model::StandardModel, rxn::Reaction) =
+    reaction_atom_balance(model, rxn.id)
+
 """
-    stoichiometry_string(rxn_dict)
+    reaction_mass_balanced(model::StandardModel, rxn)
 
-Return the reaction equation as a string.
+Checks if `rxn` is atom balanced. Returns a boolean for whether the reaction is balanced,
+and the associated balance of atoms for convenience (useful if not balanced). Calls
+`reaction_atom_balance` internally.
 
-# Example
-```
-julia> req = Dict("coa_c" => -1, "for_c" => 1, "accoa_c" => 1, "pyr_c" => -1);
-julia> stoichiometry_string(req)
-"coa_c + pyr_c = for_c + accoa_c"
-```
+See also: [`check_duplicate_reaction`](@ref), [`reaction_atom_balance`](@ref)
 """
-function stoichiometry_string(req)
-    replace_one(n) = abs(n) == 1 ? "" : string(abs(n))
-    substrates =
-        join(strip.(["$(replace_one(n)) $met" for (met, n) in req if n < 0]), " + ")
-    products = join(strip.(["$(replace_one(n)) $met" for (met, n) in req if n >= 0]), " + ")
-    return substrates * " = " * products
-end
+reaction_mass_balanced(model::StandardModel, rxn_id::String) =
+    all(values(reaction_atom_balance(model, rxn_id)) .== 0)
+
+reaction_mass_balanced(model::StandardModel, rxn::Reaction) =
+    reaction_mass_balanced(model, rxn.id)
+
+reaction_mass_balanced(model::StandardModel, reaction_dict::Dict{String,Float64}) =
+    all(values(reaction_atom_balance(model, reaction_dict)) .== 0)

src/base/utils/StandardModel.jl

@@ -48,22 +48,3 @@ Base.copy(m::Metabolite) = Metabolite(
 Shallow copy of a [`Gene`](@ref)
 """
 Base.copy(g::Gene) = Gene(g.id; notes = g.notes, annotations = g.annotations)
-
-"""
-    atom_exchange(model::StandardModel, rxn_id::String)
-
-Specialized version of [`atom_exchange`](@ref) for [`StandardModel`](@ref) that
-quickly returns a dictionary mapping the flux-relative stoichiometry of atoms
-through a single reaction in `model`.
-"""
-function atom_exchange(model::StandardModel, rxn_id::String)
-    atom_flux = Dict{String,Float64}()
-    for (met, stoich_rxn) in model.reactions[rxn_id].metabolites
-        adict = metabolite_formula(model, met)
-        isnothing(adict) && continue
-        for (atom, stoich_molecule) in adict
-            atom_flux[atom] = get(atom_flux, atom, 0.0) + stoich_rxn * stoich_molecule
-        end
-    end
-    return atom_flux
-end

src/base/utils/atom_exchange.jl

-"""
-    atom_exchange(model::MetabolicModel, flux_dict::Dict{String, Float64})
-
-Return a dictionary mapping the flux of atoms across the boundary of the model
-given `flux_dict` (the solution of a constraint based analysis) of reactions in `model`.
-"""
-function atom_exchange(model::MetabolicModel, flux_dict::Dict{String,Float64})
-    S = stoichiometry(model)
-    mids = metabolites(model)
-    rids = reactions(model)
-
-    atom_flux = Dict{String,Float64}()
-    for (ridx, rid) in enumerate(rids)
-        haskey(flux_dict, rid) || continue
-        rflux = flux_dict[rid]
-        for (midx, mstoi) in zip(findnz(S[:, ridx])...)
-            atoms = metabolite_formula(model, mids[midx])
-            isnothing(atoms) && continue
-            for (atom, abundance) in atoms
-                atom_flux[atom] = get(atom_flux, atom, 0.0) + rflux * mstoi * abundance
-            end
-        end
-    end
-    return atom_flux
-end

src/base/utils/fluxes.jl

@@ -32,3 +32,39 @@ function metabolite_fluxes(model::MetabolicModel, flux_dict::Dict{String,Float64
     end
     return consuming, producing
 end
+
+"""
+    atom_fluxes(model::MetabolicModel, reaction_fluxes::Dict{String, Float64})
+
+Return a dictionary mapping the flux of atoms across a flux solution given by
+`reaction_fluxes` using the reactions in `model` to determine the appropriate stoichiometry.
+
+Note, this function ignores metabolites with no formula assigned to them, no error message
+will be displayed.
+
+Note, if a model is mass balanced there should be not net flux of any atom. By removing
+reactions from the flux_solution you can investigate how that impacts the mass balances.
+
+# Example
+```
+# Find flux of Carbon through all metabolic reactions except the biomass reaction
+delete!(fluxes, "BIOMASS_Ecoli_core_w_GAM")
+atom_fluxes(model, fluxes)["C"]
+```
+"""
+function atom_fluxes(model::MetabolicModel, reaction_fluxes::Dict{String,Float64})
+    rids = reactions(model)
+    atom_flux = Dict{String,Float64}()
+    for (ridx, rid) in enumerate(rids)
+        haskey(reaction_fluxes, rid) || continue
+        rflux = reaction_fluxes[rid]
+        for (mid, mstoi) in reaction_stoichiometry(model, rid)
+            atoms = metabolite_formula(model, mid)
+            isnothing(atoms) && continue # missing formulas are ignored
+            for (atom, abundance) in atoms
+                atom_flux[atom] = get(atom_flux, atom, 0.0) + rflux * mstoi * abundance
+            end
+        end
+    end
+    return atom_flux
+end

test/base/types/Metabolite.jl

@@ -22,9 +22,6 @@
     m3.formula = "X"
     m3.annotations = Dict("sboterm" => ["sbo"], "kegg.compound" => ["ad2s", "asds"])
 
-    ats = get_atoms(m1)
-    @test ats["C"] == 6 && ats["N"] == 1
-
     m4 = Metabolite("met4")
     m4.formula = "X"
     m4.annotations = Dict("sboterm" => ["sbo"], "kegg.compound" => ["adxxx2s", "asdxxxs"])

test/base/utils/StandardModel.jl

@@ -8,14 +8,6 @@
         modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
     )
 
-    # single-reaction atom exchanges
-    @test atom_exchange(model, "FBA")["C"] == 0.0
-    @test isapprox(
-        atom_exchange(model, "BIOMASS_Ecoli_core_w_GAM")["C"],
-        -42.5555;
-        atol = TEST_TOLERANCE,
-    )
-
     # bounds setting
     cbm = make_optimization_model(model, Tulip.Optimizer)
     ubs = cbm[:ubs]

test/base/utils/atom_exchange.jl

-@testset "atom exchanges" begin
-    model = load_model(StandardModel, model_paths["e_coli_core.json"])
-
-    fluxes = flux_balance_analysis_dict(
-        model,
-        Tulip.Optimizer;
-        modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
-    )
-
-    # remove the biomass production from the fluxes, so that there's some atom
-    # disbalance that can be measured
-    delete!(fluxes, "BIOMASS_Ecoli_core_w_GAM")
-
-    # atom tracker
-    atom_fluxes = atom_exchange(model, fluxes)
-    @test isapprox(atom_fluxes["C"], 37.190166489763214; atol = TEST_TOLERANCE)
-    @test isapprox(atom_fluxes["O"], 41.663071522672226; atol = TEST_TOLERANCE)
-    @test isapprox(atom_fluxes["N"], 4.765319167566247; atol = TEST_TOLERANCE)
-end

test/base/utils/fluxes.jl

-@testset "flux utilities" begin
+@testset "Flux utilities" begin
     model = load_model(StandardModel, model_paths["e_coli_core.json"])
 
     fluxes = flux_balance_analysis_dict(
@@ -10,4 +10,14 @@
     consuming, producing = metabolite_fluxes(model, fluxes)
     @test isapprox(consuming["atp_c"]["PFK"], -7.47738; atol = TEST_TOLERANCE)
     @test isapprox(producing["atp_c"]["PYK"], 1.75818; atol = TEST_TOLERANCE)
+
+    # remove the biomass production from the fluxes, so that there's some atom
+    # disbalance that can be measured
+    delete!(fluxes, "BIOMASS_Ecoli_core_w_GAM")
+
+    # atom tracker
+    atom_fluxes_out = atom_fluxes(model, fluxes)
+    @test isapprox(atom_fluxes_out["C"], 37.190166489763214; atol = TEST_TOLERANCE)
+    @test isapprox(atom_fluxes_out["O"], 41.663071522672226; atol = TEST_TOLERANCE)
+    @test isapprox(atom_fluxes_out["N"], 4.765319167566247; atol = TEST_TOLERANCE)
 end

test/base/utils/reaction.jl

 @testset "Reaction utilities" begin
-    # test if reaction equation can be built back into a sensible reaction string
-    req = Dict("coa_c" => -1, "for_c" => 1, "accoa_c" => 1, "pyr_c" => -1)
-    rstr_out = stoichiometry_string(req)
-    @test occursin("coa_c", split(rstr_out, " = ")[1])
-    @test occursin("for", split(rstr_out, " = ")[2])
+    model = load_model(StandardModel, model_paths["e_coli_core.json"])
+
+    # FBA
+    fluxes = flux_balance_analysis_dict(
+        model,
+        Tulip.Optimizer;
+        modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
+    )
+
+    # test if reaction is balanced
+    @test reaction_mass_balanced(model, "PFL")
+    @test !reaction_mass_balanced(model, "BIOMASS_Ecoli_core_w_GAM")
+    @test reaction_mass_balanced(model, model.reactions["PFL"])
+    @test !reaction_mass_balanced(model, model.reactions["BIOMASS_Ecoli_core_w_GAM"])
+    @test reaction_mass_balanced(model, Dict("h_c" => -1.0, "h_e" => 1.0))
+    @test !reaction_mass_balanced(model, Dict("h_c" => -1.0, "h2o_c" => 1.0))
+
+    # test if a reaction is a boundary reaction
+    @test !is_boundary(model.reactions["PFL"])
+    @test is_boundary(model.reactions["EX_glc__D_e"])
+    @test !is_boundary(model, "PFL")
+    @test is_boundary(model, "EX_glc__D_e")
+    @test !is_boundary(model, model.reactions["PFL"])
+    @test is_boundary(model, model.reactions["EX_glc__D_e"])
+
+    # single-reaction atom balance
+    @test reaction_atom_balance(model, "FBA")["C"] == 0.0
+    @test isapprox(
+        reaction_atom_balance(model, "BIOMASS_Ecoli_core_w_GAM")["C"],
+        -42.5555;
+        atol = TEST_TOLERANCE,
+    )
+    @test reaction_atom_balance(model, model.reactions["FBA"])["C"] == 0.0
+    @test isapprox(
+        reaction_atom_balance(model, model.reactions["BIOMASS_Ecoli_core_w_GAM"])["C"],
+        -42.5555;
+        atol = TEST_TOLERANCE,
+    )
+    @test reaction_atom_balance(model, Dict("h_c" => -1.0, "h2o_c" => 1.0))["H"] == 1.0
 end