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

begin doc

parent c60830fd
......@@ -2,4 +2,5 @@
*.code-workspace
*.jld
*.png
*.brenda.txt
\ No newline at end of file
*.brenda.txt
/src/depr_functions.jl
\ No newline at end of file
......@@ -8,6 +8,7 @@ Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
......
......@@ -2,7 +2,7 @@
This is package aims to provide constraint based reconstruction and analysis (COBRA) tools in the Julia environment.
This package provides basic convenience functions, e.g. FBA, pFBA, sampling, model construction, etc.
More importantly, it also exposes the user to the core machinary 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 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).
## Installation
......@@ -26,12 +26,14 @@ To install this package: `] add CobraTools`.
- [X] pFBA
- [ ] MOMA
- [ ] FVA
- [x] Implement sampling
- [x] Implement sampling (hit and run)
- [x] Implement sampling (achr - kind of?)
- [ ] Single gene knockouts
- [ ] Double gene knockout
- [ ] Model construction tools
- [x] Gibbs integration
- [x] Brenda integration (basic)
- [x] Reaction construction
- [ ] Model modifications
- [ ] Distributed analysis (COBRA.jl integration?)
## Usage
......
......@@ -28,7 +28,7 @@ include("construction_tools.jl")
include("basic_analysis.jl")
include("gibbs_tools.jl")
include("brenda_tools.jl")
include("sampling.jl")
include("sampling_tools.jl")
# export
= Metabolite("∅") # for exchange reactions
......
function ACHRSampler(model, warmupPoints, fileName, nFiles, pointsPerFile, stepsPerPoint, initPoint, fileBaseNo, maxTime, printLevel, minFlux, maxFlux)
% Artificial Centering Hit-and-Run sampler
%
% USAGE:
%
% ACHRSampler(model, warmupPoints, fileName, nFiles, pointsPerFile, stepsPerPoint, initPoint, fileBaseNo, maxTime)
%
% INPUTS:
% model: Model structure
% warmupPoints: Warmup points
% fileName: Base `fileName` for saving results
% nFiles: Number of files created
% pointsPerFile: Number of points per file saved
% stepsPerPoint: Number of sampler steps per point saved
% printLevel: 0 silent / 1 active
%
% OPTIONAL INPUTS:
% initPoint: Initial point (Default = centerPoint)
% fileBaseNo: Base file number for continuing previous sampler run
% (Default = 0)
% maxTime: Maximum time limit (Default = 36000 s)
% minFlux: Lower bound of the solution space given by FVA
% maxFlux: Upper bound of the solution space given by FVA
%
% .. Authors: -
% - Markus Herrgard 4/14/06
% - Gregory Hannum 4/14/06
% - Ines Thiele 4/14/06
% - Nathan Price 4/14/06
warning off MATLAB:divideByZero;
if (nargin < 8) || isempty(fileBaseNo)
fileBaseNo = 0;
end
if (nargin < 9) || isempty(maxTime)
maxTime = 10*3600;
end
if nargin<10
printLevel=1;
end
if nargin<12
minFlux = model.lb;
maxFlux = model.ub;
end
N = null(full(model.S));
% Minimum allowed distance to the closest constraint
maxMinTol = 1e-9;
% Ignore directions where u is really small
uTol = 1e-9;
% Project out of directions that are too close to the boundary
dTol = 1e-14;
% Number of warmup points
[nRxns,nWrmup] = size(warmupPoints);
% Find the center of the space
centerPoint = mean(warmupPoints,2);
% Set the start point
if (nargin < 7) || isempty(initPoint)
prevPoint = centerPoint;
else
prevPoint = initPoint;
end
fidErr = fopen('ACHRerror.txt','w');
totalStepCount = 0;
showprogress(0,'ACHR sampling in progress ...');
totalCount = nFiles*pointsPerFile*stepsPerPoint;
t0 = cputime;
fprintf('File #\tPoint #\tStep #\tTime\t#Time left\n');
for i = 1:nFiles
% Allocate memory for all points
points = zeros(nRxns,pointsPerFile);
pointCount = 1;
while (pointCount <= pointsPerFile)
% Create the random step size vector
randVector = rand(stepsPerPoint,1);
stepCount = 1;
while (stepCount <= stepsPerPoint)
% Pick a random warmup point
randPointID = ceil(nWrmup*rand);
randPoint = warmupPoints(:,randPointID);
% Get a direction from the center point to the warmup point
u = (randPoint-centerPoint);
u = u/norm(u);
% Figure out the distances to upper and lower bounds
distUb = (maxFlux - prevPoint);
distLb = (prevPoint - minFlux);
% Figure out if we are too close to a boundary
validDir = ((distUb > dTol) & (distLb > dTol));
%model.rxns(~validDir)
% Zero out the directions that would bring us too close to the
% boundary. This may cause problems.
%u(~validDir) = 0;
% Figure out positive and negative directions
posDirn = find(u(validDir) > uTol);
negDirn = find(u(validDir) < -uTol);
% Figure out all the possible maximum and minimum step sizes
maxStepTemp = distUb(validDir)./u(validDir);
minStepTemp = -distLb(validDir)./u(validDir);
maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn)];
minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn)];
% Figure out the true max & min step sizes
maxStep = min(maxStepVec);
minStep = max(minStepVec);
%fprintf('%f\t%f\n',minStep,maxStep);
% Find new direction if we're getting too close to a constraint
if (abs(minStep) < maxMinTol & abs(maxStep) < maxMinTol) | (minStep > maxStep)
fprintf('Warning %f %f\n',minStep,maxStep);
continue;
end
% Pick a rand out of list_of_rands and use it to get a random
% step distance
stepDist = randVector(stepCount)*(maxStep-minStep)+minStep;
% Advance to the next point
curPoint = prevPoint + stepDist*u;
% Reproject the current point and go to the next step
if mod(totalStepCount,10) == 0
if (full(max(max(abs(model.S*curPoint)))) > 1e-9)
curPoint = N*(N'*curPoint);
end
end
% Print out errors
if (mod(totalStepCount,2000)==0)
fprintf(fidErr,'%10.8f\t%10.8f\t',max(curPoint-maxFlux),max(minFlux-curPoint));
end
timeElapsed = cputime-t0;
% Print step information
if (mod(totalStepCount,5000)==0)
timePerStep = timeElapsed/totalStepCount;
if printLevel
fprintf('%d\t%d\t%d\t%8.2f\t%8.2f\n',i,pointCount,totalStepCount,timeElapsed/60,(totalCount-totalStepCount)*timePerStep/60);
end
end
overInd = find(curPoint > maxFlux);
underInd = find(curPoint < minFlux);
if (any((maxFlux-curPoint) < 0) | any((curPoint-minFlux) ...
< 0))
curPoint(overInd) = maxFlux(overInd);
curPoint(underInd) = minFlux(underInd);
end
if (mod(totalStepCount,2000)==0)
if printLevel
fprintf(fidErr,'%10.8f\n',full(max(max(abs(model.S*curPoint)))));
end
end
prevPoint = curPoint;
stepCount = stepCount + 1;
% Count the total number of steps
totalStepCount = totalStepCount + 1;
if printLevel
showprogress(totalStepCount/totalCount);
end
%recalculate the center point
centerPoint = ((nWrmup+totalStepCount)*centerPoint + curPoint)/(nWrmup+totalStepCount+1);
% Exit if time exceeded
if (timeElapsed > maxTime)
points(:,pointCount) = curPoint;
file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
save (file,'points');
save ACHR_last_point.mat curPoint
return;
end
end % Steps per point
% Add the current point to points
points(:,pointCount) = curPoint;
pointCount = pointCount + 1;
end % Points per cycle
% Save current points to a file
file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
save (file,'points');
end
......@@ -7,7 +7,7 @@ This is useful if you want to write your own optimization problem.
cbmodel is the JuMP model. v are the fluxes, mb is S*v == 0, and lbs <= v <= ubs.
"""
function CBM(model::Model)
function build_cbm(model::Model)
S, b, ubs, lbs = get_core_model(model) # Construct S, b, lbs, ubs from model
cbmodel = JuMP.Model()
nvars = size(S, 2) # number of variables in model
......@@ -25,7 +25,7 @@ Run flux balance analysis (FBA) on the model using objective_rxn(s) and optional
Optimiser can also be set here. Uses the constraints implied by the model object.
"""
function fba(model::Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}; weights=Float64[], optimizer="gurobi")
cbm, _, _, _, _ = CBM(model) # get the base constraint based model
cbm, _, _, _, _ = build_cbm(model) # get the base constraint based model
# ensure that an array of objective indices are fed in
if typeof(objective_rxns) == Reaction
......@@ -88,7 +88,7 @@ Optimiser can also be set here. Uses the constraints implied by the model object
"""
function pfba(model::Model, objective_rxns::Union{Reaction, Array{Reaction, 1}}; weights=Float64[], optimizer="gurobi")
## FBA ################################################
cbm, _, _, _, _ = CBM(model) # get the base constraint based model
cbm, _, _, _, _ = build_cbm(model) # get the base constraint based model
# ensure that an array of objective indices are fed in
if typeof(objective_rxns) == Reaction
......@@ -169,27 +169,7 @@ function pfba(model::Model, objective_rxns::Union{Reaction, Array{Reaction, 1}};
end
"""
atom_balance_dict = atom_exchange(fluxdict::Dict{String, Float64}, model::Model)
Return the composition of atoms consumed or produced by the model according to fluxdict.
"""
function atom_exchange(fluxdict::Dict{String, Float64}, model::Model)
exrxns = [k for k in keys(fluxdict) if startswith(k, "EX_")] # get exchange reactions
atom_balance = Dict{String, Float64}()
for exrxn in exrxns
rxn = findfirst(model.rxns, exrxn)
for (met, w) in rxn.metabolites
for (atom, stoich) in get_atoms(met)
atom_balance[atom] = get(atom_balance, atom, 0.0) + stoich*w*fluxdict[exrxn]
end
end
end
return atom_balance
end
"""
map_fluxes(v, model::Model)
rxn_flux_dict = map_fluxes(v, model::Model)
Map fluxes from an optimization problem (v) to rxns in a model. v can be a JuMP object (fluxes) or an array of Float64 fluxes.
Assumes they are in order, which they should be since they are constructed from model.
......@@ -231,7 +211,7 @@ get_exchanges(rxndict::Dict{String, Float64}; topN=8, ignorebound=1000)
Display the topN producing and consuming exchange fluxes. Ignores infinite (problem upper/lower bound) fluxes (set with ignorebound).
"""
function get_exchanges(rxndict::Dict{String, Float64}; topN=8, ignorebound=1000)
function display_exchange_reactions(rxndict::Dict{String, Float64}; topN=8, ignorebound=1000)
fluxes = Float64[]
rxns = String[]
for (k, v) in rxndict
......@@ -245,10 +225,22 @@ function get_exchanges(rxndict::Dict{String, Float64}; topN=8, ignorebound=1000)
println("Consuming fluxes:")
for i in 1:topN
if rxndict[rxns[inds_cons[i]]] > 0
continue
end
println(rxns[inds_cons[i]], " = ", round(rxndict[rxns[inds_cons[i]]], digits=4))
end
println("Producing fluxes:")
for i in 1:topN
if rxndict[rxns[inds_cons[i]]] < 0
continue
end
println(rxns[inds_prod[i]], " = ", round(rxndict[rxns[inds_prod[i]]], digits=4))
end
end
"""
"""
function display_metabolite_fluxes(fluxdict::Dict{String, Float64}, model::Model)
end
......@@ -103,3 +103,219 @@
# end
# return total_ΔG, missing_flux/(missing_flux+found_flux) # units J/gDW/h
# end
function ACHRSampler(model, warmupPoints, fileName, nFiles, pointsPerFile, stepsPerPoint, initPoint, fileBaseNo, maxTime, printLevel, minFlux, maxFlux)
% Artificial Centering Hit-and-Run sampler
%
% USAGE:
%
% ACHRSampler(model, warmupPoints, fileName, nFiles, pointsPerFile, stepsPerPoint, initPoint, fileBaseNo, maxTime)
%
% INPUTS:
% model: Model structure
% warmupPoints: Warmup points
% fileName: Base `fileName` for saving results
% nFiles: Number of files created
% pointsPerFile: Number of points per file saved
% stepsPerPoint: Number of sampler steps per point saved
% printLevel: 0 silent / 1 active
%
% OPTIONAL INPUTS:
% initPoint: Initial point (Default = centerPoint)
% fileBaseNo: Base file number for continuing previous sampler run
% (Default = 0)
% maxTime: Maximum time limit (Default = 36000 s)
% minFlux: Lower bound of the solution space given by FVA
% maxFlux: Upper bound of the solution space given by FVA
%
% .. Authors: -
% - Markus Herrgard 4/14/06
% - Gregory Hannum 4/14/06
% - Ines Thiele 4/14/06
% - Nathan Price 4/14/06
warning off MATLAB:divideByZero;
if (nargin < 8) || isempty(fileBaseNo)
fileBaseNo = 0;
end
if (nargin < 9) || isempty(maxTime)
maxTime = 10*3600;
end
if nargin<10
printLevel=1;
end
if nargin<12
minFlux = model.lb;
maxFlux = model.ub;
end
N = null(full(model.S));
% Minimum allowed distance to the closest constraint
maxMinTol = 1e-9;
% Ignore directions where u is really small
uTol = 1e-9;
% Project out of directions that are too close to the boundary
dTol = 1e-14;
% Number of warmup points
[nRxns,nWrmup] = size(warmupPoints);
% Find the center of the space
centerPoint = mean(warmupPoints,2);
% Set the start point
if (nargin < 7) || isempty(initPoint)
prevPoint = centerPoint;
else
prevPoint = initPoint;
end
fidErr = fopen('ACHRerror.txt','w');
totalStepCount = 0;
showprogress(0,'ACHR sampling in progress ...');
totalCount = nFiles*pointsPerFile*stepsPerPoint;
t0 = cputime;
fprintf('File #\tPoint #\tStep #\tTime\t#Time left\n');
for i = 1:nFiles
% Allocate memory for all points
points = zeros(nRxns,pointsPerFile);
pointCount = 1;
while (pointCount <= pointsPerFile)
% Create the random step size vector
randVector = rand(stepsPerPoint,1);
stepCount = 1;
while (stepCount <= stepsPerPoint)
% Pick a random warmup point
randPointID = ceil(nWrmup*rand);
randPoint = warmupPoints(:,randPointID);
% Get a direction from the center point to the warmup point
u = (randPoint-centerPoint);
u = u/norm(u);
% Figure out the distances to upper and lower bounds
distUb = (maxFlux - prevPoint);
distLb = (prevPoint - minFlux);
% Figure out if we are too close to a boundary
validDir = ((distUb > dTol) & (distLb > dTol));
%model.rxns(~validDir)
% Zero out the directions that would bring us too close to the
% boundary. This may cause problems.
%u(~validDir) = 0;
% Figure out positive and negative directions
posDirn = find(u(validDir) > uTol);
negDirn = find(u(validDir) < -uTol);
% Figure out all the possible maximum and minimum step sizes
maxStepTemp = distUb(validDir)./u(validDir);
minStepTemp = -distLb(validDir)./u(validDir);
maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn)];
minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn)];
% Figure out the true max & min step sizes
maxStep = min(maxStepVec);
minStep = max(minStepVec);
%fprintf('%f\t%f\n',minStep,maxStep);
% Find new direction if we're getting too close to a constraint
if (abs(minStep) < maxMinTol & abs(maxStep) < maxMinTol) | (minStep > maxStep)
fprintf('Warning %f %f\n',minStep,maxStep);
continue;
end
% Pick a rand out of list_of_rands and use it to get a random
% step distance
stepDist = randVector(stepCount)*(maxStep-minStep)+minStep;
% Advance to the next point
curPoint = prevPoint + stepDist*u;
% Reproject the current point and go to the next step
if mod(totalStepCount,10) == 0
if (full(max(max(abs(model.S*curPoint)))) > 1e-9)
curPoint = N*(N'*curPoint);
end
end
% Print out errors
if (mod(totalStepCount,2000)==0)
fprintf(fidErr,'%10.8f\t%10.8f\t',max(curPoint-maxFlux),max(minFlux-curPoint));
end
timeElapsed = cputime-t0;
% Print step information
if (mod(totalStepCount,5000)==0)
timePerStep = timeElapsed/totalStepCount;
if printLevel
fprintf('%d\t%d\t%d\t%8.2f\t%8.2f\n',i,pointCount,totalStepCount,timeElapsed/60,(totalCount-totalStepCount)*timePerStep/60);
end
end
overInd = find(curPoint > maxFlux);
underInd = find(curPoint < minFlux);
if (any((maxFlux-curPoint) < 0) | any((curPoint-minFlux) ...
< 0))
curPoint(overInd) = maxFlux(overInd);
curPoint(underInd) = minFlux(underInd);
end
if (mod(totalStepCount,2000)==0)
if printLevel
fprintf(fidErr,'%10.8f\n',full(max(max(abs(model.S*curPoint)))));
end
end
prevPoint = curPoint;
stepCount = stepCount + 1;
% Count the total number of steps
totalStepCount = totalStepCount + 1;
if printLevel
showprogress(totalStepCount/totalCount);
end
%recalculate the center point
centerPoint = ((nWrmup+totalStepCount)*centerPoint + curPoint)/(nWrmup+totalStepCount+1);
% Exit if time exceeded
if (timeElapsed > maxTime)
points(:,pointCount) = curPoint;
file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
save (file,'points');
save ACHR_last_point.mat curPoint
return;
end
end % Steps per point
% Add the current point to points
points(:,pointCount) = curPoint;
pointCount = pointCount + 1;
end % Points per cycle
% Save current points to a file
file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
save (file,'points');
end
\ No newline at end of file
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