readsbml.jl 15.6 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
    function readSBML(fn::String, sbml_conversion = model->nothing)::SBML.Model
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
76

anand jain's avatar
anand jain committed
77
Read the SBML from a XML file in `fn` and return the contained `SBML.Model`.
78
79
80
81

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
82
[`set_level_and_version`](@ref), [`libsbml_convert`](@ref), and
83
84
85
86
87
[`convert_simplify_math`](@ref).

# Example
```
m = readSBML("my_model.xml", doc -> begin
88
    set_level_and_version(3, 1)(doc)
89
90
91
    convert_simplify_math(doc)
end)
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
92
"""
93
function readSBML(fn::String, sbml_conversion = document -> nothing)::SBML.Model
94
95
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
96
        get_error_messages(doc, AssertionError("Opening SBML document has reported errors"))
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
97

98
99
        sbml_conversion(doc)

100
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
101
            throw(AssertionError("SBML document contains no model"))
102
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
103

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

106
107
108
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
109
    end
110
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
111

112
113
get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
114

115
"""
116
117
118
119
120
121
122
123
124
125
    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
126
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
    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
146
    function extractModel(mdl::VPtr)::SBML.Model
147
148

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

155
    # get the parameters
St. Elmo's avatar
St. Elmo committed
156
157
158
    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)
159
        id = get_string(p, :Parameter_getId)
St. Elmo's avatar
St. Elmo committed
160
161
162
163
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

164
    # parse out the unit definitions
anand jain's avatar
anand jain committed
165
    units = Dict{String,Vector{SBML.UnitPart}}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
166
167
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
168
        id = get_string(ud, :UnitDefinition_getId)
169
170
        units[id] = [
            begin
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
171
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
anand jain's avatar
anand jain committed
172
                SBML.UnitPart(
173
174
175
176
177
178
179
180
181
182
183
184
                    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
185
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
186
        ]
187
188
    end

189
    # parse out compartment names
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    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
208

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

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

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

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

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

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

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

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

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

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

335
        # extract stoichiometry
336
        stoi = Dict{String,Float64}()
337
338
339
340
341
342
343
        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
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