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

finished community stuff

parent 0d629d65
"""
add_objective!(cmodel, objective_mets, objective_weights, objective_column)
add_objective!(community::CoreModel, objective_mets::Vector{String}; objective_weights=Float64[], objective_column_index=0)
Add an objective to the `community` model. Supply the string names of the
objective metabolites in `objective_mets`. Optionally specify the weight to
assign each metabolite in the objective function, if unassigned then equal
weight is assumed. Also, optionally specify whether the objective already exists
in the model by assigning `objective_column_index`. If unassigned then an
objective column will be added, otherwise the column at `objective_column_index`
will be updated.
Note, the weights are negated inside the function so that the objective metabolites
are seen as reagents/substrates, not products in the reaction equation.
# Example
```
add_objective!(model, ["met1", "met2"]) # adds a new column with weights = [1,1]
add_objective!(model, ["met1", "met2"]; objective_weights=[0.1, 0.9]) # adds a new column
add_objective!(model, ["met1", "met2"]; objective_weights=[0.1, 0.9], objective_column_index=10) # updates column 10
```
"""
function add_objective!()
function add_objective!(community::CoreModel, objective_mets::Vector{String}; objective_weights=Float64[], objective_column_index=0)
obj_inds = indexin(objective_mets, metabolites(community))
if isempty(objective_weights)
objective_weights = repeat([1.0], inner=length(objective_mets))
end
if objective_column_index == 0 # needs to be created
nr, _ = size(community.S)
objcol = spzeros(nr)
objcol[obj_inds] .= -objective_weights
# extend model by one reaction
community.S = hcat(community.S, objcol)
community.xl = [community.xl; 0.0]
community.xu = [community.xu; 1000.0]
community.rxns = [community.rxns; "community_biomass"]
community.c = spzeros(size(community.S, 2))
community.c[end] = 1.0
else # only update
community.S[:, objective_column_index] .= 0.0 # reset
community.S[obj_inds, objective_column_index] .= -objective_weights
dropzeros!(community.S)
community.c = spzeros(size(community.S, 2))
community.c[objective_column_index] = 1.0
community.xl[objective_column_index] = 0.0
community.xu[objective_column_index] = 1000.0
end
end
add_objective!(model::CoreModel, objective_met::String, objective_weight::Float64) = add_objective!(model, [objective_met], [objective_weight])
"""
push!()
add_model(
community::CoreModel,
model::M,
exchange_rxn_ids::Vector{String},
exchange_met_ids::Vector{String};
species_name="",
biomass_id=""
) where {M<:MetabolicModel}
Add `model` to `community`, which is a pre-existing community model with
`exchange_rxn_ids` and `exchange_met_ids`. The `species_name` is appended to
each reaction and metabolite, see [`join`](@ref). If `biomass_id` is specified
then a biomass metabolite for `model` is also added to the resulting model. The
column corresponding to the `biomass_id` reaction then produces this new biomass
metabolite with unit coefficient. Note, `exchange_rxn_ids` and
`exchange_met_ids` must already exist in the `community` model.
Append `model` to `cmodel` where `cmodel` is a pre-existing community model with `exchange_rxn_ids` and
`exchange_met_ids`. If an objective function has already been assigned then supply its column index in `objective_col`
and the metabolites used by the objective in `objective_rows` as well as the weight to assign the new
# Example
```
community = add_model(community, model, exchange_rxn_ids, exchange_met_ids; species_name="species_2", biomass_id="BIOMASS_Ecoli_core_w_GAM")
```
"""
function Base.push!(
cmodel::CoreModel,
function add_model(
community::CoreModel,
model::M,
exchange_met_ids::Vector{String},
exchange_rxn_ids::Vector{String},
species_name,
has_biomass_objective;
biomass_id = "",
) where {M<:MetabolicModel}
exchange_met_ids::Vector{String};
species_name="",
biomass_id=""
) where {M<:MetabolicModel}
if has_biomass_objective && biomass_id == ""
exchange_met_community_inds = indexin(exchange_met_ids, metabolites(community))
exchange_rxn_community_inds = indexin(exchange_rxn_ids, reactions(community))
if any(isnothing.(exchange_met_community_inds)) || any(isnothing.(exchange_rxn_community_inds))
throw(
DomainError(
"Argument required.",
"The community uses a biomass objective function, please supply the objective id of the model you want to add.",
"exchange metabolite/reaction not found.",
"Exchange metabolites/reactions must already be contained in the community model.",
),
)
end
n_cmodel_rows, n_cmodel_cols = size(stoichiometry(cmodel))
n_cmodel_rows, n_cmodel_cols = size(stoichiometry(community))
n_model_rows, n_model_cols = size(stoichiometry(model))
Iadd, Jadd, Vadd = findnz(stoichiometry(model))
add_row = has_biomass_objective ? 1 : 0
n_metabolites_total =
n_reactions_total =
# shift to fit into bigger model
Iadd .+= n_cmodel_rows
Jadd .+= n_cmodel_cols
exchange_met_community_inds = indexin(exchange_met_ids, metabolites(cmodel))
if any(isnothing.(exchange_met_community_inds))
throw(
DomainError(
"Exchange metabolite not found.",
"Exchange metabolite not found in community model.",
),
)
end
exchange_rxn_model_inds = indexin(exchange_rxn_ids, reactions(model))
# shift to fit into community
Iadd .+= n_cmodel_rows
Jadd .+= n_cmodel_cols
# when adding a single model not that many reactions, push! okay?
for i = 1:length(exchange_met_community_inds)
exchange_rxn_model_inds = indexin(exchange_rxn_ids, reactions(model))
for i = 1:length(exchange_rxn_ids)
isnothing(exchange_rxn_model_inds[i]) && continue
push!(Iadd, n_cmodel_rows + exchange_met_community_inds[i])
push!(Jadd, n_cmodel_cols + exchange_rxn_model_inds[i])
push!(Iadd, exchange_met_community_inds[i]) # already present ex met in community model
push!(Jadd, n_cmodel_cols + exchange_rxn_model_inds[i]) # new column hence the offset
push!(Vadd, 1.0)
end
# if has_biomass_objective
# biomass_ind = first(indexin([biomass_id]), reactions(model))
biomass_met = 0.0
if biomass_id != "" # add biomass metabolite
biomass_rxn = first(indexin([biomass_id], reactions(model)))
push!(Iadd, n_model_rows + n_cmodel_rows + 1)
push!(Jadd, biomass_rxn + n_cmodel_cols)
push!(Vadd, 1.0)
biomass_met = 1
end
# push!(Iadd, n_cmodel_rows + )
# push!(Jadd, n_cmodel_cols + )
# push!(Vadd, 1.0)
# end
n_metabolites_total = n_model_rows + n_cmodel_rows + biomass_met
n_reactions_total = n_cmodel_cols + n_model_cols
I, J, V = findnz(stoichiometry(cmodel))
I, J, V = findnz(stoichiometry(community))
I = [I; Iadd]
J = [J; Jadd]
V = [V; Vadd]
S = sparse(I, J, V, n_metabolites_total, n_reactions_total)
lbsadd, ubsadd = bounds(model)
lbs, ubs = bounds(community)
lbs = [lbs; lbsadd]
ubs = [lbs; ubsadd]
ubs = [ubs; ubsadd]
rxnsadd = "$(species_name)_".reactions(model)
if has_biomass_objective
rxnsadd = "$(species_name)_".*reactions(model)
if biomass_id != ""
metsadd =
["$(species_name)_" .* metabolites(model); "$(species_name)_" * biomass_id]
else
metsadd = "$(species_name)_".metabolites(model)
metsadd = "$(species_name)_".*metabolites(model)
end
rxns = [rxns; rxnsadd]
mets = [mets; metsadd]
rxns = [reactions(community); rxnsadd]
mets = [metabolites(community); metsadd]
return CoreModel(S, spzeros(size(S, 1)), spzeros(size(S, 2)), lbs, ubs, rxns, mets)
# adds to the original community data, could possibly reset?
I, V = findnz(balance(community))
b = sparsevec(I, V, n_metabolites_total)
I, V = findnz(objective(community))
c = sparsevec(I, V, n_reactions_total)
return CoreModel(S, b, c, lbs, ubs, rxns, mets)
end
"""
......@@ -176,6 +234,15 @@ function Base.join(
)
end
if length(exchange_met_ids) != length(exchange_rxn_ids)
throw(
DomainError(
"Exchange identifiers are misspecified",
"The lenghts of the exchange metabolite and reactions are different.",
),
)
end
# get offsets to construct community S
reaction_lengths = [n_reactions(model) for model in models]
metabolite_lengths = [n_metabolites(model) for model in models]
......
......@@ -27,15 +27,8 @@
community.xl[env_ex_inds] .= m2.xl[m2_ex_inds]
community.xu[env_ex_inds] .= m2.xu[m2_ex_inds]
biomass_metabolite_inds = indexin(
["species_1_BIOMASS_Ecoli_core_w_GAM", "species_2_BIOMASS_Ecoli_core_w_GAM"],
metabolites(community),
)
community.S[biomass_metabolite_inds, end] .= -1.0
community.c[end] = 1.0
community.xl[end] = 0.0
community.xu[end] = 1000.0
biomass_ids = ["species_1_BIOMASS_Ecoli_core_w_GAM", "species_2_BIOMASS_Ecoli_core_w_GAM"]
add_objective!(community, biomass_ids; objective_column_index=first(indexin(["community_biomass"], reactions(community))))
d = flux_balance_analysis_dict(community, Tulip.Optimizer)
@test size(stoichiometry(community)) == (166, 211)
......@@ -105,3 +98,42 @@ end
@test size(stoichiometry(community)) == (2203, 3003)
@test isapprox(d["community_biomass"], 0.8739215069675402, atol = TEST_TOLERANCE)
end
@testset "Community model modifications" begin
model_path = download_data_file(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
joinpath("data", "e_coli_core.json"),
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)
m1 = load_model(CoreModel, model_path)
boundary_rxn_ids, boundary_met_ids = all_boundaries(m1)
exchange_rxn_ids = filter(startswith("EX_"), boundary_rxn_ids)
exchange_met_ids = filter(endswith("_e"), boundary_met_ids)
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM"]
community = COBREXA.join(
[m1],
exchange_rxn_ids,
exchange_met_ids;
add_biomass_objective = true,
biomass_ids = biomass_ids,
)
env_ex_inds = indexin(exchange_rxn_ids, reactions(community))
m1_ex_inds = indexin(exchange_rxn_ids, reactions(m1))
community.xl[env_ex_inds] .= m1.xl[m1_ex_inds]
community.xu[env_ex_inds] .= m1.xu[m1_ex_inds]
m2 = load_model(CoreModel, model_path)
community = add_model(community, m2, exchange_rxn_ids, exchange_met_ids; species_name="species_2", biomass_id="BIOMASS_Ecoli_core_w_GAM")
biomass_ids = ["species_1_BIOMASS_Ecoli_core_w_GAM", "species_2_BIOMASS_Ecoli_core_w_GAM"]
add_objective!(community, biomass_ids; objective_column_index=first(indexin(["community_biomass"], reactions(community))))
d = flux_balance_analysis_dict(community, Tulip.Optimizer)
@test size(stoichiometry(community)) == (166, 211)
@test isapprox(d["community_biomass"], 0.41559777495618294, atol = TEST_TOLERANCE)
end
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