-
Notifications
You must be signed in to change notification settings - Fork 8
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Added loopless fba #60
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nice, and thanks for adding in tests!
2a98b5a
to
a90f7cb
Compare
a90f7cb
to
1cf4854
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See comments :)
|
||
# Is there a way to get the objective indices if we set the | ||
# objective via `change_objective`? then we could recycle the flux_balance_analysis | ||
# method in a better way |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, JumP.objective_function
, look at line 37 of parimonious_flux_balance_analysis.jl
for a use case.
objective_func::Union{Reaction,Array{Reaction,1}}, | ||
weights=[], |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can be wrapped into the modifications? Ideally all the analysis functions would have the same arguments where possible
objective_indices = _get_reaction_indices(model, objective_func) | ||
opt_weights = _get_objective_weights(model, objective_indices, weights) | ||
change_objective(objective_func, weights=weights)(model, opt_model) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not sure this is necessary if you have the callbacks.
function _set_solver_attibutes!(opt_model, solver_attributes) | ||
if !isempty(solver_attributes) # set other attributes | ||
for (k, val) in solver_attributes | ||
set_optimizer_attribute(opt_model, k, val) | ||
end | ||
end | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Already have similar functionality in src/base/utilities.jl
as change_solver_attributes
.
function _set_additional_constraints!(model, opt_model, constraints) | ||
for (rxnid, con) in constraints | ||
ind = model.reactions[findfirst(model.reactions, rxnid)] | ||
set_bound(ind, opt_model; lb=con[1], ub=con[2]) | ||
end | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See previous comment and change_constraint
function _get_reaction_indices(model, objective_func) | ||
if typeof(objective_func) == Reaction | ||
objective_indices = [model[objective_func]] | ||
else | ||
objective_indices = [model[rxn] for rxn in objective_func] | ||
end | ||
return objective_indices | ||
end | ||
|
||
function _get_objective_weights(model, objective_indices, weights) | ||
if isempty(weights) | ||
weights = ones(length(objective_indices)) | ||
end | ||
opt_weights = zeros(length(model.reactions)) | ||
|
||
# update the objective function tracker | ||
wcounter = 1 | ||
for i in eachindex(model.reactions) | ||
if i in objective_indices | ||
opt_weights[i] = weights[wcounter] | ||
wcounter += 1 | ||
end | ||
end | ||
return opt_weights | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could probably read this out using JuMP, see comment below.
function _constrain_objective_value!(opt_model, opt_weights, objective_indices) | ||
λ = objective_value(opt_model) | ||
v = opt_model[:x] | ||
@constraint( | ||
opt_model, | ||
λ <= sum(opt_weights[i] * v[i] for i in objective_indices) <= λ | ||
) | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not sure if necessary to make this a function. Also look at parsimonious_flux_balance_analysis.jl
to see how to make this more robust, specifically lines 51 - 58.
opt_weights, | ||
objective_indices, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
hopefully not necessary after comments
return loopless_solution( | ||
model, | ||
opt_model, | ||
opt_weights, | ||
objective_indices, | ||
fluxes, | ||
) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nice!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Very nice! Just a couple of tiny changes and then it'll be ready :)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] || | ||
return nothing | ||
return Dict(zip(reactions(model), value.(opt_model[:x]))) | ||
function add_loopless(model::M, opt_model, fluxes::Dict{String,Float64})::JuMP.Model where {M <: MetabolicModel} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we need COBREXA.JuMP.Model
here
""" | ||
loopless_solution(args...)::Union{Dict{String,Float64},Nothing} | ||
|
||
Convert an existing FBA solution to a loopless one. Removes as many loops as possible usingthe method from [1] | ||
|
||
References: | ||
[1] CycleFreeFlux: efficient removal of thermodynamically infeasible | ||
loops from flux distributions. Desouki AA, Jarre F, Gelius-Dietrich | ||
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi: | ||
10.1093/bioinformatics/btv096. | ||
""" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Don't get rid of this!
function constrain_objective_value() | ||
return (model, opt_model) -> begin | ||
λ = objective_value(opt_model) | ||
old_objective = objective_function(opt_model) | ||
@constraint( | ||
opt_model, | ||
λ <= sum(old_objective) <= λ | ||
) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe move this to src/utilities.jl
since it is now a callback, it could maybe be used in other functions. But give it an optimum_bound
argument so that I can directly use it in src/analysis/flux_variability_analysis.jl
on line 159 - 166 :)
Ey, can we separate out the "fixes" of CobraModel to StandardModel in In particular, I'm doing that in parallel now, so I was kinda expected to see the conflict here :] Otherwise it looks cool, will give it a deeper dive later today |
@marvinvanaalst and @exaexa can we wrap this up? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
👍 from me (if the formatting doesn't trigger formatcheck failure)
EDIT: apparently it doesn't trigger any failure, feel free to merge. I'll eventually move my cleaning monsterpullrequest onto this.
@@ -97,7 +97,7 @@ function flux_balance_analysis_dict( | |||
model::M, | |||
args...; | |||
kwargs..., | |||
)::Union{Dict{String,Float64},Nothing} where {M<:MetabolicModel} | |||
)::Union{Dict{String,Float64},Nothing} where {M <: MetabolicModel} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this might cause a tiny problem later, Formatter formats the typenames without spaces around <:
for me. Does this come from the formatter?
G, Lercher MJ. Bioinformatics. 2015 Jul 1;31(13):2159-65. doi: | ||
10.1093/bioinformatics/btv096. | ||
""" | ||
function add_loopless_vec(model::M, opt_model, fluxes::Dict{String,Float64})::Union{Array{Float64,1},Nothing} where {M <: MetabolicModel} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this looks nice, I'll see if we can unify the code with the main FBA later.
top_n = 8, | ||
ignorebound = 1000.0, | ||
verbose = true, | ||
top_n=8, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is this from JuliaFormatter?
Added loopless FBA
Two interfaces:
It should be straight-forward to replace parts of the functionality, as most of it was pulled out into separate functions. This should make merging it with other developments in the FBA routine easier, since it heavily depends on it