Unverified Commit 5c6c4f92 authored by St. Elmo's avatar St. Elmo
Browse files

fixed gecko

parent c240484a
Pipeline #55274 failed with stages
in 9 minutes and 49 seconds
......@@ -139,11 +139,11 @@ function make_gecko_model(
mg_gid_lookup[mg] = [gid]
end
end
coupling_row_mass_group = Vector{Tuple{Vector{Int}, Vector{Float64}, Float64}}()
coupling_row_mass_group = Vector{Tuple{String, Vector{Int}, Vector{Float64}, Float64}}()
for (grp, gs) in mg_gid_lookup
idxs = Int.(indexin(gs, gids))
idxs = [gene_row_lookup[x] for x in Int.(indexin(gs, gids))]
mms = gpmm_.(gs)
push!(coupling_row_mass_group, (idxs, mms, gmgb_(grp)))
push!(coupling_row_mass_group, (grp, idxs, mms, gmgb_(grp)))
end
gm = GeckoModel(
......
......@@ -55,7 +55,7 @@ struct GeckoModel <: ModelWrapper
columns::Vector{_gecko_column}
coupling_row_reaction::Vector{Int}
coupling_row_gene_product::Vector{Tuple{Int,Tuple{Float64,Float64}}}
coupling_row_mass_group::Vector{Tuple{Vector{Int}, Vector{Float64}, Float64}}
coupling_row_mass_group::Vector{Tuple{String, Vector{Int}, Vector{Float64}, Float64}}
inner::MetabolicModel
end
......@@ -91,11 +91,10 @@ objective(model::GeckoModel) = model.objective
Returns the internal reactions in a [`GeckoModel`](@ref) (these may be split
to forward- and reverse-only parts with different isozyme indexes; reactions
IDs are mangled accordingly with suffixes), as well as the genes associated
with enzymatic reactions. In the context of a GeckoModel, this is better described
as the variables in the problem.
with enzymatic reactions.
"""
function reactions(model::GeckoModel)
rxnnames = let inner_reactions = reactions(model.inner)
reactions(model::GeckoModel) =
let inner_reactions = reactions(model.inner)
[
_gecko_reaction_name(
inner_reactions[col.reaction_idx],
......@@ -104,18 +103,14 @@ function reactions(model::GeckoModel)
) for col in model.columns
]
end
[rxnnames; genes(model)]
end
"""
n_reactions(model::GeckoModel)
Returns the number of all irreversible reactions in `model` as well as the number of gene products
that take part in enzymatic reactions. In the context of a GeckoModel, this is better described
as the number of variables in the problem.
that take part in enzymatic reactions.
"""
n_reactions(model::GeckoModel) =
length(model.columns) + length(model.coupling_row_gene_product)
n_reactions(model::GeckoModel) = length(reactions(model))
"""
bounds(model::GeckoModel)
......@@ -140,14 +135,7 @@ end
Get the mapping of the reaction rates in [`GeckoModel`](@ref) to the original
fluxes in the wrapped model.
"""
function reaction_flux(model::GeckoModel)
gecko_mat = _gecko_column_reactions(model)'
inner_mat = reaction_flux(model.inner)
[
gecko_mat*inner_mat spzeros(size(gecko_mat, 1), n_genes(model))
spzeros(n_genes(model), size(inner_mat, 2)) I(n_genes(model))
]
end
reaction_flux(model::GeckoModel) = _gecko_column_reactions(model)' * reaction_flux(model.inner)
"""
coupling(model::GeckoModel)
......@@ -196,7 +184,7 @@ function coupling_bounds(model::GeckoModel)
vcat(
icub,
iub[model.coupling_row_reaction],
[ub for (_, _, ub) in model.coupling_row_mass_group],
[ub for (_, _, _, ub) in model.coupling_row_mass_group],
),
)
end
......@@ -214,7 +202,7 @@ balance(model::GeckoModel) = [balance(model.inner); spzeros(length(model.couplin
Return the number of genes that have enzymatic constraints associated with them.
"""
n_genes(model::GeckoModel) = length(model.coupling_row_gene_product)
n_genes(model::GeckoModel) = length(genes(model))
"""
genes(model::GeckoModel)
......@@ -224,9 +212,15 @@ Return the gene ids of genes that have enzymatic constraints associated with the
genes(model::GeckoModel) = genes(model.inner)[[idx for (idx, _) in model.coupling_row_gene_product]]
"""
fluxes(model::GeckoModel)
metabolites(model::GeckoModel)
Return the ids of all metabolites, both real and pseudo, for a [`GeckoModel`](@ref).
"""
metabolites(model::GeckoModel) = [metabolites(model.inner); genes(model).*"#supply"]
"""
n_metabolites(model::GeckoModel)
Return the original reaction ids and the gene ids that were variables
in a [`GeckoModel`](@ref).
Return the number of metabolites, both real and pseudo, for a [`GeckoModel`](@ref).
"""
fluxes(model::GeckoModel) = [reactions(model.inner); genes(model)]
\ No newline at end of file
n_metabolites(model::GeckoModel) = length(metabolites(model))
\ No newline at end of file
......@@ -6,12 +6,8 @@ argument `opt_model` is a solved optimization problem, typically returned by
[`flux_balance_analysis`](@ref).
"""
protein_dict(model::GeckoModel, opt_model) =
let gids = genes(model)
is_solved(opt_model) ?
Dict(
[gids[gidx] for (gidx, _) in model.coupling_row_gene_product] .=> _gecko_gene_product_coupling(model) * value.(opt_model[:x]),
) : nothing
end
is_solved(opt_model) ?
Dict(genes(model) .=> value.(opt_model[:x])[(n_reactions(model)+1):end]) : nothing
"""
protein_dict(model::GeckoModel)
......@@ -28,8 +24,8 @@ Extract the mass utilization in mass groups from a solved [`GeckoModel`](@ref).
protein_mass_group_dict(model::GeckoModel, opt_model) =
is_solved(opt_model) ?
Dict(
(group for (group, _) in model.coupling_row_mass_group) .=>
_gecko_mass_group_coupling(model) * value.(opt_model[:x]),
grp[1] => dot(value.(opt_model[:x])[n_reactions(model) .+ grp[2]], grp[3]) for
grp in model.coupling_row_mass_group
) : nothing
"""
......@@ -39,7 +35,6 @@ A pipe-able variant of [`mass_group_dict`](@ref).
"""
protein_mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
"""
protein_mass(model::SMomentModel)
......
......@@ -73,7 +73,7 @@ mass of each group of gene products.
function _gecko_mass_group_coupling(model::GeckoModel)
tmp = [
(i, j, mm) for (i, mg) in enumerate(model.coupling_row_mass_group) for
(j, mm) in zip(mg[1], mg[2])
(j, mm) in zip(mg[2], mg[3])
]
sparse(
[i for (i, _, _) in tmp],
......@@ -82,4 +82,14 @@ function _gecko_mass_group_coupling(model::GeckoModel)
length(model.coupling_row_mass_group),
n_genes(model),
)
end
\ No newline at end of file
end
"""
flux_dict(model::GeckoModel, opt_model)
Returns the fluxes (not gene product concentrations) of the model as a
reaction-keyed dictionary, if solved.
"""
flux_dict(model::GeckoModel, opt_model) =
is_solved(opt_model) ?
Dict(fluxes(model) .=> reaction_flux(model)' * value.(opt_model[:x])[1:n_reactions(model)] ) : nothing
\ No newline at end of file
......@@ -49,3 +49,107 @@
@test isapprox(prot_mass, total_protein_mass, atol = TEST_TOLERANCE)
end
@testset "GECKO small model" begin
#=
Implement the small model found in the supplment of the
original GECKO paper. This model is nice to troubleshoot with,
because the stoich matrix is small.
=#
m = StandardModel("gecko")
m1 = Metabolite("m1")
m2 = Metabolite("m2")
m3 = Metabolite("m3")
m4 = Metabolite("m4")
@add_reactions! m begin
"r1", nothing m1, -100, 100
"r2", nothing m2, -100, 100
"r3", m1 + m2 m3, 0, 100
"r4", m3 m4, -100, 100 # make reversible instead
"r5", m2 m4, -100, 100
"r6", nothing m4, 0, 100
end
gs = [Gene("g$i") for i in 1:4]
m.reactions["r3"].grr = [["g1"]]
m.reactions["r4"].grr = [["g1"], ["g2"]]
m.reactions["r5"].grr = [["g3", "g4"]]
m.reactions["r4"].objective_coefficient = 1.0
add_genes!(m, gs)
add_metabolites!(m, [m1, m2, m3, m4])
reaction_isozymes = Dict(
"r3" => [
Isozyme(
Dict("g1" => 1),
1.0,
1.0,
),
],
"r4" => [
Isozyme(
Dict("g1" => 1),
2.0,
2.0,
),
Isozyme(
Dict("g2" => 1),
3.0,
3.0,
),
],
"r5" => [
Isozyme(
Dict("g3" => 1, "g4" => 2),
5.0,
5.0,
),
],
)
gene_product_bounds = Dict(
"g1" => (0.0, 10.0),
"g2" => (0.0, 10.0),
"g3" => (0.0, 10.0),
"g4" => (0.0, 10.0),
)
gene_product_molar_mass = Dict(
"g1" => 1.0,
"g2" => 2.0,
"g3" => 3.0,
"g4" => 4.0,
)
gene_mass_group_bound = Dict("uncategorized" => 0.5)
gm = make_gecko_model(
m;
reaction_isozymes,
gene_product_bounds,
gene_product_molar_mass,
gene_mass_group_bound,
)
S = stoichiometry(gm)
l, u = bounds(gm)
coupling(gm)
cl, cu = coupling_bounds(gm)
opt_model = flux_balance_analysis(
gm,
Tulip.Optimizer;
modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
rxn_fluxes = flux_dict(gm, opt_model)
gene_products = protein_dict(gm, opt_model)
mass_groups = protein_mass_group_dict(gm, opt_model)
@test isapprox(rxn_fluxes["r4"], 0.142857, atol = TEST_TOLERANCE)
@test isapprox(gene_products["g3"], 0.0285714, atol = TEST_TOLERANCE)
@test isapprox(mass_groups["uncategorized"], 0.5, atol = TEST_TOLERANCE)
end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment