Commit a24bab78 authored by St. Elmo's avatar St. Elmo
Browse files

added lots of docs, tests and code

parent 78d19a4d
......@@ -79,28 +79,24 @@ sol = map_fluxes(v, model)
More funcionality is described in the documention, e.g. model construction and analysis in pure Julia.
## Progress
- [ ] Read SBML models
- [ ] Write SBML models
- [ ] Single gene knockouts
- [ ] Thermodynamic FBA (and related functions)
- [x] Read JSON "Cobrapy" models
- [x] Read Matlab models
- [ ] Read SBML models
- [x] Write JSON models
- [x] Write Matlab models
- [ ] Write SBML models
- [x] FBA
- [X] pFBA
- [ ] FVA
- [x] pFBA
- [x] FVA
- [x] Implement sampling (hit and run)
- [x] Implement sampling (achr - kind of, fix constraint issue...)
- [ ] Single gene knockouts
- [ ] Double gene knockout
- [x] Equilibrator integration
- [x] Brenda integration (basic)
- [x] Reaction construction
- [x] Model construction
- [x] Model modifications
- [ ] Thermodynamic FBA (and related functions)
### Citations
1) Ebrahim, A., Lerman, J.A., Palsson, B.O. & Hyduke, D. R. (2013). COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Systems Biology, 7(74). https://doi.org/10.1186/1752-0509-7-74
......
......@@ -13,11 +13,9 @@ makedocs(
"Model Structure" => "model_structure.md",
"Model IO" => "io.md",
"Model Construction" => "model_construction.md",
"Optimization Based Analysis" => "basic_analysis.md",
"Optimization Based Analysis Tools" => "basic_analysis.md",
"Sampling Tools" => "sampling_tools.md",
"Equilibrator Interface" => "thermo_tools.md",
"Brenda Interface" => "brenda_tools.md",
"Thermodynamic Analysis" => "thermodynamics.md"
"External Tools" => "external_tools.md",
]
)
......
# Optimization Based Analysis
A selection of standard COBRA functions have been implemented to make basic model analysis more convenient.
Additionally, `CobraTools.jl` allows you to easily formulate your own optimization problems using the structure of a constraint based model.
This makes it easy to experiment with custom algorithms etc.
## Flux balance analysis (FBA)
Flux balance analysis solves the linear program,
......@@ -40,6 +42,18 @@ open("fluxes.json", "w") do io
JSON.print(io, sol)
end
```
## Solution inspection
Sometimes it is useful to investigate which reactions consume or produce a certain metabolite, or which exchange reactions are active.
This functionality is exposed via `metabolite_fluxes` and `exchange_reactions`.
```@docs
metabolite_fluxes
exchange_reactions
```
```@example fba
consuming, producing = metabolite_fluxes(sol, model)
consuming["atp_c"]
```
## Parsimonious FBA
Parsimonious FBA (pFBA) solves a two stage optimization problem. First, a classic FBA problem is solved to identify the unique maximum of the objective.
However, it should be noted that the fluxes from FBA are not unique (i.e. many fluxes may yield the objective optimum).
......@@ -47,15 +61,15 @@ To yield a unique set of fluxes, and remove internal futile cycles, a secondary
Suppose that FBA has found the optimum of ``v_\mu = \mu``, pFBA then solves,
```math
\begin{aligned}
& \underset{v}{\text{max}}
& \underset{v}{\text{min}}
& \sum_{i} {v_i^2} \\
& \text{s. t.}
& Sv = 0 \\
& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\
& v_\mu = \mu
& & v_\mu = \mu
\end{aligned}
```
again using any JuMP compatible solver(s). If multiple solvers are given, the first solver is used to solve the LP, and the second solver the QP, otherwise the same solver is used to solve both problems.
again using any JuMP compatible solver(s). In the `CobraTools.jl` implementation of pFBA, both the FBA and QP problem are solved internally in `pfba`, using similar input arguments as in `fba`. If multiple solvers are given, the first solver is used to solve the LP, and the second solver the QP, otherwise the same solver is used to solve both problems.
This is useful if the QP solver does not handle the LP problem well, as with OSQP.
An alternative, related formulation of this idea exists, called "CycleFreeFlux".
......@@ -64,7 +78,7 @@ really have much benefit beyond only solving two linear programs. See [Building
```@docs
pfba
```
Here, we use `Tulip.jl` followed by `OSQP.jl`, with the `pfba` function from `CobraTools.jl`. Note that `OSQP.jl` has iffy performance, and is only included here because it is open source. We strongly recommend that a commercial solver, e.g. `Gubobi.jl` be used to simplify your user experience.
Here, we use `Tulip.jl` followed by `OSQP.jl`, with the `pfba` function from `CobraTools.jl`. Note that `OSQP.jl` has iffy performance, and is only included here because it is open source. We recommend that a commercial solver, e.g. `Gubobi.jl`, be used to simplify your user experience.
```@example fba
using OSQP
......@@ -72,10 +86,43 @@ atts = Dict("eps_abs" => 5e-4,"eps_rel" => 5e-4, "max_iter" => 100_000, "verbose
sol = pfba(model, biomass, [Tulip.Optimizer, OSQP.Optimizer]; solver_attributes=Dict("opt1" => Dict{Any, Any}(), "opt2" => atts))
```
## Flux variability analysis (FVA)
Flux variability analysis can also be used to investigate the degeneracy associated with flux balance analysis derived solutions (see also [Sampling Tools](@ref)).
`CobraTools.jl` exposes `fva` that sequentially maximizes and minimizes each reaction in a model subject to the constraint that each optimization problem also satisfies an initial FBA type objective optimum, below denoted by ``v_{\mu}=\mu``,
```math
\begin{aligned}
& \underset{v}{\text{max or min}}
& v_i \\
& \text{s. t.}
& Sv = 0 \\
& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\
& & v_{\mu} = \mu \\
\end{aligned}
```
```@docs
fva
```
## Building your own optimization analysis script
CobraTools.jl makes it simple to access the
`CobraTools.jl` also makes it simple to construct customized optimization problems by making judicious use of [JuMP](https://jump.dev/).
Convenience functions make optimization problem construction, modification and data extraction from JuMP result objects easy.
```@docs
map_fluxes
```
\ No newline at end of file
get_core_model
build_cbm
set_bound
map_fluxes(::Array{Float64,1}, ::CobraTools.Model)
```
```@example fba
using CobraTools
using JuMP
using Tulip
model = read_model(model_location)
cbm, v, mb, ubs, lbs = build_cbm(model)
glucose_index = model[findfirst(model.reactions, "EX_glc__D_e")]
set_bound(glucose_index, ubs, lbs; ub=-12.0, lb=-12.0)
set_optimizer(cbm, Tulip.Optimizer)
@objective(cbm, Max, v[model[biomass]])
optimize!(cbm)
sol = map_fluxes(v, model)
```
# Brenda Interface
\ No newline at end of file
# External Tools
Numerous external databases exist that can add functionality to constraint based models.
Currently, `CobraTools.jl` tries to make it easier to incorporate two such databases in analyses.
The first is Brenda, a database that maps enzymes to kinetic data.
The second is Equilibrator, a database that allows one to calculate thermodynamic information, e.g. ``\Delta G``, from reactions.
`CobraTools.jl` includes lightweight interfaces to these sources.
## Brenda Interface
To use the Brenda interface functions, you will need to download the database as a txt file [available here](https://www.brenda-enzymes.org/download_brenda_without_registration.php) (~250 MB).
Once the database has been downloaded, it can be parsed by `parse_brenda`.
```@docs
parse_brenda
```
This function returns an array of `BrendaEntry` structs, which are composed of `EnzymeParams` for each field extracted from the Brenda database.
Currently, only the ID (=EC number), TN (=turn over number), KM (=Michaelis-Menten constant, ``K_M``), KI (=Inhibition term for Michaelis-Menten kintics), and KMM (=ratio of TN/KM) numbers for each enzyme class (ID, or EC number) are extracted. All the structs have pretty printing enabled.
```@docs
CobraTools.BrendaEntry
CobraTools.EnzymeParams
```
```@setup brenda
brenda_loc = model_location = joinpath("..","..", "test", "data", "small_brenda.txt")
```
```@example brenda
using CobraTools
brenda_data = parse_brenda(brenda_loc)
```
```@example brenda
brenda_data[1]
```
```@example brenda
brenda_data[1].TN
```
## Equilibrator Interface
The Equilibrator interface requires that the Equilibrator-API has been installed and can be accessed through Julia's PyCall package. Refer to the [Equilibrator-API website](https://gitlab.com/equilibrator/equilibrator-api) for installation instructions. Within Julia, if you can call `pyimport("equilibrator_api")` successfully, then you will be able to use the functions exposed here. To actually use the functions insert `using PyCall` in your main level script (before or after `using CobraTools`).
```@docs
map_gibbs_rxns
map_gibbs_external
map_gibbs_internal
```
# Sampling Tools
\ No newline at end of file
# Sampling Tools
It is well known that the FBA does not yield a unique solution, i.e. many flux distributions are capable of satisfying the system constraints as well as optimizing the imposed objective function.
Let the feasible space be defined by ``\mathcal{P} = \left\{ v : Sv = 0 \cap v_{\text{min}} \leq v \leq v_{\text{max}} \right\}``.
Sampling methods have been developed to uniformly sample from this feasible solution space.
`CobraTools.jl` implements both `hit_and_run` and `achr` to sample from ``\mathcal{P}``.
```@docs
hit_and_run
achr
```
```@setup sample
model_location = joinpath("..","..", "models", "e_coli_core.json")
```
```@example sample
using CobraTools
using JuMP
using Tulip
model = CobraTools.read_model(model_location)
optimizer = Tulip.Optimizer
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
cons = Dict("EX_glc__D_e" => (-12.0, -12.0))
sol = fba(model, biomass, optimizer, constraints=cons) # classic flux balance analysis
cons["BIOMASS_Ecoli_core_w_GAM"] = (sol["BIOMASS_Ecoli_core_w_GAM"], sol["BIOMASS_Ecoli_core_w_GAM"]*0.99)
samples = CobraTools.hit_and_run(100_000, model, optimizer; keepevery=10, samplesize=5000, constraints=cons)
```
\ No newline at end of file
# Equilibrator Interface
The standard (temperature 25 °C, pressue 1 bar, concentration 1M) Gibbs energies of biochemical reactions at various pH levels (ionic strength = 0.1M) mapped to the KEGG database are made available from [Equilibrator](http://equilibrator.weizmann.ac.il/download).
\ No newline at end of file
# Thermodynamic Analysis
\ No newline at end of file
# Thermodynamic Analysis
TBD.
\ No newline at end of file
......@@ -25,9 +25,9 @@ include(joinpath("io", "io_tools.jl"))
include(joinpath("construction", "construction_overloading.jl"))
include(joinpath("construction", "model_manipulations.jl"))
include(joinpath("optimization_analysis", "basic_analysis.jl"))
# include("sampling_tools.jl")
# include("brenda_tools.jl")
# include("equilibrator_tools.jl")
include(joinpath("sampling", "sampling_tools.jl"))
include(joinpath("external", "brenda_tools.jl"))
include(joinpath("external", "equilibrator_tools.jl"))
= Metabolite("∅") # for exchange reactions
......@@ -44,7 +44,13 @@ export
read_model, save_model,
# optimization_analysis
get_core_model, build_cbm, fba, map_fluxes, set_bound, pfba, atom_exchange, exchange_reactions, metabolite_fluxes
get_core_model, build_cbm, fba, map_fluxes, set_bound, pfba, atom_exchange, exchange_reactions, metabolite_fluxes, fva,
# sampling
hit_and_run, test_samples, achr,
# external
parse_brenda, map_gibbs_rxns, map_gibbs_external, map_gibbs_internal
# Initialization functions
include("init_functions.jl")
......
"""
EnzymeParams
This struct has `val`, `substrate`, `ph`, and `temperature` fields, extrated from the database.
"""
struct EnzymeParams
val :: Union{Nothing, Float64}
......@@ -17,6 +19,8 @@ end
"""
BrendaEntry
This struct has ID, TN, KI, KM, KMM fields, corresponding to the fields in the database.
"""
struct BrendaEntry
ID :: String # EC number
......@@ -70,10 +74,11 @@ function Base.show(io::IO, m::BrendaEntry)
end
"""
brenda_data = parse_brenda(brenda_txt_file_path)
parse_brenda(brenda_txt_file_path)
Parse the Brenda txt file line by line.
Download it from https://www.brenda-enzymes.org/download_brenda_without_registration.php.
Returns an array with Brenda
"""
function parse_brenda(brenda_loc)
brenda_data = Array{BrendaEntry, 1}()
......
......@@ -37,6 +37,11 @@ end
map_gibbs_rxns(rxns::Array{Reaction, 1}; dgtype="zero", ph=7.0, ionic_str="100 mM", usekegg=true)
Return a dict of rxn.id => ΔG of the specified dgtype.
Takes as inputs an array of reactions, `rxns`, and optional keyword arguments `dgtype`, which specifies which type of ΔG calculation should be returned.
Valid options for `dgtype` are "phys" and "prime" (anything else is "zero").
Ionic strength can be set through `ionic_str` which takes a string input, e.g. "150 mM".
Only BIGG and KEGG metabolite identifiers are supported, i.e. the reaction needs to have a KEGG or BIGG `id` listed in the annotation field in the reaction struct.
By default KEGG annotations are used to build the reaction strings that are fed to Equilibrator. Note that the first metabolite `id` is used.
"""
function map_gibbs_rxns(rxns::Array{Reaction, 1}; dgtype="zero", ph=7.0, ionic_str="100 mM", usekegg=true)
if usekegg
......@@ -115,33 +120,4 @@ function map_gibbs_internal(fluxres::Dict{String, Float64}, gibbs, biomassid="BI
end
end
return total_ΔG, missing_flux/(missing_flux+found_flux) # units J/gDW/h
end
"""
save_gibbs(path, gibbs)
Save gibbs dict. Saved as String => [mag, err]
"""
function save_gibbs(path, gibbs)
decomp = Dict{String, Array{Float64,1}}()
for (k, v) in gibbs
decomp[k] = [Measurements.value(v), Measurements.uncertainty(v)]
end
open(path, "w") do io
JSON.print(io, decomp)
end
end
"""
load_Gibbs(path)
Load Gibbs dict. Loads String => [mag, err]
"""
function load_gibbs(path)
decomp = JSON.parsefile(path)
gibbs = Dict{String, Measurement{Float64}}()
for (k, vs) in decomp
gibbs[k] = vs[1] ± vs[2]
end
return gibbs
end
\ No newline at end of file
......@@ -3,9 +3,9 @@
Return stoichiometrix matrix (S), mass balance right hand side (b), upper (ubs), and lower bounds (lbs) of constraint based `model`.
That is, S*v=b with lbs ≤ v ≤ ubs where v is the flux vector.
This is useful if you want to construct your own optimization problem.
This is useful if you want to construct your own optimization problem and just want the raw data.
Returns: S, b, ubs, lbs
Returns: S, b, ubs, lbs. All these data are arrays.
"""
function get_core_model(model::CobraTools.Model)
ubs = [rxn.ub for rxn in model.reactions]
......@@ -30,10 +30,10 @@ end
Initialize a constraint based `model` using `JuMP`.
Creates a model that satisfies the mass balance and flux constraints but no objective or optimizer is set.
Returns the `JuMP` model.
This is useful if you want to write your own optimization problem.
This is useful if you want to write your own optimization problem, but want `CobraTools.jl` to construct the basic optimization problem for you.
Returns: `cbmodel`, `v`, `mb`, `ubs`, and `lbs`, where `cbmodel` is the JuMP model, `v` are the fluxes, `mb` is S*v == 0, and lbs <= v <= ubs.
All these variables are JuMP objects (not arrays).
"""
function build_cbm(model::CobraTools.Model)
S, b, ubs, lbs = get_core_model(model) # Construct S, b, lbs, ubs from model
......@@ -47,9 +47,10 @@ function build_cbm(model::CobraTools.Model)
end
@doc raw"""
fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
Run flux balance analysis (FBA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction).
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
......@@ -66,15 +67,22 @@ biomass = findfirst(model.reactions, "BIOMASS\_Ec\_iJO1366\_WT\_53p95M")
sol = fba(model, biomass, optimizer; solver_attributes=atts)
"""
function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
cbm, _, _, _, _ = build_cbm(model) # get the base constraint based model
function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
cbm, _, _, ubcons, lbcons = build_cbm(model) # get the base constraint based model
set_optimizer(cbm, optimizer) # choose optimizer
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
set_optimizer_attribute(cbm, k, v)
for (k, val) in solver_attributes
set_optimizer_attribute(cbm, k, val)
end
end
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
set_bound(ind, ubcons, lbcons; ub=con[1], lb=con[2])
end
# ensure that an array of objective indices are fed in
if typeof(objective_rxns) == Reaction
objective_indices = [model[objective_rxns]]
......@@ -114,10 +122,11 @@ function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reac
end
@doc raw"""
pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
Run parsimonious flux balance analysis (pFBA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction) for the initial FBA problem.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.
When `optimizer` is an array of optimizers, e.g. `[opt1, opt2]`, then `opt1` is used to solve the FBA problem, and `opt2` is used to solve the QP problem.
This strategy is useful when the QP solver is not good at solving the LP problem.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, k, v)`.
......@@ -136,10 +145,10 @@ biomass = findfirst(model.reactions, "BIOMASS\_Ec\_iJO1366\_WT\_53p95M")
sol = pfba(model, biomass, optimizer; solver_attributes=atts)
"""
function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
## FBA ################################################
cbm, _, _, _, _ = build_cbm(model) # get the base constraint based model
cbm, _, _, ubcons, lbcons = build_cbm(model) # get the base constraint based model
if typeof(optimizer) <: AbstractArray # choose optimizer
set_optimizer(cbm, optimizer[1])
if !isempty(solver_attributes["opt1"]) # set other attributes
......@@ -156,6 +165,12 @@ function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Rea
end
end
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
set_bound(ind, ubcons, lbcons; ub=con[1], lb=con[2])
end
# ensure that an array of objective indices are fed in
if typeof(objective_rxns) == Reaction
objective_indices = [model[objective_rxns]]
......@@ -235,30 +250,28 @@ Map fluxes from an optimization problem (`v`) to rxns in a model.
`v` can be a JuMP object (fluxes) or an array of Float64 fluxes.
Assumes they are in order of `model.reactions`, which they should be since the optimization problem is constructed from the model.
"""
function map_fluxes(v::Array{VariableRef,1}, model::CobraTools.Model)
function map_fluxes(v::Array{Float64,1}, model::CobraTools.Model)
rxndict = Dict{String, Float64}()
for i in eachindex(model.reactions)
rxndict[model.reactions[i].id] = value(v[i])
rxndict[model.reactions[i].id] = v[i]
end
return rxndict
end
function map_fluxes(v::Array{Float64,1}, model::CobraTools.Model)
function map_fluxes(v::Array{VariableRef,1}, model::CobraTools.Model)
rxndict = Dict{String, Float64}()
for i in eachindex(model.reactions)
rxndict[model.reactions[i].id] = v[i]
rxndict[model.reactions[i].id] = value(v[i])
end
return rxndict
end
"""
setbound(index, ubconstaintref, lbconstaintref; ub=1000, lb=-1000)
set_bound(index, ubconstaintref, lbconstaintref; 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
The JuMP `set_normalized_rhs` function is a little confusing, so this function simplifies setting
constraints. Just supply the constraint `index` and the bound objects (`ubconstaintref`, `lbconstaintref`) and they will be set to `ub` and `lb`.
See also: [`set_normalized_rhs`](@ref)
"""
function set_bound(vind, ubs, lbs; ub=1000, lb=-1000)
if lb <= 0
......@@ -338,7 +351,7 @@ end
"""
metabolite_fluxes(fluxdict::Dict{String, Float64}, model::CobraTools.Model)
Return two dictionaries of metabolite `id`s mapped to reactions that consume or produce them.
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::CobraTools.Model)
S, _, _, _ = get_core_model(model)
......@@ -370,3 +383,107 @@ function metabolite_fluxes(fluxdict::Dict{String, Float64}, model::CobraTools.Mo
end
return consuming, producing
end
@doc raw"""
fva(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; optimum_bound=0.9999, weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
Run flux variability analysis (FVA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction).
It runs fba on the model once to determine the optimum of the objective.
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.
The model is then constrained to produce objective flux bounded by `optimum_bound` from below (set to slightly less than 1.0 for stability) and each flux in the model sequentially minimized and maximized.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns two dictionaries (fva_max and fva_min) that each reaction `id`s to dictionaries of the resultant flux distributions (if solved successfully) when that `id` is optimized.
# Example
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS\_Ec\_iJO1366\_WT\_53p95M")
fva\_max, fva\_min = fva(model, biomass, optimizer; solver_attributes=atts)
"""
function fva(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; optimum_bound=0.9999, weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
cbm, _, _, ubcons, lbcons = build_cbm(model) # get the base constraint based model
set_optimizer(cbm, optimizer) # choose optimizer
if !isempty(solver_attributes) # set other attributes
for (k, v) in solver_attributes
set_optimizer_attribute(cbm, k, v)
end
end
# set additional constraints
for (rxnid, con) in constraints
ind = model.reactions[findfirst(model.reactions, rxnid)]
set_bound(ind, ubcons, lbcons; ub=con[1], lb=con[2])
end
# ensure that an array of objective indices are fed in
if typeof(objective_rxns) == Reaction
objective_indices = [model[objective_rxns]]
else
objective_indices = [model[rxn] for rxn in objective_rxns]
end
if isempty(weights)
weights = ones(length(objective_indices))
end
opt_weights = zeros(length(model.reactions))
# update the objective function tracker
wcounter = 1
for i in eachindex(model.reactions)
if i in objective_indices
model.reactions[i].objective_coefficient = weights[wcounter]
opt_weights[i] = weights[wcounter]
wcounter += 1
else
model.reactions[i].objective_coefficient = 0.0
end
end
v = all_variables(cbm)
@objective(cbm, Max, sum(opt_weights[i]*v[i] for i in objective_indices))
optimize!(cbm)
status = (termination_status(cbm) == MOI.OPTIMAL || termination_status(cbm) == MOI.LOCALLY_SOLVED)
fva_min = Dict{String, Dict{String, Float64}}()
fva_max = Dict{String, Dict{String, Float64}}()
if !status
@warn "Error getting the initial optimum, aborting "
return fva_max, fva_min
end
λ = objective_value(cbm) # objective value
@constraint(cbm, optimum_bound*λ <= sum(opt_weights[i]*v[i] for i in objective_indices) <= λ) # for stability
for i in 1:length(v)
@objective(cbm, Max, v[i])
optimize!(cbm)
status = (termination_status(cbm) == MOI.OPTIMAL || termination_status(cbm) == MOI.LOCALLY_SOLVED)
if status
fva_max[model.reactions[i].id] = map_fluxes(v, model)
else
@warn "Error maximizing index: $i with error $(termination_status(cbm))"
end
@objective(cbm, Min, v[i])
optimize!(cbm)
status = (termination_status(cbm) == MOI.OPTIMAL || termination_status(cbm) == MOI.LOCALLY_SOLVED)
if status
fva_min[model.reactions[i].id] = map_fluxes(v, model)
else
@warn "Error minimizing index: $i with error $(termination_status(cbm))"
end
end
return fva_max, fva_min
end
\ No newline at end of file
"""
get_warmup_points(cbmodel, v, mb, ubs, lbs; randomobjective=false, numstop=1e10)
get_warmup_points(cbmodel, v, mb, ubs, lbs; random_objective=false, numstop=1e10)
Generate warmup points for all the reactions on the model that
are not fixed. Assumes you feed in a JuMP model that is already
......@@ -8,7 +8,7 @@ Note, extra constraints applied to ubs and lbs will have no effect.
numstop = 2*number of warmup points - to reduce the time this takes
"""
function get_warmup_points(cbmodel, v, ubs, lbs; randomobjective=false, numstop=1e10)
function get_warmup_points(cbmodel, v, ubs, lbs; random_objective=false, numstop=1e10)
# determine which rxns should be max/min-ized
fixed_rxns = Int64[]
for i in eachindex(v)
......@@ -29,7 +29,7 @@ function get_warmup_points(cbmodel, v, ubs, lbs; randomobjective=false, numstop=
for (i, ii) in zip(1:length(var_rxn_inds), 1:2:(2*length(var_rxn_inds)))
i > NN && break
if randomobjective
if random_objective
@objective(cbmodel, Max, sum(rand()*v[iii] for iii in var_rxn_inds))