readsbml.jl 12.2 KB
Newer Older
1

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

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
4
5
6
7
8
"""
    function readSBML(fn::String)::Model

Read the SBML from a XML file in `fn` and return the contained `Model`.
"""
9
10
11
12
function readSBML(fn::String)::Model
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
13
14
15
16
17
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
            msg = unsafe_string(ccall(sbml(:XMLError_getMessage), Cstring, (VPtr,), err))
            @warn "SBML reported error: $msg"
        end
18
        if n_errs > 0
19
            throw(AssertionError("Opening SBML document has reported errors"))
20
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
21

22
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
23
            throw(AssertionError("SBML document contains no model"))
24
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
25

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

28
29
30
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31
    end
32
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
33

34
35
36
37
38
39
40
"""
    function getOptionalString(x::VPtr, fn_sym)::Maybe{String}

C-call the SBML function `fn_sym` with a single parameter `x`, interpret the result as a nullable string pointer and return appropriately.

This is used to get notes and annotations and several other things (see `getNotes`, `getAnnotations`)
"""
41
42
function getOptionalString(x::VPtr, fn_sym)::Maybe{String}
    str = ccall(sbml(fn_sym), Cstring, (VPtr,), x)
43
44
45
46
47
48
49
    if str != C_NULL
        return unsafe_string(str)
    else
        return nothing
    end
end

50
51
52
getNotes(x::VPtr)::Maybe{String} = getOptionalString(x, :SBase_getNotesString)
getAnnotation(x::VPtr)::Maybe{String} = getOptionalString(x, :SBase_getAnnotationString)

53
"""
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
    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
        return GPARef(
            unsafe_string(ccall(sbml(:GeneProductRef_getGeneProduct), Cstring, (VPtr,), x)),
        )
    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


""""
86
87
88
89
90
    function extractModel(mdl::VPtr)::Model

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

95
    # get the parameters
St. Elmo's avatar
St. Elmo committed
96
97
98
99
100
101
102
103
    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)
        id = unsafe_string(ccall(sbml(:Parameter_getId), Cstring, (VPtr,), p))
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

104
    # parse out the unit definitions
105
    units = Dict{String,Vector{UnitPart}}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
106
107
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
108
109
110
        id = unsafe_string(ccall(sbml(:UnitDefinition_getId), Cstring, (VPtr,), ud))
        units[id] = [
            begin
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
111
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
112
113
114
115
116
117
118
119
120
121
122
123
124
                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
125
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
126
        ]
127
128
    end

129
    # parse out compartment names
130
131
132
133
134
135
    compartments = [
        unsafe_string(
            ccall(
                sbml(:Compartment_getId),
                Cstring,
                (VPtr,),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
136
                ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1),
137
            ),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
138
        ) for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
139
140
    ]

141
    # parse out species
142
    species = Dict{String,Species}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
143
144
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
145
        sp_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), sp, "fbc") # FbcSpeciesPlugin_t
146
        formula = nothing
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
147
148
        charge = nothing
        if sp_fbc != C_NULL
149
            # if the FBC plugin is present, try to get the chemical formula and charge
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
150
151
152
153
154
155
156
157
158
159
160
161
162
163
            if 0 !=
               ccall(sbml(:FbcSpeciesPlugin_isSetChemicalFormula), Cint, (VPtr,), sp_fbc)
                formula = unsafe_string(
                    ccall(
                        sbml(:FbcSpeciesPlugin_getChemicalFormula),
                        Cstring,
                        (VPtr,),
                        sp_fbc,
                    ),
                )
            end
            if 0 != ccall(sbml(:FbcSpeciesPlugin_isSetCharge), Cint, (VPtr,), sp_fbc)
                charge = ccall(sbml(:FbcSpeciesPlugin_getCharge), Cint, (VPtr,), sp_fbc)
            end
164
        end
165
166
167
        species[unsafe_string(ccall(sbml(:Species_getId), Cstring, (VPtr,), sp))] = Species(
            unsafe_string(ccall(sbml(:Species_getName), Cstring, (VPtr,), sp)),
            unsafe_string(ccall(sbml(:Species_getCompartment), Cstring, (VPtr,), sp)),
168
            formula,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
169
            charge,
170
            getNotes(sp),
171
            getAnnotation(sp),
172
        )
173
174
    end

175
176
    # 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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    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)
                objectives_fbc[unsafe_string(
                    ccall(sbml(:FluxObjective_getReaction), Cstring, (VPtr,), fo),
                )] = ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
195

196
    # reactions!
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
197
    reactions = Dict{String,Reaction}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
198
199
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
200
        lb = (-Inf, "") # (bound value, unit id)
201
202
        ub = (Inf, "")
        oc = 0.0
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
203

204
        # kinetic laws store a second version of the bounds and objectives
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
        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)
                id = unsafe_string(ccall(sbml(:Parameter_getId), Cstring, (VPtr,), p))
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
                punit =
                    () ->
                        unsafe_string(ccall(sbml(:Parameter_getUnits), Cstring, (VPtr,), p))
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
221
            end
222
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
223

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
224
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
        # 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
            fbcb =
                ccall(sbml(:FbcReactionPlugin_getLowerFluxBound), Cstring, (VPtr,), re_fbc)
            if fbcb != C_NULL && haskey(parameters, unsafe_string(fbcb))
                lb = (parameters[unsafe_string(fbcb)], "[fbc]")
            end
            fbcb =
                ccall(sbml(:FbcReactionPlugin_getUpperFluxBound), Cstring, (VPtr,), re_fbc)
            if fbcb != C_NULL && haskey(parameters, unsafe_string(fbcb))
                ub = (parameters[unsafe_string(fbcb)], "[fbc]")
            end
        end

244
        # extract stoichiometry
245
246
247
248
249
250
251
252
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
                stoi[unsafe_string(
                    ccall(sbml(:SpeciesReference_getSpecies), Cstring, (VPtr,), sr),
                )] =
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
253

254
        # reactants and products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
255
256
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
257
258
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
259

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
260
261
        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
262
263
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
264

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
        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)
                @info "assoc:" association
            end
        end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
281
        reid = unsafe_string(ccall(sbml(:Reaction_getId), Cstring, (VPtr,), re))
282
283
284
285
286
287
        reactions[reid] = Reaction(
            stoi,
            lb,
            ub,
            haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc,
            getNotes(re),
288
            getAnnotation(re),
289
        )
290
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
291

292
    # extract gene products
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    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,
            )

            id = getOptionalString(gp, :GeneProduct_getId) # IDs don't need to be set

            if id != nothing
                gene_products[id] = GeneProduct(
                    getOptionalString(gp, :GeneProduct_getName),
                    getOptionalString(gp, :GeneProduct_getLabel),
                    getNotes(gp),
                    getAnnotation(gp),
                )
            end
        end
    end

317
318
319
320
321
322
    return Model(
        parameters,
        units,
        compartments,
        species,
        reactions,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
323
        gene_products,
324
325
326
        getNotes(mdl),
        getAnnotation(mdl),
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
327
end