diff --git a/Project.toml b/Project.toml index 8d9395a..cea88ae 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SBML" uuid = "e5567a89-2604-4b09-9718-f5f78e97c3bb" authors = ["Mirek Kratochvil ", "LCSB R3 team "] -version = "1.0.1" +version = "1.1.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" diff --git a/src/readsbml.jl b/src/readsbml.jl index 419806d..580d2ff 100644 --- a/src/readsbml.jl +++ b/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) diff --git a/src/structs.jl b/src/structs.jl index a0ba7d1..f0084a2 100644 --- a/src/structs.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 88884ae..0edf193 100644 --- a/src/utils.jl +++ b/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)) diff --git a/src/writesbml.jl b/src/writesbml.jl index b2227d7..90b6994 100644 --- a/src/writesbml.jl +++ b/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 = diff --git a/test/loadmodels.jl b/test/loadmodels.jl index eee0e94..3b878ea 100644 --- a/test/loadmodels.jl +++ b/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