From ab8dfda8bc3959aad9b4138408444069589a90b4 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 11:52:43 +0100 Subject: [PATCH 01/26] ignore models --- .gitignore | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.gitignore b/.gitignore index 46e71d415..caac77fde 100644 --- a/.gitignore +++ b/.gitignore @@ -20,3 +20,8 @@ Manifest.toml # ignore container files *.sif + +# ignore models +*.xml +*.json +*.mat \ No newline at end of file From 7d99991d2e4cf56c58c6796089e0f796a1ed02c9 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 12:15:35 +0100 Subject: [PATCH 02/26] get cobrexa to load --- Project.toml | 2 + docs/Project.toml | 9 +- docs/src/examples/02-flux-balance-analysis.jl | 4 + src/COBREXA.jl | 3 + src/analysis/flux_balance_analysis.jl | 54 +++++++++++ .../parsimonious_flux_balance_analysis.jl | 92 +++++++++++++++++++ src/builders/genes.jl | 21 ++--- src/builders/objectives.jl | 8 +- 8 files changed, 172 insertions(+), 21 deletions(-) create mode 100644 src/analysis/flux_balance_analysis.jl create mode 100644 src/analysis/parsimonious_flux_balance_analysis.jl diff --git a/Project.toml b/Project.toml index dc4027e93..24409afba 100644 --- a/Project.toml +++ b/Project.toml @@ -32,6 +32,8 @@ Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" [targets] test = ["Aqua", "Clarabel", "GLPK", "Test", "Tulip"] diff --git a/docs/Project.toml b/docs/Project.toml index 1f523ccd8..528ff1adf 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,15 +1,10 @@ [deps] +AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c" COBREXA = "babc4406-5200-4a30-9033-bf5ae714c842" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" -Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" -ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -Escher = "8cc96de1-1b23-48cb-9272-618d67962629" -GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index 678930970..8d36fd1a9 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -2,5 +2,9 @@ # # Flux balance analysis using COBREXA +import AbstractFBCModels as A +import JSONFBCModels as J # TODO: run FBA on a FBC model + +model diff --git a/src/COBREXA.jl b/src/COBREXA.jl index 5baaa18cc..a223a25fb 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -16,4 +16,7 @@ include("builders/core.jl") include("builders/genes.jl") include("builders/objectives.jl") +include("analysis/flux_balance_analysis.jl") +include("analysis/parsimonious_flux_balance_analysis.jl") + end # module COBREXA diff --git a/src/analysis/flux_balance_analysis.jl b/src/analysis/flux_balance_analysis.jl new file mode 100644 index 000000000..4e48c9a19 --- /dev/null +++ b/src/analysis/flux_balance_analysis.jl @@ -0,0 +1,54 @@ +""" +$(TYPEDSIGNATURES) + +Run flux balance analysis (FBA) on the `model` optionally specifying +`modifications` to the problem. Basically, FBA solves this optimization problem: +``` +max cᵀx +s.t. S x = b + xₗ ≤ x ≤ xᵤ +``` +See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat +Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more +information. + +The `optimizer` must be set to a `JuMP`-compatible optimizer, such as +`GLPK.Optimizer` or `Tulip.Optimizer` + +Optionally, you may specify one or more modifications to be applied to the +model before the analysis, such as [`modify_optimizer_attribute`](@ref), +[`change_objective`](@ref), and [`modify_sense`](@ref). + +Returns an optimized `JuMP` model. + +# Example +``` +model = load_model("e_coli_core.json") +solution = flux_balance_analysis(model, GLPK.optimizer) +value.(solution[:x]) # extract flux steady state from the optimizer + +biomass_reaction_id = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM") + +modified_solution = flux_balance_analysis(model, GLPK.optimizer; + modifications=[modify_objective(biomass_reaction_id)]) +``` +""" +function flux_balance_analysis(model::C.ConstraintTree, optimizer; modifications = []) + opt_model = make_optimization_model(model, optimizer) + + for mod in modifications + mod(model, opt_model) + end + + optimize!(opt_model) + + ModelWithResult(model, opt_model) +end + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant of [`flux_balance_analysis`](@ref). +""" +flux_balance_analysis(optimizer; modifications = []) = + model -> flux_balance_analysis(model, optimizer; modifications) diff --git a/src/analysis/parsimonious_flux_balance_analysis.jl b/src/analysis/parsimonious_flux_balance_analysis.jl new file mode 100644 index 000000000..9e5405fe9 --- /dev/null +++ b/src/analysis/parsimonious_flux_balance_analysis.jl @@ -0,0 +1,92 @@ +""" +$(TYPEDSIGNATURES) + +Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA +runs two consecutive optimization problems. The first is traditional FBA: +``` +max cᵀx = μ +s.t. S x = b + xₗ ≤ x ≤ xᵤ +``` +And the second is a quadratic optimization problem: +``` +min Σᵢ xᵢ² +s.t. S x = b + xₗ ≤ x ≤ xᵤ + μ = μ⁰ +``` +Where the optimal solution of the FBA problem, μ⁰, has been added as an +additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M, +Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N, +Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils, +Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data +from evolved E. coli are consistent with computed optimal growth from +genome-scale models. Molecular Systems Biology, 6. 390. doi: +accession:10.1038/msb.2010.47" for more details. + +pFBA gets the model optimum by standard FBA (using +[`flux_balance_analysis`](@ref) with `optimizer` and `modifications`), then +finds a minimal total flux through the model that still satisfies the (slightly +relaxed) optimum. This is done using a quadratic problem optimizer. If the +original optimizer does not support quadratic optimization, it can be changed +using the callback in `qp_modifications`, which are applied after the FBA. See +the documentation of [`flux_balance_analysis`](@ref) for usage examples of +modifications. + +The optimum relaxation sequence can be specified in `relax` parameter, it +defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original +bound. + +Returns an optimized model that contains the pFBA solution (or an unsolved model +if something went wrong). + +# Performance + +This implementation attempts to save time by executing all pFBA steps on a +single instance of the optimization model problem, trading off possible +flexibility. For slightly less performant but much more flexible use, one can +construct parsimonious models directly using +[`with_parsimonious_objective`](@ref). + +# Example +``` +model = load_model("e_coli_core.json") +parsimonious_flux_balance_analysis(model, biomass, Gurobi.Optimizer) |> values_vec +``` +""" +function parsimonious_flux_balance_analysis( + model::C.ConstraintTree, + optimizer; + modifications = [], + qp_modifications = [], + relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99], +) + # Run FBA + opt_model = flux_balance_analysis(model, optimizer; modifications) + J.is_solved(opt_model) || return nothing # FBA failed + + # get the objective + Z = J.objective_value(opt_model) + original_objective = J.objective_function(opt_model) + + # prepare the model for pFBA + for mod in qp_modifications + mod(model, opt_model) + end + + # add the minimization constraint for total flux + v = opt_model[:x] # fluxes + J.@objective(opt_model, Min, sum(dot(v, v))) + + for rb in relax_bounds + # lb, ub = objective_bounds(rb)(Z) + J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub) + + J.optimize!(opt_model) + J.is_solved(opt_model) && break + + J.delete(opt_model, pfba_constraint) + J.unregister(opt_model, :pfba_constraint) + end + +end diff --git a/src/builders/genes.jl b/src/builders/genes.jl index 2996a08f9..17263d4ad 100644 --- a/src/builders/genes.jl +++ b/src/builders/genes.jl @@ -1,19 +1,18 @@ -knockout_constraints(; fluxes::C.ConstraintTree, knockout_test) = C.ConstraintTree( - rxn => C.Constraint(C.value(fluxes[rxn]), 0.0) for - rxn in keys(fluxes) if knockout_test(rxn) -) +knockout_constraints(;fluxes::C.ConstraintTree, knockout_test) = + C.ConstraintTree( + rxn => C.Constraint(C.value(fluxes[rxn]), 0.0) for rxn in keys(fluxes) if knockout_test(rxn) + ) export knockout_constraints -fbc_gene_knockout_constraints(; - fluxes::C.ConstraintTree, - genes, - fbc_model::A.AbstractFBCModel, -) = knockout_constraints(; +fbc_gene_knockout_constraints(;fluxes::C.ConstraintTree, genes, fbc_model::A.AbstractFBCModel) = knockout_constraints(; fluxes, - knockout_test = rxn -> - !A.reaction_gene_products_available(rxn, g -> not(g in genes)), + knockout_test = + rxn -> !A.reaction_gene_products_available( + rxn, + g -> not(g in genes) + ), ) export fbc_gene_knockout_constraints diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index 93cca18ce..b6d5f7994 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -1,6 +1,7 @@ -sum_objectve(x) = C.Constraint(sum(C.value.(x), init = zero(C.LinearValue))) - +sum_objectve(x) = + C.Constraint(sum(C.value.(x), init = zero(C.LinearValue))) + sum_objective(x::C.ConstraintTree) = squared_error_objective(values(x)) squared_sum_objective(x) = @@ -9,7 +10,8 @@ squared_sum_objective(x) = squared_sum_objective(x::C.ConstraintTree) = squared_sum_objective(values(x)) squared_error_objective(constraints::C.ConstraintTree, target) = C.Constraint( - sum(let tmp = (C.value(c) - target[k]) + sum( + let tmp = (C.value(c) - target[k]) tmp * tmp end for (k, c) in constraints if haskey(target, k)), ) From 7e82e479df4620bf8620a2f00349df6fef5057a1 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 12:23:12 +0100 Subject: [PATCH 03/26] temporarily add JSONFBCModels to make coding easier --- Project.toml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 24409afba..54e6789c9 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ConstraintTrees = "5515826b-29c3-47a5-8849-8513ac836620" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" DistributedData = "f6a0035f-c5ac-4ad0-b410-ad102ced35df" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -29,11 +30,11 @@ julia = "1.5" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" +SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" -Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" [targets] test = ["Aqua", "Clarabel", "GLPK", "Test", "Tulip"] From 5090be89cb7cc76ca0b15adbdeabda2874ff6379 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 12:26:06 +0100 Subject: [PATCH 04/26] load model and create ctmodel --- docs/src/examples/02-flux-balance-analysis.jl | 9 ++++++--- src/COBREXA.jl | 1 + src/builders/core.jl | 19 ++++++++----------- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index 8d36fd1a9..b21fb525c 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -1,10 +1,13 @@ # # Flux balance analysis -using COBREXA +import COBREXA as X import AbstractFBCModels as A import JSONFBCModels as J -# TODO: run FBA on a FBC model +model = A.load(J.JSONFBCModel, "e_coli_core.json") + + +ctmodel = X.fbc_model_constraints(model) + -model diff --git a/src/COBREXA.jl b/src/COBREXA.jl index a223a25fb..cb0bab37a 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -8,6 +8,7 @@ using DocStringExtensions import AbstractFBCModels as A import ConstraintTrees as C import JuMP as J +import SparseArrays: sparse include("types.jl") include("solver.jl") diff --git a/src/builders/core.jl b/src/builders/core.jl index d85e1d71c..221b8ee4a 100644 --- a/src/builders/core.jl +++ b/src/builders/core.jl @@ -1,20 +1,17 @@ -import AbstractFBCModels as F -import SparseArrays: sparse - """ $(TYPEDSIGNATURES) A constraint tree that models the content of the given instance of `AbstractFBCModel`. """ -function fbc_model_constraints(model::F.AbstractFBCModel) - rxns = Symbol.(F.reactions(model)) - mets = Symbol.(F.metabolites(model)) - lbs, ubs = F.bounds(model) - stoi = F.stoichiometry(model) - bal = F.balance(model) - obj = F.objective(model) +function fbc_model_constraints(model::A.AbstractFBCModel) + rxns = Symbol.(A.reactions(model)) + mets = Symbol.(A.metabolites(model)) + lbs, ubs = A.bounds(model) + stoi = A.stoichiometry(model) + bals = A.balance(model) + obj = A.objective(model) #TODO: is sparse() required below? return C.ConstraintTree( @@ -23,7 +20,7 @@ function fbc_model_constraints(model::F.AbstractFBCModel) m => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for (m, row, b) in zip(mets, eachrow(stoi), bals) ), - :objective => C.Constraint(value = C.Value(sparse(obj))), + :objective => C.Constraint(value = C.LinearValue(sparse(obj))), ) end From e877a1b51299050bf9a10cb57ea09618c0e0c74a Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 12:32:43 +0100 Subject: [PATCH 05/26] small fixes to core --- src/builders/core.jl | 50 ++++++++++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/src/builders/core.jl b/src/builders/core.jl index 221b8ee4a..67a7179eb 100644 --- a/src/builders/core.jl +++ b/src/builders/core.jl @@ -1,26 +1,29 @@ +import AbstractFBCModels as F +import SparseArrays: sparse + """ $(TYPEDSIGNATURES) A constraint tree that models the content of the given instance of `AbstractFBCModel`. """ -function fbc_model_constraints(model::A.AbstractFBCModel) - rxns = Symbol.(A.reactions(model)) - mets = Symbol.(A.metabolites(model)) - lbs, ubs = A.bounds(model) - stoi = A.stoichiometry(model) - bals = A.balance(model) - obj = A.objective(model) +function fbc_model_constraints(model::F.AbstractFBCModel) + rxns = Symbol.(F.reactions(model)) + mets = Symbol.(F.metabolites(model)) + lbs, ubs = F.bounds(model) + stoi = F.stoichiometry(model) + bal = F.balance(model) + obj = F.objective(model) #TODO: is sparse() required below? return C.ConstraintTree( - :fluxes => C.variables(keys = rxns, bounds = zip(lbs, ubs)), - :balances => C.ConstraintTree( - m => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for - (m, row, b) in zip(mets, eachrow(stoi), bals) - ), - :objective => C.Constraint(value = C.LinearValue(sparse(obj))), + :fluxes^C.variables(keys = rxns, bounds = zip(lbs, ubs)) * + :flux_stoichiometry^C.ConstraintTree( + met => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for + (met, row, b) in zip(mets, eachrow(stoi), bal) + ) * + :objective^C.Constraint(C.LinearValue(sparse(obj))), ) end @@ -67,11 +70,26 @@ sign_split_constraints(; signed::C.ConstraintTree, ) = C.ConstraintTree( k => C.Constraint( - value = s + (haskey(negative, k) ? C.value(negative[k]) : zero(C.Value)) - - (haskey(positive, k) ? C.value(positive[k]) : zero(C.Value)), + value = s.value + + (haskey(negative, k) ? negative[k].value : zero(typeof(s.value))) - + (haskey(positive, k) ? positive[k].value : zero(typeof(s.value))), bound = 0.0, - ) for (k, s) in C.elems(signed) + ) for (k, s) in signed ) #TODO the example above might as well go to docs export sign_split_constraints + +function fluxes_in_direction(fluxes::C.ConstraintTree, direction = :forward) + keys = Symbol[] + for (id, flux) in fluxes + if direction == :forward + last(flux.bound) > 0 && push!(keys, id) + else + first(flux.bound) < 0 && push!(keys, id) + end + end + C.variables(; keys, bounds = Ref((0.0, Inf))) +end + +export fluxes_in_direction From 612d5f40771acba6b96aa35a5c4148ea0b561575 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 12:34:14 +0100 Subject: [PATCH 06/26] update deps --- Project.toml | 4 ++-- docs/src/examples/02-flux-balance-analysis.jl | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 54e6789c9..068dddaf7 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ ConstraintTrees = "5515826b-29c3-47a5-8849-8513ac836620" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" DistributedData = "f6a0035f-c5ac-4ad0-b410-ad102ced35df" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -35,6 +34,7 @@ GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" +JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8" [targets] -test = ["Aqua", "Clarabel", "GLPK", "Test", "Tulip"] +test = ["Aqua", "Clarabel", "Downloads", "GLPK", "SHA", "Test", "Tulip", "JSONFBCModels"] diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index b21fb525c..f14be4bdc 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -7,7 +7,6 @@ import JSONFBCModels as J model = A.load(J.JSONFBCModel, "e_coli_core.json") - ctmodel = X.fbc_model_constraints(model) From ab3696d4249e15c36bb73393196e6163985e9612 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 13:20:23 +0100 Subject: [PATCH 07/26] fba working --- docs/src/examples/02-flux-balance-analysis.jl | 6 +++ src/COBREXA.jl | 1 + src/analysis/flux_balance_analysis.jl | 40 ++++++++++++------- src/analysis/modifications/optimizer.jl | 36 +++++++++++++++++ 4 files changed, 69 insertions(+), 14 deletions(-) create mode 100644 src/analysis/modifications/optimizer.jl diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index f14be4bdc..f4d811fab 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -4,9 +4,15 @@ import COBREXA as X import AbstractFBCModels as A import JSONFBCModels as J +import ConstraintTrees as C +using Gurobi model = A.load(J.JSONFBCModel, "e_coli_core.json") ctmodel = X.fbc_model_constraints(model) +vt = C.ValueTree(ctmodel, value.(opt_model[:x])) +vt = X.flux_balance_analysis(model, Gurobi.Optimizer) + +vt = X.flux_balance_analysis(ctmodel, Gurobi.Optimizer; modifications=[X.silence]) diff --git a/src/COBREXA.jl b/src/COBREXA.jl index cb0bab37a..934ca14b3 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -17,6 +17,7 @@ include("builders/core.jl") include("builders/genes.jl") include("builders/objectives.jl") +include("analysis/modifications/optimizer.jl") include("analysis/flux_balance_analysis.jl") include("analysis/parsimonious_flux_balance_analysis.jl") diff --git a/src/analysis/flux_balance_analysis.jl b/src/analysis/flux_balance_analysis.jl index 4e48c9a19..fcfdb34ae 100644 --- a/src/analysis/flux_balance_analysis.jl +++ b/src/analysis/flux_balance_analysis.jl @@ -1,7 +1,7 @@ """ $(TYPEDSIGNATURES) -Run flux balance analysis (FBA) on the `model` optionally specifying +Run flux balance analysis (FBA) on the `model`, optionally specifying `modifications` to the problem. Basically, FBA solves this optimization problem: ``` max cᵀx @@ -13,36 +13,46 @@ Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more information. The `optimizer` must be set to a `JuMP`-compatible optimizer, such as -`GLPK.Optimizer` or `Tulip.Optimizer` +`GLPK.Optimizer` or `Tulip.Optimizer`. Optionally, you may specify one or more modifications to be applied to the model before the analysis, such as [`modify_optimizer_attribute`](@ref), [`change_objective`](@ref), and [`modify_sense`](@ref). -Returns an optimized `JuMP` model. +Returns a [`C.ValueTree`](@ref). # Example ``` model = load_model("e_coli_core.json") solution = flux_balance_analysis(model, GLPK.optimizer) -value.(solution[:x]) # extract flux steady state from the optimizer +``` +""" +function flux_balance_analysis(model::A.AbstractFBCModel, optimizer; modifications = []) + ctmodel = fbc_model_constraints(model) + flux_balance_analysis(ctmodel, optimizer; modifications) +end -biomass_reaction_id = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM") +""" +$(TYPEDSIGNATURES) -modified_solution = flux_balance_analysis(model, GLPK.optimizer; - modifications=[modify_objective(biomass_reaction_id)]) -``` +A variant of [`flux_balance_analysis`](@ref) that takes in a +[`C.ConstraintTree`](@ref) as the model to optimize. The objective is inferred +from the field `objective` in `ctmodel`. All other arguments are forwarded. """ -function flux_balance_analysis(model::C.ConstraintTree, optimizer; modifications = []) - opt_model = make_optimization_model(model, optimizer) +function flux_balance_analysis(ctmodel::C.ConstraintTree, optimizer; modifications = []) + opt_model = optimization_model( + ctmodel; + objective = ctmodel.objective.value, + optimizer, + ) for mod in modifications - mod(model, opt_model) + mod(ctmodel, opt_model) end - optimize!(opt_model) + J.optimize!(opt_model) - ModelWithResult(model, opt_model) + C.ValueTree(ctmodel, J.value.(opt_model[:x])) end """ @@ -51,4 +61,6 @@ $(TYPEDSIGNATURES) Pipe-able variant of [`flux_balance_analysis`](@ref). """ flux_balance_analysis(optimizer; modifications = []) = - model -> flux_balance_analysis(model, optimizer; modifications) + m -> flux_balance_analysis(m, optimizer; modifications) + +export flux_balance_analysis \ No newline at end of file diff --git a/src/analysis/modifications/optimizer.jl b/src/analysis/modifications/optimizer.jl new file mode 100644 index 000000000..c29b7f7c0 --- /dev/null +++ b/src/analysis/modifications/optimizer.jl @@ -0,0 +1,36 @@ + +""" +$(TYPEDSIGNATURES) + +Change the objective sense of optimization. Possible arguments are +`JuMP.MAX_SENSE` and `JuMP.MIN_SENSE`. +""" +modify_sense(objective_sense) = + (_, opt_model) -> set_objective_sense(opt_model, objective_sense) + +""" +$(TYPEDSIGNATURES) + +Change the JuMP optimizer used to run the optimization. +""" +modify_optimizer(optimizer) = (_, opt_model) -> J.set_optimizer(opt_model, optimizer) + +""" +$(TYPEDSIGNATURES) + +Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer +to the JuMP documentation and the documentation of the specific optimizer for +usable keys and values. +""" +modify_optimizer_attribute(attribute_key, value) = + (_, opt_model) -> J.set_optimizer_attribute(opt_model, attribute_key, value) + +""" + silence + +Modification that disable all output from the JuMP optimizer (shortcut for +`set_silent` from JuMP). +""" +const silence = (_, opt_model) -> J.set_silent(opt_model) + +export modify_sense, modify_optimizer, modify_optimizer_attribute, silence From 0015fa5240f74b110cc2ba7fdbebbba1b4d41a26 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 22:00:09 +0100 Subject: [PATCH 08/26] finish docs and tests for fba --- docs/make.jl | 4 + docs/src/examples/02-flux-balance-analysis.jl | 77 ++++++++++++++++--- src/COBREXA.jl | 2 +- src/analysis/flux_balance_analysis.jl | 12 ++- .../{optimizer.jl => optimizer_settings.jl} | 10 +-- src/solver.jl | 32 -------- 6 files changed, 86 insertions(+), 51 deletions(-) rename src/analysis/modifications/{optimizer.jl => optimizer_settings.jl} (68%) diff --git a/docs/make.jl b/docs/make.jl index 0b98cc0c3..26526817c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,6 +2,10 @@ using Documenter using Literate, JSON using COBREXA +# testing constants +const TEST_TOLERANCE = 1e-3 +const QP_TEST_TOLERANCE = 1e-2 # for Clarabel + # build the examples examples_path = joinpath(@__DIR__, "src", "examples") examples_basenames = sort(filter(x -> endswith(x, ".jl"), readdir(examples_path))) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index f4d811fab..cce5e6bd5 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -1,18 +1,77 @@ -# # Flux balance analysis +# # Flux balance analysis (FBA) + +# We will use [`flux_balance_analysis`](@ref) and several related functions to +# find the optimal flux distribution in the *E. coli* "core" model. + +# If it is not already present, download the model and load the package: +import Downloads: download + +!isfile("e_coli_core.json") && + download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") + +# next, load the necessary packages import COBREXA as X -import AbstractFBCModels as A -import JSONFBCModels as J -import ConstraintTrees as C -using Gurobi +import AbstractFBCModels as A # for the accessors +import JSONFBCModels as J # for the model type +import Tulip as T # use any JuMP supported optimizer +import GLPK as G + +model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model -model = A.load(J.JSONFBCModel, "e_coli_core.json") +# run FBA on the model using default settings + +vt = X.flux_balance_analysis(model, T.Optimizer) + +@test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src + +# Alternatively, a constraint tree can be passed in as well ctmodel = X.fbc_model_constraints(model) -vt = C.ValueTree(ctmodel, value.(opt_model[:x])) +# We can also pass some modifications to the optimizer +# Except for `X.silence`, all other optimizer modifications +# are the same as those in JuMP. +vt = X.flux_balance_analysis( + ctmodel, + G.Optimizer; + modifications = [ + X.silence + X.set_objective_sense(X.J.MAX_SENSE) # JuMP is called J inside COBREXA + X.set_optimizer(T.Optimizer) # change to Tulip from GLPK + X.set_optimizer_attribute("IPM_IterationsLimit", 110) # Tulip specific setting + ], +) + +@test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src + +# We can also modify the model. The most explicit way to do this is +# to make a new constraint tree representation of the model. + +import ConstraintTrees as C + +fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value + +forced_mixed_fermentation = + ctmodel * :fermentation^C.Constraint(fermentation, (10.0, 1000.0)) # new modified model is created + +vt = X.flux_balance_analysis(forced_mixed_fermentation, T.Optimizer; modifications = [X.silence]) + +@test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src + +# Models that cannot be solved return `nothing`. In the example below, the +# underlying model is modified. + +ctmodel.fluxes.ATPM.bound = (1000.0, 10000.0) # TODO make mutable + +vt = X.flux_balance_analysis(ctmodel, T.Optimizer; modifications = [X.silence]) + +@test isnothing(vt) #src + +# Models can also be piped into the analysis functions -vt = X.flux_balance_analysis(model, Gurobi.Optimizer) +ctmodel.fluxes.ATPM.bound = (8.39, 10000.0) # revert +vt = ctmodel |> X.flux_balance_analysis(T.Optimizer; modifications = [X.silence]) -vt = X.flux_balance_analysis(ctmodel, Gurobi.Optimizer; modifications=[X.silence]) +@test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src diff --git a/src/COBREXA.jl b/src/COBREXA.jl index 934ca14b3..5ce9d2c6a 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -17,7 +17,7 @@ include("builders/core.jl") include("builders/genes.jl") include("builders/objectives.jl") -include("analysis/modifications/optimizer.jl") +include("analysis/modifications/optimizer_settings.jl") include("analysis/flux_balance_analysis.jl") include("analysis/parsimonious_flux_balance_analysis.jl") diff --git a/src/analysis/flux_balance_analysis.jl b/src/analysis/flux_balance_analysis.jl index fcfdb34ae..2de61855c 100644 --- a/src/analysis/flux_balance_analysis.jl +++ b/src/analysis/flux_balance_analysis.jl @@ -2,7 +2,8 @@ $(TYPEDSIGNATURES) Run flux balance analysis (FBA) on the `model`, optionally specifying -`modifications` to the problem. Basically, FBA solves this optimization problem: +`modifications` to the problem. Basically, FBA solves this optimization +problem: ``` max cᵀx s.t. S x = b @@ -15,9 +16,10 @@ information. The `optimizer` must be set to a `JuMP`-compatible optimizer, such as `GLPK.Optimizer` or `Tulip.Optimizer`. -Optionally, you may specify one or more modifications to be applied to the -model before the analysis, such as [`modify_optimizer_attribute`](@ref), -[`change_objective`](@ref), and [`modify_sense`](@ref). +Optionally, you may specify one or more modifications to be applied to the model +before the analysis, such as [`set_objective_sense`](@ref), +[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and +[`silence`](@ref). Returns a [`C.ValueTree`](@ref). @@ -51,6 +53,8 @@ function flux_balance_analysis(ctmodel::C.ConstraintTree, optimizer; modificatio end J.optimize!(opt_model) + + is_solved(opt_model) || return nothing C.ValueTree(ctmodel, J.value.(opt_model[:x])) end diff --git a/src/analysis/modifications/optimizer.jl b/src/analysis/modifications/optimizer_settings.jl similarity index 68% rename from src/analysis/modifications/optimizer.jl rename to src/analysis/modifications/optimizer_settings.jl index c29b7f7c0..ed8c8a664 100644 --- a/src/analysis/modifications/optimizer.jl +++ b/src/analysis/modifications/optimizer_settings.jl @@ -5,15 +5,15 @@ $(TYPEDSIGNATURES) Change the objective sense of optimization. Possible arguments are `JuMP.MAX_SENSE` and `JuMP.MIN_SENSE`. """ -modify_sense(objective_sense) = - (_, opt_model) -> set_objective_sense(opt_model, objective_sense) +set_objective_sense(objective_sense) = + (_, opt_model) -> J.set_objective_sense(opt_model, objective_sense) """ $(TYPEDSIGNATURES) Change the JuMP optimizer used to run the optimization. """ -modify_optimizer(optimizer) = (_, opt_model) -> J.set_optimizer(opt_model, optimizer) +set_optimizer(optimizer) = (_, opt_model) -> J.set_optimizer(opt_model, optimizer) """ $(TYPEDSIGNATURES) @@ -22,7 +22,7 @@ Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer to the JuMP documentation and the documentation of the specific optimizer for usable keys and values. """ -modify_optimizer_attribute(attribute_key, value) = +set_optimizer_attribute(attribute_key, value) = (_, opt_model) -> J.set_optimizer_attribute(opt_model, attribute_key, value) """ @@ -33,4 +33,4 @@ Modification that disable all output from the JuMP optimizer (shortcut for """ const silence = (_, opt_model) -> J.set_silent(opt_model) -export modify_sense, modify_optimizer, modify_optimizer_attribute, silence +export set_objective_sense, set_optimizer, set_optimizer_attribute, silence diff --git a/src/solver.jl b/src/solver.jl index 7cd7d1027..ff6fcb7b2 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -48,35 +48,3 @@ is_solved(opt_model::J.Model) = export is_solved -""" -$(TYPEDSIGNATURES) - -The optimized objective value of a JuMP model, if solved. -""" -optimized_objective_value(opt_model::J.Model)::Maybe{Float64} = - is_solved(opt_model) ? J.objective_value(opt_model) : nothing - -export optimized_objective_value - -""" -$(TYPEDSIGNATURES) - -The optimized variable assignment of a JuMP model, if solved. -""" -optimized_variable_assignment(opt_model::J.Model)::Maybe{Vector{Float64}} = - is_solved(opt_model) ? J.value.(opt_model[:x]) : nothing - -export optimized_variable_assignment - -""" -$(TYPEDSIGNATURES) - -Annotate a `ConstraintTree` with the values given by the optimization model, -producing a `ValueTree` (if solved). -""" -solution(c::C.ConstraintTree, opt_model::J.Model)::Maybe{C.ValueTree} = - let vars = optimized_variable_assignment(opt_model) - isnothing(vars) ? nothing : C.ValueTree(c, vars) - end - -export solution From 7eda0e958cd1af62b50a6b95c1a6ed9f81f4d60a Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Fri, 10 Nov 2023 22:29:36 +0100 Subject: [PATCH 09/26] add knockout functionality and doctest --- docs/src/examples/02-flux-balance-analysis.jl | 13 ++++- src/builders/genes.jl | 50 ++++++++++++++----- 2 files changed, 49 insertions(+), 14 deletions(-) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index cce5e6bd5..079ad534a 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -56,7 +56,11 @@ fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value forced_mixed_fermentation = ctmodel * :fermentation^C.Constraint(fermentation, (10.0, 1000.0)) # new modified model is created -vt = X.flux_balance_analysis(forced_mixed_fermentation, T.Optimizer; modifications = [X.silence]) +vt = X.flux_balance_analysis( + forced_mixed_fermentation, + T.Optimizer; + modifications = [X.silence], +) @test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src @@ -75,3 +79,10 @@ ctmodel.fluxes.ATPM.bound = (8.39, 10000.0) # revert vt = ctmodel |> X.flux_balance_analysis(T.Optimizer; modifications = [X.silence]) @test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src + +# Gene knockouts can be done with ease making use of the piping functionality. +# Here oxidative phosphorylation is knocked out. + +vt = ctmodel |> X.knockout!(["b0979", "b0734"], model) |> X.flux_balance_analysis(T.Optimizer; modifications = [X.silence]) + +@test isapprox(vt.objective, 0.21166, atol = TEST_TOLERANCE) #src diff --git a/src/builders/genes.jl b/src/builders/genes.jl index 17263d4ad..df4dfcd9d 100644 --- a/src/builders/genes.jl +++ b/src/builders/genes.jl @@ -1,18 +1,42 @@ -knockout_constraints(;fluxes::C.ConstraintTree, knockout_test) = - C.ConstraintTree( - rxn => C.Constraint(C.value(fluxes[rxn]), 0.0) for rxn in keys(fluxes) if knockout_test(rxn) - ) - -export knockout_constraints +""" +$(TYPEDSIGNATURES) +""" +knockout_constraints(; fluxes::C.ConstraintTree, knockout_test) = C.ConstraintTree( + rxn => C.Constraint(C.value(fluxes[rxn]), 0.0) for + rxn in keys(fluxes) if knockout_test(rxn) +) -fbc_gene_knockout_constraints(;fluxes::C.ConstraintTree, genes, fbc_model::A.AbstractFBCModel) = knockout_constraints(; +""" +$(TYPEDSIGNATURES) +""" +gene_knockouts(; + fluxes::C.ConstraintTree, + ko_genes::Vector{String}, + model::A.AbstractFBCModel, +) = knockout_constraints(; fluxes, - knockout_test = - rxn -> !A.reaction_gene_products_available( - rxn, - g -> not(g in genes) - ), + knockout_test = rxn -> begin + maybe_avail = A.reaction_gene_products_available( + model, + string(rxn), + g -> !(g in ko_genes), # not available if knocked out + ) + isnothing(maybe_avail) ? false : !maybe_avail # negate here because of knockout_constraints + end, ) -export fbc_gene_knockout_constraints +""" +$(TYPEDSIGNATURES) +""" +knockout!(ctmodel::C.ConstraintTree, ko_genes::Vector{String}, model::A.AbstractFBCModel) = + ctmodel * :gene_knockouts^gene_knockouts(; fluxes = ctmodel.fluxes, ko_genes, model) + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant. +""" +knockout!(ko_genes::Vector{String}, model::A.AbstractFBCModel) = + m -> knockout!(m, ko_genes, model) + From 6e673416cefe9fdf6c8607b137f1821ef773011b Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Sat, 11 Nov 2023 09:21:57 +0100 Subject: [PATCH 10/26] basic pFBA and MOMA --- docs/src/examples/03-qp-problems.jl | 88 +++++++++++++++++++ src/COBREXA.jl | 2 + .../minimize_metabolic_adjustment_analysis.jl | 83 +++++++++++++++++ .../parsimonious_flux_balance_analysis.jl | 52 +++++++---- src/builders/objectives.jl | 29 +++--- 5 files changed, 222 insertions(+), 32 deletions(-) create mode 100644 docs/src/examples/03-qp-problems.jl create mode 100644 src/analysis/minimize_metabolic_adjustment_analysis.jl diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl new file mode 100644 index 000000000..a835bae0b --- /dev/null +++ b/docs/src/examples/03-qp-problems.jl @@ -0,0 +1,88 @@ + +# # Quandratic objective flux balance analysis type problems + +# We will use [`parsimonious_flux_balance_analysis`](@ref) and +# [`minimize_metabolic_adjustment_analysis`](@ref) to find the optimal flux +# distribution in the *E. coli* "core" model. + +# If it is not already present, download the model and load the package: +import Downloads: download + +!isfile("e_coli_core.json") && + download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") + +# next, load the necessary packages + +import COBREXA as X +import AbstractFBCModels as A # for the accessors +import JSONFBCModels as J # for the model type +import Clarabel # can solve QPs + +model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model + +# Use the convenience function to run standard pFBA + +vt = X.parsimonious_flux_balance_analysis(model, Clarabel.Optimizer) + +# Or use the piping functionality + +model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; modifications = [X.silence]) + +@test isapprox(vt.objective, 0.87392; atol=TEST_TOLERANCE) #src +@test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol=QP_TEST_TOLERANCE) #src + +# Alternatively, you can construct your own constraint tree model with +# the quadratic objective (this approach is much more flexible). + +ctmodel = X.fbc_model_constraints(model) +ctmodel *= :l2objective ^ X.squared_sum_objective(ctmodel.fluxes) +ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks + +opt_model = X.optimization_model( + ctmodel; + objective = ctmodel.:l2objective.value, + optimizer = Clarabel.Optimizer, + sense = X.J.MIN_SENSE, +) + +X.J.optimize!(opt_model) # JuMP is called J in COBREXA + +X.is_solved(opt_model) # check if solved + +vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA + +@test isapprox(vt.l2objective, ?; atol=QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds + +# It is likewise as simple to run MOMA using the convenience functions. + +ref_sol = Dict("ATPS4r" => 33.0, "CYTBD" => 22.0) + +vt = X.minimize_metabolic_adjustment_analysis(model, ref_sol, Gurobi.Optimizer) + +# Or use the piping functionality + +model |> X.minimize_metabolic_adjustment_analysis(ref_sol, Clarabel.Optimizer; modifications = [X.silence]) + +@test isapprox(vt.:momaobjective, 0.81580806; atol=TEST_TOLERANCE) #src + +# Alternatively, you can construct your own constraint tree model with +# the quadratic objective (this approach is much more flexible). + +ctmodel = X.fbc_model_constraints(model) +ctmodel *= :minoxphospho ^ X.squared_sum_error_objective(ctmodel.fluxes, Dict(:ATPS4r => 33.0, :CYTBD => 22.0,)) +ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks + +opt_model = X.optimization_model( + ctmodel; + objective = ctmodel.minoxphospho.value, + optimizer = Clarabel.Optimizer, + sense = X.J.MIN_SENSE, +) + +X.J.optimize!(opt_model) # JuMP is called J in COBREXA + +X.is_solved(opt_model) # check if solved + +vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA + +@test isapprox(vt.l2objective, ?; atol=QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds diff --git a/src/COBREXA.jl b/src/COBREXA.jl index 5ce9d2c6a..fcdc70e8c 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -9,6 +9,7 @@ import AbstractFBCModels as A import ConstraintTrees as C import JuMP as J import SparseArrays: sparse +import LinearAlgebra: dot include("types.jl") include("solver.jl") @@ -20,5 +21,6 @@ include("builders/objectives.jl") include("analysis/modifications/optimizer_settings.jl") include("analysis/flux_balance_analysis.jl") include("analysis/parsimonious_flux_balance_analysis.jl") +include("analysis/minimize_metabolic_adjustment_analysis.jl") end # module COBREXA diff --git a/src/analysis/minimize_metabolic_adjustment_analysis.jl b/src/analysis/minimize_metabolic_adjustment_analysis.jl new file mode 100644 index 000000000..3f9370075 --- /dev/null +++ b/src/analysis/minimize_metabolic_adjustment_analysis.jl @@ -0,0 +1,83 @@ + +""" +$(TYPEDSIGNATURES) + +Run minimization of metabolic adjustment (MOMA) on `model` with respect to +`reference_solution`, which is a dictionary of fluxes. MOMA finds the shortest +Euclidian distance between `reference_solution` and `model` with `modifications`: +``` +min Σᵢ (xᵢ - flux_refᵢ)² +s.t. S x = b + xₗ ≤ x ≤ xᵤ +``` +Because the problem has a quadratic objective, a QP solver is required. See +"Daniel, Vitkup & Church, Analysis of Optimality in Natural and Perturbed +Metabolic Networks, Proceedings of the National Academy of Sciences, 2002" for +more details. + +Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be +found. + +# Example +``` +model = load_model("e_coli_core.json") + +``` +""" +function minimize_metabolic_adjustment_analysis( + ctmodel::C.ConstraintTree, + reference_solution::Dict{String,Float64}, + optimizer; + modifications = [], +) + _ctmodel = + ctmodel * + :momaobjective^squared_sum_error_objective( + ctmodel.fluxes, + Dict(Symbol(k) => float(v) for (k, v) in reference_solution), + ) + + opt_model = optimization_model( + _ctmodel; + objective = _ctmodel.momaobjective.value, + optimizer, + sense = J.MIN_SENSE, + ) + + for mod in modifications + mod(ctmodel, opt_model) + end + + J.optimize!(opt_model) + + is_solved(opt_model) || return nothing + + C.ValueTree(_ctmodel, J.value.(opt_model[:x])) +end + +""" +$(TYPEDSIGNATURES) + +Variant that takes an [`A.AbstractFBCModel`](@ref) as input. All other arguments are forwarded. +""" +function minimize_metabolic_adjustment_analysis( + model::A.AbstractFBCModel, + args...; + kwargs..., +) + ctmodel = fbc_model_constraints(model) + minimize_metabolic_adjustment_analysis(ctmodel, args...; kwargs...) +end + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant of [`minimize_metabolic_adjustment_analysis`](@ref). +""" +minimize_metabolic_adjustment_analysis( + reference_solution::Dict{String,Float64}, + optimizer; + kwargs..., +) = m -> minimize_metabolic_adjustment_analysis(m, reference_solution, optimizer; kwargs...) + +export minimize_metabolic_adjustment_analysis diff --git a/src/analysis/parsimonious_flux_balance_analysis.jl b/src/analysis/parsimonious_flux_balance_analysis.jl index 9e5405fe9..f4c3e9777 100644 --- a/src/analysis/parsimonious_flux_balance_analysis.jl +++ b/src/analysis/parsimonious_flux_balance_analysis.jl @@ -37,33 +37,31 @@ The optimum relaxation sequence can be specified in `relax` parameter, it defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original bound. -Returns an optimized model that contains the pFBA solution (or an unsolved model -if something went wrong). - -# Performance - -This implementation attempts to save time by executing all pFBA steps on a -single instance of the optimization model problem, trading off possible -flexibility. For slightly less performant but much more flexible use, one can -construct parsimonious models directly using -[`with_parsimonious_objective`](@ref). +Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be found. # Example ``` model = load_model("e_coli_core.json") -parsimonious_flux_balance_analysis(model, biomass, Gurobi.Optimizer) |> values_vec +parsimonious_flux_balance_analysis(model, Gurobi.Optimizer) ``` """ function parsimonious_flux_balance_analysis( - model::C.ConstraintTree, + ctmodel::C.ConstraintTree, optimizer; modifications = [], qp_modifications = [], relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99], ) # Run FBA - opt_model = flux_balance_analysis(model, optimizer; modifications) - J.is_solved(opt_model) || return nothing # FBA failed + opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) + + for mod in modifications + mod(ctmodel, opt_model) + end + + J.optimize!(opt_model) + + is_solved(opt_model) || return nothing # get the objective Z = J.objective_value(opt_model) @@ -79,14 +77,36 @@ function parsimonious_flux_balance_analysis( J.@objective(opt_model, Min, sum(dot(v, v))) for rb in relax_bounds - # lb, ub = objective_bounds(rb)(Z) + lb, ub = objective_bounds(rb)(Z) J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub) J.optimize!(opt_model) - J.is_solved(opt_model) && break + is_solved(opt_model) && break + @warn("Relaxing pFBA objective bound!") J.delete(opt_model, pfba_constraint) J.unregister(opt_model, :pfba_constraint) end + C.ValueTree(ctmodel, J.value.(opt_model[:x])) +end + +""" +$(TYPEDSIGNATURES) + +Variant that takes an [`A.AbstractFBCModel`](@ref) as input. All other arguments are forwarded. +""" +function parsimonious_flux_balance_analysis(model::A.AbstractFBCModel, optimizer; kwargs...) + ctmodel = fbc_model_constraints(model) + parsimonious_flux_balance_analysis(ctmodel, optimizer; kwargs...) end + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant of [`parsimonious_flux_balance_analysis`](@ref). +""" +parsimonious_flux_balance_analysis(optimizer; modifications = []) = + m -> parsimonious_flux_balance_analysis(m, optimizer; modifications) + +export parsimonious_flux_balance_analysis diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index b6d5f7994..7bc9075f5 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -1,19 +1,16 @@ -sum_objectve(x) = - C.Constraint(sum(C.value.(x), init = zero(C.LinearValue))) - -sum_objective(x::C.ConstraintTree) = squared_error_objective(values(x)) +squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = + C.Constraint( + sum( + (C.value(c) - target[k]) * (C.value(c) - target[k]) for + (k, c) in constraints if haskey(target, k) + ), + ) -squared_sum_objective(x) = - C.Constraint(sum(C.squared.(C.value.(x)), init = zero(C.QuadraticValue))) +squared_sum_objective(x::C.ConstraintTree) = + squared_sum_error_objective(x, Dict(keys(x) .=> 0.0)) -squared_sum_objective(x::C.ConstraintTree) = squared_sum_objective(values(x)) - -squared_error_objective(constraints::C.ConstraintTree, target) = C.Constraint( - sum( - let tmp = (C.value(c) - target[k]) - tmp * tmp - end for (k, c) in constraints if haskey(target, k)), -) - -# TODO use `mergewith` to do this reasonably (add it to ConstraintTrees) +objective_bounds(tolerance) = z -> begin + vs = (z * tolerance, z / tolerance) + (minimum(vs), maximum(vs)) +end From 9f886af765d83eaac588aeb3f93e1f14025634d9 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Sat, 11 Nov 2023 09:30:29 +0100 Subject: [PATCH 11/26] make pFBA more in line with MOMA and CTs --- .../parsimonious_flux_balance_analysis.jl | 41 ++++++++----------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/src/analysis/parsimonious_flux_balance_analysis.jl b/src/analysis/parsimonious_flux_balance_analysis.jl index f4c3e9777..162d0e97a 100644 --- a/src/analysis/parsimonious_flux_balance_analysis.jl +++ b/src/analysis/parsimonious_flux_balance_analysis.jl @@ -33,11 +33,11 @@ using the callback in `qp_modifications`, which are applied after the FBA. See the documentation of [`flux_balance_analysis`](@ref) for usage examples of modifications. -The optimum relaxation sequence can be specified in `relax` parameter, it -defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original -bound. +The `relax` parameter, which defaults to `0.999`, relaxes the objective bound of +the original bound to prevent numerical issues. -Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be found. +Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be +found. # Example ``` @@ -50,7 +50,7 @@ function parsimonious_flux_balance_analysis( optimizer; modifications = [], qp_modifications = [], - relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99], + relax = 0.999, ) # Run FBA opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) @@ -65,30 +65,25 @@ function parsimonious_flux_balance_analysis( # get the objective Z = J.objective_value(opt_model) - original_objective = J.objective_function(opt_model) + + _ctmodel = ctmodel * :pfbaobjective ^ squared_sum_objective(ctmodel.fluxes) + _ctmodel.objective.bound = relax == 1.0 ? Z : objective_bounds(relax)(Z) # TODO currently breaks + + opt_model = X.optimization_model( + _ctmodel; + objective = ctmodel.pfbaobjective.value, + optimizer, + sense = X.J.MIN_SENSE, + ) # prepare the model for pFBA for mod in qp_modifications mod(model, opt_model) end + + J.optimize!(opt_model) - # add the minimization constraint for total flux - v = opt_model[:x] # fluxes - J.@objective(opt_model, Min, sum(dot(v, v))) - - for rb in relax_bounds - lb, ub = objective_bounds(rb)(Z) - J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub) - - J.optimize!(opt_model) - is_solved(opt_model) && break - @warn("Relaxing pFBA objective bound!") - - J.delete(opt_model, pfba_constraint) - J.unregister(opt_model, :pfba_constraint) - end - - C.ValueTree(ctmodel, J.value.(opt_model[:x])) + C.ValueTree(_ctmodel, J.value.(opt_model[:x])) end """ From cc7cd5e7fb01ef37ededaff44da9c6f1dd6a75b7 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Fri, 8 Dec 2023 16:00:01 +0100 Subject: [PATCH 12/26] format --- docs/src/examples/03-qp-problems.jl | 28 ++++++++++++------- src/analysis/flux_balance_analysis.jl | 10 ++----- .../parsimonious_flux_balance_analysis.jl | 6 ++-- src/builders/genes.jl | 2 +- 4 files changed, 25 insertions(+), 21 deletions(-) diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index a835bae0b..ea8592c4a 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -28,14 +28,14 @@ vt = X.parsimonious_flux_balance_analysis(model, Clarabel.Optimizer) model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; modifications = [X.silence]) -@test isapprox(vt.objective, 0.87392; atol=TEST_TOLERANCE) #src -@test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol=QP_TEST_TOLERANCE) #src +@test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src +@test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol = QP_TEST_TOLERANCE) #src # Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). ctmodel = X.fbc_model_constraints(model) -ctmodel *= :l2objective ^ X.squared_sum_objective(ctmodel.fluxes) +ctmodel *= :l2objective^X.squared_sum_objective(ctmodel.fluxes) ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks opt_model = X.optimization_model( @@ -46,12 +46,12 @@ opt_model = X.optimization_model( ) X.J.optimize!(opt_model) # JuMP is called J in COBREXA - + X.is_solved(opt_model) # check if solved vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA -@test isapprox(vt.l2objective, ?; atol=QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds +@test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds # It is likewise as simple to run MOMA using the convenience functions. @@ -61,15 +61,23 @@ vt = X.minimize_metabolic_adjustment_analysis(model, ref_sol, Gurobi.Optimizer) # Or use the piping functionality -model |> X.minimize_metabolic_adjustment_analysis(ref_sol, Clarabel.Optimizer; modifications = [X.silence]) +model |> X.minimize_metabolic_adjustment_analysis( + ref_sol, + Clarabel.Optimizer; + modifications = [X.silence], +) -@test isapprox(vt.:momaobjective, 0.81580806; atol=TEST_TOLERANCE) #src +@test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src # Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). ctmodel = X.fbc_model_constraints(model) -ctmodel *= :minoxphospho ^ X.squared_sum_error_objective(ctmodel.fluxes, Dict(:ATPS4r => 33.0, :CYTBD => 22.0,)) +ctmodel *= + :minoxphospho^X.squared_sum_error_objective( + ctmodel.fluxes, + Dict(:ATPS4r => 33.0, :CYTBD => 22.0), + ) ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks opt_model = X.optimization_model( @@ -80,9 +88,9 @@ opt_model = X.optimization_model( ) X.J.optimize!(opt_model) # JuMP is called J in COBREXA - + X.is_solved(opt_model) # check if solved vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA -@test isapprox(vt.l2objective, ?; atol=QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds +@test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds diff --git a/src/analysis/flux_balance_analysis.jl b/src/analysis/flux_balance_analysis.jl index 2de61855c..343537002 100644 --- a/src/analysis/flux_balance_analysis.jl +++ b/src/analysis/flux_balance_analysis.jl @@ -42,18 +42,14 @@ A variant of [`flux_balance_analysis`](@ref) that takes in a from the field `objective` in `ctmodel`. All other arguments are forwarded. """ function flux_balance_analysis(ctmodel::C.ConstraintTree, optimizer; modifications = []) - opt_model = optimization_model( - ctmodel; - objective = ctmodel.objective.value, - optimizer, - ) + opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) for mod in modifications mod(ctmodel, opt_model) end J.optimize!(opt_model) - + is_solved(opt_model) || return nothing C.ValueTree(ctmodel, J.value.(opt_model[:x])) @@ -67,4 +63,4 @@ Pipe-able variant of [`flux_balance_analysis`](@ref). flux_balance_analysis(optimizer; modifications = []) = m -> flux_balance_analysis(m, optimizer; modifications) -export flux_balance_analysis \ No newline at end of file +export flux_balance_analysis diff --git a/src/analysis/parsimonious_flux_balance_analysis.jl b/src/analysis/parsimonious_flux_balance_analysis.jl index 162d0e97a..e52e412a7 100644 --- a/src/analysis/parsimonious_flux_balance_analysis.jl +++ b/src/analysis/parsimonious_flux_balance_analysis.jl @@ -65,8 +65,8 @@ function parsimonious_flux_balance_analysis( # get the objective Z = J.objective_value(opt_model) - - _ctmodel = ctmodel * :pfbaobjective ^ squared_sum_objective(ctmodel.fluxes) + + _ctmodel = ctmodel * :pfbaobjective^squared_sum_objective(ctmodel.fluxes) _ctmodel.objective.bound = relax == 1.0 ? Z : objective_bounds(relax)(Z) # TODO currently breaks opt_model = X.optimization_model( @@ -80,7 +80,7 @@ function parsimonious_flux_balance_analysis( for mod in qp_modifications mod(model, opt_model) end - + J.optimize!(opt_model) C.ValueTree(_ctmodel, J.value.(opt_model[:x])) diff --git a/src/builders/genes.jl b/src/builders/genes.jl index ae0e27368..8cd1587dd 100644 --- a/src/builders/genes.jl +++ b/src/builders/genes.jl @@ -16,7 +16,7 @@ gene_knockouts(; model::A.AbstractFBCModel, ) = knockout_constraints(; fluxes, - knockout_test = rxn -> begin + knockout_test = rxn -> begin maybe_avail = A.reaction_gene_products_available( model, string(rxn), From 800c4718743fe995c47d89bb9f00d22d09ae5cc4 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Fri, 8 Dec 2023 19:35:29 +0100 Subject: [PATCH 13/26] simplify stuff --- src/COBREXA.jl | 5 +- src/analysis/flux_balance.jl | 34 ++---- src/analysis/flux_balance_analysis.jl | 66 ----------- ...sis.jl => minimal_metabolic_adjustment.jl} | 41 +++---- .../modifications/optimizer_settings.jl | 36 ------ .../parsimonious_flux_balance_analysis.jl | 107 ------------------ src/builders/objectives.jl | 21 ++-- src/misc/bounds.jl | 20 ++++ src/solver.jl | 23 ++++ 9 files changed, 81 insertions(+), 272 deletions(-) delete mode 100644 src/analysis/flux_balance_analysis.jl rename src/analysis/{minimize_metabolic_adjustment_analysis.jl => minimal_metabolic_adjustment.jl} (59%) delete mode 100644 src/analysis/modifications/optimizer_settings.jl delete mode 100644 src/analysis/parsimonious_flux_balance_analysis.jl create mode 100644 src/misc/bounds.jl diff --git a/src/COBREXA.jl b/src/COBREXA.jl index f99f7358f..b03902860 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -9,7 +9,6 @@ import AbstractFBCModels as A import ConstraintTrees as C import JuMP as J import SparseArrays: sparse -import LinearAlgebra: dot include("types.jl") include("io.jl") @@ -22,6 +21,8 @@ include("builders/objectives.jl") include("analysis/modifications.jl") include("analysis/flux_balance.jl") include("analysis/parsimonious_flux_balance.jl") -include("analysis/minimize_metabolic_adjustment.jl") +include("analysis/minimal_metabolic_adjustment.jl") + +include("misc/bounds.jl") end # module COBREXA diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index f755b45fc..bb58d8dbf 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -29,36 +29,20 @@ model = load_model("e_coli_core.json") solution = flux_balance(model, GLPK.optimizer) ``` """ -function flux_balance(model::A.AbstractFBCModel, optimizer; modifications = []) - ctmodel = fbc_model_constraints(model) - flux_balance(ctmodel, optimizer; modifications) +function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...) + constraints = fbc_model_constraints(model) + optimize_constraints( + constraints; + objective = constraints.objective, + optimizer, + kwargs..., + ) end """ $(TYPEDSIGNATURES) -A variant of [`flux_balance`](@ref) that takes in a -[`C.ConstraintTree`](@ref) as the model to optimize. The objective is inferred -from the field `objective` in `ctmodel`. All other arguments are forwarded. -""" -function flux_balance(ctmodel::C.ConstraintTree, optimizer; modifications = []) - opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) - - for mod in modifications - mod(opt_model) - end - - J.optimize!(opt_model) - - is_solved(opt_model) || return nothing - - C.ValueTree(ctmodel, J.value.(opt_model[:x])) -end - -""" -$(TYPEDSIGNATURES) - -Pipe-able variant of [`flux_balance`](@ref). +Pipe-able overload of [`flux_balance`](@ref). """ flux_balance(optimizer; modifications = []) = m -> flux_balance(m, optimizer; modifications) diff --git a/src/analysis/flux_balance_analysis.jl b/src/analysis/flux_balance_analysis.jl deleted file mode 100644 index 343537002..000000000 --- a/src/analysis/flux_balance_analysis.jl +++ /dev/null @@ -1,66 +0,0 @@ -""" -$(TYPEDSIGNATURES) - -Run flux balance analysis (FBA) on the `model`, optionally specifying -`modifications` to the problem. Basically, FBA solves this optimization -problem: -``` -max cᵀx -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat -Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more -information. - -The `optimizer` must be set to a `JuMP`-compatible optimizer, such as -`GLPK.Optimizer` or `Tulip.Optimizer`. - -Optionally, you may specify one or more modifications to be applied to the model -before the analysis, such as [`set_objective_sense`](@ref), -[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and -[`silence`](@ref). - -Returns a [`C.ValueTree`](@ref). - -# Example -``` -model = load_model("e_coli_core.json") -solution = flux_balance_analysis(model, GLPK.optimizer) -``` -""" -function flux_balance_analysis(model::A.AbstractFBCModel, optimizer; modifications = []) - ctmodel = fbc_model_constraints(model) - flux_balance_analysis(ctmodel, optimizer; modifications) -end - -""" -$(TYPEDSIGNATURES) - -A variant of [`flux_balance_analysis`](@ref) that takes in a -[`C.ConstraintTree`](@ref) as the model to optimize. The objective is inferred -from the field `objective` in `ctmodel`. All other arguments are forwarded. -""" -function flux_balance_analysis(ctmodel::C.ConstraintTree, optimizer; modifications = []) - opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) - - for mod in modifications - mod(ctmodel, opt_model) - end - - J.optimize!(opt_model) - - is_solved(opt_model) || return nothing - - C.ValueTree(ctmodel, J.value.(opt_model[:x])) -end - -""" -$(TYPEDSIGNATURES) - -Pipe-able variant of [`flux_balance_analysis`](@ref). -""" -flux_balance_analysis(optimizer; modifications = []) = - m -> flux_balance_analysis(m, optimizer; modifications) - -export flux_balance_analysis diff --git a/src/analysis/minimize_metabolic_adjustment_analysis.jl b/src/analysis/minimal_metabolic_adjustment.jl similarity index 59% rename from src/analysis/minimize_metabolic_adjustment_analysis.jl rename to src/analysis/minimal_metabolic_adjustment.jl index 3f9370075..d586a3aa8 100644 --- a/src/analysis/minimize_metabolic_adjustment_analysis.jl +++ b/src/analysis/minimal_metabolic_adjustment.jl @@ -24,35 +24,35 @@ model = load_model("e_coli_core.json") ``` """ -function minimize_metabolic_adjustment_analysis( - ctmodel::C.ConstraintTree, +function minimal_metabolic_adjustment( + constraints::C.ConstraintTree, reference_solution::Dict{String,Float64}, optimizer; modifications = [], ) - _ctmodel = - ctmodel * + _constraints = + constraints * :momaobjective^squared_sum_error_objective( - ctmodel.fluxes, + constraints.fluxes, Dict(Symbol(k) => float(v) for (k, v) in reference_solution), ) opt_model = optimization_model( - _ctmodel; - objective = _ctmodel.momaobjective.value, + _constraints; + objective = _constraints.momaobjective.value, optimizer, sense = J.MIN_SENSE, ) for mod in modifications - mod(ctmodel, opt_model) + mod(constraints, opt_model) end J.optimize!(opt_model) is_solved(opt_model) || return nothing - C.ValueTree(_ctmodel, J.value.(opt_model[:x])) + C.ValueTree(_constraints, J.value.(opt_model[:x])) end """ @@ -60,24 +60,9 @@ $(TYPEDSIGNATURES) Variant that takes an [`A.AbstractFBCModel`](@ref) as input. All other arguments are forwarded. """ -function minimize_metabolic_adjustment_analysis( - model::A.AbstractFBCModel, - args...; - kwargs..., -) - ctmodel = fbc_model_constraints(model) - minimize_metabolic_adjustment_analysis(ctmodel, args...; kwargs...) +function minimal_metabolic_adjustment(model::A.AbstractFBCModel, args...; kwargs...) + constraints = fbc_model_constraints(model) + minimal_metabolic_adjustment(constraints, args...; kwargs...) end -""" -$(TYPEDSIGNATURES) - -Pipe-able variant of [`minimize_metabolic_adjustment_analysis`](@ref). -""" -minimize_metabolic_adjustment_analysis( - reference_solution::Dict{String,Float64}, - optimizer; - kwargs..., -) = m -> minimize_metabolic_adjustment_analysis(m, reference_solution, optimizer; kwargs...) - -export minimize_metabolic_adjustment_analysis +export minimal_metabolic_adjustment diff --git a/src/analysis/modifications/optimizer_settings.jl b/src/analysis/modifications/optimizer_settings.jl deleted file mode 100644 index ed8c8a664..000000000 --- a/src/analysis/modifications/optimizer_settings.jl +++ /dev/null @@ -1,36 +0,0 @@ - -""" -$(TYPEDSIGNATURES) - -Change the objective sense of optimization. Possible arguments are -`JuMP.MAX_SENSE` and `JuMP.MIN_SENSE`. -""" -set_objective_sense(objective_sense) = - (_, opt_model) -> J.set_objective_sense(opt_model, objective_sense) - -""" -$(TYPEDSIGNATURES) - -Change the JuMP optimizer used to run the optimization. -""" -set_optimizer(optimizer) = (_, opt_model) -> J.set_optimizer(opt_model, optimizer) - -""" -$(TYPEDSIGNATURES) - -Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer -to the JuMP documentation and the documentation of the specific optimizer for -usable keys and values. -""" -set_optimizer_attribute(attribute_key, value) = - (_, opt_model) -> J.set_optimizer_attribute(opt_model, attribute_key, value) - -""" - silence - -Modification that disable all output from the JuMP optimizer (shortcut for -`set_silent` from JuMP). -""" -const silence = (_, opt_model) -> J.set_silent(opt_model) - -export set_objective_sense, set_optimizer, set_optimizer_attribute, silence diff --git a/src/analysis/parsimonious_flux_balance_analysis.jl b/src/analysis/parsimonious_flux_balance_analysis.jl deleted file mode 100644 index e52e412a7..000000000 --- a/src/analysis/parsimonious_flux_balance_analysis.jl +++ /dev/null @@ -1,107 +0,0 @@ -""" -$(TYPEDSIGNATURES) - -Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA -runs two consecutive optimization problems. The first is traditional FBA: -``` -max cᵀx = μ -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -And the second is a quadratic optimization problem: -``` -min Σᵢ xᵢ² -s.t. S x = b - xₗ ≤ x ≤ xᵤ - μ = μ⁰ -``` -Where the optimal solution of the FBA problem, μ⁰, has been added as an -additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M, -Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N, -Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils, -Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data -from evolved E. coli are consistent with computed optimal growth from -genome-scale models. Molecular Systems Biology, 6. 390. doi: -accession:10.1038/msb.2010.47" for more details. - -pFBA gets the model optimum by standard FBA (using -[`flux_balance_analysis`](@ref) with `optimizer` and `modifications`), then -finds a minimal total flux through the model that still satisfies the (slightly -relaxed) optimum. This is done using a quadratic problem optimizer. If the -original optimizer does not support quadratic optimization, it can be changed -using the callback in `qp_modifications`, which are applied after the FBA. See -the documentation of [`flux_balance_analysis`](@ref) for usage examples of -modifications. - -The `relax` parameter, which defaults to `0.999`, relaxes the objective bound of -the original bound to prevent numerical issues. - -Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be -found. - -# Example -``` -model = load_model("e_coli_core.json") -parsimonious_flux_balance_analysis(model, Gurobi.Optimizer) -``` -""" -function parsimonious_flux_balance_analysis( - ctmodel::C.ConstraintTree, - optimizer; - modifications = [], - qp_modifications = [], - relax = 0.999, -) - # Run FBA - opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer) - - for mod in modifications - mod(ctmodel, opt_model) - end - - J.optimize!(opt_model) - - is_solved(opt_model) || return nothing - - # get the objective - Z = J.objective_value(opt_model) - - _ctmodel = ctmodel * :pfbaobjective^squared_sum_objective(ctmodel.fluxes) - _ctmodel.objective.bound = relax == 1.0 ? Z : objective_bounds(relax)(Z) # TODO currently breaks - - opt_model = X.optimization_model( - _ctmodel; - objective = ctmodel.pfbaobjective.value, - optimizer, - sense = X.J.MIN_SENSE, - ) - - # prepare the model for pFBA - for mod in qp_modifications - mod(model, opt_model) - end - - J.optimize!(opt_model) - - C.ValueTree(_ctmodel, J.value.(opt_model[:x])) -end - -""" -$(TYPEDSIGNATURES) - -Variant that takes an [`A.AbstractFBCModel`](@ref) as input. All other arguments are forwarded. -""" -function parsimonious_flux_balance_analysis(model::A.AbstractFBCModel, optimizer; kwargs...) - ctmodel = fbc_model_constraints(model) - parsimonious_flux_balance_analysis(ctmodel, optimizer; kwargs...) -end - -""" -$(TYPEDSIGNATURES) - -Pipe-able variant of [`parsimonious_flux_balance_analysis`](@ref). -""" -parsimonious_flux_balance_analysis(optimizer; modifications = []) = - m -> parsimonious_flux_balance_analysis(m, optimizer; modifications) - -export parsimonious_flux_balance_analysis diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index 7bc9075f5..c27ca5538 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -1,4 +1,17 @@ +""" +$(TYPEDSIGNATURES) + +TODO +""" +squared_sum_objective(x::C.ConstraintTree) = + squared_sum_error_objective(x, Dict(keys(x) .=> 0.0)) + +""" +$(TYPEDSIGNATURES) + +TODO +""" squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = C.Constraint( sum( @@ -6,11 +19,3 @@ squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,F (k, c) in constraints if haskey(target, k) ), ) - -squared_sum_objective(x::C.ConstraintTree) = - squared_sum_error_objective(x, Dict(keys(x) .=> 0.0)) - -objective_bounds(tolerance) = z -> begin - vs = (z * tolerance, z / tolerance) - (minimum(vs), maximum(vs)) -end diff --git a/src/misc/bounds.jl b/src/misc/bounds.jl new file mode 100644 index 000000000..5ecf68cc3 --- /dev/null +++ b/src/misc/bounds.jl @@ -0,0 +1,20 @@ + +""" +$(TYPEDSIGNATURES) + +TODO +""" +absolute_tolerance_bound(tolerance) = x -> begin + bound = (x - tolerance, x + tolerance) + (minimum(bound), maximum(bound)) +end + +""" +$(TYPEDSIGNATURES) + +TODO +""" +relative_tolerance_bound(tolerance) = x -> begin + bound = (x * tolerance, x / tolerance) + (minimum(bound), maximum(bound)) +end diff --git a/src/solver.jl b/src/solver.jl index c1a3b2b7b..ed0255f17 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -47,3 +47,26 @@ is_solved(opt_model::J.Model) = J.termination_status(opt_model) in [J.MOI.OPTIMAL, J.MOI.LOCALLY_SOLVED] export is_solved + +""" +$(TYPEDSIGNATURES) + +Make an JuMP model out of `constraints` using [`optimization_model`](@ref) +(most arguments are forwarded there), then apply the modifications, optimize +the model, and return either `nothing` if the optimization failed, or `output` +substituted with the solved values (`output` defaults to `constraints`. +""" +function optimize_constraints( + constraints::C.ConstraintTree, + args...; + modifications = [], + output = constraints, + kwargs..., +) + om = optimization_model(constraints, args..., kwargs...) + for m in modifications + m(om) + end + J.optimize!(om) + is_solved(om) ? C.ValueTree(output, J.value.(opt_model[:x])) : nothing +end From a1c2edf47aac71bfd6695acf3af9d332cad53f07 Mon Sep 17 00:00:00 2001 From: "St. Elmo Wilken" Date: Sat, 9 Dec 2023 11:34:06 +0100 Subject: [PATCH 14/26] update for CT 0.6 --- Project.toml | 2 +- src/analysis/flux_balance.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 14daa42be..ee0651e56 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" AbstractFBCModels = "0.2" Aqua = "0.7" Clarabel = "0.3" -ConstraintTrees = "0.5" +ConstraintTrees = "0.6" DistributedData = "0.2" DocStringExtensions = "0.8, 0.9" Downloads = "1" diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index f755b45fc..38eeec3dc 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -21,7 +21,7 @@ before the analysis, such as [`set_objective_sense`](@ref), [`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and [`silence`](@ref). -Returns a [`C.ValueTree`](@ref). +Returns a solved tree. # Example ``` @@ -52,7 +52,7 @@ function flux_balance(ctmodel::C.ConstraintTree, optimizer; modifications = []) is_solved(opt_model) || return nothing - C.ValueTree(ctmodel, J.value.(opt_model[:x])) + C.constraint_values(ctmodel, J.value.(opt_model[:x])) end """ From a30eea65fd9b2773d83c35e0a36e593749b61844 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 14:55:20 +0100 Subject: [PATCH 15/26] new functionality in CTs 0.6 --- Project.toml | 3 +-- docs/src/examples/03-qp-problems.jl | 4 ++-- src/analysis/flux_balance.jl | 5 +++-- src/solver.jl | 10 +++++----- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 544f065b4..40fdc022c 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" AbstractFBCModels = "0.2" Aqua = "0.7" Clarabel = "0.3" -ConstraintTrees = "0.5" +ConstraintTrees = "0.6" DistributedData = "0.2" DocStringExtensions = "0.8, 0.9" Downloads = "1" @@ -44,7 +44,6 @@ SBMLFBCModels = "3e8f9d1a-ffc1-486d-82d6-6c7276635980" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" -JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8" [targets] test = ["Aqua", "Clarabel", "Downloads", "GLPK", "JSONFBCModels", "SBMLFBCModels", "SHA", "Test", "Tulip"] diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index ea8592c4a..7109652e7 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -49,7 +49,7 @@ X.J.optimize!(opt_model) # JuMP is called J in COBREXA X.is_solved(opt_model) # check if solved -vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA +vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA @test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds @@ -91,6 +91,6 @@ X.J.optimize!(opt_model) # JuMP is called J in COBREXA X.is_solved(opt_model) # check if solved -vt = X.C.ValueTree(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA +vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA @test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index bb58d8dbf..0b9ad1d21 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -21,7 +21,8 @@ before the analysis, such as [`set_objective_sense`](@ref), [`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and [`silence`](@ref). -Returns a [`C.ValueTree`](@ref). +Returns a tree with the optimization solution of the same shape as the model +defined by [`fbc_model_constraints`](@ref). # Example ``` @@ -33,7 +34,7 @@ function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...) constraints = fbc_model_constraints(model) optimize_constraints( constraints; - objective = constraints.objective, + objective = constraints.objective.value, optimizer, kwargs..., ) diff --git a/src/solver.jl b/src/solver.jl index ed0255f17..cdd496ce2 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -7,8 +7,8 @@ JuMP `Model` created for solving in `optimizer`, with a given optional `objective` and optimization `sense`. """ function optimization_model( - cs::C.ConstraintTree; - objective::Union{Nothing,C.LinearValue,C.QuadraticValue} = nothing, + cs::C.ConstraintTreeElem; + objective::Union{Nothing,C.Value} = nothing, optimizer, sense = J.MAX_SENSE, ) @@ -57,16 +57,16 @@ the model, and return either `nothing` if the optimization failed, or `output` substituted with the solved values (`output` defaults to `constraints`. """ function optimize_constraints( - constraints::C.ConstraintTree, + constraints::C.ConstraintTreeElem, args...; modifications = [], output = constraints, kwargs..., ) - om = optimization_model(constraints, args..., kwargs...) + om = optimization_model(constraints, args...; kwargs...) for m in modifications m(om) end J.optimize!(om) - is_solved(om) ? C.ValueTree(output, J.value.(opt_model[:x])) : nothing + is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing end From f1884dbbb5e3a69cc44ed8e782120a8ce8deeac6 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 16:34:07 +0100 Subject: [PATCH 16/26] isolate errors in tested doc examples --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1bdc67ade..88b05b85b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,7 +21,9 @@ end function run_doc_examples() for dir in filter(endswith(".jl"), readdir("../docs/src/examples", join = true)) - run_test_file(dir) + @testset "docs/$(basename(dir))" begin + run_test_file(dir) + end end end From a985ffa8e4ca0e8a3a6697b5e311a8877af0d056 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 16:35:19 +0100 Subject: [PATCH 17/26] start cleaning out the various names for optimization --- .../examples/02c-constraint-modifications.jl | 23 ++++++++++++++++--- docs/src/examples/03-qp-problems.jl | 17 ++++++-------- src/analysis/flux_balance.jl | 2 +- src/solver.jl | 4 +++- 4 files changed, 31 insertions(+), 15 deletions(-) diff --git a/docs/src/examples/02c-constraint-modifications.jl b/docs/src/examples/02c-constraint-modifications.jl index 584ba25ac..3ff451473 100644 --- a/docs/src/examples/02c-constraint-modifications.jl +++ b/docs/src/examples/02c-constraint-modifications.jl @@ -43,7 +43,12 @@ fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value forced_mixed_fermentation = ctmodel * :fermentation^C.Constraint(fermentation, (10.0, 1000.0)) # new modified model is created -vt = flux_balance(forced_mixed_fermentation, Tulip.Optimizer; modifications = [silence]) +vt = optimized_constraints( + forced_mixed_fermentation, + objective = forced_mixed_fermentation.objective.value, + optimizer = Tulip.Optimizer, + modifications = [silence], +) @test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src @@ -52,13 +57,25 @@ vt = flux_balance(forced_mixed_fermentation, Tulip.Optimizer; modifications = [s ctmodel.fluxes.ATPM.bound = (1000.0, 10000.0) -vt = flux_balance(ctmodel, Tulip.Optimizer; modifications = [silence]) +#TODO explicitly show here how false sharing looks like + +vt = optimized_constraints( + ctmodel, + objective = ctmodel.objective.value, + optimizer = Tulip.Optimizer, + modifications = [silence], +) @test isnothing(vt) #src # Models can also be piped into the analysis functions ctmodel.fluxes.ATPM.bound = (8.39, 10000.0) # revert -vt = ctmodel |> flux_balance(Tulip.Optimizer; modifications = [silence]) +vt = optimized_constraints( + ctmodel, + objective = ctmodel.objective.value, + optimizer = Tulip.Optimizer, + modifications = [silence], +) @test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index 7109652e7..329088f81 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -1,8 +1,8 @@ # # Quandratic objective flux balance analysis type problems -# We will use [`parsimonious_flux_balance_analysis`](@ref) and -# [`minimize_metabolic_adjustment_analysis`](@ref) to find the optimal flux +# We will use [`parsimonious_flux_balance`](@ref) and +# [`minimize_metabolic_adjustment`](@ref) to find the optimal flux # distribution in the *E. coli* "core" model. # If it is not already present, download the model and load the package: @@ -22,11 +22,11 @@ model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model # Use the convenience function to run standard pFBA -vt = X.parsimonious_flux_balance_analysis(model, Clarabel.Optimizer) +vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer) # Or use the piping functionality -model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; modifications = [X.silence]) +model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [X.silence]) @test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src @test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol = QP_TEST_TOLERANCE) #src @@ -57,15 +57,12 @@ vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees ref_sol = Dict("ATPS4r" => 33.0, "CYTBD" => 22.0) -vt = X.minimize_metabolic_adjustment_analysis(model, ref_sol, Gurobi.Optimizer) +vt = X.minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer) # Or use the piping functionality -model |> X.minimize_metabolic_adjustment_analysis( - ref_sol, - Clarabel.Optimizer; - modifications = [X.silence], -) +model |> +X.minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [X.silence]) @test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index 0b9ad1d21..fd3583ebe 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -32,7 +32,7 @@ solution = flux_balance(model, GLPK.optimizer) """ function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...) constraints = fbc_model_constraints(model) - optimize_constraints( + optimized_constraints( constraints; objective = constraints.objective.value, optimizer, diff --git a/src/solver.jl b/src/solver.jl index cdd496ce2..4a6e41aef 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -56,7 +56,7 @@ Make an JuMP model out of `constraints` using [`optimization_model`](@ref) the model, and return either `nothing` if the optimization failed, or `output` substituted with the solved values (`output` defaults to `constraints`. """ -function optimize_constraints( +function optimized_constraints( constraints::C.ConstraintTreeElem, args...; modifications = [], @@ -70,3 +70,5 @@ function optimize_constraints( J.optimize!(om) is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing end + +export optimized_constraints From 2e0ec3c58d6b74eff04fe26de2c79b12ceffdb12 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 16:39:05 +0100 Subject: [PATCH 18/26] kill trailing spaces --- docs/src/examples/03-qp-problems.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index 329088f81..b0f1c5147 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -31,7 +31,7 @@ model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [X.silenc @test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src @test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol = QP_TEST_TOLERANCE) #src -# Alternatively, you can construct your own constraint tree model with +# Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). ctmodel = X.fbc_model_constraints(model) @@ -53,7 +53,7 @@ vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees @test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds -# It is likewise as simple to run MOMA using the convenience functions. +# It is likewise as simple to run MOMA using the convenience functions. ref_sol = Dict("ATPS4r" => 33.0, "CYTBD" => 22.0) @@ -66,7 +66,7 @@ X.minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [X. @test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src -# Alternatively, you can construct your own constraint tree model with +# Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). ctmodel = X.fbc_model_constraints(model) From 48cbc6ed0540d8b75504882c7e0d8a9242c5d425 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 19:55:51 +0100 Subject: [PATCH 19/26] split pFBA into CT and easy-modeling part --- docs/src/examples/03-qp-problems.jl | 2 +- src/analysis/flux_balance.jl | 62 ++++---- src/analysis/parsimonious_flux_balance.jl | 171 ++++++++++++---------- src/builders/objectives.jl | 8 +- src/solver.jl | 25 ---- 5 files changed, 133 insertions(+), 135 deletions(-) diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index b0f1c5147..19bdc56d5 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -22,7 +22,7 @@ model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model # Use the convenience function to run standard pFBA -vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer) +vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [X.silence]) # Or use the piping functionality diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index fd3583ebe..745c99896 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -1,34 +1,40 @@ """ $(TYPEDSIGNATURES) -Run flux balance analysis (FBA) on the `model`, optionally specifying -`modifications` to the problem. Basically, FBA solves this optimization -problem: -``` -max cᵀx -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat -Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more -information. - -The `optimizer` must be set to a `JuMP`-compatible optimizer, such as -`GLPK.Optimizer` or `Tulip.Optimizer`. - -Optionally, you may specify one or more modifications to be applied to the model -before the analysis, such as [`set_objective_sense`](@ref), -[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and -[`silence`](@ref). - -Returns a tree with the optimization solution of the same shape as the model -defined by [`fbc_model_constraints`](@ref). - -# Example -``` -model = load_model("e_coli_core.json") -solution = flux_balance(model, GLPK.optimizer) -``` +Make an JuMP model out of `constraints` using [`optimization_model`](@ref) +(most arguments are forwarded there), then apply the `modifications`, optimize +the model, and return either `nothing` if the optimization failed, or `output` +substituted with the solved values (`output` defaults to `constraints`. + +For a "nice" version for simpler finding of metabolic model optima, use +[`flux_balance`](@ref). +""" +function optimized_constraints( + constraints::C.ConstraintTreeElem, + args...; + modifications = [], + output = constraints, + kwargs..., +) + om = optimization_model(constraints, args...; kwargs...) + for m in modifications + m(om) + end + J.optimize!(om) + is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing +end + +export optimized_constraints + +""" +$(TYPEDSIGNATURES) + +Compute an optimal objective-optimizing solution of the given `model`. + +Most arguments are forwarded to [`optimized_constraints`](@ref). + +Returns a tree with the optimization solution of the same shape as +given by [`fbc_model_constraints`](@ref). """ function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...) constraints = fbc_model_constraints(model) diff --git a/src/analysis/parsimonious_flux_balance.jl b/src/analysis/parsimonious_flux_balance.jl index 26d842a6e..0658b0d73 100644 --- a/src/analysis/parsimonious_flux_balance.jl +++ b/src/analysis/parsimonious_flux_balance.jl @@ -1,92 +1,111 @@ + """ $(TYPEDSIGNATURES) -Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA -runs two consecutive optimization problems. The first is traditional FBA: -``` -max cᵀx = μ -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -And the second is a quadratic optimization problem: -``` -min Σᵢ xᵢ² -s.t. S x = b - xₗ ≤ x ≤ xᵤ - μ = μ⁰ -``` -Where the optimal solution of the FBA problem, μ⁰, has been added as an -additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M, -Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N, -Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils, -Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data -from evolved E. coli are consistent with computed optimal growth from -genome-scale models. Molecular Systems Biology, 6. 390. doi: -accession:10.1038/msb.2010.47" for more details. - -pFBA gets the model optimum by standard FBA (using -[`flux_balance`](@ref) with `optimizer` and `modifications`), then -finds a minimal total flux through the model that still satisfies the (slightly -relaxed) optimum. This is done using a quadratic problem optimizer. If the -original optimizer does not support quadratic optimization, it can be changed -using the callback in `qp_modifications`, which are applied after the FBA. See -the documentation of [`flux_balance`](@ref) for usage examples of -modifications. - -The optimum relaxation sequence can be specified in `relax` parameter, it -defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original -bound. - -Returns an optimized model that contains the pFBA solution (or an unsolved model -if something went wrong). - -# Performance - -This implementation attempts to save time by executing all pFBA steps on a -single instance of the optimization model problem, trading off possible -flexibility. For slightly less performant but much more flexible use, one can -construct parsimonious models directly using -[`with_parsimonious_objective`](@ref). - -# Example -``` -model = load_model("e_coli_core.json") -parsimonious_flux_balance(model, biomass, Gurobi.Optimizer) |> values_vec -``` +Optimize the system of `constraints` to get the optimal `objective` value. Then +try to find a "parsimonious" solution with the same `objective` value, which +optimizes the `parsimonious_objective` (possibly also switching optimization +sense, optimizer, and adding more modifications). + +For efficiency, everything is performed on a single instance of JuMP model. + +A simpler version suitable for direct work with metabolic models is available +in [`parsimonious_flux_balance`](@ref). """ -function parsimonious_flux_balance( - model::C.ConstraintTree, - optimizer; +function parsimonious_optimized_constraints( + constraints::C.ConstraintTreeElem, + args...; + objective::C.Value, modifications = [], - qp_modifications = [], - relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99], + parsimonious_objective::C.Value, + parsimonious_optimizer = nothing, + parsimonious_sense = J.MIN_SENSE, + parsimonious_modifications = [], + tolerances = [absolute_tolerance_bound(0)], + output = constraints, + kwargs..., ) - # Run FBA - opt_model = flux_balance(model, optimizer; modifications) - J.is_solved(opt_model) || return nothing # FBA failed - # get the objective - Z = J.objective_value(opt_model) - original_objective = J.objective_function(opt_model) + # first solve the optimization problem with the original objective + om = optimization_model(constraints, args...; kwargs...) + for m in modifications + m(om) + end + J.optimize!(om) + is_solved(om) || return nothing + + target_objective_value = J.objective_value(om) - # prepare the model for pFBA - for mod in qp_modifications - mod(model, opt_model) + # switch to parsimonizing the solution w.r.t. to the objective value + isnothing(parsimonious_optimizer) || J.set_optimizer(om, parsimonious_optimizer) + for m in parsimonious_modifications + m(om) end - # add the minimization constraint for total flux - v = opt_model[:x] # fluxes - J.@objective(opt_model, Min, sum(dot(v, v))) + J.@objective(om, J.MIN_SENSE, C.substitute(parsimonious_objective, om[:x])) - for rb in relax_bounds - # lb, ub = objective_bounds(rb)(Z) - J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub) + # try all admissible tolerances + for tolerance in tolerances + (lb, ub) = tolerance(target_objective_value) + J.@constraint( + om, + pfba_tolerance_constraint, + lb <= C.substitute(objective, om[:x]) <= ub + ) - J.optimize!(opt_model) - J.is_solved(opt_model) && break + J.optimize!(om) + is_solved(om) && return C.constraint_values(output, J.value.(om[:x])) - J.delete(opt_model, pfba_constraint) - J.unregister(opt_model, :pfba_constraint) + J.delete(om, pfba_tolerance_constraint) + J.unregister(om, :pfba_tolerance_constraint) end + # all tolerances failed + return nothing end + +export parsimonious_optimized_constraints + +""" +$(TYPEDSIGNATURES) + +Compute a parsimonious flux solution for the given `model`. In short, the +objective value of the parsimonious solution should be the same as the one from +[`flux_balance`](@ref), except the squared sum of reaction fluxes is minimized. +If there are multiple possible fluxes that achieve a given objective value, +parsimonious flux thus represents the "minimum energy" one, thus arguably more +realistic. + +Most arguments are forwarded to [`parsimonious_optimized_constraints`](@ref), +with some (objectives) filled in automatically to fit the common processing of +FBC models, and some (`tolerances`) provided with more practical defaults. + +Similarly to the [`flux_balance`](@ref), returns a tree with the optimization +solutions of the shape as given by [`fbc_model_constraints`](@ref). +""" +function parsimonious_flux_balance( + model::A.AbstractFBCModel, + optimizer; + tolerances = relative_tolerance_bound.(1 .- [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]), + kwargs..., +) + constraints = fbc_model_constraints(model) + parsimonious_optimized_constraints( + constraints; + optimizer, + objective = constraints.objective.value, + parsimonious_objective = squared_sum_objective(constraints.fluxes), + tolerances, + kwargs..., + ) +end + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant of [`parsimonious_flux_balance`](@ref). +""" +parsimonious_flux_balance(optimizer; kwargs...) = + model -> parsimonious_flux_balance(model, optimizer; kwargs...) + +export parsimonious_flux_balance diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index c27ca5538..e6f7a8822 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -13,9 +13,7 @@ $(TYPEDSIGNATURES) TODO """ squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = - C.Constraint( - sum( - (C.value(c) - target[k]) * (C.value(c) - target[k]) for - (k, c) in constraints if haskey(target, k) - ), + sum( + (C.squared(C.value(c) - target[k]) for (k, c) in constraints if haskey(target, k)), + init = zero(C.LinearValue), ) diff --git a/src/solver.jl b/src/solver.jl index 4a6e41aef..b0d24fd52 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -47,28 +47,3 @@ is_solved(opt_model::J.Model) = J.termination_status(opt_model) in [J.MOI.OPTIMAL, J.MOI.LOCALLY_SOLVED] export is_solved - -""" -$(TYPEDSIGNATURES) - -Make an JuMP model out of `constraints` using [`optimization_model`](@ref) -(most arguments are forwarded there), then apply the modifications, optimize -the model, and return either `nothing` if the optimization failed, or `output` -substituted with the solved values (`output` defaults to `constraints`. -""" -function optimized_constraints( - constraints::C.ConstraintTreeElem, - args...; - modifications = [], - output = constraints, - kwargs..., -) - om = optimization_model(constraints, args...; kwargs...) - for m in modifications - m(om) - end - J.optimize!(om) - is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing -end - -export optimized_constraints From a2c121e0e4dd546d8bbfaa671c291efb6305c5e3 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 23:00:53 +0100 Subject: [PATCH 20/26] fixes --- Project.toml | 6 +- ...ems.jl => 03-parsimonious-flux-balance.jl} | 56 ++++++++++--------- ...04-minimization-of-metabolic-adjustment.jl | 23 ++++++++ src/analysis/parsimonious_flux_balance.jl | 2 +- 4 files changed, 58 insertions(+), 29 deletions(-) rename docs/src/examples/{03-qp-problems.jl => 03-parsimonious-flux-balance.jl} (56%) create mode 100644 docs/src/examples/04-minimization-of-metabolic-adjustment.jl diff --git a/Project.toml b/Project.toml index 40fdc022c..1ec193418 100644 --- a/Project.toml +++ b/Project.toml @@ -17,10 +17,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" [compat] -AbstractFBCModels = "0.2" +AbstractFBCModels = "0.2.2" Aqua = "0.7" -Clarabel = "0.3" -ConstraintTrees = "0.6" +Clarabel = "0.6" +ConstraintTrees = "0.6.1" DistributedData = "0.2" DocStringExtensions = "0.8, 0.9" Downloads = "1" diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-parsimonious-flux-balance.jl similarity index 56% rename from docs/src/examples/03-qp-problems.jl rename to docs/src/examples/03-parsimonious-flux-balance.jl index 19bdc56d5..924141f41 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-parsimonious-flux-balance.jl @@ -1,9 +1,11 @@ -# # Quandratic objective flux balance analysis type problems +# # Parsimonious flux balance analysis # We will use [`parsimonious_flux_balance`](@ref) and # [`minimize_metabolic_adjustment`](@ref) to find the optimal flux # distribution in the *E. coli* "core" model. +# +# TODO pFBA citation # If it is not already present, download the model and load the package: import Downloads: download @@ -13,43 +15,45 @@ import Downloads: download # next, load the necessary packages -import COBREXA as X -import AbstractFBCModels as A # for the accessors -import JSONFBCModels as J # for the model type +using COBREXA + +import JSONFBCModels import Clarabel # can solve QPs -model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model +model = load_model("e_coli_core.json") # load the model -# Use the convenience function to run standard pFBA +# Use the convenience function to run standard pFBA on -vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [X.silence]) +vt = parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [silence]) # Or use the piping functionality -model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [X.silence]) +model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [silence]) @test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src -@test isapprox(sum(x^2 for x in values(vt.fluxes)), 11414.2119; atol = QP_TEST_TOLERANCE) #src +@test sum(x^2 for x in values(vt.fluxes)) < 15000 #src + +#= # Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). -ctmodel = X.fbc_model_constraints(model) -ctmodel *= :l2objective^X.squared_sum_objective(ctmodel.fluxes) +ctmodel = fbc_model_constraints(model) +ctmodel *= :l2objective^squared_sum_objective(ctmodel.fluxes) ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks -opt_model = X.optimization_model( +opt_model = optimization_model( ctmodel; objective = ctmodel.:l2objective.value, optimizer = Clarabel.Optimizer, - sense = X.J.MIN_SENSE, + sense = J.MIN_SENSE, ) -X.J.optimize!(opt_model) # JuMP is called J in COBREXA +J.optimize!(opt_model) # JuMP is called J in COBREXA -X.is_solved(opt_model) # check if solved +is_solved(opt_model) # check if solved -vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA +vt = C.constraint_values(ctmodel, J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA @test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds @@ -57,37 +61,39 @@ vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees ref_sol = Dict("ATPS4r" => 33.0, "CYTBD" => 22.0) -vt = X.minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer) +vt = minimize_metabolic_adjustment(model, ref_sol, Gurobi.Optimizer) # Or use the piping functionality model |> -X.minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [X.silence]) +minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [silence]) @test isapprox(vt.:momaobjective, 0.81580806; atol = TEST_TOLERANCE) #src # Alternatively, you can construct your own constraint tree model with # the quadratic objective (this approach is much more flexible). -ctmodel = X.fbc_model_constraints(model) +ctmodel = fbc_model_constraints(model) ctmodel *= - :minoxphospho^X.squared_sum_error_objective( + :minoxphospho^squared_sum_error_objective( ctmodel.fluxes, Dict(:ATPS4r => 33.0, :CYTBD => 22.0), ) ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks -opt_model = X.optimization_model( +opt_model = optimization_model( ctmodel; objective = ctmodel.minoxphospho.value, optimizer = Clarabel.Optimizer, - sense = X.J.MIN_SENSE, + sense = J.MIN_SENSE, ) -X.J.optimize!(opt_model) # JuMP is called J in COBREXA +J.optimize!(opt_model) # JuMP is called J in COBREXA -X.is_solved(opt_model) # check if solved +is_solved(opt_model) # check if solved -vt = X.C.constraint_values(ctmodel, X.J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA +vt = C.constraint_values(ctmodel, J.value.(opt_model[:x])) # ConstraintTrees.jl is called C in COBREXA @test isapprox(vt.l2objective, ?; atol = QP_TEST_TOLERANCE) #src # TODO will break until mutable bounds + +=# diff --git a/docs/src/examples/04-minimization-of-metabolic-adjustment.jl b/docs/src/examples/04-minimization-of-metabolic-adjustment.jl new file mode 100644 index 000000000..0cfdd4050 --- /dev/null +++ b/docs/src/examples/04-minimization-of-metabolic-adjustment.jl @@ -0,0 +1,23 @@ + +# # Minimization of metabolic adjustment + +# TODO MOMA citation + +import Downloads: download + +!isfile("e_coli_core.json") && + download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") + +using COBREXA +import AbstractFBCModels.CanonicalModel as CM +import JSONFBCModels +import Clarabel + +# TODO this might do the convert immediately as with the old cobrexa... +# probably better have an actual output-type argument tho rather than force the +# guessing. +model = convert(CM.Model, load_model("e_coli_core.json")) + +reference_fluxes = parsimonious_flux_balance(model, Clarabel.Optimizer).fluxes + +# TODO MOMA from here diff --git a/src/analysis/parsimonious_flux_balance.jl b/src/analysis/parsimonious_flux_balance.jl index 0658b0d73..c231ce5cd 100644 --- a/src/analysis/parsimonious_flux_balance.jl +++ b/src/analysis/parsimonious_flux_balance.jl @@ -27,7 +27,7 @@ function parsimonious_optimized_constraints( ) # first solve the optimization problem with the original objective - om = optimization_model(constraints, args...; kwargs...) + om = optimization_model(constraints, args...; objective, kwargs...) for m in modifications m(om) end From 53cc6b75fd7e8b01e9e97fec99746dcc4b3f2f40 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 12:14:25 +0100 Subject: [PATCH 21/26] clean up moma --- .../examples/03-parsimonious-flux-balance.jl | 4 +- src/analysis/flux_balance.jl | 5 +- src/analysis/minimal_metabolic_adjustment.jl | 101 ++++++++++-------- src/analysis/modifications.jl | 4 +- src/analysis/parsimonious_flux_balance.jl | 13 +-- src/builders/core.jl | 5 + src/solver.jl | 29 ++++- 7 files changed, 101 insertions(+), 60 deletions(-) diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl index 924141f41..2e1a9c813 100644 --- a/docs/src/examples/03-parsimonious-flux-balance.jl +++ b/docs/src/examples/03-parsimonious-flux-balance.jl @@ -46,7 +46,7 @@ opt_model = optimization_model( ctmodel; objective = ctmodel.:l2objective.value, optimizer = Clarabel.Optimizer, - sense = J.MIN_SENSE, + sense = Minimal, ) J.optimize!(opt_model) # JuMP is called J in COBREXA @@ -85,7 +85,7 @@ opt_model = optimization_model( ctmodel; objective = ctmodel.minoxphospho.value, optimizer = Clarabel.Optimizer, - sense = J.MIN_SENSE, + sense = Minimal, ) J.optimize!(opt_model) # JuMP is called J in COBREXA diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index 745c99896..5f9d06182 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -10,13 +10,12 @@ For a "nice" version for simpler finding of metabolic model optima, use [`flux_balance`](@ref). """ function optimized_constraints( - constraints::C.ConstraintTreeElem, - args...; + constraints::C.ConstraintTreeElem; modifications = [], output = constraints, kwargs..., ) - om = optimization_model(constraints, args...; kwargs...) + om = optimization_model(constraints; kwargs...) for m in modifications m(om) end diff --git a/src/analysis/minimal_metabolic_adjustment.jl b/src/analysis/minimal_metabolic_adjustment.jl index d586a3aa8..d7dcc29c5 100644 --- a/src/analysis/minimal_metabolic_adjustment.jl +++ b/src/analysis/minimal_metabolic_adjustment.jl @@ -2,67 +2,76 @@ """ $(TYPEDSIGNATURES) -Run minimization of metabolic adjustment (MOMA) on `model` with respect to -`reference_solution`, which is a dictionary of fluxes. MOMA finds the shortest -Euclidian distance between `reference_solution` and `model` with `modifications`: -``` -min Σᵢ (xᵢ - flux_refᵢ)² -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -Because the problem has a quadratic objective, a QP solver is required. See -"Daniel, Vitkup & Church, Analysis of Optimality in Natural and Perturbed -Metabolic Networks, Proceedings of the National Academy of Sciences, 2002" for -more details. +Find a feasible solution of the "minimal metabolic adjustment analysis" (MOMA) +for the `model`, which is the "closest" feasible solution to the given +`reference_fluxes`, in the sense of squared-sum error distance. The minimized +squared distance (the objective) is present in the result tree as +`minimal_adjustment_objective`. -Returns a [`C.ValueTree`](@ref), or `nothing` if the solution could not be -found. +This is often used for models with smaller feasible region than the reference +models (typically handicapped by a knockout, nutritional deficiency or a +similar perturbation). MOMA solution then gives an expectable "easiest" +adjustment of the organism towards a somewhat working state. -# Example -``` -model = load_model("e_coli_core.json") +Reference fluxes that do not exist in the model are ignored (internally, the +objective is constructed via [`squared_sum_error_objective`](@ref)). -``` +Additional parameters are forwarded to [`optimized_constraints`](@ref). """ function minimal_metabolic_adjustment( - constraints::C.ConstraintTree, - reference_solution::Dict{String,Float64}, + model::A.AbstractFBCModel, + reference_fluxes::Dict{Symbol,Float64}, optimizer; - modifications = [], + kwargs..., ) - _constraints = - constraints * - :momaobjective^squared_sum_error_objective( - constraints.fluxes, - Dict(Symbol(k) => float(v) for (k, v) in reference_solution), - ) - - opt_model = optimization_model( - _constraints; - objective = _constraints.momaobjective.value, + constraints = fbc_model_constraints(model) + objective = squared_sum_error_objective(constraints.fluxes, reference_fluxes) + optimized_constraints( + constraints * :minimal_adjustment_objective^C.Constraint(objective); optimizer, - sense = J.MIN_SENSE, + objective, + sense = Minimal, + kwargs..., ) - - for mod in modifications - mod(constraints, opt_model) - end - - J.optimize!(opt_model) - - is_solved(opt_model) || return nothing - - C.ValueTree(_constraints, J.value.(opt_model[:x])) end """ $(TYPEDSIGNATURES) -Variant that takes an [`A.AbstractFBCModel`](@ref) as input. All other arguments are forwarded. +A slightly easier-to-use version of [`minimal_metabolic_adjustment`](@ref) that +computes the reference flux as the optimal solution of the +[`reference_model`](@ref). The reference flux is calculated using +`reference_optimizer` and `reference_modifications`, which default to the +`optimizer` and `modifications`. + +Leftover arguments are passed to the overload of +[`minimal_metabolic_adjustment`](@ref) that accepts the reference flux +dictionary. """ -function minimal_metabolic_adjustment(model::A.AbstractFBCModel, args...; kwargs...) - constraints = fbc_model_constraints(model) - minimal_metabolic_adjustment(constraints, args...; kwargs...) +function minimal_metabolic_adjustment( + model::A.AbstractFBCModel, + reference_model::A.AbstractFBCModel, + optimizer; + reference_optimizer = optimizer, + modifications = [], + reference_modifications = modifications, + kwargs..., +) + reference_constraints = fbc_model_constraints(reference_model) + reference_fluxes = optimized_constraints( + reference_constraints; + optimizer = reference_optimizer, + modifications = reference_modifications, + output = reference_constraints.fluxes, + ) + isnothing(reference_fluxes) && return nothing + minimal_metabolic_adjustment( + model, + reference_fluxes, + optimizer; + modifications, + kwargs..., + ) end export minimal_metabolic_adjustment diff --git a/src/analysis/modifications.jl b/src/analysis/modifications.jl index 585077480..dccec943e 100644 --- a/src/analysis/modifications.jl +++ b/src/analysis/modifications.jl @@ -4,8 +4,8 @@ """ $(TYPEDSIGNATURES) -Change the objective sense of optimization. Possible arguments are -`JuMP.MAX_SENSE` and `JuMP.MIN_SENSE`. +Change the objective sense of optimization. Accepted arguments include +[`Minimal`](@ref), [`Maximal`](@ref), and [`Feasible`](@ref). """ set_objective_sense(objective_sense) = opt_model -> J.set_objective_sense(opt_model, objective_sense) diff --git a/src/analysis/parsimonious_flux_balance.jl b/src/analysis/parsimonious_flux_balance.jl index c231ce5cd..114cb4652 100644 --- a/src/analysis/parsimonious_flux_balance.jl +++ b/src/analysis/parsimonious_flux_balance.jl @@ -13,8 +13,7 @@ A simpler version suitable for direct work with metabolic models is available in [`parsimonious_flux_balance`](@ref). """ function parsimonious_optimized_constraints( - constraints::C.ConstraintTreeElem, - args...; + constraints::C.ConstraintTreeElem; objective::C.Value, modifications = [], parsimonious_objective::C.Value, @@ -27,7 +26,7 @@ function parsimonious_optimized_constraints( ) # first solve the optimization problem with the original objective - om = optimization_model(constraints, args...; objective, kwargs...) + om = optimization_model(constraints; objective, kwargs...) for m in modifications m(om) end @@ -74,7 +73,8 @@ objective value of the parsimonious solution should be the same as the one from [`flux_balance`](@ref), except the squared sum of reaction fluxes is minimized. If there are multiple possible fluxes that achieve a given objective value, parsimonious flux thus represents the "minimum energy" one, thus arguably more -realistic. +realistic. The optimized squared distance is present in the result as +`parsimonious_objective`. Most arguments are forwarded to [`parsimonious_optimized_constraints`](@ref), with some (objectives) filled in automatically to fit the common processing of @@ -90,11 +90,12 @@ function parsimonious_flux_balance( kwargs..., ) constraints = fbc_model_constraints(model) + parsimonious_objective = squared_sum_objective(constraints.fluxes) parsimonious_optimized_constraints( - constraints; + constraints * :parsimonious_objective^C.Constraint(parsimonious_objective); optimizer, objective = constraints.objective.value, - parsimonious_objective = squared_sum_objective(constraints.fluxes), + parsimonious_objective, tolerances, kwargs..., ) diff --git a/src/builders/core.jl b/src/builders/core.jl index 97cf6d62a..0d0624113 100644 --- a/src/builders/core.jl +++ b/src/builders/core.jl @@ -6,6 +6,11 @@ $(TYPEDSIGNATURES) A constraint tree that models the content of the given instance of `AbstractFBCModel`. + +The constructed tree contains subtrees `fluxes` (with the reaction-defining +"variables") and `flux_stoichiometry` (with the metabolite-balance-defining +constraints), and a single constraint `objective` thad describes the objective +function of the model. """ function fbc_model_constraints(model::A.AbstractFBCModel) rxns = Symbol.(A.reactions(model)) diff --git a/src/solver.jl b/src/solver.jl index b0d24fd52..2c48c0ccd 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -10,7 +10,7 @@ function optimization_model( cs::C.ConstraintTreeElem; objective::Union{Nothing,C.Value} = nothing, optimizer, - sense = J.MAX_SENSE, + sense = Maximal, ) model = J.Model(optimizer) @@ -47,3 +47,30 @@ is_solved(opt_model::J.Model) = J.termination_status(opt_model) in [J.MOI.OPTIMAL, J.MOI.LOCALLY_SOLVED] export is_solved + +""" + Minimal + +Objective sense for finding the minimal value of the objective. + +Same as `JuMP.MIN_SENSE`. +""" +const Minimal = J.MIN_SENSE + +""" + Maximal + +Objective sense for finding the maximal value of the objective. + +Same as `JuMP.MAX_SENSE`. +""" +const Maximal = J.MAX_SENSE + +""" + Maximal + +Objective sense for finding the any feasible value of the objective. + +Same as `JuMP.FEASIBILITY_SENSE`. +""" +const Feasible = J.FEASIBILITY_SENSE From c6b8af470bb40dbe766b173c2b65877789622166 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 14:29:40 +0100 Subject: [PATCH 22/26] some work on the docs --- docs/src/examples/02-flux-balance-analysis.jl | 4 + docs/src/examples/02a-optimizer-parameters.jl | 5 +- docs/src/examples/02b-model-modifications.jl | 130 +++++++++++++++++- 3 files changed, 130 insertions(+), 9 deletions(-) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index 64642dbbe..e7ceb8fcc 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -47,3 +47,7 @@ solution.fluxes.PFK # ...or make a "table" of all fluxes through all reactions: collect(solution.fluxes) + +# ## Advanced: Finding flux balance via the low-level interface + +# TODO ConstraintTrees (maybe put this into a separate example?) diff --git a/docs/src/examples/02a-optimizer-parameters.jl b/docs/src/examples/02a-optimizer-parameters.jl index 766d6c46a..2dd2705c0 100644 --- a/docs/src/examples/02a-optimizer-parameters.jl +++ b/docs/src/examples/02a-optimizer-parameters.jl @@ -39,12 +39,13 @@ solution = flux_balance( # To see some of the effects of the configuration changes, you may e.g. # deliberately cripple the optimizer's possibilities to a few iterations, which -# will cause it to fail and return no solution: +# will cause it to fail, return no solution, and verbosely describe what +# happened: solution = flux_balance( model, Tulip.Optimizer; - modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 2)], + modifications = [set_optimizer_attribute("IPM_IterationsLimit", 2)], ) println(solution) diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index 22307fe7c..c401a0944 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -8,13 +8,13 @@ # # With COBREXA, there are 2 different approaches that one can take: # 1. We can change the model structure and use the changed metabolic model. -# This is better for doing simple and small but systematic modifications, such -# as removing metabolites, adding reactions, etc. +# This is better for doing simple and small but systematic modifications, +# such as removing metabolites, adding reactions, etc. # 2. We can intercept the pipeline that converts the metabolic model to -# constraints and then to the optimizer representation, and make small -# modifications along that way. This is better for various technical model -# adjustments, such as using combined objectives or adding reaction-coupling -# constraints. +# constraints and then to the optimizer representation, and make small +# modifications along that way. This is better for various technical model +# adjustments, such as using combined objectives or adding reaction-coupling +# constraints. # # Here we demonstrate the first, "modelling" approach. The main advantage of # that approach is that the modified model is still a FBC model, and you can @@ -25,4 +25,120 @@ # constraints](02c-constraint-modifications.md) for a closer look on how to # modify even such complex constructions. # -# TODO here. :) +# ## Getting the base model + +using COBREXA + +download_model( + "http://bigg.ucsd.edu/static/models/e_coli_core.json", + "e_coli_core.json", + "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8", +) + +import JSONFBCModels + +# For applying the modifications, we will use the canonical model as exported +# from package `AbstractFBCModels`. There are other possibilities, but the +# canonical one is easiest to use for common purposes. + +import AbstractFBCModels.CanonicalModel as CM + +# We can now load the model: + +model = convert(CM.Model, load_model("e_coli_core.json")) + +# The canonical model is quite easy to work with, made basically of the most +# accessible Julia structures possible. For example, you can observe a reaction +# as such: + +model.reactions["PFK"] + +# + +model.reactions["CS"].stoichiometry + +# ## Running FBA on modified models +# +# Since the canonical model is completely mutable, you can change it in any way you like and feed it directly into [`flux_balance`](@ref). Let's first find a "original" solution, so that we have a base solution for comparing: + +import GLPK + +base_solution = flux_balance(model, GLPK.Optimizer) +base_solution.objective + +# Now, for example, we can limit the intake of glucose by the model: + +model.reactions["EX_glc__D_e"] + +# Since the original intake limit is 10 units, let's try limiting that to 5: + +model.reactions["EX_glc__D_e"].lower_bound = -5.0 + +# ...and solve the modified model: +# +low_glucose_solution = flux_balance(model, GLPK.Optimizer) +low_glucose_solution.objective + +@test isapprox(low_glucose_solution.objective, 0.41559777, atol=TEST_TOLERANCE) #src + +# ## Preventing reference-based sharing problems with `deepcopy` +# +# People often want to try different perturbations with a single base model. It +# would therefore look feasible to save retain the "unmodified" model in a +# single variable, and make copies of that with the modifications applied. +# Let's observe what happens: + +base_model = convert(CM.Model, load_model("e_coli_core.json")) # load the base + +modified_model = base_model # copy for modification + +modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 # modify the glucose intake limit + +# Surprisingly, the base model got modified too! + +base_model.reactions["EX_glc__D_e"] + +# This is because Julia uses reference-based sharing whenever anything mutable +# is copied using the `=` operator. While this is extremely useful in many +# scenarios for data processing efficiency and computational speed, it +# unfortunately breaks this simple use-case. +# +# To fix the situation, you should always ensure to make an actual copy of the +# model data by either carefully copying the changed parts with `copy()`, or +# simply by copying the whole model structure as is with `deepcopy()`. Let's +# try again: + +base_model = convert(CM.Model, load_model("e_coli_core.json")) +modified_model = deepcopy(base_model) # this forces an actual copy of the data +modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 + +# With `deepcopy`, the result works as intended: + +(modified_model.reactions["EX_glc__D_e"].lower_bound, base_model.reactions["EX_glc__D_e"].lower_bound) + +@test modified_model.reactions["EX_glc__D_e"].lower_bound != base_model.reactions["EX_glc__D_e"].lower_bound #src + +#md # !!! danger "Avoid overwriting base models when using in-place modifications" +#md # Whenever you are changing a copy of the model, make sure that you are not changing it by a reference. Always use some copy mechanism such as `copy` or `deepcopy` to prevent the default reference-based sharing. + +# ## Observing the differences +# +# We already have a `base_solution` and `low_glucose_solution` from above. What +# is the easiest way to see what has changed? We can quite easily compute +# squared distance between all dictionary entries using Julia function for +# merging dictionaries (called `mergewith`). + +# With that, we can extract the plain difference in fluxes: +flux_differences = mergewith(-, base_solution.fluxes, low_glucose_solution.fluxes) + +# ...and see what were the biggest directional differences: +sort(flux_differences, by = last) + +# ...or compute the squared distance, to see the "absolute" changes: +flux_changes = mergewith((x,y) -> (x-y)^2, base_solution.fluxes, low_glucose_solution.fluxes) + +# ...and again see what changed most: +sort(flux_changes, by = last) + +#md # !!! tip "For realistic comparisons always find an unique flux solution" +#md # Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. [parsimonious FBA](03-parsimonious-flux-balance) to obtain uniquely determined and safely comparable flux solutions. From a39106146eb1b029f44efffd4a7fd56f4e652696 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 14:31:55 +0100 Subject: [PATCH 23/26] format --- docs/src/examples/02b-model-modifications.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index c401a0944..cc49a4b8a 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -79,7 +79,7 @@ model.reactions["EX_glc__D_e"].lower_bound = -5.0 low_glucose_solution = flux_balance(model, GLPK.Optimizer) low_glucose_solution.objective -@test isapprox(low_glucose_solution.objective, 0.41559777, atol=TEST_TOLERANCE) #src +@test isapprox(low_glucose_solution.objective, 0.41559777, atol = TEST_TOLERANCE) #src # ## Preventing reference-based sharing problems with `deepcopy` # @@ -114,9 +114,13 @@ modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 # With `deepcopy`, the result works as intended: -(modified_model.reactions["EX_glc__D_e"].lower_bound, base_model.reactions["EX_glc__D_e"].lower_bound) +( + modified_model.reactions["EX_glc__D_e"].lower_bound, + base_model.reactions["EX_glc__D_e"].lower_bound, +) -@test modified_model.reactions["EX_glc__D_e"].lower_bound != base_model.reactions["EX_glc__D_e"].lower_bound #src +@test modified_model.reactions["EX_glc__D_e"].lower_bound != + base_model.reactions["EX_glc__D_e"].lower_bound #src #md # !!! danger "Avoid overwriting base models when using in-place modifications" #md # Whenever you are changing a copy of the model, make sure that you are not changing it by a reference. Always use some copy mechanism such as `copy` or `deepcopy` to prevent the default reference-based sharing. @@ -135,7 +139,8 @@ flux_differences = mergewith(-, base_solution.fluxes, low_glucose_solution.fluxe sort(flux_differences, by = last) # ...or compute the squared distance, to see the "absolute" changes: -flux_changes = mergewith((x,y) -> (x-y)^2, base_solution.fluxes, low_glucose_solution.fluxes) +flux_changes = + mergewith((x, y) -> (x - y)^2, base_solution.fluxes, low_glucose_solution.fluxes) # ...and again see what changed most: sort(flux_changes, by = last) From be71bf90506dafadfb91f52e71be3e3e2c6b8121 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 15:10:59 +0100 Subject: [PATCH 24/26] some documentation fixes --- docs/src/examples/02-flux-balance-analysis.jl | 2 ++ docs/src/examples/02b-model-modifications.jl | 6 +++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index e7ceb8fcc..5e55cea23 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -5,6 +5,8 @@ # find an optimal flux in the *E. coli* "core" model. We will need the model, # which we can download using [`download_model`](@ref): +using COBREXA + download_model( "http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json", diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index cc49a4b8a..15a385b2e 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -119,7 +119,7 @@ modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 base_model.reactions["EX_glc__D_e"].lower_bound, ) -@test modified_model.reactions["EX_glc__D_e"].lower_bound != +@test modified_model.reactions["EX_glc__D_e"].lower_bound != #src base_model.reactions["EX_glc__D_e"].lower_bound #src #md # !!! danger "Avoid overwriting base models when using in-place modifications" @@ -136,14 +136,14 @@ modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 flux_differences = mergewith(-, base_solution.fluxes, low_glucose_solution.fluxes) # ...and see what were the biggest directional differences: -sort(flux_differences, by = last) +sort(collect(flux_differences), by = last) # ...or compute the squared distance, to see the "absolute" changes: flux_changes = mergewith((x, y) -> (x - y)^2, base_solution.fluxes, low_glucose_solution.fluxes) # ...and again see what changed most: -sort(flux_changes, by = last) +sort(collect(flux_changes), by = last) #md # !!! tip "For realistic comparisons always find an unique flux solution" #md # Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. [parsimonious FBA](03-parsimonious-flux-balance) to obtain uniquely determined and safely comparable flux solutions. From ce818557e1ec150e3af4dda12a69f611f2d2d95f Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 15:23:28 +0100 Subject: [PATCH 25/26] more small doc fixes --- docs/make.jl | 2 +- docs/src/examples/02b-model-modifications.jl | 2 +- docs/src/reference.md | 7 +++++++ src/COBREXA.jl | 18 +++++++++++++++++- 4 files changed, 26 insertions(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index f488b18ed..0c1d35ef0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -29,7 +29,6 @@ find_mds(path) = filter(x -> endswith(x, ".md"), readdir(joinpath(@__DIR__, "src", path))), ) -#TODO migrate this to Documenter-1, and make all checks strict # build the docs makedocs( modules = [COBREXA], @@ -42,6 +41,7 @@ makedocs( ), authors = "The developers of COBREXA.jl", linkcheck = !("skiplinks" in ARGS), + warnonly = true, # TODO: remove later pages = [ "Home" => "index.md", "Examples" => [ diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index 15a385b2e..36ac943eb 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -146,4 +146,4 @@ flux_changes = sort(collect(flux_changes), by = last) #md # !!! tip "For realistic comparisons always find an unique flux solution" -#md # Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. [parsimonious FBA](03-parsimonious-flux-balance) to obtain uniquely determined and safely comparable flux solutions. +#md # Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. [parsimonious FBA](03-parsimonious-flux-balance.md) to obtain uniquely determined and safely comparable flux solutions. diff --git a/docs/src/reference.md b/docs/src/reference.md index 51abbb38c..d791bd478 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -42,6 +42,13 @@ Modules = [COBREXA] Pages = ["src/builders/objectives.jl"] ``` +### Bounds&tolerances helpers + +```@autodocs +Modules = [COBREXA] +Pages = ["src/misc/bounds.jl"] +``` + ## Analysis functions ```@autodocs diff --git a/src/COBREXA.jl b/src/COBREXA.jl index b03902860..f6efb9c24 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -1,5 +1,21 @@ """ -$(README) + module COBREXA + +COnstraint Based Reconstruction and EXascale Analysis. COBREXA provides +functions for construction, modification, simulation and analysis of +constraint-based metabolic models that follows the COBRA methodology. + +COBREXA is built as a front-end for the combination of `AbstractFBCModels.jl` +(provides the model I/O), `ConstraintTrees.jl` (provides the constraint system +organization), `Distributed.jl` (provides HPC execution capability), and +`JuMP.jl` (provides the solvers). + +See the online documentation for a complete description of functionality aided +by copy-pastable examples. + +To start quickly, load your favorite JuMP-compatible solver, use +[`load_model`](@ref) to read a metabolic model from the disk, and solve it with +[`flux_balance`](@ref). """ module COBREXA From 39cb1e925fc72feb84bfeb321de6d4c55b165b4a Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sun, 10 Dec 2023 15:31:24 +0100 Subject: [PATCH 26/26] wording --- docs/src/examples/02b-model-modifications.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index 36ac943eb..0ee5942d0 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -145,5 +145,5 @@ flux_changes = # ...and again see what changed most: sort(collect(flux_changes), by = last) -#md # !!! tip "For realistic comparisons always find an unique flux solution" +#md # !!! tip "For realistic comparisons always use a uniquely defined flux solution" #md # Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. [parsimonious FBA](03-parsimonious-flux-balance.md) to obtain uniquely determined and safely comparable flux solutions.