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

Merge pull request #331 from LCSB-BioCore/mk-317

make atom_exchange and metabolite_fluxes a bit more generic
parents c60bda6d 7f22afea
Pipeline #43058 failed with stages
in 10 minutes and 2 seconds
......@@ -41,7 +41,7 @@ model = load_model(StandardModel, "e_coli_core.json") # we specifically want to
using Tulip
dict_sol = flux_balance_analysis_dict(
fluxes = flux_balance_analysis_dict(
model,
Tulip.Optimizer;
modifications = [
......@@ -59,18 +59,24 @@ dict_sol = flux_balance_analysis_dict(
# 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`.
# [`atom_exchange`](@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:
atom_exchange(dict_sol, model) # flux of individual atoms entering and leaving the system through boundary reactions (e.g. exchange reactions) based on flux_dict
fluxes_without_biomass = copy(fluxes);
delete!(fluxes_without_biomass, "BIOMASS_Ecoli_core_w_GAM");
atom_exchange(model, fluxes_without_biomass)
# ## Inspecting the flux solution: `metabolite_fluxes`
# Another useful flux result analysis function is `metabolite_fluxes`. This
# function keeps track of reactions consuming and producing each metabolite.
# Another useful flux result analysis function is [`metabolite_fluxes`](@ref).
# This function gives an overview of reactions consuming and producing each
# metabolite.
consuming, producing = metabolite_fluxes(dict_sol, model)
consuming, producing = metabolite_fluxes(model, fluxes)
consuming["atp_c"] # reactions consuming atp_c
consuming["atp_c"] # reactions consuming `atp_c`
# ## Internals of `StandardModel`
......@@ -170,8 +176,8 @@ check_duplicate_reaction(pgm_duplicate, model.reactions; only_metabolites = fals
# ## Checking the internals of `StandardModel`s: `is_mass_balanced`
# Finally, `is_mass_balanced` can be used to check if a reaction is mass
# Finally, [`is_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(pgm_duplicate, model) # extra_atoms shows which atoms are in excess/deficit
is_bal, extra_atoms = is_mass_balanced(model, pgm_duplicate) # extra_atoms shows which atoms are in excess/deficit
......@@ -50,14 +50,14 @@ function is_boundary(rxn::Reaction)::Bool
end
"""
is_mass_balanced(rxn::Reaction, model::StandardModel)
is_mass_balanced(model::StandardModel, rxn::Reaction)
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).
See also: [`get_atoms`](@ref), [`check_duplicate_reaction`](@ref)
"""
function is_mass_balanced(rxn::Reaction, model::StandardModel)
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)
......
......@@ -52,34 +52,13 @@ Shallow copy of a [`Gene`](@ref)
Base.copy(g::Gene) = Gene(g.id; name = g.name, notes = g.notes, annotations = g.annotations)
"""
atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)
atom_exchange(model::StandardModel, rxn_id::String)
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`.
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(flux_dict::Dict{String,Float64}, model::StandardModel)
atom_flux = Dict{String,Float64}()
for (rxn_id, flux) in flux_dict
if is_boundary(model.reactions[rxn_id])
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) + flux * stoich_rxn * stoich_molecule
end
end
end
end
return atom_flux
end
"""
atom_exchange(rxn_id::String, model::StandardModel)
Return a dictionary mapping the flux of atoms through a reaction in `model`.
"""
function atom_exchange(rxn_id::String, model::StandardModel)
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)
......@@ -90,38 +69,3 @@ function atom_exchange(rxn_id::String, model::StandardModel)
end
return atom_flux
end
"""
metabolite_fluxes(flux_dict::Dict{String, Float64}, model::StandardModel)
Return two dictionaries of metabolite `id`s mapped to reactions that consume or
produce them given the flux distribution supplied in `fluxdict`.
"""
function metabolite_fluxes(flux_dict::Dict{String,Float64}, model::StandardModel)
S = stoichiometry(model)
rxnids = reactions(model)
metids = metabolites(model)
producing = Dict{String,Dict{String,Float64}}()
consuming = Dict{String,Dict{String,Float64}}()
for (row, metid) in enumerate(metids)
for (col, rxnid) in enumerate(rxnids)
mf = flux_dict[rxnid] * S[row, col]
# ignore zero flux
if mf < -_constants.tolerance # consuming rxn
if haskey(consuming, metid)
consuming[metid][rxnid] = mf
else
consuming[metid] = Dict(rxnid => mf)
end
elseif mf > _constants.tolerance
if haskey(producing, metid)
producing[metid][rxnid] = mf
else
producing[metid] = Dict(rxnid => mf)
end
end
end
end
return consuming, producing
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
"""
metabolite_fluxes(model::MetabolicModel, flux_dict::Dict{String, Float64})
Return two dictionaries of metabolite `id`s mapped to reactions that consume or
produce them, given the flux distribution supplied in `flux_dict`.
"""
function metabolite_fluxes(model::MetabolicModel, flux_dict::Dict{String,Float64})
S = stoichiometry(model)
rids = reactions(model)
mids = metabolites(model)
producing = Dict{String,Dict{String,Float64}}()
consuming = Dict{String,Dict{String,Float64}}()
for (row, mid) in enumerate(mids)
for (col, rid) in enumerate(rids)
mf = flux_dict[rid] * S[row, col]
# ignore zero flux
if mf < -_constants.tolerance # consuming rxn
if haskey(consuming, mid)
consuming[mid][rid] = mf
else
consuming[mid] = Dict(rid => mf)
end
elseif mf > _constants.tolerance
if haskey(producing, mid)
producing[mid][rid] = mf
else
producing[mid] = Dict(rid => mf)
end
end
end
end
return consuming, producing
end
......@@ -2,31 +2,22 @@
model = load_model(StandardModel, model_paths["e_coli_core.json"])
# FBA
glc = model.reactions["EX_glc__D_e"]
optimizer = Tulip.Optimizer # quiet by default
sol = flux_balance_analysis_dict(
fluxes = flux_balance_analysis_dict(
model,
optimizer;
Tulip.Optimizer;
modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
)
# atom tracker
atom_fluxes = atom_exchange(sol, model)
@test isapprox(atom_fluxes["C"], 37.19016648975907; atol = TEST_TOLERANCE)
@test atom_exchange("FBA", model)["C"] == 0.0
# single-reaction atom exchanges
@test atom_exchange(model, "FBA")["C"] == 0.0
@test isapprox(
atom_exchange("BIOMASS_Ecoli_core_w_GAM", model)["C"],
atom_exchange(model, "BIOMASS_Ecoli_core_w_GAM")["C"],
-42.5555;
atol = TEST_TOLERANCE,
)
# metabolite trackers
consuming, producing = metabolite_fluxes(sol, model)
@test isapprox(consuming["atp_c"]["PFK"], -7.47738; atol = TEST_TOLERANCE)
@test isapprox(producing["atp_c"]["PYK"], 1.75818; atol = TEST_TOLERANCE)
# set bounds
cbm = make_optimization_model(model, optimizer)
# bounds setting
cbm = make_optimization_model(model, Tulip.Optimizer)
ubs = cbm[:ubs]
lbs = cbm[:lbs]
glucose_index = first(indexin(["EX_glc__D_e"], reactions(model)))
......
@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
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")],
)
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)
end
Supports Markdown
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