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

moved scripts around

parent 8e041b1c
......@@ -5,7 +5,7 @@
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 machinery used in this type of analysis, e.g. the stoichiometric matrix, so that custom optimization routines can be written as painlessly as possible (due in large part to JuMP).
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).
## Installation
......@@ -43,4 +43,5 @@ TODO.
### Citations
1) Ebrahim, A., Lerman, J.A., Palsson, B.O. & Hyduke, D. R. (2013). COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Systems Biology, 7(74). https://doi.org/10.1186/1752-0509-7-74
2) Heirendt, L., Arreckx, S., Pfau, T. et al. (2019). Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat Protoc 14, 639–702. https://doi.org/10.1038/s41596-018-0098-2
3) Noor, E., Bar-Even, A., Flamholz, A., Lubling, Y., Davidi, D., & Milo, R. (2012). An integrated open framework for thermodynamics of reactions that combines accuracy and coverage. Bioinformatics, 28(15), 2037–2044. https://doi.org/10.1093/bioinformatics/bts317
\ No newline at end of file
3) Noor, E., Bar-Even, A., Flamholz, A., Lubling, Y., Davidi, D., & Milo, R. (2012). An integrated open framework for thermodynamics of reactions that combines accuracy and coverage. Bioinformatics, 28(15), 2037–2044. https://doi.org/10.1093/bioinformatics/bts317
4) Cite Brenda
\ No newline at end of file
......@@ -19,9 +19,8 @@ include("model_tools.jl")
include("io_tools.jl")
include("construction_tools.jl")
# include("basic_analysis.jl")
# include("equilibrator_tools.jl")
include("brenda_tools.jl")
# include("basic_analysis.jl")
# include("sampling_tools.jl")
# export
......
using CobraTools
# using SBML
modelpath = joinpath("models", "iMM904.xml") # doesn' work
# modelpath = joinpath("models", "iJO1366.xml") # doesn't work
# modelpath = joinpath("models", "e_coli_core.xml") # doesn't work
# modelpath = joinpath("models", "Ec_core_flux1.xml") # works
# model = readSBML(modelpath)
# modelpath = joinpath("models", "iAF1260.xml")
# modelpath = joinpath("models", "yeastGEM.xml")
model = CobraTools.read_model(modelpath)
using CobraTools
using Measurements
using LinearAlgebra
using Statistics
using JLD
using JuMP
using Gurobi
using Plots
pyplot()
# E. coli model
modelpath = joinpath("models", "iML1515.json")
model = CobraTools.read_model(modelpath)
###### Set constraints
cbmodel, v, mb, ubs, lbs = CobraTools.CBM(model)
set_optimizer(cbmodel, Gurobi.Optimizer)
set_optimizer_attribute(cbmodel, "OutputFlag", 0) # quiet
biomass_index = model[findfirst(model.rxns, "BIOMASS_Ec_iML1515_core_75p37M")]
glucose_index = model[findfirst(model.rxns, "EX_glc__D_e")]
o2_index = model[findfirst(model.rxns, "EX_o2_e")]
atpm_index = model[findfirst(model.rxns, "ATPM")]
etoh_index = model[findfirst(model.rxns, "EX_etoh_e")]
CobraTools.set_bound(glucose_index, ubs, lbs; ub=-1.0, lb=-1.0)
CobraTools.set_bound(o2_index, ubs, lbs; ub=1000.0, lb=0.0)
CobraTools.set_bound(atpm_index, ubs, lbs; ub=1000.0, lb=0.0)
# Do FBA
@objective(cbmodel, Max, v[biomass_index])
optimize!(cbmodel)
termination_status(cbmodel) != MOI.OPTIMAL && @warn "Optimization issue..."
μ_max = round(objective_value(cbmodel), digits=6)
CobraTools.set_bound(biomass_index, ubs, lbs; ub=μ_max, lb=0.9*μ_max)
##################################
# Get warmup points
wpoints = CobraTools.get_warmup_points(cbmodel, v, ubs, lbs, numstop=1000) # very slow
# sample!
samples = @time CobraTools.hit_and_run(100_000, wpoints, ubs, lbs; keepevery=10, samplesize=5000, W=1000)
###########################
violation_inds = CobraTools.test_samples(samples, model, ubs, lbs)
plot(samples[etoh_index, :])
using CobraTools
using JuMP
using Gurobi # use your favourite solver
using Measurements
using LinearAlgebra
using MCMCChains
using HypothesisTests
using StatsBase
using Turing
using Distributions
using DataFrames
using StatsPlots
pyplot()
# E. coli model
modelpath = joinpath("models", "e_coli_core.json")
model = CobraTools.read_model(modelpath)
#### set bounds
cbmodel, v, mb, ubs, lbs = CobraTools.CBM(model)
set_optimizer(cbmodel, Gurobi.Optimizer)
set_optimizer_attribute(cbmodel, "OutputFlag", 0) # quiet
biomass_index = model[findfirst(model.rxns, "BIOMASS_Ecoli_core_w_GAM")]
glucose_index = model[findfirst(model.rxns, "EX_glc__D_e")]
o2_index = model[findfirst(model.rxns, "EX_o2_e")]
atpm_index = model[findfirst(model.rxns, "ATPM")]
etoh_index = model[findfirst(model.rxns, "EX_etoh_e")]
# Fix glucose use 1.0 then normalization is easy. NB - if not 1 then change normalization!!
CobraTools.set_bound(glucose_index, ubs, lbs; ub=-1.0, lb=-1.01)
# Aerobic
# CobraTools.set_bound(o2_index, ubs, lbs; ub=1000.0, lb=-1000.0)
# Anaerobic
CobraTools.set_bound(o2_index, ubs, lbs; ub=1000.0, lb=0.0)
# No free ATP generation
CobraTools.set_bound(atpm_index, ubs, lbs; ub=1000.0, lb=0.0)
@objective(cbmodel, Max, v[biomass_index])
optimize!(cbmodel)
termination_status(cbmodel) != MOI.OPTIMAL && @warn "Optimization issue..."
μ = objective_value(cbmodel)
### Fix biomass as a constraint
CobraTools.set_bound(biomass_index, ubs, lbs; ub=μ, lb=0.99*μ)
ubvec, lbvec = CobraTools.get_bound_vectors(ubs, lbs)
S, _, _, _ = CobraTools.get_core_model(model)
###########################################
#### Do inference
S = [1 -1.0 0 0 0 0 0;0 1 -1 -1 -1 0 0;0 0 1 1 1 -1 0;0 0 0 0 0 1 -1]
Sin = S[:, 2:end]
@model function turingtest(S)
Σ ~ InverseWishart(6 + 1, Matrix{Float64}(I, 6, 6))
flux ~ MvNormal(zeros(6), Symmetric(Σ))
θ ~ Uniform(0.0, 100.0)
b = sum(abs, S*flux - [1.0, 0.0, 0.0, 0.0]) # input is r₀ = -1
0.0 ~ truncated(Normal(b, θ), 0.0, Inf)
end
tmodel = turingtest(Sin)
chain = sample(tmodel, NUTS(0.65), 1_000)
################
df = DataFrame(x1 = rand(Normal(), 100), x2 = rand(Normal(), 100), x3 = rand(Normal(), 100))
beta = [5, 4, 5, 6]
df[!, :y] = beta[1] .+ beta[2]*df[!, :x1] + beta[3]*df[!, :x2] + beta[4]*df[!, :x3] + rand(Normal(), 100)
@model BayesReg(y, x1, x2, x3) = begin
n_obs = length(y)
σ₂ ~ InverseGamma(0.001, 0.001)
Sigma ~ InverseWishart(
4 + 1,
Matrix{Float64}(I, 4, 4)
)
Mu ~ MvNormal(
zeros(4),
Matrix{Float64}(I, 4, 4)
)
beta_hat ~ MvNormal(Mu, Symmetric(Sigma))
mu = beta_hat[1] .+ beta_hat[2]*x1 .+ beta_hat[3]*x2 .+ beta_hat[4]*x3
for i = 1:n_obs
y[i] ~ Normal(mu[i], σ₂)
end
end
model = BayesReg(df[!, :y], df[!, :x1], df[!, :x2], df[!, :x3])
chain = sample(model, NUTS(0.65), 1000)
Supports Markdown
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