readsbml.jl 15.2 KB
Newer Older
1
"""
2
    get_string(x::VPtr, fn_sym)::Maybe{String}
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

C-call the SBML function `fn_sym` with a single parameter `x`, interpret the
result as a string and return it, or throw exception in case the pointer is
NULL.
"""
function get_string(x::VPtr, fn_sym)::String
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        throw(DomainError(x, "Calling $fn_sym returned NULL, valid string expected."))
    end
end

"""
18
    get_optional_string(x::VPtr, fn_sym)::Maybe{String}
19

20
Like [`get_string`](@ref), but returns `nothing` instead of throwing an
21
22
23
exception.

This is used to get notes and annotations and several other things (see
24
`get_notes`, `get_annotations`)
25
26
27
28
29
30
31
32
33
34
"""
function get_optional_string(x::VPtr, fn_sym)::Maybe{String}
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
    if str != C_NULL
        return unsafe_string(str)
    else
        return nothing
    end
end

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
"""
    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

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
74
"""
75
76
77
78
79
    readSBML(
        fn::String,
        sbml_conversion = document -> nothing;
        report_severities = ["Fatal", "Error"],
    )::SBML.Model
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
80

anand jain's avatar
anand jain committed
81
Read the SBML from a XML file in `fn` and return the contained `SBML.Model`.
82
83
84
85

The `sbml_conversion` is a function that does an in-place modification of the
single parameter, which is the C pointer to the loaded SBML document (C type
`SBMLDocument*`). Several functions for doing that are prepared, including
86
[`set_level_and_version`](@ref), [`libsbml_convert`](@ref), and
87
88
[`convert_simplify_math`](@ref).

89
90
91
`report_severities` switches on and off reporting of certain errors; see the
documentation of [`get_error_messages`](@ref) for details.

92
93
94
# Example
```
m = readSBML("my_model.xml", doc -> begin
95
    set_level_and_version(3, 1)(doc)
96
97
98
    convert_simplify_math(doc)
end)
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
99
"""
100
101
102
103
104
function readSBML(
    fn::String,
    sbml_conversion = document -> nothing;
    report_severities = ["Fatal", "Error"],
)::SBML.Model
105
106
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
107
108
109
110
111
        get_error_messages(
            doc,
            AssertionError("Opening SBML document has reported errors"),
            report_severities,
        )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
112

113
114
        sbml_conversion(doc)

115
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
116
            throw(AssertionError("SBML document contains no model"))
117
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
118

119
        model = ccall(sbml(:SBMLDocument_getModel), VPtr, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
120

121
122
123
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
124
    end
125
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
126

127
128
get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
129

130
"""
131
132
133
134
135
136
137
138
139
140
    function getAssociation(x::VPtr)::GeneProductAssociation

Convert a pointer to SBML `FbcAssociation_t` to the `GeneProductAssociation`
tree structure.
"""
function getAssociation(x::VPtr)::GeneProductAssociation
    # libsbml C API is currently missing functions to check this in a normal
    # way, so we use a bit of a hack.
    typecode = ccall(sbml(:SBase_getTypeCode), Cint, (VPtr,), x)
    if typecode == 808 # SBML_FBC_GENEPRODUCTREF
141
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    elseif typecode == 809 # SBML_FBC_AND
        return GPAAnd([
            getAssociation(
                ccall(sbml(:FbcAnd_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcAnd_getNumAssociations), Cuint, (VPtr,), x)
        ])
    elseif typecode == 810 # SBML_FBC_OR
        return GPAOr([
            getAssociation(
                ccall(sbml(:FbcOr_getAssociation), VPtr, (VPtr, Cuint), x, i - 1),
            ) for i = 1:ccall(sbml(:FbcOr_getNumAssociations), Cuint, (VPtr,), x)
        ])
    else
        throw(ErrorException("Unsupported FbcAssociation type"))
    end
end


""""
anand jain's avatar
anand jain committed
161
    function extractModel(mdl::VPtr)::SBML.Model
162
163

Take the `SBMLModel_t` pointer and extract all information required to make a
anand jain's avatar
anand jain committed
164
valid [`SBML.Model`](@ref) structure.
165
"""
anand jain's avatar
anand jain committed
166
function extractModel(mdl::VPtr)::SBML.Model
167
    # get the FBC plugin pointer (FbcModelPlugin_t)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
168
169
    mdl_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), mdl, "fbc")

170
    # get the parameters
St. Elmo's avatar
St. Elmo committed
171
172
173
    parameters = Dict{String,Float64}()
    for i = 1:ccall(sbml(:Model_getNumParameters), Cuint, (VPtr,), mdl)
        p = ccall(sbml(:Model_getParameter), VPtr, (VPtr, Cuint), mdl, i - 1)
174
        id = get_string(p, :Parameter_getId)
St. Elmo's avatar
St. Elmo committed
175
176
177
178
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

179
    # parse out the unit definitions
180
    units = Dict{String,Quantity}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
181
182
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
183
        id = get_string(ud, :UnitDefinition_getId)
184
        units[id] = get_units(ud)
185
    end
186
    # parse out compartment names
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
    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
205

206
    # parse out species
207
    species = Dict{String,Species}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
208
209
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
210
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t
211

212
        formula = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
213
214
        charge = nothing
        if sp_fbc != C_NULL
215
            # if the FBC plugin is present, try to get the chemical formula and charge
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
216
217
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
218
                formula = get_string(sp_fbc, :FbcSpeciesPlugin_getChemicalFormula)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
219
220
221
222
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
223
        end
224
225
226
227
228

        ia = nothing
        if ccall(sbml(:Species_isSetInitialAmount), Cint, (VPtr,), sp) != 0
            ia = (
                ccall(sbml(:Species_getInitialAmount), Cdouble, (VPtr,), sp),
229
230
231
232
233
234
235
                get_optional_string(sp, :Species_getSubstanceUnits),
            )
        end

        ic = nothing
        if ccall(sbml(:Species_isSetInitialConcentration), Cint, (VPtr,), sp) != 0
            ic = (
paulflang's avatar
paulflang committed
236
                ccall(sbml(:Species_getInitialConcentration), Cdouble, (VPtr,), sp),
237
                get_optional_string(sp, :Species_getSubstanceUnits),
238
239
240
            )
        end

241
        species[get_string(sp, :Species_getId)] = Species(
242
            get_optional_string(sp, :Species_getName),
243
            get_string(sp, :Species_getCompartment),
244
245
246
247
248
            get_optional_bool(
                sp,
                :Species_isSetBoundaryCondition,
                :Species_getBoundaryCondition,
            ),
249
            formula,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
250
            charge,
251
            ia,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
252
            ic,
253
254
255
256
257
            get_optional_bool(
                sp,
                :Species_isSetHasOnlySubstanceUnits,
                :Species_getHasOnlySubstanceUnits,
            ),
258
259
            get_notes(sp),
            get_annotation(sp),
260
        )
261
262
    end

263
264
    # parse out the flux objectives (these are complementary to the objectives
    # that appear in the reactions, see comments lower)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
265
266
267
268
269
270
271
272
273
274
275
276
    objectives_fbc = Dict{String,Float64}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumObjectives), Cuint, (VPtr,), mdl_fbc)
            o = ccall(
                sbml(:FbcModelPlugin_getObjective),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )
            for j = 1:ccall(sbml(:Objective_getNumFluxObjectives), Cuint, (VPtr,), o)
                fo = ccall(sbml(:Objective_getFluxObjective), VPtr, (VPtr, Cuint), o, j - 1)
277
278
                objectives_fbc[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
279
280
281
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
282

283
    # reactions!
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
284
    reactions = Dict{String,Reaction}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
285
286
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
287
        lb = (-Inf, "") # (bound value, unit id)
288
289
        ub = (Inf, "")
        oc = 0.0
290
        math = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
291

292
        # kinetic laws store a second version of the bounds and objectives
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
293
294
295
296
        kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
        if kl != C_NULL
            for j = 1:ccall(sbml(:KineticLaw_getNumParameters), Cuint, (VPtr,), kl)
                p = ccall(sbml(:KineticLaw_getParameter), VPtr, (VPtr, Cuint), kl, j - 1)
297
                id = get_string(p, :Parameter_getId)
298
299
                pval() = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
                punit() = get_string(p, :Parameter_getUnits)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
300
301
302
303
304
305
306
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
307
            end
308
309

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
310
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
311
            end
312
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
313

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
314
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
315
316
317
318
319
320
321
        # someone specifies both the above and below formats of the flux bounds
        # for one reaction. Notably, these do not really specify much
        # interaction with units. In this case, we'll just set a special
        # "[fbc]" unit that has no specification in `units`, and hope the users
        # can make something out of it.
        re_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), re, "fbc")
        if re_fbc != C_NULL
322
323
324
            fbcb = get_optional_string(re_fbc, :FbcReactionPlugin_getLowerFluxBound)
            if !isnothing(fbcb) && haskey(parameters, fbcb)
                lb = (parameters[fbcb], "[fbc]")
St. Elmo's avatar
St. Elmo committed
325
            end
326
327
328
            fbcb = get_optional_string(re_fbc, :FbcReactionPlugin_getUpperFluxBound)
            if !isnothing(fbcb) && haskey(parameters, fbcb)
                ub = (parameters[fbcb], "[fbc]")
St. Elmo's avatar
St. Elmo committed
329
330
331
            end
        end

332
        # extract stoichiometry
333
        stoi = Dict{String,Float64}()
334
335
336
337
338
339
340
        add_stoi(sr, factor) = begin
            s = get_string(sr, :SpeciesReference_getSpecies)
            stoi[s] =
                get(stoi, s, 0) +
                ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                factor
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
341

342
        # reactants and products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
343
344
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
345
346
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
347

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
348
349
        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
350
351
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
352

353
        # gene product associations
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
        association = nothing
        if re_fbc != C_NULL
            gpa = ccall(
                sbml(:FbcReactionPlugin_getGeneProductAssociation),
                VPtr,
                (VPtr,),
                re_fbc,
            )
            if gpa != C_NULL
                a = ccall(sbml(:GeneProductAssociation_getAssociation), VPtr, (VPtr,), gpa)
                a != C_NULL
                association = getAssociation(a)
            end
        end

paulflang's avatar
paulflang committed
369
370
371
        # explicit reversible flag (defaults to true in SBML)
        reversible = Bool(ccall(sbml(:Reaction_getReversible), Cint, (VPtr,), re))

372
        reid = get_string(re, :Reaction_getId)
373
374
375
376
377
        reactions[reid] = Reaction(
            stoi,
            lb,
            ub,
            haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
378
            association,
379
            math,
paulflang's avatar
paulflang committed
380
            reversible,
381
382
            get_notes(re),
            get_annotation(re),
383
        )
384
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
385

386
    # extract gene products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
387
388
389
390
391
392
393
394
395
396
397
    gene_products = Dict{String,GeneProduct}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumGeneProducts), Cuint, (VPtr,), mdl_fbc)
            gp = ccall(
                sbml(:FbcModelPlugin_getGeneProduct),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )

398
            id = get_optional_string(gp, :GeneProduct_getId) # IDs don't need to be set
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
399
400
401

            if id != nothing
                gene_products[id] = GeneProduct(
402
403
                    get_optional_string(gp, :GeneProduct_getName),
                    get_optional_string(gp, :GeneProduct_getLabel),
404
405
                    get_notes(gp),
                    get_annotation(gp),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
406
407
408
409
410
                )
            end
        end
    end

411
412
413
414
415
416
417
418
419
420
421
422
    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,
423
424
                get_notes(fd),
                get_annotation(fd),
425
426
427
            )
    end

428
429
430
431
432
433
    return Model(
        parameters,
        units,
        compartments,
        species,
        reactions,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
434
        gene_products,
435
436
437
        function_definitions,
        get_notes(mdl),
        get_annotation(mdl),
438
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
439
end