Commit e5eb9304 authored by St. Elmo's avatar St. Elmo
Browse files

update tests

parent 29821912
......@@ -22,8 +22,8 @@
| [![docs-img]][docs-url] | [![CI][ci-img]][ci-url] | [![codecov][cov-img]][cov-url] |
This is package aims to provide constraint based reconstruction and analysis (COBRA) tools in the Julia environment, similar to Cobrapy in Python and the Cobra Toolbox in Matlab.
This package provides basic convenience functions, e.g. FBA, pFBA, sampling, model construction, etc.
More importantly, it also exposes the user to the core structures used in COBRA, e.g. the stoichiometric matrix, etc., so that custom optimization routines can be written as painlessly as possible (due in large part to JuMP). An alternative, [COBRA.jl](https://github.com/opencobra/COBRA.jl), but its scope is more restricted than `CobraTools.jl`.
This package provides basic convenience functions, e.g. model IO, construction, modification, FBA, pFBA, sampling, etc. It can also be used to interface with [Equilibrator](http://equilibrator.weizmann.ac.il/) and [BRENDA](https://www.brenda-enzymes.org/).
More importantly, it also exposes the user to the core structures used in COBRA, e.g. the stoichiometric matrix, etc., so that custom optimization routines can be written as painlessly as possible (due in large part to JuMP). An alternative Julia package hosted by OpenCobra, [COBRA.jl](https://github.com/opencobra/COBRA.jl), also exists, but its scope is more restricted than `CobraTools.jl`.
## Installation
......@@ -46,8 +46,23 @@ biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
# Use convenience functions
sol = fba(model, biomass, GLPK.Optimizer; solver_attributes=Dict("msg_lev" => GLPK.GLP_MSG_OFF)) # classic flux balance analysis
# DIY
# 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.
# Really DIY - manually construct a JuMP model from constraint based model
S, b, ubvec, lbvec = get_core_model(model) # 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
......
......@@ -66,13 +66,13 @@ function test_gene()
g.notes = Dict("notes"=>["blah", "blah"])
g.annotation = Dict("sboterm" => "sbo", "ncbigene" => ["ads", "asds"])
g # test IO
repr(show(g)) # test IO
g2 = Gene("gene2")
genes = [g, g2]
genes # test IO
repr(show(genes)) # test IO
ind = genes[g]
if ind != 1
......@@ -108,7 +108,7 @@ function test_metabolite()
m1.notes = Dict("notes"=>["blah", "blah"])
m1.annotation = Dict("sboterm" => "sbo", "kegg.compound" => ["ads", "asds"])
m1 # test IO
repr(show(m1)) # test IO
m2 = Metabolite("met2")
m2.formula = "C6H12O6N"
......@@ -119,7 +119,7 @@ function test_metabolite()
mets = [m1, m2, m3]
mets # test IO
repr(show(mets)) # test IO
ind = mets[m2]
if ind != 2
......@@ -176,7 +176,7 @@ function test_reaction()
r1.annotation = Dict("sboterm" => "sbo", "biocyc" => ["ads", "asds"])
r1.objective_coefficient = 1.0
r1 # test IO
repr(show(r1)) # test IO
r2 = Reaction("r2", Dict(m1 => -2.0, m4 => 1.0), "rev")
if r2.lb != -1000.0 && r2.ub != 0.0
......@@ -190,7 +190,7 @@ function test_reaction()
rxns = [r1, r2, r3]
rxns # test IO
repr(show(rxns)) # test IO
ind = rxns[r3]
if ind != 3
......@@ -268,7 +268,7 @@ function test_model()
model.metabolites = mets
model.genes = genes
model # test IO
repr(show(model)) # test IO
ind = model[r2]
if ind != 2
......
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