diff --git a/.gitattributes b/.gitattributes index 3e0ca4ab..7f4d5667 100644 --- a/.gitattributes +++ b/.gitattributes @@ -10,3 +10,4 @@ mafft binary .gitattributes export-ignore .gitignore export-ignore .github export-ignore +*.mat filter=lfs diff=lfs merge=lfs -text diff --git a/INIT/runINIT.m b/INIT/runINIT.m index ca064f50..f20fd586 100644 --- a/INIT/runINIT.m +++ b/INIT/runINIT.m @@ -1,4 +1,4 @@ -function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops) +function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params) % runINIT % Generates a model using the INIT algorithm, based on proteomics and/or % transcriptomics and/or metabolomics and/or metabolic tasks. This is the @@ -38,6 +38,7 @@ % problem significantly more computationally intensive to % solve (two more integer constraints per reversible reaction) % (optional, default false) +% params parameter structure for use by optimizeProb % % outModel the resulting model structure % deletedRxns reactions which were deleted by the algorithm diff --git a/core/FSEOF.m b/core/FSEOF.m index 585abf1f..e60740ce 100755 --- a/core/FSEOF.m +++ b/core/FSEOF.m @@ -119,7 +119,12 @@ A2=char(model.rxns(num)); %enzyme ID A3=char(model.rxnNames(num)); %enzyme Name if isfield(model,'subSystems') && ~isempty(model.subSystems{num}); - A4=char(strjoin(model.subSystems{num,1},';')); %Subsystems + if ~any(cellfun(@(x) iscell(x), model.subSystems)); + subSys = cellfun(@(x) {x}, model.subSystems, 'uni', 0); + else + subSys = model.subSystems; + end + A4=char(strjoin(subSys{num},';')); %Subsystems else A4=''; end diff --git a/core/addExchangeRxns.m b/core/addExchangeRxns.m index 31df9b80..831e888b 100755 --- a/core/addExchangeRxns.m +++ b/core/addExchangeRxns.m @@ -63,8 +63,8 @@ end if isfield(model,'subSystems') fillerSub = filler; - if iscell(model.subSystems(1,1)) - fillerSub = repmat({fillerSub},numel(J),1); + if any(cellfun(@(x) iscell(x), model.subSystems)) + fillerSub(:)={{''}}; end model.subSystems=[model.subSystems;fillerSub]; end diff --git a/core/addRxns.m b/core/addRxns.m index f23222e1..3ce890c6 100755 --- a/core/addRxns.m +++ b/core/addRxns.m @@ -227,10 +227,12 @@ nOldRxns=numel(model.rxns); filler=cell(nRxns,1); filler(:)={''}; -cellfiller=cellfun(@(x) cell(0,0),filler,'UniformOutput',false); +cellfiller=cell(nRxns,1); +cellfiller(:)={{''}}; largeFiller=cell(nOldRxns,1); largeFiller(:)={''}; -celllargefiller=cellfun(@(x) cell(0,0),largeFiller,'UniformOutput',false); +largeCellFiller=cell(nOldRxns,1); +largeCellFiller(:)={{''}}; %***Add everything to the model except for the equations. if numel(rxnsToAdd.equations)~=nRxns @@ -350,24 +352,45 @@ end if isfield(rxnsToAdd,'subSystems') - if numel(rxnsToAdd.subSystems)~=nRxns - EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns'; - dispEM(EM); + % Has to be cell array + if ischar(rxnsToAdd.subSystems) + rxnsToAdd.subSystems = {rxnsToAdd.subSystems}; + end + % If all nested cells are 1x1, then unnest + if all(cellfun(@(x) iscell(x) && isscalar(x), rxnsToAdd.subSystems)) + rxnsToAdd.subSystems = transpose([rxnsToAdd.subSystems{:}]); end - for i=1:numel(rxnsToAdd.subSystems) - if ischar(rxnsToAdd.subSystems{i}) - rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i); + % Cell array should now be as simple as possible. Check if it is nested + subSysRxnsNested = any(cellfun(@(x) iscell(x), rxnsToAdd.subSystems)); + if isfield(newModel,'subSystems') + subSysModelNested = any(cellfun(@(x) iscell(x), newModel.subSystems)); + else + subSysModelNested = subSysRxnsNested; + if subSysRxnsNested + newModel.subSystems=largeCellFiller; + else + newModel.subSystems=largeFiller; end end - %Fill with standard if it doesn't exist - if ~isfield(newModel,'subSystems') - newModel.subSystems=celllargefiller; + if subSysRxnsNested && ~subSysModelNested + % Make all existing subSystems nested + newModel.subSystems = cellfun(@(x) {x}, newModel.subSystems, 'uni', 0); + elseif ~subSysRxnsNested && subSysModelNested + rxnsToAdd.subSystems = cellfun(@(x) {x}, rxnsToAdd.subSystems, 'uni', 0); + end + if numel(rxnsToAdd.subSystems)~=nRxns + EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns'; + dispEM(EM); end newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)]; else - %Fill with standard if it doesn't exist + %Fill with standard if it does not exist if isfield(newModel,'subSystems') - newModel.subSystems=[newModel.subSystems;cellfiller]; + if any(cellfun(@(x) iscell(x), newModel.subSystems)) + newModel.subSystems=[newModel.subSystems;cellfiller]; + else + newModel.subSystems=[newModel.subSystems;filler]; + end end end if isfield(rxnsToAdd,'rxnMiriams') diff --git a/core/addTransport.m b/core/addTransport.m index 160f9fd2..4b88f022 100755 --- a/core/addTransport.m +++ b/core/addTransport.m @@ -134,11 +134,15 @@ model.eccodes=[model.eccodes;filler]; end if isfield(model,'subSystems') - ssFiller=filler; + isNested = any(cellfun(@(x) iscell(x), model.subSystems)); + ssFiller = filler; if isRev==1 - ssFiller(:)={{['Transport between ' fromComp ' and ' toComps{i}]}}; + ssFiller(:) = {['Transport between ' fromComp ' and ' toComps{i}]}; else - ssFiller(:)={{['Transport from ' fromComp ' to ' toComps{i}]}}; + ssFiller(:) = {['Transport from ' fromComp ' to ' toComps{i}]}; + end + if isNested + ssFiller = cellfun(@(x) {x}, ssFiller, 'uni', 0); end model.subSystems=[model.subSystems;ssFiller]; end diff --git a/core/checkModelStruct.m b/core/checkModelStruct.m index 3fb1e15f..eb9f8e6d 100755 --- a/core/checkModelStruct.m +++ b/core/checkModelStruct.m @@ -175,11 +175,11 @@ function checkModelStruct(model,throwErrors,trimWarnings) end end if isfield(model,'subSystems') - for i=1:numel(model.subSystems) - if ~iscell(model.subSystems{i,1}) - EM='The "subSystems" field must be a cell array'; + isNested = any(cellfun(@(x) iscell(x), model.subSystems)); + isCellStr = any(cellfun(@(x) ischar(x), model.subSystems)); + if ~xor(isNested,isCellStr) + EM='The "subSystems" field must be a cell array of chars, *or* a cell array of cell arrays of chars'; dispEM(EM,throwErrors); - end end end if isfield(model,'eccodes') diff --git a/core/compareMultipleModels.m b/core/compareMultipleModels.m index bb34169f..596d0957 100755 --- a/core/compareMultipleModels.m +++ b/core/compareMultipleModels.m @@ -92,16 +92,6 @@ fprintf('*** Done \n\n') -%% Flatten models' subSystems field -% Convert from cell array of cells to cell array of strings -% NOTE: this function currently only recognizes one subSystem per reaction; -% additional subSystems will be ignored! -for i = 1:numel(models) - cells = cellfun(@iscell,models{i}.subSystems); - models{i}.subSystems(cells) = cellfun(@(s) s{1}, models{i}.subSystems(cells), 'UniformOutput', false); -end - - %% Compare models structure & function based on high-dimensional methods % Compare number of reactions in each subsystem in each model using a heatmap field = 'subSystems'; @@ -110,6 +100,16 @@ fprintf('\nWARNING: At least one model does not contain the field "subSystems". \n') fprintf(' Skipping subsystem comparison. \n\n') else + %% Flatten model subSystems field + % Convert from cell array of cells to cell array of strings + % NOTE: this function currently only recognizes one subSystem per reaction; + % additional subSystems will be ignored! + for i = 1:numel(models) + if any(cellfun(@(x) iscell(x), models{i}.subSystems)); + cells = cellfun(@iscell,models{i}.subSystems); + models{i}.subSystems(cells) = cellfun(@(s) s{1}, models{i}.subSystems(cells), 'UniformOutput', false); + end + end [id,compMat] = compareModelField(models,field); compStruct.subsystems.ID = id; compStruct.subsystems.matrix = compMat; diff --git a/core/copyToComps.m b/core/copyToComps.m index a285e650..042b4c0b 100755 --- a/core/copyToComps.m +++ b/core/copyToComps.m @@ -25,24 +25,26 @@ % % Usage: model=copyToComps(model,toComps,rxns,deleteOriginal,compNames,compOutside) -if nargin<3 - rxns=model.rxns; -elseif ~islogical(rxns) && ~isnumeric(rxns) - rxns=convertCharArray(rxns); +arguments + model (1,1) struct + toComps {emptyOrTextOrCellOfText} + rxns = model.rxns + deleteOriginal {emptyOrLogicalScalar} = false + compNames {emptyOrTextOrCellOfText} = toComps + compOutside {emptyOrTextOrCellOfText} = ''; end -if nargin<4 - deleteOriginal=false; + +if nargin >= 3 && ~islogical(rxns) && ~isnumeric(rxns) + rxns = convertCharArray(rxns); end -if nargin<5 - compNames=toComps; -else +if nargin >= 5 compNames=convertCharArray(compNames); end -if nargin<6 - compOutside=cell(numel(toComps),1); - compOutside(:)={''}; -else +if nargin >= 6 compOutside=convertCharArray(compOutside); + if length(compOutside) ~= length(compNames) + error('compOutside and compNames should be of equal size.'); + end end originalID=model.id; @@ -79,15 +81,20 @@ modelToAdd.compMiriams=modelToAdd.compMiriams(J); end modelToAdd.metComps=ones(numel(modelToAdd.mets),1); - + if isfield(modelToAdd,'metFrom') + modelToAdd = rmfield(modelToAdd,'metFrom'); + end + if isfield(modelToAdd,'rxnFrom') + modelToAdd = rmfield(modelToAdd,'rxnFrom'); + end + if isfield(modelToAdd,'geneFrom') + modelToAdd = rmfield(modelToAdd,'geneFrom'); + end + %Merge the models - model=mergeModels({model;modelToAdd},'metNames'); + model=mergeModels({model;modelToAdd},'metNames',[],true); end -model=rmfield(model,'rxnFrom'); -model=rmfield(model,'metFrom'); -model=rmfield(model,'geneFrom'); - if deleteOriginal==true model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments end diff --git a/core/findRAVENroot.m b/core/findRAVENroot.m index 10c28ce4..17e4b878 100755 --- a/core/findRAVENroot.m +++ b/core/findRAVENroot.m @@ -1,23 +1,27 @@ function [ravenPath, prevDir] = findRAVENroot() % findRAVENroot -% Finds the root of the RAVEN directory, by searching for the path to +% Finds the root of the RAVEN directory, first by by searching for the path to % RAVEN2.png. Can also record the current directory, in case a function will % use the ravenPath to navigate to a precise folder, and it should return to % the previous directory afterwards. See e.g. optimizeProb calling glpk. ST=dbstack('-completenames'); prevDir = pwd(); -ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file; -rootFound = 0; -while rootFound == 0 - isRoot = isfile(fullfile(ravenPath,'RAVEN2.png')); - if isRoot - rootFound = 1; - else - ravenPathOld = ravenPath; - ravenPath = fileparts(ravenPath); - if strcmp(ravenPathOld,ravenPath) - error('Cannot find the RAVEN root directory. Make sure you have not removed the RAVEN2.png file from your RAVEN installation.') +if ispref('RAVEN','ravenPath') + ravenPath = getpref('RAVEN','ravenPath'); +else + ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file; + rootFound = 0; + while rootFound == 0 + isRoot = isfile(fullfile(ravenPath,'RAVEN2.png')); + if isRoot + rootFound = 1; + else + ravenPathOld = ravenPath; + ravenPath = fileparts(ravenPath); + if strcmp(ravenPathOld,ravenPath) + error('Cannot find the RAVEN root directory. Make sure you have not removed the RAVEN2.png file from your RAVEN installation.') + end end end end \ No newline at end of file diff --git a/core/fitTasks.m b/core/fitTasks.m index 9963f9e7..f823d459 100755 --- a/core/fitTasks.m +++ b/core/fitTasks.m @@ -34,7 +34,7 @@ % the result. % % Usage: [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,... -% rxnScores,taskStructure,params) +% rxnScores,taskStructure) if nargin<4 printOutput=true; @@ -48,9 +48,6 @@ if nargin<6 taskStructure=[]; end -if nargin<7 - params=[]; -end if isempty(taskStructure) && ~isfile(inputFile) error('Task file %s cannot be found',string(inputFile)); @@ -266,9 +263,9 @@ %Only do gap-filling if it cannot be solved failed=false; try - [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params); + [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores); if exitFlag==-2 - EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n']; + EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality.\n']; dispEM(EM,false); end catch diff --git a/core/mergeModels.m b/core/mergeModels.m index 29faff96..d7ef3669 100755 --- a/core/mergeModels.m +++ b/core/mergeModels.m @@ -1,37 +1,76 @@ -function model=mergeModels(models,metParam,supressWarnings) +function model=mergeModels(models,metParam,supressWarnings,copyToComps) % mergeModels % Merges models into one model structure. Reactions are added without any % checks, so duplicate reactions might appear. Metabolites are matched by % their name and compartment (metaboliteName[comp]), while genes are % matched by their name. % +% Input: % models a cell array with model structures -% metParam string specifying whether to refer to metabolite name -% (metNames) or ID (mets) for matching (default, metNames) -% supressWarnings true if warnings should be supressed (optional, default -% false) +% metParam string metabolite name ('metNames') or ID ('mets') are +% used for matching (optional, default 'metNames') +% supressWarnings logical whether warnings should be supressed (optional, +% default false) +% copyToComps logical whether mergeModels is run via copyToComps +% (optional, default false) % -% model a model structure with the merged model. Follows the structure -% of normal models but also has 'rxnFrom/metFrom/geneFrom' fields -% to indicate from which model each reaction/metabolite/gene was -% taken +% Output: +% model a model structure with the merged model. Follows the +% structure of normal models but also has 'rxnFrom/ +% metFrom/geneFrom' fields to indicate from which model +% each reaction/metabolite/gene was taken. If the model +% already has 'rxnFrom/metFrom/geneFrom' fields, then +% these fields are not modified. % % Usage: model=mergeModels(models) +arguments + models; + metParam {emptyOrTextScalar} = "metNames" + supressWarnings {emptyOrLogicalScalar} = false + copyToComps {emptyOrLogicalScalar} = false +end + +metParam = char(metParam); + %Just return the model if numel(models)<=1 model=models{1}; return; end -if nargin<2 - metParam='metNames'; -else - metParam=char(metParam); -end +hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models); +hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models); +hasRxnFrom = cellfun(@(s) isfield(s,'rxnFrom'), models); -if nargin<3 - supressWarnings=false; +for i = 1:numel(models) + if copyToComps + if hasMetFrom(1) + models{2}.metFrom = repmat({''},numel(models{i}.mets),1); + end + elseif ~any(hasMetFrom) + models{i}.metFrom = repmat({models{i}.id},numel(models{i}.mets),1); + elseif ~hasMetFrom(i) + models{i}.metFrom = repmat({''},numel(models{i}.mets),1); + end + if copyToComps + if hasRxnFrom(1) + models{2}.rxnFrom = repmat({''},numel(models{i}.rxns),1); + end + elseif ~any(hasRxnFrom) + models{i}.rxnFrom = repmat({models{i}.id},numel(models{i}.rxns),1); + elseif ~hasRxnFrom(i) + models{i}.rxnFrom = repmat({''},numel(models{i}.rxns),1); + end + if copyToComps + if hasGeneFrom(1) + models{2}.geneFrom = repmat({''},numel(models{i}.genes),1); + end + elseif ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models)) + models{i}.geneFrom = repmat({models{i}.id},numel(models{i}.genes),1); + elseif ~hasGeneFrom(i) + models{i}.geneFrom = repmat({''},numel(models{i}.genes),1); + end end %Add new functionality in the order specified in models @@ -39,19 +78,15 @@ model.id='MERGED'; model.name=''; -model.rxnFrom=cell(numel(models{1}.rxns),1); -model.rxnFrom(:)={models{1}.id}; -model.metFrom=cell(numel(models{1}.mets),1); -model.metFrom(:)={models{1}.id}; -if isfield(models{1},'genes') - model.geneFrom=cell(numel(models{1}.genes),1); - model.geneFrom(:)={models{1}.id}; -end - if isfield(model,'equations') model=rmfield(model,'equations'); end +for i=1:numel(models) + if isfield(models{i},'subSystems') + models{i}.subSystems = cellfun(@(x) {x}, models{i}.subSystems, 'uni', 0); + end +end for i=2:numel(models) %Add the model id to the rxn id id it already exists in the model (id %have to be unique) This is because it makes a '[]' string if no new @@ -78,15 +113,15 @@ end %Add all static stuff - rxnFrom=cell(numel(models{i}.rxns),1); - rxnFrom(:)={models{i}.id}; - model.rxnFrom=[model.rxnFrom;rxnFrom]; - model.rxns=[model.rxns;models{i}.rxns]; - model.rxnNames=[model.rxnNames;models{i}.rxnNames]; - model.lb=[model.lb;models{i}.lb]; - model.ub=[model.ub;models{i}.ub]; - model.c=[model.c;models{i}.c]; - model.rev=[model.rev;models{i}.rev]; + if any(hasRxnFrom) || (~copyToComps && ~any(hasRxnFrom)) + model.rxnFrom = [model.rxnFrom; models{i}.rxnFrom]; + end + model.rxns = [model.rxns; models{i}.rxns]; + model.rxnNames = [model.rxnNames; models{i}.rxnNames]; + model.lb = [model.lb; models{i}.lb]; + model.ub = [model.ub; models{i}.ub]; + model.c = [model.c; models{i}.c]; + model.rev = [model.rev; models{i}.rev]; if isfield(models{i},'subSystems') if isfield(model,'subSystems') @@ -253,7 +288,6 @@ if strcmpi(metParam,'mets') %Get the new metabolites from matching the models. Metabolites are matched by metabolite ID (model.mets). - oldMetComps=model.comps(model.metComps); oldMets=model.mets; if ~isempty(models{i}.mets) @@ -287,12 +321,12 @@ end %Add static info on the metabolites - metFrom=cell(numel(metsToAdd),1); - metFrom(:)={models{i}.id}; - model.metFrom=[model.metFrom;metFrom]; - model.mets=[model.mets;models{i}.mets(metsToAdd)]; - model.metNames=[model.metNames;models{i}.metNames(metsToAdd)]; - model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))]; + if any(hasMetFrom) + model.metFrom = [model.metFrom; models{i}.metFrom(metsToAdd)]; + end + model.mets = [model.mets; models{i}.mets(metsToAdd)]; + model.metNames = [model.metNames; models{i}.metNames(metsToAdd)]; + model.b = [model.b; zeros(numel(metsToAdd),size(model.b,2))]; if isfield(model,'unconstrained') if isfield(models{i},'unconstrained') @@ -481,13 +515,13 @@ if isfield(models{i},'genes') if ~isfield(model,'genes') %If there was no gene info before - model.genes=models{i}.genes; - model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat]; - emptyGene=cell(numel(model.rxns),1); - emptyGene(:)={''}; - model.grRules=[emptyGene;models{i}.grRules]; - model.geneFrom=cell(numel(models{i}.genes),1); - model.geneFrom(:)={models{i}.id}; + model.genes = models{i}.genes; + model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat]; + emptyGene = repmat({''},numel(model.rxns),1); + model.grRules = [emptyGene;models{i}.grRules]; + if any(hasGeneFrom) + model.geneFrom = models{i}.geneFrom; + end if isfield(models{i},'geneShortNames') model.geneShortNames=models{i}.geneShortNames; @@ -513,11 +547,9 @@ %Only add extra gene info on new genes. This might not be %correct and should be changed later... if ~isempty(genesToAdd) - model.genes=[model.genes;models{i}.genes(genesToAdd)]; - emptyGene=cell(numel(genesToAdd),1); - emptyGene(:)={models{i}.id}; - model.geneFrom=[model.geneFrom;emptyGene]; - model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))]; + model.genes = [model.genes; models{i}.genes(genesToAdd)]; + model.geneFrom = [model.geneFrom; models{i}.geneFrom(genesToAdd)]; + model.rxnGeneMat = [model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))]; if isfield(models{i},'geneShortNames') if isfield(model,'geneShortNames') @@ -585,9 +617,9 @@ end %Remap the genes from the new model. The same thing as with - %mets; this is a wasteful way to do it but I don't care right + %mets; this is a wasteful way to do it but I do not care right %now - [a, b]=ismember(models{i}.genes,model.genes); + a = ismember(models{i}.genes,model.genes); %Just a check if ~all(a) @@ -609,4 +641,10 @@ [grRules,rxnGeneMat] = standardizeGrRules(model,true); model.grRules = grRules; model.rxnGeneMat = rxnGeneMat; +%Flatten subSystems if possible +if isfield(model,'subSystems') + if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems)) + model.subSystems = transpose([model.subSystems{:}]); + end +end end diff --git a/core/parseTaskList.m b/core/parseTaskList.m index 614710c5..e6cafca0 100755 --- a/core/parseTaskList.m +++ b/core/parseTaskList.m @@ -2,10 +2,11 @@ % parseTaskList % Parses a task list file. % -% inputFile a task list in Excel format. The file must contain a -% sheet named TASKS, which in turn may contain the -% following column headers (note, all rows starting with -% a non-empty cell are removed. The first row after that +% inputFile a task list in either Excel (*.xlsx, with a sheet named +% TASKS with all relevant content) or tab-delimited +% (*.txt) format. The file may contain the following +% column headers (note, all rows starting with a +% non-empty cell are removed. The first row after that % is considered the headers): % ID % the only required header. Each task must have a diff --git a/core/predictLocalization.m b/core/predictLocalization.m index 2bfff833..1908caa6 100755 --- a/core/predictLocalization.m +++ b/core/predictLocalization.m @@ -621,7 +621,7 @@ outModel.compNames(2)=GSS.compartments(1); end end -outModel.compNames=[outModel.compNames;GSS.compartments(2:end)']; +outModel.compNames=[outModel.compNames;GSS.compartments(2:end)]; %Ugly little loop for i=1:numel(GSS.compartments)-1 @@ -726,7 +726,11 @@ outModel.rxnMiriams=[outModel.rxnMiriams;{[]}]; end if isfield(outModel,'subSystems') - outModel.subSystems=[outModel.subSystems;{{'Inferred transport reactions'}}]; + if any(cellfun(@(x) iscell(x), outModel.subSystems)) + outModel.subSystems=[outModel.subSystems;{{'Inferred transport reactions'}}]; + else + outModel.subSystems=[outModel.subSystems;{'Inferred transport reactions'}]; + end end if isfield(outModel,'eccodes') outModel.eccodes=[outModel.eccodes;{''}]; diff --git a/core/setParam.m b/core/setParam.m index 734853b7..70705f3f 100755 --- a/core/setParam.m +++ b/core/setParam.m @@ -68,7 +68,7 @@ params(~Lia)=[]; indexes(~Lia)=[]; paramType(~Lia)=[]; - dispEM('Reactions not present in model, will be ignored:',false,rxnLise(~Lia)); + dispEM('Reactions not present in model, will be ignored:',false,rxnList(~Lia)); end %Change the parameters diff --git a/core/sortModel.m b/core/sortModel.m index 8d329fbc..7dc00d7e 100755 --- a/core/sortModel.m +++ b/core/sortModel.m @@ -65,6 +65,9 @@ subsystemsUnique=''; subsystemsConcatenated=''; for i=1:numel(model.subSystems) + if ~iscell(model.subSystems{i}) + model.subSystems{i} = {model.subSystems{i}}; + end subsystemsConcatenated{i,1}=strjoin(model.subSystems{i,1},';'); if ~isempty(model.subSystems{i,1}) for j=1:numel(model.subSystems{i,1}) diff --git a/doc/INIT/getINITModel.html b/doc/INIT/getINITModel.html index 739882c2..d0913510 100644 --- a/doc/INIT/getINITModel.html +++ b/doc/INIT/getINITModel.html @@ -133,7 +133,7 @@



runINIT @@ -66,6 +66,7 @@DESCRIPTION
-CROSS-REFERENCE INFORMATION
SOURCE CODE
0001 function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops) ++0041 % params parameter structure for use by optimizeProb +0042 % +0043 % outModel the resulting model structure +0044 % deletedRxns reactions which were deleted by the algorithm +0045 % metProduction array that indicates which of the +0046 % metabolites in presentMets that could be +0047 % produced +0048 % -2: metabolite name not found in model +0049 % -1: metabolite found, but it could not be produced +0050 % 1: metabolite could be produced +0051 % fValue objective value (sum of (the negative of) +0052 % reaction scores for the included reactions and +0053 % prodWeight*number of produced metabolites) +0054 % +0055 % This function is the actual implementation of the algorithm. See +0056 % getINITModel for a higher-level function for model reconstruction. See +0057 % PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the +0058 % implementation. +0059 % +0060 % Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,... +0061 % rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,... +0062 % noRevLoops,params) +0063 +0064 if nargin<2 +0065 rxnScores=zeros(numel(model.rxns),1); +0066 end +0067 if isempty(rxnScores) +0068 rxnScores=zeros(numel(model.rxns),1); +0069 end +0070 if nargin<3 || isempty(presentMets) +0071 presentMets={}; +0072 else +0073 presentMets=convertCharArray(presentMets); +0074 end +0075 if nargin<4 || isempty(essentialRxns) +0076 essentialRxns={}; +0077 else +0078 essentialRxns=convertCharArray(essentialRxns); +0079 end +0080 if nargin<5 || isempty(prodWeight) +0081 prodWeight=0.5; +0082 end +0083 if nargin<6 +0084 allowExcretion=false; +0085 end +0086 if nargin<7 +0087 noRevLoops=false; +0088 end +0089 if nargin<8 +0090 params=[]; +0091 end +0092 +0093 if numel(presentMets)~=numel(unique(presentMets)) +0094 EM='Duplicate metabolite names in presentMets'; +0095 dispEM(EM); +0096 end +0097 +0098 %Default is that the metabolites cannot be produced +0099 if ~isempty(presentMets) +0100 metProduction=ones(numel(presentMets),1)*-2; +0101 presentMets=upper(presentMets); +0102 pmIndexes=find(ismember(presentMets,upper(model.metNames))); +0103 metProduction(pmIndexes)=-1; %Then set that they are at least found +0104 else +0105 metProduction=[]; +0106 pmIndexes=[]; +0107 end +0108 +0109 %The model should be in the reversible format and all relevant exchange +0110 %reactions should be open +0111 if isfield(model,'unconstrained') +0112 EM='Exchange metabolites are still present in the model. Use simplifyModel if this is not intended'; +0113 dispEM(EM,false); +0114 end +0115 +0116 %The irreversible reactions that are essential must have a flux and are +0117 %therefore not optimized for using MILP, which reduces the problem size. +0118 %However, reversible reactions must have a flux in one direction, so they +0119 %have to stay in the problem. The essentiality constraint on reversible +0120 %reactions is implemented in the same manner as for reversible reactions +0121 %when noRevLoops==true, but with the additional constraint that C ub=-1. +0122 %This forces one of the directions to be active. +0123 revRxns=find(model.rev~=0); +0124 essentialReversible=find(ismember(model.rxns(revRxns),essentialRxns)); +0125 essentialRxns=intersect(essentialRxns,model.rxns(model.rev==0)); +0126 +0127 %Convert the model to irreversible +0128 irrevModel=convertToIrrev(model); +0129 rxnScores=[rxnScores;rxnScores(model.rev==1)]; +0130 %These are used if noRevLoops is true +0131 if noRevLoops==true +0132 forwardIndexes=find(model.rev~=0); +0133 backwardIndexes=(numel(model.rxns)+1:numel(irrevModel.rxns))'; +0134 else +0135 %Then they should only be used for essential reversible reactions +0136 forwardIndexes=revRxns(essentialReversible); +0137 backwardIndexes=essentialReversible+numel(model.rxns); +0138 end +0139 +0140 %Get the indexes of the essential reactions and remove them from the +0141 %scoring vector +0142 essentialIndex=find(ismember(irrevModel.rxns,essentialRxns)); +0143 rxnScores(essentialIndex)=[]; +0144 +0145 %Go through each of the presentMets (if they exist) and modify the S matrix +0146 %so that each reaction which produces any of them also produces a +0147 %corresponding fake metabolite and the opposite in the reverse direction. +0148 +0149 %This is to deal with the fact that there is no compartment info regarding +0150 %the presentMets. This modifies the irrevModel structure, but that is fine +0151 %since it's the model structure that is returned. +0152 if any(pmIndexes) +0153 irrevModel.metNames=upper(irrevModel.metNames); +0154 metsToAdd.mets=strcat({'FAKEFORPM'},num2str(pmIndexes)); +0155 metsToAdd.metNames=metsToAdd.mets; +0156 metsToAdd.compartments=irrevModel.comps{1}; +0157 +0158 %There is no constraints on the metabolites yet, since maybe not all of +0159 %them could be produced +0160 irrevModel=addMets(irrevModel,metsToAdd); +0161 end +0162 +0163 %Modify the matrix +0164 for i=1:numel(pmIndexes) +0165 %Get the matching mets +0166 I=ismember(irrevModel.metNames,presentMets(pmIndexes(i))); +0167 +0168 %Find the reactions where any of them are used. +0169 [~, K, L]=find(irrevModel.S(I,:)); +0170 +0171 %This ugly loop is to avoid problems if a metabolite occurs several +0172 %times in one reaction +0173 KK=unique(K); +0174 LL=zeros(numel(KK),1); +0175 for j=1:numel(KK) +0176 LL(j)=sum(L(K==KK(j))); +0177 end +0178 irrevModel.S(numel(irrevModel.mets)-numel(pmIndexes)+i,KK)=LL; +0179 end +0180 +0181 %Some nice to have numbers +0182 nMets=numel(irrevModel.mets); +0183 nRxns=numel(irrevModel.rxns); +0184 nEssential=numel(essentialIndex); +0185 nNonEssential=nRxns-nEssential; +0186 nonEssentialIndex=setdiff(1:nRxns,essentialIndex); +0187 S=irrevModel.S; +0188 +0189 %Add so that each non-essential reaction produces one unit of a fake +0190 %metabolite +0191 temp=sparse(1:nRxns,1:nRxns,1); +0192 temp(essentialIndex,:)=[]; +0193 S=[S;temp]; +0194 +0195 %Add another set of reactions (will be binary) which also produce these +0196 %fake metabolites, but with a stoichiometry of 1000 +0197 temp=sparse(1:nNonEssential,1:nNonEssential,1000); +0198 temp=[sparse(nMets,nNonEssential);temp]; +0199 S=[S temp]; +0200 +0201 %Add reactions for net-production of (real) metabolites +0202 if prodWeight~=0 +0203 temp=[speye(nMets-numel(pmIndexes))*-1;sparse(nNonEssential+numel(pmIndexes),nMets-numel(pmIndexes))]; +0204 S=[S temp]; +0205 %To keep the number of reactions added like this +0206 nNetProd=nMets-numel(pmIndexes); +0207 else +0208 nNetProd=0; +0209 end +0210 +0211 %Add constraints so that reversible reactions can only be used in one +0212 %direction. This is done by adding the fake metabolites A, B, C for each +0213 %reversible reaction in the following manner +0214 % forward: A + .. => ... backwards: B + ... => ... int1: C => 1000 A int2: +0215 % C => 1000 B A ub=999.9 B ub=999.9 C lb=-1 int1 and int2 are binary +0216 if any(forwardIndexes) +0217 nRevBounds=numel(forwardIndexes); +0218 +0219 %Add the A metabolites for the forward reactions and the B metabolites +0220 %for the reverse reactions +0221 I=speye(numel(irrevModel.rxns))*-1; +0222 temp=[I(forwardIndexes,:);I(backwardIndexes,:)]; +0223 +0224 %Padding +0225 temp=[temp sparse(size(temp,1),size(S,2)-numel(irrevModel.rxns))]; +0226 +0227 %Add the int1 & int2 reactions that produce A and B +0228 temp=[temp speye(nRevBounds*2)*1000]; +0229 +0230 %And add that they also consume C +0231 temp=[temp;[sparse(nRevBounds,size(S,2)) speye(nRevBounds)*-1 speye(nRevBounds)*-1]]; +0232 +0233 %Add the new reactions and metabolites +0234 S=[S sparse(size(S,1),nRevBounds*2)]; +0235 S=[S;temp]; +0236 else +0237 nRevBounds=0; +0238 end +0239 +0240 %Add so that the essential reactions must have a small flux and that the +0241 %binary ones (and net-production reactions) may have zero flux. The integer +0242 %reactions for reversible reactions have [0 1] +0243 prob.blx=[irrevModel.lb;zeros(nNonEssential+nNetProd+nRevBounds*2,1)]; +0244 prob.blx(essentialIndex)=max(0.1,prob.blx(essentialIndex)); +0245 +0246 %Add so that the binary ones and net-production reactions can have at the +0247 %most flux 1.0 +0248 prob.bux=[irrevModel.ub;ones(nNonEssential+nNetProd+nRevBounds*2,1)]; +0249 +0250 %Add that the fake metabolites must be produced in a small amount and that +0251 %the A and B metabolites for reversible reactions can be [0 999.9] and C +0252 %metabolites [-1 0] +0253 prob.blc=[irrevModel.b(:,1);ones(nNonEssential,1);zeros(nRevBounds*2,1);ones(nRevBounds,1)*-1]; +0254 +0255 %Add that normal metabolites can be freely excreted if +0256 %allowExcretion==true, and that the fake ones can be excreted 1000 units at +0257 %most. C metabolites for essential reversible reactions should have an +0258 %upper bound of -1. If noRevLoops is false, then add this constraint for +0259 %all the reactions instead. +0260 if noRevLoops==true +0261 revUB=zeros(nRevBounds,1); +0262 revUB(essentialReversible)=-1; +0263 else +0264 revUB=ones(nRevBounds,1)*-1; +0265 end +0266 if allowExcretion==true +0267 metUB=inf(nMets,1); +0268 else +0269 metUB=irrevModel.b(:,min(size(irrevModel.b,2),2)); +0270 end +0271 prob.buc=[metUB;ones(nNonEssential,1)*1000;ones(nRevBounds*2,1)*999.9;revUB]; +0272 +0273 %Add objective coefficients for the binary reactions. The negative is used +0274 %since we're minimizing. The negative is taken for the prodWeight as well, +0275 %in order to be consistent with the syntax that positive scores are good +0276 prob.c=[zeros(nRxns,1);rxnScores;ones(nNetProd,1)*prodWeight*-1;zeros(nRevBounds*2,1)]; +0277 prob.a=S; +0278 +0279 % adapt problem structure for cobra-style solver +0280 prob.c=[prob.c;zeros(size(prob.a,1),1)]; +0281 prob.A=[prob.a -speye(size(prob.a,1))]; +0282 prob.b=zeros(size(prob.a,1), 1); +0283 prob.ub=[prob.bux; prob.buc]; +0284 prob.osense=1; +0285 prob.csense=char(zeros(1,size(prob.a,1))); +0286 prob.csense(:)='E'; +0287 +0288 %We still don't know which of the presentMets that can be produced. Go +0289 %through them, force production, and see if the problem can be solved +0290 for i=1:numel(pmIndexes) +0291 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=1; +0292 prob.lb=[prob.blx; prob.blc]; +0293 res=optimizeProb(prob,params); +0294 isFeasible=checkSolution(res); +0295 if ~isFeasible +0296 %Reset the constraint again +0297 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=0; +0298 else +0299 %Metabolite produced +0300 metProduction(pmIndexes(i))=1; +0301 end +0302 end +0303 prob.lb=[prob.blx; prob.blc]; +0304 +0305 %Add that the binary reactions may only take integer values. +0306 prob.vartype = repmat('C', 1, size(prob.A, 2)); +0307 allInt=[(nRxns+1):(nRxns+nNonEssential) size(S,2)-nRevBounds*2+1:size(S,2)]; +0308 prob.vartype(allInt) = 'B'; +0309 +0310 % solve problem +0311 res=optimizeProb(prob,params); +0312 +0313 %Problem should not be infeasible, but it is possible that the time limit +0314 %was reached before finding any solutions. +0315 if ~checkSolution(res) +0316 if strcmp(res.origStat, 'TIME_LIMIT') +0317 EM='Time limit reached without finding a solution. Try increasing the TimeLimit parameter.'; +0318 else +0319 EM='The problem is infeasible'; +0320 end +0321 dispEM(EM); +0322 end +0323 +0324 fValue=res.obj; +0325 +0326 %Get all reactions used in the irreversible model +0327 usedRxns=(nonEssentialIndex(res.full(nRxns+1:nRxns+nNonEssential)<0.1))'; +0328 +0329 %Map to reversible model IDs +0330 usedRxns=[usedRxns(usedRxns<=numel(model.rxns));revRxns(usedRxns(usedRxns>numel(model.rxns))-numel(model.rxns))]; +0331 +0332 %Then get the ones that are not used in either direction or is essential +0333 I=true(numel(model.rxns),1); +0334 I(usedRxns)=false; +0335 I(essentialIndex)=false; +0336 deletedRxns=model.rxns(I); +0337 outModel=removeReactions(model,I,true,true); +0338 end0001 function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params) 0002 % runINIT 0003 % Generates a model using the INIT algorithm, based on proteomics and/or 0004 % transcriptomics and/or metabolomics and/or metabolic tasks. This is the @@ -141,303 +142,304 @@SOURCE CODE
% problem significantly more computationally intensive to 0039 % solve (two more integer constraints per reversible reaction) 0040 % (optional, default false) -0041 % -0042 % outModel the resulting model structure -0043 % deletedRxns reactions which were deleted by the algorithm -0044 % metProduction array that indicates which of the -0045 % metabolites in presentMets that could be -0046 % produced -0047 % -2: metabolite name not found in model -0048 % -1: metabolite found, but it could not be produced -0049 % 1: metabolite could be produced -0050 % fValue objective value (sum of (the negative of) -0051 % reaction scores for the included reactions and -0052 % prodWeight*number of produced metabolites) -0053 % -0054 % This function is the actual implementation of the algorithm. See -0055 % getINITModel for a higher-level function for model reconstruction. See -0056 % PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the -0057 % implementation. -0058 % -0059 % Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,... -0060 % rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,... -0061 % noRevLoops,params) -0062 -0063 if nargin<2 -0064 rxnScores=zeros(numel(model.rxns),1); -0065 end -0066 if isempty(rxnScores) -0067 rxnScores=zeros(numel(model.rxns),1); -0068 end -0069 if nargin<3 || isempty(presentMets) -0070 presentMets={}; -0071 else -0072 presentMets=convertCharArray(presentMets); -0073 end -0074 if nargin<4 || isempty(essentialRxns) -0075 essentialRxns={}; -0076 else -0077 essentialRxns=convertCharArray(essentialRxns); -0078 end -0079 if nargin<5 || isempty(prodWeight) -0080 prodWeight=0.5; -0081 end -0082 if nargin<6 -0083 allowExcretion=false; -0084 end -0085 if nargin<7 -0086 noRevLoops=false; -0087 end -0088 if nargin<8 -0089 params=[]; -0090 end -0091 -0092 if numel(presentMets)~=numel(unique(presentMets)) -0093 EM='Duplicate metabolite names in presentMets'; -0094 dispEM(EM); -0095 end -0096 -0097 %Default is that the metabolites cannot be produced -0098 if ~isempty(presentMets) -0099 metProduction=ones(numel(presentMets),1)*-2; -0100 presentMets=upper(presentMets); -0101 pmIndexes=find(ismember(presentMets,upper(model.metNames))); -0102 metProduction(pmIndexes)=-1; %Then set that they are at least found -0103 else -0104 metProduction=[]; -0105 pmIndexes=[]; -0106 end -0107 -0108 %The model should be in the reversible format and all relevant exchange -0109 %reactions should be open -0110 if isfield(model,'unconstrained') -0111 EM='Exchange metabolites are still present in the model. Use simplifyModel if this is not intended'; -0112 dispEM(EM,false); -0113 end -0114 -0115 %The irreversible reactions that are essential must have a flux and are -0116 %therefore not optimized for using MILP, which reduces the problem size. -0117 %However, reversible reactions must have a flux in one direction, so they -0118 %have to stay in the problem. The essentiality constraint on reversible -0119 %reactions is implemented in the same manner as for reversible reactions -0120 %when noRevLoops==true, but with the additional constraint that C ub=-1. -0121 %This forces one of the directions to be active. -0122 revRxns=find(model.rev~=0); -0123 essentialReversible=find(ismember(model.rxns(revRxns),essentialRxns)); -0124 essentialRxns=intersect(essentialRxns,model.rxns(model.rev==0)); -0125 -0126 %Convert the model to irreversible -0127 irrevModel=convertToIrrev(model); -0128 rxnScores=[rxnScores;rxnScores(model.rev==1)]; -0129 %These are used if noRevLoops is true -0130 if noRevLoops==true -0131 forwardIndexes=find(model.rev~=0); -0132 backwardIndexes=(numel(model.rxns)+1:numel(irrevModel.rxns))'; -0133 else -0134 %Then they should only be used for essential reversible reactions -0135 forwardIndexes=revRxns(essentialReversible); -0136 backwardIndexes=essentialReversible+numel(model.rxns); -0137 end -0138 -0139 %Get the indexes of the essential reactions and remove them from the -0140 %scoring vector -0141 essentialIndex=find(ismember(irrevModel.rxns,essentialRxns)); -0142 rxnScores(essentialIndex)=[]; -0143 -0144 %Go through each of the presentMets (if they exist) and modify the S matrix -0145 %so that each reaction which produces any of them also produces a -0146 %corresponding fake metabolite and the opposite in the reverse direction. -0147 -0148 %This is to deal with the fact that there is no compartment info regarding -0149 %the presentMets. This modifies the irrevModel structure, but that is fine -0150 %since it's the model structure that is returned. -0151 if any(pmIndexes) -0152 irrevModel.metNames=upper(irrevModel.metNames); -0153 metsToAdd.mets=strcat({'FAKEFORPM'},num2str(pmIndexes)); -0154 metsToAdd.metNames=metsToAdd.mets; -0155 metsToAdd.compartments=irrevModel.comps{1}; -0156 -0157 %There is no constraints on the metabolites yet, since maybe not all of -0158 %them could be produced -0159 irrevModel=addMets(irrevModel,metsToAdd); -0160 end -0161 -0162 %Modify the matrix -0163 for i=1:numel(pmIndexes) -0164 %Get the matching mets -0165 I=ismember(irrevModel.metNames,presentMets(pmIndexes(i))); -0166 -0167 %Find the reactions where any of them are used. -0168 [~, K, L]=find(irrevModel.S(I,:)); -0169 -0170 %This ugly loop is to avoid problems if a metabolite occurs several -0171 %times in one reaction -0172 KK=unique(K); -0173 LL=zeros(numel(KK),1); -0174 for j=1:numel(KK) -0175 LL(j)=sum(L(K==KK(j))); -0176 end -0177 irrevModel.S(numel(irrevModel.mets)-numel(pmIndexes)+i,KK)=LL; -0178 end -0179 -0180 %Some nice to have numbers -0181 nMets=numel(irrevModel.mets); -0182 nRxns=numel(irrevModel.rxns); -0183 nEssential=numel(essentialIndex); -0184 nNonEssential=nRxns-nEssential; -0185 nonEssentialIndex=setdiff(1:nRxns,essentialIndex); -0186 S=irrevModel.S; -0187 -0188 %Add so that each non-essential reaction produces one unit of a fake -0189 %metabolite -0190 temp=sparse(1:nRxns,1:nRxns,1); -0191 temp(essentialIndex,:)=[]; -0192 S=[S;temp]; -0193 -0194 %Add another set of reactions (will be binary) which also produce these -0195 %fake metabolites, but with a stoichiometry of 1000 -0196 temp=sparse(1:nNonEssential,1:nNonEssential,1000); -0197 temp=[sparse(nMets,nNonEssential);temp]; -0198 S=[S temp]; -0199 -0200 %Add reactions for net-production of (real) metabolites -0201 if prodWeight~=0 -0202 temp=[speye(nMets-numel(pmIndexes))*-1;sparse(nNonEssential+numel(pmIndexes),nMets-numel(pmIndexes))]; -0203 S=[S temp]; -0204 %To keep the number of reactions added like this -0205 nNetProd=nMets-numel(pmIndexes); -0206 else -0207 nNetProd=0; -0208 end -0209 -0210 %Add constraints so that reversible reactions can only be used in one -0211 %direction. This is done by adding the fake metabolites A, B, C for each -0212 %reversible reaction in the following manner -0213 % forward: A + .. => ... backwards: B + ... => ... int1: C => 1000 A int2: -0214 % C => 1000 B A ub=999.9 B ub=999.9 C lb=-1 int1 and int2 are binary -0215 if any(forwardIndexes) -0216 nRevBounds=numel(forwardIndexes); -0217 -0218 %Add the A metabolites for the forward reactions and the B metabolites -0219 %for the reverse reactions -0220 I=speye(numel(irrevModel.rxns))*-1; -0221 temp=[I(forwardIndexes,:);I(backwardIndexes,:)]; -0222 -0223 %Padding -0224 temp=[temp sparse(size(temp,1),size(S,2)-numel(irrevModel.rxns))]; -0225 -0226 %Add the int1 & int2 reactions that produce A and B -0227 temp=[temp speye(nRevBounds*2)*1000]; -0228 -0229 %And add that they also consume C -0230 temp=[temp;[sparse(nRevBounds,size(S,2)) speye(nRevBounds)*-1 speye(nRevBounds)*-1]]; -0231 -0232 %Add the new reactions and metabolites -0233 S=[S sparse(size(S,1),nRevBounds*2)]; -0234 S=[S;temp]; -0235 else -0236 nRevBounds=0; -0237 end -0238 -0239 %Add so that the essential reactions must have a small flux and that the -0240 %binary ones (and net-production reactions) may have zero flux. The integer -0241 %reactions for reversible reactions have [0 1] -0242 prob.blx=[irrevModel.lb;zeros(nNonEssential+nNetProd+nRevBounds*2,1)]; -0243 prob.blx(essentialIndex)=max(0.1,prob.blx(essentialIndex)); -0244 -0245 %Add so that the binary ones and net-production reactions can have at the -0246 %most flux 1.0 -0247 prob.bux=[irrevModel.ub;ones(nNonEssential+nNetProd+nRevBounds*2,1)]; -0248 -0249 %Add that the fake metabolites must be produced in a small amount and that -0250 %the A and B metabolites for reversible reactions can be [0 999.9] and C -0251 %metabolites [-1 0] -0252 prob.blc=[irrevModel.b(:,1);ones(nNonEssential,1);zeros(nRevBounds*2,1);ones(nRevBounds,1)*-1]; -0253 -0254 %Add that normal metabolites can be freely excreted if -0255 %allowExcretion==true, and that the fake ones can be excreted 1000 units at -0256 %most. C metabolites for essential reversible reactions should have an -0257 %upper bound of -1. If noRevLoops is false, then add this constraint for -0258 %all the reactions instead. -0259 if noRevLoops==true -0260 revUB=zeros(nRevBounds,1); -0261 revUB(essentialReversible)=-1; -0262 else -0263 revUB=ones(nRevBounds,1)*-1; -0264 end -0265 if allowExcretion==true -0266 metUB=inf(nMets,1); -0267 else -0268 metUB=irrevModel.b(:,min(size(irrevModel.b,2),2)); -0269 end -0270 prob.buc=[metUB;ones(nNonEssential,1)*1000;ones(nRevBounds*2,1)*999.9;revUB]; -0271 -0272 %Add objective coefficients for the binary reactions. The negative is used -0273 %since we're minimizing. The negative is taken for the prodWeight as well, -0274 %in order to be consistent with the syntax that positive scores are good -0275 prob.c=[zeros(nRxns,1);rxnScores;ones(nNetProd,1)*prodWeight*-1;zeros(nRevBounds*2,1)]; -0276 prob.a=S; -0277 -0278 % adapt problem structure for cobra-style solver -0279 prob.c=[prob.c;zeros(size(prob.a,1),1)]; -0280 prob.A=[prob.a -speye(size(prob.a,1))]; -0281 prob.b=zeros(size(prob.a,1), 1); -0282 prob.ub=[prob.bux; prob.buc]; -0283 prob.osense=1; -0284 prob.csense=char(zeros(1,size(prob.a,1))); -0285 prob.csense(:)='E'; -0286 -0287 %We still don't know which of the presentMets that can be produced. Go -0288 %through them, force production, and see if the problem can be solved -0289 for i=1:numel(pmIndexes) -0290 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=1; -0291 prob.lb=[prob.blx; prob.blc]; -0292 res=optimizeProb(prob,params); -0293 isFeasible=checkSolution(res); -0294 if ~isFeasible -0295 %Reset the constraint again -0296 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=0; -0297 else -0298 %Metabolite produced -0299 metProduction(pmIndexes(i))=1; -0300 end -0301 end -0302 prob.lb=[prob.blx; prob.blc]; -0303 -0304 %Add that the binary reactions may only take integer values. -0305 prob.vartype = repmat('C', 1, size(prob.A, 2)); -0306 allInt=[(nRxns+1):(nRxns+nNonEssential) size(S,2)-nRevBounds*2+1:size(S,2)]; -0307 prob.vartype(allInt) = 'B'; -0308 -0309 % solve problem -0310 res=optimizeProb(prob,params); -0311 -0312 %Problem should not be infeasible, but it is possible that the time limit -0313 %was reached before finding any solutions. -0314 if ~checkSolution(res) -0315 if strcmp(res.origStat, 'TIME_LIMIT') -0316 EM='Time limit reached without finding a solution. Try increasing the TimeLimit parameter.'; -0317 else -0318 EM='The problem is infeasible'; -0319 end -0320 dispEM(EM); -0321 end -0322 -0323 fValue=res.obj; -0324 -0325 %Get all reactions used in the irreversible model -0326 usedRxns=(nonEssentialIndex(res.full(nRxns+1:nRxns+nNonEssential)<0.1))'; -0327 -0328 %Map to reversible model IDs -0329 usedRxns=[usedRxns(usedRxns<=numel(model.rxns));revRxns(usedRxns(usedRxns>numel(model.rxns))-numel(model.rxns))]; -0330 -0331 %Then get the ones that are not used in either direction or is essential -0332 I=true(numel(model.rxns),1); -0333 I(usedRxns)=false; -0334 I(essentialIndex)=false; -0335 deletedRxns=model.rxns(I); -0336 outModel=removeReactions(model,I,true,true); -0337 end
Generated by m2html © 2005