StandardModel.jl 12 KB
Newer Older
St. Elmo's avatar
St. Elmo committed
1
"""
St. Elmo's avatar
St. Elmo committed
2
3
    mutable struct StandardModel

St. Elmo's avatar
St. Elmo committed
4
5
6
`StandardModel` is used to store a constraint based metabolic model with meta-information.
Meta-information is defined as annotation details, which include gene-reaction-rules, formulas, etc.

7
8
This model type seeks to keep as much meta-information as possible, as opposed to `CoreModel` and `CoreModelCoupled`,
which keep the bare neccessities only.
St. Elmo's avatar
St. Elmo committed
9
When merging models and keeping meta-information is important, use this as the model type. 
St. Elmo's avatar
St. Elmo committed
10
11
12
If meta-information is not important, use the more efficient core model types. 
See [`CoreModel`](@ref) and [`CoreModelCoupled`](@ref) for comparison.

St. Elmo's avatar
St. Elmo committed
13
14
In this model, reactions, metabolites, and genes are stored in ordered dictionaries indexed by each struct's `id` field.
For example, `model.reactions["rxn1_id"]` returns a `Reaction` where the field `id` equals `"rxn1_id"`.
St. Elmo's avatar
St. Elmo committed
15
16
This makes adding and removing reactions efficient.   

St. Elmo's avatar
St. Elmo committed
17
18
Note that the stoichiometric matrix (or any other core data, e.g. flux bounds) is not stored directly as in `CoreModel`. 
When this model type is used in analysis functions, these core data structures are built from scratch each time an analysis function is called.
St. Elmo's avatar
St. Elmo committed
19
20
21
22
This can cause performance issues if you run many small analysis functions sequentially. 
Consider using the core model types if performance is critical.

See also: [`Reaction`](@ref), [`Metabolite`](@ref), [`Gene`](@ref)
St. Elmo's avatar
St. Elmo committed
23
24

# Fields
St. Elmo's avatar
St. Elmo committed
25
```
St. Elmo's avatar
St. Elmo committed
26
id :: String
St. Elmo's avatar
St. Elmo committed
27
28
29
reactions :: OrderedDict{String, Reaction}
metabolites :: OrderedDict{String, Metabolite}
genes :: OrderedDict{String, Gene}
St. Elmo's avatar
St. Elmo committed
30
```
St. Elmo's avatar
St. Elmo committed
31
32
33
34
35

# Example
```
model = load_model(StandardModel, "model_location")
```
St. Elmo's avatar
St. Elmo committed
36
"""
St. Elmo's avatar
St. Elmo committed
37
mutable struct StandardModel <: MetabolicModel
St. Elmo's avatar
St. Elmo committed
38
    id::String
St. Elmo's avatar
St. Elmo committed
39
40
41
    reactions::OrderedDict{String,Reaction}
    metabolites::OrderedDict{String,Metabolite}
    genes::OrderedDict{String,Gene}
42
43

    StandardModel(
St. Elmo's avatar
St. Elmo committed
44
        id = "";
St. Elmo's avatar
St. Elmo committed
45
46
47
        reactions = OrderedDict{String,Reaction}(),
        metabolites = OrderedDict{String,Metabolite}(),
        genes = OrderedDict{String,Gene}(),
48
    ) = new(id, reactions, metabolites, genes)
St. Elmo's avatar
St. Elmo committed
49
50
end

51
# MetabolicModel interface follows
52
"""
St. Elmo's avatar
St. Elmo committed
53
    reactions(model::StandardModel)
54
55
56
57
58

Return a vector of reaction id strings contained in `model`.
The order of reaction ids returned here matches the order used to construct the
stoichiometric matrix.
"""
59
reactions(model::StandardModel)::StringVecType = collect(keys(model.reactions))
60
61
62
63
64
65

"""
    n_reactions(model::StandardModel)

Return the number of reactions contained in `model`.
"""
66
n_reactions(model::StandardModel)::Int = length(model.reactions)
St. Elmo's avatar
St. Elmo committed
67

68

69
70
71
72
73
74
75
"""
    metabolites(model::StandardModel)

Return a vector of metabolite id strings contained in `model`.
The order of metabolite strings returned here matches the order used to construct
the stoichiometric matrix.
"""
76
metabolites(model::StandardModel)::StringVecType = collect(keys(model.metabolites))
77
78
79
80
81
82

"""
n_metabolites(model::StandardModel)

Return the number of metabolites in `model`.
"""
83
n_metabolites(model::StandardModel)::Int = length(model.metabolites)
84

85
86
87
88
89
"""
    genes(model::StandardModel)

Return a vector of gene id strings in `model`.
"""
90
genes(model::StandardModel)::StringVecType = collect(keys(model.genes))
91
92
93
94
95
96

"""
    n_genes(model::StandardModel)

Return the number of genes in `model`.
"""
St. Elmo's avatar
St. Elmo committed
97
98
n_genes(model::StandardModel)::Int = length(model.genes)

99
100
101
102
103
"""
    stoichiometry(model::StandardModel)

Return the stoichiometric matrix associated with `model` in sparse format.
"""
St. Elmo's avatar
St. Elmo committed
104
105
function stoichiometry(model::StandardModel)::SparseMat
    S = SparseArrays.spzeros(length(model.metabolites), length(model.reactions))
St. Elmo's avatar
St. Elmo committed
106
    met_ids = metabolites(model) # vector of metabolite ids
St. Elmo's avatar
St. Elmo committed
107
108
    rxn_ids = reactions(model)
    for (i, rxn_id) in enumerate(rxn_ids) # column, in order
St. Elmo's avatar
St. Elmo committed
109
        for (met_id, coeff) in model.reactions[rxn_id].metabolites
110
            j = findfirst(==(met_id), met_ids) # row
St. Elmo's avatar
fix    
St. Elmo committed
111
112
113
114
115
116
117
118
            if isnothing(j)
                throw(
                    DomainError(
                        met_id,
                        "Metabolite not found in model but occurs in stoichiometry of $(rxn_id). Perhaps it was deleted?",
                    ),
                )
            end
119
120
121
122
123
124
            S[j, i] = coeff
        end
    end
    return S
end

125
126
127
128
129
"""
    lower_bounds(model::StandardModel)

Return the lower bounds for all reactions in `model` in sparse format.
"""
St. Elmo's avatar
St. Elmo committed
130
function lower_bounds(model::StandardModel)::SparseVec
131
    sparse([model.reactions[rxn].lb for rxn in reactions(model)])
St. Elmo's avatar
St. Elmo committed
132
133
end

134
135
136
137
"""
    upper_bounds(model::StandardModel)

Return the upper bounds for all reactions in `model` in sparse format.
St. Elmo's avatar
St. Elmo committed
138
Order matches that of the reaction ids returned in `reactions()`.
139
"""
St. Elmo's avatar
St. Elmo committed
140
function upper_bounds(model::StandardModel)::SparseVec
141
    sparse([model.reactions[rxn].ub for rxn in reactions(model)])
St. Elmo's avatar
St. Elmo committed
142
143
end

144
145
146
147
"""
    bounds(model::StandardModel)

Return the lower and upper bounds, respectively, for reactions in `model`.
St. Elmo's avatar
St. Elmo committed
148
Order matches that of the reaction ids returned in `reactions()`.
149
"""
cylon-x's avatar
cylon-x committed
150
function bounds(model::StandardModel)::Tuple{SparseVec,SparseVec}
St. Elmo's avatar
St. Elmo committed
151
152
    ubs = upper_bounds(model)
    lbs = lower_bounds(model)
153
154
155
    return lbs, ubs
end

156
157
158
159
160
161
"""
    balance(model::StandardModel)

Return the balance of the linear problem, i.e. b in Sv = 0 where S is the stoichiometric matrix
and v is the flux vector.
"""
162
balance(model::StandardModel)::SparseVec = spzeros(length(model.metabolites))
163

164
165
166
167
168
"""
    objective(model::StandardModel)

Return sparse objective vector for `model`.
"""
St. Elmo's avatar
St. Elmo committed
169
function objective(model::StandardModel)::SparseVec
170
171
172
    return sparse([
        model.reactions[rid].objective_coefficient for rid in keys(model.reactions)
    ])
173
end
St. Elmo's avatar
St. Elmo committed
174

175
"""
176
    reaction_gene_association(model::StandardModel, id::String)
177

178
179
Return the gene reaction rule in string format for reaction with `id` in `model`.
Return `nothing` if not available.
180
"""
St. Elmo's avatar
St. Elmo committed
181
function reaction_gene_association(model::StandardModel, id::String)::Maybe{GeneAssociation}
182
    _maybemap(identity, model.reactions[id].grr)
St. Elmo's avatar
St. Elmo committed
183
184
end

185
"""
186
    metabolite_formula(model::StandardModel, id::String)
187
188

Return the formula of reaction `id` in `model`.
189
Return `nothing` if not present.
190
"""
191
function metabolite_formula(model::StandardModel, id::String)::Maybe{MetaboliteFormula}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
192
    _maybemap(_parse_formula, model.metabolites[id].formula)
St. Elmo's avatar
St. Elmo committed
193
194
end

195
"""
196
    metabolite_charge(model::StandardModel, id::String)
197
198

Return the charge associated with metabolite `id` in `model`.
199
Return nothing if not present.
200
"""
201
function metabolite_charge(model::StandardModel, id::String)::Maybe{Int}
202
    model.metabolites[id].charge
St. Elmo's avatar
St. Elmo committed
203
204
end

205
"""
206
    metabolite_compartment(model::StandardModel, id::String)
207
208

Return compartment associated with metabolite `id` in `model`.
209
Return `nothing` if not present.
210
"""
211
function metabolite_compartment(model::StandardModel, id::String)::Maybe{String}
212
    model.metabolites[id].compartment
St. Elmo's avatar
St. Elmo committed
213
214
end

215
216
217
218
"""
    reaction_subsystem(id::String, model::StandardModel)

Return the subsystem associated with reaction `id` in `model`.
219
Return `nothing` if not present.
220
"""
221
function reaction_subsystem(model::StandardModel, id::String)::Maybe{String}
222
    model.reactions[id].subsystem
St. Elmo's avatar
St. Elmo committed
223
224
end

225
"""
226
    metabolite_notes(model::StandardModel, id::String)::Notes
227
228

Return the notes associated with metabolite `id` in `model`.
229
Return an empty Dict if not present.
230
"""
231
function metabolite_notes(model::StandardModel, id::String)::Maybe{Notes}
232
    model.metabolites[id].notes
St. Elmo's avatar
St. Elmo committed
233
234
end

235
"""
236
    metabolite_annotations(model::StandardModel, id::String)::Annotations
237
238

Return the annotation associated with metabolite `id` in `model`.
239
Return an empty Dict if not present.
240
"""
241
242
function metabolite_annotations(model::StandardModel, id::String)::Maybe{Annotations}
    model.metabolites[id].annotations
St. Elmo's avatar
St. Elmo committed
243
244
end

245
"""
246
    gene_notes(model::StandardModel, id::String)::Notes
247
248

Return the notes associated with gene `id` in `model`.
249
Return an empty Dict if not present.
250
"""
251
function gene_notes(model::StandardModel, id::String)::Maybe{Notes}
St. Elmo's avatar
St. Elmo committed
252
    model.genes[id].notes
St. Elmo's avatar
St. Elmo committed
253
254
end

255
"""
256
    gene_annotations(model::StandardModel, id::String)::Annotations
257
258

Return the annotation associated with gene `id` in `model`.
259
Return an empty Dict if not present.
260
"""
261
262
function gene_annotations(model::StandardModel, id::String)::Maybe{Annotations}
    model.genes[id].annotations
St. Elmo's avatar
St. Elmo committed
263
264
end

265
"""
266
    reaction_notes(model::StandardModel, id::String)::Notes
267
268

Return the notes associated with reaction `id` in `model`.
269
Return an empty Dict if not present.
270
"""
271
function reaction_notes(model::StandardModel, id::String)::Maybe{Notes}
St. Elmo's avatar
St. Elmo committed
272
    model.reactions[id].notes
St. Elmo's avatar
St. Elmo committed
273
274
end

275
"""
276
    reaction_annotations(model::StandardModel, id::String)::Annotations
277
278

Return the annotation associated with reaction `id` in `model`.
279
Return an empty Dict if not present.
280
"""
281
282
function reaction_annotations(model::StandardModel, id::String)::Maybe{Annotations}
    model.reactions[id].annotations
St. Elmo's avatar
St. Elmo committed
283
end
St. Elmo's avatar
St. Elmo committed
284

285
286
287
288
289
290
291
"""
Base.convert(::Type{StandardModel}, model::MetabolicModel)

Convert any `MetabolicModel` into a `StandardModel`.
Note, some data loss may occur since only the generic interface is used during
the conversion process.
"""
St. Elmo's avatar
St. Elmo committed
292
function Base.convert(::Type{StandardModel}, model::MetabolicModel)
293
294
295
296
    if typeof(model) == StandardModel
        return model
    end

St. Elmo's avatar
St. Elmo committed
297
298
299
300
301
302
303
304
305
306
    timer = TimerOutput()

    @timeit timer "instantiate dicts" begin
        id = "" # TODO: add accessor to get model ID
        modelreactions = OrderedDict{String,Reaction}()
        modelmetabolites = OrderedDict{String,Metabolite}()
        modelgenes = OrderedDict{String,Gene}()
    end

    @timeit timer "access model details" begin
St. Elmo's avatar
St. Elmo committed
307
        gids = genes(model)
St. Elmo's avatar
St. Elmo committed
308
309
310
311
312
313
        metids = metabolites(model)
        rxnids = reactions(model)
        S = stoichiometry(model)
        lbs, ubs = bounds(model)
        ocs = objective(model)
    end
cylon-x's avatar
cylon-x committed
314

St. Elmo's avatar
St. Elmo committed
315
    @timeit timer "assign genes" begin
cylon-x's avatar
cylon-x committed
316
        # gtask = Base.Threads.@spawn begin
St. Elmo's avatar
St. Elmo committed
317
        for gid in gids
St. Elmo's avatar
St. Elmo committed
318
319
            notes = @timeit timer "gene_notes" gene_notes(model, gid)
            annotations = @timeit timer "gene_annotations" gene_annotations(model, gid)
cylon-x's avatar
cylon-x committed
320
            modelgenes[gid] = Gene(gid; notes = notes, annotations = annotations) # TODO: add name accessor
St. Elmo's avatar
St. Elmo committed
321
        end
St. Elmo's avatar
St. Elmo committed
322
    end
St. Elmo's avatar
St. Elmo committed
323

St. Elmo's avatar
St. Elmo committed
324
    @timeit timer "assign metabolites" begin
cylon-x's avatar
cylon-x committed
325
        # mtask = Base.Threads.@spawn begin
St. Elmo's avatar
St. Elmo committed
326
        for mid in metids
St. Elmo's avatar
St. Elmo committed
327
            charge = @timeit timer "charge" metabolite_charge(model, mid)
cylon-x's avatar
cylon-x committed
328
329
330
331
            formula = @timeit timer "formula" _maybemap(
                _unparse_formula,
                metabolite_formula(model, mid),
            )
St. Elmo's avatar
St. Elmo committed
332
333
334
            compartment = @timeit timer "compartment" metabolite_compartment(model, mid)
            notes = @timeit timer "notes" metabolite_notes(model, mid)
            annotations = @timeit timer "annotations" metabolite_annotations(model, mid)
St. Elmo's avatar
St. Elmo committed
335
336
            modelmetabolites[mid] = Metabolite(
                mid;
St. Elmo's avatar
St. Elmo committed
337
338
339
340
341
                charge = charge,
                formula = formula,
                compartment = compartment,
                notes = notes,
                annotations = annotations,
St. Elmo's avatar
St. Elmo committed
342
            )
St. Elmo's avatar
St. Elmo committed
343
        end
St. Elmo's avatar
St. Elmo committed
344
    end
St. Elmo's avatar
St. Elmo committed
345

St. Elmo's avatar
St. Elmo committed
346
    @timeit timer "assign reactions" begin
cylon-x's avatar
cylon-x committed
347
        # rtask = Base.Threads.@spawn begin
St. Elmo's avatar
St. Elmo committed
348

St. Elmo's avatar
St. Elmo committed
349
        for (i, rid) in enumerate(rxnids)
St. Elmo's avatar
St. Elmo committed
350
351
352
353
354
355

            grr = @timeit timer "grr" reaction_gene_association(model, rid)
            notes = @timeit timer "notes" reaction_notes(model, rid)
            annotations = @timeit timer "annotations" reaction_annotations(model, rid)
            subsys = @timeit timer "subsystem" reaction_subsystem(model, rid)
            @timeit timer "get reaction equation" begin
cylon-x's avatar
cylon-x committed
356
357
358
359
                rmets = Dict{String,Float64}()
                for (j, stoich) in zip(findnz(S[:, i])...)
                    rmets[metids[j]] = stoich
                end
St. Elmo's avatar
St. Elmo committed
360
            end
St. Elmo's avatar
St. Elmo committed
361
362
363
364
365
            modelreactions[rid] = Reaction(
                rid;
                metabolites = rmets,
                lb = lbs[i],
                ub = ubs[i],
St. Elmo's avatar
St. Elmo committed
366
                grr = grr,
St. Elmo's avatar
St. Elmo committed
367
                objective_coefficient = ocs[i],
St. Elmo's avatar
St. Elmo committed
368
369
370
                notes = notes,
                annotations = annotations,
                subsystem = subsys,
St. Elmo's avatar
St. Elmo committed
371
            ) # TODO: add name accessor
St. Elmo's avatar
St. Elmo committed
372
373
        end
    end
cylon-x's avatar
cylon-x committed
374

St. Elmo's avatar
St. Elmo committed
375
376
377
    # wait(rtask)
    # wait(mtask)
    # wait(gtask)
cylon-x's avatar
cylon-x committed
378

St. Elmo's avatar
St. Elmo committed
379
380
    show(timer)
    println()
cylon-x's avatar
cylon-x committed
381
382
    return StandardModel(
        id;
383
384
385
        reactions = modelreactions,
        metabolites = modelmetabolites,
        genes = modelgenes,
cylon-x's avatar
cylon-x committed
386
    )
St. Elmo's avatar
St. Elmo committed
387
end