Unverified Commit f5b7f450 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil Committed by GitHub
Browse files

Merge pull request #414 from LCSB-BioCore/mo-update-docstrings-atom-exchange

馃獡 some code and streamline utilities for reaction/metabolite balancing
parents c73f5155 3df683c3
Pipeline #45075 passed with stages
in 33 minutes and 46 seconds
......@@ -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`.
"""
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)
......@@ -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)
......@@ -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
"""
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
......@@ -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
......@@ -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"])
......
......@@ -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]
......
@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
@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
@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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment