fba.jl 5.6 KB
Newer Older
1
"""
2
    flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4

Flux balance analysis solves the following problem for the input `model`:
Sylvain Arreckx's avatar
Sylvain Arreckx committed
5
6
7
8
9
```
max cᵀx
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
10

11
Returns a solved model from [`optimize_model`](@ref).
Sylvain Arreckx's avatar
Sylvain Arreckx committed
12
"""
13
14
flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel} =
    optimize_model(model, optimizer; sense = MOI.MAX_SENSE)
15

16
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
17
    flux_balance_analysis_vec(args...)::Union{Vector{Float64},Nothing}
18
19
20
21
22

A variant of FBA that returns a vector of fluxes in the same order as reactions
of the model, if the solution is found.
Arguments are passed to [`flux_balance_analysis`](@ref).
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
23
function flux_balance_analysis_vec(args...)::Union{Vector{Float64},Nothing}
St. Elmo's avatar
St. Elmo committed
24
    optmodel = flux_balance_analysis(args...)
25
26
    vars = optmodel[:x]

27
    termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || return nothing
28
    value.(vars)
29
30
31
end

"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
32
    flux_balance_analysis_dict(model::M, args...)::Union{Dict{String, Float64},Nothing} where {M <: MetabolicModel}
33
34
35
36
37
38
39

A variant of FBA that returns a dictionary assigning fluxes to reactions, if
the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref).
"""
function flux_balance_analysis_dict(
    model::M,
    args...,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
40
)::Union{Dict{String,Float64},Nothing} where {M<:MetabolicModel}
41
42
43
44
45
    v = flux_balance_analysis_vec(model, args...)
    isnothing(v) && return nothing
    Dict(zip(reactions(model), v))
end

46
"""
47
    fba(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
48

St. Elmo's avatar
St. Elmo committed
49
Run flux balance analysis (FBA) on the `model` optionally specifying `objective_rxn(s)` and their `weights` (empty `weights` mean equal weighting per reaction).
50
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.
51
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
52
53
54
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.
55

56
57
58
59
60
61
62
63
# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = fba(model, biomass, optimizer; solver_attributes=atts)
```
64
"""
65
function fba(
66
67
    model::CobraModel,
    optimizer;
68
    modifications = [(model, opt_model)->nothing]
69
)
70
71
72
73
74
75

objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[]
weights = Float64[]
solver_attributes = Dict{Any,Any}()
sense = MOI.MAX_SENSE

St. Elmo's avatar
St. Elmo committed
76
    # get core optimization problem
77
78
79
80
81
82
    cbm = make_optimization_model(model, optimizer, sense = sense)
    v = cbm[:x] # fluxes

    # apply callbacks
    for mod in modifications
        mod(model, cbm)
St. Elmo's avatar
St. Elmo committed
83
    end
St. Elmo's avatar
St. Elmo committed
84

85
86
87
88
89
90
    # modify core optimization problem according to user specifications
    if !isempty(solver_attributes) # set other attributes
        for (k, val) in solver_attributes
            set_optimizer_attribute(cbm, k, val)
        end
    end
91

St. Elmo's avatar
St. Elmo committed
92
<<<<<<< Updated upstream
93
94
95
96
97
98
99
100
    # if an objective function is supplied, modify the default objective
    if typeof(objective_func) == Reaction || !isempty(objective_func)
        # ensure that an array of objective indices are fed in
        if typeof(objective_func) == Reaction
            objective_indices = [model[objective_func]]
        else
            objective_indices = [model[rxn] for rxn in objective_func]
        end
101

102
103
104
105
        if isempty(weights)
            weights = ones(length(objective_indices))
        end
        opt_weights = zeros(length(model.reactions))
St. Elmo's avatar
St. Elmo committed
106

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
        # update the objective function tracker
        # don't update model objective function - silly thing to do
        wcounter = 1
        for i in eachindex(model.reactions)
            if i in objective_indices
                # model.reactions[i].objective_coefficient = weights[wcounter]
                opt_weights[i] = weights[wcounter]
                wcounter += 1
                # else
                # model.reactions[i].objective_coefficient = 0.0
            end
        end

        @objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
    else # use default objective
        # automatically assigned by make_optimization_model
    end

    optimize!(cbm)

St. Elmo's avatar
St. Elmo committed
127
128
129
130
=======
function flux_balance_analysis_vec(args)
    cbm = flux_balance_analysis(args...)
    
131
132
133
134
    status = (
        termination_status(cbm) == MOI.OPTIMAL ||
        termination_status(cbm) == MOI.LOCALLY_SOLVED
    )
St. Elmo's avatar
St. Elmo committed
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    
    if status
        return value.(cbm[:x])
    else
        @warn "Optimization issues occurred."
        return Vector{Float64}[]
    end    
end

function flux_balance_analysis_dict(model, args)
    cbm = flux_balance_analysis(model, args...)
    
>>>>>>> Stashed changes
    status = (
        termination_status(cbm) == MOI.OPTIMAL ||
        termination_status(cbm) == MOI.LOCALLY_SOLVED
    )
<<<<<<< Updated upstream
St. Elmo's avatar
St. Elmo committed
153

154
155
156
157
158
159
    if status
        return map_fluxes(v, model)
    else
        @warn "Optimization issues occurred."
        return Dict{String,Float64}()
    end
St. Elmo's avatar
St. Elmo committed
160
161
162
163
164
165
166
167
168
=======
    
    if status
        return Dict(zip(reactions(model), value.(cbm[:x])))
    else
        @warn "Optimization issues occurred."
        return Dict{String, Float64}()
    end   
>>>>>>> Stashed changes
169
end