Merge pull request #608 from LCSB-BioCore/develop

develop -> master for v1.3.0

Project.toml

@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "1.2.3"
+version = "1.3.0"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -21,9 +21,9 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-DistributedData = "0.1.4"
+DistributedData = "0.1.4, 0.2"
 JSON = "0.21"
-JuMP = "0.21.0, 0.22.0"
+JuMP = "0.21.0, 0.22.0, 0.23, 1"
 MAT = "0.10"
 MacroTools = "0.5.6"
 OSQP = "0.6"

README.md

@@ -128,20 +128,20 @@ some reactions were disabled independently:
 
 ```julia
 # convert to a model type that is efficient to modify
-m = convert(StandardModel, m)
+m = convert(StandardModel, model)
 
 # find the model objective value if oxygen or carbon dioxide transports are disabled
 screen(m, # the base model
  variants=[ # this specifies how to generate the desired model variants
  [], # one with no modifications, i.e. the base case
- [with_changed_bound("O2t", lower=0.0, upper=0.0)], # disable oxygen
- [with_changed_bound("CO2t", lower=0.0, upper=0.0)], # disable CO2
- [with_changed_bound("O2t", lower=0.0, upper=0.0),
- with_changed_bound("CO2t", lower=0.0, upper=0.0)], # disable both
+ [with_changed_bound("R_O2t", lower=0.0, upper=0.0)], # disable oxygen
+ [with_changed_bound("R_CO2t", lower=0.0, upper=0.0)], # disable CO2
+ [with_changed_bound("R_O2t", lower=0.0, upper=0.0),
+ with_changed_bound("R_CO2t", lower=0.0, upper=0.0)], # disable both
  ],
  # this specifies what to do with the model variants (received as the argument `x`)
  analysis = x ->
- flux_balance_analysis_dict(x, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+ flux_balance_analysis_dict(x, Tulip.Optimizer)["R_BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 You should receive a result showing that missing oxygen transport makes the

docs/src/advanced.md

@@ -5,4 +5,3 @@
 Pages = joinpath.("advanced", filter(x -> endswith(x, ".md"), readdir("advanced")))
 Depth = 2
 ```
-

docs/src/tutorials.md

@@ -5,4 +5,3 @@
 Pages = joinpath.("tutorials", filter(x -> endswith(x, ".md"), readdir("tutorials")))
 Depth = 2
 ```
-

src/COBREXA.jl

@@ -18,9 +18,7 @@ import Base: findfirst, getindex, show
 import Pkg
 import SBML # conflict with Reaction struct name
 
-
 include("banner.jl")
-_print_banner()
 
 # autoloading
 const _inc(path...) = include(joinpath(path...))
@@ -34,6 +32,7 @@ _inc_all.(
  joinpath("base", "logging"),
  joinpath("base", "macros"),
  joinpath("base", "types"),
+ joinpath("base", "types", "wrappers"),
  "base",
  "io",
  joinpath("io", "show"),

src/analysis/flux_balance_analysis.jl

@@ -1,5 +1,5 @@
 """
- flux_balance_analysis_vec(args...)::Maybe{Vector{Float64}}
+ flux_balance_analysis_vec(model::MetabolicModel, args...)::Maybe{Vector{Float64}}
 
 A variant of FBA that returns a vector of fluxes in the same order as reactions
 of the model, if the solution is found.
@@ -9,8 +9,12 @@ Arguments are passed to [`flux_balance_analysis`](@ref).
 This function is kept for backwards compatibility, use [`flux_vector`](@ref)
 instead.
 """
-flux_balance_analysis_vec(args...; kwargs...)::Maybe{Vector{Float64}} =
- flux_vector(flux_balance_analysis(args...; kwargs...))
+flux_balance_analysis_vec(
+ model::MetabolicModel,
+ args...;
+ kwargs...,
+)::Maybe{Vector{Float64}} =
+ flux_vector(model, flux_balance_analysis(model, args...; kwargs...))
 
 """
  flux_balance_analysis_dict(model::MetabolicModel, args...)::Maybe{Dict{String, Float64}}
@@ -66,7 +70,6 @@ biomass_reaction_id = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
 modified_solution = flux_balance_analysis(model, GLPK.optimizer;
  modifications=[change_objective(biomass_reaction_id)])
 ```
-
 """
 function flux_balance_analysis(
  model::M,

src/analysis/flux_variability_analysis.jl

@@ -1,7 +1,7 @@
 """
  flux_variability_analysis(
  model::MetabolicModel,
- reactions::Vector{Int},
+ fluxes::Vector{Int},
  optimizer;
  modifications = [],
  workers = [myid()],
@@ -11,9 +11,9 @@
  )::Matrix{Float64}
 
 Flux variability analysis solves a pair of optimization problems in `model` for
-each flux listed in `reactions`:
+each flux `f` described in `fluxes`:
 ```
- min,max xᵢ
+min,max fᵀxᵢ
 s.t. S x = b
  xₗ ≤ x ≤ xᵤ
  cᵀx ≥ bounds(Z₀)[1]
@@ -38,13 +38,13 @@ multithreaded variability computation can be used to improve resource
 allocation efficiency in many common use-cases.
 
 `ret` is a function used to extract results from optimized JuMP models of the
-individual reactions. By default, it calls and returns the value of
+individual fluxes. By default, it calls and returns the value of
 `JuMP.objective_value`. More information can be extracted e.g. by setting it to
 a function that returns a more elaborate data structure; such as `m ->
 (JuMP.objective_value(m), JuMP.value.(m[:x]))`.
 
 Returns a matrix of extracted `ret` values for minima and maxima, of total size
-(`length(reactions)`,2). The optimizer result status is checked with
+(`size(fluxes,2)`,2). The optimizer result status is checked with
 [`is_solved`](@ref); `nothing` is returned if the optimization failed for any
 reason.
 
@@ -56,16 +56,21 @@ flux_variability_analysis(model, [1, 2, 3, 42], GLPK.optimizer)
 """
 function flux_variability_analysis(
  model::MetabolicModel,
- reactions::Vector{Int},
+ fluxes::SparseMat,
  optimizer;
  modifications = [],
  workers = [myid()],
  optimal_objective_value = nothing,
  bounds = z -> (z, Inf),
  ret = objective_value,
 )
- if any(reactions .< 1) || any(reactions .> n_reactions(model))
- throw(DomainError(reactions, "Index exceeds number of reactions."))
+ if size(fluxes, 1) != n_reactions(model)
+ throw(
+ DomainError(
+ size(fluxes, 1),
+ "Flux matrix size is not compatible with model reaction count.",
+ ),
+ )
  end
 
  Z = bounds(
@@ -75,6 +80,8 @@ function flux_variability_analysis(
  ) : optimal_objective_value,
  )
 
+ flux_vector = [fluxes[:, i] for i = 1:size(fluxes, 2)]
+
  return screen_optmodel_modifications(
  model,
  optimizer;
@@ -93,12 +100,36 @@ function flux_variability_analysis(
  end,
  ],
  ),
- args = Tuple.([-reactions reactions]),
- analysis = (_, opt_model, ridx) -> _max_variability_flux(opt_model, ridx, ret),
+ args = tuple.([flux_vector flux_vector], [MIN_SENSE MAX_SENSE]),
+ analysis = (_, opt_model, flux, sense) ->
+ _max_variability_flux(opt_model, flux, sense, ret),
  workers = workers,
  )
 end
 
+"""
+ flux_variability_analysis(model::MetabolicModel, flux_indexes::Vector{Int}, optimizer; kwargs...)
+
+An overload of [`flux_variability_analysis`](@ref) that explores the fluxes specified by integer indexes
+"""
+function flux_variability_analysis(
+ model::MetabolicModel,
+ flux_indexes::Vector{Int},
+ optimizer;
+ kwargs...,
+)
+ if any((flux_indexes .< 1) .| (flux_indexes .> n_fluxes(model)))
+ throw(DomainError(flux_indexes, "Flux index out of range"))
+ end
+
+ flux_variability_analysis(
+ model,
+ reaction_flux(model)[:, flux_indexes],
+ optimizer;
+ kwargs...,
+ )
+end
+
 """
  flux_variability_analysis(
  model::MetabolicModel,
@@ -107,12 +138,10 @@ end
  )
 
 A simpler version of [`flux_variability_analysis`](@ref) that maximizes and
-minimizes all reactions in the model. Arguments are forwarded.
+minimizes all declared fluxes in the model. Arguments are forwarded.
 """
-function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
- n = n_reactions(model)
- return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
-end
+flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...) =
+ flux_variability_analysis(model, reaction_flux(model), optimizer; kwargs...)
 
 """
  flux_variability_analysis_dict(
@@ -122,8 +151,8 @@ end
  )
 
 A variant of [`flux_variability_analysis`](@ref) that returns the individual
-maximized and minimized fluxes of all reactions as two dictionaries (of
-dictionaries). All keyword arguments except `ret` are passed through.
+maximized and minimized fluxes as two dictionaries (of dictionaries). All
+keyword arguments except `ret` are passed through.
 
 # Example
 ```
@@ -140,24 +169,71 @@ mins, maxs = flux_variability_analysis_dict(
 ```
 """
 function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
- fluxes = flux_variability_analysis(model, optimizer; kwargs..., ret = flux_vector)
- rxns = reactions(model)
- dicts = zip.(Ref(rxns), fluxes)
+ vs = flux_variability_analysis(
+ model,
+ optimizer;
+ kwargs...,
+ ret = sol -> flux_vector(model, sol),
+ )
+ flxs = fluxes(model)
+ dicts = zip.(Ref(flxs), vs)
 
- return (Dict(rxns .=> Dict.(dicts[:, 1])), Dict(rxns .=> Dict.(dicts[:, 2])))
+ return (Dict(flxs .=> Dict.(dicts[:, 1])), Dict(flxs .=> Dict.(dicts[:, 2])))
 end
 
 """
- _max_variability_flux(opt_model, rid, ret)
+ _max_variability_flux(opt_model, flux, sense, ret)
 
 Internal helper for maximizing reactions in optimization model.
 """
-function _max_variability_flux(opt_model, rid, ret)
- sense = rid > 0 ? MAX_SENSE : MIN_SENSE
- var = all_variables(opt_model)[abs(rid)]
-
- @objective(opt_model, sense, var)
+function _max_variability_flux(opt_model, flux, sense, ret)
+ @objective(opt_model, sense, sum(flux .* opt_model[:x]))
  optimize!(opt_model)
 
  is_solved(opt_model) ? ret(opt_model) : nothing
 end
+
+"""
+ reaction_variability_analysis(model::MetabolicModel, reaction_indexes::Vector{Int}, optimizer; kwargs...)
+
+A variant for [`flux_variability_analysis`](@ref) that examines actual
+reactions (selected by their indexes in `reactions` argument) instead of whole
+fluxes. This may be useful for models where the sets of reactions and fluxes
+differ.
+"""
+function reaction_variability_analysis(
+ model::MetabolicModel,
+ reaction_indexes::Vector{Int},
+ optimizer;
+ kwargs...,
+)
+ if any((reaction_indexes .< 1) .| (reaction_indexes .> n_reactions(model)))
+ throw(DomainError(reaction_indexes, "Flux index out of range"))
+ end
+
+ flux_variability_analysis(
+ model,
+ sparse(
+ reaction_indexes,
+ 1:length(reaction_indexes),
+ 1.0,
+ n_reactions(model),
+ length(reaction_indexes),
+ ),
+ optimizer;
+ kwargs...,
+ )
+end
+
+"""
+ reaction_variability_analysis( model::MetabolicModel, optimizer; kwargs...)
+
+Shortcut for [`reaction_variability_analysis`](@ref) that examines all reactions.
+"""
+reaction_variability_analysis(model::MetabolicModel, optimizer; kwargs...) =
+ reaction_variability_analysis(
+ model,
+ collect(1:n_reactions(model)),
+ optimizer;
+ kwargs...,
+ )