Skip to content

Commit

Permalink
Merge pull request #60 from LCSB-BioCore/mva-loopless-fba
Browse files Browse the repository at this point in the history
Added loopless fba
  • Loading branch information
laurentheirendt committed Apr 12, 2021
2 parents 469ec56 + dbdfa9d commit 3db6eb6
Show file tree
Hide file tree
Showing 9 changed files with 275 additions and 44 deletions.
8 changes: 6 additions & 2 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -23,5 +23,9 @@ docs/site/
# environment.
Manifest.toml

# Ignore temporary script for testing functions
temp.jl
# Ignore temporary files for testing functions
temp.*

# Ignore jupyter notebook stuff
.ipynb_checkpoints

2 changes: 1 addition & 1 deletion docs/src/basic_analysis.md
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ data extraction from JuMP result objects easy.
get_core_model
build_cbm
set_bound
map_fluxes(::Array{Float64,1}, ::CobraModel)
map_fluxes(::Array{Float64,1}, ::StandardModel)
```

```@example fba
Expand Down
2 changes: 1 addition & 1 deletion docs/src/io.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ model # pretty print the model
Currently, JSON and Matlab models can be exported.

```@docs
save_model(model::CobraModel, file_location::String)
save_model(model::StandardModel, file_location::String)
```

```@example ioexample
Expand Down
30 changes: 15 additions & 15 deletions docs/src/model_construction.md
Original file line number Diff line number Diff line change
Expand Up @@ -184,10 +184,10 @@ Once you have defined some metabolites, genes, and reactions, you can construct
a model! This is most simply done by using the empty model constructor:

```@docs
CobraModel()
StandardModel()
```

The fields of `CobraModel` can then be assigned as usual.
The fields of `StandardModel` can then be assigned as usual.

```@example
using COBREXA
Expand Down Expand Up @@ -218,7 +218,7 @@ mets = [atp, adp, pp, h2o, glc, lac]
genes = [g1, g2, g3, g4]
rxns = [catabolism, anabolism, lac_ex, glc_ex]
model = CobraModel()
model = StandardModel()
model.id = "Test model"
add!(model, mets)
add!(model, rxns)
Expand All @@ -231,9 +231,9 @@ Here the `add` functions were used to add the reactions, metabolites and genes
to the model.

```@docs
add!(model::CobraModel, rxns::Array{Reaction, 1})
add!(model::CobraModel, mets::Array{Metabolite, 1})
add!(model::CobraModel, genes::Array{Gene, 1})
add!(model::StandardModel, rxns::Array{Reaction, 1})
add!(model::StandardModel, mets::Array{Metabolite, 1})
add!(model::StandardModel, genes::Array{Gene, 1})
```

Checking for duplicates of genes, metabolites or reactions can also be done.
Expand Down Expand Up @@ -281,9 +281,9 @@ Similar functionality exists for genes and reactions. Duplicate reactions,
metabolites or genes can be removed using `rm!`.

```@docs
rm!(model::CobraModel, rxns::Union{Array{Reaction, 1}, Reaction})
rm!(model::CobraModel, mets::Union{Array{Metabolite, 1}, Metabolite})
rm!(model::CobraModel, genes::Union{Array{Gene, 1}, Gene})
rm!(model::StandardModel, rxns::Union{Array{Reaction, 1}, Reaction})
rm!(model::StandardModel, mets::Union{Array{Metabolite, 1}, Metabolite})
rm!(model::StandardModel, genes::Union{Array{Gene, 1}, Gene})
```

```@example
Expand All @@ -293,7 +293,7 @@ atp2 = Metabolite("atp2")
mets = [atp, atp2]
model = CobraModel()
model = StandardModel()
add!(model, mets)
rm!(model, atp2)
Expand All @@ -304,7 +304,7 @@ relative to the reactions. `fix_model` also ensures that no extra metabolites
are present.

```@docs
fix_model!(model::CobraModel)
fix_model!(model::StandardModel)
```

```@example
Expand All @@ -318,7 +318,7 @@ anabolism.id = "anabolism"
mets = [atp]
rxns = [anabolism]
model = CobraModel()
model = StandardModel()
model.id = "Test model"
add!(model, mets) # missing adp
add!(model, rxns)
Expand All @@ -332,7 +332,7 @@ Helper functions from Base have also been overwritten to make accessing
reactions, metabolites and genes easy from a model.

```@docs
getindex(model::CobraModel, rxn::Reaction)
getindex(model::CobraModel, rxn::Metabolite)
getindex(model::CobraModel, rxn::Gene)
getindex(model::StandardModel, rxn::Reaction)
getindex(model::StandardModel, rxn::Metabolite)
getindex(model::StandardModel, rxn::Gene)
```
4 changes: 2 additions & 2 deletions docs/src/model_structure.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@
Before reading, writing, or building models, it is important to understand how
they are represented internally in `COBREXA.jl`.

Each model is a struct of the type `CobraModel`, which is composed of a model
Each model is a struct of the type `StandardModel`, which is composed of a model
`id`, and arrays of `Reaction`s, `Metabolite`s and `Gene`s.

```@docs
CobraModel
StandardModel
```

The fields of `Reaction`, `Metabolite`, `Gene` types are shown below. When
Expand Down
6 changes: 3 additions & 3 deletions src/analysis/flux_balance_analysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ function flux_balance_analysis_dict(
model::M,
args...;
kwargs...,
)::Union{Dict{String,Float64},Nothing} where {M<:MetabolicModel}
)::Union{Dict{String,Float64},Nothing} where {M <: MetabolicModel}
v = flux_balance_analysis_vec(model, args...; kwargs...)
isnothing(v) && return nothing
Dict(zip(reactions(model), v))
Expand Down Expand Up @@ -141,8 +141,8 @@ solved_model = fba(model, optimizer; modifications=[modify_objective(biomass)])
function flux_balance_analysis(
model::M,
optimizer;
modifications = [(model, opt_model) -> nothing],
) where {M<:MetabolicModel}
modifications=[(model, opt_model) -> nothing],
) where {M <: MetabolicModel}

opt_model = make_optimization_model(model, optimizer)

Expand Down
159 changes: 158 additions & 1 deletion src/analysis/loopless_flux_balance_analysis.jl
Original file line number Diff line number Diff line change
@@ -1 +1,158 @@
# TODO
function _get_boundary_reaction_ids(model::StandardModel)::Array{String,1}
return [i.id for i in model.reactions if length(i.metabolites) == 1]
end

function add_cycle_free(fluxes::Dict{String,Float64})
return (model, opt_model) -> begin
v = opt_model[:x]
old_objective = objective_function(opt_model)
boundary_ids = _get_boundary_reaction_ids(model)
min_objectives = zeros(Int64, 1, length(v))
for (i, reaction) in enumerate(model.reactions)
id = reaction.id
if v[i] in old_objective.terms.keys
continue
end
flux = fluxes[id]
if id in boundary_ids
lb = flux
ub = flux
elseif flux >= 0
lb = max(0, reaction.lb)
ub = max(flux, reaction.ub)
min_objectives[i] = 1
else
lb = min(flux, reaction.lb)
ub = min(0, reaction.ub)
min_objectives[i] = -1
end
set_bound(i, opt_model; ub=ub, lb=lb)
end
@objective(opt_model, Min, dot(min_objectives, v))
end
end


"""
add_loopless(args...)::Union{COBREXA.JuMP.Model,Nothing}
Convert an existing FBA solution to a loopless one. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function add_loopless(model::M, opt_model, fluxes::Dict{String,Float64})::Union{COBREXA.JuMP.Model,Nothing} where {M <: MetabolicModel}
constrain_objective_value(1)(model, opt_model)
add_cycle_free(fluxes)(model, opt_model)
COBREXA.JuMP.optimize!(opt_model)
termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return opt_model
end

"""
add_loopless_vec(args...)::Union{Array{Float64,1},Nothing}
Convert an existing FBA solution to a loopless one. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function add_loopless_vec(model::M, opt_model, fluxes::Dict{String,Float64})::Union{Array{Float64,1},Nothing} where {M <: MetabolicModel}
opt_model = add_loopless(model, opt_model, fluxes)
termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return value.(opt_model[:x])
end

"""
add_loopless_dict(args...)::Union{Dict{String,Float64},Nothing}
Convert an existing FBA solution to a loopless one. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function add_loopless_dict(model::M, opt_model, fluxes::Dict{String,Float64})::Union{Dict{String,Float64},Nothing} where {M <: MetabolicModel}
v = add_loopless_vec(model, opt_model, fluxes)
isnothing(v) && return nothing
return Dict(zip(reactions(model), v))
end

"""
loopless_flux_balance_analysis(args...)::Union{COBREXA.JuMP.Model,Nothing}
Run loopless FBA. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function loopless_flux_balance_analysis(
model::M,
optimizer;
modifications=[(model, opt_model) -> nothing],

)::Union{COBREXA.JuMP.Model,Nothing} where {M <: MetabolicModel}
opt_model = flux_balance_analysis(model, optimizer, modifications=modifications)
termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
fluxes = Dict(zip(reactions(model), value.(opt_model[:x])))
return add_loopless(
model,
opt_model,
fluxes,
)
end

"""
loopless_flux_balance_analysis_vec(args...)::Union{Array{Float64,1},Nothing}
Run loopless FBA. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function loopless_flux_balance_analysis_vec(
model::M,
optimizer;
modifications=[(model, opt_model) -> nothing],

)::Union{Array{Float64,1},Nothing} where {M <: MetabolicModel}
opt_model = loopless_flux_balance_analysis(model, optimizer, modifications=modifications)
termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
return value.(opt_model[:x])
end

"""
loopless_flux_balance_analysis_dict(args...)::Union{Dict{String,Float64},Nothing}
Run loopless FBA. Removes as many loops as possible using the method from [1]
References:
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi:
10.1093/bioinformatics/btv096.
"""
function loopless_flux_balance_analysis_dict(
model::M,
optimizer;
modifications=[(model, opt_model) -> nothing],

)::Union{Dict{String,Float64},Nothing} where {M <: MetabolicModel}
v = loopless_flux_balance_analysis_vec(model, optimizer, modifications=modifications)
isnothing(v) && return nothing
return Dict(zip(reactions(model), v))
end
Loading

0 comments on commit 3db6eb6

Please sign in to comment.