readsbml.jl 15.6 KB
Newer Older
1

2
const VPtr = Ptr{Cvoid}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3

4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
"""
    function get_string(x::VPtr, fn_sym)::Maybe{String}

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

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

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
27
`get_notes`, `get_annotations`)
28
29
30
31
32
33
34
35
36
37
"""
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

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
74
75
76
"""
    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
77
"""
78
    function readSBML(fn::String, sbml_conversion = model->nothing)::SBML.Model
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
79

anand jain's avatar
anand jain committed
80
Read the SBML from a XML file in `fn` and return the contained `SBML.Model`.
81
82
83
84
85
86
87
88
89
90
91
92
93
94

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
[`convert_level_and_version`](@ref), [`libsbml_convert`](@ref), and
[`convert_simplify_math`](@ref).

# Example
```
m = readSBML("my_model.xml", doc -> begin
    convert_level_and_version(3, 1)(doc)
    convert_simplify_math(doc)
end)
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
95
"""
96
function readSBML(fn::String, sbml_conversion = document -> nothing)::SBML.Model
97
98
99
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
100
101
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
102
            msg = get_string(err, :XMLError_getMessage)
103
104
            @warn "SBML reported error: $msg"
        end
105
        if n_errs > 0
106
            throw(AssertionError("Opening SBML document has reported errors"))
107
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
108

109
110
        sbml_conversion(doc)

111
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
112
            throw(AssertionError("SBML document contains no model"))
113
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
114

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

117
118
119
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
120
    end
121
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
122

123
124
get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
125

126
"""
127
128
129
130
131
132
133
134
135
136
    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
137
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    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
157
    function extractModel(mdl::VPtr)::SBML.Model
158
159

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

166
    # get the parameters
St. Elmo's avatar
St. Elmo committed
167
168
169
    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)
170
        id = get_string(p, :Parameter_getId)
St. Elmo's avatar
St. Elmo committed
171
172
173
174
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

175
    # parse out the unit definitions
anand jain's avatar
anand jain committed
176
    units = Dict{String,Vector{SBML.UnitPart}}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
177
178
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
179
        id = get_string(ud, :UnitDefinition_getId)
180
181
        units[id] = [
            begin
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
182
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
anand jain's avatar
anand jain committed
183
                SBML.UnitPart(
184
185
186
187
188
189
190
191
192
193
194
195
                    unsafe_string(
                        ccall(
                            sbml(:UnitKind_toString),
                            Cstring,
                            (Cint,),
                            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
                        ),
                    ),
                    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getScale), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u),
                )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
196
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
197
        ]
198
199
    end

200
    # parse out compartment names
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    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
219

220
    # parse out species
221
    species = Dict{String,Species}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
222
223
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
224
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t
225

226
        formula = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
227
228
        charge = nothing
        if sp_fbc != C_NULL
229
            # if the FBC plugin is present, try to get the chemical formula and charge
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
230
231
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
232
                formula = get_string(sp_fbc, :FbcSpeciesPlugin_getChemicalFormula)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
233
234
235
236
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
237
        end
238
239
240
241
242
243
244
245
246

        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

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

268
269
    # 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
270
271
272
273
274
275
276
277
278
279
280
281
    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)
282
283
                objectives_fbc[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
284
285
286
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
287

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

297
        # kinetic laws store a second version of the bounds and objectives
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
298
299
300
301
        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)
302
                id = get_string(p, :Parameter_getId)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
303
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
304
                punit = () -> get_string(p, :Parameter_getUnits)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
305
306
307
308
309
310
311
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
312
            end
313
314

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
315
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
316
            end
317
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
318

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
319
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
320
321
322
323
324
325
326
        # 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
327
328
329
            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
330
            end
331
332
333
            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
334
335
336
            end
        end

337
        # extract stoichiometry
338
339
340
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
341
                stoi[get_string(sr, :SpeciesReference_getSpecies)] =
342
343
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
344

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

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

356
        # gene product associations
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
        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
372
373
374
        # explicit reversible flag (defaults to true in SBML)
        reversible = Bool(ccall(sbml(:Reaction_getReversible), Cint, (VPtr,), re))

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

389
    # extract gene products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
390
391
392
393
394
395
396
397
398
399
400
    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,
            )

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

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

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

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