Commit a9677a09 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

switch gecko to gene-based mass groups (seems simpler)

parent 54d9bed3
Pipeline #55040 failed with stages
in 5 minutes and 18 seconds
function make_gecko_model(
model::StandardModel;
reaction_isozymes::Function,
reaction_isozyme_masses::Function,
gene_product_mass::Function,
gene_product_limit::Function,
reaction_mass_group::Function = _ -> "uncategorized",
gene_mass_group::Function = _ -> "uncategorized",
mass_fraction_limit::Function,
)
......@@ -22,27 +22,15 @@ function make_gecko_model(
for i = 1:n_reactions(model)
isozymes = reaction_isozymes(rids[i])
if isempty(isozymes)
push!(columns, _gecko_column(i, 0, 0, 0, lbs[i], ubs[i], [], 0, 0))
push!(columns, _gecko_column(i, 0, 0, 0, lbs[i], ubs[i], [], []))
continue
end
group = reaction_mass_group(rids[i])
mass_group_row =
isnothing(group) ? 0 :
haskey(mass_group_lookup, group) ? mass_group_lookup[group] :
begin
push!(coupling_row_mass_group, group)
mass_group_lookup[group] = length(coupling_row_mass_group)
end
push!(coupling_row_reaction, i)
reaction_coupling_row = length(coupling_row_reaction)
masses = mass_group_row > 0 ? reaction_isozyme_masses(rids[i]) : zeros(length(isozymes))
for (iidx, isozyme) in enumerate(isozymes)
if min(lbs[i], ubs[i]) < 0 && isozyme.kcat_reverse > _constants.tolerance
# reaction can run in reverse
push!(
columns,
_gecko_column(
......@@ -59,13 +47,18 @@ function make_gecko_model(
gene_row_lookup,
coupling_row_gene_product,
),
mass_group_row,
masses[iidx] / isozyme.kcat_reverse,
_gecko_make_mass_group_coupling(
isozyme.gene_product_count,
isozyme.kcat_reverse,
gene_mass_group,
gene_product_mass,
mass_group_lookup,
coupling_row_mass_group,
),
),
)
end
if max(lbs[i], ubs[i]) > 0 && isozyme.kcat_forward > _constants.tolerance
# reaction can run forward
push!(
columns,
_gecko_column(
......@@ -82,8 +75,14 @@ function make_gecko_model(
gene_row_lookup,
coupling_row_gene_product,
),
mass_group_row,
masses[iidx] / isozyme.kcat_forward,
_gecko_make_mass_group_coupling(
isozyme.gene_product_count,
isozyme.kcat_forward,
gene_mass_group,
gene_product_mass,
mass_group_lookup,
coupling_row_mass_group,
),
),
)
end
......@@ -96,7 +95,9 @@ function make_gecko_model(
collect(
zip(coupling_row_gene_product, gene_product_limit.(coupling_row_gene_product)),
),
collect( zip(coupling_row_mass_group, mass_fraction_limit.(coupling_row_mass_group))),
collect(
zip(coupling_row_mass_group, mass_fraction_limit.(coupling_row_mass_group)),
),
model,
)
end
......@@ -119,3 +120,34 @@ _gecko_make_gene_product_coupling(
(row_idx, 1 / kcat)
end for (gene, count) in gene_product_count if haskey(name_lookup, gene)
)
function _gecko_make_mass_group_coupling(
gene_product_count::Dict{String,Int},
kcat::Float64,
gene_row::Function,
gene_product_mass::Function,
row_lookup::Dict{String,Int},
rows::Vector{String},
)
gp_groups = gene_row.(keys(gene_product_count))
gp_mass = gene_product_mass.(keys(gene_product_count))
groups = unique(filter(!isnothing, gp_groups))
group_idx = Dict(groups .=> 1:length(groups))
vals = [0.0 for _ in groups]
for (gpg, mass) in zip(gp_groups, gp_mass)
if !isnothing(gpg)
vals[group_idx[gpg]] += mass / kcat
end
end
collect(
isnothing(group) ? 0 : begin
if !haskey(row_lookup, group)
push!(rows, group)
row_lookup[group] = length(rows)
end
(row_lookup[group], val)
end for (group, val) in zip(groups, vals)
)
end
......@@ -7,8 +7,7 @@ struct _gecko_column
lb::Float64
ub::Float64
gene_product_coupling::Vector{Tuple{Int,Float64}}
mass_group_row::Int
mass_required::Float64
mass_group_coupling::Vector{Tuple{Int,Float64}}
end
struct GeckoModel <: ModelWrapper
......
"""
protein_dict(model::GeckoModel, opt_model)
Return a dictionary mapping protein concentrations to their ids. The argument
`opt_model` is a solved optimization problem, typically returned by
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).
"""
protein_dict(model::GeckoModel, opt_model) = let gids = genes(model)
is_solved(opt_model) ?
Dict(
[gids[gidx] for (gidx,_) = model.coupling_row_gene_product] .=> _gecko_gene_product_coupling(model) * value.(opt_model[:x]),
) : nothing
end
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
"""
protein_dict(model::GeckoModel)
......@@ -21,19 +22,20 @@ protein_dict(model::GeckoModel) = x -> protein_dict(model, x)
"""
mass_group_dict(model::GeckoModel, opt_model)
protein_mass_group_dict(model::GeckoModel, opt_model)
Extract the mass utilization in mass groups from a solved [`GeckoModel`](@ref).
"""
mass_group_dict(model::GeckoModel, opt_model) =
protein_mass_group_dict(model::GeckoModel, opt_model) =
is_solved(opt_model) ?
Dict(
(group for (group,_)=model.coupling_row_mass_group) .=> _gecko_mass_group_coupling(model) * value.(opt_model[:x])
(group for (group, _) in model.coupling_row_mass_group) .=>
_gecko_mass_group_coupling(model) * value.(opt_model[:x]),
) : nothing
"""
mass_group_dict(model::GeckoModel)
protein_mass_group_dict(model::GeckoModel)
A pipe-able variant of [`mass_group_dict`](@ref).
"""
mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
protein_mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
......@@ -68,12 +68,13 @@ _gecko_gene_product_coupling(model::GeckoModel) =
_gecko_mass_group_coupling(model::GeckoModel)
Compute the part of the coupling for [`GeckoModel`](@ref) that limits the total
mass of enzymes available in each group of reactions.
mass of each group of gene products.
"""
_gecko_mass_group_coupling(model::GeckoModel) =
let tmp = [
(col.mass_group_row, i, col.mass_required) for
(i, col) = enumerate(model.columns) if col.direction != 0
let
tmp = [
(row, i, val) for (i, col) = enumerate(model.columns) if col.direction != 0
for (row, val) in col.mass_group_coupling
]
sparse(
[row for (row, _, _) in tmp],
......
......@@ -11,15 +11,7 @@
) for (i, grr) in enumerate(reaction_gene_association(model, rid))
) : Isozyme[]
get_reaction_isozyme_masses =
rid ->
haskey(ecoli_core_protein_stoichiometry, rid) ?
[
sum(
counts .*
get.(Ref(ecoli_core_protein_masses), gids, 0.0),
) for (gids, counts) in zip(reaction_gene_association(model, rid), ecoli_core_protein_stoichiometry[rid])
] : []
get_gene_product_mass = gid -> get(ecoli_core_protein_masses, gid, 0.0)
total_protein_mass = 100.0
......@@ -32,8 +24,8 @@
) |>
with_gecko(
reaction_isozymes = get_reaction_isozymes,
reaction_isozyme_masses = get_reaction_isozyme_masses,
gene_product_limit = g -> g == "b2779" ? (0.01, 0.06) : (0.0, 1.0),
gene_product_mass = get_gene_product_mass,
mass_fraction_limit = _ -> total_protein_mass,
)
......
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