Commit 1799b794 authored by St. Elmo's avatar St. Elmo
Browse files

updated tests

parent fe33d56b
......@@ -28,46 +28,52 @@ More importantly, it also exposes the user to the core structures used in COBRA,
## Installation
To install this package: `] add CobraTools`. See the documentation for more information.
To install this package: `] add https://github.com/stelmo/CobraTools.jl`. See the documentation for more information.
## Quick Example
Let's use `CobraTools` to perform classic flux balance analysis on an *E. coli* genome-scale metabolic model.
```julia
using CobraTools
using JuMP
using GLPK # pick any solver supported by JuMP
using Tulip # pick any solver supported by JuMP
# import E. coli model
model = read_model("iJO1366.json") # models have pretty printing
# Import E. coli model (models have pretty printing)
model = read_model("iJO1366.json")
# choose objective to maximize
# Choose objective to maximize (biomass is a reaction struct, which also has pretty printing)
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
# FBA - use convenience functions
sol = fba(model, biomass, GLPK.Optimizer; solver_attributes=Dict("msg_lev" => GLPK.GLP_MSG_OFF)) # classic flux balance analysis
sol = fba(model, biomass, Tulip.Optimizer))
```
# FBA DIY - automatically construct basic JuMP model from a constraint based model
cbm, v, mb, ubs, lbs = build_cbm(model) # get the constraint based model (cbm) in JuMP format: S*v=b (mb: mass balance constraints) with lbs <= v <= ubs.
set_optimizer(cbm, GLPK.Optimizer) # use JuMP functions to set optimizer
set_optimizer_attribute(cbm, "msg_lev", GLPK.GLP_MSG_OFF) # use JuMP functions to set optimizer attributes
@objective(cbm, Max, v[model[biomass]]) # use index notation to get biomass equation index
optimize!(cbm)
sol = map_fluxes(v, model) # map fluxes to reaction ids.
# FBA really DIY - manually construct a JuMP model from constraint based model
S, b, ubvec, lbvec = get_core_model(model) # get S*v = b with lbvec <= v <= ubvec from model
cbm = JuMP.Model()
nvars = size(S, 2) # number of variables in model
v = @variable(cbmodel, v[1:nvars]) # flux variables
mb = @constraint(cbmodel, mb, S*v .== b) # mass balance
lbs = @constraint(cbmodel, lbs, lbs .<= v) # lower bounds
ubs = @constraint(cbmodel, ubs, v .<= ubs) # upper bounds
set_optimizer(cbm, GLPK.Optimizer) # use JuMP functions to set optimizer
set_optimizer_attribute(cbm, "msg_lev", GLPK.GLP_MSG_OFF) # use JuMP functions to set optimizer attributes
@objective(cbm, Max, v[model[biomass]]) # use index notation to get biomass equation index
If you are feeling more adventurous you can perform the optimization yourself using `JuMP`.
```julia
# Get the constraint based model (cbm) in JuMP format: S*v=b (mb: mass balance constraints) with lbs <= v <= ubs
cbm, v, mb, ubs, lbs = build_cbm(model)
# Use JuMP functions to optimize the constraint based model
set_optimizer(cbm, Tulip.Optimizer)
# Use index notation to get biomass equation index
@objective(cbm, Max, v[model[biomass]])
optimize!(cbm)
sol = map_fluxes(v, model) # map fluxes to reaction ids.
# Map fluxes to reaction ids
sol = map_fluxes(v, model)
```
If you are feeling even more adventurous you can do everything yourself!
```julia
# Get S*v = b with lbvec <= v <= ubvec from model
S, b, ubvec, lbvec = get_core_model(model)
# Manually define the model in JuMP
cbm = JuMP.Model(Tulip.Optimizer)
nvars = size(S, 2)
v = @variable(cbmodel, v[1:nvars])
mb = @constraint(cbmodel, mb, S*v .== b)
lbs = @constraint(cbmodel, lbs, lbs .<= v)
ubs = @constraint(cbmodel, ubs, v .<= ubs)
@objective(cbm, Max, v[model[biomass]])
optimize!(cbm)
sol = map_fluxes(v, model)
```
More funcionality is described in the documention.
......
......@@ -17,7 +17,7 @@
## Installation
To install this package: `] add CobraTools`.
To install this package: `] add https://github.com/stelmo/CobraTools.jl`.
Some of the optional features used in this package require external programs and/or data to be available. These are described below:
......
......@@ -5,16 +5,15 @@
g.notes = Dict("notes"=>["blah", "blah"])
g.annotation = Dict("sboterm" => "sbo", "ncbigene" => ["ads", "asds"])
# @test repr("text/plain", g) == "Gene ID: gene1\nGene name: gene_name\n"
@test sprint(show, MIME("text/plain"), g) == "Gene ID: gene1\nGene name: gene_name\n"
g2 = Gene("gene2")
genes = [g, g2]
@test repr("text/plain", genes) == "Gene set of length: 2\n"
@test sprint(show, MIME("text/plain"), genes) == "Gene set of length: 2\n"
gene_list = [[g], [g2]]
@test repr("text/plain", gene_list) == "(gene1) or (gene2)\n"
@test sprint(show, MIME("text/plain"), gene_list) == "(gene1) or (gene2)\n"
@test genes[g] == 1
......@@ -26,4 +25,11 @@
dup, ind = check_duplicate_annotations(genes, g3)
@test dup && ind == 1
@test isnothing(findfirst(genes, "nope"))
g4 = Gene("g4")
g4.annotation = Dict("ncbigene" => "sbo", "ncbigene" => ["ads22", "asd22s"])
dup, ind = check_duplicate_annotations(genes, g4)
@test !dup && ind == -1
end
......@@ -8,7 +8,7 @@
m1.notes = Dict("notes"=>["blah", "blah"])
m1.annotation = Dict("sboterm" => "sbo", "kegg.compound" => ["ads", "asds"])
@test repr("text/plain", m1) == "Metabolite ID: met1\nMetabolite name: metabolite 1\nFormula: C6H12O6N\nCharge: 1\n"
@test sprint(show, MIME("text/plain"), m1) == "Metabolite ID: met1\nMetabolite name: metabolite 1\nFormula: C6H12O6N\nCharge: 1\n"
m2 = Metabolite("met2")
m2.formula = "C6H12O6N"
......@@ -19,7 +19,7 @@
mets = [m1, m2, m3]
@test repr("text/plain", mets) == "Metabolite set of length: 3\n"
@test sprint(show, MIME("text/plain"), mets) == "Metabolite set of length: 3\n"
@test mets[m2] == 2
......@@ -34,4 +34,13 @@
ats = get_atoms(mms[1])
@test ats["C"] == 6 && ats["N"] == 1
@test isnothing(findfirst(mets, "nope"))
m4 = Metabolite("met4")
m4.formula = "X"
m4.annotation = Dict("sboterm" => "sbo", "kegg.compound" => ["adxxx2s", "asdxxxs"])
dup, ind = check_duplicate_annotations(mets, m4)
@test !dup && ind == -1
end
\ No newline at end of file
......@@ -37,7 +37,7 @@
model.metabolites = mets
model.genes = genes
@test repr("text/plain", model) == "Constraint based model: model\nNumber of reactions: 4\nNumber of metabolites: 4\nNumber of genes: 3\n"
@test sprint(show, MIME("text/plain"), model) == "Constraint based model: model\nNumber of reactions: 4\nNumber of metabolites: 4\nNumber of genes: 3\n"
@test model[r2] == 2
......
......@@ -5,6 +5,14 @@
m2.formula = "H3C2"
m3 = Metabolite("m3")
m4 = Metabolite("m4")
m5 = Metabolite("m5")
m6 = Metabolite("m6")
m7 = Metabolite("m7")
m8 = Metabolite("m8")
m9 = Metabolite("m9")
m10 = Metabolite("m10")
m11 = Metabolite("m11")
m12 = Metabolite("m12")
g1 = Gene("g1")
g2 = Gene("g2")
......@@ -22,8 +30,16 @@
r1.annotation = Dict("sboterm" => "sbo", "biocyc" => ["ads", "asds"])
r1.objective_coefficient = 1.0
@test repr("text/plain", r1) == "Reaction ID: r1\nReaction name: reaction 1\nReaction subsystem: glycolysis\n1.0 m1 ⟷ 1.0 m2\nLower bound: -100.0\nUpper bound: 100.0\nGenes: (g1 and g2) or (g3)\nE.C. number: \n"
@test sprint(show, MIME("text/plain"), r1) == "Reaction ID: r1\nReaction name: reaction 1\nReaction subsystem: glycolysis\n1.0 m1 ⟷ 1.0 m2\nLower bound: -100.0\nUpper bound: 100.0\nGenes: (g1 and g2) or (g3)\nE.C. number: \n"
rlongfor = Reaction("rlongfor", Dict(m1 => -1.0, m2 => -1.0, m3 => -1.0, m4 => -1.0, m5 => -1.0, m6 => -1.0,
m7 => 1.0, m8 => 1.0, m9 => 1.0, m10 => 1.0, m11 => 1.0, m12 => 1.0,), "for")
@test occursin("...", sprint(show, MIME("text/plain"), rlongfor)) # from dictionaries so order is not consistent.
rlongrev = Reaction("rlongrev", Dict(m1 => -1.0, m2 => -1.0, m3 => -1.0, m4 => -1.0, m5 => -1.0, m6 => -1.0,
m7 => 1.0, m8 => 1.0, m9 => 1.0, m10 => 1.0, m11 => 1.0, m12 => 1.0,), "rev")
@test occursin("...", sprint(show, MIME("text/plain"), rlongrev))
r2 = Reaction("r2", Dict(m1 => -2.0, m4 => 1.0), "rev")
@test r2.lb == -1000.0 && r2.ub == 0.0
......@@ -32,7 +48,7 @@
rxns = [r1, r2, r3]
@test repr("text/plain", rxns) == "Reaction set of length: 3\n"
@test sprint(show, MIME("text/plain"), rxns) == "Reaction set of length: 3\n"
@test rxns[r3] == 3
......@@ -58,4 +74,8 @@
bal, d = is_mass_balanced(r1)
@test bal
@test isnothing(findfirst(rxns, "nope"))
end
......@@ -30,4 +30,10 @@
rxn = 1.0*nadh + 4.0*h_c + 1.0*q8 1.0*q8h2 + 1.0*nad + 3.0*h_p
@test prod(values(rxn.metabolites)) == -12
@test ("q8h2_c" in [x.id for x in keys(rxn.metabolites)])
rxn = nadh + 4.0*h_c + 1.0*q8 1.0*q8h2 + 1.0*nad + 3.0*h_p
@test rxn.lb == 0.0 && rxn.ub > 0.0
@test length(h_p + h_p) == 2
@test length(h_p + h_p + h_p) == 3
end
\ No newline at end of file
......@@ -6,6 +6,10 @@
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
optimizer = Tulip.Optimizer # quiet by default
sol = fba(model, biomass, optimizer)
flux_vec = [sol[rxn.id] for rxn in model.reactions]
sol_mapped = map_fluxes(flux_vec, model)
@test isapprox(sol_mapped["BIOMASS_Ecoli_core_w_GAM"], 0.8739215022678488, atol=1e-6)
@test isapprox(sol["BIOMASS_Ecoli_core_w_GAM"], 0.8739215022678488, atol=1e-6)
# exchange trackers
......
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