gapfill_minimum_reactions.jl 6.52 KB
Newer Older
St. Elmo's avatar
St. Elmo committed
1
"""
2
    function gapfill_minimum_reactions(
St. Elmo's avatar
St. Elmo committed
3
4
5
        model::MetabolicModel,
        universal_reactions::Vector{Reaction},
        optimizer;
6
7
        objective_bounds = (_constants.tolerance, _constants.default_reaction_bound),
        maximum_new_reactions = 5,
St. Elmo's avatar
St. Elmo committed
8
        weights = fill(1.0, length(universal_reactions)),
9
        modifications = [],
St. Elmo's avatar
St. Elmo committed
10
    )
11
12
13
14
15
16
17
18
19
20
21
22
23

Find a minimal set of reactions from `universal_reactions` that should be added
to `model` so that the model has a feasible solution with bounds on its
objective function given in `objective_bounds`.  Weights of the added reactions
may be specified in `weights` to prefer adding reactions with lower weights.

Internally, this builds and solves a mixed integer program, following
the method of Reed et al. (Reed, Jennifer L., et al. "Systems approach to
refining genome annotation." *Proceedings of the National Academy of Sciences*
(2006)).

The function returns a solved JuMP optimization model, with the boolean
reaction inclusion indicators in variable vector `y`. Use
24
[`gapfilled_mask`](@ref) or [`gapfilled_rids`](@ref) to collect the reaction
25
26
27
28
29
information in Julia datatypes.

To reduce the uncertainty in the MILP solver (and likely reduce the
complexity), you may put a limit on the size of the added reaction set in
`maximum_new_reactions`.
St. Elmo's avatar
St. Elmo committed
30
31
32
33
34
"""
function gapfill_minimum_reactions(
    model::MetabolicModel,
    universal_reactions::Vector{Reaction},
    optimizer;
35
36
    objective_bounds = (_constants.tolerance, _constants.default_reaction_bound),
    maximum_new_reactions = length(universal_reactions),
St. Elmo's avatar
St. Elmo committed
37
    weights = fill(1.0, length(universal_reactions)),
38
    modifications = [],
St. Elmo's avatar
St. Elmo committed
39
)
40
    precache!(model)
St. Elmo's avatar
St. Elmo committed
41
42

    # constraints from universal reactions that can fill gaps
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
    univs = _universal_stoichiometry(universal_reactions, metabolites(model))

    # add space for additional metabolites and glue with the universal reaction
    # stoichiometry
    extended_stoichiometry = [[
        stoichiometry(model)
        spzeros(length(univs.new_mids), n_reactions(model))
    ] univs.stoichiometry]

    # make the model anew (we can't really use make_optimization_model because
    # we need the balances and several other things completely removed. Could
    # be solved either by parametrizing make_optimization_model or by making a
    # tiny temporary wrapper for this.
    # keep this in sync with src/base/solver.jl, except for adding balances.
    opt_model = Model(optimizer)
    @variable(opt_model, x[1:n_reactions(model)])
    xl, xu = bounds(model)
    @constraint(opt_model, lbs, xl .<= x)
    @constraint(opt_model, ubs, x .<= xu)

    C = coupling(model)
    isempty(C) || begin
        cl, cu = coupling_bounds(model)
        @constraint(opt_model, c_lbs, cl .<= C * x)
        @constraint(opt_model, c_ubs, C * x .<= cu)
    end
St. Elmo's avatar
St. Elmo committed
69

70
71
72
    # add the variables for new stuff
    @variable(opt_model, ux[1:length(universal_reactions)]) # fluxes from universal reactions
    @variable(opt_model, y[1:length(universal_reactions)], Bin) # indicators
St. Elmo's avatar
St. Elmo committed
73

74
    # combined metabolite balances
St. Elmo's avatar
St. Elmo committed
75
76
    @constraint(
        opt_model,
77
78
        extended_stoichiometry * [x; ux] .==
        [balance(model); zeros(length(univs.new_mids))]
St. Elmo's avatar
St. Elmo committed
79
80
    )

81
82
83
84
    # objective bounds
    @constraint(opt_model, objective_bounds[1] <= objective(model)' * x)
    @constraint(opt_model, objective_bounds[2] >= objective(model)' * x)

St. Elmo's avatar
St. Elmo committed
85
    # flux bounds of universal reactions with indicators
86
87
    @constraint(opt_model, ulb, univs.lbs .* y .<= ux)
    @constraint(opt_model, uub, univs.ubs .* y .>= ux)
St. Elmo's avatar
St. Elmo committed
88

89
90
    # minimize the total number of indicated reactions
    @objective(opt_model, Min, weights' * y)
St. Elmo's avatar
St. Elmo committed
91

92
93
94
    # limit the number of indicated reactions
    # (prevents the solver from exploring too far)
    @constraint(opt_model, sum(y) <= maximum_new_reactions)
St. Elmo's avatar
St. Elmo committed
95

96
    # apply all modifications
St. Elmo's avatar
St. Elmo committed
97
98
99
100
101
102
    for mod in modifications
        mod(opt_model, model)
    end

    optimize!(opt_model)

103
    return opt_model
St. Elmo's avatar
St. Elmo committed
104
105
end

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
"""
    gapfilled_mask(opt_model::BitVector)

Get a `BitVector` of added reactions from the model solved by
[`gapfill_minimum_reactions`](@ref). The bit indexes correspond to the indexes
of `universal_reactions` given to the gapfilling function. In case the model is
not solved, this returns `nothing`.

# Example

    gapfill_minimum_reactions(myModel, myReactions, Tulip.Optimizer) |> gapfilled_mask
"""
gapfilled_mask(opt_model)::BitVector =
    is_solved(opt_model) ? value.(opt_model[:y]) .> 0 : nothing

"""
    gapfilled_rids(opt_model, universal_reactions::Vector{Reaction})::Vector{String}

Utility to extract a short vector of IDs of the reactions added by the
gapfilling algorithm. Use with `opt_model` returned from
[`gapfill_minimum_reactions`](@ref).
"""
gapfilled_rids(opt_model, universal_reactions::Vector{Reaction}) =
    let v = gapfilled_mask(opt_model)
        isnothing(v) ? nothing : [rxn.id for rxn in universal_reactions[v]]
    end

"""
    gapfilled_rids(universal_reactions::Vector{Reaction})

Overload of [`gapfilled_rids`](@ref) that can be piped easily.

# Example

    gapfill_minimum_reactions(myModel, myReactions, Tulip.Optimizer) |> gapfilled_rids(myReactions)
"""
gapfilled_rids(universal_reactions::Vector{Reaction}) =
    opt_model -> gapfilled_rids(opt_model, universal_reactions)

St. Elmo's avatar
St. Elmo committed
145
146
147
"""
    _universal_stoichiometry(
        universal_reactions::Vector{Reaction},
148
        mids,
St. Elmo's avatar
St. Elmo committed
149
150
151
152
    )

A helper function that constructs the stoichiometric matrix of a set of
`universal_reactions`. The order of the metabolites is determined with
153
`mids`, so that this stoichiometric matrix can be combined with
St. Elmo's avatar
St. Elmo committed
154
155
another one.
"""
156
157
158
159
160
161
162
function _universal_stoichiometry(urxns::Vector{Reaction}, mids::Vector{String})

    # traversal over all elements in stoichiometry of universal_reactions
    stoiMap(f) = [
        f(ridx, mid, stoi) for (ridx, rxn) in enumerate(urxns) for
        (mid, stoi) in rxn.metabolites
    ]
St. Elmo's avatar
St. Elmo committed
163

164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
    # make an index and find new metabolites
    met_id_lookup = Dict(mids .=> eachindex(mids))

    new_mids =
        collect(Set(filter(x -> !haskey(met_id_lookup, x), stoiMap((_, mid, _) -> mid))))
    all_mids = vcat(mids, new_mids)

    # remake the index with all metabolites
    met_id_lookup = Dict(all_mids .=> eachindex(all_mids))

    # build the result
    return (
        stoichiometry = float.(
            sparse(
                stoiMap((_, mid, _) -> met_id_lookup[mid]),
                stoiMap((ridx, _, _) -> ridx),
                stoiMap((_, _, stoi) -> stoi),
                length(all_mids),
                length(urxns),
            ),
        ),
        lbs = [rxn.lb for rxn in urxns],
        ubs = [rxn.ub for rxn in urxns],
        new_mids = new_mids,
    )
St. Elmo's avatar
St. Elmo committed
189
end