diff --git a/docs/src/notebooks/7_community_model.jl b/docs/src/notebooks/7_community_model.jl index a5ece1c56103d2921b18e7310533c26350ecb777..3ac595469d75db2a3ea78420102d553daed46fca 100644 --- a/docs/src/notebooks/7_community_model.jl +++ b/docs/src/notebooks/7_community_model.jl @@ -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", ) diff --git a/src/reconstruction/community.jl b/src/reconstruction/community.jl index c5fcb2415d6135ff56bbf68b8687320cf676dd29..4295e841a2385b89a8d38db182472881f30b3197 100644 --- a/src/reconstruction/community.jl +++ b/src/reconstruction/community.jl @@ -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 diff --git a/test/reconstruction/community.jl b/test/reconstruction/community.jl index 449513eebd2501efcb94336e4ddcdd9edbf090be..a2f18a4a830822a2a891e33de22f500b2de5331c 100644 --- a/test/reconstruction/community.jl +++ b/test/reconstruction/community.jl @@ -1,11 +1,9 @@ @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", )