Commit 6663c2f7 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

works partially

parent 308190aa
......@@ -22,14 +22,14 @@ 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, 0, 0))
push!(columns, _gecko_column(i, 0, 0, 0, lbs[i], ubs[i], [], 0, 0))
continue
end
group = reaction_mass_group(rids[i])
mass_group_row =
isnothing(group) ? 0 :
haskey(massgroup_lookup, group) ? mass_group_lookup[group] :
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)
......@@ -38,7 +38,7 @@ function make_gecko_model(
push!(coupling_row_reaction, i)
reaction_coupling_row = length(coupling_row_reaction)
masses = group > 0 ? reaction_isozyme_masses(rids[i]) : zeros(length(isozymes))
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
......@@ -95,9 +95,9 @@ function make_gecko_model(
coupling_row_reaction,
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)),
)model,
),
collect( zip(coupling_row_mass_group, mass_fraction_limit.(coupling_row_mass_group))),
model,
)
end
......
......@@ -14,8 +14,8 @@ end
struct GeckoModel <: ModelWrapper
columns::Vector{_gecko_column}
coupling_row_reaction::Vector{Int}
coupling_row_gene_product::Vector{Tuple{Int,Float64}}
coupling_row_mass_group::Vector{Tuple{String,Float64}} #TODO add to matrices
coupling_row_gene_product::Vector{Tuple{Int,Tuple{Float64,Float64}}}
coupling_row_mass_group::Vector{Tuple{String,Float64}}
inner::MetabolicModel
end
......@@ -121,14 +121,14 @@ function coupling_bounds(model::GeckoModel)
vcat(
iclb,
ilb[model.coupling_row_reaction],
[0.0 for _ in model.coupling_row_gene_product],
[lb for (_, (lb, _)) in model.coupling_row_gene_product],
[0.0 for _ in model.coupling_row_mass_group],
),
vcat(
icub,
rub[model.coupling_row_reaction],
[c for (i, c) in model.coupling_row_gene_product],
[c for (i, c) in model.coupling_row_mass_group],
iub[model.coupling_row_reaction],
[ub for (_, (_, ub)) in model.coupling_row_gene_product],
[ub for (_, ub) in model.coupling_row_mass_group],
),
)
end
......@@ -5,69 +5,35 @@ Return a dictionary mapping protein 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) =
protein_dict(model::GeckoModel, opt_model) = let gids = genes(model)
is_solved(opt_model) ?
Dict(
model.gene_ids .=> value.(opt_model[:x][(length(model.irrev_reaction_ids)+1):end]),
[gids[gidx] for (gidx,_) = model.coupling_row_gene_product] .=> _gecko_gene_product_coupling(model) * value.(opt_model[:x]),
) : nothing
end
"""
protein_dict(model::GeckoModel)
A pipe-able variant of `protein_dict`.
A pipe-able variant of [`protein_dict`](@ref).
"""
protein_dict(model::GeckoModel) = x -> protein_dict(model, x)
"""
get_genes_with_kcats(rid_isozymes::Dict{String, Vector{Isozyme}})
Return all protein (gene ids) that have a kcat from `model` based on `reaction_kcats` field.
Assume that if a reaction has a kcat then each isozyme has a kcat.
"""
function get_genes_with_kcats(rid_isozymes::Dict{String,Vector{Isozyme}})
gids = String[]
for isozymes in values(rid_isozymes)
for isozyme in isozymes
append!(gids, collect(keys(isozyme.stoichiometry)))
end
end
return unique(gids)
end
mass_group_dict(model::GeckoModel, opt_model)
Extract the mass utilization in mass groups from a solved [`GeckoModel`](@ref).
"""
remove_low_expressed_isozymes!(
model::StandardModel,
rid_isozymes = Dict{String, Vector{Isozyme}}()
gid_measurements = Dict(),
)
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])
) : nothing
Modify `rid_isozymes` in place by keeping only the highest expressed isozyme.
"""
function remove_low_expressed_isozymes!(
model::StandardModel,
rid_isozymes = Dict{String,Vector{Isozyme}}(),
gid_measurements = Dict(),
)
for (rid, isozymes) in rid_isozymes
measured_proteins = Float64[]
for isozyme in isozymes
gid_stoich = isozyme.stoichiometry
push!(
measured_proteins,
sum(
map(
*,
collect(values(gid_stoich)),
[get(gid_measurements, gid, 0.0) for gid in keys(gid_stoich)],
[model.genes[gid].molar_mass for gid in keys(gid_stoich)],
),
),
)
end
idx = argmax(measured_proteins)
rid_isozymes[rid] = [rid_isozymes[rid][idx]]
end
mass_group_dict(model::GeckoModel)
return nothing
end
A pipe-able variant of [`mass_group_dict`](@ref).
"""
mass_group_dict(model::GeckoModel) = x -> mass_group_dict(model, x)
......@@ -16,9 +16,9 @@
haskey(ecoli_core_protein_stoichiometry, rid) ?
[
sum(
values(counts) .*
get.(Ref(ecoli_core_protein_masses), keys(counts), 0.0),
) for (iidx, counts) in enumerate(ecoli_core_protein_stoichiometry[rid])
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])
] : []
total_protein_mass = 100.0
......@@ -26,14 +26,14 @@
gm =
model |>
with_changed_bounds(
["EX_glc__D_e", "b2779", "GLCpts"];
lower = [-1000.0, 0.01, -1.0],
upper = [nothing, 0.06, 12.0],
["EX_glc__D_e", "GLCpts"];
lower = [-1000.0, -1.0],
upper = [nothing, 12.0],
) |>
with_gecko(
reaction_isozymes = get_reaction_isozymes,
reaction_isozyme_masses = get_reaction_oisozyme_masses,
gene_product_limit = _ -> 1.0,
reaction_isozyme_masses = get_reaction_isozyme_masses,
gene_product_limit = g -> g == "b2779" ? (0.01, 0.06) : (0.0, 1.0),
mass_fraction_limit = _ -> total_protein_mass,
)
......
......@@ -37,16 +37,16 @@ end
# set up the workers for Distributed, so that the tests that require more
# workers do not unnecessarily load the stuff multiple times
W = addprocs(2)
t = @elapsed @everywhere using COBREXA, Tulip, JuMP
print_timing("import of packages", t)
t = @elapsed @everywhere begin
model = Model(Tulip.Optimizer)
@variable(model, 0 <= x <= 1)
@objective(model, Max, x)
optimize!(model)
end
print_timing("JuMP+Tulip code warmup", t)
#W = addprocs(2)
#t = @elapsed @everywhere using COBREXA, Tulip, JuMP
#print_timing("import of packages", t)
#t = @elapsed @everywhere begin
#model = Model(Tulip.Optimizer)
#@variable(model, 0 <= x <= 1)
#@objective(model, Max, x)
#optimize!(model)
#end
#print_timing("JuMP+Tulip code warmup", t)
# make sure there's a directory for temporary data
tmpdir = "tmpfiles"
......@@ -59,16 +59,17 @@ run_test_file("data_downloaded.jl")
# import base files
@testset "COBREXA test suite" begin
run_test_dir(joinpath("base", "types", "abstract"), "Abstract types")
run_test_dir(joinpath("base", "types"), "Base model types")
run_test_dir(joinpath("base", "logging"), "Logging")
run_test_dir("base", "Base functionality")
run_test_dir(joinpath("base", "utils"), "Utilities")
run_test_dir("io", "I/O functions")
run_test_dir("reconstruction")
run_test_dir("analysis")
run_test_dir(joinpath("analysis", "sampling"), "Sampling")
run_test_file("aqua.jl")
#run_test_dir(joinpath("base", "types", "abstract"), "Abstract types")
#run_test_dir(joinpath("base", "types"), "Base model types")
#run_test_dir(joinpath("base", "logging"), "Logging")
#run_test_dir("base", "Base functionality")
#run_test_dir(joinpath("base", "utils"), "Utilities")
#run_test_dir("io", "I/O functions")
#run_test_dir("reconstruction")
#run_test_dir("analysis")
#run_test_dir(joinpath("analysis", "sampling"), "Sampling")
#run_test_file("analysis/smoment.jl")
run_test_file("analysis/gecko.jl")
end
rmprocs(W)
#rmprocs(W)
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