readsbml.jl 15 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`.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
81
"""
82
function readSBML(fn::String, sbml_conversion = document -> nothing)::SBML.Model
83
84
85
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
86
87
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
88
            msg = get_string(err, :XMLError_getMessage)
89
90
            @warn "SBML reported error: $msg"
        end
91
        if n_errs > 0
92
            throw(AssertionError("Opening SBML document has reported errors"))
93
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
94

95
96
        sbml_conversion(doc)

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

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

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

109
110
get_notes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
get_annotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
111

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

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

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

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

        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

233
        species[get_string(sp, :Species_getId)] = Species(
234
            get_optional_string(sp, :Species_getName),
235
            get_string(sp, :Species_getCompartment),
236
237
238
239
240
            get_optional_bool(
                sp,
                :Species_isSetBoundaryCondition,
                :Species_getBoundaryCondition,
            ),
241
            formula,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
242
            charge,
243
            ia,
244
245
246
247
248
            get_optional_bool(
                sp,
                :Species_isSetHasOnlySubstanceUnits,
                :Species_getHasOnlySubstanceUnits,
            ),
249
250
            get_notes(sp),
            get_annotation(sp),
251
        )
252
253
    end

254
255
    # 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
256
257
258
259
260
261
262
263
264
265
266
267
    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)
268
269
                objectives_fbc[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
270
271
272
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
273

274
    # reactions!
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
275
    reactions = Dict{String,Reaction}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
276
277
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
278
        lb = (-Inf, "") # (bound value, unit id)
279
280
        ub = (Inf, "")
        oc = 0.0
281
        math = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
282

283
        # kinetic laws store a second version of the bounds and objectives
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
284
285
286
287
        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)
288
                id = get_string(p, :Parameter_getId)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
289
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
290
                punit = () -> get_string(p, :Parameter_getUnits)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
291
292
293
294
295
296
297
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
298
            end
299
300

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
301
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
302
            end
303
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
304

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
305
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
306
307
308
309
310
311
312
        # 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
313
314
315
            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
316
            end
317
318
319
            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
320
321
322
            end
        end

323
        # extract stoichiometry
324
325
326
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
327
                stoi[get_string(sr, :SpeciesReference_getSpecies)] =
328
329
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
330

331
        # reactants and products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
332
333
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
334
335
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
336

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
337
338
        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
339
340
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
341

342
        # gene product associations
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
        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

358
        reid = get_string(re, :Reaction_getId)
359
360
361
362
363
        reactions[reid] = Reaction(
            stoi,
            lb,
            ub,
            haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
364
            association,
365
            math,
366
367
            get_notes(re),
            get_annotation(re),
368
        )
369
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
370

371
    # extract gene products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
372
373
374
375
376
377
378
379
380
381
382
    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,
            )

383
            id = get_optional_string(gp, :GeneProduct_getId) # IDs don't need to be set
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
384
385
386

            if id != nothing
                gene_products[id] = GeneProduct(
387
388
                    get_optional_string(gp, :GeneProduct_getName),
                    get_optional_string(gp, :GeneProduct_getLabel),
389
390
                    get_notes(gp),
                    get_annotation(gp),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
391
392
393
394
395
                )
            end
        end
    end

396
397
398
399
400
401
402
403
404
405
406
407
    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,
408
409
                get_notes(fd),
                get_annotation(fd),
410
411
412
            )
    end

413
414
415
416
417
418
    return Model(
        parameters,
        units,
        compartments,
        species,
        reactions,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
419
        gene_products,
420
421
422
        function_definitions,
        get_notes(mdl),
        get_annotation(mdl),
423
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
424
end