Skip to content
Snippets Groups Projects
Unverified Commit 6c97da9b authored by St. Elmo's avatar St. Elmo
Browse files

Update community construction functions to take a dictionary

parent 8705af8d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
docs/src/notebooks/7_community_model.jl
+ 9
− 10
View file @ 6c97da9b
......@@ -40,15 +40,15 @@ sol["BIOMASS_Ecoli_core_w_GAM"] # Enolase knockout μ, cannot grow by itself
# ## Build a community model of the cytochrome oxidase knockout and the ATP synthase knockout models
ex_rxns = filter(looks_like_exchange_reaction, reactions(base_model)) # identify exchange reactions heuristically
ex_mets = [first(keys(reaction_stoichiometry(base_model, ex_rxn))) for ex_rxn in ex_rxns] # identify exchange metabolites IN THE SAME ORDER as ex_rxns
[ex_rxns ex_mets]
ex_rxn_mets = Dict(
ex_rxn => first(keys(reaction_stoichiometry(base_model, ex_rxn))) for
ex_rxn in filter(looks_like_exchange_reaction, reactions(base_model))
) # identify exchange reactions heuristically
#
model_names = ["cytbd_ko", "atps4r_ko"]
community_model = join_with_exchanges(
[cytbd_knockout_model, atps4r_knockout_model],
ex_rxns,
ex_mets;
ex_rxn_mets;
add_biomass_objective = true,
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"],
model_names = model_names,
......@@ -56,9 +56,9 @@ community_model = join_with_exchanges(
# ## Set exchange reaction bounds of community model based on the bounds of the individual models
env_ex_rxn_idxs = indexin(ex_rxns, reactions(community_model)) # identify the global (environmental exchange reactions)
cytbd_ex_rxn_idxs = indexin(ex_rxns, reactions(cytbd_knockout_model)) # identify the indices of the corresponding exchange reactions in the original models
atps4r_ex_rxn_idxs = indexin(ex_rxns, reactions(atps4r_knockout_model))
env_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(community_model)) # identify the global (environmental exchange reactions)
cytbd_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(cytbd_knockout_model)) # identify the indices of the corresponding exchange reactions in the original models
atps4r_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(atps4r_knockout_model))
# In case some exchange reactions are not present in both models, set
# environmental exchange bound to the sum of the individual exchange bounds
......@@ -96,8 +96,7 @@ d["community_biomass"] # community μ
community_model = add_model_with_exchanges(
community_model,
eno_knockout_model,
ex_rxns,
ex_mets;
ex_rxn_mets;
model_name = "eno_ko",
biomass_id = "BIOMASS_Ecoli_core_w_GAM",
)
......
src/reconstruction/community.jl
+ 52
− 52
View file @ 6c97da9b
......@@ -63,35 +63,38 @@ end
add_model_with_exchanges(
community::CoreModel,
model::MetabolicModel,
exchange_rxn_ids::Vector{String},
exchange_met_ids::Vector{String};
exchange_rxn_mets::Dict{String,String};
model_name = "unknown_species",
biomass_id = nothing,
)::CoreModel
Add `model` to `community`, which is a pre-existing community model with
`exchange_rxn_ids` and `exchange_met_ids`. The `model_name` is appended to
each reaction and metabolite, see [`join_with_exchanges`](@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.
Add `model` to `community`, which is a pre-existing community model with exchange reactions
and metabolites in the dictionary `exchange_rxn_mets`. The `model_name` is appended to each
reaction and metabolite, see [`join_with_exchanges`](@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.
# Example
```
community = add_model_with_exchanges(community, model, exchange_rxn_ids, exchange_met_ids; model_name="species_2", biomass_id="BIOMASS_Ecoli_core_w_GAM")
community = add_model_with_exchanges(community,
model,
exchange_rxn_mets;
model_name="species_2",
biomass_id="BIOMASS_Ecoli_core_w_GAM")
```
"""
function add_model_with_exchanges(
community::CoreModel,
model::MetabolicModel,
exchange_rxn_ids::Vector{String},
exchange_met_ids::Vector{String};
exchange_rxn_mets::Dict{String,String};
model_name = "unknown_species",
biomass_id = nothing,
)::CoreModel
exchange_rxn_ids = keys(exchange_rxn_mets)
exchange_met_ids = values(exchange_rxn_mets)
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)) ||
......@@ -171,36 +174,33 @@ end
"""
join_with_exchanges(models::Vector{M},
exchange_rxn_ids::Vector{String},
exchange_met_ids::Vector{String};
exchange_rxn_mets::Dict{String, String};
add_biomass_objective=false,
biomass_ids::Vector{String},
model_names=String[]
)
Return a `CoreModel` representing the community model of `models` joined through
their `exchange_rxn_ids` and `exchange_met_ids`. These exchange reactions and
metabolites link to environmental metabolites and reactions.
Optionally specify `model_names` to append a specific name to each reaction
and metabolite of an organism for easier reference (default is `species_i` for
each model index i in `models`). Note, the bounds of the environmental variables
are all set to zero. Thus, to run a simulation you need to constrain them
appropriately. All the other bounds are inherited from the models used to
construct the community model.
If `add_biomass_objective` is true then `biomass_ids` needs to be supplied as
well. This creates a model with an extra reaction added to the end of the
stoichiometric matrix (last column) that can be assigned as the objective
reaction. It also creates biomass "metabolites" that can be used in this
objective reaction. Note, this reaction is unspecified, further action needs to
be taken to specify it, e.g. assign weights to the last column of the
stoichiometric matrix in the rows corresponding to the biomass metabolites.
To further clarify how this `join` works. Suppose you have 2 organisms with
stoichiometric matrices S₁ and S₂ and you want to link them with
`exchange_met_ids = [em₁, em₂, em₃, ...]` and `exchange_rxn_ids = [er₁, er₂,
er₃, ...]`. Then a new community stoichiometric matrix is constructed that looks
like this:
Return a `CoreModel` representing the community model of `models` joined through their
exchange reactions and metabolites in the dictionary `exchange_rxn_mets`, which maps
exchange reactions to their associated metabolite. These exchange reactions and metabolites
link model metabolites to environmental metabolites and reactions. Optionally specify
`model_names` to append a specific name to each reaction and metabolite of an organism for
easier reference (default is `species_i` for each model index i in `models`). Note, the
bounds of the environmental variables are all set to zero. Thus, to run a simulation you
need to constrain them appropriately. All the other bounds are inherited from the models
used to construct the community model.
If `add_biomass_objective` is true then `biomass_ids` needs to be supplied as well. This
creates a model with an extra reaction added to the end of the stoichiometric matrix (last
column) that can be assigned as the objective reaction. It also creates biomass
"metabolites" that can be used in this objective reaction. Note, this reaction is
unspecified, further action needs to be taken to specify it, e.g. assign weights to the last
column of the stoichiometric matrix in the rows corresponding to the biomass metabolites.
To further clarify how this `join` works. Suppose you have 2 organisms with stoichiometric
matrices S₁ and S₂ and you want to link them with `exchange_rxn_mets = Dict(er₁ => em₁, er₂
=> em₂, er₃ => em₃, ...)`. Then a new community stoichiometric matrix is constructed that
looks like this:
```
_ er₁ er₂ er₃ ... b_
|S₁ |
......@@ -213,13 +213,12 @@ S = em₂| |
bm₂|_ _|
```
The exchange reactions in each model get linked to environmental metabolites,
`emᵢ`, and these get linked to environmental exchanges, `erᵢ`. These `erᵢ`
behave like normal single organism exchange reactions. When
`add_biomass_objective` is true each model's biomass becomes a pseudo-metabolite
(`bmᵢ`). These can be weighted in column `b`, called the `community_biomass`
reaction in the community model, if desired. Refer to the tutorial if this is
unclear.
The exchange reactions in each model get linked to environmental metabolites, `emᵢ`, and
these get linked to environmental exchanges, `erᵢ`. These `erᵢ` behave like normal single
organism exchange reactions. When `add_biomass_objective` is true each model's biomass
becomes a pseudo-metabolite (`bmᵢ`). These can be weighted in column `b`, called the
`community_biomass` reaction in the community model, if desired. Refer to the tutorial if
this is unclear.
# Example
```
......@@ -227,15 +226,14 @@ m1 = load_model(core_model_path)
m2 = load_model(CoreModel, core_model_path)
# need to list ALL the exchanges that will form part of the entire model
exchange_rxn_ids = filter(looks_like_exchange_reaction, reactions(m1))
exchange_met_ids = [first(keys(reaction_stoichiometry(m1, ex_rxn))) for ex_rxn in exchange_rxn_ids]
exchange_rxn_mets = Dict(k => first(keys(reaction_stoichiometry(m1, ex_rxn)))
for filter(looks_like_exchange_reaction, reactions(m1)))
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"]
community = join_with_exchanges(
[m1, m2],
exchange_rxn_ids,
exchange_met_ids;
exchange_rxn_mets;
add_biomass_objective = true,
biomass_ids = biomass_ids,
)
......@@ -243,18 +241,20 @@ community = join_with_exchanges(
"""
function join_with_exchanges(
models::Vector{M},
exchange_rxn_ids::Vector{String},
exchange_met_ids::Vector{String};
exchange_rxn_mets::Dict{String,String};
add_biomass_objective = true,
biomass_ids = String[],
model_names = String[],
)::CoreModel where {M<:MetabolicModel}
exchange_rxn_ids = keys(exchange_rxn_mets)
exchange_met_ids = values(exchange_rxn_mets)
if add_biomass_objective && isempty(biomass_ids)
throw(
DomainError(
"biomass_ids",
"Please add supply the string ids of the biomass functions when `add_biomass_objective` is true.",
"Please add supply the string ids of the biomass functions when
`add_biomass_objective` is true.",
),
)
end
......
test/reconstruction/community.jl
+ 26
− 30
View file @ 6c97da9b
@testset "Detailed community stoichiometrix matrix check" begin
m1 = test_toyModel()
m2 = test_toyModel()
ex_rxn_ids = ["EX_m1(e)", "EX_m3(e)"]
ex_met_ids = ["m1[e]", "m3[e]"]
ex_rxn_mets = Dict("EX_m1(e)" => "m1[e]", "EX_m3(e)" => "m3[e]")
c1 =
join_with_exchanges([m1, m2], ex_rxn_ids, ex_met_ids; add_biomass_objective = false)
c1 = join_with_exchanges([m1, m2], ex_rxn_mets; add_biomass_objective = false)
# test of stoichs are the same
@test all(c1.S[1:6, 1:7] .== c1.S[7:12, 8:14])
......@@ -28,8 +26,7 @@
c2 = join_with_exchanges(
[m1, m2],
ex_rxn_ids,
ex_met_ids;
ex_rxn_mets;
add_biomass_objective = true,
biomass_ids = ["biomass1", "biomass1"],
)
......@@ -54,22 +51,22 @@ end
m1 = load_model(model_paths["e_coli_core.json"])
m2 = load_model(CoreModel, model_paths["e_coli_core.json"])
exchange_rxn_ids = filter(looks_like_exchange_reaction, reactions(m2))
exchange_met_ids =
[first(keys(reaction_stoichiometry(m2, ex_rxn))) for ex_rxn in exchange_rxn_ids]
exchange_rxn_mets = Dict(
ex_rxn => first(keys(reaction_stoichiometry(m2, ex_rxn))) for
ex_rxn in filter(looks_like_exchange_reaction, reactions(m2))
)
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"]
community = join_with_exchanges(
[m1, m2],
exchange_rxn_ids,
exchange_met_ids;
exchange_rxn_mets;
add_biomass_objective = true,
biomass_ids = biomass_ids,
)
env_ex_inds = indexin(exchange_rxn_ids, reactions(community))
m2_ex_inds = indexin(exchange_rxn_ids, reactions(m2))
env_ex_inds = indexin(keys(exchange_rxn_mets), reactions(community))
m2_ex_inds = indexin(keys(exchange_rxn_mets), reactions(m2))
community.xl[env_ex_inds] .= m2.xl[m2_ex_inds]
community.xu[env_ex_inds] .= m2.xu[m2_ex_inds]
......@@ -92,23 +89,23 @@ end
m1 = load_model(CoreModel, model_paths["e_coli_core.json"])
m2 = load_model(CoreModel, model_paths["iJO1366.mat"])
exchange_rxn_ids = filter(looks_like_exchange_reaction, reactions(m2))
exchange_met_ids =
[first(keys(reaction_stoichiometry(m2, ex_rxn))) for ex_rxn in exchange_rxn_ids]
exchange_rxn_mets = Dict(
ex_rxn => first(keys(reaction_stoichiometry(m2, ex_rxn))) for
ex_rxn in filter(looks_like_exchange_reaction, reactions(m2))
)
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ec_iJO1366_core_53p95M"]
community = join_with_exchanges(
[m1, m2],
exchange_rxn_ids,
exchange_met_ids;
exchange_rxn_mets;
add_biomass_objective = true,
biomass_ids = biomass_ids,
)
env_ex_inds = indexin(exchange_rxn_ids, reactions(community))
m2_ex_inds = indexin(exchange_rxn_ids, reactions(m2))
m1_ex_inds = indexin(exchange_rxn_ids, reactions(m1))
env_ex_inds = indexin(keys(exchange_rxn_mets), reactions(community))
m2_ex_inds = indexin(keys(exchange_rxn_mets), reactions(m2))
m1_ex_inds = indexin(keys(exchange_rxn_mets), reactions(m1))
for (env_ex, m2_ex, m1_ex) in zip(env_ex_inds, m2_ex_inds, m1_ex_inds)
m2lb = isnothing(m2_ex) ? 0.0 : m2.xl[m2_ex]
......@@ -144,22 +141,22 @@ end
@testset "Community model modifications" begin
m1 = load_model(CoreModel, model_paths["e_coli_core.json"])
exchange_rxn_ids = filter(looks_like_exchange_reaction, reactions(m1))
exchange_met_ids =
[first(keys(reaction_stoichiometry(m1, ex_rxn))) for ex_rxn in exchange_rxn_ids]
exchange_rxn_mets = Dict(
ex_rxn => first(keys(reaction_stoichiometry(m1, ex_rxn))) for
ex_rxn in filter(looks_like_exchange_reaction, reactions(m1))
)
biomass_ids = ["BIOMASS_Ecoli_core_w_GAM"]
community = join_with_exchanges(
[m1],
exchange_rxn_ids,
exchange_met_ids;
exchange_rxn_mets;
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))
env_ex_inds = indexin(keys(exchange_rxn_mets), reactions(community))
m1_ex_inds = indexin(keys(exchange_rxn_mets), reactions(m1))
community.xl[env_ex_inds] .= m1.xl[m1_ex_inds]
community.xu[env_ex_inds] .= m1.xu[m1_ex_inds]
......@@ -168,8 +165,7 @@ end
community = add_model_with_exchanges(
community,
m2,
exchange_rxn_ids,
exchange_met_ids;
exchange_rxn_mets;
model_name = "species_2",
biomass_id = "BIOMASS_Ecoli_core_w_GAM",
)
......
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