Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Gap filling produces nonempty result when filling a model that is already feasible #78

Closed
oxinabox opened this issue Nov 14, 2024 · 4 comments · Fixed by #80
Closed

Gap filling produces nonempty result when filling a model that is already feasible #78

oxinabox opened this issue Nov 14, 2024 · 4 comments · Fixed by #80

Comments

@oxinabox
Copy link

I was showing someone that if you give gapfilling a model that is feasible already,
even if you change the objective to be something that can be improved by knocking in some genes, it still won't say to insert them.
To do this i found 10 genes i could KO that would decrease CO2 production,
and then set objective to maximise CO2 production (which is kinda nonsense to do but anyway)
and ran gap-filling
to my surprise it actually decided to fill in something.
But it wasn't even something I had knocked out (though it was something that had fixed bounds, so knocking in an extra copy accelerates it)

In contrast, when I run this exact same investigation in Cobrapy (I can provide code for that if you want) this does not occur and as expected it return an empty set.

There may be some interaction with having basically identical model and universal model,
as I did not see this when i used metanetx 's largest universal model (which probably has this reaction).

Minimal code example to reproduce the problem

using my fav yeast-9:
https://raw.githubusercontent.com/SysBioChalmers/yeast-GEM/refs/heads/main/model/yeast-GEM.xml

using COBREXA
using SBMLFBCModels
using HiGHS
using AbstractFBCModels:AbstractFBCModels, reaction_gene_products_available, reactions


function ko_model!(model, knockout_genes)
    for reaction_id in reactions(model)
        can_run = reaction_gene_products_available(model, reaction_id, !in(knockout_genes))
        #@show can_run
        if !something(can_run, true)
            reaction = model.reactions[reaction_id]
            @info "KOing" reaction.name
            reaction.lower_bound = reaction.upper_bound = 0.0
        end
    end
    return model
end


function set_product_objective!(model)
    for reaction in values(model.reactions)
        uris = get(reaction.annotations, "RESOURCE_URI", String[]) 
        reaction.objective_coefficient = "https://identifiers.org/bigg.reaction/EX_co2_e" ∈ uris
    end
    return model
end

data_dir = joinpath(dirname(@__DIR__), "data")
yeast_file = joinpath(data_dir, "yeast-GEM.xml")
const original_model = load_model(yeast_file)
const model = convert(AbstractFBCModels.CanonicalModel.Model, original_model)  # so we can mutate it directly

ko_genes = ["YLR028C", "YBR111C", "YOR125C", "YNL038W", "YJL134W", "YGR202C", "YGR055W", "YER152C", "YJR040W", "YGR243W"]
ko_model!(model, ko_genes)
set_product_objective!(model)

x = gap_filling_analysis(
    model,
    original_model,
    0.05,
    optimizer = HiGHS.Optimizer,
)

Expected result

It should not try and fill anything, because the model is already feasible

ConstraintTrees.Tree{Float64} with 6 elements:
  :fill_flags            => ConstraintTrees.Tree{Float64}(#= 4130 elements =#)
  :n_filled              => 0.0
  :....

Actual behavior

It tries to fill 1 reaction

ConstraintTrees.Tree{Float64} with 6 elements:
  :fill_flags            => ConstraintTrees.Tree{Float64}(#= 4130 elements =#)
  :n_filled              => 1.0
  :....

you can see which with:

@show any(!=(0), x.fill_flags)  # this should be false
cuni = convert(AbstractFBCModels.CanonicalModel.Model, original_model)
r = only([cuni.reactions[string(k)] for (k,v) in x.fill_flags if v!=0])

What was the unexpected thing that happened instead?

The extra surprising thing is it didn't even put in one i deleted.
It put in one that was already there, and which doesn't have a gene associated

AbstractFBCModels.CanonicalModel.Reaction(
  name = "non-growth associated maintenance reaction",
  lower_bound = 0.7,
  upper_bound = 0.7,
  stoichiometry = Dict("M_s_0803" => -1.0, "M_s_0434" => -1.0, "M_s_0394" => 1.0, "M_s_0794" => 1.0, "M_s_1322" => 1.0),
  objective_coefficient = 0.0,
  gene_association_dnf = nothing,
  annotations = Dict("sbo" => ["SBO:0000630"], "RESOURCE_URI" => ["https://identifiers.org/bigg.reaction/ATPM"]),
  notes = Dict("" => ["<notes>\n  <body xmlns=\"http://www.w3.org/1999/xhtml\">\n    <p>Confidence Level: 0</p>\n  </body>\n</notes>"]),
)
@exaexa
Copy link
Member

exaexa commented Nov 14, 2024

Hi! Thanks for the detailed reproducer, this indeed looks weird.

If I get it right, the reaction that gets gapfilled is the ATP maintenance or something similar? If that is so, it gets gapfilled literally because it is constrainted to strictly nonzero flux, and if you wouldn't include it in the fill, it would have a zero value in the model, making the model infeasible. This is the outcome of how the problem is encoded into a MILP, and of the fact that original bounds in the universal reactions are retained.

If the above is the case: maintenance pseudoreactions should not be included in universal reactions (they are more like extra constraints encoded as reactions), and removing them will solve the issue.

I recall we have this removal of ATPM in at least one of the examples. If it is not in the gapfilling one, we should certainly add it there. Maybe together with a mild warning about this possible issue.

If the reason is different then please do elaborate, because that would be very weird.

@exaexa
Copy link
Member

exaexa commented Nov 14, 2024

PS re why cobrapy works: they might either auto-ignore these reactions, or use a different encoding into MILP that also selectively ignores the constraints, or fix/remove the bounds. Will have a look.

@oxinabox
Copy link
Author

This is indeed ATPM.
It is mentioned removing it in example but no details as to why.
Maybe we need to add those details to the docs or have a function or a fag to remove this kind automatically?

This was referenced Nov 22, 2024
Merged
Open
@exaexa
Copy link
Member

exaexa commented Nov 22, 2024

automated cleaning is kinda desirable so I'll continue with this within a greater scope of #81

@exaexa exaexa closed this as completed in a4dd3b3 Nov 27, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

add an explicit warning about what happens if people gapfill pseudoreactions COBREXA/COBREXA.jl
2 participants
@exaexa @oxinabox