Unverified Commit dc8b02ba authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

bulldoze FBA and FVA

parent 4fc2d47c
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
A variant of FBA that returns a vector of fluxes in the same order as reactions A variant of FBA that returns a vector of fluxes in the same order as reactions
of the model, if the solution is found. of the model, if the solution is found.
Arguments are passed to [`flux_balance_analysis`](@ref). Arguments are passed to [`flux_balance_analysis`](@ref).
""" """
function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},Nothing} function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},Nothing}
...@@ -14,16 +15,16 @@ function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},No ...@@ -14,16 +15,16 @@ function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},No
end end
""" """
flux_balance_analysis_dict(model::M, args...)::Union{Dict{String, Float64},Nothing} where {M <: MetabolicModel} flux_balance_analysis_dict(model::MetabolicModel, args...)::Union{Dict{String, Float64},Nothing}
A variant of FBA that returns a dictionary assigning fluxes to reactions, if A variant of FBA that returns a dictionary assigning fluxes to reactions, if
the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref). the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref).
""" """
function flux_balance_analysis_dict( function flux_balance_analysis_dict(
model::M, model::MetabolicModel,
args...; args...;
kwargs..., kwargs...,
)::Union{Dict{String,Float64},Nothing} where {M<:MetabolicModel} )::Union{Dict{String,Float64},Nothing}
v = flux_balance_analysis_vec(model, args...; kwargs...) v = flux_balance_analysis_vec(model, args...; kwargs...)
isnothing(v) && return nothing isnothing(v) && return nothing
Dict(zip(reactions(model), v)) Dict(zip(reactions(model), v))
...@@ -33,51 +34,45 @@ end ...@@ -33,51 +34,45 @@ end
flux_balance_analysis( flux_balance_analysis(
model::M, model::M,
optimizer; optimizer;
modifications = [(model, opt_model) -> nothing], modifications = [],
) where {M<:MetabolicModel} ) where {M<:MetabolicModel}
Run flux balance analysis (FBA) on the `model` optionally specifying Run flux balance analysis (FBA) on the `model` optionally specifying
`modifications` to the problem. `modifications` to the problem. Basically, FBA solves this optimization problem:
Effectively solves the optimization problem:
``` ```
max cᵀx max cᵀx
s.t. S x = b s.t. S x = b
xₗ ≤ x ≤ xᵤ xₗ ≤ x ≤ xᵤ
``` ```
Optionally, you may specify one or more "modifications" to be applied to the models. The `optimizer` must be set to a `JuMP`-compatible optimizer, such as
[`change_solver_attribute`](@ref),[`change_objective`](@ref), or `GLPK.Optimizer` or `Tulip.Optimizer`
[`change_sense`](@ref) for examples of modifications.
The `optimizer` must be set to perform the analysis, any JuMP solver will work. Optionally, you may specify one or more modifications to be applied to the
model before the analysis, such as
[`change_solver_attribute`](@ref),[`change_objective`](@ref), and
[`change_sense`](@ref).
Returns a solved JuMP model from [`optimize_model`](@ref). Returns an optimized `JuMP` model.
# Example # Example
``` ```
optimizer = GLPK.Optimizer
model = load_model(StandardModel, "e_coli_core.json") model = load_model(StandardModel, "e_coli_core.json")
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM") biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
solved_model = fba(model, optimizer; modifications=[change_objective(biomass)]) solution = flux_balance_analysis(model, GLPK.optimizer; modifications=[change_objective(biomass)])
``` ```
""" """
function flux_balance_analysis( function flux_balance_analysis(
model::M, model::M,
optimizer; optimizer;
modifications = [(model, opt_model) -> nothing], modifications = [],
) where {M<:MetabolicModel} ) where {M<:MetabolicModel}
opt_model = make_optimization_model(model, optimizer) opt_model = make_optimization_model(model, optimizer)
# support for multiple modification, fallback to single one for mod in modifications
if typeof(modifications) <: AbstractVector mod(model, opt_model)
for mod in modifications
mod(model, opt_model)
end
else
modifications(model, opt_model)
end end
COBREXA.JuMP.optimize!(opt_model) COBREXA.JuMP.optimize!(opt_model)
......
""" """
flux_variability_analysis( flux_variability_analysis(
model::LM, model::MetabolicModel,
reactions::Vector{Int}, reactions::Vector{Int},
optimizer, optimizer;
workers = [myid()]; modifications = [],
gamma::AbstractFloat = 1.0, workers = [myid()],
)::Matrix{Float64} where {LM<:MetabolicModel} bounds = z -> (z,z),
ret = objective_value,
# Flux variability analysis (FVA) )::Matrix{Float64}
FVA solves the pair of optimization problems in `model` for each flux xᵢ listed Flux variability analysis solves a pair of optimization problems in `model` for
in `reactions` each flux listed in `reactions`:
``` ```
min/max xᵢ min,max xᵢ
s.t. S x = b s.t. S x = b
xₗ ≤ x ≤ xᵤ xₗ ≤ x ≤ xᵤ
cᵀx ≥ γ Z₀ cᵀx ≥ bounds(Z₀)[1]
cᵀx ≤ bounds(Z₀)[2]
``` ```
where Z₀:= cᵀx₀ is the objective value of an optimal solution to the associated where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated
FBA problem. FBA problem (see [`flux_balance_analysis`](@ref)).
Internally uses the specified JuMP-compatible `optimizer`, and the work is The `bounds` is a user-supplied function that specifies the objective bounds
scheduled in parallel on `workers`. for the variability optimizations, by default it restricts the flux objective
value to the precise optimum reached in FBA. It can return `-Inf` and `Inf` in
Returns a matrix of minima and maxima of size (length(reactions),2). first and second pair to remove the limit. Use [`gamma_bound`](@ref) and
[`optimum_bound`](@ref) for simple bounds.
`optimizer` must be set to a `JuMP`-compatible optimizer. The computation of
the individual optimization problems is transparently distributed to `workers`
(see `Distributed.workers()`).
`ret` is a function used to extract results from optimized JuMP models of the
individual reactions. More detailed information can be extracted e.g. by
setting it to `m -> (JuMP.objective_value(m), JuMP.value.(m[:x]))`.
Returns a matrix of extracted `ret` values for minima and maxima, of total size
`length(reactions)`×2. The optimizer result status is not checked by default,
instead `ret` function can access the `JuMP.termination_status` of the model
and react accordingly, depending on user decision.
""" """
function flux_variability_analysis( function flux_variability_analysis(
model::LM, model::MetabolicModel,
reactions::Vector{Int}, reactions::Vector{Int},
optimizer, optimizer;
workers = [myid()]; modifications = [],
gamma::AbstractFloat = 1.0, workers = [myid()],
)::Matrix{Float64} where {LM<:MetabolicModel} bounds = z -> (z, z),
ret = objective_value,
)
if any(reactions .< 1) || any(reactions .> n_reactions(model)) if any(reactions .< 1) || any(reactions .> n_reactions(model))
throw(DomainError(reactions, "Index exceeds number of reactions.")) throw(DomainError(reactions, "Index exceeds number of reactions."))
end end
optimization_model = flux_balance_analysis(model, optimizer) Z = bounds(
Z0 = COBREXA.JuMP.objective_value(optimization_model) objective_value(
optimization_model = nothing # we won't need this one anymore, so free the memory flux_balance_analysis(model, optimizer; modifications = modifications),
),
)
# store a JuMP optimization model at all workers # store a JuMP optimization model at all workers
save_model = :( save_model = :(
begin begin
optmodel = $COBREXA.make_optimization_model($model, $optimizer) optmodel = $COBREXA.make_optimization_model($model, $optimizer)
$COBREXA._FVA_add_constraint( $COBREXA._FVA_add_constraint(optmodel, $(objective(model)), optmodel[:x], $Z)
optmodel,
$(objective(model)),
optmodel[:x],
$Z0,
$gamma,
)
optmodel optmodel
end end
) )
...@@ -60,7 +72,7 @@ function flux_variability_analysis( ...@@ -60,7 +72,7 @@ function flux_variability_analysis(
# schedule FVA parts parallely using pmap # schedule FVA parts parallely using pmap
fluxes = dpmap( fluxes = dpmap(
rid -> :($COBREXA._FVA_optimize_reaction(cobrexa_parfva_model, $rid)), rid -> :($COBREXA._FVA_optimize_reaction(cobrexa_parfva_model, $rid, $ret)),
CachingPool(workers), CachingPool(workers),
[-reactions reactions], [-reactions reactions],
) )
...@@ -73,125 +85,95 @@ end ...@@ -73,125 +85,95 @@ end
""" """
flux_variability_analysis( flux_variability_analysis(
model::LM, model::MetabolicModel,
optimizer; optimizer;
gamma::AbstractFloat = 1.0, kwargs...
) where {LM<:MetabolicModel} )
A simpler version of FVA that maximizes and minimizes all reactions in the model. A simpler version of [`flux_variability_analysis`](@ref) that maximizes and minimizes all reactions in the model. Arguments are forwarded.
""" """
function flux_variability_analysis( function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
model::LM,
optimizer;
gamma::AbstractFloat = 1.0,
) where {LM<:MetabolicModel}
n = n_reactions(model) n = n_reactions(model)
return flux_variability_analysis(model, collect(1:n), optimizer; gamma = gamma) return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
end
"""
_FVA_add_constraint(model, c, x, Z0, gamma)
Internal helper function for adding constraints to a model. Exists mainly
because for avoiding namespace problems on remote workers.
"""
function _FVA_add_constraint(model, c, x, Z0, gamma)
@constraint(model, c' * x gamma * Z0)
end end
""" """
_FVA_get_opt(model, rid) flux_variability_analysis_dict(
model::MetabolicModel,
optimizer;
kwargs...
)
Helper for creating the optimized model on a remote worker, for avoiding A variant of [`flux_variability_analysis`](@ref) that returns the individual
namespace problems. maximized and minimized fluxes of all reactions as two dictionaries (of
dictionaries). All keyword arguments except `ret` are passed through.
""" """
function _FVA_optimize_reaction(model, rid) function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
sense = rid > 0 ? MOI.MAX_SENSE : MOI.MIN_SENSE vs = flux_variability_analysis(
var = all_variables(model)[abs(rid)] model,
optimizer;
kwargs...,
ret = m -> (objective_value(m), JuMP.value.(m[:x])),
)
rxns = reactions(model)
@objective(model, sense, var) return (
optimize!(model) Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for (_, fluxes) in vs[:, 1]])),
return objective_value(model) Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for (_, fluxes) in vs[:, 2]])),
)
end end
""" """
fva(model::StandardModel, optimizer; optimum_bound=1.0-_constants.tolerance, modifications) gamma_bounds(gamma)
Run flux variability analysis (FVA) on the `model` (of type `StandardModel`).
Optionally specifying problem modifications like in [`flux_balance_analysis`](@ref).
This algorithm runs FBA on the model to determine the optimum of the objective.
This optimum then constrains subsequent problems, where `optimum_bound` can be used to
relax this constraint as a fraction of the FBA optimum, e.g.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
Note, this function only runs serially.
Consider using a different model type for parallel implementations.
Returns two dictionaries (`fva_max` and `fva_min`) that map each reaction `id` to dictionaries of the resultant flux distributions when that `id` is optimized.
See also: [`CoreModel`](@ref)
# Example A bounds-generating function for [`flux_variability_analysis`](@ref) that
limits the objective value to be at least `gamma*Z₀`, as usual in COBRA
packages. Use as the `bounds` argument:
``` ```
optimizer = Gurobi.Optimizer flux_variability_analysis(model, some_optimizer; bounds = gamma_bounds(0.9))
model = load_model(StandardModel, "iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
``` ```
""" """
function flux_variability_analysis( gamma_bounds(gamma) = z -> (gamma * z, Inf)
model::StandardModel,
optimizer;
optimum_bound = 1.0 - _constants.tolerance,
modifications = [(model, opt_model) -> nothing],
)
# Run FBA
opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || """
(return nothing, nothing) (tolerance) = z -> begin
fva_min = Dict{String,Union{Nothing,Dict{String,Float64}}}() A bounds-generating function for [`flux_variability_analysis`](@ref) that
fva_max = Dict{String,Union{Nothing,Dict{String,Float64}}}() limits the objective value to a small multiple of Z₀. Use as `bounds` argument,
similarly to [`gamma_bounds`](@ref).
"""
objective_bounds(tolerance) = z -> begin
vs = (z * tolerance, z / tolerance)
(minimum(vs), maximum(vs))
end
# Now do FVA """
v = opt_model[:x] _FVA_add_constraint(model, c, x, Z)
λ = COBREXA.JuMP.objective_value(opt_model) # objective value Internal helper function for adding constraints to a model. Exists mainly
λmin = min(optimum_bound * λ, λ * 1.0 / optimum_bound) because for avoiding namespace problems on remote workers.
λmax = max(optimum_bound * λ, λ * 1.0 / optimum_bound) """
function _FVA_add_constraint(model, c, x, Z)
Z[1] > -Inf && @constraint(model, c' * x >= Z[1])
Z[2] < Inf && @constraint(model, c' * x <= Z[2])
end
COBREXA.JuMP.@constraint( """
opt_model, _FVA_get_opt(model, rid)
λmin <= COBREXA.JuMP.objective_function(opt_model) <= λmax # in case there is a negative bound
)
for (i, r_id) in enumerate(reactions(model)) Internal helper for creating the optimized model on a remote worker, for
COBREXA.JuMP.@objective(opt_model, Max, v[i]) avoiding namespace problems.
COBREXA.JuMP.optimize!(opt_model) """
status = ( function _FVA_optimize_reaction(model, rid, ret)
COBREXA.JuMP.termination_status(opt_model) == MOI.OPTIMAL || sense = rid > 0 ? MOI.MAX_SENSE : MOI.MIN_SENSE
COBREXA.JuMP.termination_status(opt_model) == MOI.LOCALLY_SOLVED var = all_variables(model)[abs(rid)]
)
if status
fva_max[r_id] = Dict(zip(reactions(model), value.(opt_model[:x])))
else
@warn "Error maximizing index: $i with error $(termination_status(opt_model))"
fva_max[r_id] = nothing
end
@objective(opt_model, Min, v[i]) @objective(model, sense, var)
optimize!(opt_model) optimize!(model)
status = (
termination_status(opt_model) == MOI.OPTIMAL ||
termination_status(opt_model) == MOI.LOCALLY_SOLVED
)
if status
fva_min[r_id] = Dict(zip(reactions(model), value.(opt_model[:x])))
else
@warn "Error minimizing index: $i with error $(termination_status(opt_model))"
fva_min[r_id] = nothing
end
end
return fva_max, fva_min if termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
return ret(model)
else
return nothing
end
end end
...@@ -71,10 +71,10 @@ end ...@@ -71,10 +71,10 @@ end
sol_multi = flux_balance_analysis_dict( sol_multi = flux_balance_analysis_dict(
model, model,
Tulip.Optimizer; Tulip.Optimizer;
modifications = change_objective( modifications = [change_objective(
["BIOMASS_Ecoli_core_w_GAM", "PFL"]; ["BIOMASS_Ecoli_core_w_GAM", "PFL"];
weights = [biomass_frac, pfl_frac], weights = [biomass_frac, pfl_frac],
), )],
) )
@test isapprox( @test isapprox(
biomass_frac * sol_multi["BIOMASS_Ecoli_core_w_GAM"] + pfl_frac * sol_multi["PFL"], biomass_frac * sol_multi["BIOMASS_Ecoli_core_w_GAM"] + pfl_frac * sol_multi["PFL"],
......
...@@ -30,13 +30,13 @@ ...@@ -30,13 +30,13 @@
0.0 0.0 0.0 0.0
-1.0 -1.0 -1.0 -1.0
] ]
fluxes = flux_variability_analysis(cp, optimizer; gamma = 0.5) fluxes = flux_variability_analysis(cp, optimizer; bounds = gamma_bounds(0.5))
@test fluxes [ @test fluxes [
0.5 1.0 0.5 1.0
0.0 0.5 0.0 0.5
-1.0 -0.5 -1.0 -0.5
] ]
fluxes = flux_variability_analysis(cp, optimizer; gamma = 0.0) fluxes = flux_variability_analysis(cp, optimizer; bounds = _ -> (0, Inf))
@test fluxes [ @test fluxes [
0.0 1.0 0.0 1.0
0.0 1.0 0.0 1.0
...@@ -52,7 +52,7 @@ end ...@@ -52,7 +52,7 @@ end
cp = test_simpleLP() cp = test_simpleLP()
pids = addprocs(2, topology = :master_worker) pids = addprocs(2, topology = :master_worker)
@everywhere using COBREXA, GLPK @everywhere using COBREXA, GLPK
fluxes = flux_variability_analysis(cp, [1, 2], GLPK.Optimizer, pids) fluxes = flux_variability_analysis(cp, [1, 2], GLPK.Optimizer; workers=pids)
@test fluxes [ @test fluxes [
1.0 1.0 1.0 1.0
2.0 2.0 2.0 2.0
...@@ -68,10 +68,10 @@ end ...@@ -68,10 +68,10 @@ end
) )
model = load_model(StandardModel, model_path) model = load_model(StandardModel, model_path)
fva_max, fva_min = flux_variability_analysis( mins, maxs = flux_variability_analysis_dict(
model, model,
Tulip.Optimizer; Tulip.Optimizer;
optimum_bound = 0.99, bounds=objective_bounds(0.99),
modifications = [ modifications = [
change_solver_attribute("IPM_IterationsLimit", 500), change_solver_attribute("IPM_IterationsLimit", 500),
change_constraint("EX_glc__D_e", -10, -10), change_constraint("EX_glc__D_e", -10, -10),
...@@ -79,6 +79,6 @@ end ...@@ -79,6 +79,6 @@ end
], ],
) )
@test isapprox(fva_max["EX_ac_e"]["EX_ac_e"], 8.518549434876208, atol = TEST_TOLERANCE) @test isapprox(maxs["EX_ac_e"]["EX_ac_e"], 8.518549434876208, atol = TEST_TOLERANCE)
@test isapprox(fva_min["EX_ac_e"]["EX_ac_e"], 7.448388738973361, atol = TEST_TOLERANCE) @test isapprox(mins["EX_ac_e"]["EX_ac_e"], 7.448388738973361, atol = TEST_TOLERANCE)
end 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