Merge remote-tracking branch 'origin/develop' into next

Project.toml

@@ -32,7 +32,7 @@ JuMP = "1"
 MAT = "0.10"
 MacroTools = "0.5.6"
 OrderedCollections = "1.4"
-SBML = "~1.1, ~1.2"
+SBML = "~1.3"
 StableRNGs = "1.0"
 Tulip = "0.7.0, 0.8.0, 0.9.2"
 julia = "1.5"

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 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"
@@ -10,7 +11,6 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [compat]

docs/src/examples/05b_fba_mods.jl

@@ -10,7 +10,7 @@
 !isfile("e_coli_core.xml") &&
     download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")
 
-using COBREXA, GLPK, Tulip
+using COBREXA, GLPK, Tulip, JuMP
 
 model = load_model("e_coli_core.xml")
 
@@ -33,6 +33,6 @@ fluxes = flux_balance_analysis_dict(
         knockout(["b0978", "b0734"]), # knock out two genes
         change_optimizer(Tulip.Optimizer), # ignore the above optimizer and switch to Tulip
         change_optimizer_attribute("IPM_IterationsLimit", 1000), # customize Tulip
-        change_sense(MAX_SENSE), # explicitly tell Tulip to maximize the objective
+        change_sense(JuMP.MAX_SENSE), # explicitly tell Tulip to maximize the objective
     ],
 )

docs/src/examples/08_pfba.jl

@@ -43,7 +43,6 @@ fluxes = parsimonious_flux_balance_analysis_dict(
     Clarabel.Optimizer;
     modifications = [
         silence, # optionally silence the optimizer (Clarabel is very verbose by default)
-        change_optimizer_attribute("polish", true), # tell Clarabel to invest time into improving the precision of the solution
         change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12), # fix glucose consumption rate
     ],
 )
@@ -70,7 +69,6 @@ flux_vector = parsimonious_flux_balance_analysis_vec(
     ],
     qp_modifications = [
         change_optimizer(Clarabel.Optimizer), # now switch to Clarabel (Tulip wouldn't be able to finish the computation)
-        change_optimizer_attribute("polish", true), # get an accurate solution, see Clarabel's documentation
         silence, # and make it quiet.
     ],
 )

docs/src/examples/13_moma.jl

@@ -25,11 +25,8 @@ using Clarabel
 
 # We will need a reference solution, which represents the original state of the
 # organism before the change.
-reference_flux = flux_balance_analysis_dict(
-    model,
-    Clarabel.Optimizer;
-    modifications = [silence, change_optimizer_attribute("polish", true)],
-)
+reference_flux =
+    flux_balance_analysis_dict(model, Clarabel.Optimizer; modifications = [silence])
 
 # As the change here, we manually knock out CYTBD reaction:
 changed_model = change_bound(model, "R_CYTBD", lower = 0.0, upper = 0.0);
@@ -41,7 +38,7 @@ flux_summary(
         changed_model,
         reference_flux,
         Clarabel.Optimizer;
-        modifications = [silence, change_optimizer_attribute("polish", true)],
+        modifications = [silence],
     ),
 )
 
@@ -52,6 +49,6 @@ flux_summary(
     flux_balance_analysis_dict(
         changed_model,
         Clarabel.Optimizer;
-        modifications = [silence, change_optimizer_attribute("polish", true)],
+        modifications = [silence],
     ),
 )

src/misc/constants.jl

@@ -26,7 +26,7 @@ const constants = (
         ids = ["id", "description"],
         metformulas = ["metFormula", "metFormulas"],
         metcharges = ["metCharge", "metCharges"],
-        metcompartments = ["metCompartment", "metCompartments"],
+        metcompartments = ["metCompartment", "metCompartments", "metComp", "metComps"],
         rxnnames = ["rxnNames"],
         metnames = ["metNames"],
     ),

src/types/models/MATModel.jl

@@ -155,12 +155,12 @@ $(TYPEDSIGNATURES)
 
 Extract metabolite charge from `metCharge` or `metCharges`.
 """
-function Accessors.metabolite_charge(m::MATModel, mid::String)
+function Accessors.metabolite_charge(m::MATModel, mid::String)::Maybe{Int}
     met_charge = maybemap(
         x -> x[findfirst(==(mid), metabolites(m))],
         gets(m.mat, nothing, constants.keynames.metcharges),
     )
-    isnan(met_charge) ? 0 : met_charge
+    _maybemap(Int, isnan(met_charge) ? nothing : met_charge)
 end
 
 """

src/types/models/SBMLModel.jl

@@ -123,24 +123,66 @@ Accessors.metabolite_formula(model::SBMLModel, mid::String)::Maybe{MetaboliteFor
 """
 $(TYPEDSIGNATURES)
 
+Get the compartment of a chosen metabolite from [`SBMLModel`](@ref).
+"""
+metabolite_compartment(model::SBMLModel, mid::String) = model.sbml.species[mid].compartment
+
+"""
+$(TYPEDSIGNATURES)
+
 Get charge of a chosen metabolite from [`SBMLModel`](@ref).
 """
 Accessors.metabolite_charge(model::SBMLModel, mid::String)::Maybe{Int} =
     model.sbml.species[mid].charge
 
-function _sbml_export_annotation(annotation)::Maybe{String}
-    if isnothing(annotation) || isempty(annotation)
-        nothing
-    elseif length(annotation) != 1 || first(annotation).first != ""
-        @io_log @warn "Possible data loss: multiple annotations converted to text for SBML" annotation
-        join(["$k: $v" for (k, v) in annotation], "\n")
-    else
-        @io_log @warn "Possible data loss: trying to represent annotation in SBML is unlikely to work " annotation
-        first(annotation).second
+function _parse_sbml_identifiers_org_uri(uri::String)::Tuple{String,String}
+    m = match(r"^http://identifiers.org/([^/]+)/(.*)$", uri)
+    isnothing(m) ? ("RESOURCE_URI", uri) : (m[1], m[2])
+end
+
+function _sbml_import_cvterms(sbo::Maybe{String}, cvs::Vector{SBML.CVTerm})::Annotations
+    res = Annotations()
+    isnothing(sbo) || (res["sbo"] = [sbo])
+    for cv in cvs
+        cv.biological_qualifier == :is || continue
+        for (id, val) in _parse_sbml_identifiers_org_uri.(cv.resource_uris)
+            push!(get!(res, id, []), val)
+        end
+    end
+    return res
+end
+
+function _sbml_export_cvterms(annotations::Annotations)::Vector{SBML.CVTerm}
+    isempty(annotations) && return []
+    length(annotations) == 1 && haskey(annotations, "sbo") && return []
+    [
+        SBML.CVTerm(
+            biological_qualifier = :is,
+            resource_uris = [
+                id == "RESOURCE_URI" ? val : "http://identifiers.org/$id/$val" for
+                (id, vals) = annotations if id != "sbo" for val in vals
+            ],
+        ),
+    ]
+end
+
+function _sbml_export_sbo(annotations::Annotations)::Maybe{String}
+    haskey(annotations, "sbo") || return nothing
+    if length(annotations["sbo"]) != 1
+        @_io_log @error "Data loss: SBO term is not unique for SBML export" annotations["sbo"]
+        return
     end
+    return annotations["sbo"][1]
+end
+
+function _sbml_import_notes(notes::Maybe{String})::Notes
+    isnothing(notes) ? Notes() : Notes("" => [notes])
 end
 
-const _sbml_export_notes = _sbml_export_annotation
+function _sbml_export_notes(notes::Notes)::Maybe{String}
+    isempty(notes) || @_io_log @error "Data loss: notes not exported to SBML" notes
+    nothing
+end
 
 """
 $(TYPEDSIGNATURES)
@@ -183,6 +225,56 @@ Accessors.gene_name(model::SBMLModel, gid::String) = model.sbml.gene_products[gi
 """
 $(TYPEDSIGNATURES)
 
+Return the annotations of reaction with ID `rid`.
+"""
+reaction_annotations(model::SBMLModel, rid::String) =
+    _sbml_import_cvterms(model.sbml.reactions[rid].sbo, model.sbml.reactions[rid].cv_terms)
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the annotations of metabolite with ID `mid`.
+"""
+metabolite_annotations(model::SBMLModel, mid::String) =
+    _sbml_import_cvterms(model.sbml.species[mid].sbo, model.sbml.species[mid].cv_terms)
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the annotations of gene with ID `gid`.
+"""
+gene_annotations(model::SBMLModel, gid::String) = _sbml_import_cvterms(
+    model.sbml.gene_products[gid].sbo,
+    model.sbml.gene_products[gid].cv_terms,
+)
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the notes about reaction with ID `rid`.
+"""
+reaction_notes(model::SBMLModel, rid::String) =
+    _sbml_import_notes(model.sbml.reactions[rid].notes)
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the notes about metabolite with ID `mid`.
+"""
+metabolite_notes(model::SBMLModel, mid::String) =
+    _sbml_import_notes(model.sbml.species[mid].notes)
+
+"""
+$(TYPEDSIGNATURES)
+
+Return the notes about gene with ID `gid`.
+"""
+gene_notes(model::SBMLModel, gid::String) =
+    _sbml_import_notes(model.sbml.gene_products[gid].notes)
+
+"""
+$(TYPEDSIGNATURES)
+
 Convert any metabolic model to [`SBMLModel`](@ref).
 """
 function Base.convert(::Type{SBMLModel}, mm::MetabolicModel)
@@ -197,38 +289,45 @@ function Base.convert(::Type{SBMLModel}, mm::MetabolicModel)
     comps = default.("compartment", metabolite_compartment.(Ref(mm), mets))
     compss = Set(comps)
 
+    metid(x) = startswith(x, "M_") ? x : "M_$x"
+    rxnid(x) = startswith(x, "R_") ? x : "R_$x"
+    gprid(x) = startswith(x, "G_") ? x : "G_$x"
+
     return SBMLModel(
         SBML.Model(
             compartments = Dict(
                 comp => SBML.Compartment(constant = true) for comp in compss
             ),
             species = Dict(
-                mid => SBML.Species(
+                metid(mid) => SBML.Species(
                     name = metabolite_name(mm, mid),
                     compartment = default("compartment", comps[mi]),
                     formula = metabolite_formula(mm, mid),
+                    formula = maybemap(unparse_formula, metabolite_formula(mm, mid)),
                     charge = metabolite_charge(mm, mid),
                     constant = false,
                     boundary_condition = false,
                     only_substance_units = false,
+                    sbo = _sbml_export_sbo(metabolite_annotations(mm, mid)),
                     notes = _sbml_export_notes(metabolite_notes(mm, mid)),
-                    annotation = _sbml_export_annotation(metabolite_annotations(mm, mid)),
+                    metaid = metid(mid),
+                    cv_terms = _sbml_export_cvterms(metabolite_annotations(mm, mid)),
                 ) for (mi, mid) in enumerate(mets)
             ),
             reactions = Dict(
-                rid => SBML.Reaction(
+                rxnid(rid) => SBML.Reaction(
                     name = reaction_name(mm, rid),
                     reactants = [
                         SBML.SpeciesReference(
-                            species = mets[i],
+                            species = metid(mets[i]),
                             stoichiometry = -stoi[i, ri],
                             constant = true,
                         ) for
                         i in SparseArrays.nonzeroinds(stoi[:, ri]) if stoi[i, ri] <= 0
                     ],
                     products = [
                         SBML.SpeciesReference(
-                            species = mets[i],
+                            species = metid(mets[i]),
                             stoichiometry = stoi[i, ri],
                             constant = true,
                         ) for
@@ -245,16 +344,20 @@ function Base.convert(::Type{SBMLModel}, mm::MetabolicModel)
                         reaction_gene_association(mm, rid),
                     ),
                     reversible = true,
+                    sbo = _sbml_export_sbo(reaction_annotations(mm, rid)),
                     notes = _sbml_export_notes(reaction_notes(mm, rid)),
-                    annotation = _sbml_export_annotation(reaction_annotations(mm, rid)),
+                    metaid = rxnid(rid),
+                    cv_terms = _sbml_export_cvterms(reaction_annotations(mm, rid)),
                 ) for (ri, rid) in enumerate(rxns)
             ),
             gene_products = Dict(
-                gid => SBML.GeneProduct(
+                gprid(gid) => SBML.GeneProduct(
                     label = gid,
                     name = gene_name(mm, gid),
+                    sbo = _sbml_export_sbo(gene_annotations(mm, gid)),
                     notes = _sbml_export_notes(gene_notes(mm, gid)),
-                    annotation = _sbml_export_annotation(gene_annotations(mm, gid)),
+                    metaid = gprid(gid),
+                    cv_terms = _sbml_export_cvterms(gene_annotations(mm, gid)),
                 ) for gid in genes(mm)
             ),
             active_objective = "objective",