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

implement reviews

parent 0127dfc4
Pipeline #55313 failed with stages
in 10 minutes and 10 seconds
......@@ -4,8 +4,8 @@
reaction_isozymes::Union{Function,Dict{String,Vector{Isozyme}}}
gene_product_bounds::Union{Function,Dict{String,Tuple{Float64,Float64}}},
gene_product_molar_mass::Union{Function,Dict{String,Float64}},
gene_mass_group::Union{Function,Dict{String,String}} = _ -> "uncategorized",
gene_mass_group_bound::Union{Function,Dict{String,Float64}},
gene_product_mass_group::Union{Function,Dict{String,String}} = _ -> "uncategorized",
gene_product_mass_group_bound::Union{Function,Dict{String,Float64}},
)
Wrap a model into a [`GeckoModel`](@ref), following the structure given by
......@@ -20,11 +20,11 @@ GECKO algorithm (see [`GeckoModel`](@ref) documentation for details).
`reaction_isozymes`), as `Tuple{Float64,Float64}`.
- `gene_product_molar_mass` is a function that returns a numeric molar mass of
a given gene product specified by string gene ID.
- `gene_mass_group` is a function that returns a string group identifier for a
- `gene_product_mass_group` is a function that returns a string group identifier for a
given gene product, again specified by string gene ID. By default, all gene
products belong to group `"uncategorized"` which is the behavior of original
GECKO.
- `gene_mass_group_bound` is a function that returns the maximum mass for a given
- `gene_product_mass_group_bound` is a function that returns the maximum mass for a given
mass group.
Alternatively, all function arguments may be also supplied as dictionaries that
......@@ -35,8 +35,8 @@ function make_gecko_model(
reaction_isozymes::Union{Function,Dict{String,Vector{Isozyme}}},
gene_product_bounds::Union{Function,Dict{String,Tuple{Float64,Float64}}},
gene_product_molar_mass::Union{Function,Dict{String,Float64}},
gene_mass_group::Union{Function,Dict{String,String}} = _ -> "uncategorized",
gene_mass_group_bound::Union{Function,Dict{String,Float64}},
gene_product_mass_group::Union{Function,Dict{String,String}} = _ -> "uncategorized",
gene_product_mass_group_bound::Union{Function,Dict{String,Float64}},
)
ris_ =
reaction_isozymes isa Function ? reaction_isozymes :
......@@ -47,10 +47,12 @@ function make_gecko_model(
gpmm_ =
gene_product_molar_mass isa Function ? gene_product_molar_mass :
(gid -> gene_product_molar_mass[gid])
gmg_ = gene_mass_group isa Function ? gene_mass_group : (gid -> gene_mass_group[gid])
gmg_ =
gene_product_mass_group isa Function ? gene_product_mass_group :
(gid -> gene_product_mass_group[gid])
gmgb_ =
gene_mass_group_bound isa Function ? gene_mass_group_bound :
(grp -> gene_mass_group_bound[grp])
gene_product_mass_group_bound isa Function ? gene_product_mass_group_bound :
(grp -> gene_product_mass_group_bound[grp])
# ...it would be nicer to have an overload for this, but kwargs can't be used for dispatch
columns = Vector{_gecko_column}()
......@@ -138,11 +140,11 @@ function make_gecko_model(
mg_gid_lookup[mg] = [gid]
end
end
coupling_row_mass_group = Vector{Tuple{String,Vector{Int},Vector{Float64},Float64}}()
coupling_row_mass_group = Vector{_gecko_capacity}()
for (grp, gs) in mg_gid_lookup
idxs = [gene_row_lookup[x] for x in Int.(indexin(gs, gids))]
mms = gpmm_.(gs)
push!(coupling_row_mass_group, (grp, idxs, mms, gmgb_(grp)))
push!(coupling_row_mass_group, _gecko_capacity(grp, idxs, mms, gmgb_(grp)))
end
# create model with dummy objective
......
......@@ -13,6 +13,19 @@ struct _gecko_column
gene_product_coupling::Vector{Tuple{Int,Float64}}
end
"""
struct _gecko_capacity
A helper struct that contains the gene product capacity terms organized by
the grouping type, e.g. metabolic or membrane groups etc.
"""
struct _gecko_capacity
group_id::String
gene_product_idxs::Vector{Int}
gene_product_molar_masses::Vector{Float64}
group_upper_bound::Float64
end
"""
struct GeckoModel <: ModelWrapper
......@@ -55,7 +68,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{String,Vector{Int},Vector{Float64},Float64}}
coupling_row_mass_group::Vector{_gecko_capacity}
inner::MetabolicModel
end
......@@ -113,7 +126,7 @@ 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.
"""
n_reactions(model::GeckoModel) = length(reactions(model))
n_reactions(model::GeckoModel) = length(model.columns)
"""
bounds(model::GeckoModel)
......@@ -138,8 +151,13 @@ end
Get the mapping of the reaction rates in [`GeckoModel`](@ref) to the original
fluxes in the wrapped model.
"""
reaction_flux(model::GeckoModel) =
_gecko_column_reactions(model)' * reaction_flux(model.inner)
function reaction_flux(model::GeckoModel)
rxnmat = _gecko_column_reactions(model)' * reaction_flux(model.inner)
[
rxnmat
spzeros(n_genes(model), size(rxnmat, 2))
]
end
"""
coupling(model::GeckoModel)
......@@ -188,7 +206,7 @@ function coupling_bounds(model::GeckoModel)
vcat(
icub,
iub[model.coupling_row_reaction],
[ub for (_, _, _, ub) in model.coupling_row_mass_group],
[grp.group_upper_bound for grp in model.coupling_row_mass_group],
),
)
end
......@@ -207,7 +225,7 @@ balance(model::GeckoModel) =
Return the number of genes that have enzymatic constraints associated with them.
"""
n_genes(model::GeckoModel) = length(genes(model))
n_genes(model::GeckoModel) = length(model.coupling_row_gene_product)
"""
genes(model::GeckoModel)
......@@ -229,4 +247,4 @@ metabolites(model::GeckoModel) = [metabolites(model.inner); genes(model) .* "#su
Return the number of metabolites, both real and pseudo, for a [`GeckoModel`](@ref).
"""
n_metabolites(model::GeckoModel) = length(metabolites(model))
n_metabolites(model::GeckoModel) = n_metabolites(model.inner) + n_genes(model)
"""
protein_dict(model::GeckoModel, opt_model)
gene_product_dict(model::GeckoModel, opt_model)
Return a dictionary mapping protein molar concentrations to their ids. The
argument `opt_model` is a solved optimization problem, typically returned by
[`flux_balance_analysis`](@ref). See [`flux_dict`](@ref) for the corresponding
function that returns a dictionary of solved fluxes.
"""
protein_dict(model::GeckoModel, opt_model) =
gene_product_dict(model::GeckoModel, opt_model) =
is_solved(opt_model) ?
Dict(genes(model) .=> value.(opt_model[:x])[(n_reactions(model)+1):end]) : nothing
"""
protein_dict(model::GeckoModel)
gene_product_dict(model::GeckoModel)
A pipe-able variant of [`protein_dict`](@ref).
A pipe-able variant of [`gene_product_dict`](@ref).
"""
protein_dict(model::GeckoModel) = x -> protein_dict(model, x)
gene_product_dict(model::GeckoModel) = x -> gene_product_dict(model, x)
"""
protein_mass_group_dict(model::GeckoModel, opt_model)
gene_product_mass_group_dict(model::GeckoModel, opt_model)
Extract the mass utilization in mass groups from a solved [`GeckoModel`](@ref).
"""
protein_mass_group_dict(model::GeckoModel, opt_model) =
gene_product_mass_group_dict(model::GeckoModel, opt_model) =
is_solved(opt_model) ?
Dict(
grp[1] => dot(value.(opt_model[:x])[n_reactions(model).+grp[2]], grp[3]) for
grp in model.coupling_row_mass_group
grp.group_id => dot(
value.(opt_model[:x])[n_reactions(model).+grp.gene_product_idxs],
grp.gene_product_molar_masses,
) for grp in model.coupling_row_mass_group
) : nothing
"""
protein_mass_group_dict(model::GeckoModel)
gene_product_mass_group_dict(model::GeckoModel)
A pipe-able variant of [`mass_group_dict`](@ref).
"""
protein_mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
gene_product_mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
"""
protein_mass(model::SMomentModel)
gene_product_mass(model::SMomentModel)
Extract the total mass utilization in a solved [`SMomentModel`](@ref).
"""
protein_mass(model::SMomentModel, opt_model) =
gene_product_mass(model::SMomentModel, opt_model) =
is_solved(opt_model) ?
sum((col.capacity_required for col in model.columns) .* value.(opt_model[:x])) : nothing
"""
protein_mass(model::SMomentModel)
gene_product_mass(model::SMomentModel)
A pipe-able variant of [`protein_mass`](@ref).
A pipe-able variant of [`gene_product_mass`](@ref).
"""
protein_mass(model::SMomentModel) = x -> protein_mass(model, x)
gene_product_mass(model::SMomentModel) = x -> gene_product_mass(model, x)
......@@ -73,7 +73,7 @@ mass of each group of gene products.
function _gecko_mass_group_coupling(model::GeckoModel)
tmp = [ # mm = molar mass, mg = mass group, i = row idx, j = col idx
(i, j, mm) for (i, mg) in enumerate(model.coupling_row_mass_group) for
(j, mm) in zip(mg[2], mg[3])
(j, mm) in zip(mg.gene_product_idxs, mg.gene_product_molar_masses)
]
sparse(
[i for (i, _, _) in tmp],
......@@ -83,17 +83,3 @@ function _gecko_mass_group_coupling(model::GeckoModel)
n_genes(model),
)
end
"""
flux_dict(model::GeckoModel, opt_model)
Returns the fluxes (not gene product concentrations) of the model as a
reaction-keyed dictionary, if solved. See [`protein_dict`](@ref) for a
function to get the gene product concentrations.
"""
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
......@@ -11,9 +11,9 @@
) for (i, grr) in enumerate(reaction_gene_association(model, rid))
) : Isozyme[]
get_gene_product_mass = gid -> get(ecoli_core_protein_masses, gid, 0.0)
get_gene_product_mass = gid -> get(ecoli_core_gene_product_masses, gid, 0.0)
total_protein_mass = 100.0
total_gene_product_mass = 100.0
bounded_model =
model |> with_changed_bounds(
......@@ -27,7 +27,7 @@
reaction_isozymes = get_reaction_isozymes,
gene_product_bounds = g -> g == "b2779" ? (0.01, 0.06) : (0.0, 1.0),
gene_product_molar_mass = get_gene_product_mass,
gene_mass_group_bound = _ -> total_protein_mass,
gene_product_mass_group_bound = _ -> total_gene_product_mass,
)
opt_model = flux_balance_analysis(
......@@ -37,7 +37,7 @@
)
rxn_fluxes = flux_dict(gm, opt_model)
prot_concens = protein_dict(gm, opt_model)
prot_concens = gene_product_dict(gm, opt_model)
@test isapprox(
rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"],
......@@ -45,9 +45,9 @@
atol = TEST_TOLERANCE,
)
prot_mass = sum(ecoli_core_protein_masses[gid] * c for (gid, c) in prot_concens)
prot_mass = sum(ecoli_core_gene_product_masses[gid] * c for (gid, c) in prot_concens)
@test isapprox(prot_mass, total_protein_mass, atol = TEST_TOLERANCE)
@test isapprox(prot_mass, total_gene_product_mass, atol = TEST_TOLERANCE)
end
@testset "GECKO small model" begin
......@@ -95,14 +95,14 @@ end
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)
gene_product_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,
gene_product_mass_group_bound,
)
opt_model = flux_balance_analysis(
......@@ -112,8 +112,8 @@ end
)
rxn_fluxes = flux_dict(gm, opt_model)
gene_products = protein_dict(gm, opt_model)
mass_groups = protein_mass_group_dict(gm, opt_model)
gene_products = gene_product_dict(gm, opt_model)
mass_groups = gene_product_mass_group_dict(gm, opt_model)
@test isapprox(rxn_fluxes["r6"], 3.181818181753438, atol = TEST_TOLERANCE)
@test isapprox(gene_products["g4"], 0.09090909090607537, atol = TEST_TOLERANCE)
......
@testset "SMOMENT" begin
model = load_model(StandardModel, model_paths["e_coli_core.json"])
get_gene_product_mass = gid -> get(ecoli_core_protein_masses, gid, 0.0)
get_gene_product_mass = gid -> get(ecoli_core_gene_product_masses, gid, 0.0)
get_reaction_isozyme =
rid ->
......
Markdown is supported
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