readsbml.jl 12.4 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
27
28
29
30
31
32
33
34
35
36
37
"""
    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
`getNotes`, `getAnnotations`)
"""
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

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
38
39
40
41
42
"""
    function readSBML(fn::String)::Model

Read the SBML from a XML file in `fn` and return the contained `Model`.
"""
43
44
45
46
function readSBML(fn::String)::Model
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
47
48
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
49
            msg = get_string(err, :XMLError_getMessage)
50
51
            @warn "SBML reported error: $msg"
        end
52
        if n_errs > 0
53
            throw(AssertionError("Opening SBML document has reported errors"))
54
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
55

56
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
57
            throw(AssertionError("SBML document contains no model"))
58
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
59

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

62
63
64
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
65
    end
66
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
67

68
69
getNotes(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getNotesString)
getAnnotation(x::VPtr)::Maybe{String} = get_optional_string(x, :SBase_getAnnotationString)
70

71
"""
72
73
74
75
76
77
78
79
80
81
    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
82
        return GPARef(get_string(x, :GeneProductRef_getGeneProduct))
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
    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


""""
102
103
104
105
106
    function extractModel(mdl::VPtr)::Model

Take the `SBMLModel_t` pointer and extract all information required to make a
valid [`Model`](@ref) structure.
"""
107
function extractModel(mdl::VPtr)::Model
108
    # get the FBC plugin pointer (FbcModelPlugin_t)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
109
110
    mdl_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), mdl, "fbc")

111
    # get the parameters
St. Elmo's avatar
St. Elmo committed
112
113
114
    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)
115
        id = get_string(p, :Parameter_getId)
St. Elmo's avatar
St. Elmo committed
116
117
118
119
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

120
    # parse out the unit definitions
121
    units = Dict{String,Vector{UnitPart}}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
122
123
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
124
        id = get_string(ud, :UnitDefinition_getId)
125
126
        units[id] = [
            begin
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
127
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
128
129
130
131
132
133
134
135
136
137
138
139
140
                UnitPart(
                    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
141
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
142
        ]
143
144
    end

145
    # parse out compartment names
146
    compartments = [
147
148
149
        get_string(
            ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1),
            :Compartment_getId,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
150
        ) for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
151
152
    ]

153
    # parse out species
154
    species = Dict{String,Species}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
155
156
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
157
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t
158

159
        formula = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
160
161
        charge = nothing
        if sp_fbc != C_NULL
162
            # if the FBC plugin is present, try to get the chemical formula and charge
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
163
164
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
165
                formula = get_string(sp_fbc, :FbcSpeciesPlugin_getChemicalFormula)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
166
167
168
169
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
170
        end
171
172
173
174
175
176
177
178
179

        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

180
        species[get_string(sp, :Species_getId)] = Species(
181
            get_optional_string(sp, :Species_getName),
182
            get_string(sp, :Species_getCompartment),
183
            formula,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
184
            charge,
185
            ia,
186
            ccall(sbml(:Species_getHasOnlySubstanceUnits), Cint, (VPtr,), sp) != 0,
187
            getNotes(sp),
188
            getAnnotation(sp),
189
        )
190
191
    end

192
193
    # 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
194
195
196
197
198
199
200
201
202
203
204
205
    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)
206
207
                objectives_fbc[get_string(fo, :FluxObjective_getReaction)] =
                    ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
208
209
210
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
211

212
    # reactions!
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
213
    reactions = Dict{String,Reaction}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
214
215
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
216
        lb = (-Inf, "") # (bound value, unit id)
217
218
        ub = (Inf, "")
        oc = 0.0
219
        math = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
220

221
        # kinetic laws store a second version of the bounds and objectives
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
222
223
224
225
        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)
226
                id = get_string(p, :Parameter_getId)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
227
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
228
                punit = () -> get_string(p, :Parameter_getUnits)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
229
230
231
232
233
234
235
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
236
            end
237
238

            if ccall(sbml(:KineticLaw_isSetMath), Cint, (VPtr,), kl) != 0
239
                math = parse_math(ccall(sbml(:KineticLaw_getMath), VPtr, (VPtr,), kl))
240
            end
241
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
242

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
243
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
244
245
246
247
248
249
250
        # 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
251
252
253
            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
254
            end
255
256
257
            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
258
259
260
            end
        end

261
        # extract stoichiometry
262
263
264
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
265
                stoi[get_string(sr, :SpeciesReference_getSpecies)] =
266
267
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
268

269
        # reactants and products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
270
271
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
272
273
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
274

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
275
276
        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
277
278
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
279

280
        # gene product associations
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
        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

296
        reid = get_string(re, :Reaction_getId)
297
298
299
300
301
        reactions[reid] = Reaction(
            stoi,
            lb,
            ub,
            haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
302
            association,
303
            math,
304
            getNotes(re),
305
            getAnnotation(re),
306
        )
307
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
308

309
    # extract gene products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
310
311
312
313
314
315
316
317
318
319
320
    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,
            )

321
            id = get_optional_string(gp, :GeneProduct_getId) # IDs don't need to be set
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
322
323
324

            if id != nothing
                gene_products[id] = GeneProduct(
325
326
                    get_optional_string(gp, :GeneProduct_getName),
                    get_optional_string(gp, :GeneProduct_getLabel),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
327
328
329
330
331
332
333
                    getNotes(gp),
                    getAnnotation(gp),
                )
            end
        end
    end

334
335
336
337
338
339
    return Model(
        parameters,
        units,
        compartments,
        species,
        reactions,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
340
        gene_products,
341
342
343
        getNotes(mdl),
        getAnnotation(mdl),
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
344
end