Gapfilling not giving any solutions

[HettieC]

A model made by removing reactions from iML1515 and then gapfilled from iML1515 to do the reaction

Dict("atp_c"=>-1,
"adp_c"=>1,
"pi_c"=>1,
"h2o_c"=>-1,
"h_c"=>1)

doesn't give a solution. test.txt

Here is the gapfilling code:

using COBREXA, Tulip, GLPK, JuMP

ecoli_rxns = collect(values(ecoli.reactions))
for r in ecoli_rxns
    r.lb = -1000
    r.ub = 1000
end

# first see if model can produce atp
mets = Dict("atp_c"=>-1,
"adp_c"=>1,
"pi_c"=>1,
"h2o_c"=>-1,
"h_c"=>1)

biomass_atp = Reaction("biomass_atp";
metabolites=mets,
objective_coefficient=1.0)

add_reaction!(model,biomass_atp)

a = flux_balance_analysis_dict(model,Tulip.Optimizer)
### this gives zero

#### gapfill for atp production with all iML1515 reactions

m = gapfill_atp = gapfill_minimum_reactions(
    model,
    ecoli_rxns,
    GLPK.Optimizer,
    maximum_new_reactions=1000,
)

gapfilled_rids(m,ecoli_rxns) # this is nothing

for r in values(ecoli.reactions)
    if r.objective_coefficient != 0.0
        r.objective_coefficient = 0
    end
end

add_reaction!(ecoli,biomass_atp)
b = flux_balance_analysis_dict(ecoli,Tulip.Optimizer)
## ecoli model can definitely do this reaction as flux_summary(b) is non-zero.

[exaexa]

Hi! just to be sure, what's ecoli and what is model ?

I assume

model = ecoli = load_model(StandardModel, "iML1515.json")

and the iML1515.json is from BIGG?

[exaexa]

OK, as an observation the a actually is feasible in my case, but with really weird biomass production (something like 999 or so). Please can you send a standalone script to avoid uncertainty in what model and ecoli is?

With the model and ecoli loaded as above, I get:

julia> a
Dict{String, Float64} with 2713 entries:
  "PACOAT"        => 18.2651
  "Zn2tex"        => 0.142092
  "GUI1"          => -25.4034
  "DXYLK"         => 0.0
  "FE3DCITtonex"  => 9.96505
  "METSOXR1"      => -32.3627
  "FACOAL180t2pp" => -50.5778
  "CBL1tonex"     => 10.7797
  "LIPOtex"       => -0.0159199
  "NTD11"         => -159.944
  "GLUNpp"        => 134.854
  "ORNDC"         => -55.5699
  "UAGCVT"        => -107.138
  "ALAGLUE"       => -6.79532
  "I2FE2ST"       => 10.8841
  "BMOCOS"        => 3.69732e-16
  "6D6SPA"        => -46.1326
  "GLCURt2rpp"    => 26.0364
  "EX_g3pg_e"     => -34.1899
  "PGSA141"       => -31.9436
  "GMPS2"         => 35.7147
  "RNDR2b"        => -7.69883
  "NMNPtpp"       => 42.8984
  "EX_gua_e"      => -10.9186
  "ACGK"          => -59.7721
  ⋮               => ⋮

julia> a["BIOMASS_Ec_iML1515_core_75p37M"]
999.9999979361254

julia> gapfilled_rids(m,ecoli_rxns)
String[]

[HettieC]

Hey! ecoli = load_model(StandardModel, "iML1515.json") but model = load_model(StandardModel, "test.json") .

I couldn't upload a .json file to the github issue so it is saved as a .txt, hopefully you can still load this into cobrexa... if not i can send the file directly to you as .json.

So model is what I made by blasting and then taking reactions from ecoli that were hits in the bidirectional blast.

[exaexa]

Ah I see, good, thanks, I'll try asap.

[exaexa]

OK so after some verification in slack it turns out the model is already feasible, so the gapfilling doesn't need to add anything.

Solution: use objective_bounds= parameter for the gapfilling to force it to actually add something.

I'll mark this bug as a documentation bug because it's not at all obvious from the current doc.

Maybe we might print a small explanation in gapfilled_rids or so, that tells the user that the model might have been feasible already.

[exaexa]

closed by #662