Merge pull request #39 from LCSB-BioCore/mk-species-units

Extract kinetics-relevant information

Project.toml

 name = "SBML"
 uuid = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
 authors = ["Mirek Kratochvil <miroslav.kratochvil@uni.lu>", "LCSB R3 team <lcsb-r3@uni.lu>"]
-version = "0.3.2"
+version = "0.4.0"
 
 [deps]
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

src/SBML.jl

@@ -6,6 +6,7 @@ using SparseArrays
 include("structs.jl")
 include("version.jl")
 include("readsbml.jl")
+include("math.jl")
 include("utils.jl")
 
 sbml = (sym::Symbol) -> dlsym(SBML_jll.libsbml_handle, sym)

src/math.jl

+
+"""
+    ast_is(ast::VPtr, what::Symbol)::Bool
+
+Helper for quickly recognizing kinds of ASTs
+"""
+ast_is(ast::VPtr, what::Symbol)::Bool = ccall(sbml(what), Cint, (VPtr,), ast) != 0
+
+"""
+    parse_math_children(ast::VPtr)::Vector{Math}
+
+Recursively parse all children of an AST node.
+"""
+parse_math_children(ast::VPtr)::Vector{Math} = [
+    parse_math(ccall(sbml(:ASTNode_getChild), VPtr, (VPtr, Cuint), ast, i - 1)) for
+    i = 1:ccall(sbml(:ASTNode_getNumChildren), Cuint, (VPtr,), ast)
+]
+
+"""
+    parse_math(ast::VPtr)::Math
+
+This attempts to parse out a decent Julia-esque ([`Math`](@ref) AST from a
+pointer to `ASTNode_t`.
+"""
+function parse_math(ast::VPtr)::Math
+    if ast_is(ast, :ASTNode_isName) || ast_is(ast, :ASTNode_isConstant)
+        return MathIdent(get_string(ast, :ASTNode_getName))
+    elseif ast_is(ast, :ASTNode_isInteger)
+        return MathVal(ccall(sbml(:ASTNode_getInteger), Cint, (VPtr,), ast))
+    elseif ast_is(ast, :ASTNode_isReal)
+        return MathVal(ccall(sbml(:ASTNode_getReal), Cdouble, (VPtr,), ast))
+    elseif ast_is(ast, :ASTNode_isFunction)
+        return MathApply(get_string(ast, :ASTNode_getName), parse_math_children(ast))
+    elseif ast_is(ast, :ASTNode_isOperator)
+        return MathApply(
+            string(Char(ccall(sbml(:ASTNode_getCharacter), Cchar, (VPtr,), ast))),
+            parse_math_children(ast),
+        )
+    elseif ast_is(ast, :ASTNode_isLambda)
+        children = parse_math_children(ast)
+        if !isempty(children)
+            body = pop!(children)
+            return MathLambda(broadcast((x::MathIdent) -> x.id, children), body)
+        else
+            @warn "invalid function definition found"
+            return MathIdent("?invalid?")
+        end
+    else
+        @warn "unsupported math element found"
+        return MathIdent("?unsupported?")
+    end
+end

src/readsbml.jl

@@ -24,7 +24,7 @@ Like [`get_string`](@Ref), but returns `nothing` instead of throwing an
 exception.
 
 This is used to get notes and annotations and several other things (see
-`getNotes`, `getAnnotations`)
+`get_notes`, `get_annotations`)
 """
 function get_optional_string(x::VPtr, fn_sym)::Maybe{String}
     str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
@@ -35,6 +35,45 @@ function get_optional_string(x::VPtr, fn_sym)::Maybe{String}
     end
 end
 
+"""
+    get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}
+
+Helper for getting out boolean flags.
+"""
+function get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}
+    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
+        return ccall(sbml(get_sym), Cint, (VPtr,), x) != 0
+    else
+        return nothing
+    end
+end
+
+"""
+    get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{UInt}
+
+Helper for getting out unsigned integers.
+"""
+function get_optional_int(x::VPtr, is_sym, get_sym)::Maybe{Int}
+    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
+        return ccall(sbml(get_sym), Cint, (VPtr,), x)
+    else
+        return nothing
+    end
+end
+
+"""
+    get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}
+
+Helper for getting out C doubles aka Float64s.
+"""
+function get_optional_double(x::VPtr, is_sym, get_sym)::Maybe{Float64}
+    if ccall(sbml(is_sym), Cint, (VPtr,), x) != 0
+        return ccall(sbml(get_sym), Cdouble, (VPtr,), x)
+    else
+        return nothing
+    end
+end
+
 """
     function readSBML(fn::String)::Model
 
@@ -65,8 +104,8 @@ function readSBML(fn::String)::Model
     end
 end
 
-getNotes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
-getAnnotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
+get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
+get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
 
 """
     function getAssociation(x::VPtr)::GeneProductAssociation
@@ -143,18 +182,31 @@ function extractModel(mdl::VPtr)::Model
     end
 
     # parse out compartment names
-    compartments = [
-        get_string(
-            ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1),
-            :Compartment_getId,
-        ) for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
-    ]
+    compartments = Dict{String,Compartment}()
+    for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
+        co = ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1)
+
+        compartments[get_string(co, :Compartment_getId)] = Compartment(
+            get_optional_string(co, :Compartment_getName),
+            get_optional_bool(co, :Compartment_isSetConstant, :Compartment_getConstant),
+            get_optional_int(
+                co,
+                :Compartment_isSetSpatialDimensions,
+                :Compartment_getSpatialDimensions,
+            ),
+            get_optional_double(co, :Compartment_isSetSize, :Compartment_getSize),
+            get_optional_string(co, :Compartment_getUnits),
+            get_notes(co),
+            get_annotation(co),
+        )
+    end
 
     # parse out species
     species = Dict{String,Species}()
     for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
         sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
         sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t
+
         formula = nothing
         charge = nothing
         if sp_fbc != C_NULL
@@ -167,13 +219,33 @@ function extractModel(mdl::VPtr)::Model
                 charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
             end
         end
+
+        ia = nothing
+        if ccall(sbml(:Species_isSetInitialAmount), Cint, (VPtr,), sp) != 0
+            ia = (
+                ccall(sbml(:Species_getInitialAmount), Cdouble, (VPtr,), sp),
+                get_string(sp, :Species_getSubstanceUnits),
+            )
+        end
+
         species[get_string(sp, :Species_getId)] = Species(
             get_optional_string(sp, :Species_getName),
             get_string(sp, :Species_getCompartment),
+            get_optional_bool(
+                sp,
+                :Species_isSetBoundaryCondition,
+                :Species_getBoundaryCondition,
+            ),
             formula,
             charge,
-            getNotes(sp),
-            getAnnotation(sp),
+            ia,
+            get_optional_bool(
+                sp,
+                :Species_isSetHasOnlySubstanceUnits,
+                :Species_getHasOnlySubstanceUnits,
+            ),
+            get_notes(sp),
+            get_annotation(sp),
         )
     end
 
@@ -204,6 +276,7 @@ function extractModel(mdl::VPtr)::Model
         lb = (-Inf, "") # (bound value, unit id)
         ub = (Inf, "")
         oc = 0.0
+        math = nothing
 
         # kinetic laws store a second version of the bounds and objectives
         kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
@@ -221,6 +294,10 @@ function extractModel(mdl::VPtr)::Model
                     oc = pval()
                 end
             end
+
+            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
+                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
+            end
         end
 
         # TRICKY: SBML spec is completely silent about what should happen if
@@ -260,6 +337,7 @@ function extractModel(mdl::VPtr)::Model
             add_stoi(sr, 1)
         end
 
+        # gene product associations
         association = nothing
         if re_fbc != C_NULL
             gpa = ccall(
@@ -282,8 +360,9 @@ function extractModel(mdl::VPtr)::Model
             ub,
             haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
             association,
-            getNotes(re),
-            getAnnotation(re),
+            math,
+            get_notes(re),
+            get_annotation(re),
         )
     end
 
@@ -305,13 +384,30 @@ function extractModel(mdl::VPtr)::Model
                 gene_products[id] = GeneProduct(
                     get_optional_string(gp, :GeneProduct_getName),
                     get_optional_string(gp, :GeneProduct_getLabel),
-                    getNotes(gp),
-                    getAnnotation(gp),
+                    get_notes(gp),
+                    get_annotation(gp),
                 )
             end
         end
     end
 
+    function_definitions = Dict{String,FunctionDefinition}()
+    for i = 1:ccall(sbml(:Model_getNumFunctionDefinitions), Cuint, (VPtr,), mdl)
+        fd = ccall(sbml(:Model_getFunctionDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
+        def = nothing
+        if ccall(sbml(:FunctionDefinition_isSetMath), Cint, (VPtr,), fd) != 0
+            def = parse_math(ccall(sbml(:FunctionDefinition_getMath), VPtr, (VPtr,), fd))
+        end
+
+        function_definitions[get_string(fd, :FunctionDefinition_getId)] =
+            FunctionDefinition(
+                get_optional_string(fd, :FunctionDefinition_getName),
+                def,
+                get_notes(fd),
+                get_annotation(fd),
+            )
+    end
+
     return Model(
         parameters,
         units,
@@ -319,7 +415,8 @@ function extractModel(mdl::VPtr)::Model
         species,
         reactions,
         gene_products,
-        getNotes(mdl),
-        getAnnotation(mdl),
+        function_definitions,
+        get_notes(mdl),
+        get_annotation(mdl),
     )
 end

src/structs.jl

@@ -55,6 +55,54 @@ struct GPAOr <: GeneProductAssociation
     terms::Vector{GeneProductAssociation}
 end
 
+"""
+A simplified representation of MathML-specified math AST
+"""
+abstract type Math end
+
+"""
+A literal value (usually a numeric constant) in mathematical expression
+"""
+struct MathVal{T} <: Math where {T}
+    val::T
+end
+
+"""
+An identifier (usually a variable name) in mathematical expression
+"""
+struct MathIdent <: Math
+    id::String
+end
+
+"""
+Function application ("call by name", no tricks allowed) in mathematical expression
+"""
+struct MathApply <: Math
+    fn::String
+    args::Vector{Math}
+end
+
+"""
+Function definition (aka "lambda") in mathematical expression
+"""
+struct MathLambda <: Math
+    args::Vector{String}
+    body::Math
+end
+
+"""
+SBML Compartment with sizing information.
+"""
+struct Compartment
+    name::Maybe{String}
+    constant::Maybe{Bool}
+    spatial_dimensions::Maybe{Int}
+    size::Maybe{Float64}
+    units::Maybe{String}
+    notes::Maybe{String}
+    annotation::Maybe{String}
+    Compartment(na, c, sd, s, u, no = nothing, an = nothing) = new(na, c, sd, s, u, no, an)
+end
 
 """
 Reaction with stoichiometry that assigns reactants and products their relative
@@ -68,9 +116,10 @@ struct Reaction
     ub::Tuple{Float64,String}
     oc::Float64
     gene_product_association::Maybe{GeneProductAssociation}
+    kinetic_math::Maybe{Math}
     notes::Maybe{String}
     annotation::Maybe{String}
-    Reaction(s, l, u, o, as, n = nothing, an = nothing) = new(s, l, u, o, as, n, an)
+    Reaction(s, l, u, o, as, km, n = nothing, an = nothing) = new(s, l, u, o, as, km, n, an)
 end
 
 """
@@ -80,11 +129,15 @@ identifier, `formula`, `charge`, and additional `notes` and `annotation`.
 struct Species
     name::Maybe{String}
     compartment::String
+    boundary_condition::Maybe{Bool}
     formula::Maybe{String}
     charge::Maybe{Int}
+    initial_amount::Maybe{Tuple{Float64,String}}
+    only_substance_units::Maybe{Bool}
     notes::Maybe{String}
     annotation::Maybe{String}
-    Species(na, co, f, ch, no = nothing, a = nothing) = new(na, co, f, ch, no, a)
+    Species(na, co, b, f, ch, ia, osu, no = nothing, a = nothing) =
+        new(na, co, b, f, ch, ia, osu, no, a)
 end
 
 """
@@ -98,6 +151,17 @@ struct GeneProduct
     GeneProduct(na, l, no = nothing, a = nothing) = new(na, l, no, a)
 end
 
+"""
+Custom function definition.
+"""
+struct FunctionDefinition
+    name::Maybe{String}
+    body::Maybe{Math}
+    notes::Maybe{String}
+    annotation::Maybe{String}
+    FunctionDefinition(na, b, no = nothing, a = nothing) = new(na, b, no, a)
+end
+
 """
 Structure that collects the model-related data. Contains `parameters`, `units`,
 `compartments`, `species` and `reactions` and `gene_products`, and additional
@@ -108,11 +172,12 @@ objects.
 struct Model
     parameters::Dict{String,Float64}
     units::Dict{String,Vector{UnitPart}}
-    compartments::Vector{String}
+    compartments::Dict{String,Compartment}
     species::Dict{String,Species}
     reactions::Dict{String,Reaction}
     gene_products::Dict{String,GeneProduct}
+    function_definitions::Dict{String,FunctionDefinition}
     notes::Maybe{String}
     annotation::Maybe{String}
-    Model(p, u, c, s, r, g, n = nothing, a = nothing) = new(p, u, c, s, r, g, n, a)
+    Model(p, u, c, s, r, g, f, n = nothing, a = nothing) = new(p, u, c, s, r, g, f, n, a)
 end