flux_variability_analysis.jl 5.18 KB
Newer Older
Sylvain Arreckx's avatar
Sylvain Arreckx committed
1
"""
2
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
        model::MetabolicModel,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
4
        reactions::Vector{Int},
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
5
6
7
8
9
10
11
12
13
        optimizer;
        modifications = [],
        workers = [myid()],
        bounds = z -> (z,z),
        ret = objective_value,
    )::Matrix{Float64}

Flux variability analysis solves a pair of optimization problems in `model` for
each flux listed in `reactions`:
Sylvain Arreckx's avatar
Sylvain Arreckx committed
14
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
15
 min,max xᵢ
Sylvain Arreckx's avatar
Sylvain Arreckx committed
16
s.t. S x = b
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
17
18
19
    xₗ ≤ x ≤ xᵤ
     cᵀx ≥ bounds(Z₀)[1]
     cᵀx ≤ bounds(Z₀)[2]
Sylvain Arreckx's avatar
Sylvain Arreckx committed
20
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
21
where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated
St. Elmo's avatar
St. Elmo committed
22
23
24
FBA problem (see [`flux_balance_analysis`](@ref)). See "Gudmundsson, S., Thiele,
I. Computationally efficient flux variability analysis. BMC Bioinformatics 11,
489 (2010). https://doi.org/10.1186/1471-2105-11-489" for more information.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
25
26
27
28

The `bounds` is a user-supplied function that specifies the objective bounds
for the variability optimizations, by default it restricts the flux objective
value to the precise optimum reached in FBA. It can return `-Inf` and `Inf` in
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
29
30
first and second pair to remove the limit. Use [`gamma_bounds`](@ref) and
[`objective_bounds`](@ref) for simple bounds.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31
32
33
34
35
36

`optimizer` must be set to a `JuMP`-compatible optimizer. The computation of
the individual optimization problems is transparently distributed to `workers`
(see `Distributed.workers()`).

`ret` is a function used to extract results from optimized JuMP models of the
37
38
39
40
individual reactions. By default, it calls and returns the value of
`JuMP.objective_value`. More information can be extracted e.g. by setting it to
a function that returns a more elaborate data structure; such as `m ->
(JuMP.objective_value(m), JuMP.value.(m[:x]))`.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
41
42

Returns a matrix of extracted `ret` values for minima and maxima, of total size
43
44
45
(`length(reactions)`,2). The optimizer result status is checked with
[`is_solved`](@ref); `nothing` is returned if the optimization failed for any
reason.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
46
"""
47
function flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
48
    model::MetabolicModel,
49
    reactions::Vector{Int},
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
50
51
52
53
54
55
    optimizer;
    modifications = [],
    workers = [myid()],
    bounds = z -> (z, z),
    ret = objective_value,
)
56
    if any(reactions .< 1) || any(reactions .> n_reactions(model))
57
        throw(DomainError(reactions, "Index exceeds number of reactions."))
58
59
    end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
60
61
62
63
64
    Z = bounds(
        objective_value(
            flux_balance_analysis(model, optimizer; modifications = modifications),
        ),
    )
65

66
67
68
    # store a JuMP optimization model at all workers
    save_model = :(
        begin
69
70
71
72
73
            model = $model
            opt_model = $COBREXA.make_optimization_model(model, $optimizer)
            for mod in $modifications
                mod(model, opt_model)
            end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
74
75
            $COBREXA._FVA_add_constraint(opt_model, $(objective(model)), opt_model[:x], $Z)
            opt_model
76
        end
77
    )
78
79
    map(fetch, save_at.(workers, :cobrexa_parfva_model, Ref(save_model)))
    save_model = nothing # this has some volume, free it again
80

81
    # schedule FVA parts parallely using pmap
82
    fluxes = dpmap(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
83
        rid -> :($COBREXA._FVA_optimize_reaction(cobrexa_parfva_model, $rid, $ret)),
84
        CachingPool(workers),
85
86
        [-reactions reactions],
    )
87

88
    # free the data on workers
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
89
    map(fetch, remove_from.(workers, :cobrexa_parfva_model))
90

91
92
    return fluxes
end
93

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
94
"""
95
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
96
        model::MetabolicModel,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
97
        optimizer;
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
98
99
        kwargs...
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
100

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
101
A simpler version of [`flux_variability_analysis`](@ref) that maximizes and minimizes all reactions in the model. Arguments are forwarded.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
102
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
103
function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
104
    n = n_reactions(model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
105
    return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
106
107
end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
108
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
109
110
111
112
113
    flux_variability_analysis_dict(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
114

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
115
116
117
A variant of [`flux_variability_analysis`](@ref) that returns the individual
maximized and minimized fluxes of all reactions as two dictionaries (of
dictionaries). All keyword arguments except `ret` are passed through.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
118
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
119
120
121
122
123
function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
    vs = flux_variability_analysis(
        model,
        optimizer;
        kwargs...,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
124
        ret = m -> JuMP.value.(m[:x]),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
125
126
    )
    rxns = reactions(model)
127

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
128
    return (
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
129
130
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 1]])),
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 2]])),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
131
    )
132
end
133

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
134
135
"""
    _FVA_add_constraint(model, c, x, Z)
St. Elmo's avatar
St. Elmo committed
136

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
137
138
139
140
141
142
143
Internal helper function for adding constraints to a model. Exists mainly
because for avoiding namespace problems on remote workers.
"""
function _FVA_add_constraint(model, c, x, Z)
    Z[1] > -Inf && @constraint(model, c' * x >= Z[1])
    Z[2] < Inf && @constraint(model, c' * x <= Z[2])
end
St. Elmo's avatar
St. Elmo committed
144

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
145
146
"""
    _FVA_get_opt(model, rid)
147

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
148
149
150
151
Internal helper for creating the optimized model on a remote worker, for
avoiding namespace problems.
"""
function _FVA_optimize_reaction(model, rid, ret)
St. Elmo's avatar
St. Elmo committed
152
    sense = rid > 0 ? MAX_SENSE : MIN_SENSE
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
153
    var = all_variables(model)[abs(rid)]
154

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
155
156
    @objective(model, sense, var)
    optimize!(model)
157

St. Elmo's avatar
St. Elmo committed
158
    if is_solved(model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
159
160
161
162
        return ret(model)
    else
        return nothing
    end
163
end