Develop → master merge for 1.4.3

Project.toml

@@ -1,7 +1,7 @@
 name = "COBREXA"
 uuid = "babc4406-5200-4a30-9033-bf5ae714c842"
 authors = ["The developers of COBREXA.jl"]
-version = "1.4.2"
+version = "1.4.3"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"

src/analysis/flux_variability_analysis.jl

@@ -42,7 +42,7 @@ reason.
 # Example
 ```
 model = load_model("e_coli_core.json")
-flux_variability_analysis(model, [1, 2, 3, 42], GLPK.optimizer)
+flux_variability_analysis(model, GLPK.optimizer)
 ```
 """
 function flux_variability_analysis(

src/base/types/SBMLModel.jl

@@ -9,44 +9,103 @@ $(TYPEDFIELDS)
 """
 struct SBMLModel <: MetabolicModel
     sbml::SBML.Model
+    reaction_ids::Vector{String}
+    reaction_idx::Dict{String,Int}
+    metabolite_ids::Vector{String}
+    metabolite_idx::Dict{String,Int}
+    gene_ids::Vector{String}
+    active_objective::String
+end
+
+"""
+$(TYPEDEF)
+
+Construct the SBML model and add the necessary cached indexes, possibly choosing an active objective.
+"""
+function SBMLModel(sbml::SBML.Model, active_objective::String = "")
+    rxns = sort(collect(keys(sbml.reactions)))
+    mets = sort(collect(keys(sbml.species)))
+    genes = sort(collect(keys(sbml.gene_products)))
+
+    SBMLModel(
+        sbml,
+        rxns,
+        Dict(rxns .=> eachindex(rxns)),
+        mets,
+        Dict(mets .=> eachindex(mets)),
+        genes,
+        active_objective,
+    )
 end
 
 """
 $(TYPEDSIGNATURES)
 
 Get reactions from a [`SBMLModel`](@ref).
 """
-reactions(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.reactions)]
+reactions(model::SBMLModel)::Vector{String} = model.reaction_ids
 
 """
 $(TYPEDSIGNATURES)
 
 Get metabolites from a [`SBMLModel`](@ref).
 """
-metabolites(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.species)]
+metabolites(model::SBMLModel)::Vector{String} = model.metabolite_ids
 
 """
 $(TYPEDSIGNATURES)
 
 Efficient counting of reactions in [`SBMLModel`](@ref).
 """
-n_reactions(model::SBMLModel)::Int = length(model.sbml.reactions)
+n_reactions(model::SBMLModel)::Int = length(model.reaction_ids)
 
 """
 $(TYPEDSIGNATURES)
 
 Efficient counting of metabolites in [`SBMLModel`](@ref).
 """
-n_metabolites(model::SBMLModel)::Int = length(model.sbml.species)
+n_metabolites(model::SBMLModel)::Int = length(model.metabolite_ids)
 
 """
 $(TYPEDSIGNATURES)
 
 Recreate the stoichiometry matrix from the [`SBMLModel`](@ref).
 """
 function stoichiometry(model::SBMLModel)::SparseMat
-    _, _, S = SBML.stoichiometry_matrix(model.sbml)
-    return S
+
+    # find the vector size for preallocation
+    nnz = 0
+    for (_, r) in model.sbml.reactions
+        for _ in r.reactants
+            nnz += 1
+        end
+        for _ in r.products
+            nnz += 1
+        end
+    end
+
+    Rows = Int[]
+    Cols = Int[]
+    Vals = Float64[]
+    sizehint!(Rows, nnz)
+    sizehint!(Cols, nnz)
+    sizehint!(Vals, nnz)
+
+    row_idx = Dict(k => i for (i, k) in enumerate(model.metabolite_ids))
+    for (ridx, rid) in enumerate(model.reaction_ids)
+        r = model.sbml.reactions[rid]
+        for sr in r.reactants
+            push!(Rows, model.metabolite_idx[sr.species])
+            push!(Cols, ridx)
+            push!(Vals, isnothing(sr.stoichiometry) ? -1.0 : -sr.stoichiometry)
+        end
+        for sr in r.products
+            push!(Rows, model.metabolite_idx[sr.species])
+            push!(Cols, ridx)
+            push!(Vals, isnothing(sr.stoichiometry) ? 1.0 : sr.stoichiometry)
+        end
+    end
+    return sparse(Rows, Cols, Vals, n_metabolites(model), n_reactions(model))
 end
 
 """
@@ -56,21 +115,56 @@ Get the lower and upper flux bounds of model [`SBMLModel`](@ref). Throws `Domain
 case if the SBML contains mismatching units.
 """
 function bounds(model::SBMLModel)::Tuple{Vector{Float64},Vector{Float64}}
-    lbu, ubu = SBML.flux_bounds(model.sbml)
-
-    unit = lbu[1][2]
-    getvalue = (val, _)::Tuple -> val
-    getunit = (_, unit)::Tuple -> unit
+    # There are multiple ways in SBML to specify a lower/upper bound. There are
+    # the "global" model bounds that we completely ignore now because no one
+    # uses them. In reaction, you can specify the bounds using "LOWER_BOUND"
+    # and "UPPER_BOUND" parameters, but also there may be a FBC plugged-in
+    # parameter name that refers to the parameters.  We extract these, using
+    # the units from the parameters. For unbounded reactions we use -Inf or Inf
+    # as a default.
+
+    common_unit = ""
+
+    function get_bound(rid, fld, param, default)
+        rxn = model.sbml.reactions[rid]
+        param_name = SBML.mayfirst(getfield(rxn, fld), param)
+        param = get(
+            rxn.kinetic_parameters,
+            param_name,
+            get(model.sbml.parameters, param_name, default),
+        )
+        unit = SBML.mayfirst(param.units, "")
+        if unit != ""
+            if common_unit != ""
+                if unit != common_unit
+                    throw(
+                        DomainError(
+                            units_in_sbml,
+                            "The SBML file uses multiple units; loading would need conversion",
+                        ),
+                    )
+                end
+            else
+                common_unit = unit
+            end
+        end
+        return param.value
+    end
 
-    units_in_sbml = filter(!isempty, unique([getunit.(lbu) getunit.(ubu)]))
-    length(units_in_sbml) <= 1 || throw(
-        DomainError(
-            units_in_sbml,
-            "The SBML file uses multiple units; loading would need conversion",
+    return (
+        get_bound.(
+            model.reaction_ids,
+            :lower_bound,
+            "LOWER_BOUND",
+            Ref(SBML.Parameter(value = -Inf)),
+        ),
+        get_bound.(
+            model.reaction_ids,
+            :upper_bound,
+            "UPPER_BOUND",
+            Ref(SBML.Parameter(value = Inf)),
         ),
     )
-
-    return (getvalue.(lbu), getvalue.(ubu))
 end
 
 """
@@ -85,21 +179,41 @@ $(TYPEDSIGNATURES)
 
 Objective of the [`SBMLModel`](@ref).
 """
-objective(model::SBMLModel)::SparseVec = SBML.flux_objective(model.sbml)
+function objective(model::SBMLModel)::SparseVec
+    res = sparsevec([], [], n_reactions(model))
+
+    objective = get(model.sbml.objectives, model.active_objective, nothing)
+    if isnothing(objective) && length(model.sbml.objectives) == 1
+        objective = first(values(model.sbml.objectives))
+    end
+    if !isnothing(objective)
+        direction = objective.type == "maximize" ? 1.0 : -1.0
+        for (rid, coef) in objective.flux_objectives
+            res[model.reaction_idx[rid]] = float(direction * coef)
+        end
+    else
+        # old-style objectives
+        for (rid, r) in model.sbml.reactions
+            oc = get(r.kinetic_parameters, "OBJECTIVE_COEFFICIENT", nothing)
+            isnothing(oc) || (res[model.reaction_idx[rid]] = float(oc.value))
+        end
+    end
+    return res
+end
 
 """
 $(TYPEDSIGNATURES)
 
 Get genes of a [`SBMLModel`](@ref).
 """
-genes(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.gene_products)]
+genes(model::SBMLModel)::Vector{String} = model.gene_ids
 
 """
 $(TYPEDSIGNATURES)
 
 Get number of genes in [`SBMLModel`](@ref).
 """
-n_genes(model::SBMLModel)::Int = length(model.sbml.gene_products)
+n_genes(model::SBMLModel)::Int = length(model.gene_ids)
 
 """
 $(TYPEDSIGNATURES)

test/analysis/envelopes.jl

@@ -5,15 +5,15 @@
     rxns = [1, 2, 3]
 
     lat = collect.(envelope_lattice(m, rxns; samples = 3))
-    @test lat == [[0, 500, 1000], [-1000, 0, 1000], [-1000, 0, 1000]]
+    @test lat == [[-1000.0, 0.0, 1000.0], [-1000.0, 0.0, 1000.0], [-1000.0, 0.0, 1000.0]]
     @test lat == collect.(envelope_lattice(m, reactions(m)[rxns]; samples = 3))
 
     vals =
         objective_envelope(
             m,
             reactions(m)[rxns],
             Tulip.Optimizer;
-            lattice_args = (samples = 3, ranges = [(-5, 5), (-5, 5), (-5, 5)]),
+            lattice_args = (samples = 3, ranges = [(-5, 0), (-5, 0), (-5, 5)]),
             workers = W,
         ).values
 
@@ -22,18 +22,18 @@
     @test isapprox(
         filter(!isnothing, vals),
         [
-            0.0,
-            0.0,
-            0.3809552834584839,
-            0.0,
-            0.0,
-            0.0,
-            0.7619105669322938,
-            0.38095528320252253,
-            0.0,
-            0.0,
-            1.0542729305608431,
-            0.7619105589042587,
+            0.5833914451564178,
+            0.5722033617617296,
+            0.673404874265479,
+            0.5610152783118744,
+            0.6606736594947443,
+            0.7593405739802911,
+            0.7020501075121626,
+            0.689318892439569,
+            0.787985806919745,
+            0.6765876776826879,
+            0.7752545924613183,
+            0.8739215069575214,
         ],
         atol = TEST_TOLERANCE,
     )

test/io/sbml.jl

@@ -10,8 +10,8 @@
     @test length.(bounds(m)) == (95, 95)
     @test all([length(m.xl), length(m.xu), length(m.c)] .== 95)
 
-    @test metabolites(m)[1:3] == ["M_succoa_c", "M_ac_c", "M_fru_b"]
-    @test reactions(m)[1:3] == ["R_EX_fum_e", "R_ACONTb", "R_GLNS"]
+    @test metabolites(m)[1:3] == ["M_13dpg_c", "M_2pg_c", "M_3pg_c"]
+    @test reactions(m)[1:3] == ["R_ACALD", "R_ACALDt", "R_ACKr"]
 
     cm = convert(CoreModelCoupled, sbmlm)
     @test n_coupling_constraints(cm) == 0