Commit 35a73fab authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

make atom_exchange and metabolite_fluxes a bit more generic

Closes #317
parent 131adbf3
......@@ -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,21 @@ 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 entering
# and leaving the system through boundary reactions (e.g. exchange reactions),
# based on `fluxes`:
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
atom_exchange(model, fluxes)
# ## 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`
......
......@@ -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::Dlicict{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)
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) + flux * 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,30 +2,21 @@
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
# bounds setting
cbm = make_optimization_model(model, optimizer)
ubs = cbm[:ubs]
lbs = cbm[:lbs]
......
@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")],
)
# atom tracker
atom_fluxes = atom_exchange(fluxes, model)
@test isapprox(atom_fluxes["C"], 37.19016648975907; atol = TEST_TOLERANCE)
@test isapprox(atom_fluxes["N"], 37.19016648975907; 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