Unverified Commit 3b355fb1 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil Committed by GitHub
Browse files

Merge pull request #617 from LCSB-BioCore/mo-gapfil

Add a gap filling algorithm
parents f6f086c2 fa0dc910
Pipeline #56586 passed with stages
in 21 minutes and 5 seconds
......@@ -19,16 +19,18 @@ function make_optimization_model(model::MetabolicModel, optimizer; sense = MAX_S
xl, xu = bounds(model)
optimization_model = Model(optimizer)
@variable(optimization_model, x[i = 1:n])
@variable(optimization_model, x[1:n])
@objective(optimization_model, sense, objective(model)' * x)
@constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
@constraint(optimization_model, lbs, xl .<= x) # lower bounds
@constraint(optimization_model, ubs, x .<= xu) # upper bounds
C = coupling(model) # empty if no coupling
cl, cu = coupling_bounds(model)
isempty(C) || @constraint(optimization_model, c_lbs, cl .<= C * x) # coupling lower bounds
isempty(C) || @constraint(optimization_model, c_ubs, C * x .<= cu) # coupling upper bounds
isempty(C) || begin
cl, cu = coupling_bounds(model)
@constraint(optimization_model, c_lbs, cl .<= C * x) # coupling lower bounds
@constraint(optimization_model, c_ubs, C * x .<= cu) # coupling upper bounds
end
return optimization_model
end
......
"""
function gapfill_minimum_reactions(
model::MetabolicModel,
universal_reactions::Vector{Reaction},
optimizer;
objective_bounds = (_constants.tolerance, _constants.default_reaction_bound),
maximum_new_reactions = 5,
weights = fill(1.0, length(universal_reactions)),
modifications = [],
)
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
[`mask`](@ref) or [`gapfilled_rids`](@ref) to collect the reaction
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`.
"""
function gapfill_minimum_reactions(
model::MetabolicModel,
universal_reactions::Vector{Reaction},
optimizer;
objective_bounds = (_constants.tolerance, _constants.default_reaction_bound),
maximum_new_reactions = length(universal_reactions),
weights = fill(1.0, length(universal_reactions)),
modifications = [],
)
precache!(model)
# constraints from universal reactions that can fill gaps
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
# 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
# combined metabolite balances
@constraint(
opt_model,
extended_stoichiometry * [x; ux] .==
[balance(model); zeros(length(univs.new_mids))]
)
# objective bounds
@constraint(opt_model, objective_bounds[1] <= objective(model)' * x)
@constraint(opt_model, objective_bounds[2] >= objective(model)' * x)
# flux bounds of universal reactions with indicators
@constraint(opt_model, ulb, univs.lbs .* y .<= ux)
@constraint(opt_model, uub, univs.ubs .* y .>= ux)
# minimize the total number of indicated reactions
@objective(opt_model, Min, weights' * y)
# limit the number of indicated reactions
# (prevents the solver from exploring too far)
@constraint(opt_model, sum(y) <= maximum_new_reactions)
# apply all modifications
for mod in modifications
mod(opt_model, model)
end
optimize!(opt_model)
return opt_model
end
"""
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)
"""
_universal_stoichiometry(
universal_reactions::Vector{Reaction},
mids,
)
A helper function that constructs the stoichiometric matrix of a set of
`universal_reactions`. The order of the metabolites is determined with
`mids`, so that this stoichiometric matrix can be combined with
another one.
"""
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
]
# 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,
)
end
@testset "Gap fill with minimum reactions" begin
#=
Implement the small model that should be gapfilled.
=#
model = StandardModel("partial model")
(m1, m2, m3, m4, m5, m6, m7, m8) = Metabolite.("m$i" for i = 1:8)
@add_reactions! model begin
"r1", nothing m1, 0, 1
"r2", m1 m2, -10, 100
"r3", m1 m3, 0, 100
"r4", m2 m4, 0, 100
# "r5", m3 → m4, 0, 100
"r6", m4 nothing, 0, 100
# "r7", m2 → m7 + m6, 0, 100
"r8", m7 m8, 0, 100
"r9", m8 nothing, 0, 100
# "r10", m6 → nothing, 0, 100
"r11", m2 + m3 + m7 nothing, 0, 100
"r12", m3 m5, -10, 10
end
model.reactions["r11"].objective_coefficient = 1.0
add_metabolites!(model, [m1, m2, m3, m4, m5, m7, m8])
r5 = Reaction("r5", Dict("m3" => -1, "m4" => 1), :forward)
r7 = Reaction("r7", Dict("m2" => -1, "m7" => 1, "m6" => 1), :forward)
r10 = Reaction("r10", Dict("m6" => -1), :forward)
rA = Reaction("rA", Dict("m1" => -1, "m2" => 1, "m3" => 1), :forward)
rB = Reaction("rB", Dict("m2" => -1, "m9" => 1), :forward)
rC = Reaction("rC", Dict("m9" => -1, "m10" => 1), :bidirectional)
rD = Reaction("rC", Dict("m10" => -1), :reverse)
universal_reactions = [r5, r7, r10, rA, rB, rC, rD]
rxns =
gapfill_minimum_reactions(
model,
universal_reactions,
GLPK.Optimizer;
objective_bounds = (0.1, 1000.0),
) |> gapfilled_rids(universal_reactions)
@test issetequal(["r7", "r10"], rxns)
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment