From db2140dc392bf4105f012394afcf2557afa87c60 Mon Sep 17 00:00:00 2001 From: lvalcarcel Date: Wed, 13 Feb 2019 13:08:06 +0100 Subject: [PATCH 1/5] first ROOM implementation --- src/analysis/ROOM/ROOM.m | 130 ++++++++++++++++++ src/analysis/ROOM/linearROOM.m | 125 +++++++++++++++++ .../analysis/testROOM/testROOM.m | 79 +++++++++++ 3 files changed, 334 insertions(+) create mode 100644 src/analysis/ROOM/ROOM.m create mode 100644 src/analysis/ROOM/linearROOM.m create mode 100644 test/verifiedTests/analysis/testROOM/testROOM.m diff --git a/src/analysis/ROOM/ROOM.m b/src/analysis/ROOM/ROOM.m new file mode 100644 index 0000000000..5ee62fd190 --- /dev/null +++ b/src/analysis/ROOM/ROOM.m @@ -0,0 +1,130 @@ +function [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, fluxWT, rxnKO, varargin) +% Performs a MILP version of the ROOM (Regulatory on/off minimization of +% metabolic flux changes) approach +% +% USAGE: +% +% [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, WTflux, rxnKO, delta, epsilon, printLevel) +% +% INPUTS: +% model: Metabolic model +% fluxWT: Numeric array with flux distribution of wild type +% rxnKO: List of perturbations performed to the model +% (reactions that are eliminated) +% +% OPTIONAL INPUTS: +% delta: Multiplicative tol for flux change (Default = 0.03) +% epsilon: Additive tolerance for flux change (Default = 0.001) +% printLevel: Verbose output (Default = 1) +% +% OUTPUTS: +% fluxROOM: Flux distribution after ROOM calculation +% solutionROOM: Solution structure +% totalFluxDiff: Euclidean distance of ROOM objective, i.e. +% :math:`\sum (v_{wt}-v_{del})^2` +% +% Solve the following problem: +% +% .. math:: +% min ~&~ \sum y_{i} \\ +% ~&~ S_{del}v_{del} = 0 \\ +% ~&~ lb_{del} \leq v_{del} \leq ub_{del} \\ +% ~&~ for i=1:nRxns\\ +% ~&~ v_{i} - y_{i}(v_{max,i}-w_{wt,i}^u) \leq w_{wt,i}^u \\ +% ~&~ v_{i} - y_{i}(v_{min,i}-w_{wt,i}^l) \geq w_{wt,i}^l \\ +% ~&~ y_{i} \in {0,1} \\ +% ~&~ w_{wt,i}^u = w_{wt,i} + \delta |w_{wt,i}| + \epsilon \\ +% ~&~ w_{wt,i}^l = w_{wt,i} - \delta |w_{wt,i}| - \epsilon \\ +% +% NOTE:: +% +% The code here has been based on: +% Shlomi, T., Berkman, O., & Ruppin, E. (2005). Regulatory on/off +% minimization of metabolic flux changes after genetic perturbations. +% Proceedings of the National Academy of Sciences, 102(21), 7695-7700 +% +% .. Authors: +% - Luis V. Valcarcel, 23/01/2019, University of Navarra, CIMA & TECNUN School of Engineering. + + +p = inputParser; +% check required arguments +addRequired(p, 'model'); +addRequired(p, 'WTflux', @(x)isnumeric(x)&&isvector(x)); +addRequired(p, 'rxnKO', @(x)iscell(x)); +% Check optional arguments +addOptional(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x)); +addOptional(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x)); +addOptional(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x)); +% extract variables from parser +parse(p, model, fluxWT, rxnKO, varargin{:}); +delta = p.Results.delta; +epsilon = p.Results.epsilon; +printLevel = p.Results.printLevel; + +% LP solution tolerance +global CBT_LP_PARAMS +if (exist('CBT_LP_PARAMS', 'var')) + if isfield(CBT_LP_PARAMS, 'objTol') + tol = CBT_LP_PARAMS.objTol; + else + tol = 1e-6; + end +else + tol = 1e-6; +end + + +% Check the inputs +fluxWT = reshape(fluxWT,[],1); % reshape flux vector +assert(numel(model.rxns)==numel(fluxWT), 'This flux distribution has different number of reactions than the model') +assert(norm(model.S * fluxWT)< 10*tol, 'This flux distribution cannot exist in this model') +assert(all(ismember(rxnKO, model.rxns)), 'Some reactions are not in the model') + +% Eliminate almost-zero fluxes +fluxWT(abs(fluxWT) B', '2 B -> C + byp', '2 B + cof -> D + byp',... +'D -> E + cof', 'C + cof -> D', 'C -> E', ' -> A', 'E -> ', 'byp -> '}; +ReactionNames = {'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'b1', 'b2', 'b3'}; +lowerbounds = [0 0 0 0 0 0 0 0 0]; +upperbounds = [20, 20, 20, 20, 20, 20, 20, 20, 20]; +model = createModel(ReactionNames, ReactionNames, ReactionFormulas,... +'lowerBoundList', lowerbounds, 'upperBoundList', upperbounds); + +% Generate reference flux +WT_flux = [10, 5, 0, 0, 0, 5, 10, 5, 5]; + +% Check errors when missing argument +assert(verifyCobraFunctionError('ROOM', 'inputs', {model})) +assert(verifyCobraFunctionError('ROOM', 'inputs', {model, WT_flux})) +assert(verifyCobraFunctionError('ROOM', 'inputs', {model, WT_flux, 'a'})) + +% define the solver packages to be used to run this test +solverPkgs = solversToUse.MILP; + +predictedROOMsol = [10 5 0 5 5 0 10 5 5]'; + +% Test solving for different solvers +for k = 1:length(solverPkgs) + fprintf(' -- Running testROOM using the solver interface: %s ... ', solverPkgs{k}); + + solverOK = changeCobraSolver(solverPkgs{k}, 'all', 0); + + if solverOK + % Calculate ROOM and check solutions + fluxROOM = ROOM(model, WT_flux, {'v6'},1e-5,1e-5); + assert(norm(fluxROOM-predictedROOMsol)<1e-2) + + [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, WT_flux, {'v6'},1e-6,1e-6); + assert(norm(fluxROOM-predictedROOMsol)<1e-2) + assert(abs(totalFluxDiff^2-75)<1e-2) + + [fluxROOM2, solutionROOM2, totalFluxDiff2] = linearROOM(model, WT_flux, {'v6'},1e-6,1e-6); + assert(norm(fluxROOM-predictedROOMsol)<1e-2) + assert(abs(totalFluxDiff^2-75)<1e-2) + + else + warning('The test testROOM cannot run using the solver interface: %s. The solver interface is not installed or not configured properly.\n', solverPkgs{k}); + end + + % output a success message + fprintf('Done.\n'); +end + +% remove the file created during the test +if exist('CobraLPSolver.log','file') + delete('CobraLPSolver.log'); +end + +% Set seed to default value +rng('default') +% change the directory +cd(currentDir) From 042edc29e5d14b5c9103c362aa025a430a14a0c8 Mon Sep 17 00:00:00 2001 From: lvalcarcel Date: Thu, 14 Feb 2019 11:09:46 +0100 Subject: [PATCH 2/5] change from optional to paramenter --- src/analysis/ROOM/ROOM.m | 6 +++--- src/analysis/ROOM/linearROOM.m | 6 +++--- test/verifiedTests/analysis/testROOM/testROOM.m | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/analysis/ROOM/ROOM.m b/src/analysis/ROOM/ROOM.m index 5ee62fd190..0584971b34 100644 --- a/src/analysis/ROOM/ROOM.m +++ b/src/analysis/ROOM/ROOM.m @@ -53,9 +53,9 @@ addRequired(p, 'WTflux', @(x)isnumeric(x)&&isvector(x)); addRequired(p, 'rxnKO', @(x)iscell(x)); % Check optional arguments -addOptional(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x)); -addOptional(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x)); -addOptional(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x)); % extract variables from parser parse(p, model, fluxWT, rxnKO, varargin{:}); delta = p.Results.delta; diff --git a/src/analysis/ROOM/linearROOM.m b/src/analysis/ROOM/linearROOM.m index 9610aa5d86..10a37d4c4e 100644 --- a/src/analysis/ROOM/linearROOM.m +++ b/src/analysis/ROOM/linearROOM.m @@ -53,9 +53,9 @@ addRequired(p, 'WTflux', @(x)isnumeric(x)&&isvector(x)); addRequired(p, 'rxnKO', @(x)iscell(x)); % Check optional arguments -addOptional(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x)); -addOptional(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x)); -addOptional(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x)); +addParameter(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x)); % extract variables from parser parse(p, model, fluxWT, rxnKO, varargin{:}); delta = p.Results.delta; diff --git a/test/verifiedTests/analysis/testROOM/testROOM.m b/test/verifiedTests/analysis/testROOM/testROOM.m index 1e25b65240..e2bdb33386 100644 --- a/test/verifiedTests/analysis/testROOM/testROOM.m +++ b/test/verifiedTests/analysis/testROOM/testROOM.m @@ -49,14 +49,14 @@ if solverOK % Calculate ROOM and check solutions - fluxROOM = ROOM(model, WT_flux, {'v6'},1e-5,1e-5); + fluxROOM = ROOM(model, WT_flux, {'v6'},'delta', 1e-5, 'epsilon', 1e-5); assert(norm(fluxROOM-predictedROOMsol)<1e-2) - [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, WT_flux, {'v6'},1e-6,1e-6); + [fluxROOM, solutionROOM, totalFluxDiff] = ROOM(model, WT_flux, {'v6'},'delta',1e-6,'epsilon',1e-6); assert(norm(fluxROOM-predictedROOMsol)<1e-2) assert(abs(totalFluxDiff^2-75)<1e-2) - [fluxROOM2, solutionROOM2, totalFluxDiff2] = linearROOM(model, WT_flux, {'v6'},1e-6,1e-6); + [fluxROOM2, solutionROOM2, totalFluxDiff2] = linearROOM(model, WT_flux, {'v6'},'delta',1e-6,'epsilon',1e-6); assert(norm(fluxROOM-predictedROOMsol)<1e-2) assert(abs(totalFluxDiff^2-75)<1e-2) From 0752c2d1a9ed14dfb2ae0248d7cdb0ebb08a5a94 Mon Sep 17 00:00:00 2001 From: lvalcarcel Date: Mon, 11 Mar 2019 10:13:09 +0100 Subject: [PATCH 3/5] eliminate startsWith --- src/analysis/ROOM/ROOM.m | 2 +- src/analysis/ROOM/linearROOM.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/analysis/ROOM/ROOM.m b/src/analysis/ROOM/ROOM.m index 0584971b34..afe02f4b4f 100644 --- a/src/analysis/ROOM/ROOM.m +++ b/src/analysis/ROOM/ROOM.m @@ -109,7 +109,7 @@ ROOMmodel = changeRxnBounds(ROOMmodel, rxnKO, 0, 'b'); % change objective function ROOMmodel.c(:) = 0; -ROOMmodel.evarc(startsWith(ROOMmodel.evars,'binary_')) = 1; +ROOMmodel.evarc(cellfun(@length,regexp(ROOMmodel.evars,'^binary_'))==1) = 1; ROOMmodel.osenseStr = 'min'; MILPproblem = buildLPproblemFromModel(ROOMmodel); diff --git a/src/analysis/ROOM/linearROOM.m b/src/analysis/ROOM/linearROOM.m index 10a37d4c4e..9c931c5112 100644 --- a/src/analysis/ROOM/linearROOM.m +++ b/src/analysis/ROOM/linearROOM.m @@ -109,7 +109,7 @@ ROOMmodel = changeRxnBounds(ROOMmodel, rxnKO, 0, 'b'); % change objective function ROOMmodel.c(:) = 0; -ROOMmodel.evarc(:) = 1; +ROOMmodel.evarc(cellfun(@length,regexp(ROOMmodel.evars,'^binary_'))==1) = 1; ROOMmodel.osenseStr = 'min'; solutionROOM = optimizeCbModel(ROOMmodel); From fcbadf57031bb3a5c724af0aff48cde7e884d604 Mon Sep 17 00:00:00 2001 From: laurentheirendt Date: Fri, 15 Mar 2019 09:54:26 +0100 Subject: [PATCH 4/5] essential changes for special characters --- src/analysis/gMCS/GPR2models.m | 30 +++++----- src/analysis/gMCS/buildGmatrix.m | 30 +++++----- src/analysis/gMCS/calculateGeneMCS.m | 64 +++++++++++----------- src/analysis/gMCS/calculateMCS.m | 46 ++++++++-------- src/base/solvers/getAvailableGAMSSolvers.m | 18 +++--- 5 files changed, 94 insertions(+), 94 deletions(-) diff --git a/src/analysis/gMCS/GPR2models.m b/src/analysis/gMCS/GPR2models.m index 446accd96f..43bc9aded5 100644 --- a/src/analysis/gMCS/GPR2models.m +++ b/src/analysis/gMCS/GPR2models.m @@ -21,7 +21,7 @@ % * 1 - sequential % * 2+ - parallel % printLevel: show the reactions created in models. -% * 0 - shows nothing +% * 0 - shows nothing % * 1 - shows progress by reactions (default) % * 2+ - shows everything (reaction and network generation) % @@ -31,10 +31,10 @@ % rule for each reaction. % % .. Authors: -% - Iņigo Apaolaza, Aug 2017, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, Aug 2017, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, Aug 2017, University of Navarra, TECNUN School of Engineering. % - Francisco J. Planes, Aug 2017, University of Navarra, TECNUN School of Engineering. -% - Iņigo Apaolaza, April 2018, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, April 2018, University of Navarra, TECNUN School of Engineering. if (nargin < 5 || isempty(printLevel)) printLevel = 1; % Default is show progress @@ -100,13 +100,13 @@ % include metabolic reaction as objective metabolite model.mets{1} = RXN; % Add objective function - model.rxns{1} = RXN; + model.rxns{1} = RXN; model.S(1,1) = -1; - + % Check if there is any gene related to this reaction and add them genes = metabolic_model.genes(metabolic_model.rxnGeneMat(selected_rxns(i),:)>0); if length(genes)>1 - model.mets = vertcat(RXN, genes); % Include genes in the model + model.mets = vertcat(RXN, genes); % Include genes in the model % model = addMetabolite(model,genes); model.genes = genes; fp = FormulaParser(); @@ -130,12 +130,12 @@ model.rxns(demand_idx) = strcat('DM_',genes); model.rxns = model.rxns'; model.S(gene_idx,demand_idx) = eye(length(genes)); - + % converge isoforms (if neccesary) if ~isempty(separate_transcript) && length(genes)>1 model = convergeTranscripts2Gene(model, separate_transcript); end - + % Generate rest of fields idx_rxns = (1:length(model.rxns))'; idx_mets = (1:length(model.mets))'; @@ -152,7 +152,7 @@ model.lb = model.lb'; model.ub = model.ub'; model.rules = model.rules'; - + % Store the model networks{i} = model; end @@ -177,7 +177,7 @@ % g1 & g2 & g3 is one layer in the model, RXN = g1 + g2 + g3 % (g1 | g2) & g3 is a two layer model: RXN = (g1 | g2) + g3 % (g1 | g2) = g1 -% (g1 | g2) = g2 +% (g1 | g2) = g2 @@ -241,7 +241,7 @@ [~, rxn_idx] = find((model.S(node_idx(i),:)<0)); rxn_name = model.rxns(rxn_idx); % in this reaction, change (node -> rule) for (gene -> rule) - model = changeRxnMets(model, node_name, genes(gen), rxn_name); + model = changeRxnMets(model, node_name, genes(gen), rxn_name); end % remove stoichimetri of intermediate nodes model.S(node_idx,:)=0; @@ -263,10 +263,10 @@ % Some COBRA models have transcripts instead of genes. Most of the % transcripts have the same functionallity, so we can converge them as the % same gene. -% Ej: gene 10005.1 -% gene 10005.2 ==> gene 10005 +% Ej: gene 10005.1 +% gene 10005.2 ==> gene 10005 % gene 10005.3 -% +% % This is done for every gene of the model at exchange reaction level % seach demand reactions @@ -313,4 +313,4 @@ model.S(gene_idx,demand_idx) = eye(length(genes_unique)); -end \ No newline at end of file +end diff --git a/src/analysis/gMCS/buildGmatrix.m b/src/analysis/gMCS/buildGmatrix.m index aa1b07844c..40b8317dae 100644 --- a/src/analysis/gMCS/buildGmatrix.m +++ b/src/analysis/gMCS/buildGmatrix.m @@ -1,42 +1,42 @@ function [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix(model_name, model, separate_isoform, numWorkers, printLevel) % Build the G matrix required for the calculation of genetic Minimal Cut % Sets (gMCSs). -% +% % USAGE: -% +% % [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix(model_name, model_struct, separate_isoform, numWorkers, printLevel) -% +% % INPUTS: % model_name: Name of the metabolic model under study. % model_struct: Metabolic model structure (COBRA Toolbox format). % separate_isoform: Character used to discriminate different isoforms of a gene. -% +% % OPTIONAL INPUTS: % numWorkers: Maximum number of workers % * 0 - maximum provided by the system (automatic) % * 1 - sequential % * 2+ - parallel % printLevel: show the reactions created in models. -% * 0 - shows nothing +% * 0 - shows nothing % * 1 - shows progress by reactions (default) % * 2+ - shows everything (reaction and network generation) -% +% % OUTPUTS: % G: G matrix. % G_ind: Gene knockouts associated with each row in the G matrix. % related: Relationships among rows of the G matrix. % n_genes_KO: Number of genes whose knockout is added by each row of the G matrix. % G_time: Calculation times for each of the steps. -% +% % EXAMPLE: -% +% % [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix('Recon2.v04', modelR204, '.') -% +% % .. Authors: -% - Iņigo Apaolaza, 16/11/2017, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 16/11/2017, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, 19/11/2017, University of Navarra, TECNUN School of Engineering. % - Francisco J. Planes, 20/11/2017, University of Navarra, TECNUN School of Engineering. -% - Iņigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. time_a = tic; % Generate name for temporary folder @@ -194,7 +194,7 @@ for i = 1:n_rxns_or_and load([tmpFolderName filesep 'rxn_level_gMCSs_by_rxn' filesep 'rxn_level_gMCSs_' model_name '_rxn' num2str(pos_rxns_or_and(i)) '.mat']); n_act_mcs = length(act_mcs); - + for j = 1:n_act_mcs act_G_ind = act_mcs{j}; if ~iscell(act_G_ind) && isnan(act_G_ind) @@ -255,7 +255,7 @@ else act_G_ind = sort(tmp_G_ind{i}); end - + pos_equal = cellfun(@isequal, G_ind, repmat({act_G_ind}, length(G_ind), 1)); if sum(pos_equal) > 0 G(pos_equal, :) = G(pos_equal, :) + tmp_G(i, :); @@ -284,7 +284,7 @@ if sum(ismember(G_ind{pos(j)}, act_G_ind)) == n_act_G_ind G(pos(j), :) = G(pos(j), :) + G(i, :); end - end + end end G = double(G>0); @@ -336,4 +336,4 @@ disp('The G Matrix has been successfully calculated'); end rmdir(tmpFolderName, 's'); -end \ No newline at end of file +end diff --git a/src/analysis/gMCS/calculateGeneMCS.m b/src/analysis/gMCS/calculateGeneMCS.m index d266975f90..8cf8454f51 100644 --- a/src/analysis/gMCS/calculateGeneMCS.m +++ b/src/analysis/gMCS/calculateGeneMCS.m @@ -3,23 +3,23 @@ % available in CPLEX, namely cplex.populate(), with or without selecting a % given knockout, among all the genes included in the model or a given % subset of them. Apaolaza et al., 2017 (Nature Communications). -% +% % USAGE: -% +% % [gmcs, gmcs_time] = calculateGeneMCS(model_name, model_struct, n_gmcs, max_len_gmcs, options) -% +% % INPUTS: % model_name: Name of the metabolic model under study (in order to % identify the G matrix). % model_struct: Metabolic model structure (COBRA Toolbox format). % n_gmcs: Number of gMCSs to calculate. % max_len_gmcs: Number of genes in the largest gMCS to be calculated. -% +% % OPTIONAL INPUT: % options: Structure with fields: -% +% % * .KO - Selected gene knockout. Default: []. -% * .gene_set - Cell array containing the set of +% * .gene_set - Cell array containing the set of % genes among which the gMCSs are wanted to be calculated. % Default: [] (all genes are included). % * .timelimit - Time limit for the calculation of gMCSs @@ -28,13 +28,13 @@ % task to be disrupted. Default: 1e-3; % * .separate_transcript - Character used to discriminate % different transcripts of a gene. Default: ''. -% Example: separate_transcript = '' -% gene 10005.1 ==> gene 10005.1 +% Example: separate_transcript = '' +% gene 10005.1 ==> gene 10005.1 % gene 10005.2 ==> gene 10005.2 -% gene 10005.3 ==> gene 10005.3 -% separate_transcript = '.' -% gene 10005.1 -% gene 10005.2 ==> gene 10005 +% gene 10005.3 ==> gene 10005.3 +% separate_transcript = '.' +% gene 10005.1 +% gene 10005.2 ==> gene 10005 % gene 10005.3 % * .forceLength - 1 if the constraint limiting the % length of the gMCSs is to be active (recommended for @@ -45,12 +45,12 @@ % sequential, >1 = parallel. Default = 0; % * .printLevel - 1 if the process is wanted to be % shown on the screen, 0 otherwise. Default: 1. -% +% % OUTPUTS: % gmcs: Cell array containing the calculated gMCSs. -% gmcs_time: Calculation times of the different processes in +% gmcs_time: Calculation times of the different processes in % the algorithm. -% +% % EXAMPLE: % %With optional values % [gmcs, gmcs_time] = calculateGeneMCS('Recon2.v04', modelR204, 100, 10, options) @@ -62,16 +62,16 @@ % %options.separate_transcript = '.'; % %options.forceLength = 0 % %options.printLevel = 0 -% -% %Without optional values +% +% %Without optional values % [gmcs, gmcs_time] = calculateGeneMCS('ecoli_core_model', model, 100, 10) -% +% % .. Authors: -% - Iņigo Apaolaza, 30/01/2017, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 30/01/2017, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, 19/11/2017, University of Navarra, TECNUN School of Engineering. % - Francisco J. Planes, 20/11/2017, University of Navarra, TECNUN School of Engineering. % .. Revisions: -% - Iņigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, 17/04/2018, University of Navarra, TECNUN School of Engineering. % Check the installation of cplex @@ -188,12 +188,12 @@ if n_relations > 0 cell_related_1 = mat2cell(related(:, 1), ones(size(related, 1), 1), 1); cell_pos_set = mat2cell(pos_set, ones(n_poss_KO, 1), 1); - cell_related_1 = cellfun(@num2str, cell_related_1, 'UniformOutput', false); + cell_related_1 = cellfun(@num2str, cell_related_1, 'UniformOutput', false); cell_pos_set = cellfun(@num2str, cell_pos_set, 'UniformOutput', false); tmp_related_1 = cellfun(@ismember, cell_related_1, repmat({cell_pos_set}, size(related, 1), 1), 'UniformOutput', false); tmp_related_1 = logical(cell2mat(tmp_related_1)); related = related(tmp_related_1, :); - n_relations = size(related, 1); + n_relations = size(related, 1); pos_set(:, 2) = 1:length(pos_set); tmp_related = related(:); n_tmp_related = length(tmp_related); @@ -236,7 +236,7 @@ cons.forceLength = cons.forceBioCons(end)+1:cons.forceBioCons(end)+1; end n_cons = cons.forceLength(end); - + % Cplex - A matrix A = sparse(zeros(n_cons, n_vars)); A(cons.Ndual, var.u) = S'; @@ -297,7 +297,7 @@ % Cplex - sense of the optimization sense = 'minimize'; - + % Cplex - Introduce all data in a Cplex structure cplex = Cplex('geneMCS'); Model = struct(); @@ -322,7 +322,7 @@ % Cplex Parameters [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50); - [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); + [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ... sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ... sP.preprocessing.dual, sP.preprocessing.fill,... @@ -388,7 +388,7 @@ % CALCULATE gMCSs WITH A GIVEN KNOCKOUT % Select the row(s) in G_ind related to the KO under study n_G_ind = length(G_ind); - tmp = repmat({KO}, n_G_ind, 1); + tmp = repmat({KO}, n_G_ind, 1); dp = cellfun(@ismember, tmp, G_ind); % Define variables @@ -416,9 +416,9 @@ if n_relations > 0 cons.relations = cons.linearComb(end)+1:cons.linearComb(end)+n_relations; cons.forceLength = cons.relations(end)+1:cons.relations(end)+1; - else + else cons.forceLength = cons.linearComb(end)+1:cons.linearComb(end)+1; - end + end n_cons = cons.forceLength(end); % Cplex - A matrix @@ -459,7 +459,7 @@ try lhs(cons.relations) = 0; end if forceLength == 1 lhs(cons.forceLength) = 1; - end + end % Cplex - ub and lb vectors ub(var.u, 1) = inf; @@ -540,7 +540,7 @@ a(var_group.eps(ivar)) = 1; cplex.addIndicators(var_group.z(ivar), 0, a, 'L', 0); end - + % Cplex Indicators % z = 0 --> epsilon <= M for ivar = 1:length(var_group.z) @@ -552,7 +552,7 @@ % Cplex Parameters sP = struct(); [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50); - [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); + [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ... sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ... sP.preprocessing.dual, sP.preprocessing.fill,... @@ -618,4 +618,4 @@ n_time = size(gmcs_time, 1); gmcs_time{n_time+1, 1} = 'TOTAL gMCSs'; gmcs_time{n_time+1, 2} = toc(time_aa); -end \ No newline at end of file +end diff --git a/src/analysis/gMCS/calculateMCS.m b/src/analysis/gMCS/calculateMCS.m index 60119d0a59..2b9b2ecba3 100644 --- a/src/analysis/gMCS/calculateMCS.m +++ b/src/analysis/gMCS/calculateMCS.m @@ -4,22 +4,22 @@ % knockout, among all the reactions included in the model or a given subset % of them. Tobalina et al., 2016 (Bioinformatics); von Kamp and Klamt, 2014 % (PLoS Computational Biology). -% +% % USAGE: -% +% % [mcs, mcs_time] = populateMCS(model_struct, n_mcs, max_len_mcs, options) -% +% % INPUTS: % model_struct: Metabolic model structure (COBRA Toolbox format). % n_mcs: Number of MCSs to calculate. % max_len_mcs: Number of reactions in the largest MCS to be calculated. -% +% % OPTIONAL INPUT: % options: Structure with fields: -% +% % * .KO - Selected reaction knockout. Default: []. -% * .rxn_set - Cell array containing the set of -% reactions ('rxns') among which the MCSs are +% * .rxn_set - Cell array containing the set of +% reactions ('rxns') among which the MCSs are % wanted to be calculated. Default: [] (all reactions % are included). % * .timelimit - Time limit for the calculation of MCSs @@ -35,12 +35,12 @@ % sequential, >1 = parallel. Default = 0; % * .printLevel - 1 if the process is wanted to be % shown on the screen, 0 otherwise. Default: 1. -% +% % OUTPUTS: % mcs: Cell array containing the calculated MCSs. -% mcs_time: Calculation times of the different processes in +% mcs_time: Calculation times of the different processes in % the algorithm. -% +% % EXAMPLE: % %With optional values % [mcs, mcs_time] = calculateMCS(modelR204, 100, 10, options) @@ -51,16 +51,16 @@ % %options.target_b = 1e-4 % %options.forceLength = 0 % %options.printLevel = 0 -% -% %Without optional values +% +% %Without optional values % [mcs, mcs_time] = calculateMCS(model, 100, 10) -% +% % .. Authors: -% - Iņigo Apaolaza, 30/01/2017, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 30/01/2017, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, 19/11/2017, University of Navarra, TECNUN School of Engineering. % - Francisco J. Planes, 20/11/2017, University of Navarra, TECNUN School of Engineering. % .. Revisions: -% - Iņigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. +% - Inigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering. % - Luis V. Valcarcel, 17/04/2018, University of Navarra, TECNUN School of Engineering. % Check the installation of cplex @@ -227,7 +227,7 @@ % Cplex - sense of the optimization sense = 'minimize'; - + % Cplex - Introduce all data in a Cplex structure cplex = Cplex('MCS'); Model = struct(); @@ -253,7 +253,7 @@ % Cplex Parameters sP = struct(); [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50); - [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); + [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ... sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ... sP.preprocessing.dual, sP.preprocessing.fill,... @@ -261,7 +261,7 @@ sP.preprocessing.presolve, sP.preprocessing.reduce,..., ... sP.preprocessing.relax, sP.preprocessing.symmetry] = deal(50, 1, 2, 1, 1, 50, 1, 50, 1, 3, 1, 1); cplex = setCplexParam(cplex, sP); - + if printLevel == 0 cplex.DisplayFunc = []; end @@ -311,7 +311,7 @@ mcs_time{n_time+1, 1} = 'TOTAL MCSs'; mcs_time{n_time+1, 2} = toc(time_aa); return; - end + end end try save('tmp.mat', 'mcs', 'mcs_time'); end try largest_mcs = max(cellfun(@length, mcs)); end @@ -458,7 +458,7 @@ a(var_group.eps(ivar)) = 1; cplex.addIndicators(var_group.z(ivar), 0, a, 'L', 0); end - + % Cplex Indicators % z = 0 --> epsilon <= M for ivar = 1:length(var_group.z) @@ -470,7 +470,7 @@ % Cplex Parameters sP = struct(); [sP.mip.tolerances.integrality, sP.mip.strategy.heuristicfreq, sP.mip.strategy.rinsheur] = deal(integrality_tolerance, 1000, 50); - [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); + [sP.emphasis.mip, sP.output.clonelog, sP.timelimit, sP.threads] = deal(4, -1, max(10, timelimit), numWorkers); [sP.preprocessing.aggregator, sP.preprocessing.boundstrength, ... sP.preprocessing.coeffreduce, sP.preprocessing.dependency, ... sP.preprocessing.dual, sP.preprocessing.fill,... @@ -527,7 +527,7 @@ mcs_time{n_time+1, 1} = 'TOTAL MCSs'; mcs_time{n_time+1, 2} = toc(time_aa); return; - end + end end try save('tmp.mat', 'mcs', 'mcs_time'); end try largest_mcs = max(cellfun(@length, mcs)); end @@ -536,4 +536,4 @@ n_time = size(mcs_time, 1); mcs_time{n_time+1, 1} = 'TOTAL MCSs'; mcs_time{n_time+1, 2} = toc(time_aa); -end \ No newline at end of file +end diff --git a/src/base/solvers/getAvailableGAMSSolvers.m b/src/base/solvers/getAvailableGAMSSolvers.m index bea395def1..5c325c6dfa 100644 --- a/src/base/solvers/getAvailableGAMSSolvers.m +++ b/src/base/solvers/getAvailableGAMSSolvers.m @@ -2,8 +2,8 @@ % This function return the GAMS solvers which are available in your system % according to the license. % -% USAGE: -% +% USAGE: +% % [summaryTable, booleanTable, problemTypes, solvers] = getAvailableGAMSSolvers % % OUTPUTS: @@ -26,7 +26,7 @@ % solvers Type: cell array of strings % Description: list of solvers available in GAMS % -% .. Author: - Sebastián Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl +% .. Author: - Sebastiian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl gamsPath = which('gams'); % verify that gams is installed if isempty(gamsPath) @@ -48,10 +48,10 @@ %run licememo.gms run = system('gams licememo'); -%if licememo.gms ran correctly. +%if licememo.gms ran correctly. if run == 0 % read the content of the file generated by licememo.gms to find the - % GAMS solvers available + % GAMS solvers available fileID = fopen('licememo.txt', 'r'); summaryTable = cell(2, 2); booleanTable = zeros(120, 30); @@ -68,7 +68,7 @@ end end nProblemType = length(problemTypes); - + nSolvers = 0; tline = fgetl(fileID); while ischar(tline) @@ -93,14 +93,14 @@ end tline = fgetl(fileID); end - + booleanTable = booleanTable(1:nSolvers, 1:nProblemType); summaryTable = [{''}, problemTypes; [solvers', summaryTable]]; fclose(fileID); %delete autogenerated file delete('licememo.txt') else - + summaryTable ={}; booleanTable = 0; problemTypes = {}; @@ -110,4 +110,4 @@ delete('licememo.gms') delete('licememo.lst') -end \ No newline at end of file +end From 18f9051a2f8362d7ba68725a8a1d10bc8a6e8d8b Mon Sep 17 00:00:00 2001 From: laurentheirendt Date: Fri, 15 Mar 2019 10:05:10 +0100 Subject: [PATCH 5/5] adding index files for ROOM --- docs/source/modules/analysis/ROOM/index.rst | 9 +++++++++ docs/source/modules/analysis/index.rst | 1 + 2 files changed, 10 insertions(+) create mode 100644 docs/source/modules/analysis/ROOM/index.rst diff --git a/docs/source/modules/analysis/ROOM/index.rst b/docs/source/modules/analysis/ROOM/index.rst new file mode 100644 index 0000000000..e4b44f69a2 --- /dev/null +++ b/docs/source/modules/analysis/ROOM/index.rst @@ -0,0 +1,9 @@ + +.. _ROOM: + + +ROOM +---- + +.. automodule:: src.analysis.ROOM + :members: diff --git a/docs/source/modules/analysis/index.rst b/docs/source/modules/analysis/index.rst index df2385d3ad..86f2f58438 100644 --- a/docs/source/modules/analysis/index.rst +++ b/docs/source/modules/analysis/index.rst @@ -42,6 +42,7 @@ Analysis rFBA/index rMTA/index robustnessAnalysis/index + ROOM/index rumba/index sampling/index sparseFBA/index