Unverified Commit 0128a96b authored by St. Elmo's avatar St. Elmo
Browse files

implement reviews

parent 6ea1e8f6
......@@ -12,6 +12,7 @@ DistributedData = "f6a0035f-c5ac-4ad0-b410-ad102ced35df"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
......
......@@ -25,6 +25,9 @@ import SBML # conflict with Reaction struct name
include("banner.jl")
_print_banner()
# Constants
const DEFAULT_FVA_TOL = 1e-6 # for numerical stability when doing FVA
# autoloading
const _inc(path...) = include(joinpath(path...))
const _inc_all(dir) = _inc.(joinpath.(dir, filter(fn -> endswith(fn, ".jl"), readdir(dir))))
......
......@@ -82,7 +82,8 @@ Arguments are passed to [`flux_balance_analysis`](@ref).
function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},Nothing}
optmodel = flux_balance_analysis(args...; kwargs...)
JuMP.termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
COBREXA.JuMP.termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing
value.(optmodel[:x])
end
......@@ -154,6 +155,6 @@ function flux_balance_analysis(
modifications(model, opt_model)
end
JuMP.optimize!(opt_model)
COBREXA.JuMP.optimize!(opt_model)
return opt_model
end
......@@ -116,17 +116,20 @@ function _FVA_optimize_reaction(model, rid)
end
"""
fva(model::StandardModel, optimizer; optimum_bound=0.9999, modifications)
fva(model::StandardModel, optimizer; optimum_bound=1.0-DEFAULT_FVA_TOL, modifications)
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.
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: [`LinearModel`](@ref)
# Example
```
optimizer = Gurobi.Optimizer
......@@ -138,50 +141,58 @@ fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
function flux_variability_analysis(
model::StandardModel,
optimizer;
optimum_bound = 0.9999,
optimum_bound = 1.0-DEFAULT_FVA_TOL,
modifications = [(model, opt_model) -> nothing],
)
# Run FBA
cbm = flux_balance_analysis(model, optimizer; modifications=modifications)
opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || (return nothing, nothing)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
(return nothing, nothing)
fva_min = Dict{String,Union{Nothing, Dict{String,Float64}}}()
fva_max = Dict{String,Union{Nothing, Dict{String,Float64}}}()
fva_min = Dict{String,Union{Nothing,Dict{String,Float64}}}()
fva_max = Dict{String,Union{Nothing,Dict{String,Float64}}}()
# Now do FVA
v = cbm[:x]
v = opt_model[:x]
λ = COBREXA.JuMP.objective_value(opt_model) # objective value
λmin = min(optimum_bound*λ, λ * 1.0 / optimum_bound)
λmax = max(optimum_bound*λ, λ * 1.0 / optimum_bound)
λ = JuMP.objective_value(cbm) # objective value
JuMP.@constraint(
cbm,
optimum_bound * λ <= JuMP.objective_function(cbm) <= λ * 1.0/optimum_bound # in case there is a negative bound and the problem has min sense
) # for stability
COBREXA.JuMP.@constraint(
opt_model,
λmin <=
COBREXA.JuMP.objective_function(opt_model) <=
λmax # in case there is a negative bound
)
for i = 1:length(v)
JuMP.@objective(cbm, Max, v[i])
JuMP.optimize!(cbm)
COBREXA.JuMP.@objective(opt_model, Max, v[i])
COBREXA.JuMP.optimize!(opt_model)
status = (
JuMP.termination_status(cbm) == MOI.OPTIMAL ||
JuMP.termination_status(cbm) == MOI.LOCALLY_SOLVED
COBREXA.JuMP.termination_status(opt_model) == MOI.OPTIMAL ||
COBREXA.JuMP.termination_status(opt_model) == MOI.LOCALLY_SOLVED
)
if status
fva_max[model.reactions[i].id] = Dict(zip(reactions(model), value.(cbm[:x])))
fva_max[model.reactions[i].id] =
Dict(zip(reactions(model), value.(opt_model[:x]) ))
else
@warn "Error maximizing index: $i with error $(termination_status(cbm))"
@warn "Error maximizing index: $i with error $(termination_status(opt_model))"
fva_max[model.reactions[i].id] = nothing
end
@objective(cbm, Min, v[i])
optimize!(cbm)
@objective(opt_model, Min, v[i])
optimize!(opt_model)
status = (
termination_status(cbm) == MOI.OPTIMAL ||
termination_status(cbm) == MOI.LOCALLY_SOLVED
termination_status(opt_model) == MOI.OPTIMAL ||
termination_status(opt_model) == MOI.LOCALLY_SOLVED
)
if status
fva_min[model.reactions[i].id] = Dict(zip(reactions(model), value.(cbm[:x])))
fva_min[model.reactions[i].id] =
Dict(zip(reactions(model), value.(opt_model[:x]) ))
else
@warn "Error minimizing index: $i with error $(termination_status(cbm))"
@warn "Error minimizing index: $i with error $(termination_status(opt_model))"
fva_min[model.reactions[i].id] = nothing
end
end
......
......@@ -20,41 +20,52 @@ biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = pfba(model, biomass, optimizer; solver_attributes=atts)
```
"""
function parsimonious_flux_balance_analysis(model::MetabolicModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
function parsimonious_flux_balance_analysis(
model::MetabolicModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
)
# Run FBA
cbm = flux_balance_analysis(model, optimizer)#; modifications=modifications)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing # FBA failed
opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing # FBA failed
# Run pFBA
λ = objective_value(cbm)
old_objective = JuMP.objective_function(cbm)
λ = objective_value(opt_model)
old_objective = COBREXA.JuMP.objective_function(opt_model)
qp_solver(model, cbm) # change the solver if specified, otherwise default argument does nothing
if typeof(qp_solver_attributes) <: Vector # many modifications
qp_solver(model, opt_model) # change the solver if specified, otherwise default argument does nothing
if typeof(qp_solver_attributes) <: AbstractVector # many modifications
for mod in qp_solver_attributes
mod(model, cbm)
mod(model, opt_model)
end
else # single modification
qp_solver_attributes(model, cbm)
qp_solver_attributes(model, opt_model)
end
v = cbm[:x] # fluxes
@objective(cbm, Min, sum(dot(v, v)))
v = opt_model[:x] # fluxes
@objective(opt_model, Min, sum(dot(v, v)))
for lbconval in [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99] # sequentially relax bound for stability
λmin = min(lbconval*λ, λ * 1.0 / lbconval)
λmax = max(lbconval*λ, λ * 1.0 / lbconval)
@constraint(
cbm,
opt_model,
pfbacon,
lbconval * λ <= old_objective <= λ * 1.0/lbconval # in case of negative objective and minimized fba
λmin <= old_objective <= λmax # in case of negative constraints
)
optimize!(cbm)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] && break
COBREXA.JuMP.delete(cbm, pfbacon)
optimize!(opt_model)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] &&
break
COBREXA.JuMP.delete(opt_model, pfbacon)
end
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing # pFBA failed
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing # pFBA failed
return cbm
return opt_model
end
"""
......@@ -64,11 +75,24 @@ Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a vector of fluxes in the same order as the reactions in `model`.
Calls [`parsimonious_flux_balance_analysis_dict`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_vec(model::StandardModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
cbm = parsimonious_flux_balance_analysis(model, optimizer; modifications=modifications, qp_solver=qp_solver, qp_solver_attributes=qp_solver_attributes)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return value.(cbm[:x])
function parsimonious_flux_balance_analysis_vec(
model::StandardModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
)
opt_model = parsimonious_flux_balance_analysis(
model,
optimizer;
modifications = modifications,
qp_solver = qp_solver,
qp_solver_attributes = qp_solver_attributes,
)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing
return value.(opt_model[:x])
end
"""
......@@ -78,9 +102,22 @@ Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping reaction `id`s to fluxes.
Calls [`parsimonious_flux_balance_analysis_dict`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_dict(model::StandardModel, optimizer; modifications = [(model, opt_model) -> nothing], qp_solver = [(model, opt_model) -> nothing], qp_solver_attributes=[(model, opt_model) -> nothing])
cbm = parsimonious_flux_balance_analysis(model, optimizer; modifications=modifications, qp_solver=qp_solver, qp_solver_attributes=qp_solver_attributes)
JuMP.termination_status(cbm) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return Dict(zip(reactions(model), value.(cbm[:x])))
function parsimonious_flux_balance_analysis_dict(
model::StandardModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
)
opt_model = parsimonious_flux_balance_analysis(
model,
optimizer;
modifications = modifications,
qp_solver = qp_solver,
qp_solver_attributes = qp_solver_attributes,
)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing
return Dict(zip(reactions(model), value.(opt_model[:x])))
end
"""
atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)
Return a dictionary mapping the flux of atoms across the boundary of the model given `flux_dict` of reactions in `model`.
Here `flux_dict` is a mapping of reaction `id`s to fluxes, e.g. from FBA.
"""
function atom_exchange(flux_dict::Dict{String,Float64}, model::StandardModel)
atom_flux = Dict{String,Float64}()
for (rxnid, flux) in flux_dict
if startswith(rxnid, "EX_") || startswith(rxnid, "DM_") # exchange, demand reaction
for (met, stoich) in findfirst(model.reactions, rxnid).metabolites
adict = get_atoms(met)
for (atom, stoich) in adict
atom_flux[atom] = get(atom_flux, atom, 0.0) + flux * stoich
end
end
end
end
return atom_flux
end
"""
get_exchanges(rxndict::Dict{String, Float64}; top_n=8, ignorebound=1000.0, verbose=true)
Display the top_n producing and consuming exchange fluxes.
Set top_n to a large number to get all the consuming/producing fluxes.
Ignores infinite (problem upper/lower bound) fluxes (set with ignorebound).
When `verbose` is false, the output is not printed out.
Return these reactions in two dictionaries: `consuming`, `producing`
"""
function exchange_reactions(
rxndict::Dict{String,Float64};
top_n = 8,
ignorebound = 1000.0,
verbose = true,
)
fluxes = Float64[]
rxns = String[]
for (k, v) in rxndict
if startswith(k, "EX_") && abs(v) < ignorebound
push!(rxns, k)
push!(fluxes, v)
end
end
inds_prod = sortperm(fluxes, rev = true)
inds_cons = sortperm(fluxes)
consuming = Dict{String,Float64}()
producing = Dict{String,Float64}()
verbose && println("Consuming fluxes:")
for i = 1:min(top_n, length(rxndict))
if rxndict[rxns[inds_cons[i]]] < -eps()
verbose && println(
rxns[inds_cons[i]],
" = ",
round(rxndict[rxns[inds_cons[i]]], digits = 4),
)
consuming[rxns[inds_cons[i]]] = rxndict[rxns[inds_cons[i]]]
else
continue
end
end
verbose && println("Producing fluxes:")
for i = 1:min(top_n, length(rxndict))
if rxndict[rxns[inds_prod[i]]] > eps()
verbose && println(
rxns[inds_prod[i]],
" = ",
round(rxndict[rxns[inds_prod[i]]], digits = 4),
)
producing[rxns[inds_prod[i]]] = rxndict[rxns[inds_prod[i]]]
else
continue
end
end
return consuming, producing
end
"""
metabolite_fluxes(fluxdict::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(fluxdict::Dict{String,Float64}, model::StandardModel)
S = Array(stoichiometry(model))
met_flux = Dict{String,Float64}()
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 = fluxdict[rxnid] * S[row, col]
# ignore zero flux
if mf < -eps() # consuming rxn
if haskey(consuming, metid)
consuming[metid][rxnid] = mf
else
consuming[metid] = Dict(rxnid => mf)
end
elseif mf > eps()
if haskey(producing, metid)
producing[metid][rxnid] = mf
else
producing[metid] = Dict(rxnid => mf)
end
end
end
end
return consuming, producing
end
"""
set_bound(index, optimization_model; ub=1000, lb=-1000)
Helper function to set the bounds of variables.
The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
constraints. Just supply the constraint `index` and the model and they will be set to `ub` and `lb`.
"""
function set_bound(vind, opt_model; ub = 1000, lb = -1000)
if lb <= 0
set_normalized_rhs(opt_model[:lbs][vind], abs(lb))
else
set_normalized_rhs(opt_model[:lbs][vind], -abs(lb))
end
set_normalized_rhs(opt_model[:ubs][vind], ub)
end
"""
modify_constraint(reaction::Reaction, lb, ub)
Modify constraints of model reaction.
"""
function modify_constraint(reaction::Reaction, lb, ub)
(model, opt_model) -> begin
ind = model.reactions[reaction]
set_bound(ind, opt_model, lb = lb, ub = ub)
end
end
"""
modify_solver_attribute(option_key, option_val)
Modify a solver attribute. These attributes are solver specific,
refer the either JuMP or the solver you are using's documentation.
"""
function modify_solver_attribute(option_key, option_val)
(model, opt_model) -> begin
JuMP.set_optimizer_attribute(opt_model, option_key, option_val)
end
end
"""
modify_sense(objective_sense)
Modify the objective sense.
Possible arguments are `MOI.MAX_SENSE` and `MOI.MIN_SENSE`.
Note, [`modify_objective`](@ref) sets the sense of the objective,
so it doesn't make sense to use this function AND [`modify_objective`](@ref) simultaneously.
"""
function modify_sense(objective_sense)
(model, opt_model) -> begin
JuMP.set_objective_sense(opt_model, objective_sense)
end
end
"""
modify_objective(objective_functions::Union{Reaction, Array{Reaction, 1}}; weights=[], sense=MOI.MAX_SENSE)
Callback function to modify the objective function used in a constraint based analysis function.
`objective_functions` can be a single reaction or an array of reactions (of type `Reaction`).
Optionally specify their `weights` and the sense of the new objective (`MOI.MAX_SENSE`, `MOI.MIN_SENSE`).
Note, this function sets the sense of the objective.
"""
function modify_objective(
objective_functions::Union{Reaction,Array{Reaction,1}};
weights = [],
sense = MOI.MAX_SENSE,
)
(model, opt_model) -> begin
# Construct objective_indices array
if typeof(objective_functions) == Reaction
objective_indices = [model[objective_functions]]
else
objective_indices = [model[rxn] for rxn in objective_functions]
end
# Initialize weights
opt_weights = zeros(length(model.reactions))
isempty(weights) && (weights = ones(length(objective_indices))) # equal weights
wcounter = 1
for i in eachindex(model.reactions)
if i in objective_indices
opt_weights[i] = weights[wcounter]
wcounter += 1
end
end
v = opt_model[:x]
@objective(opt_model, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
end
end
"""
modify_solver(optimizer)
Modify the solver used to solve the model.
Typically the solver is specified as a required argument.
This function is useful if the problem has multiple subparts that require different solvers.
See
"""
function modify_solver(optimizer)
(model, opt_model) -> begin
JuMP.set_optimizer(opt_model, optimizer)
end
end
......@@ -10,7 +10,7 @@ function make_optimization_model(
m, n = size(stoichiometry(model))
xl, xu = bounds(model)
optimization_model = JuMP.Model(optimizer)
optimization_model = COBREXA.JuMP.Model(optimizer)
@variable(optimization_model, x[i = 1:n])
@objective(optimization_model, sense, objective(model)' * x)
@constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
......@@ -29,6 +29,6 @@ function optimize_model(
sense = MOI.MIN_SENSE,
) where {LM<:MetabolicModel}
optimization_model = make_optimization_model(model, optimizer; sense = sense)
JuMP.optimize!(optimization_model)
COBREXA.JuMP.optimize!(optimization_model)
return optimization_model
end
......@@ -18,3 +18,291 @@ Base.isequal(model1::CoupledLinearModel, model2::CoupledLinearModel) =
Base.copy(model::CoupledLinearModel) =
CoupledLinearModel(model.lm, model.C, model.cl, model.cu)
"""
atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)
Return a dictionary mapping the flux of atoms across the boundary of the model given `flux_dict` of reactions in `model`.
Here `flux_dict` is a mapping of reaction `id`s to fluxes, e.g. from FBA.
"""
function atom_exchange(flux_dict::Dict{String,Float64}, model::StandardModel)
atom_flux = Dict{String,Float64}()
for (rxnid, flux) in flux_dict
if startswith(rxnid, "EX_") || startswith(rxnid, "DM_") # exchange, demand reaction
for (met, stoich) in findfirst(model.reactions, rxnid).metabolites
adict = get_atoms(met)
for (atom, stoich) in adict
atom_flux[atom] = get(atom_flux, atom, 0.0) + flux * stoich
end
end
end
end
return atom_flux
end
"""
get_exchanges(rxndict::Dict{String, Float64}; top_n=8, ignorebound=1000.0, verbose=true)
Display the top_n producing and consuming exchange fluxes.
Set top_n to a large number to get all the consuming/producing fluxes.
Ignores infinite (problem upper/lower bound) fluxes (set with ignorebound).
When `verbose` is false, the output is not printed out.
Return these reactions in two dictionaries: `consuming`, `producing`
"""
function exchange_reactions(
rxndict::Dict{String,Float64};
top_n = 8,
ignorebound = 1000.0,
verbose = true,
)
fluxes = Float64[]
rxns = String[]
for (k, v) in rxndict
if startswith(k, "EX_") && abs(v) < ignorebound
push!(rxns, k)
push!(fluxes, v)
end
end
inds_prod = sortperm(fluxes, rev = true)
inds_cons = sortperm(fluxes)
consuming = Dict{String,Float64}()
producing = Dict{String,Float64}()
verbose && println("Consuming fluxes:")
for i = 1:min(top_n, length(rxndict))
if rxndict[rxns[inds_cons[i]]] < -eps()
verbose && println(
rxns[inds_cons[i]],
" = ",
round(rxndict[rxns[inds_cons[i]]], digits = 4),
)
consuming[rxns[inds_cons[i]]] = rxndict[rxns[inds_cons[i]]]
else
continue
end
end
verbose && println("Producing fluxes:")
for i = 1:min(top_n, length(rxndict))
if rxndict[rxns[inds_prod[i]]] > eps()
verbose && println(
rxns[inds_prod[i]],
" = ",
round(rxndict[rxns[inds_prod[i]]], digits = 4),
)
producing[rxns[inds_prod[i]]] = rxndict[rxns[inds_prod[i]]]
else
continue
end
end
return consuming, producing
end
"""
metabolite_fluxes(fluxdict::Dict{String, Float64}, model::StandardModel)