parsimonious_flux_balance_analysis.jl 4.71 KB
Newer Older
St. Elmo's avatar
St. Elmo committed
1
"""
St. Elmo's avatar
St. Elmo committed
2
    pfba(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
St. Elmo's avatar
St. Elmo committed
3
4
5
6
7

Run parsimonious flux balance analysis (pFBA) on the `model`.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
Optionally, specify problem `modifications` as in [`flux_balance_analysis`](@ref).
Also, `qp_solver` can be set to be different from `optimizer`, where the latter is then the LP optimizer only.
St. Elmo's avatar
St. Elmo committed
8
Also note that `qp_solver_attributes` is meant to set the attributes for the `qp_solver`.
St. Elmo's avatar
St. Elmo committed
9
10
11
12
13
14
15
16
17
Any solver attributes changed in `modifications` will only affect he LP solver.
This function automatically relaxes the constraint that the FBA solution objective matches the pFBA solution. 
This is iteratively relaxed like 1.0, 0.999999, 0.9999, etc. of the bound until 0.99 when the function fails and returns nothing.
Return a solved pFBA model.

# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
18
model = load_model(StandardModel, "iJO1366.json")
St. Elmo's avatar
St. Elmo committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = pfba(model, biomass, optimizer; solver_attributes=atts)
```
"""
function parsimonious_flux_balance_analysis(
    model::MetabolicModel,
    optimizer;
    modifications = [(model, opt_model) -> nothing],
    qp_solver = (model, opt_model) -> nothing,
    qp_solver_attributes = [(model, opt_model) -> nothing],
)
    # Run FBA
    opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
    COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
        return nothing # FBA failed

    # Run pFBA
    λ = objective_value(opt_model)
    old_objective = COBREXA.JuMP.objective_function(opt_model)

    qp_solver(model, opt_model) # change the solver if specified, otherwise default argument does nothing
    if typeof(qp_solver_attributes) <: AbstractVector # many modifications
        for mod in qp_solver_attributes
            mod(model, opt_model)
        end
    else # single modification
        qp_solver_attributes(model, opt_model)
    end

    v = opt_model[:x] # fluxes
    @objective(opt_model, Min, sum(dot(v, v)))

    for lbconval in [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99] # sequentially relax bound for stability
St. Elmo's avatar
St. Elmo committed
52
53
        λmin = min(lbconval * λ, λ * 1.0 / lbconval)
        λmax = max(lbconval * λ, λ * 1.0 / lbconval)
St. Elmo's avatar
St. Elmo committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
        @constraint(
            opt_model,
            pfbacon,
            λmin <= old_objective <= λmax # in case of negative constraints
        )
        optimize!(opt_model)
        COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] &&
            break
        COBREXA.JuMP.delete(opt_model, pfbacon)
    end

    COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
        return nothing # pFBA failed

    return opt_model
end

"""
St. Elmo's avatar
St. Elmo committed
72
    parsimonious_flux_balance_analysis_dict(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
St. Elmo's avatar
St. Elmo committed
73
74
75

Perform parsimonious flux balance analysis on `model` using `optimizer`. 
Returns a vector of fluxes in the same order as the reactions in `model`. 
St. Elmo's avatar
St. Elmo committed
76
Calls [`parsimonious_flux_balance_analysis`](@ref) internally.
St. Elmo's avatar
St. Elmo committed
77
78
"""
function parsimonious_flux_balance_analysis_vec(
St. Elmo's avatar
St. Elmo committed
79
    model::MetabolicModel,
St. Elmo's avatar
St. Elmo committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    optimizer;
    modifications = [(model, opt_model) -> nothing],
    qp_solver = (model, opt_model) -> nothing,
    qp_solver_attributes = [(model, opt_model) -> nothing],
)
    opt_model = parsimonious_flux_balance_analysis(
        model,
        optimizer;
        modifications = modifications,
        qp_solver = qp_solver,
        qp_solver_attributes = qp_solver_attributes,
    )
    COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
        return nothing

    return value.(opt_model[:x])
end

"""
St. Elmo's avatar
St. Elmo committed
99
    parsimonious_flux_balance_analysis_dict(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
St. Elmo's avatar
St. Elmo committed
100
101
102

Perform parsimonious flux balance analysis on `model` using `optimizer`. 
Returns a dictionary mapping reaction `id`s to fluxes. 
St. Elmo's avatar
St. Elmo committed
103
Calls [`parsimonious_flux_balance_analysis`](@ref) internally.
St. Elmo's avatar
St. Elmo committed
104
105
"""
function parsimonious_flux_balance_analysis_dict(
St. Elmo's avatar
St. Elmo committed
106
    model::MetabolicModel,
St. Elmo's avatar
St. Elmo committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    optimizer;
    modifications = [(model, opt_model) -> nothing],
    qp_solver = (model, opt_model) -> nothing,
    qp_solver_attributes = [(model, opt_model) -> nothing],
)
    opt_model = parsimonious_flux_balance_analysis(
        model,
        optimizer;
        modifications = modifications,
        qp_solver = qp_solver,
        qp_solver_attributes = qp_solver_attributes,
    )
    COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
        return nothing

    return Dict(zip(reactions(model), value.(opt_model[:x])))
end