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
96
97
        get_error_messages(doc,
                           AssertionError("Opening SBML document has reported errors"))
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
98

99
100
        sbml_conversion(doc)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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