Commit 78d19a4d authored by St. Elmo's avatar St. Elmo
Browse files

added documentation for fba and pfba

parent 02284cca
# CobraTools.jl
*CobraTools is a Julia package for constraint based reconstruction and analysis of metabolic models.*
*CobraTools.jl is a Julia package for constraint based reconstruction and analysis of metabolic models.*
[docs-img]:https://img.shields.io/badge/docs-latest-blue.svg
[docs-url]: https://stelmo.github.io/CobraTools.jl/dev
......@@ -31,7 +31,7 @@ More importantly, it also exposes the user to the core structures used in COBRA,
To install this package: `] add https://github.com/stelmo/CobraTools.jl`. See the documentation for more information.
## Quick Example
Let's use `CobraTools` to perform classic flux balance analysis on an *E. coli* genome-scale metabolic model.
Let's use `CobraTools.jl` to perform classic flux balance analysis on an *E. coli* genome-scale metabolic model.
```julia
using CobraTools
using JuMP
......@@ -88,7 +88,6 @@ More funcionality is described in the documention, e.g. model construction and a
- [ ] Write SBML models
- [x] FBA
- [X] pFBA
- [ ] MOMA
- [ ] FVA
- [x] Implement sampling (hit and run)
- [x] Implement sampling (achr - kind of, fix constraint issue...)
......
# Optimization Based Analysis
## Convenience functions
A selection of standard COBRA functions have been implemented to make basic model analysis more convenient.
### FBA
## Flux balance analysis (FBA)
Flux balance analysis solves the linear program,
```math
\begin{aligned}
& \underset{v}{\text{max}}
& \sum_{i}^I {w_i}{v_i} \\
& \text{s. t.}
& Sv = 0 \\
& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\
\end{aligned}
```
using any `JuMP` compatible solver. Typically ``I`` is a singleton set that only includes the index of the biomass objective function, ``v_\mu``, and weight, ``w_\mu=1``.
```@docs
fba
```
Here, we use `Tulip.jl`, a pure Julia interior point linear program solver, with the `fba` function from `CobraTools.jl`.
```@setup fba
model_location = joinpath("..","..", "models", "e_coli_core.json")
```
```@example fba
using CobraTools
using JuMP
using Tulip
model = read_model(model_location)
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
optimizer = Tulip.Optimizer
sol = fba(model, biomass, optimizer)
```
Note that the result of the optimization problem, `sol` maps fluxes to reaction `id`s in the model, to simplify down stream analysis.
Since the fluxes are mapped to a dictionary, it makes it simple to export them as a JSON file for visualization in, e.g. [Escher](https://escher.github.io/#/).
```@example fba
using JSON
### pFBA
open("fluxes.json", "w") do io
JSON.print(io, sol)
end
```
## 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).
To yield a unique set of fluxes, and remove internal futile cycles, a secondary quadratic problem is imposed *ad hoc*.
Suppose that FBA has found the optimum of ``v_\mu = \mu``, pFBA then solves,
```math
\begin{aligned}
& \underset{v}{\text{max}}
& \sum_{i} {v_i^2} \\
& \text{s. t.}
& Sv = 0 \\
& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\
& 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.
This is useful if the QP solver does not handle the LP problem well, as with OSQP.
### Escher integration
An alternative, related formulation of this idea exists, called "CycleFreeFlux".
This replaces the quadratic formulation with an L1 (taxi cab) norm. While this should also remove futile cycles, it doesn't have the same uniqueness qualities and doesn't
really have much benefit beyond only solving two linear programs. See [Building your own optimization analysis script](@ref) for ideas about how to implement this yourself if you really need it.
```@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.
```@example fba
using OSQP
atts = Dict("eps_abs" => 5e-4,"eps_rel" => 5e-4, "max_iter" => 100_000, "verbose"=>false) # set solver attributes for QSQP
sol = pfba(model, biomass, [Tulip.Optimizer, OSQP.Optimizer]; solver_attributes=Dict("opt1" => Dict{Any, Any}(), "opt2" => atts))
```
## Flux variability analysis (FVA)
## Building your own optimization analysis script
CobraTools.jl makes it simple to access the
\ No newline at end of file
CobraTools.jl makes it simple to access the
```@docs
map_fluxes
```
\ No newline at end of file
# CobraTools.jl
*CobraTools is a Julia package for constraint based reconstruction and analysis of metabolic models.*
*CobraTools.jl is a Julia package for constraint based reconstruction and analysis of metabolic models.*
## Contents
```@contents
......
# Model IO
## Reading constraint based models
Currently, SBML, JSON and Matlab formatted models can be imported.
Currently, JSON and Matlab formatted models can be imported.
```@docs
read_model(file_location::String)
......
# Model Construction
## Defining genes
Genes are represented by the `Gene` type in `CobraTools`, see [Model Structure](@ref) for details.
Genes are represented by the `Gene` type in `CobraTools.jl`, see [Model Structure](@ref) for details.
`Gene`s can be constructed using either an empty constructor, or a constructor taking only
the string `id` of the gene.
```@docs
......@@ -21,7 +21,7 @@ findfirst(::Array{Gene, 1}, ::String)
getindex(::Array{Gene, 1}, ::Gene)
```
## Defining metabolites
Metabolites are represented by the `Metabolite` type in `CobraTools`, see [Model Structure](@ref) for details.
Metabolites are represented by the `Metabolite` type in `CobraTools.jl`, see [Model Structure](@ref) for details.
The simplest way to define a new metabolite is by using the empty constructor `Metabolite()`.
Alternatively, `Metabolite(id::String)` can be used to assign only the `id` field of the `Metabolite`.
```@docs
......@@ -60,7 +60,7 @@ findfirst(mets::Array{Metabolite, 1}, metid::String)
getindex(mets::Array{Metabolite, 1}, met::Metabolite)
```
## Defining reactions
Reactions are represented by the `Reaction` type in `CobraTools`, see [Model Structure](@ref) for details.
Reactions are represented by the `Reaction` type in `CobraTools.jl`, see [Model Structure](@ref) for details.
The simplest way to define a new reaction is by using the empty constructor `Reaction()`.
All the other fields still need to be assigned.
```@docs
......
......@@ -5,7 +5,7 @@ Each model is a struct of the type `CobraTools.Model`, which is composed of a mo
Model
```
The fields of `Reaction`, `Metabolite`, `Gene` types are shown below.
When reading, writing, building or analysing models, these fields are what is used by `CobraTools`.
When reading, writing, building or analysing models, these fields are what is used by `CobraTools.jl`.
```@docs
Reaction
```
......
......@@ -46,20 +46,24 @@ function build_cbm(model::CobraTools.Model)
return cbmodel, v, mb, ubs, lbs
end
"""
@doc raw"""
fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
Run flux balance analysis (FBA) on the `model` with `objective_rxn(s)` and optionally specifying their `weights` (empty `weights` mean equal weighting per reaction).
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, k, v)`.
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 a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.
# Example
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
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}())
......@@ -78,16 +82,17 @@ function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reac
objective_indices = [model[rxn] for rxn in objective_rxns]
end
# ensure corrects weights are given to objective
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
......@@ -95,7 +100,7 @@ function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reac
end
v = all_variables(cbm)
@objective(cbm, Max, sum(v[i] for i in objective_indices))
@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)
......@@ -108,12 +113,12 @@ function fba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reac
end
end
"""
@doc raw"""
pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}, optimizer; weights=Float64[], solver_attributes=Dict{Any, Any}())
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.
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.
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)`.
If more than one solver is specified in `optimizer`, then `solver_attributes` must be a dictionary of dictionaries with keys "opt1" and "opt2", e.g. Dict("opt1" => Dict{Any, Any}(),"opt2" => Dict{Any, Any}()).
......@@ -122,9 +127,13 @@ Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully,
# Example
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
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}())
......@@ -154,16 +163,17 @@ function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Rea
objective_indices = [model[rxn] for rxn in objective_rxns]
end
# ensure corrects weights are given to objective
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
......@@ -171,7 +181,7 @@ function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Rea
end
v = all_variables(cbm)
@objective(cbm, Max, sum(v[i] for i in objective_indices))
@objective(cbm, Max, sum(opt_weights[i]*v[i] for i in objective_indices))
optimize!(cbm)
fba_status = (termination_status(cbm) == MOI.OPTIMAL || termination_status(cbm) == MOI.LOCALLY_SOLVED)
......@@ -188,7 +198,7 @@ function pfba(model::CobraTools.Model, objective_rxns::Union{Reaction, Array{Rea
end
end
@constraint(cbm, pfbacon, 0.999999*λ <= sum(v[i] for i in objective_indices) <= λ) # constrain model - 0.9999 should be close enough?
@constraint(cbm, pfbacon, 0.999999*λ <= sum(opt_weights[i]*v[i] for i in objective_indices) <= λ) # constrain model - 0.9999 should be close enough?
@objective(cbm, Min, sum(dot(v, v)))
optimize!(cbm)
......
Markdown is supported
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