MATModel.jl 7.17 KB
Newer Older
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
1
2
3
4
5
6
"""
    struct MATModel

Wrapper around the models loaded in dictionaries from the MATLAB representation.
"""
struct MATModel <: MetabolicModel
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
7
    mat::Dict{String,Any}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
8
9
end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
10
11
n_metabolites(m::MATModel)::Int = size(m.mat["S"], 1)
n_reactions(m::MATModel)::Int = size(m.mat["S"], 2)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
12

13
14
15
16
17
"""
    reactions(m::MATModel)::Vector{String}

Extracts reaction names from `rxns` key in the MAT file.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
18
function reactions(m::MATModel)::Vector{String}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
19
20
    if haskey(m.mat, "rxns")
        reshape(m.mat["rxns"], n_reactions(m))
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
21
22
23
24
25
    else
        "rxn" .* string.(1:n_reactions(m))
    end
end

26
27
28
29
30
31
32
33
34
"""
    _mat_has_squashed_coupling(mat)

Guesses whether C in the MAT file is stored in A=[S;C].
"""
_mat_has_squashed_coupling(mat) =
    haskey(mat, "A") && haskey(mat, "b") && length(mat["b"]) == size(mat["A"], 1)


35
36
37
38
39
"""
    metabolites(m::MATModel)::Vector{String}

Extracts metabolite names from `mets` key in the MAT file.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
40
function metabolites(m::MATModel)::Vector{String}
41
    nm = n_metabolites(m)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
42
    if haskey(m.mat, "mets")
43
        reshape(m.mat["mets"], length(m.mat["mets"]))[begin:nm]
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
44
45
46
47
48
    else
        "met" .* string.(1:n_metabolites(m))
    end
end

49
50
51
52
53
"""
    stoichiometry(m::MATModel)

Extract the stoichiometry matrix, stored under key `S`.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
54
stoichiometry(m::MATModel) = sparse(m.mat["S"])
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
55

56
57
58
59
60
"""
    bounds(m::MATModel)

Extracts bounds from the MAT file, saved under `lb` and `ub`.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
61
bounds(m::MATModel) = (
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
62
63
    sparse(reshape(get(m.mat, "lb", fill(-Inf, n_reactions(m), 1)), n_reactions(m))),
    sparse(reshape(get(m.mat, "ub", fill(Inf, n_reactions(m), 1)), n_reactions(m))),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
64
65
)

66
67
68
69
70
"""
    balance(m::MATModel)

Extracts balance from the MAT model, defaulting to zeroes if not present.
"""
71
72
73
74
75
76
77
function balance(m::MATModel)
    b = get(m.mat, "b", spzeros(n_metabolites(m), 1))
    if _mat_has_squashed_coupling(m.mat)
        b = b[1:n_metabolites(m), :]
    end
    sparse(reshape(b, n_metabolites(m)))
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
78

79
80
81
82
83
"""
    objective(m::MATModel)

Extracts the objective from the MAT model (defaults to zeroes).
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
84
85
objective(m::MATModel) =
    sparse(reshape(get(m.mat, "c", zeros(n_reactions(m), 1)), n_reactions(m)))
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
86

87
88
89
90
91
"""
    coupling(m::MATModel)

Extract coupling matrix stored, in `C` key.
"""
92
93
94
coupling(m::MATModel) =
    _mat_has_squashed_coupling(m.mat) ? sparse(m.mat["A"][n_reactions(m)+1:end, :]) :
    sparse(get(m.mat, "C", zeros(0, n_reactions(m))))
95
96
97
98
99
100
101
102

"""
    coupling_bounds(m::MATModel)

Extracts the coupling constraints. Currently, there are several accepted ways to store these in MATLAB models; this takes the constraints from vectors `cl` and `cu`.
"""
function coupling_bounds(m::MATModel)
    nc = n_coupling_constraints(m)
103
104
105
106
107
108
109
110
111
112
113
    if _mat_has_squashed_coupling(m.mat)
        (
            sparse(fill(-Inf, nc)),
            sparse(reshape(m.mat["b"], length(m.mat["b"]))[n_reactions(m)+1:end]),
        )
    else
        (
            sparse(reshape(get(m.mat, "cl", fill(-Inf, nc, 1)), nc)),
            sparse(reshape(get(m.mat, "cu", fill(Inf, nc, 1)), nc)),
        )
    end
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
end

"""
    genes(m::MATModel)

Extracts the possible gene list from `genes` key.
"""
function genes(m::MATModel)
    x = get(m.mat, "genes", [])
    reshape(x, length(x))
end

"""
    reaction_gene_association(m::MATModel, rid::String)

Extracts the associations from `grRules` key, if present.
"""
function reaction_gene_association(m::MATModel, rid::String)
    if haskey(m.mat, "grRules")
133
134
        grr = m.mat["grRules"][findfirst(==(rid), reactions(m))]
        typeof(grr) == String ? _parse_grr(grr) : nothing
135
136
137
138
139
140
141
142
143
144
    else
        nothing
    end
end

"""
    metabolite_formula(m::MATModel, mid::String)

Extract metabolite formula from key `metFormula` or `metFormulas`.
"""
145
metabolite_formula(m::MATModel, mid::String) = _maybemap(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
146
    x -> _parse_formula(x[findfirst(==(mid), metabolites(m))]),
147
148
149
150
151
152
153
154
    get(m.mat, "metFormula", get(m.mat, "metFormulas", nothing)),
)

"""
    metabolite_charge(m::MATModel, mid::String)

Extract metabolite charge from `metCharge` or `metCharges`.
"""
155
metabolite_charge(m::MATModel, mid::String) = _maybemap(
156
157
158
159
160
161
162
163
164
    x -> x[findfirst(==(mid), metabolites(m))],
    get(m.mat, "metCharge", get(m.mat, "metCharges", nothing)),
)

"""
    metabolite_compartment(m::MATModel, mid::String)

Extract metabolite compartment from `metCompartment` or `metCompartments`.
"""
165
metabolite_compartment(m::MATModel, mid::String) = _maybemap(
166
167
168
169
    x -> x[findfirst(==(mid), metabolites(m))],
    get(m.mat, "metCompartment", get(m.mat, "metCompartments", nothing)),
)

St. Elmo's avatar
St. Elmo committed
170
171

"""
St. Elmo's avatar
St. Elmo committed
172
    reaction_stoichiometry(model::MATModel, rxn_id::String)::Dict{String, Float64}
St. Elmo's avatar
St. Elmo committed
173
174
175
176

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
177
function reaction_stoichiometry(m::MATModel, rxn_id::String)::Dict{String,Float64}
St. Elmo's avatar
St. Elmo committed
178
    rxn_ind = first(indexin([rxn_id], m.mat["rxns"]))
St. Elmo's avatar
St. Elmo committed
179
    reaction_stoichiometry(m, rxn_ind)
St. Elmo's avatar
St. Elmo committed
180
181
182
end

"""
St. Elmo's avatar
St. Elmo committed
183
    reaction_stoichiometry(model::MATModel, rxn_index)::Dict{String, Float64}
St. Elmo's avatar
St. Elmo committed
184

St. Elmo's avatar
St. Elmo committed
185
186
187
Return the reaction equation of reaction with index `rxn_index` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients. Note, `rxn_index` can 
be any suitable type that can index into an array.
St. Elmo's avatar
St. Elmo committed
188
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
189
function reaction_stoichiometry(m::MATModel, rxn_ind)::Dict{String,Float64}
St. Elmo's avatar
St. Elmo committed
190
191
192
193
    met_inds = findall(m.mat["S"][:, rxn_ind] .!= 0.0)
    Dict(m.mat["mets"][met_ind] => m.mat["S"][met_ind, rxn_ind] for met_ind in met_inds)
end

194
195
196
197
198
199
200
# NOTE: There's no useful standard on how and where to store notes and
# annotations in MATLAB models. We therefore leave it very open for the users,
# who can easily support any annotation scheme using a custom wrapper.
# Even the (simple) assumptions about grRules, formulas and charges that we use
# here are very likely completely incompatible with >50% of the MATLAB models
# out there.

201
202
"""
    Base.convert(::Type{MATModel}, m::MetabolicModel)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
203

204
205
Convert any metabolic model to `MATModel`.
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
206
function Base.convert(::Type{MATModel}, m::MetabolicModel)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
207
208
209
210
    if typeof(m) == MATModel
        return m
    end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
211
    lb, ub = bounds(m)
212
    cl, cu = coupling_bounds(m)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
213
214
215
216
217
    nr = n_reactions(m)
    nm = n_metabolites(m)
    return MATModel(
        Dict(
            "S" => stoichiometry(m),
218
219
220
221
222
223
            "rxns" => reactions(m),
            "mets" => metabolites(m),
            "lb" => Vector(lb),
            "ub" => Vector(ub),
            "b" => Vector(balance(m)),
            "c" => Vector(objective(m)),
224
225
226
227
228
            "C" => coupling(m),
            "cl" => Vector(cl),
            "cu" => Vector(cu),
            "genes" => genes(m),
            "grRules" =>
229
                _default.(
230
                    "",
231
                    _maybemap.(
cylon-x's avatar
cylon-x committed
232
                        x -> _unparse_grr(String, x),
233
234
235
236
                        reaction_gene_association.(Ref(m), reactions(m)),
                    ),
                ),
            "metFormulas" =>
237
                _default.(
238
                    "",
239
                    _maybemap.(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
240
                        _unparse_formula,
241
242
243
                        metabolite_formula.(Ref(m), metabolites(m)),
                    ),
                ),
244
            "metCharges" => _default.(0, metabolite_charge.(Ref(m), metabolites(m))),
245
            "metCompartments" =>
246
                _default.("", metabolite_compartment.(Ref(m), metabolites(m))),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
247
248
249
        ),
    )
end