Unverified Commit d9dab2a1 authored by Laurent Heirendt's avatar Laurent Heirendt Committed by GitHub
Browse files

Merge pull request #111 from LCSB-BioCore/mo-fix-standardmodel-interface

Add generic interface for StandardModels
parents 8e20d36d c8ab3c18
......@@ -21,6 +21,7 @@ MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79"
OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
......@@ -34,6 +35,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
[compat]
julia = "1"
DistributedData = ">=0.1.3"
SBML = ">=0.4"
julia = "1"
......@@ -16,6 +16,7 @@ using Measurements
using Statistics
using Random
using MacroTools # for DSL :)
using OrderedCollections
import Base: findfirst, getindex, show
import Pkg
......
......@@ -165,7 +165,7 @@ function flux_variability_analysis(
λmin <= COBREXA.JuMP.objective_function(opt_model) <= λmax # in case there is a negative bound
)
for i = 1:length(v)
for (i, r_id) in enumerate(reactions(model))
COBREXA.JuMP.@objective(opt_model, Max, v[i])
COBREXA.JuMP.optimize!(opt_model)
status = (
......@@ -173,11 +173,10 @@ function flux_variability_analysis(
COBREXA.JuMP.termination_status(opt_model) == MOI.LOCALLY_SOLVED
)
if status
fva_max[model.reactions[i].id] =
Dict(zip(reactions(model), value.(opt_model[:x])))
fva_max[r_id] = Dict(zip(reactions(model), value.(opt_model[:x])))
else
@warn "Error maximizing index: $i with error $(termination_status(opt_model))"
fva_max[model.reactions[i].id] = nothing
fva_max[r_id] = nothing
end
@objective(opt_model, Min, v[i])
......@@ -187,11 +186,10 @@ function flux_variability_analysis(
termination_status(opt_model) == MOI.LOCALLY_SOLVED
)
if status
fva_min[model.reactions[i].id] =
Dict(zip(reactions(model), value.(opt_model[:x])))
fva_min[r_id] = Dict(zip(reactions(model), value.(opt_model[:x])))
else
@warn "Error minimizing index: $i with error $(termination_status(opt_model))"
fva_min[model.reactions[i].id] = nothing
fva_min[r_id] = nothing
end
end
......
function _get_boundary_reaction_ids(model::StandardModel)::Array{String,1}
return [i.id for i in model.reactions if length(i.metabolites) == 1]
return [r.id for r in values(model.reactions) if is_boundary(r)]
end
function add_cycle_free(fluxes::Dict{String,Float64})
......@@ -8,7 +8,7 @@ function add_cycle_free(fluxes::Dict{String,Float64})
old_objective = objective_function(opt_model)
boundary_ids = _get_boundary_reaction_ids(model)
min_objectives = zeros(Int64, 1, length(v))
for (i, reaction) in enumerate(model.reactions)
for (i, reaction) in enumerate(values(model.reactions))
id = reaction.id
if v[i] in old_objective.terms.keys
continue
......
......@@ -8,7 +8,7 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
modelid = modeldict["id"]
mets = Metabolite[]
mets = OrderedDict{String,Metabolite}()
for i in modeldict["metabolites"]
met = Metabolite()
......@@ -28,7 +28,7 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
elseif k == "annotation"
for (kk, vv) in v
if typeof(vv) == String
met.annotation[kk] = vv
met.annotation[kk] = [vv]
else
met.annotation[kk] = convert(Vector{String}, vv)
end
......@@ -37,10 +37,10 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
@warn "Unrecognized metabolite field: $k"
end
end
push!(mets, met)
mets[met.id] = met
end
genes = Gene[]
genes = OrderedDict{String,Gene}()
for i in modeldict["genes"]
gene = Gene()
for (k, v) in i
......@@ -55,7 +55,7 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
elseif k == "annotation"
for (kk, vv) in v
if typeof(vv) == String
gene.annotation[kk] = vv
gene.annotation[kk] = [vv]
else
gene.annotation[kk] = convert(Vector{String}, vv)
end
......@@ -64,10 +64,10 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
@warn "Unrecognized gene field: $k"
end
end
push!(genes, gene)
genes[gene.id] = gene
end
rxns = Reaction[]
rxns = OrderedDict{String,Reaction}()
for i in modeldict["reactions"]
rxn = Reaction()
for (k, v) in i
......@@ -77,20 +77,14 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
rxn.name = v
elseif k == "metabolites"
for (kk, vv) in v
ind = findfirst(x -> x.id == kk, mets)
if isnothing(ind)
@warn "Metabolite $kk not found in reaction assignment."
continue
else
rxn.metabolites[mets[ind]] = vv
end
rxn.metabolites[mets[kk].id] = vv
end
elseif k == "lower_bound"
rxn.lb = v
elseif k == "upper_bound"
rxn.ub = v
elseif k == "gene_reaction_rule"
rxn.grr = _parse_grr(v, genes)
rxn.grr = _parse_grr(v)
elseif k == "subsystem"
rxn.subsystem = v
elseif k == "notes"
......@@ -100,7 +94,7 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
elseif k == "annotation"
for (kk, vv) in v
if typeof(vv) == String
rxn.annotation[kk] = vv
rxn.annotation[kk] = [vv]
else
rxn.annotation[kk] = convert(Vector{String}, vv)
end
......@@ -111,8 +105,8 @@ function _read_model(file_location::String, ::Type{JSONFile}, ::Type{StandardMod
@warn "Unrecognized reaction field: $k"
end
end
push!(rxns, rxn)
rxns[rxn.id] = rxn
end
return StandardModel(modelid, rxns, mets, genes)
return StandardModel(modelid; reactions = rxns, metabolites = mets, genes = genes)
end
......@@ -14,7 +14,7 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
model_id = haskey(modeldict, "description") ? modeldict["description"] : model_name
model_id = haskey(modeldict, "modelName") ? modeldict["modelName"] : model_name # more specific
mets = Metabolite[]
mets = OrderedDict{String,Metabolite}()
for i in eachindex(modeldict["mets"])
met = Metabolite()
......@@ -58,17 +58,17 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
end
if haskey(modeldict, "metSBOTerms")
met.annotation["sbo"] = string(modeldict["metSBOTerms"][i])
met.annotation["sbo"] = [string(modeldict["metSBOTerms"][i])]
end
if haskey(modeldict, "metNotes")
met.notes["note"] = string.(split(string(modeldict["metNotes"][i]), "; "))
end
push!(mets, met)
mets[met.id] = met
end
genes = Gene[]
genes = OrderedDict{String,Gene}()
for i in eachindex(modeldict["genes"])
gene = Gene()
if haskey(modeldict, "genes")
......@@ -80,25 +80,23 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
# name and other fields don't generally exist in the matlab models,
# ignoring them for now
push!(genes, gene)
genes[gene.id] = gene
end
rxns = Reaction[]
rxns = OrderedDict{String,Reaction}()
for i in eachindex(modeldict["rxns"])
rxn = Reaction()
if haskey(modeldict, "rxns")
rxn.id = modeldict["rxns"][i]
else
continue
end
rxn.id = modeldict["rxns"][i]
if haskey(modeldict, "rxnNames")
rxn.name = modeldict["rxnNames"][i]
end
metinds = findall(x -> x .!= 0.0, modeldict["S"][:, i])
rxn.metabolites =
Dict{Metabolite,Float64}(mets[j] => modeldict["S"][j, i] for j in metinds)
rxn.metabolites = Dict{String,Float64}(
mets[sj].id => modeldict["S"][j, i] for
(sj, j) in zip(modeldict["mets"][metinds], metinds)
)
if haskey(modeldict, "lb")
rxn.lb = modeldict["lb"][i]
......@@ -110,7 +108,7 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
if haskey(modeldict, "grRules")
grr_string = modeldict["grRules"][i]
rxn.grr = _parse_grr(grr_string, genes)
rxn.grr = _parse_grr(grr_string)
end
rxn.subsystem = join(modeldict["subSystems"][i], "; ")
......@@ -133,7 +131,7 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
end
if haskey(modeldict, "rxnSBOTerms")
rxn.annotation["sbo"] = string(modeldict["rxnSBOTerms"][i])
rxn.annotation["sbo"] = [string(modeldict["rxnSBOTerms"][i])]
end
# look for some notes
......@@ -141,10 +139,10 @@ function _read_model(file_location::String, ::Type{MFile}, ::Type{StandardModel}
rxn.notes["note"] = string.(split(string(modeldict["rxnNotes"][i]), "; "))
end
push!(rxns, rxn)
rxns[rxn.id] = rxn
end
return StandardModel(model_id, rxns, mets, genes)
return StandardModel(model_id; reactions = rxns, metabolites = mets, genes = genes)
end
function _read_model(file_location::String, ::Type{MFile}, ::Type{CoreModel})
......
......@@ -2,8 +2,8 @@
Pretty printing of `Gene`.
"""
function Base.show(io::IO, ::MIME"text/plain", g::Gene)
println(io, "Gene ID: ", g.id)
println(io, "Gene name: ", g.name)
println(io, "Gene ID: $(g.id)")
println(io, "Gene name: $(g.name)")
end
"""
......@@ -12,14 +12,3 @@ Pretty printing of `Vector{Gene}`.
function Base.show(io::IO, ::MIME"text/plain", gs::Vector{Gene})
println(io, "Gene set of length: ", length(gs))
end
"""
Pretty printing of gene reaction rules in type `Vector{Vector{Gene}}`.
"""
function Base.show(io::IO, ::MIME"text/plain", grr::Vector{Vector{Gene}})
grr_strings = String[]
for gr in grr
push!(grr_strings, "(" * join([g.id for g in gr], " and ") * ")")
end
println(io, join(grr_strings, " or "))
end
......@@ -5,17 +5,7 @@ Pretty printing of metabolite::Metabolite.
function Base.show(io::IO, ::MIME"text/plain", m::Metabolite)
println(
io,
"Metabolite ID: ",
m.id,
"\n",
"Metabolite name: ",
m.name,
"\n",
"Formula: ",
m.formula,
"\n",
"Charge: ",
m.charge,
"Metabolite ID: $(m.id)\nMetabolite name: $(m.name)\nFormula: $(m.formula)\nCharge: $(m.charge)",
)
end
......
"""
Pretty printing of reaction::Reaction.
"""
......@@ -16,17 +15,17 @@ function Base.show(io::IO, ::MIME"text/plain", r::Reaction)
products = String[]
for (k, v) in r.metabolites
if v < 0.0
push!(substrates, string(abs(v)) * " " * k.id)
push!(substrates, string(abs(v)) * " " * k)
else
push!(products, string(abs(v)) * " " * k.id)
push!(products, string(abs(v)) * " " * k)
end
end
isempty(substrates) && (substrates = "∅")
isempty(products) && (products = "∅")
println(io, "Reaction ID: ", r.id)
println(io, "Reaction name: ", r.name)
println(io, "Reaction subsystem: ", r.subsystem)
println(io, "Reaction ID: $(r.id)")
println(io, "Reaction name: $(r.name)")
println(io, "Reaction subsystem: $(r.subsystem)")
if length(substrates) > 5 && length(products) > 5
sp = substrates[1] * " + ... + " * substrates[end]
pp = products[1] * " + ... + " * products[end]
......@@ -40,16 +39,16 @@ function Base.show(io::IO, ::MIME"text/plain", r::Reaction)
else
println(io, join(substrates, " + ") * arrow * join(products, " + "))
end
println(io, "Lower bound: ", r.lb)
println(io, "Upper bound: ", r.ub)
println(io, "Lower bound: $(r.lb)")
println(io, "Upper bound: $(r.ub)")
grr_strings = String[]
for gr in r.grr
push!(grr_strings, "(" * join([g.id for g in gr], " and ") * ")")
push!(grr_strings, "(" * join([g for g in gr], " and ") * ")")
end
grr_string = join(grr_strings, " or ")
(isnothing(grr_string) || grr_string == "") && (grr_string = "")
println(io, "Genes: ", grr_string)
println(io, "Genes: $grr_string")
println(io, "E.C. number: ", join(get(r.annotation, "ec-code", [""]), " or "))
end
......
......@@ -5,16 +5,6 @@ Pretty printing of model::StandardModel.
function Base.show(io::IO, ::MIME"text/plain", m::StandardModel)
println(
io,
"Constraint based model: ",
m.id,
"\n",
"Number of reactions: ",
length(m.reactions),
"\n",
"Number of metabolites: ",
length(m.metabolites),
"\n",
"Number of genes: ",
length(m.genes),
"Constraint based model: $(m.id)\nNumber of reactions: $(length(m.reactions))\nNumber of metabolites: $(length(m.metabolites))\nNumber of genes: $(length(m.genes))",
)
end
"""
change_constraint(reaction::Reaction, lb, ub)
change_constraint(id::String, lb, ub)
Change the lower and upper bounds (`lb` and `ub` respectively) of `reaction`.
Change the lower and upper bounds (`lb` and `ub` respectively) of reaction `id`.
"""
function change_constraint(reaction::Reaction, lb, ub)
function change_constraint(id::String, lb, ub)
(model, opt_model) -> begin
ind = model.reactions[reaction]
ind = first(indexin([id], reactions(model)))
set_bound(ind, opt_model, lb = lb, ub = ub)
end
end
"""
change_objective(objective_functions::Union{Reaction, Vector{Reaction}}; weights=[], sense=MOI.MAX_SENSE)
change_objective(objective_functions::Union{String,Vector{String}}; weights=[], sense=MOI.MAX_SENSE)
Callback function to change the objective function used in a constraint based analysis function.
`objective_functions` can be a single reaction or an array of reactions (of type `Reaction`).
`objective_functions` can be a single reaction or an array of reactions (input their string `id`s).
Optionally specify their `weights` and the sense of the new objective (`MOI.MAX_SENSE`, `MOI.MIN_SENSE`).
Note, this function sets the sense of the objective to `MOI.MAX_SENSE` by default if not specified.
"""
function change_objective(
objective_functions::Union{Reaction,Vector{Reaction}};
objective_functions::Union{String,Vector{String}};
weights = [],
sense = MOI.MAX_SENSE,
)
(model, opt_model) -> begin
# Construct objective_indices array
if typeof(objective_functions) == Reaction
objective_indices = [model[objective_functions]]
if typeof(objective_functions) == String
objective_indices = [first(indexin([objective_functions], reactions(model)))]
else
objective_indices = [model[rxn] for rxn in objective_functions]
objective_indices = [
first(indexin([rxnid], reactions(model))) for rxnid in objective_functions
]
end
# Initialize weights
opt_weights = zeros(length(model.reactions))
opt_weights = spzeros(length(model.reactions))
isempty(weights) && (weights = ones(length(objective_indices))) # equal weights
wcounter = 1
for i in eachindex(model.reactions)
if i in objective_indices
opt_weights[i] = weights[wcounter]
wcounter += 1
end
for (j, i) in enumerate(objective_indices)
opt_weights[i] = weights[j]
end
v = opt_model[:x]
......
......@@ -31,34 +31,34 @@ function Base.:+(
return push!(m1, m2)
end
function mkrxn(substrates, products)
metdict = Dict{Metabolite,Float64}()
function _mkrxn(substrates, products)
metdict = Dict{String,Float64}()
if typeof(substrates) == Metabolite
substrates != && (metdict[substrates] = get(metdict, substrates, 0.0) - 1.0)
substrates != && (metdict[substrates.id] = get(metdict, substrates.id, 0.0) - 1.0)
elseif typeof(substrates) == MetaboliteWithCoefficient
substrates.metabolite != && (
metdict[substrates.metabolite] =
get(metdict, substrates.metabolite, 0.0) - 1.0 * abs(substrates.coeff)
metdict[substrates.metabolite.id] =
get(metdict, substrates.metabolite.id, 0.0) - 1.0 * abs(substrates.coeff)
)
else
for mwc in substrates
metdict[mwc.metabolite] =
get(metdict, mwc.metabolite, 0.0) - 1.0 * abs(mwc.coeff)
metdict[mwc.metabolite.id] =
get(metdict, mwc.metabolite.id, 0.0) - 1.0 * abs(mwc.coeff)
end
end
if typeof(products) == Metabolite
products != && (metdict[products] = get(metdict, products, 0.0) + 1.0)
products != && (metdict[products.id] = get(metdict, products.id, 0.0) + 1.0)
elseif typeof(products) == MetaboliteWithCoefficient
products.metabolite != && (
metdict[products.metabolite] =
get(metdict, products.metabolite, 0.0) + abs(products.coeff)
metdict[products.metabolite.id] =
get(metdict, products.metabolite.id, 0.0) + abs(products.coeff)
)
else
for mwc in products
metdict[mwc.metabolite] =
get(metdict, mwc.metabolite, 0.0) + 1.0 * abs(mwc.coeff)
metdict[mwc.metabolite.id] =
get(metdict, mwc.metabolite.id, 0.0) + 1.0 * abs(mwc.coeff)
end
end
......@@ -76,7 +76,7 @@ function ⟶(
},
products::Union{Metabolite,MetaboliteWithCoefficient,Vector{MetaboliteWithCoefficient}},
)
metdict = mkrxn(substrates, products)
metdict = _mkrxn(substrates, products)
return Reaction("", metdict, :forward)
end
const =
......@@ -92,7 +92,7 @@ function ⟵(
},
products::Union{Metabolite,MetaboliteWithCoefficient,Vector{MetaboliteWithCoefficient}},
)
metdict = mkrxn(substrates, products)
metdict = _mkrxn(substrates, products)
return Reaction("", metdict, :reverse)
end
const =
......@@ -108,7 +108,7 @@ function ⟷(
},
products::Union{Metabolite,MetaboliteWithCoefficient,Vector{MetaboliteWithCoefficient}},
)
metdict = mkrxn(substrates, products)
metdict = _mkrxn(substrates, products)
return Reaction("", metdict, :bidirectional)
end
const =
......@@ -10,12 +10,7 @@ function add!(model::StandardModel, rxns::Vector{Reaction})
end
function add!(model::StandardModel, rxn::Reaction)
if model[rxn] == -1
push!(model.reactions, rxn)
else
@warn "$(rxn.id) already present in model."
end
return nothing
model.reactions[rxn.id] = rxn
end
"""
......@@ -30,12 +25,7 @@ function add!(model::StandardModel, mets::Vector{Metabolite})
end
function add!(model::StandardModel, met::Metabolite)
if model[met] == -1
push!(model.metabolites, met)
else
@warn "$(met.id) already present in model."
end
return nothing
model.metabolites[met.id] = met
end
"""
......@@ -50,12 +40,7 @@ function add!(model::StandardModel, genes::Vector{Gene})
end
function add!(model::StandardModel, gene::Gene)
if model[gene] == -1
push!(model.genes, gene)
else
@warn "$(gene.id) already present in model."
end
return nothing
model.genes[gene.id] = gene
end
"""
......@@ -115,122 +100,58 @@ macro add_reactions!(model::Symbol, ex::Expr)
end
"""
rm!(model::StandardModel, rxns::Union{Vector{Reaction}, Reaction})
rm!(::Type{Reaction}, model::StandardModel, ids::Vector{String})
Remove all `rxn(s)` from `model` if the `id`s match those in `rxns`.
Remove all reactions with `ids` from `model`.
# Example
rm!(Reaction, model, ["EX_glc__D_e", "fba"])
rm!(Reaction, model, "EX_glc__D_e")
"""
function rm!(model::StandardModel, rxns::Union{Vector{Reaction},Reaction})
new_rxn_list = Reaction[]
for r in model.reactions
if typeof(rxns) == Reaction
if rxns.id != r.id
push!(new_rxn_list, r)
end
else
if !(r.id in [rr.id for rr in rxns])
push!(new_rxn_list, r)
end
end
function rm!(::Type{Reaction}, model::StandardModel, ids::Vector{String})
for id in ids
rm!(Reaction, model, id)
end
model.reactions = new_rxn_list
return nothing
end
"""
rm!(model::StandardModel, mets::Union{Vector{Metabolite}, Metabolite})
Remove `met(s)` from `model` based on metabolite `id`.
"""
function rm!(model::StandardModel, mets::Union{Vector{Metabolite},Metabolite})
new_met_list = Metabolite[]
for m in model.metabolites
if typeof(mets) == Metabolite
if mets.id != m.id
push!(new_met_list, m)