Merge pull request #224 from LCSB-BioCore/mk-species-reference

implement full species references

Project.toml

@@ -1,7 +1,7 @@
 name = "SBML"
 uuid = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
 authors = ["Mirek Kratochvil <miroslav.kratochvil@uni.lu>", "LCSB R3 team <lcsb-r3@uni.lu>"]
-version = "1.0.1"
+version = "1.1.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/readsbml.jl

@@ -411,14 +411,34 @@ function get_model(mdl::VPtr)::SBML.Model
             )
         end
 
-        # extract stoichiometry
-        reactants = Dict{String,Float64}()
-        products = Dict{String,Float64}()
+        # extract species references
+        reactants = SpeciesReference[]
+        products = SpeciesReference[]
 
         add_stoi(sr, coll) = begin
-            s = get_string(sr, :SpeciesReference_getSpecies)
-            coll[s] =
-                ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr)
+            stoichiometry =
+                if ccall(sbml(:SpeciesReference_isSetStoichiometry), Cint, (VPtr,), sr) != 0
+                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr)
+                end
+
+            constant =
+                if ccall(sbml(:SpeciesReference_isSetConstant), Cint, (VPtr,), sr) != 0
+                    ccall(sbml(:SpeciesReference_getConstant), Cint, (VPtr,), sr) != 0
+                end
+
+            push!(
+                coll,
+                SpeciesReference(
+                    id = get_optional_string(
+                        sr,
+                        :SpeciesReference_isSetId,
+                        :SpeciesReference_getId,
+                    ),
+                    species = get_string(sr, :SpeciesReference_getSpecies);
+                    stoichiometry,
+                    constant,
+                ),
+            )
         end
 
         for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)

src/structs.jl

@@ -199,6 +199,21 @@ end
 """
 $(TYPEDEF)
 
+SBML SpeciesReference.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+Base.@kwdef struct SpeciesReference
+    id::Maybe{String} = nothing
+    species::String
+    stoichiometry::Maybe{Float64} = nothing
+    constant::Maybe{Bool} = nothing
+end
+
+"""
+$(TYPEDEF)
+
 Reaction with stoichiometry that assigns reactants and products their relative
 consumption/production rates, lower/upper bounds (in tuples `lb` and `ub`, with
 unit names), and objective coefficient (`oc`). Also may contains `notes` and
@@ -209,8 +224,8 @@ $(TYPEDFIELDS)
 """
 Base.@kwdef struct Reaction
     name::Maybe{String} = nothing
-    reactants::Dict{String,Float64} = Dict()
-    products::Dict{String,Float64} = Dict()
+    reactants::Vector{SpeciesReference} = []
+    products::Vector{SpeciesReference} = []
     kinetic_parameters::Dict{String,Parameter} = Dict()
     lower_bound::Maybe{String} = nothing
     upper_bound::Maybe{String} = nothing

src/utils.jl

@@ -64,15 +64,15 @@ function stoichiometry_matrix(m::SBML.Model)
 
     for (rid, r) in m.reactions
         ridx = col_idx[rid]
-        for (sid, stoi) in r.reactants
-            push!(SI, row_idx[sid])
+        for sr in r.reactants
+            push!(SI, row_idx[sr.species])
             push!(RI, ridx)
-            push!(SV, -stoi)
+            push!(SV, isnothing(sr.stoichiometry) ? -1.0 : -sr.stoichiometry)
         end
-        for (sid, stoi) in r.products
-            push!(SI, row_idx[sid])
+        for sr in r.products
+            push!(SI, row_idx[sr.species])
             push!(RI, ridx)
-            push!(SV, stoi)
+            push!(SV, isnothing(sr.stoichiometry) ? 1.0 : sr.stoichiometry)
         end
     end
     return rows, cols, SparseArrays.sparse(SI, RI, SV, length(rows), length(cols))

src/writesbml.jl

@@ -325,56 +325,41 @@ function model_to_sbml!(doc::VPtr, mdl::Model)::VPtr
             reaction_ptr,
             reaction.name,
         )
-        for (species, stoichiometry) in reaction.reactants
-            reactant_ptr =
-                ccall(sbml(:Reaction_createReactant), VPtr, (VPtr,), reaction_ptr)
-            ccall(
-                sbml(:SpeciesReference_setSpecies),
-                Cint,
-                (VPtr, Cstring),
-                reactant_ptr,
-                species,
-            )
-            ccall(
-                sbml(:SpeciesReference_setStoichiometry),
-                Cint,
-                (VPtr, Cdouble),
-                reactant_ptr,
-                stoichiometry,
-            )
-            # Assume constant reactant for the time being
-            ccall(
-                sbml(:SpeciesReference_setConstant),
-                Cint,
-                (VPtr, Cint),
-                reactant_ptr,
-                true,
-            )
-        end
-        for (species, stoichiometry) in reaction.products
-            product_ptr = ccall(sbml(:Reaction_createProduct), VPtr, (VPtr,), reaction_ptr)
-            ccall(
-                sbml(:SpeciesReference_setSpecies),
-                Cint,
-                (VPtr, Cstring),
-                product_ptr,
-                species,
-            )
-            ccall(
-                sbml(:SpeciesReference_setStoichiometry),
-                Cint,
-                (VPtr, Cdouble),
-                product_ptr,
-                stoichiometry,
-            )
-            # Assume constant product for the time being
-            ccall(
-                sbml(:SpeciesReference_setConstant),
-                Cint,
-                (VPtr, Cint),
-                product_ptr,
-                true,
-            )
+        for (sr_create, srs) in [
+            :Reaction_createReactant => reaction.reactants,
+            :Reaction_createProduct => reaction.products,
+        ]
+            for sr in srs
+                reactant_ptr = ccall(sbml(sr_create), VPtr, (VPtr,), reaction_ptr)
+                ccall(
+                    sbml(:SpeciesReference_setSpecies),
+                    Cint,
+                    (VPtr, Cstring),
+                    reactant_ptr,
+                    sr.species,
+                )
+                isnothing(sr.id) || ccall(
+                    sbml(:SpeciesReference_setId),
+                    Cint,
+                    (VPtr, Cstring),
+                    reactant_ptr,
+                    sr.id,
+                )
+                isnothing(sr.stoichiometry) || ccall(
+                    sbml(:SpeciesReference_setStoichiometry),
+                    Cint,
+                    (VPtr, Cdouble),
+                    reactant_ptr,
+                    sr.stoichiometry,
+                )
+                isnothing(sr.constant) || ccall(
+                    sbml(:SpeciesReference_setConstant),
+                    Cint,
+                    (VPtr, Cint),
+                    reactant_ptr,
+                    sr.constant,
+                )
+            end
         end
         if !isempty(reaction.kinetic_parameters) || !isnothing(reaction.kinetic_math)
             kinetic_law_ptr =

test/loadmodels.jl

@@ -406,8 +406,25 @@ end
     @test m.parameters["S1conv"].constant == true
 end
 
+function fix_constant!(model::SBML.Model)
+    # We only write SBML L3v2. If a model is L2 or less, the `constant`
+    # attributes may be missing (which is true for some models). We add the
+    # main problematic ones (in speciesReferences) here, to make sure the
+    # round trip has a chance to finish.
+    _clean(sr::SBML.SpeciesReference) = SBML.SpeciesReference(
+        species = sr.species,
+        stoichiometry = sr.stoichiometry,
+        constant = isnothing(sr.constant) ? true : sr.constant,
+    )
+    for (_, r) in model.reactions
+        r.reactants .= map(_clean, r.reactants)
+        r.products .= map(_clean, r.products)
+    end
+end
+
 @testset "writeSBML" begin
     model = readSBML(joinpath(@__DIR__, "data", "Dasgupta2020.xml"))
+    fix_constant!(model)
     expected = read(joinpath(@__DIR__, "data", "Dasgupta2020-written.xml"), String)
     # Remove carriage returns, if any
     expected = replace(expected, '\r' => "")
@@ -424,6 +441,7 @@ end
     # with the original model.
     @testset "Round-trip - $(basename(file))" for file in first.(sbmlfiles)
         model = readSBML(file)
+        fix_constant!(model)
         round_trip_model = readSBMLFromString(@test_logs(writeSBML(model)))
         @test model.parameters == round_trip_model.parameters
         @test model.units == round_trip_model.units